Modeling Gamma-Ray Light Curves with More Realistic Pulsar Magnetospheres

The launch of the Fermi Gamma-Ray Space Telescope has opened a new era in the study of pulsar physics. In the first year of operation of the Large Area Telescope (LAT), more than 40 new gamma-ray pulsars were discovered (Abdo et al. 2010). To date, more than 150 gamma-ray pulsars have been detected by the LAT (Ackermann et al. 2015), 117 of which are included in the Second Fermi Pulsar Catalog (Abdo et al. 2013). They are divided into three groups: millisecond pulsars, young radio-loud pulsars, and young radio-quiet pulsars. High-quality light curves, phase-averaged spectra, and phase-resolved spectra obtained by the LAT observations provide an opportunity for understanding the emission sites and magnetic field geometries of these gamma-ray pulsars.

The realistic structures of the pulsar magnetosphere still remain uncertain. Knowledge about the pulsar magnetosphere structures can be used to identify the potential sites of particle acceleration and gamma-ray emission. A vacuum dipole field is generally adopted in the early study of pulsar emission, because it has an exact analytical solution given by Deutsch (1955). Based on this field structure, different theoretical models have been developed to explain the observed pulsar emission. In these models, it is widely believed that particles are accelerated in the gap region where an accelerating electric field is created because of the deficit of charges. Gamma-ray emission is produced by the curvature or inverse-Compton radiation from high-energy particles accelerated in these gaps. Due to different emission zone locations, standard pulsar radiation models include the polar cap (PC; e.g., Ruderman & Sutherland 1975; Daugherty & Harding 1982), the slot-gap (SG) (e.g., Dyks & Rudak 2003; Dyks et al. 2004; Muslimov & Harding 2004), and the outer gap (OG; e.g., Cheng et al. 1986, 2000; Zhang & Cheng 1997, 2001; Zhang et al. 2004) models. These gap models have achieved great successes in explaining pulsar high-energy emissions and light curves (e.g., Watters et al. 2009; Romani & Watters 2010). However, the vacuum solution has no plasma; it is not able to reproduce any pulsar phenomena. It is well known that the pulsar magnetosphere should be filled with plasma (Goldreich & Julian 1969). In the presence of abundant plasma, all accelerating electric fields can be efficiently screened to form a force-free (FF) magnetosphere. The FF solution for an aligned rotator was first obtained by Contopoulos et al. (1999). The CKF solution consists of a closed field line region extending to the light cylinder (LC), an open field line region, and an equatorial current sheet beyond the LC. Moreover, the time-dependent simulations for the FF axisymmetric rotator also confirmed the closed–open CKF solution (e.g., Komissarov 2006; McKinney 2006; Timokhin 2006; Yu 2011; Parfrey et al. 2012; Cao et al. 2016a; Etienne et al. 2017).

Three-dimensional (3D) structures of the FF pulsar magnetosphere have been available in recent years (e.g., Spitkovsky 2006; Kalapotharakos & Contopoulos 2009; Pétri 2012). Moreover, these studies have been also extended to the general-relativistic FF regime that takes time–space curvature and frame-dragging effects into account (Pétri 2016; Carrasco et al. 2018) and the full magnetohydrodynamic regime that takes plasma pressures and inertial into account (Tchekhovskoy et al. 2013). The 3D FF solutions are similar to the CKF solution with an equatorial current sheet outside the LC. These simulations gave an impetus to the study of pulsar high-energy radiation with a more realistic field geometry instead of the vacuum dipole field. Gamma-ray light curves have been modeled using the FF magnetic field (Bai & Spitkovsky 2010b; Bogovalov et al. 2018). The FF solutions provide different pulse profiles due to the increased PC size compared to that of the vacuum dipole. It should be noted that the FF solutions cannot allow any accelerating electric fields. Therefore, they cannot provide any information about the sites of particle acceleration and radiation.

