Three-dimensional Kinetic Pulsar Magnetosphere Models: Connecting to Gamma-Ray Observations

It has been half a century since the first pulsar was observed (Hewish et al. 1968). Pulsars are identified as rapidly rotating Neutron Stars (NSs) with huge surface magnetic fields B_⋆ that reach up to ≈10¹⁴ G. However, they continuously lose rotational energy while they radiate over almost the entire spectrum (from radio to γ-rays). The corresponding brightness temperature of radio emission is extremely high (10²³–10²⁵ K), indicating its coherent nature. The exact mechanism of radio emission remains unknown, although it is believed that it is related to pair-production processes.

However, the radio emission is energetically unimportant in comparison to the high-energy (e.g., γ-rays) emission. For decades, we had very limited information about the pulsar γ-ray emission, and it was only after the launch of Fermi in 2008 when we gradually gained access to a plethora of observational data. Thus, we now have more than 205 detected γ-ray pulsars (117 of them are compiled in the second pulsar catalog (2PC); Abdo et al. 2013), which has led to the derivation of meaningful statistical trends and correlations. This ushered in the Fermi era in pulsar research, significantly affecting the theoretical modeling of pulsar high-energy emission.

Pulsars are considered to be spherical perfect conductors that rotate within their own magnetic fields. This rotation polarizes the charges, supporting an electric field inside the conductor star,

$$\begin{eqnarray}&&{\boldsymbol{E}}=-\displaystyle \frac{1}{c}({\boldsymbol{\Omega }}\times r)\times {\boldsymbol{B}},\end{eqnarray} \tag{ 1 }$$

where Ω is the angular frequency of the star. Assuming a dipole magnetic field and no charges outside the stellar surface, the field structure is provided by the analytic Vacuum-retarded Dipole (VRD) solution (Deutsch 1955).

Soon after the discovery of pulsars, it became evident that the high voltages due to the accelerating electric field components ${{\boldsymbol{E}}}_{\mathrm{acc}}$, which exist in VRD solutions, are able to initiate pair cascade processes (Sturrock 1971). Actually, it was believed that the high pair-creation efficiency (especially of the B–γ process) would short-out ${{\boldsymbol{E}}}_{\mathrm{acc}}$ everywhere in the magnetosphere, leading to the so-called force-free (FF) regime. Thus, pulsars' radiation can be interpreted as the result of their (not totally successful) efforts to construct perfect (conductive) wires using their huge fields and microphysical processes.

Even though the main principles of the FF solutions had been described in the late 1960s and early 1970s (e.g., Goldreich & Julian 1969; Scharlemann & Wagoner 1973), the first realistic FF solution for pulsar magnetospheres was presented only two decades ago.

Contopoulos et al. (1999) presented the first FF solution with a dipolar magnetic field for an aligned rotator (α = 0°, where α is the inclination angle; see also Gruzinov 2005; Komissarov 2006; McKinney 2006; Timokhin 2006; Parfrey et al. 2012; Cao et al. 2016a), while Spitkovsky (2006) presented that for an oblique rotator (see also Kalapotharakos & Contopoulos 2009; Pétri 2012a; Tchekhovskoy et al. 2013). These solutions provided not only the field structure, but also the corresponding current and charge densities, and revealed that the current closure is achieved by the return current, the main part of it flowing along the equatorial current sheet (ECS) outside the light cylinder (LC). The ECS is stable up to large distances (Kalapotharakos et al. 2012a), while its 3D shape is close to the solution for the current sheet of an FF split monopole solution (Bogovalov 1999).

The FF solutions reveal the currents and field structure that make ${{\boldsymbol{E}}}_{\mathrm{acc}}=0$ everywhere in the magnetosphere. However, due to their ideal character, they do not provide any information about the acceleration regions, the corresponding emission, and apparently the related ${{\boldsymbol{E}}}_{\mathrm{acc}}$. The only way to study the high-energy emission in FF models was to place arbitrarily the emission regions based on the pulsar radiation models. Fermi observations of the high-energy cutoff in the spectrum of the Vela pulsar (Abdo et al. 2009) ruled out polar cap emission models (Arons & Scharlemann 1979; Daugherty & Harding 1982, 1996), which would produce a super-exponential spectral cutoff from magnetic pair attenuation. Thus, the candidate accelerating regions were either (i) along the last open field line, like the Slot Gap (SG; Arons 1983; Muslimov & Harding 2003, 2004), or (ii) near the Outer Gap (OG) region (Cheng et al. 1986; Romani 1996; Hirotani & Shibata 2001), or (iii) in the ECS (Lyubarskii 1996; Bai & Spitkovsky 2010; Contopoulos & Kalapotharakos 2010; Pétri 2012b; Arka & Dubus 2013), and also (iv) in the striped pulsar wind (Pétri & Kirk 2005).

Kalapotharakos et al. (2012c), Li et al. (2012), and more recently, Cao et al. (2016b), motivated by the aforementioned limitations of the FF solutions, started exploring resistive/dissipative solutions. These models, using rather simple macroscopic prescriptions for the current density J, regulate the ${{\boldsymbol{E}}}_{\mathrm{acc}}$ through a macroscopic conductivity σ. Varying σ from 0 to $\infty$, an entire spectrum of solutions from VRD to FF is produced, respectively. The physical meaning of σ is supposedly related to the local pair multiplicity and its efficiency in "killing" ${{\boldsymbol{E}}}_{\mathrm{acc}}$.

These dissipative solutions provided the spatial distribution of ${{\boldsymbol{E}}}_{\mathrm{acc}}$ and therefore, allowed for the integration of test particle trajectories in these models and consequently the study of the corresponding emission. Thus, assuming emission due to curvature radiation (CR), Kalapotharakos et al. (2012b) produced the first γ-ray light curves based on dissipative solutions, while Kalapotharakos et al. (2014) found that a hybrid conductivity model that has infinite σ within the LC and high and finite σ outside the LC reproduced the radio lag, δ, versus peak separation, Δ, correlation shown in 2PC. In the so-called FIDO (FF Inside, Dissipative Outside) models, the emission is produced near the ECS, though the emissivity is not uniformly distributed.

The successful fit of the 2PC δ–Δ correlation gave confidence that the FIDO macroscopic models provide viable γ-ray emission geometries. However, a successful high-energy emission model should be able to reproduce the observed spectral properties. Thus, in Brambilla et al. (2015), we used the FIDO models with approximate ${E}_{\mathrm{acc}}$ values at low σ in order to fit eight bright pulsars that have published phase-resolved spectra. This study indicated an increase of σ with the spin-down power $\dot{{ \mathcal E }}$ and a decrease with the pulsar age.

Guided by Fermi observations, in Kalapotharakos et al. (2017), we expanded our studies and refined the FIDO models by restricting the dissipative regions outside the LC only near the ECS. Running series of models that covered the entire range of the observed $\dot{{ \mathcal E }}$ of Young Pulsars (YPs) and Millisecond Pulsars (MPs), and matching the corresponding cutoff energies ${\epsilon }_{\mathrm{cut}}$, we revealed a dependence of σ with $\dot{{ \mathcal E }}$. We found that σ increases with $\dot{{ \mathcal E }}$ for high $\dot{{ \mathcal E }}$, while it saturates for low $\dot{{ \mathcal E }}$. We also found clear indications that the size of the dissipative zones increases toward low $\dot{{ \mathcal E }}$ and that the multiplicity of the emitting particles increases with $\dot{{ \mathcal E }}$ and α.

Recently, a new approach to modeling pulsar magnetospheres has been explored. A few groups have attempted to build kinetic models of pulsar magnetospheres using the Particle-in-cell (PIC) technique, wherein charged particle trajectories and electromagnetic fields are treated self-consistently (Chen & Beloborodov 2014; Philippov & Spitkovsky 2014; Belyaev 2015a, 2015b; Cerutti et al. 2015, 2016; Philippov et al. 2015a, 2015b; Brambilla et al. 2018). Some of these studies (Chen & Beloborodov 2014; Belyaev 2015b; Philippov et al. 2015a) tried to simulate the pair-creation physics but failed to produce magnetospheres of low inclination angles filled with plasma close to the FF ones. One of the problems in these studies is the poor numerical resolution inherent to any global model which might prevent an accurate modeling of the microphysics of particle acceleration and pair production. Cerutti et al. (2016) studied the high-energy emission from their PIC simulations and found that the main source is synchrotron radiation but with fixed (with $\dot{{ \mathcal E }}$) γ-ray efficiency (not totally consistent with observations; 2PC). More recently, Philippov & Spitkovsky (2018) agreed that higher particle multiplicities (i.e., higher σ) from pair creation in the outer magnetosphere is needed at higher $\dot{{ \mathcal E }}$ values to explain the dependence of the total γ-ray luminosity on $\dot{{ \mathcal E }}$ (i.e., ${L}_{\gamma }\propto {\dot{{ \mathcal E }}}^{1/2}$).

In the present paper, we present our PIC code, and we take the first step by applying ad hoc (not physically motivated) particle injections that provide at least field configurations and particle distributions that are consistent with each other. The main goal of this study is to explore the role of the global particle injection rate ${ \mathcal F }$ in both the development of the field structures toward the FF regime and the characteristics of the corresponding high-energy emission. More specifically, taking into account the constraints imposed by Fermi data, we reveal a relation between ${ \mathcal F }$ and $\dot{{ \mathcal E }}$ that is able to reproduce a broad spectrum of the observed phenomenology (i.e., γ-ray light curves and spectra).

