overestimated the effect of photon splitting; see Section 3.2). Baring & Harding (2007) model the process of resonant in- verse Compton scattering of photons from the neutron star surface (with the blackbody temperature T = 6 × 10 6 K) in the same field range, but only for single scattering events (see also Dermer 1990). The magnetosphere acceleration zone in the superstrong, twisted field regime of magnetars is investigated analytically by Beloborodov & Thompson (2007) for cascades occurring in the closed field line region of the magnetosphere and by Thompson (2008a,b) in the open field line region. 2002, 2008; Lyutikov 2007), while the high-energy pho- tons emitted in pair cascade models can reproduce the ob- served pulse profiles and phase-resolved spectra of gamma- ray pulsars once the three-dimensional emission geometry is taken into account (see, e.g., Romani & Yadigaroglu 1995; Cheng, Ruderman, & Zhang 2000; Dyks & Rudak 2003; Harding et al. 2008; Bai & Spitkovsky 2009). We note in passing that the dense pair plasma generated by this cas- cade also plays an important role in models of pulsar wind nebulae (see Arons 2007 for a review). 2 Zach Medin and Dong Lai The behavior of the pair cascade in the superstrong field regime (magnetic field strengths B > ∼ B$_{Q}$ ≡ 4 . 414 × 10 13 G) and its effect on emission from pulsars and magnetars is somewhat puzzling. For example, of the dozen-or-so ob- served magnetars, only two show pulsed radio emission, and it is of a completely different nature than the emis- sion from “standard” radio pulsars (e.g., the radio pulsa- tions are transient and appear to be correlated with strong X-ray outbursts from the magnetars; see Camilo et al. 2007, 2008). In contrast, several radio pulsars with in- ferred surface field strengths similar to those of magne- tars have been discovered (e.g., Kaspi & McLaughlin 2005; Vranevsevic, Manchester, & Melrose 2007). Why the stan- dard mechanism for pulsed radio emission turns off for mag- netars but not for these pulsars is unknown. In this paper we present numerical simulations of the pair cascade from onset to completion. Motivated by the lack of full cascade results for the superstrong field regime, and in light of the unexplained differences between the ob- served emission properties of high-field radio pulsars and magnetars, we run our simulations in magnetospheres with field strengths up to 10 15 G. We consider several impor- tant factors that affect high-field cascades, including pho- ton splitting, pair creation in low Landau levels, photon polarization modes ( ⊥ or ‖ to the magnetic field direc- tion), and resonant inverse Compton scattering. We use our simulations to generate spectra of the high-energy pho- tons and the electron-positron plasma produced by the cas- cade. Additionally, we use our simulations to comment on the conditions for when the radio emission mechanism no longer operates in the neutron star magnetosphere, the so- called “pulsar death line” (e.g., Ruderman & Sutherland 1975; Chen & Ruderman 1993; Hibschman & Arons 2001a; Harding & Muslimov 2002; Harding, Muslimov, & Zhang 2002; Medin & Lai 2007). While the results of our simulation are most applicable to cascades occurring in the open field line region of the magnetosphere (since the primary particles are injected into the magnetosphere along open field lines), some of our results are also relevant to cascades occurring in the closed field line region for magnetars, e.g., the products of a cascade initiated by a photon injected into a non-dipole magnetosphere. A necessary component of any pair cascade simulation is a model of the magnetosphere acceleration zone, or “gap”, where the cascade originates. In real magnetospheres of pul- sars and magnetars, the acceleration of primary particles is coupled to the rest of the cascade (e.g., charged particles pro- duced in the cascade can screen out the acceleration poten- tial). However, there is significant uncertainty about the pre- cise nature of the acceleration gap. A number of models have been proposed for the location of the gap, from inner mag- netosphere accelerators (both “vacuum” and “space-charge- limited flow” types; see, e.g., Ruderman & Sutherland 1975; Arons & Scharlemann 1979; Muslimov & Tsygan 1992; Hibschman & Arons 2001a; Medin & Lai 2007; Thompson 2008a,b), to outer magnetosphere accelerators (e.g., Cheng et al. 1986a,b; Romani 1996; Cheng et al. 2000; Takata et al. 2006), to hybrid inner-outer magnetosphere ac- celerators (“slot” gaps and extended outer gaps; e.g., Arons 1983; Muslimov & Harding 2003, 2004; Hirotani 2006). Non-steady (oscillatory) inner gaps have also been dis- cussed recently (e.g., Sakai & Shibata 2003; Levinson et al. 2005; Beloborodov 2008; Luo & Melrose 2008). Numer- ical simulations of force-free global magnetospheres in- cluding magnetic-field twisting near the light cylinder have been performed (e.g., Contopoulos, Kazanas, & Fendt However, for superstrong magnetic fields ( B > ∼ B$_{Q}$ ≡ 4 . 414 × 10 13 G) only limited aspects of the full cascades have been studied. For example, Arendt & Eilek (2002) simulate the cascade for B$_{p}$ ≤ 10 13 G and P = 0 . 033 s (for both a pure dipole and a more complex field geometry), but with the simplification that all photons radiated by the primary particle are emitted from the surface. Baring & Harding (2001) (see also Harding, Baring, & Gonthier 1997) use this same simplification to study the effects of photon splitting on the cascade for field strengths up to B = 2 × 10 14 G (however, they assumed that both photon modes can split, and thus There have been only a few publications devoted to numerical simulations of the pair cascade in pul- sar magnetospheres. For moderate-strength magnetic fields ( B < ∼ 5 × 10 12 G), significant progress has been made. Daugherty & Harding (1982) present simulations of the cas- cade initiated by a single electron injected from the neu- tron star surface, emitting photons through curvature ra- diation, for (polar) surface field strengths B$_{p}$ up to 5 × 10 12 G and rotation periods P = 0 . 033–1 s. In a later paper (Daugherty & Harding 1996) they consider gamma ray emission from the entire open-field-line region of the magnetosphere, using a simplified acceleration model and for Vela-like pulsar parameters ( B$_{p}$ = 3 × 10 12 G and P = 0 . 089 s). Sturner, Dermer, & Michel (1995) present a similar simulation to that of Daugherty & Harding, but for cascades initiated by electrons upscattering photons through the inverse Compton process (again for Vela-like parameters). Hibschman & Arons (2001b) develop a semi- analytic model of the inner gap cascade, both for curva- ture radiation-initiated and inverse Compton scattering- initiated cascades, applicable for B < ∼ 3 × 10 12 G (see also Zhang & Harding 2000). Cascades occurring in the outer magnetosphere have also been simulated, by Romani (1996) for Vela- and Crab-like ( B$_{p}$ = 4 × 10 12 G and P = 0 . 033 s) parameters (see also Cheng, Ho, & Ruderman 1986a,b; Cheng, Ruderman, & Zhang 2000).