It is well known that the vacuum solution has no particle distributions, while the FF solution does not allow particle acceleration along magnetic field lines. Therefore, a more realistic pulsar magnetosphere should lie between the vacuum and FF limits. In fact, the plasma resistivity can produce a nonzero accelerating electric field along the magnetic field and a set of resistive solutions that smoothly bridges the gap between the vacuum and FF solutions that have been constructed based on the finite-difference time-domain approach (Kalapotharakos et al. 2012b; Li et al. 2012). Moreover, the resistive pulsar magnetospheres have been used to model the pulsar γ-ray spectra and light curves (Kalapotharakos et al. 2014; Brambilla et al. 2015). Also, particle-in-cell (PIC) methods with a self-consistent treatment between particles and fields are developed to simulate the structures of the pulsar magnetosphere (Chen & Beloborodov 2014; Philippov & Spitkovsky 2014; Belyaev 2015; Cerutti et al. 2015; Philippov et al. 2015; Brambilla et al. 2018; Kalapotharakos et al. 2018). Very recently, full PIC simulations also started to predict the pulsar light curves by including the radiation reaction (Cerutti et al. 2016; Kalapotharakos et al. 2018; Philippov & Spitkovsky 2018).

In a previous paper, we presented the structures of oblique pulsar magnetosphere with the uniform conductivity of a pseudospectral method (Cao et al. 2016b). We constructed a set of resistive solutions that smoothly bridges the gap between the vacuum and FF limits. In this paper, we focus on the the influence of the conductivity on gamma-ray light curves in dissipation pulsar magnetospheres. The light curves are produced by the geometric method and the particle trajectory method. As an application, we compare the predicted light curves with the gamma-ray light curves of the Vela pulsar observed by Fermi-LAT. The paper is organized as follows: In Section 2, we present the basic equations describing the structure of pulsar magnetosphere. In Sections 3 and 4, we present the method that we have used to produce the gamma-ray light curves. In Section 5, we apply the model to the Vela pulsar. Finally, a brief discussion and conclusions are given in Section 6.

The time-dependent Maxwell equations are given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{c}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-{\boldsymbol{\nabla }}\times {\boldsymbol{E}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{c}\displaystyle \frac{\partial {\boldsymbol{E}}}{\partial t}={\boldsymbol{\nabla }}\times {\boldsymbol{B}}-\displaystyle \frac{4\pi }{c}{\boldsymbol{J}},\end{eqnarray} \tag{ 2 }$$

with two initial conditions

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{E}}=4\pi {\rho }_{{\rm{e}}},\end{eqnarray} \tag{ 4 }$$

where ρ_e is the charge density and ${\boldsymbol{J}}$ is the current density. The structure of pulsar magnetosphere can be determined by the prescription for the current density ${\boldsymbol{J}}$. For the FF magnetosphere, the current density can be derived by the FF condition $({\boldsymbol{E}}\cdot {\boldsymbol{B}}=0)$ and the Maxwell equations as (Gruzinov 1999; Blandford 2002)

$$\begin{eqnarray}&&{\boldsymbol{J}}=c{\rho }_{{\rm{e}}}\displaystyle \frac{{\boldsymbol{E}}\times {\boldsymbol{B}}}{{B}^{2}}+\displaystyle \frac{c}{4\pi }\displaystyle \frac{({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}\times {\boldsymbol{B}}-{\boldsymbol{E}}\cdot {\boldsymbol{\nabla }}\times {\boldsymbol{E}}){\boldsymbol{B}}}{{B}^{2}},\end{eqnarray} \tag{ 5 }$$

which implies that the electric field parallel to the magnetic field component ${{\boldsymbol{E}}}_{\parallel }$ will be zero.

The FF approximations do not allow for any particle acceleration along the magnetic field, ${{\boldsymbol{E}}}_{\parallel }=0$. Therefore, they do not account for the particle acceleration and pulsed emission in the magnetosphere. A more realistic pulsar magnetosphere should allow for ${{\boldsymbol{E}}}_{\parallel }\ne 0$ in some regions of the magnetosphere. The resistive magnetosphere can allow for a nonzero parallel electric field component by involving a conductivity parameter σ. In fact, the resistive pulsar plasma is effectively collisionless. Therefore, the resistivity coming from Coulomb collisions can be practically neglected. We can define the current density for the resistive magnetosphere by a form of Ohm's law given by (Kalapotharakos et al. 2014)

$$\begin{eqnarray}&&{\boldsymbol{J}}=c{\rho }_{{\rm{e}}}\displaystyle \frac{{\boldsymbol{E}}\times {\boldsymbol{B}}}{{B}^{2}+{E}_{0}^{2}}+\sigma {{\boldsymbol{E}}}_{\parallel },\end{eqnarray} \tag{ 6 }$$

where

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

$$\begin{eqnarray}&&{E}_{0}{B}_{0}={\boldsymbol{E}}\cdot {\boldsymbol{B}},\quad {E}_{0}\geqslant 0.\end{eqnarray} \tag{ 8 }$$