The structure of this paper is as follows. In Section 2, we introduce our PIC code, discuss technical details, and present some basic tests. In Section 3, we discuss how we build our models and the adopted model parameters. In Section 4, we show a development of solutions from the VRD to the FF one as a function of the particle injection rate, ${ \mathcal F }$. In Section 5, we discuss the model interpretation and the method we use in order to scale the particle energies as well as the corresponding emitted photon energies to those of real pulsars. In Section 6, we present our main results for two series of models that cover the observed $\dot{{ \mathcal E }}$ values of YPs and MPs. For these models, we show the dependence of the spectral cutoff energies on ${ \mathcal F }$ which, by taking into account the variation of the Fermi ${\epsilon }_{\mathrm{cut}}$, leads to a dependence of ${ \mathcal F }$ on $\dot{{ \mathcal E }}$. Finally, in Section 7, we summarize our conclusions and discuss future prospects.

We have developed a Cartesian 3D PIC code for Astrophysical studies (hereafter C-3PA). Our code follows the well-known, simple, but also powerful algorithm:

• 1.  
  Integrate the time-dependent Maxwell's equations one time step using the current provided by the particle collective motions;
• 2.  
  Integrate the equations of motion one time step for all the particles, taking into account the field structure;
• 3.  
  Go to 1.

Our code is written in Fortran 90 and is fully parallelized through Message Passing Interface (MPI). The code is charge conservative, and in the current version, it uses the so-called cloud-in-cell scheme for the shape function (i.e., the shape of each macroparticle is a cube of size equal to one grid cell of the lattice; see Villasenor & Buneman 1992). For the particle mover, we implemented Vay's algorithm (Vay 2008), which, in general, provides better accuracy for the drift motion of relativistic particle trajectories (compared to the standard Boris algorithm; see Birdsall & Langdon 1991 and references therein). The forces exerted on the particles are weighted over the individual cubic particle shapes. The field solver integrates the time-dependent Maxwell's equations through a finite-difference time-domain method applied in a Yee staggered mesh. In the outer boundary, we use a Perfectly Matched Layer (PML; Berenger 1996; Kalapotharakos & Contopoulos 2009) that absorbs the outgoing electromagnetic waves and minimizes the inward reflections. PML is applied outside a central cubic domain with side ≈8.6 R_LC (where R_LC is the light-cylinder radius) and its width is ≈0.5 R_LC. The particles are removed from the simulation after they enter the outer layer. They are also removed when their entire shape enters the central spherical conductor (rotating NS). Moreover, a Gaussian kernel can be applied for the current density that mimics a higher-order shape function, considerably reducing the noise level.

Inside the stellar surface, the electric field is always defined by Equation (1). The problem with this configuration is that it does not incorporate the current closure of the charge carriers that reach the stellar surface. This leads to the charging of the surface mostly in regions the polar cap passes through (see also Spitkovsky & Arons 2002). The problem is more prominent for low α values because the corresponding polar caps do not significantly change their position on the spherical surface. Philippov & Spitkovsky (2014), in order to resolve this problem, simulated the behavior of the spherical conductor by enforcing charge motions inside the conductor that restore/sustain the electric field of Equation (1). We have implemented similar methods where the particles in a thin layer at the stellar surface are enforced to move in a way that restores the conductive electric field and the continuity of the parallel (to the stellar surface) electric field component. Nonetheless, in our simulations, we have adopted a much simpler and less computationally expensive method that provides very similar results. Namely, we integrate Maxwell's equations down to the stellar radius, while at the end of each time step, we enforce E, within a thin layer (∼3 cells wide) outside the stellar surface, to go to the conductive value (i.e., Equation (1)). This treatment "cleans" the traces of the charges that enter the stellar surface, and it actually mimics the current that maintains the electric field in the conductor. In this approach, the outer surface of this thin layer determines the effective stellar surface.

We have also incorporated into the particle equations of motion the radiation reaction forces. The complete expression reads (Landau & Lifshitz 1987)

$$\begin{eqnarray}{{\boldsymbol{F}}}_{\mathrm{rr}} & = & \displaystyle \frac{2{q}_{{\rm{e}}}^{4}}{3{m}_{{\rm{e}}}^{2}{c}^{4}}\left[{\boldsymbol{E}}\times {\boldsymbol{B}}+{\boldsymbol{B}}\times \left({\boldsymbol{B}}\times \displaystyle \frac{{\boldsymbol{v}}}{{\boldsymbol{c}}}\right)+{\boldsymbol{E}}\left(\displaystyle \frac{{\boldsymbol{v}}}{c}\cdot {\boldsymbol{E}}\right)\right]\\ & & -\displaystyle \frac{2{q}_{{\rm{e}}}^{4}{\gamma }_{{\rm{L}}}^{2}}{3{m}_{{\rm{e}}}^{2}{c}^{5}}{\boldsymbol{v}}\left[{\left({\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{v}}}{c}\times {\boldsymbol{B}}\right)}^{2}-{\left({\boldsymbol{E}}\cdot \displaystyle \frac{{\boldsymbol{v}}}{c}\right)}^{2}\right]\\ & & +\displaystyle \frac{2{q}_{{\rm{e}}}^{3}{\gamma }_{{\rm{L}}}}{3{m}_{{\rm{e}}}{c}^{3}}\left[\left(\displaystyle \frac{\partial }{\partial t}+{\boldsymbol{v}}\cdot {\rm{\nabla }}\right){\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{v}}}{c}\times \left(\displaystyle \frac{\partial }{\partial t}+{\boldsymbol{v}}\cdot {\rm{\nabla }}\right){\boldsymbol{B}}\right],\end{eqnarray} \tag{ 2 }$$

where q_e and m_e are the electron charge and mass, respectively, and v and γ_L are the corresponding particle velocity and Lorentz factor. The third term involves the convective derivatives of the fields, and besides the implementation difficulty, especially in an evolving system, we have found (by applying it to stationary magnetosphere solutions) that it is negligible compared to the first two terms (see also Tamburini et al. 2010; Cerutti et al. 2016). Thus, in C-3PA we have incorporated only the first two terms following the numerical scheme suggested in Tamburini et al. (2010).

