The Energy-dependent Gamma-Ray Light Curves and Spectra of the Vela Pulsar in the Dissipative Magnetospheres

We study the pulsar energy-dependent γ-ray light curves and spectra from curvature radiation in the dissipative magnetospheres. The dissipative magnetospheres with the combined force-free (FFE) and Aristotelian are computed by a pseudo-spectral method with a high-resolution simulation in the rotating coordinate system, which produces a near-FFE field structure with the dissipative region only near the equatorial current sheet outside the light cylinder. We use the test-particle trajectory method to compute the energy-dependent γ-ray light curves and phase-average and phase-resolved spectra by including both the accelerating electric field and radiation reaction. The predicted energy-dependent γ-ray light curves and spectra are then compared with those of the Vela pulsar observed by Fermi. Our results can generally reproduce the observed trends of the energy-dependent γ-ray light curves and spectra for the Vela pulsar.


INTRODUCTION
Pulsars are rapidly rotating neutron star with extremely strong magnetic field, which can convert a substantial fraction of rotational kinetic energy into particle acceleration and radiation.These objects can emit the broadband electromagnetic emission from the radio to very high-energy GeV, and even TeV bands.The radiation from pulsars originates from the high-energy particles accelerated by unscreened electric fields in the magnetosphere.These high-energy particles flow out along the open field lines and then produce the broadband electromagnetic spectra by the synchrotron, curvature, and inverse-Compton processes.The Fermi Gamma-Ray Space Telescope launched in 2008 has revolutionized our knowledge of pulsar physics with the detection of more than 270 γ-ray pulsars.The Second Fermi Pulsar Catalog (2PC) contained 117 γ-ray pulsars with the measurement of good-quality light curves and spectra (Abdo et al. 2010(Abdo et al. , 2013)).They are divided into young radio-loud, young radio-quiet, and millisecond pulsars.Fermi pulsars usually show the double-peak γ-ray light curves with ∼ 0.5 peak separation and the first peak typically lags behind the radio peak.An apparent feature is that the relative ratio of first to second γ-ray peaks decrease as the the photon energy increases for Vela and Geminga.The observed γ-ray spectra have the exponentially cutoff power-law shape with the cutoff energy in the range of 1-5 GeV.Good-quality γ-ray data from Fermi observation provide a wealth of information to explore the radiation mechanism and particle acceleration in the magnetosphere.However, it is still unclear where the particles in the magnetosphere are accelerated and how the emitting photons are produced.The particle acceleration and pulsed emission in the magnetosphere is connected to the structure of the global pulsar magnetosphere.It is compulsory to self-consistently solve for Maxwell equations and particle dynamics and radiation to model the global pulsar magnetosphere and explain the observed pulsed emission.
Significant progresses towards the self-consistent modeling of global pulsar magnetosphere have been made over the last decades.It is expected that the pulsar magnetospheres are filled with plasmas by the pair production (Goldreich & Julian 1969).A good first guess about plasma-filled magnetosphere is referred to as the force-free electrodynamic, which corresponds to the zeroth order approximation of global pulsar magnetosphere.Force-free solutions have recently become available by the numerical simulation pioneered by Contopoulos et al. (1999, hereafter CKF) for the aligned rotator and then by Spitkovsky (2006) for the oblique rotator.The force-free solutions are later studied by the time-dependent simulations with the spectral method (Parfrey et al. 2012;Pétri 2012;Cao et al. 2016a;Pétri 2016) and the finite-difference method (Kalapotharakos & Contopoulos 2009;Contopoulos & Kalapotharakos 2010).All these force-free simulations reveal the existence of an equatorial current sheet outside the LC, which is suggested to be an alternative site of the pulsed emission in addition to the "gap" regions (Bai & Spitkovsky 2010;Harding & Kalapotharakos 2015;Brambilla et al. 2015;Bogovalov et al. 2018;Harding et al. 2018).The force-free solutions are thought to be closest to realistic pulsar magnetospheres, but no particle acceleration and radiation is allowed in the magnetosphere.Therefore, some dissipation mechanism should be included to address the problem of particle acceleration and radiation production in the magnetosphere.The resistive magnetosphere by introducing a conductivity parameter are then developed to selfconsistently include the dissipative regions in the magnetosphere, which can switch the magnetospheric solutions from the vacuum to force-free limits (Li et al. 2012;Kalapotharakos et al. 2012a;Cao et al. 2016b).The resistive solutions have been also used to model the pulsar light curves and spectra by including the accelerating electric field from the simulations (Kalapotharakos et al. 2012b(Kalapotharakos et al. , 2014(Kalapotharakos et al. , 2017;;Cao & Yang 2019;Yang & Cao 2021).However, it is not clear about the origin and physical motivation of these resistivities by involving an arbitrary conductivity.The aforementioned fluid description offer us a good knowledge of global pulsar magnetosphere with the macroscopic current and charge densities.However, the fluid methods can not model the source of particles that provide the current and charge densities.Therefore, a full self-consistent PIC (Particle-In-Cell) method is developed to model the global pulsar magnetosphere by including the particle dynamics and radiation reactions (Philippov & Spitkovsky 2014;Chen & Beloborodov 2014;Philippov et al. 2015;Belyaev 2015;Cerutti et al. 2015;Kalapotharakos et al. 2018;Brambilla et al. 2018).The PIC methods are also connected to the observational signature by simultaneously extracting the pulsar light curves from the simulations (Cerutti et al. 2016;Philippov & Spitkovsky 2018;Kalapotharakos et al. 2018).However, the Lorentz factors from the PIC simulation are much smaller than ones in the real pulsars, which is not enough to explain the observed γ-ray photons.
Aristotelian electrodynamics is a good approximation between the resistive model and PIC model, which can include the back-reaction of emission onto particle dynamics and allow for a local dissipation where the forcefree condition is violated.The AE method has been introduced to model the structure of pulsar magnetospheres (Gruzinov 2012(Gruzinov , 2013;;Contopoulos et al. 2016).A clever method by combining FFE and AE is also proposed to model the dissipative pulsar mag- netosphere for the aligned rotator (Pétri 2020a) and for oblique rotators (Cao & Yang 2020;Pétri 2022).Recently, Cao & Yang (2022) perform the high-resolution simulations of the combined FFE and AE magnetospheres to better capture the dissipative regions near the current sheet outside the LC.The pulsar light curves and spectra are also computed by using the accelerating electric field from the simulations.We confirm that the the particle acceleration and the γ-ray radiation is produced near the current sheets outside the LC.Moreover, the predicted light curves and spectra can general reproduce the double-peak γ-ray profiles and GeV spectral cutoff energy in agreement with the Fermi observations for the high pair multiplicity.However, the light curves and spectra prediction are not directly compared with the Fermi observed γ-ray data in this study.In this paper, we compute the energy-dependent light curves, phase-average and phase-resolved spectra by expanding the study of Cao & Yang (2022) and directly compare the predicted energy-dependent light curves and spectra with the Fermi observation for the Vela pulsar.

DISSIPATIVE PULSAR MAGNETOSPHERES
The structure of pulsar magnetospheres are obtained by solving the time-dependent Maxwell equations in the rotating coordinate system (Muslimov & Harding 2005;Pétri 2020b) with the 3D spectral method (Cao & Yang 2022), where J is the current density, ρ e is the charge density, V rot = Ω × r is the corotating velocity.The current density with the combined FFE and AE is implemented to the Maxwell equations by introducing the pair multiplicity κ (Cao & Yang 2020) where B 0 and E 0 are the magnetic and electric field in the frame in which E and B are parallel.E 0 is denoted as E acc in the rest of the paper, which represents the effective accelerating electric component.|ρ e | can be interpreted as the charge density of the primary electrons, κ| ρ e | can be interpreted as the charge density of the secondary pairs from the pair cascades.Therefore, the pair multiplicity κ is physically related with the pair cascade processes in the magnetosphere.The quantities B 0 and E 0 is defined by the Lorentz invariants relations The magnetic field is initialized as an oblique vacuum dipole magnetic field with magnetic inclination angles α from 0 • to 90 • with a interval of 5 • .The inner boundary condition at the stellar surface is enforced with a rotating electric field E = −(Ω × r) × B. A non-reflecting boundary condition is implemented to prevent spurious reflection from the outer boundary.A high resolution with N r × N θ × N ϕ = 129 × 64 × 128 is used to obtain the accurate magnetospheric solutions from the stellar surface r = 0.2 r L to r = 3 r L .The pair multiplicity is set to κ = {0, 1, 3}.The dissipative pulsar magnetospheres with the combined FFE and AE is computed by applying the force-free description where E ≤ B and the AE description where E > B. It is expected that our model will locate the dissipative region at the current sheet outside the LC where E > B.
Figure 1 shows the structure of magnetic field line and the distribution of the accelerating electric field E acc for a 65 • rotator with the pair multiplicity κ = 3 in the x−z plane.We see that the combined FFE and AE magnetosphere tends to the force-free solution and the dissipative region is only confined to near the current sheet outside the LC.For the higher κ values, we can expect that the magnetospheric solutions keep the near force-free structures and the accelerating electric fields also decrease as the pair multiplicity κ increases (see Cao & Yang 2022).The numerical integration for the higher κ value become cumbersome because of the stiff nature of the pair multiplicity term.Therefore, the E acc value at κ > 3 is determined by the scaling relation where E acc,0 is the accelerating electric field at κ 0 = 3.

CURVATURE RADIATION
The particle velocity in the dissipative magnetosphere is defined by (Gruzinov 2012;Pétri 2020b;Cai et al. 2023) where the two signs correspond to positrons and electrons, they follow a differen trajectory that only depends on the local electric and magnetic field.The Lorentz factor along particle trajectory is integrated by which takes into account the influence of the local accelerating electric field and the curvature radiation losses.
The photon spectrum of the particle curvature radiation with Lorentz factor γ is given by where R C is the curvature radius of particles, x = E γ /E cur , E γ is the radiation photon energy, E cur = 3 2 cℏ γ 3 RCR is the characteristic energy of the curvature radiation photon, and the function F (x) is given by The particles are randomly injected at the stellar surface with the small Lorentz factor.We integrate the particle trajectory from the the neutron star surface up to r = 2.5 r L by equation ( 8).The Lorentz factor along each trajectory is then computed by equation ( 9).We assume the direction of the photon emission along the direction of particle motion, the direction of the photon emission and the curvature radiation photons are then computed along each trajectory.The sky maps are produced by collecting all curvature radiation photons in the (ζ,ϕ) plane.The energy-dependent light curves are obtained by accumulating all emitting photons at a given energy range for a constant viewing angle ζ.The phase-averaged spectra are obtained by averaging all emitting photons in energy over the observed phase ϕ for a constant ζ.The phase-resolved spectra are obtained by collecting all emitting photons in energy at a given ϕ range for a constant ζ.

RESULTS
The light curves can provide an diagnostic for pulsar magnetosphere geometry, while the energy spectra can both explore the geometry and location of the emission regions in the magnetosphere.The combined light curves and spectra can put a stronger constraint on the location of the particle acceleration and the geometry of the emission regions in the magnetosphere.We simultaneously fit the energy-dependent light curves, phaseaveraged spectra and phase-resolved spectra of the Vela pulsar to constrain the model parameters.We show the sky maps and the predicted γ-ray light curves at > 0.1 GeV energies for the Vela pulsar in Figure 2. We can see that the predicted sky maps show the two bright caustics that are characteristic of radiation from the current sheet.Our sky maps show that the emission pattern from the current sheet is non-uniform, which depends on the field line geometry and the E acc distribution in the current sheet.On the one hand, the emission from the current sheet is concentrated in the same region of the sky map due to the "stagnation" caustic effects in the near force-free magnetosphere (Bai & Spitkovsky 2010).On the other hand, the E acc distribution in the current sheet is non-uniform, the E acc is more confined to the current sheet in the range of 1−2 R L , and the E acc varys with radius and azimuth and become very weak beyond the 2 R L .Therefore, we can not see an extended geometric emission pattern with ± 65 • above/below the equatorial plane in sky map due to the non-uniform emission from the current sheet.In fact, our emission pattern from the current sheet is also very similar to those of the PIC model (Philippov & Spitkovsky 2018;Kalapotharakos et al. 2018).

Cao & Yang
In our model, the field line geometry and the E acc distribution determines the light curve shape, while the E acc value determines the peak energy of the γ-ray spectrum, which can be adjusted by the free parameter κ.The free parameter κ is chosen to adjust the E acc value to match the Fermi γ-ray spectra with the real pulsar parameters.The model parameters are B * = 4×10 12 G, P = 0.089 s, κ = 17, α = 65 • and ζ = 64 • .The parameter κ = 17 is only chosen to obtain the suitable E acc value to match the Fermi γ-ray spectra with the real Vela pulsar parameters of B * = 4 × 10 12 G and P = 0.089 s.The same results can also be obtained by using a simulated pair multiplicity κ = 3 and a slightly lower surface magnetic field B * = 7 × 10 11 G.The pair multiplicity κ = 17 can not be understood as the real Veal pulsar parameter in our model, which is only physically related to the pair cascade efficiency in "killing" E acc .We find that the predicted light curve can well match the Fermi observed γ-ray data for the Vela pulsar.The peak phase and the peak separation of the Vela pulsar can be well reproduced by our model.Our model can better match the γ-ray light curves of the Vela pulsar compared to those in Harding et al. (2021) and Barnard et al. (2022).They can not provide a good match to the peaks phases of the Vela pulsar and the predicted first peak phase is larger than the observed one.We show the predicted energy-dependent γ-ray light curves and compare them with those of the Vela pulsar observed by Fermi in Figure 3.We observe that the flux of the first peak relative to the second one decreases towards the higher photon energies and the peak width also narrow with increasing photon energies.Our model can well reproduce the decreasing ratio of the first peak to the second one (P1/P2) and the narrow of the peak width.Our results indicate that the highest-energy photons originate from the second peak of the Fermi light curves.The predicted general trend of the light curves with the energy evolution is in good agreement with the Fermi observed one.We show the predicted phaseaveraged spectra and the phase-resolved spectra of the first and second peaks for the Vela pulsar in Figure 4. We see that the Fermi observed phase-averaged spectra can well be explained by the curvature radiation from the current sheet outside the LC.Our model also can well reproduce the Fermi observed phase-resolved spectra of the first and second peaks for the Vela pulsar.We find that the phase-resolved spectra of the second peak have a higher flux and a larger spectral cutoff compared to the first peak one.This is because that the light curves of the second peak survive at the higher energy compared to the first peak one.We show the sky map of R C and the variation of R C with the pulsed phase in Figure 5.We see that the sky pattern of R C is similar to that of the caustic emission, because the variation of R C closely reflects the variations of the caustic emission.We also find that the curvature radius in the second peak is the systematically larger than that in the first peak.The cutoff energy E cut of the curvature radiation in the radiation reaction limit is related with the accelerating electric field and the curvature radius by which scales with R 1/2 C .Therefore, we can expect a large E cut value in the second peak due to the systematically larger curvature radius than that in the first peak.Therefore, the decreasing P1/P2 ratio with increasing energy can be explained by the systematically larger curvature radius of the second peak than the first one.A similar result is also found in the force-free magnetosphere by Harding et al. (2021) and Barnard et al. (2022).

DISCUSSION AND CONCLUSIONS
We explore the properties of the pulsar energydependent γ-ray light curves and spectra in the dissipative magnetospheres.The dissipative magnetospheres with the combined FFE and AE are constructed by applying the force-free description where E ≤ B and the AE description where E > B. We compute the dissipative magnetospheres by a pseudo-spectral method with the high-resolution simulations in the rotating coordinate system.Our simulations show that the dissipative magnetospheres have the near force-free field structure and the dissipative region is confined to only near the equatorial current sheet outside the LC.We compute the energy-dependent γ-ray light curves, phase-average and phase-resolved spectra by the test particle trajectory method in the dissipative magnetosphere, taking into account both the local accelerating electric field and radiation reaction.We can generally reproduce the Fermi observed energy-dependent γ-ray light curves and spectra of the Vela pulsar.The decreasing ratio of the first peak to the second one with increasing energy can well reproduced by our model, which is attributed to the systematically larger local curvature radius and the higher spectral cutoff in the second peak.Brambilla et al. (2015) studied the energy-dependent γ-ray light curves of the Vela pulsar by the FIDO (Forcefree Inside Dissipative Outside) model.They can reproduce the observed trends of the light curve with the energy evolution for the Vela pulsar.However, an approximate accelerating electric field is used based on the corresponding force-free solutions in their study, which is not self-consistently determined by the simulation them-selves.Moreover, the predicted phase-averaged and phase-resolved spectra is not directly compared with the observed ones.The energy-dependent light curves and spectra of the Vela pulsar is also explained by an extended slot gap and current-sheet model with the constant accelerating electric field in the force-free magnetospheres (Harding et al. 2021;Barnard et al. 2022).However, the predicted light curve with the constant accelerating electric field can not well match the peaks phases of the Vela pulsar, because the accelerating electric field is the function of the altitude and azimuth.Our results can provide a better match to the peaks phases of the Vela pulsar by using the accelerating electric field from the simulation themselves.Our MHD simulations show that the acceleration electric fields in the current sheet decrease with increasing pair multiplicity.However, the particle accelerations in the current sheet not only depend on the plasma supply from the pair multiplicity but also the reconnection rate in the current sheet.Therefore, it is necessary to further study magnetic reconnection acceleration in the current sheet by the PIC method.The observed bridge emission of the Vela pulsar can not be explained by our model.It is suggested that the the bridge emission may originate in the inner regions of the magnetosphere (Barnard et al. 2022).We will further explore the the bridge emission of the Vela pulsar by introducing the local dissipation not only in the current sheet region but also in the inner region.The study in the Fermi γ-ray band is not enough to distinguish between different radiation mechanisms and radiation locations in the pulsar magnetosphere.A future study on the multi-wavelength light curves and spectra can better constrain the radiation mechanisms and localize the photon production sites in the magnetosphere.We will use the presented dissipative magnetospheres to present the study of the multi-wavelength light curves and spectra in a forthcoming work.

Figure 1 .
Figure 1.The structure of magnetic field line and the distribution of the accelerating electric field for a 65 • rotator with the pair multiplicity κ = 3 in the x−z plane.

Figure 4 .
Figure 4.The predicted γ-ray phase-averaged spectra and the phase-resolved spectra for the Vela pulsar with the same model parameters as in Figure 2. The Fermi observed data is taken from Abdo et al. (2013).

Figure 5 .
Figure 5.The sky map of RC and the variation of RC with the pulsed phase for the same model parameters as in Figure 2.