The first term in Equation (6) is the drift current component perpendicular to ${\boldsymbol{B}}$, while the second term controls ${{\boldsymbol{E}}}_{\parallel }$ by the conductivity σ. The E₀ term in Equation (6) ensures the drift current will be subluminal. When the conductivity σ increases from σ = 0 to $\sigma \to \infty$, we expect to obtain a set of solutions that smoothly transition from the vacuum field to the FF field. It is noted that there is no unique prescription for the current density in the resistive electrodynamics. We checked the field structures by using the prescription for the current density given by Li et al. (2012). Our results are also qualitatively very similar to those of Li et al. (2012).

3.1. Magnetic Field Structures

The time-dependent Maxwell equations are solved by a spectral method in spherical coordinates. The radial coordinate r is expanded into the Chebyshev function. However, the angle coordinates θ and ϕ are expanded into the vector spherical harmonic functions. A spectral filter is used in all directions in order to ensure the stability of the algorithm and increase the convergent rate of the solution. For a detailed discussion about the pseudospectral algorithm, see Pétri (2012) and Cao et al. (2016a, 2016b). The computational domain extends from the inner boundary r_min = 0.2 r_L to the outer boundary r_max = 3 r_L. A resolution of N_r × N_θ × N_ϕ = 128 × 32 × 64 is necessary to get a good accuracy. The magnetic field is initialized to be a rotating vacuum dipolar with magnetic inclination angle α = {0°, 15°, 30°, 45°, 60°, 75°, 90°}. We impose inner boundary conditions at the stellar surface with a rotating electric field ${\boldsymbol{E}}=-({\boldsymbol{\Omega }}\times {\boldsymbol{r}})\times {\boldsymbol{B}}/c$. A characteristic compatibility method is used in order to prevent the inward reflection from the outer boundary. We let the system evolve for several rotational periods to reach a final steady solution. We perform a series of simulations with uniform conductivities σ = {0.1, 0.3, 1, 3, 5, 10, 30, 60} Ω. Magnetic field structures in the Ω–μ plane for magnetic inclination α = 60° with increasing conductivity are shown in Figure 1. The normalized Poynting flux L/L_aligned as a function of radius r for magnetic inclination α = 60° with different σ values at time t = 3 P are shown in Figure 2. We can see that the Poynting flux increases with increasing conductivity. To demonstrate the convergence of the solution, the normalized Poynting flux for σ = 30 Ω at time t = 4 P is also shown as the red dashed line. The Poynting flux is nearly constant when the time t > 3 P, the system will relax to the steady-state solution. In Figure 3, we show the normalized Poynting flux L/L_aligned as a function of radius r for magnetic inclination α = 0° with different σ values at time t = 3 P. It is found that our results still have some numerical dissipations due to the low resolution used in the presented 3D simulations. A higher resolution is necessary to resolve more fine plasmoids in the current sheet. Our simulations only present the lower limit on resolution to get convergence. We can see a smooth transition for magnetic field structures and the Poynting flux from the vacuum limit to the FF limit with increasing conductivity. For a detailed description about the structures of resistive magnetosphere, see Cao et al. (2016b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The PC shape is defined by the footprints of the last closed field lines on the stellar surface. The last closed field lines can be found by checking whether the field lines close inside or outside the LC. It has been known that the PC shape in the vacuum case is different from that in the FF case (e.g., Bai & Spitkovsky 2010b; Kalapotharakos et al. 2014). Here, a third-order Runge–Kutta integration is used to find the last closed field lines in the lab frame, and the bisection method is used to find the magnetic colatitude ${\theta }_{m}^{\mathrm{rim}}$ of the rim of the PC for the fixed magnetic azimuth ${\phi }_{m}^{\mathrm{rim}}$. Open volume coordinates are be defined by (r_ovc, ϕ_m) on the PC. r_ovc ≥ 1 represent the closed field lines and r_ovc < 1 represent the open field lines. The PC shapes on a sphere of radius 0.2 r_L for magnetic inclination α = 60° with different values of σ are shown in Figure 4. The PC shapes for the low σ value are similar to that for the vacuum field, and a notch appears on the PC rim. In the resistive solution, on the one hand, the polar shapes are larger than that of the vacuum solution because of the larger open magnetic flux. On the other hand, the PCs are shifted toward the trailing side. The notches on the PCs disappear at high σ values and the PC shapes are more circular.

Zoom In Zoom Out Reset image size

3.2. Sky Maps and Light Curves

The light curves can be computed by the geometry of pulsar emission models for a given magnetosphere structure. In the SG model, all accelerating electric fields are screened by the pair cascades above the pair-formation fronts. There is a narrow gap with ${E}_{\parallel }\ne 0$ between the last closed lines and the pair-formation fronts. The particles continue to accelerate and radiate from the neutron star surface to high altitude. The OG model is a vacuum gap that also forms between the null-charge surfaces and the last closed field lines. The outer gap extends from the null-charge surfaces to the LC, and the gap width is limited by the screening of ${E}_{\parallel }\ne 0$ by pair cascades. After the current through the gap is taken into account, it has been shown that the inner boundary of the outer gap is shifted inward from the null-charge surface (Takata et al. 2004; Hirotani 2006). Therefore, the gaps would exist inside the LC of the pulsar magnetosphere, the details of E_∥ depending on the pair cascade process that has not been taken into account in the simulation of pulsar magnetosphere so far.

To explore how the magnetic field structures and offset PCs influence γ-ray pulsar light curves, we use a geometric method to produce the γ-ray light curves based on the uniform σ model. We assume that the γ-ray emission originates in a layer of width ω = 0.1 along the last closed field lines between r_ovc = 0.9 and r_ovc = 1. The emission region extends from the neutron star surface to the LC, and the γ-ray emission is restricted to lie within a cylindrical radius ρ_max ≲ r_L. We assume the uniform emissivity along the field lines in the comoving frame (CF). The photon direction is assumed to be tangential to the magnetic field in the CF, obtained through a Lorentz transformation from the inertial observer's frame (IOF; Bai & Spitkovsky 2010a). The emitted photon is collected in the sky map in viewing angle ζ and the observed phase Φ, taking into account the aberration and time-delay effects. The light curves are then obtained by a cut through this sky map at constant ζ.

In Figure 5, we show the sky maps and the corresponding light curves from the geometric method in different inclination angles and view angles with different magnetic field configurations. We see narrow double-peak profiles at relative large α and ζ, which are very similar to the observed γ-ray light curves. We see that the light curves for the low σ values are very close to that of the vacuum dipole, but the peaks are shifted in phase. As the conductivity σ increases, the peak phase is shifted to the lager phase relative to the magnetic pole. Similar results are also found by Harding et al. (2011) and Kalapotharakos et al. (2012a). Overall, there is a significant progression in the LC shapes as σ increases. This is because that the magnetospheres with high conductivity have more sweepback field lines, which produce the larger shift of PC and thus cause a larger phase lag of the light curve peaks.

Zoom In Zoom Out Reset image size

4.1. Real Particle Trajectories

We study the contribution of the curvature radiation to the γ-ray light curves by including the accelerating electric fields provided by the solutions themselves. We define the trajectory of particles in IOF by (Kalapotharakos et al. 2014)

$$\begin{eqnarray}&&{\boldsymbol{v}}=\left(\displaystyle \frac{{\boldsymbol{E}}\times {\boldsymbol{B}}}{{B}^{2}+{E}_{0}^{2}}+f\displaystyle \frac{{\boldsymbol{B}}}{B}\right)c,\end{eqnarray} \tag{ 9 }$$

where the first term in Equation (8) is a drift velocity component, while the second term is a component parallel to the magnetic field. By requiring that the condition v ≃ c and that the motion of the particle is outward, the spatial distribution of f can be uniquely determined by Equation (9). We assume that the direction of photon emission, ${{\boldsymbol{\eta }}}_{\mathrm{em}}$, is along the direction of particle motion ${\boldsymbol{\beta }}={\boldsymbol{v}}/c$. Therefore, the direction of photon emission in the IOF is determined by

$$\begin{eqnarray}&&{\mu }_{\mathrm{em}}={\beta }_{z},\quad {\phi }_{\mathrm{em}}=\mathrm{atan}\left(\displaystyle \frac{{\beta }_{y}}{{\beta }_{x}}\right),\end{eqnarray} \tag{ 10 }$$

where the view angle $\zeta =\mathrm{acos}({\mu }_{\mathrm{em}})$.

The observe phase is determined by include the rotation and time-delay correction

$$\begin{eqnarray}&&{\rm{\Phi }}={\phi }_{\mathrm{rot}}-{\phi }_{\mathrm{em}}-{{\boldsymbol{r}}}_{\mathrm{em}}\cdot {{\boldsymbol{\eta }}}_{\mathrm{em}}/{r}_{{\rm{L}}},\end{eqnarray} \tag{ 11 }$$

where ϕ_rot = Ωdt is the rotation phase, ϕ_em is the azimuthal angle of the emitted photon, and ${{\boldsymbol{r}}}_{\mathrm{em}}$ is the location of the emitted photon.

Under the above assumptions, we can calculate the trajectories of radiating particles passing through each magnetospheric point in the computational domain, which also allows for the determination of curvature radius R_CR along each particle trajectory. We assume that the charge particles (e⁻, e⁺) are uniformly distributed on the PC. The e^± pairs are then injected from the stellar surface with a small Lorentz factor (γ < 100). The Lorentz factor γ of the radiating particle along each trajectory is calculated by including the influence of the accelerating electric field and curvature radiation loss. Specifically, the γ value along each trajectory is integrated by the expression

$$\begin{eqnarray}&&\displaystyle \frac{d\gamma }{{dt}}=f\displaystyle \frac{{q}_{{\rm{e}}}{{cE}}_{\parallel }}{{m}_{{\rm{e}}}{c}^{2}}-\displaystyle \frac{2{q}_{{\rm{e}}}^{2}{\gamma }^{4}}{3{R}_{\mathrm{CR}}^{2}{m}_{{\rm{e}}}c},\end{eqnarray} \tag{ 12 }$$

where q_e and m_e are the electron charge and rest mass, respectively. The first term in Equation (11) is the energy-gain rates of the particles due to the accelerating electric field, and the second term is the energy-loss rates due to curvature radiation. Then, we calculate the characteristic energy of curvature radiation ${E}_{{\rm{c}}}=\tfrac{3}{2}c\hslash \tfrac{{\gamma }^{3}}{{R}_{{\rm{CR}}}}$ and its bolometric luminosity ${L}_{{\rm{bol}}}=\tfrac{2}{3}{q}_{{\rm{e}}}^{2}c\tfrac{{\gamma }^{4}}{{R}_{{\rm{CR}}}^{2}}$ along each trajectory. We integrated each particle trajectory from the neutron star surface up to r = 2.5 r_L. We construct the γ-ray light curve by collecting the bolometric emission with E_c > 0.1 GeV from all the emitting particles in sky maps. For a detailed description of the trajectory and radiation of particles, see Kalapotharakos et al. (2014).

4.2. Sky Maps and Light Curves

Kalapotharakos et al. (2014) studied the γ-ray light curves using dissipative pulsar magnetospheres. It is found that the uniform σ model cannot well match the distribution of γ-ray peak separation and the radio lag observed by Fermi-LAT. Further, they found that a significant improvement can be achieved by applying an FF regime inside the LC and a dissipative regime (FIDO) outside the LC. In fact, we also find that the uniform σ model cannot explain the pulsar light curves observed by Fermi-LAT. Therefore, we adopt a model that is similar to the FIDO model introduced by Kalapotharakos et al. (2014). Kalapotharakos et al. (2014) used a σ linear approximation for E_∥ based on the FF solution at high σ values. In this work, we adopt a more accurate treatment for E_∥. We employ a very high conductivity (σ = 60 Ω) within the LC and a finite conductivity outside the LC, we then remove any E_∥ solutions within the LC. We assume that all the γ-ray emission comes from the outer magnetosphere outside the LC. In Figure 6, we shown the E_∥ distribution in the Ω–μ plane for magnetic inclination α = 60° with σ = 60 Ω within the LC and with σ = 30 Ω outside the LC. The field line structure is very similar to that of the FF solution. We see the strong parallel electric field along the equatorial current sheet outside the LC. We also note that the parallel electric distribution is very similar to that of the global PIC simulation (Kalapotharakos et al. 2018; Philippov & Spitkovsky 2018).

Zoom In Zoom Out Reset image size

In the following study, we assume the standard pulsar parameters for pulsar period P = 0.1 s and surface magnetic field B_⋆ = 10¹² G. We calculate the direction of photon emission and the corresponding bolometric emission along each particle trajectory. This information allowed us to produce γ-ray light curves by collecting all the photons with E_c > 0.1 GeV in sky maps. In Figures 7–8, we show the sky maps and the corresponding light curves in different inclination angles and observer viewing angles for the magnetosphere with a near FF regime within the LC and with σ = 0.3 Ω and 30 Ω outside the LC. For σ = 0.3 Ω outside the LC, the light curves generally display only one broad peak. We see the double-peak light curves as α and ζ increase, but the peaks seem to be broad. For very high σ values outside the LC, we see the significant double-peak profiles, except for some additional secondary peaks. There are some cases where the light curves are narrow, which are very similar to those observed by Fermi-LAT. We find that all the emission originates in a thin layer with r_ovc = 0.9–1 near the equatorial current sheet outside the LC.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The observed pulsar light curves potentially provide an important diagnostic for pulsar magnetosphere structure. As an application, we compare the light curves predicted from the two methods with those observed at >0.1 GeV energies for the Vela pulsar. The results are shown in Figures 9 and 10. For the geometric method, we find that the predicted light curve with σ = 0.3 Ω provides a better match to the Fermi observed data for the Vela pulsar. The double-peak profile with a phase separation of ∼0.43 can be better reproduced by our model with a low value of σ. In fact, the first peak phase lag of the high σ model is too large to explain the observed data for the Vela pulsar. We suggest that the magnetosphere with the low σ value is favorable for the geometric method. For the particle trajectory method, we use the measured parameters of the Vela pulsar: P = 0.089 s, B = 3.4 × 10¹² G. We find that the magnetosphere with a near FF regime within the LC and with σ = 30 Ω outside the LC can better produce the observed data of the Vela pulsar. Our results indicate that the magnetosphere with a near FF regime within the LC and a high σ value outside the LC is favorable for the particle trajectory method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we explore the γ-ray emission patterns and light curves in dissipation pulsar magnetospheres. The γ-ray light curves are produced by the geometric method and particle trajectory method. For the geometric method, we assume that the γ-ray emission comes from a finite-width layer along the last closed lines; we generate the γ-ray light curves by assuming the uniform emissivity along the field lines in the CF. We find that there is a significant progression in the γ-ray peak phase with increasing σ. For the particle trajectory method, we consider the spatial distribution of σ by assuming a very high conductivity within the LC and a finite conductivity outside the LC. Then, we use the field structure of these models to define realistic trajectories of radiated particles. Assuming that all the γ-ray emission comes from the outer magnetosphere outside the LC, we compute the Lorentz factors of radiating particles and the characteristic energy of curvature radiation along each particle trajectory under the influence of both the accelerating electric field and curvature radiation loss. We produce the γ-ray light curves by collecting the bolometric emission with E_c > 0.1 GeV from all the emitting particles. We find that the light curve shapes are very sensitive to the value of σ, and all the γ-ray emission is produced in regions near the equatorial current sheet for very high σ values.

As an application, we compare the modeling light curves from the two methods with observed light curves at >0.1 GeV energies for the Vela pulsar. For the geometric method, we find that the magnetosphere with the low σ value provides a better match to the observed data for the Vela pulsar. For the particle trajectory method, we find that the magnetosphere with a near FF regime within the LC and a high σ value outside the LC can better produce the light curves of the Vela pulsar. We note that the particle trajectory method uses the parallel electric fields that are self-consistent with the magnetic fields of the magnetosphere. Therefore, the results from the particle trajectory method are more physically motivated. In fact, the microphysical processes of the pair production and particle acceleration are not self-consistently included in the simulation of the resistive magnetosphere. However, the parallel electric field distribution from the high σ model outside the LC is only confined near the equatorial current sheet, which is very similar to that of the global PIC simulation (Kalapotharakos et al. 2018; Philippov & Spitkovsky 2018). This gives a microphysical explanation for the origin of the parallel electric field from the high σ model outside the LC. Our results from the particle trajectory are also similar to those from the PIC simulation. The particle trajectory method also allows for the derivation of an energy spectrum, which can be compared directly with the Fermi observations. Therefore, we will explore the γ-ray energetic and spectral properties in future work. Further, we will improve our code by identifying different field lines at each time step of the simulation. This will allow us to apply σ along specific magnetic field lines. This technique will be helpful for us to understand the underlying pair production process in the gaps or current sheets.

We thank the anonymous referee for valuable comments and suggestions. We would like to thank Jérôme Pétri, Gabriele Brambilla, Xuening Bai, and Li Zhang for some useful discussions. We acknowledge the financial support from the National Natural Science Foundation of China 11573060 and 11661161010, the National Science Foundation of China 11673060, the National Science Foundation of China 11871418, and Fund for Science and Technology Innovation team in universities at Yunnan. The simulation presented in this work is performed at HPC Center, Yunnan Observatories, CAS, China.