In almost every PIC code, the run time is determined mainly by the work load that corresponds to the integration of the particle equations of motion (and not the integration of Maxwell's equations). In our problem, the spatial particle distribution is quite non-uniform. The particle number density outside the stellar surface is much higher than the one at the edges of the computational domain. Assuming an MPI implementation where all of the central processing units (CPUs) "control" equal computational volumes (Figure 1(a)), the performance is slower (in several cases by more than one order of magnitude) compared to the theoretical one that corresponds to the total number of operating CPUs. This happens because the CPUs that "control" the outer parts of the domain have to wait for the inner CPUs to finish the calculations for relatively much higher numbers of particles. C-3PA takes care of this load balance issue by implementing a rather simple non-uniform volume distribution (Figure 1(b)) that, in general, increases the controlled volume with the distance from the center. The advantage of the scheme shown in Figure 1(b) is that it maintains straightforward communication between neighbor cuboids. A proper choice of the different length sizes of the rectangular cuboids (rectangles in Figure 1(b)) considerably improves the computational speed. In particular, for simulations with 16³ = 4096 CPUs, we gain a factor ∼7.

Zoom In Zoom Out Reset image size

Besides all trivial tests, we have reproduced the two-stream instability (Birdsall & Langdon 1991). It is well-known that two sets of opposite streams of charged particles are unstable. Any small charge imbalance (i.e., perturbation) creates an electric field along the motion direction that grows exponentially. In Figure 2(a) we show, for a configuration of two initially opposite streams of e⁻ within a proton background, the evolution of the amplitude of the fastest growing mode of the electric field component E_x along the motion direction x. The blue line corresponds to the simulated data, while the red line indicates the theoretical growth rate, which is equal to ω_p/2, where ω_p is the electron plasma frequency. Figures 2(b) and (c) show the phase space for the indicated times (see also the dashed light gray vertical lines in Figure 2(a)) that correspond to a snapshot during the exponential growth and to a snapshot when the instability has been fully developed. Moreover, we have checked for the Weibel type of instability. In the bottom row of Figure 2, we show the results corresponding to a 3D simulation, the initial state of which consists of uniformly distributed e⁻–e⁺ pairs that have opposite velocities along the x axis with γ_L = 15. The total net current is 0 since half of the pairs have e⁻ (e⁺) moving along the positive (negative) direction of the x axis while the other half has e⁻ (e⁺) moving along the corresponding negative (positive) direction. We note that the plasma is cold in the perpendicular direction (${\gamma }_{{{\rm{L}}}_{\perp }}=1$). In Figure 2(d), we plot, in the indicated color scale, the values of the component B_z of the magnetic field at $t\approx 22{\omega }_{p}^{-1}$, where ω_p is the fundamental (i.e., non-relativistic) plasma frequency. In Figure 2(e), the blue line shows the time evolution of the amplitude of the k_max mode of the B_z magnetic field component, where k_max is the wave number corresponding to the highest growth rate (Medvedev & Loeb 1999; Shaisultanov et al. 2012; Kumar et al. 2015). The red line denotes the theoretical growth rate for the k_max mode. We see that C-3PA reproduces not only the corresponding filamentary field structure but also the correct growth rate. We note also that saturation occurs when B reaches subequipartition values ∼0.1 (Medvedev & Loeb 1999).

Zoom In Zoom Out Reset image size

The present study focuses on exploring kinetic (PIC) pulsar magnetosphere models that ignore pair-creation microphysics. As mentioned above, the primary goal is to find how the solutions and the related properties depend on a single parameter, the global particle injection rate ${ \mathcal F }$. Thus, for the models presented throughout this paper, particle injection is made based on the following prescription: at each time step and at each cell up to r = 2.5 R_LC, one pair (e⁺, e⁻) is injected at rest (γ_L = 1; γ_L is the Lorentz factor) as long as the plasma magnetization,

$$\begin{eqnarray}&&{\sigma }_{{\rm{M}}}=\displaystyle \frac{{B}^{2}}{8\pi ({n}_{{e}^{+}}+{n}_{{e}^{-}}){m}_{{\rm{e}}}{c}^{2}},\end{eqnarray} \tag{ 3 }$$

exceeds a locally predefined value, Σ, that varies according to

$$\begin{eqnarray}{\rm{\Sigma }}=\left\{\begin{array}{ll}{{\rm{\Sigma }}}_{0}{\left(\tfrac{{r}_{{\rm{s}}}}{r}\right)}^{3} & {\rm{if}}\,{\rm{}}r\leqslant {R}_{\mathrm{LC}}\\ {{\rm{\Sigma }}}_{0}{\left(\tfrac{{r}_{{\rm{s}}}}{{R}_{\mathrm{LC}}}\right)}^{3}\tfrac{{R}_{\mathrm{LC}}}{r} & {\rm{if}}\,{\rm{}}r\gt {R}_{\mathrm{LC}}\end{array}\right.,\end{eqnarray} \tag{ 4 }$$

where m_e is the electron mass, ${n}_{{e}^{-}}$ and ${n}_{{e}^{+}}$ are the ${e}^{-}$ and e⁺ number densities, r_s is the stellar radius in the simulation, and Σ₀ is the value of Σ at r = r_s. The adopted dependence of Σ on r is similar to the B dependence in the FF solutions. Assuming n ∝ B, then σ_M ∝ B. Thus, Equation (4) provides ideally the same σ_M value everywhere. Even though this is an oversimplification (n is not proportional to B everywhere in the magnetosphere), we use this threshold prescription because it places appropriate amounts of the injected plasma (compared to the total/global injection rate) in a large area of the magnetosphere while avoiding over-injection of particles (e.g., in the closed field line regions). In that sense, Equation (4) is not essential, and a somewhat different radial profile would provide similar results.

Our code starts with some initial (estimated) Σ₀ value which is adjusted in time so that ${ \mathcal F }$ reaches the originally adopted goal value. The ${ \mathcal F }$ unit corresponds to the Goldreich–Julian flux from both polar caps, ${{ \mathcal F }}_{\mathrm{GJ}}={{ \mathcal F }}_{\mathrm{GJ}}^{0}\cos (\alpha )$,⁵ where ${{ \mathcal F }}_{\mathrm{GJ}}^{0}\,={B}_{{\rm{s}}}{\rm{\Omega }}{A}_{\mathrm{pc}}/(\pi {q}_{{\rm{e}}})$ (q_e is the elementary electric charge, B_s is the stellar surface magnetic field, and A_pc is the corresponding polar cap area) is the ${{ \mathcal F }}_{\mathrm{GJ}}$ value for α = 0°.

Starting from a configuration close to the VRD, we start injecting particles based on the aforementioned prescription, and let our simulations evolve for at least 1.5 periods while steady state is achieved after ∼1 period.

The stellar period in the simulation is P_s = 0.1 s, which places the LC at R_LC = 4.8 × 10⁸ cm. The adopted magnetic field at r_⋆ = 10 km (actual NS radius) is B_s = 10⁶ G, which allows us to achieve high magnetization values even for the highest injection rates. For the simulations with the highest particle injection rates, the magnetization is ∼500 at the stellar surface and ∼40 near the LC, and it never falls below 10 in the regions⁶ we present our results. We have also run simulations with B_s = 10⁵ G, and the results presented below remained unaffected. If not mentioned otherwise, the simulation resolution presented below is 25 grid points per R_LC. The simulated stellar radius is at r_s = 0.28 R_LC (resolved by 7 grid points) while the kernel layer thickness is 2.5 grid cells, which places the effective stellar radius at r_eff = 0.38 R_LC. We note that these stellar radii are unrealistically high. Nonetheless, many previous studies (Spitkovsky 2006; Kalapotharakos & Contopoulos 2009; Kalapotharakos et al. 2012a, 2012c; Li et al. 2012; Tchekhovskoy et al. 2013; Philippov & Spitkovsky 2014; Cerutti et al. 2015; Philippov et al. 2015a) have shown that the field and current structures on the polar caps remain practically unaffected as long as the stellar radius is not close to the LC.

The charge q_M of the individual macroparticles determines, for each ${ \mathcal F }$ value, the average number N_Mc of macroparticles per cell. The smaller the q_M (i.e., the higher N_Mc), the lower the noise level (i.e., the field and current fluctuations). However, we have found that the higher ${ \mathcal F }$ is, the lower q_M should be in order to keep the noise at relatively low levels. As we also discuss in Section 5 for solutions near the FF (i.e., high ${ \mathcal F }$), the actual accelerating electric fields become small, and the noisy fields can become comparable. Apparently, this implies that for higher ${ \mathcal F }$, we need disproportionately high numbers of macroparticles, which makes the corresponding simulations very expensive.

The time step Δt_s = 4 × 10⁻⁵ P_s guarantees that, for the densest magnetosphere regions of the models with the highest injection rates presented in this paper, (ω_pΔt_s)⁻¹ > 6, where ω_p is the fundamental plasma frequency ${\omega }_{{\rm{p}}}\,=\sqrt{4\pi ({n}_{{e}^{+}}+{n}_{{e}^{-}}){q}_{{\rm{e}}}^{2}/{m}_{{\rm{e}}}}$. In some cases, the adopted global time step does not resolve the gyromotion. In these cases, the particle equations of motion are integrated in smaller time steps. Following the orbits finely does not affect the currents that depend only on the motion within the global time step. Nonetheless, we have found that it provides more accurate particle energies, and as we discuss in Section 5, allows us to adjust the particle gyromotion.

We note that in some cases the adopted regular spatial resolution does not resolve the skin depth λ_D = c/ω_p for the low-γ_L particles in the densest magnetosphere regions. This implies that this particle population could be artificially heated to energies higher than expected physically. However, the corresponding thermal γ_L values are much smaller (∼1 order of magnitude) than that corresponding to the full potential drop across the polar cap (i.e., γ_Lmax = q_eB_sr_⋆³/m_ec² R_LC²) while the number of these particles seems to be small and does not affect the particle energy distributions. Actually, we have found that the particle energy distribution is more sensitive to the q_M and Δt_s values discussed above. Nonetheless, we have run simulations for the highest ${ \mathcal F }$ values, doubling the spatial resolution, and the particle energy distributions and all results presented in the next sections remained unaffected.

We present simulations for α = 15°, α = 45°, and α = 75°. For each of these α values, we have simulated 11 ${ \mathcal F }$ values that effectively produce an entire spectrum of solutions from close to VRD to close to FF. Table 1 shows the adopted ${ \mathcal F }$ values in ${{ \mathcal F }}_{\mathrm{GJ}}^{0}$ units (second column) and in the corresponding ${{ \mathcal F }}_{\mathrm{GJ}}$ units (third to fifth columns).

Table 1.  Models Corresponding to 11 ${ \mathcal F }$ Values

	α = 15°	α = 45°	α = 75°
${ \mathcal F }({{ \mathcal F }}_{\mathrm{GJ}}^{0})$	${ \mathcal F }({{ \mathcal F }}_{\mathrm{GJ}})$	${ \mathcal F }({{ \mathcal F }}_{\mathrm{GJ}})$	${ \mathcal F }({{ \mathcal F }}_{\mathrm{GJ}})$
0.1	0.1	0.14	0.38
0.25	0.26	0.35	0.97
0.6	0.62	0.85	2.3
1.5	1.55	2.1	5.8
3.5	3.6	4.9	13.5
6.5	6.7	9.2	25.1
8.5	8.8	12.0	32.8
11.5	11.9	16.2	44.4
15.0	15.5	21.2	58.0
20.0	20.7	28.3	77.3
28.0	29.0	39.6	108.2

Note. In the first column, the adopted ${ \mathcal F }$ values are shown in ${{ \mathcal F }}_{\mathrm{GJ}}^{0}$ units (i.e., the particle injection rate that corresponds to the GJ flux from the two polar caps of the aligned, α = 0°, rotator). The second to fourth columns show the adopted ${ \mathcal F }$ values in ${{ \mathcal F }}_{\mathrm{GJ}}$ units (i.e., the particle injection rate that corresponds to the GJ flux from the two polar caps for the indicated α value).

Download table as:  ASCIITypeset image

Figure 3 presents, for the indicated α values, the spin-down power, $\dot{{ \mathcal E }}$, versus the particle injection rate, ${ \mathcal F }$, in a log-linear plot. We see that for low and high ${ \mathcal F }$-values $\dot{{ \mathcal E }}$ asymptotically goes to the values corresponding to VRD and FF ones denoted by the dashed lines.

Zoom In Zoom Out Reset image size

In Figure 4, we see the field structure, for α = 45° and for the indicated ${ \mathcal F }$ values, on the poloidal ${\boldsymbol{\mu }}\mbox{--}{\boldsymbol{\Omega }}$ plane. The charge density ρ (left-hand column) as well as the number density of positrons ${n}_{{e}^{+}}$ (middle column) and electrons ${n}_{{e}^{-}}$ (right-hand column) are shown in the indicated color scales. As ${ \mathcal F }$ increases, features of the FF solution start developing. Thus, the magnetic field lines open gradually beyond the LC and become straight, while the ECS beyond the LC starts forming. The low-${ \mathcal F }$ magnetospheres (i.e., top row) tend to be charge-separated (similar to the dome-torus configuration; see Krause-Polstorff & Michel 1985; Spitkovsky & Arons 2002), while the high-${ \mathcal F }$ values produce magnetospheres with both kinds of charge carriers being abundant everywhere (i.e., bottom row).

Zoom In Zoom Out Reset image size

Figure 5 shows, in logarithmic scale on the ${\boldsymbol{\mu }}\mbox{--}{\boldsymbol{\Omega }}$ plane, the number density of the injected particles per unit time for the indicated α values and for a high injection rate (${ \mathcal F }=20{{ \mathcal F }}_{\mathrm{GJ}}^{0}$). The much higher densities near the stellar surface are compensated for by the much larger volumes at larger distances. We note that our injection prescription activates particle injection only when the magnetization σ_M exceeds some threshold value (Equation (4)), and so particle injection is strongly suppressed in the ECS region outside the LC where $B\to 0$. This implies that in our simulations, the particles appearing eventually in the ECS originate from the separatrix within the LC and/or reach the ECS region by drifting across magnetic field lines (see also the discussion in Section 6). In Brambilla et al. (2018), studying magnetospheres with particles injected only near the stellar surface, we saw drifting positrons contributing to the structure of the ECS.

Zoom In Zoom Out Reset image size

Figure 6 shows for the indicated α and ${ \mathcal F }$ values the normalized particle energy (i.e., γ_L) distributions for both e⁻ (blue lines) and e⁺ (red lines). Even though the particle energy range, which in reality extends to many orders of magnitude, is compressed within three orders of magnitude, our simulations are still able to provide not only the current and field structures corresponding to an entire spectrum of solutions from vacuum to FF but also the accelerating component of the particle distribution and the magnetosphere regions where the particle acceleration takes place. Thus, the particle energy distributions consist of a low-energy population that peaks at γ_L ∼ 1 and one at higher energy (γ_L ≳ 30) that extends to γ_L ∼ 10³. The first component is the non-accelerated bulk plasma, which is associated mainly with the short-scale (∼λ_D) fluctuating fields while the second component is the accelerating component, which is associated with the larger-scale unscreened fields of the magnetosphere. We note that the relative strength of the accelerating component decreases with increasing ${ \mathcal F }$. As we also discuss later, a significant part of the particle acceleration takes place at large distances at and beyond the LC, even for relatively low ${ \mathcal F }$ (i.e., ${ \mathcal F }\gtrsim 1{{ \mathcal F }}_{\mathrm{GJ}}^{0}$). However, for high-${ \mathcal F }$ values and solutions close to the FF ones, the acceleration takes place exclusively at and beyond the LC in regions near the ECS.

Zoom In Zoom Out Reset image size

In Figure 7 we plot, for different α values and for ${ \mathcal F }=28{{ \mathcal F }}_{\mathrm{GJ}}^{0}$, the average (within a computational cell) γ_L value. These simulations have field structures very close to the FF ones, and the average γ_L values indicate the existence of energetic particles only in the region close to the Y point⁷ and the ECS beyond the LC. The FF ECS for α = 0° is positively charged (i.e., ρ > 0), but as α increases, the corresponding charge density decreases, and it becomes 0 for α = 90°.⁸ This can be concluded also from the fact that the orthogonal rotator α = 90° actually lies in the middle of the aligned and anti-aligned rotator. In Figure 6, the consistent excess of positrons at high energies reflects the strong acceleration that takes place at and beyond the LC near the ECS.⁹ The effect is suppressed as we go toward higher α values.

Zoom In Zoom Out Reset image size

In the previous section, we saw that an entire spectrum of solutions from near VRD to near FF can be covered by applying arbitrary particle injection everywhere in the magnetosphere. However, a vital question that arises is whether and how these simulations can be compared to recent rich observational data (i.e., Fermi). The comparison with the observations requires a detailed modeling of the high-energy emission in PIC simulations. This implies not only the determination of the emitting regions but also the derivation of the corresponding ${\epsilon }_{\mathrm{cut}}$ and L_γ values. However, the PIC particle energies are much smaller than the energies in a real pulsar magnetosphere (see Figure 6) due to the much smaller fields (see Section 3). In addition, there is no trivial relation (e.g., proportionality) that can provide the particle energy transformation between the PIC model and the real system.

The particle energies in the real system are not only larger (than those in the PIC simulations)—they are also crucially affected by the radiation reaction forces, with the phenomenon being more prominent for high-energy emitting particles. In our PIC simulations, the impact of the actual radiation reaction forces (Equation (2)) on the particle energies is negligible (due to the correspondingly small γ_L and B values). Moreover, in principle, synchrotron and curvature losses act very differently, depending strongly on the specific regime (i.e., γ_L, B, and radius of curvature R_C of the guiding center). Furthermore, the gyroradius in PIC simulations is much larger than that in real pulsar magnetospheres (in R_LC units), and in low-field regions this can, in principle, affect magnetosphere features (see Cerutti et al. 2015).

In principle, the radiation reaction force, ${F}_{\mathrm{rr}}$, on the PIC particles can be enhanced (artificially) to the level that produces cooling times equal (in pulsar period units) to those of realistic B_⋆ and P magnetosphere values, to establish the correct loss timescales. However, the realistic cooling times are so small, especially for the case of high-$\dot{{ \mathcal E }}$ values, that the time steps needed to capture them are far too short to be implemented in our simulations. In addition, the maximum particle energies would be severely reduced relative to the maximum potential drop (because of the drastic enhancement in F_rr); this effect would be more dramatic for the high-$\dot{{ \mathcal E }}$ cases. This reduction would then bring the entire particle population (including the accelerating component) either close to or within the statistical noise level that inherently exists in these simulations, thereby invalidating any claim of realistic radiation modeling.

In real pulsar magnetospheres, particles will lose their perpendicular (to the magnetic field) momentum, and therefore their pitch angle, very quickly due to the corresponding synchrotron losses (unless there is some mechanism to sustain pitch angles). Under this regime, particles remain at the zero Landau level, thus terminating the synchrotron radiation emission. Nonetheless, there are, in principle, certain conditions that can sustain the pitch angles and make synchrotron emission important. Lyubarskii & Petrova (1998) showed that cyclotron resonant absorption of low-energy radio photons is efficient in pulsar magnetosphere environments. In regions (at high altitudes within the LC) where the resonant condition is fulfilled, the Landau states are populated, and the pitch angles of the particles increase up to the equilibrium point where the absorption energy gain balances the energy losses due to (mainly) synchrotron emission. In Harding & Kalapotharakos (2015; see also Harding et al. 2008), we showed that this process may be important in some cases (e.g., in the Crab pulsar) to explain the optical to hard X-ray observed emission. Moreover, Uzdensky & Spitkovsky (2014) showed that in the reconnection region at the ECS of a Crab-like pulsar, equipartition can trigger synchrotron emission. Nonetheless, any synchrotron emission that is triggered by cyclotron resonant absorption of low-energy radio photons contributes to energies much smaller than the cutoff energies observed by Fermi while the microphysical emission properties of the reconnection at the ECS is difficult to study properly in global magnetosphere simulations that resolve this region rather poorly.

Taking into account all of the above, we focused on the properties and efficiency of CR assuming that the contribution of the synchrotron emission is small (especially in the Fermi energy band) because of the rapid decrease of the perpendicular (to B) momentum.

However, the CR calculation depends on both the electric fields the particles encounter and the orbital radius of curvature. On the one hand, the electric fields and the corresponding potential drops depend exclusively on the magnetosphere model (i.e., on the global particle injection rate). On the other hand, an accurate calculation of CR implies that the orbital radii of curvature are close to those corresponding to asymptotic drift approximation trajectories. The latter ones are actually defined to be those of the so-called Aristotelian electrodynamics (AE; see Gruzinov 2012, 2013; Kelner et al. 2015).

For these asymptotic trajectories, the velocity flow reads

$$\begin{eqnarray}&&{\boldsymbol{v}}=\displaystyle \frac{{\boldsymbol{E}}\times {\boldsymbol{B}}\pm ({B}_{0}{\boldsymbol{B}}+{E}_{0}{\boldsymbol{E}})}{{B}^{2}+{E}_{0}^{2}},\end{eqnarray} \tag{ 5 }$$

where the two signs correspond to the two different types of charges, while the quantities E₀ and B₀ are related to the Lorentz invariants (Gruzinov 2008; Li et al. 2012)

$$\begin{eqnarray}&&{E}_{0}{B}_{0}={\boldsymbol{E}}\cdot {\boldsymbol{B}},\,\,{E}_{0}^{2}-{B}_{0}^{2}={E}^{2}-{B}^{2},\end{eqnarray} \tag{ 6 }$$

and E₀ is the electric field in the frame where E and B are parallel and is the effective accelerating electric component ${E}_{\mathrm{acc}}$, which becomes zero only when E · B = 0 and E < B. We note that particles that follow trajectories defined by Equation (5) do not emit synchrotron radiation (Kelner et al. 2015). For the rest of the paper, E₀ denotes ${E}_{\mathrm{acc}}$.

In our simulations, the high-energy particles, in general, do not have high pitch angles because on the one hand the particles are injected with zero velocities, and on the other hand, the acceleration along the accelerating electric fields tends to reduce any pitch angle. Nonetheless, in order to make the CR calculation more accurate, we used an artificially enhanced F_rr in order to reduce the pitch angles even further. For this, we increased properly the magnetic fields B_eff that enter the radiation reaction forces on the PIC particles (Equation (2)). The B_eff scaling varies in time (for each particle trajectory), ensuring always that the corresponding synchrotron cooling time

$$\begin{eqnarray}&&{t}_{\mathrm{sc}}=-\displaystyle \frac{\gamma }{\dot{\gamma }}\propto {B}_{\mathrm{eff}}^{-2}{\gamma }_{{\rm{L}}}^{-1}\end{eqnarray} \tag{ 7 }$$

is effectively small, which means it is much smaller than the light crossing time t_l = 1/Ω but always resolved by the adopted time step (at least five times). We note that the B_eff fields affect only the radiation reaction forces (i.e., Equation (2)) and not the ordinary Lorentz force, ${{\boldsymbol{F}}}_{{\rm{L}}}={q}_{{\rm{e}}}({\boldsymbol{E}}+{\boldsymbol{v}}\times {\boldsymbol{B}}/c),$ where the actual (i.e., simulated) fields enter. The resulting trajectories (especially the high-energy ones) have geometric properties closer to the real ones (i.e., the perpendicular momentum is significantly reduced). We have actually checked that the vast majority of the PIC high-energy particles have velocities and orbital radii of curvature close to those corresponding to AE (i.e., Equation (5)).

Nonetheless, the resulting particle energies are still much smaller than the particle energies in a real pulsar magnetosphere. In order to derive PIC particle energies with meaningful values, we scale up the particle energies to the values they would have for real pulsar surface fields assuming that the trajectories of the high-energy particle are geometrically correct. For the calculation of the scaled-up (i.e., realistic) particle energies, γ_R, done in parallel inside the PIC simulation, we choose realistic (P, B_⋆) values, and we integrate along each particle trajectory the expression

$$\begin{eqnarray}&&\displaystyle \frac{d{\gamma }_{{\rm{R}}}}{{dt}}=\displaystyle \frac{{q}_{{\rm{e}}}{\boldsymbol{v}}\cdot {\boldsymbol{E}}}{{m}_{{\rm{e}}}{c}^{2}}-\displaystyle \frac{2{q}_{{\rm{e}}}^{2}{\gamma }_{{\rm{R}}}^{4}}{3{R}_{{\rm{C}}}^{2}{m}_{{\rm{e}}}c},\end{eqnarray} \tag{ 8 }$$

which provides realistic γ_R values. We note that in Equation (8) R_C is the local radius of curvature of the trajectory while v is the corresponding particle velocity. The integration of Equation (8) is made in parallel with the PIC particle trajectories, and it does not affect the trajectories themselves. It is used in order to provide realistic estimations (i.e., in real pulsar environments) of particle energies, taking into account the realistic energy gains (i.e., the first term in Equation (8)) due to any accelerating electric fields and the energy losses due to CR reaction losses (i.e., the second term in Equation (8)). For the calculation of R_C, three trajectory points that are separated in time by d/c, where d is the computational cell size, are used. This implies that the R_C value is updated at intervals of d/c.¹⁰ Apparently, as discussed above, R_C is not critically affected by the (reduced) gyromotion due to the enhanced ${{\boldsymbol{F}}}_{\mathrm{rr}}$. Thus, in our adopted approach, instead of using test particles in order to estimate the effectiveness of CR, we use the particles that the kinetic simulations indicate lie in the model's strong accelerating regions. Finally, for the integration of Equation (8), we take into account only the effective accelerating electric field components, as described below.

In our PIC simulations, there are magnetosphere regions where nearly all of the corresponding accelerating electric components are screened, and magnetosphere regions with consistent accelerating electric components. The former regions are the effective FF regions where the field fluctuations dominate while the latter regions are the actual acceleration locations. We have found that the best way to identify these regions is to use the local values of the normalized effective accelerating electric field E₀/E, where E is the local total electric field (in the observer's inertial frame). In Figure 8, we plot, in $\mathrm{log}-\mathrm{log}$ scale, the distributions of E₀/E values that have been calculated at points that are uniformly distributed within the magnetosphere volume that extends from the stellar surface up to r = 2.5 R_LC. These results are for simulations with α = 45° corresponding to the different particle injection rates (${ \mathcal F };$ see Table 1) that are denoted by the different colors according to the indicated color scale. We observe the presence of a maximum at $\mathrm{log}({E}_{0}/E)\lesssim -2$ associated with a thermal-like particle distribution along with a bump that always appears in these distributions at higher values. The thermal-like maximum is related to the fluctuating (noisy) field behavior; its value depends primarily on the time step and the number of macroparticles in our simulations. The higher-value component (i.e., the bump), beyond the vertical dashed line, actually indicates the magnetosphere regions where truly high acceleration takes place. We have checked that the E₀/E position of the maximum decreases (increases) as we reduce (enhance) noise. Thus, when we increase the number of macroparticles, and/or decrease the time step, and/or increase the spatial resolution, and/or increase the effective area of the current smoothing kernel, the maximum of the distributions shown in Figure 8 moves to lower $\mathrm{log}({E}_{0}/E)$ values. In Figure 8, the gradual drift of the maximum toward higher $\mathrm{log}({E}_{0}/E)$ as ${ \mathcal F }$ increases is due to the corresponding poorer resolution for λ_D and ω_p (a.k.a. higher noise level). The increase of ${ \mathcal F }$ implies higher particle number densities and consequently smaller λ_D and higher ω_p. Keeping all other parameters (e.g., d, dt) the same, these lower λ_D and higher ω_p values are more poorly resolved. This gradually enhances the noise, leading to the increase of the position of the maximum. On the other hand, we see that the area under the nonthermal extension (i.e., bump) decreases with increasing ${ \mathcal F }$. This implies smaller magnetosphere volumes with consistent accelerating electric field components and lower actual accelerating electric fields. Actually, Figure 8 also indicates our limitations for the highest ${ \mathcal F }$ values we can use for a specific parameter set. Thus, as ${ \mathcal F }$ increases and the $\mathrm{log}({E}_{0}/E)$ maximum moves toward higher values, the gradually reducing area under the bump is eventually (at some ${ \mathcal F }$ value) buried under the corresponding noise.

Zoom In Zoom Out Reset image size

In Figure 9, we show, in the indicated linear color scale, the $\mathrm{log}({E}_{0}/E)$ values on the poloidal ${\boldsymbol{\mu }}\mbox{--}{\boldsymbol{\Omega }}$ plane for the α = 45° simulations and for the indicated ${ \mathcal F }$ values. The reddish regions denote the regions with values higher than the value corresponding to the dashed vertical line shown in Figure 8. We see that at lower ${ \mathcal F }$ values, large magnetosphere regions have consistent accelerating electric fields, while for higher ${ \mathcal F }$ values, the reddish regions shrink and finally are restricted to only near the ECS outside the LC. For all of the model results presented below, we assume that the particle acceleration that produces the high-energy emission takes place in the magnetosphere regions corresponding to values greater than the one indicated by the dashed line in Figure 8. We have decided to use the same threshold value for consistency even though the bump appears at slightly different values. However, we have checked that taking into account the regions corresponding to smaller $\mathrm{log}({E}_{0}/E)$ values (for the simulations where the bump appears earlier than the dashed line in Figure 8) does not change the results, mainly because the corresponding accelerating electric components are small. Moreover, we have seen that the Ohmic J · E dissipation takes place largely in the magnetosphere regions corresponding to $\mathrm{log}({E}_{0}/E)$ values higher than the dashed line threshold. For the rest of the magnetosphere, the total J ·  E dissipation tends to zero.

Zoom In Zoom Out Reset image size

Finally, we note that in real pulsar magnetospheres, the noise level (i.e., the maximum of the distribution in Figure 8) is expected to be much smaller relative to the maximum of E₀/E in our simulations. This allows, in principle, the existence of consistent accelerating electric fields for $\mathrm{log}({E}_{0}/E)$ values much lower than those we assume in our study. Nonetheless, these accelerating electric fields are expected to be important for the lower part of the emission spectrum (and not for energies near the ${\epsilon }_{\mathrm{cut}}$ values observed by Fermi).

In real pulsars, the direct measurements of the period P and its derivative $\dot{P}$ allow the spin-down power to be calculated. Assuming that the moment of inertia I is known, the spin-down power is given by

$$\begin{eqnarray}&&\dot{{ \mathcal E }}=\displaystyle \frac{4{\pi }^{2}I\dot{P}}{{P}^{3}},\end{eqnarray} \tag{ 9 }$$

which can indirectly lead to an estimate of the magnetic field on the stellar surface. However, the spin-down power depends on the magnetosphere regime and is given by (Deutsch 1955; Spitkovsky 2006)

$$\begin{eqnarray}&&\dot{{ \mathcal E }}=\displaystyle \frac{8{\pi }^{4}{r}_{\star }^{6}}{3{c}^{3}}\displaystyle \frac{{B}_{\star }^{2}}{{P}^{4}}{\sin }^{2}\alpha \,\,\,\mathrm{for}\ \mathrm{the}\ \mathrm{VRD}\ \mathrm{regime},\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&\dot{{ \mathcal E }}=\displaystyle \frac{4{\pi }^{4}{r}_{\star }^{6}}{{c}^{3}}\displaystyle \frac{{B}_{\star }^{2}}{{P}^{4}}(1+{\sin }^{2}\alpha )\,\,\,\mathrm{for}\ \mathrm{the}\ \mathrm{FF}\ \mathrm{regime}.\end{eqnarray} \tag{ 10b }$$

By combining Equations (9) and (10), we can obtain the B_⋆ value with an uncertainty due to the unknown regime and the unknown α value. In the left-hand (right-hand) panel of Figure 10, we plot on (P, B_⋆) diagrams all of the Fermi YPs (MPs). Each pulsar is represented by a horizontal light gray line that denotes the uncertainty in the B_⋆ estimation. We note that the left and right ends of the light gray lines denote the B_⋆ values that correspond to α = 90° at the FF and VRD regimes, respectively. The color stripes denote the indicated $\dot{{ \mathcal E }}$ values using the same B_⋆ uncertainty as the light gray lines for the Fermi pulsars.

Zoom In Zoom Out Reset image size

For the study of the high-energy emission in our simulations, we follow the approach that is discussed in the previous section. Table 2 shows the realistic (P, B_⋆) values that we have chosen. The calculations described in the previous section are computationally demanding, and therefore we have selected only 12 cases (six for YPs and six for MPs) that trace well the corresponding observed area (see the black points in Figure 10).

Table 2.  The 12 (6+6) Adopted Realistic (B_⋆, P) Value Sets for the YP and MP Models. For Each Value Set, the FF Spin-down Power Range Is Indicated (see Equation 10(b))

Model	Young Pulsars			Millisecond Pulsars
No.	B⋆	P	${\dot{{ \mathcal E }}}_{\mathrm{FF}}$	B⋆	P	${\dot{{ \mathcal E }}}_{\mathrm{FF}}$
	(1012 G)	(ms)	(erg s−1)	(108 G)	(ms)	(erg s−1)
1	0.63	302.0	(0.5–1.4) × 1033	0.7	5.1	(0.7–2.0) × 1032
2	1.26	239.9	(0.5–1.4) × 1034	2.2	5.0	(0.7–2.2) × 1033
3	2.24	177.8	(0.5–1.4) × 1035	5.0	4.1	(0.9–2.6) × 1034
4	3.16	120.2	(0.5–1.4) × 1036	3.3	2.7	(2.0–6.0) × 1034
5	3.98	75.9	(0.5–1.4) × 1037	3.5	2.0	(0.7–2.2) × 1035
6	7.94	60.3	(0.5–1.4) × 1038	10.0	1.8	(0.9–2.6) × 1036

Download table as:  ASCIITypeset image

By calculating the realistic particle ${\gamma }_{{\rm{R}}}$ values, we can calculate not only the individual particle emissivity ($\propto {\gamma }_{{\rm{R}}}^{4}{R}_{{\rm{C}}}^{-2};$ see Equation (8)) but also the entire CR spectrum of the corresponding emission. In Figure 11, we plot, in $\mathrm{log}\mbox{--}\mathrm{log}$ scale, the normalized distribution of the realistic γ_R values corresponding to a simulation with α = 45° and ${ \mathcal F }=28{F}_{\mathrm{GJ}}^{0}$. The top and bottom panels show the distributions for the different (P, B_⋆) values that correspond to YPs and MPs, respectively.

Zoom In Zoom Out Reset image size

Taking into account the emission from the entire magnetosphere (up to 2.5 R_LC), we construct the "global" spectrum, which we then fit with the model used in 2PC, namely,

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\epsilon }=A{\epsilon }^{-{\rm{\Gamma }}}\exp \left(-\displaystyle \frac{\epsilon }{{\epsilon }_{\mathrm{cut}}}\right),\end{eqnarray} \tag{ 11 }$$

where Γ is the photon index.

Figure 12 shows, in $\mathrm{log}\mbox{--}\mathrm{log}$ scale, the model ${\epsilon }_{\mathrm{cut}}$ values as a function of the particle injection rates ${ \mathcal F }$. The three different columns show the results for the indicated α values, while the top and bottom rows show the results for the YP and MP model cases, respectively. The light greenish zones in each panel of Figure 12 denote the zone of the observed (by Fermi) ${\epsilon }_{\mathrm{cut}}$ values (∼1–6 GeV).¹¹ The different colors indicate the results for the different (P, B_⋆) (i.e., $\dot{{ \mathcal E }}$) values that are shown in Table 2. First, we see that, for the same ${ \mathcal F }$ value, the higher the $\dot{{ \mathcal E }},$ the higher the corresponding ${\epsilon }_{\mathrm{cut}}$ value. This increase is due to the corresponding higher accelerating electric components E₀. We see also that for the same $\dot{{ \mathcal E }}$ value (i.e., along a line) and for high ${ \mathcal F }$ values $(\gtrsim 1{{ \mathcal F }}_{\mathrm{GJ}}^{0})$, ${\epsilon }_{\mathrm{cut}}$ decreases with ${ \mathcal F }$. This decrease is more prominent for low α and $\dot{{ \mathcal E }}$ values. On the other hand, toward low ${ \mathcal F }$ values, ${\epsilon }_{\mathrm{cut}}$ either decreases (for low α) or tends to stabilize (for high α). We have already discussed that low ${ \mathcal F }$ values produce solutions near VRD whose spin-down power (i.e., the Poynting flux) for low α values goes to very low values (see Equation 10(a)). The low $\dot{{ \mathcal E }}$ values imply a reservoir of low E fields, and therefore low accelerating electric components E₀ (≤E). This is mainly the reason why, for low α values, ${\epsilon }_{\mathrm{cut}}$ values decrease toward low ${ \mathcal F }$.

Zoom In Zoom Out Reset image size

Figure 12 provides a possible explanation for why we start observing YPs and MPs for spin-down powers ≳10³⁴ erg s⁻¹ and ≳10³³ erg s⁻¹, respectively. We see that for low ${ \mathcal F }$ values, the low $\dot{{ \mathcal E }}$ YP and MP models (see Table 2) struggle to reach the zone observed by Fermi (≈1 GeV). Taking also into account now the Fermi threshold (≈100 MeV) and the sensitivity of the instrument, we understand that YPs and MPs with $\dot{{ \mathcal E }}$ values lower than those observed are difficult to detect with Fermi. These objects, if they exist, should have spectra with ${\epsilon }_{\mathrm{cut}}$ values lower than those observed by Fermi (at the MeV levels). We note here that all MeV pulsars detected thus far (Kuiper & Hermsen 2015) have high $\dot{{ \mathcal E }}$ values. However, the light-curve characteristics of these objects indicate that the corresponding line of sight (i.e., ζ value) might not cross the core of the high-energy emission, particularly for small α and for intermediate ζ values (see Figure 15; see also Harding & Kalapotharakos 2017).

On the other hand, models with higher $\dot{{ \mathcal E }}$ values lie within the Fermi zone for a range of ${ \mathcal F }$ values adopted in the current study. We see also that the high $\dot{{ \mathcal E }}$ models (≳10³⁷ erg s⁻¹ for YPs and ≳10³⁵ erg s⁻¹ for MPs) require higher ${ \mathcal F }$ values than those we have applied. However, obtaining higher ${ \mathcal F }$ values is a computationally very demanding task, and so we get an estimation of the corresponding ${ \mathcal F }$ values through linear extrapolations (see dashed lines in Figure 12).

In Figure 13, we plot, in $\mathrm{log}\mbox{--}\mathrm{log}$ scale, the ${ \mathcal F }$ value as a function of $\dot{{ \mathcal E }}$ that produces ${\epsilon }_{\mathrm{cut}}$ values within the Fermi zone. The left-hand (right-hand) panel shows the results for the YP (MP) models. The horizontal dimension of each rectangle denotes the model $\dot{{ \mathcal E }}$ range (see Table 2) while the vertical dimension denotes the model ${ \mathcal F }$ range (for all α values) that provides ${\epsilon }_{\mathrm{cut}}$ values within the Fermi zone. We note that for the optimum ${ \mathcal F }$ range, we have taken into account the extrapolations shown in Figure 12. The empty rectangles denote that no ${ \mathcal F }$ value can take the models within the Fermi zone. In this case, the vertical ${ \mathcal F }$ range denotes the values that produce ${\epsilon }_{\mathrm{cut}}$ values closest to the Fermi zone. We see that the highest ${ \mathcal F }$ for the high $\dot{{ \mathcal E }}$ models can reach up to ${10}^{5}\mbox{--}{10}^{6}{{ \mathcal F }}_{\mathrm{GJ}}^{0}$, which are consistent with the particle multiplicities found in local pair cascade studies (Timokhin & Harding 2015) and to the multiplicities that are needed to explain the Crab Nebula spectrum (de Jager et al. 1996).

Zoom In Zoom Out Reset image size

In Figure 14, we plot the ${\epsilon }_{\mathrm{cut}}$ (left-hand column) and L_γ (right-hand column) values for all YP (top row) and MP (bottom row) models together with the Fermi data, as functions of the corresponding spin-down power, $\dot{{ \mathcal E }}$. We note that the darker the model points, the higher the corresponding ${ \mathcal F }$ value. Similarly to Figure 12, Figure 14 reveals (from a different perspective than Figure 12) the model ${ \mathcal F }$ value ranges that are consistent with both ${\epsilon }_{\mathrm{cut}}$ and L_γ Fermi values. In the right-hand column, we see that a maximum appears as ${ \mathcal F }$ varies from low to high values. This expected behavior reflects the ideal nature of the VRD and FF solutions. The former solutions have high fields but limited numbers of particles to dissipate the energy while the latter ones have high numbers of particles but limited energy available for dissipation. Taking into account the corresponding $\dot{{ \mathcal E }}$ value, we see that the maximum efficiency ranges from ≈13% to ≈35% and is higher (lower) for the lower (higher) α values. We note that Figure 14 implies that the optimum ${ \mathcal F }$ values that reproduce the 2PC L_γ values are somehow higher then the optimum ${ \mathcal F }$ values that reproduce the 2PC ${\epsilon }_{\mathrm{cut}}$ values. Nonetheless, as we have discussed in Kalapotharakos et al. (2017), the 2PC L_γ values are less reliable than the corresponding ${\epsilon }_{\mathrm{cut}}$ values. The former ones depend on the assumed beaming-factor F_b and the estimated distances, while their wide spread with $\dot{{ \mathcal E }}$ indicates that other factors (i.e., α values, variability of F_b with observer angle, ζ) play an important role in their determination. On the other hand, the range of ${\epsilon }_{\mathrm{cut}}$ values is more limited, it does not suffer from geometry or distance uncertainties, and depends weakly only on the adopted fit model. A detailed study of the model F_b values that takes into account the observationally dependent parameters (e.g., α, ζ values, distances, instrument thresholds etc.) could eventually answer whether the observed γ-ray fluxes are consistent either with the model L_γ values corresponding to the optimum ${ \mathcal F }$ values that reproduce the observed ${\epsilon }_{\mathrm{cut}}$ values or to the 2PC L_γ values, which always assume F_b = 1. Even though this study goes beyond the scope of the present paper, our preliminary calculations indicate that the average model F_b values are smaller than 1, which means, assuming model correctness, that the 2PC L_γ values are overestimated. Finally, we note that the ${\epsilon }_{\mathrm{cut}}$ of MPs for the same $\dot{{ \mathcal E }}$ and ${ \mathcal F }$ are ∼3× higher than the corresponding ones of YPs, consistent with Fermi data (Kalapotharakos et al. 2017). This results from the ${\epsilon }_{\mathrm{cut}}\propto {B}_{\star }^{-1/8}$ relation that is presented in Kalapotharakos et al. (2017) and comes from the assumption that emission takes place at the LC near the ECS due to CR at the radiation-reaction-limit regime. Taking into account that ${B}_{{\star }_{\mathrm{MP}}}\approx {10}^{-4}{B}_{{\star }_{\mathrm{YP}}}$, the previous relation leads to ${\epsilon }_{{\mathrm{cut}}_{\mathrm{MP}}}\approx 3{\epsilon }_{{\mathrm{cut}}_{\mathrm{YP}}}$.

Zoom In Zoom Out Reset image size

In Figures 15–17, we plot sky-map atlases for α = 15°, α = 45°, and α = 75°, respectively. In each figure, the different rows correspond to the different indicated particle injection rates (${ \mathcal F }$) while the different columns correspond to the different indicated YP models (i.e., different P, B_⋆ values and therefore different $\dot{{ \mathcal E }}$ values; see Table 2). Each sky map shows, in the indicated color scale, the emitted luminosity per solid angle (${dL}/{dw}={dL}/\sin (\zeta )d\zeta d{\phi }_{\mathrm{ph}}$). The horizontal axis depicts the phase, ϕ_ph, of the pulsar rotation, while the vertical axis depicts the observer inclination angle, ζ, which is the angle between the rotational axis and the line of sight. Horizontal cuts of each sky map correspond to the light curve that is observed by an observer that lies at some ζ value. We note that ϕ_ph = 0 corresponds to the observed phase of a photon that originates from the magnetic pole at $r\to 0$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The sky maps were produced by collecting all of the photons emitted by the particles from t = 1.5P to t = 2.0P (half a period). The simulations are in a steady-state equilibrium during that time interval, and therefore the photon collection from multiple time steps (during this interval) reduces the noise considerably. The emissivity of each particle is considered to be $\propto {\gamma }_{{\rm{R}}}^{4}{R}_{{\rm{C}}}^{-2}$, and the corresponding photon emission is along the particles' velocities. To determine the photons' relative time of arrival (i.e., ϕ_ph), we take into account the corresponding relativistic time delay effects. Thus, ϕ_ph reads (see also Kalapotharakos et al. 2014)

$$\begin{eqnarray}&&{\phi }_{\mathrm{ph}}=\left({\rm{\Omega }}{t}_{{\rm{S}}}-{\phi }_{{{\boldsymbol{v}}}_{{\rm{p}}}}-\displaystyle \frac{{{\boldsymbol{r}}}_{{\rm{p}}}\cdot {{\boldsymbol{v}}}_{{\rm{p}}}}{{v}_{{\rm{p}}}}\displaystyle \frac{1}{{R}_{\mathrm{LC}}}\right)\mathrm{mod}2\pi ,\end{eqnarray} \tag{ 12 }$$

where t_S is the global simulation time, v_p and r_p are the particle velocity and position vectors, and ${\phi }_{{{\boldsymbol{v}}}_{{\rm{p}}}}$ is the azimuth angle of the particle velocity v_p with respect to the magnetic axis at t_S = 0 oriented according to ${\boldsymbol{\Omega }}$. The last term in Equation (12) formulates the light travel time delay.

The sky-map patterns (especially toward high ${ \mathcal F }$ values) are actually similar to those of the FIDO model (see Figure 17 of Kalapotharakos et al. 2014). For large inclination angles, the locus of the ECS deviates considerably from the rotational equator, allowing, in principle, the emission to contribute to very low and very high ζ values. However, in our models (of high α values), similarly to the FIDO models, the main part of the emission does not come from the high latitude (with respect to the rotational equator) zones. This makes the corresponding emission more concentrated around the rotational equator (i.e., around ζ = 90°).

Actually, in Figure 18, we plot, for the indicated α values and for ${ \mathcal F }=20{{ \mathcal F }}_{\mathrm{GJ}}^{0}$, the particles that produce the highest 95% of the corresponding total L_γ value. The individual particle colors represent the corresponding emissivity as indicated in the color scale. The particle location is always near the ECS. However, we see that for low α values, the particle distribution is quite uniform, while for high α values, most of the emission is produced by particles that are concentrated mainly closer (with respect to the theoretical extent) to the rotational equator. This behavior is similar to what had been originally seen in the FIDO models (see Figure 19(b) of Kalapotharakos et al. 2014). A similar behavior is also seen in Cerutti et al. (2016). This implies that the E_acc, far from the equator, are shorted out relatively easily (i.e., the corresponding current requirement is not that strong). Figure 18 indicates also that for low α values, the emission is more concentrated near the LC while for high α values, it extends to larger distances (see the distribution of the red points in Figure 18).

Zoom In Zoom Out Reset image size

Moreover, Figures 15–17 indicate that toward high ${ \mathcal F }$ and $\dot{{ \mathcal E }}$ values, the luminous sky-map regions become narrower. This effect is further enhanced if we take into account that ${ \mathcal F }$ and $\dot{{ \mathcal E }}$ increase together (see Figure 13). This behavior is consistent with the observations, which show broader γ-ray light curves toward low $\dot{{ \mathcal E }}$ values of YPs and MPs (2PC).

In Figures 19–21, we present, as a demonstration, γ-ray light curves for the indicated ${ \mathcal F }$, α, $\dot{{ \mathcal E }}$, and ζ values. These correspond to the ζ values indicated by the dashed white horizontal lines shown in Figures 15–17. The γ-ray light curves corresponding to high ${ \mathcal F }$ values appear very similar to the observed ones (2PC). As mentioned above, in these cases, the sky maps and the corresponding γ-ray light curves are similar to the ones from the FIDO model, and so they yield radio lag values, δ, consistent with those indicated by Fermi (assuming that the radio emission comes from near the polar cap region).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, we note that there are cases (mainly for low ${ \mathcal F }$ and high α values; see the first rows of Figures 17 and 21) when the light curves seem to have more complicated features (e.g., more than two peaks, considerable interpeak and off-peak emission) not seen in the majority of the observed γ-ray light curves. This is mainly due to the adopted "uniform"¹² regulation of the particle injection rate. By reducing ${ \mathcal F }$ uniformly, we start introducing accelerating electric fields in regions (even inside the closed zone) that produce additional emission components that are not always consistent with the observed features. However, we note that peculiar and unique features in γ-ray light curves exist in some MPs that have low $\dot{{ \mathcal E }}$ values and are expected to be pair starved (i.e., low ${ \mathcal F }$ values; see Muslimov & Harding 2004; Johnson et al. 2014). Independently of this, it seems that the construction of realistic low ${ \mathcal F }$ models (i.e., weak pulsars, Gruzinov 2013; Contopoulos 2016) should be made more carefully (compared to the approach followed in this study). Our results indicate that the observed γ-ray emission is regulated mainly by the particle abundance in regions near the ECS. For the rest of the magnetosphere regions, the accelerating electric components are more easily screened, and beyond some point (i.e., some ${ \mathcal F }$ value), no considerable emission is produced. Of course, the role of these regions might be important for the lower-energy part of the spectrum (see also the related discussion in Section 5).

In Figure 22, we plot, on the poloidal ${\boldsymbol{\mu }}\,\mbox{--}\,{\boldsymbol{\Omega }}$ plane, the origin of the particles (red for e⁺ and blue for e⁻) that eventually produce the highest 95% of the observed YP emission, in the case of ${ \mathcal F }=28{{ \mathcal F }}_{\mathrm{GJ}}^{0}$. We see that the particles that produce most of the observed emission are e⁺ and originate along the separatrix. This particle population goes into the ECS, where it emits. A lower-energy population that originates mainly from the same locus self-consistently regulates the corresponding accelerating electric fields.¹³ From all of the above, it becomes apparent that an approach similar to what we followed in Kalapotharakos et al. (2017) becomes essential. In Kalapotharakos et al. (2017; see also Contopoulos et al. 2014), we assumed that the dissipative region is only near the ECS beyond the LC, and we applied the finite σ values only there. The equivalent in the PIC simulations would be to apply a rather high and constant particle injection rate in all magnetosphere regions except for the ones indicated in Figure 22. It seems that the regulation of only this population would provide results that are in full agreement with the observations for the entire range of the parameter space. This approach highlights the role of the return-current region in the high-energy emission observed in pulsars, something that, eventually, is expected to be supported by the (micro)physically motivated pair creation.

Zoom In Zoom Out Reset image size

In this paper, we present 3D kinetic global models of pulsar magnetospheres. In our study, we use C-3PA, an efficient, 3D Cartesian relativistic PIC code we have developed. In this first step, we ignore the microphysics of particle production, and so we apply ad hoc particle injection rates (${ \mathcal F }$) that at least produce pulsar magnetosphere models where the field structure and the corresponding particle distribution are consistent with each other.

The particle injection is applied up to large distances (2.5R_LC) and is regulated by the local magnetization. From low to high particle injection rates, we cover an entire spectrum of solutions from near VRD to near FF, respectively. For models with high particle injection rates ($\gt 1{{ \mathcal F }}_{\mathrm{GJ}}^{0}$), particle acceleration takes place mainly near the region of the ECS outside the LC, while the volume of this region decreases with the particle injection rate.

A main goal of this study is to uncover the properties of the model γ-ray emission and how these vary with the particle injection rate. However, in our PIC simulations, the particle Lorentz factors are (due to numerical limitations) much smaller than those in real pulsar magnetospheres. Moreover, the relative behaviors of synchrotron and curvature losses can be quite different for different particle energy regimes. In our approach, we assume that the synchrotron losses in the environment of real pulsar magnetospheres rapidly reduce the corresponding pitch angles, and so the synchrotron emission is not the main component of the observed γ-ray component. In our simulations, we use artificially scaled-up magnetic fields inside the radiation reaction force expressions in order to reduce the synchrotron cooling times. This mostly affects the geometry of the trajectories because of the suppressed gyromotion. At the same time, along the PIC particle trajectories, we derive realistic particle Lorentz factors by rescaling the fields and the length and timescales, always taking into account the corresponding energy gain and CR reaction losses. This approach allows the corresponding high-energy emission to be calculated (i.e., sky maps, light curves, spectra).

Using realistic value sets for the stellar surface magnetic field (B_⋆) and period (P_s) that cover the corresponding observed value ranges of YPs and MPs, we find the following main results:

• (a)  
  For the same particle injection rate, ${ \mathcal F }$, the cutoff energies, ${\epsilon }_{\mathrm{cut}}$, increase with the spin-down power, $\dot{{ \mathcal E }}$.
• (b)  
  For the same B_⋆ and P_s values (i.e., $\dot{{ \mathcal E }}$), the ${\epsilon }_{\mathrm{cut}}$ decrease with ${ \mathcal F }$ for high ${ \mathcal F }$ and saturate for low ${ \mathcal F }$.
• (c)  
  The optimum ${ \mathcal F }$ values that produce ${\epsilon }_{\mathrm{cut}}$ within the narrow value range observed by Fermi increase with $\dot{{ \mathcal E }}$ and can reach up to 10⁵ for both YPs and MPs. These optimum ${ \mathcal F }$ values also produce the observed L_γ behavior (${L}_{\gamma }\propto \dot{{ \mathcal E }}$ for low $\dot{{ \mathcal E }}$ and ${L}_{\gamma }\propto \sqrt{\dot{{ \mathcal E }}}$ for high $\dot{{ \mathcal E }}$).
• (d)  
  YP and MP models with $\dot{{ \mathcal E }}\lesssim {10}^{33}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and $\dot{{ \mathcal E }}\lesssim {10}^{32}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ respectively produce ${\epsilon }_{\mathrm{cut}}$ values that are always smaller than those observed by Fermi and close to the Fermi threshold, which provides a possible explanation for why we do not observe γ-ray YPs and MPs with smaller $\dot{{ \mathcal E }}$.
• (e)  
  The γ-ray efficiency shows a maximum (at ${ \mathcal F }\approx 1.5{{ \mathcal F }}_{\mathrm{GJ}}^{0}$) that depends on the inclination angle, α, and ranges from 13% to 35%. For higher ${ \mathcal F }$ values, the corresponding γ-ray efficiency drops considerably.
• (f)  
  The higher $\dot{{ \mathcal E }}$ and ${ \mathcal F }$, the narrower the peaks in the γ-ray light curves are.

For lower ${ \mathcal F }$ values (i.e., ${ \mathcal F }\lesssim 10{{ \mathcal F }}_{\mathrm{GJ}}^{0}$), additional features begin appearing (especially for high α values) in the model γ-ray light curves that are not always consistent with the observed ones. This indicates that the decrease of ${ \mathcal F }$ should not be uniform along all magnetic field lines but should be more focused along the separatrices since these regions are the main suppliers of the particles that eventually lie in the region near the ECS. Actually, this is totally consistent with the approach we have followed for the FIDO models where the decrease of the plasma conductivity σ was made only near the ECS and not in other places inside or outside the LC (Kalapotharakos et al. 2017). The above indicates that the number of particles that eventually enter the ECS region regulates the observed emission. For this reason, we have started producing kinetic models using a magnetic-field-line-dependent particle injection prescription. This approach allows much better (i.e., autonomous) control of the particle population that participates in and regulates the area at and around the ECS. We expect the γ-ray light curves of these models to be consistent with the Fermi ones even for lower particle injection rates.

In Kalapotharakos et al. (2017; FIDO models), we found a dependence of the plasma conductivity σ on $\dot{{ \mathcal E }}$; in the present study (kinetic models), we find a dependence of the global particle injection rate ${ \mathcal F }$ on $\dot{{ \mathcal E }}$ as well. Both these relations (σ versus $\dot{{ \mathcal E }}$ and ${ \mathcal F }$ versus $\dot{{ \mathcal E }}$) are imposed by the Fermi data and more specifically by the requirement that the models reproduce the observed ${\epsilon }_{\mathrm{cut}}$. Aside from the ${\epsilon }_{\mathrm{cut}}$ success, the optimum kinetic models reproduce the trends of the observed total γ-ray luminosity as a function of $\dot{{ \mathcal E }}$. Moreover, the FIDO models and the kinetic models for high particle injection rates provide very similar sky maps and γ-ray light curves that are consistent with the observed δ–Δ correlation while, as mentioned above, the γ-ray light curves for lower particle injection rates often indicate that a more careful treatment of the corresponding particle injection is needed. Nonetheless, the apparent $\sigma \,\mbox{--}\,{ \mathcal F }\,\mbox{--}\,\dot{{ \mathcal E }}$ relation uncovered by our studies connects fundamental macroscopic quantities, providing unique insight into the physical mechanisms behind the high-energy emission in pulsar magnetospheres.

Eventually, this successful description needs to be supported by the microphysical mechanisms that govern the particle creation processes, providing a complete justification for the aforementioned relations.

However, incorporating pair-creation microphysics (see Timokhin & Arons 2013) in global PIC simulations is not a trivial task because of the corresponding different length and time scales and the fact that the PIC energies are much smaller than those for real pulsars. Nonetheless, in trying to expand our studies toward this direction, we have started exploring not only global magnetosphere models with magnetic-field-line-dependent particle injection but also fully self-consistent models that incorporate a physically motivated particle injection. We will present our results in forthcoming papers.

We would like to thank an anonymous referee for helpful suggestions that improved the paper. This work is supported by the National Science Foundation under Grant No. AST-1616632, by the NASA Astrophysics Theory Program, by the NASA Astrophysics Data Analysis Program, and by Fermi Guest Investigator Program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Facility at NASA Ames Research Center and NASA Center for Climate Simulation (NCCS) at NASA Goddard Space Flight Center.