Modeling the Multiwavelength Radiation Properties in Pulsar Dissipative Magnetospheres

We explore the multiwavelength radiation properties of the light curves and energy spectra in the dissipative magnetospheres of pulsars. The dissipative magnetospheres are simulated by the pseudo-spectral method with the combined force-free and Aristotelian electrodynamics, which can produce self-consistent accelerating electric fields mainly distributed in the equatorial current sheet outside the light cylinder. The multiwavelength light curves and spectra are computed by using the multiple emission mechanisms of both the primary particles accelerated by the accelerating electric fields in the equatorial current sheet and the secondary pairs with an assumed distribution spectrum. We then compare the predicted multiwavelength light curves and spectra with the observed data from the Crab, Vela, and Geminga pulsars. Our modeling results can systematically reproduce the observed trends of the multiwavelength light curves and the spectra for these three pulsars well.


INTRODUCTION
The extreme circumstances of rapidly rotating and highly magnetized pulsars make it possible for them to be ideal accelerators in the Universe.They can emit very stable and regular electromagnetic radiation from accelerated high-energy particles in the radio up to the γ-ray band.In the past 15 years, the Fermi Large Area Telescope has provided us with abundant high-quality γ-ray light curves and spectra in the 100 MeV−300 GeV range.The Fermi γ-ray light curves generally show double-peak profiles and energy-dependent evolutionary patterns.The Fermi γ-ray spectra usually follow an exponential power-law shape with a cutoff energy around a few GeV (Abdo et al. 2013).Recently, ground-based air-Cherenkov telescopes, such as MAGIC (Aliu et al. 2008;Aleksic et al. 2011Aleksic et al. , 2012;;Ansoldi et al. 2016;Acciari et al. 2020), VERITAS, and H.E.S.S. (Abdalla et al. 2018), have detected sub-TeV γ-ray emission from a few γ-ray pulsars.The observed TeV spectra can be smoothly joined to those observed by Fermi, and perhaps indicate a transition from curvature radiation (CR) to the inverse-Compton scattering (ICS) regime.It is well known that the high-energy (HE) radiation from a pulsar is directly linked with the structure of its magnetosphere.Therefore, it is necessary to obtain a realistic structure of the magnetosphere to constrain the HE emission production sites and the emission mechanisms.
The analytic vacuum-retarded dipole derived by Deutsch (1955) was widely used as the background field in the standard gap models to study pulsar spectra and light-curve properties.In these models, the CR (Romani 1996;Hirotani & Shibata 2001;Harding et al. 2008;Tang 2008;Pétri 2019) or synchrotron curvature (Torres 2018) from the primary particles was often thought to produce the GeV radiation, and the ICS radiation (Harding et al. 2008;Aleksic et al. 2011;Pétri 2011;Du et al. 2012;Lyutikov et al. 2012;Lyutikov 2013;Mochol & Pétri 2015;Rudak & Dyks 2017) from the primary particles or the secondary pairs is considered as the origin of TeV photons.Actually, a realistic pulsar magnetosphere will deviate from the vacuum state and become a plasma-filled force-free (FF) one (Goldreich & Julian 1969).The aligned FF magnetosphere was first obtained by Contopoulos et al. (1999, CKF) by solving the pulsar equation; they obtained an open−closed magnetospheric configuration with the equatorial current sheet outside the light cylinder (LC).The CKF solution was further confirmed by performing a time-dependent FF simulation (e.g., Komissarov 2006;McKinney 2006;Timokhin 2006;Parfrey et al. 2012;Cao et al. 2016a).The three-dimensional inclined FF magnetosphere was first obtained by Spitkovsky (2006) with the finite-difference time-domain method, which was further explored by Kalapotharakos & Contopoulos (2009) with the inclusion of a perfectly matched layer, and by Pétri (2012) and Cao et al. (2016b) with the pseudospectral method.FF magnetospheres are also used to predict the radiation characteristics of pulsars (e.g., Bai & Spitkovsky 2010;Contopoulos & Kalapotharakos 2010;Harding & Kalapotharakos 2015).
Although the FF magnetosphere was thought to be closer to a realistic pulsar magnetosphere, it was inherently nondissipative, whereas realistic pulsar magnetospheres should have some dissipative regions that Yang & Cao can self-consistently produce the observed pulsar emission.Therefore, the resistive model with a conductivity was developed to produce the dissipative regions.Resistive magnetospheres also became available (e.g., Kalapotharakos et al. 2012a;Li et al. 2012;Cao et al. 2016b) and were used to predict the pulsar γ-ray emission by the test particle trajectory method (Kalapotharakos et al. 2012b(Kalapotharakos et al. , 2014(Kalapotharakos et al. , 2017;;Brambilla et al. 2015;Cao & Yang 2019;Yang & Cao 2021).These studies revealed that the γ-ray radiation mainly comes from the outer magnetosphere near the equatorial current sheet outside the LC.Resistive magnetospheres can produce a self-consistent accelerating electric field by introducing a conductivity, but they cannot explain the microscopic physical origin of the conductivity itself.Therefore, the particle-in-cell (PIC) method was developed to model the pulsar magnetosphere by self-consistently treating the interactions of the particle motions, the emitting photons, and the electromagnetic fields.Some groups attribute the synchrotron radiation (SR, (Cerutti et al. 2016;Philippov & Spitkovsky 2018)) from particles accelerated via magnetic reconnection in the current sheet to the main GeV emission mechanism, while other groups ascribe CR (Brambilla et al. 2018;Kalapotharakos et al. 2018Kalapotharakos et al. , 2022Kalapotharakos et al. , 2023) ) to the main GeV emission mechanism.However, the particle Lorentz factors obtained by the PIC method are much smaller than the real ones that can produce the observed GeV photons.
An alternative Aristotelian electrodynamics (AE), which can include the reaction of the emitting photons on the particle dynamics, has been proposed to model the pulsar magnetosphere (Gruzinov 2012(Gruzinov , 2013)).Moreover, combined FF and AE magnetospheres are developed to simulate the pulsar magnetosphere, which can constrain the dissipative region only near the current sheet (Contopoulos 2016;Cao & Yang 2020;Pétri 2020Pétri , 2022)).Cao & Yang (2022) further calculate the CR energy spectra and light curves from the primary particles in the combined FF and AE magnetospheres by the particle trajectory method.Furthermore, the Fermi energy-dependent double-peak light curves, and the phase-averaged and phase-resolved GeV spectra, can also be reproduced for a high pair multiplicity, indicating that particle acceleration and γ-ray emission mainly originate in the current sheet outside the LC (Cao & Yang 2024).However, the individual Fermi γ-ray band is not yet enough to distinguish between different emission locations and emission mechanisms in the magnetospheres.The study of combined multiwavelength light curves and spectra is thus expected to put a stronger constraint on these emission locations and mechanisms.
Recently, Harding & Kalapotharakos (2015) and Harding et al. (2018Harding et al. ( , 2021) ) explored the multiwavelength radiation properties of pulsars in the FF magnetosphere with a constant accelerating electric field distribution.However, they cannot well reproduce the peak phases of the observed multiwavelength light curves, indicating that a radially and azimuthally dependent accelerating electric field distribution is needed in modeling the pulsar multiwavelength radiation.
In the paper, we further expand the study of Cao & Yang (2022, 2024) by simultaneously computing the multiwavelength radiation properties of pulsars with multiple radiation mechanisms from both the primary particles accelerated by the self-consistent accelerating electric field in the equatorial current sheet and the secondary pairs with an assumed pair spectrum in the combined FF and AE dissipative magnetospheres.The modeled multiwavelength light curves and spectra are then directly compared with those of the Crab, Vela, and Geminga pulsars.In Section 2, we describe the dissipative magnetospheric model and the multiwavelength radiation model.In Section 3, we apply our model to explain the multiwavelength light curves and spectra of the Crab, Vela, and Geminga pulsars.Finally, the conclusions and discussions are presented in Section 4.

The combined FF and AE magnetospheres
The pulsar magnetosphere was obtained by solving the time-dependent Maxwell equations in the comoving frame where B is the magnetic field, E is the electric field, ρ is the charge density, J is the current density, and V rot is the corotating velocity.It is noted that all quantities are defined in the inertial observed frame(IOF).
The FF magnetospheres are thought to be closest to the realistic pulsar magnetosphere, but they cannot include any dissipation to produce the observed multiwavelength radiation.Recently, the combined FF and AE magnetospheres have been developed to approximate a realistic pulsar magnetosphere, which can indirectly include the reaction of the emitting photons on the particle dynamics and introduce the local dissipation where the FF condition is violated.In this paper, we use the combined FF and AE magnetospheres to ap-proximate a realistic pulsar magnetosphere and predict the pulsar multiwavelength radiation.The current density J in the combined FF and AE magnetospheres is defined by the local electromagnetic fields as (Cao & Yang 2020, 2022) where E 0 and B 0 are the electric and magnetic fields in the fluid frame where E is parallel to B. They are related by the following conditions: , where E 0 is the accelerating electric field parallel to the magnetic field.κ is a pair multiplicity phenomenologically connected with the pair cascade process, and is introduced to adjust the strength and the distribution of the accelerating electric field E 0 in the equatorial current sheet.
We simulate a series of the dissipative magnetospheres for the magnetic inclination angles χ from 0 • to 90 • in 5 • intervals, where the FF description is enforced in the E < B regions and the AE description is implemented in the E ≥ B regions.The electric and magnetic fields are discretized in the regions from r ⋆ = 0.2R LC to r max = 3R LC .The pair multiplicities are taken as κ = {0, 1, 3}.Moreover, a high resolution of N r ×N θ ×N ϕ = 129×64×128 is used to better resolve the equatorial current sheet.Lastly, a nonreflecting boundary condition is implemented to avoid inward reflection from the outer boundary so as to ensure the combined FF and AE magnetospheres reach a stationary state for several rotation periods.Figure 1.We plot the assumed normalized pair spectra for Crab, Vela, and Geminga pulsars.The pair spectra of Vela and Geminga pulsars display an exponential power-law shape, while that of the Crab should be presented as a piecewise exponential power law.(2019).The 50 hr HAWC observational sensitivity curve (Abeysekara et al. 2017) is also plotted.

The injecting and tracking of the primary particles and the secondary pairs
It is commonly recognized that the charged particles can be extracted from the pulsar PC surface and rapidly accelerated by the rotating induced electric field to extremely relativistic velocity to produce HE photons, which will be converted in bursts to the secondary pairs through the strong magnetic field in multiple pair cascade processes.This will cause the pulsar magnetosphere to reach the FF state in most regions.The secondary pairs can still produce low-energy SR photons, which may be boosted to HE photons through the ICS or synchrotron self-Compton (SSC) processes.
In our previous work, we only produce the pulsar Fermi γ-ray emission by computing the CR radiation from the primary particles (Cao & Yang 2022, 2024).In this paper, we will further produce the pulsar multiwavelength radiation by computing the SR, CR and ICS spectra from both the primary particles and the secondary pairs.We will inject a series of the primary particles and the secondary pairs from the PC surface in the specified open volume coordinates of (r ovc , ϕ m ) Yang & Cao (Dyks et al. 2004).The radial coordinate r ovc = 1 denotes the PC rim, and r ovc = 0 is the magnetic pole.The azimuthal coordinate ϕ m is the magnetic phase around the PC in the counterclockwise direction, which is divided into 360 equal portions.Primary particles with low initial Lorentz factor(γ 0 = 100) are uniformly injected from the PC surface in the range r ovc = 0.9 − 1 with 0.01 intervals, and will be accelerated by the selfconsistent accelerating electric field along their trajectories.Secondary pairs are injected with an assumed pair spectrum in the range r ovc = 0.8 − 0.9 with 0.01 intervals, but they do not undergo any acceleration along their trajectories.The trajectories of both the primary particles and secondary pairs are tracked from the neutron surface up to 2.5R LC in the IOF by where the two signs correspond to positrons and electrons, which follow different trajectories in the combined FF and AE magnetospheres.
For the primary particles, taking into account the influence of the local accelerating electric field and the radiation losses, their Lorentz factors (γ) along the trajectories are integrated by the following expression Where q e and m e are the charge of the electron and its rest mass, respectively.c is the velocity of light in free space, and R CR = dl/dθ is the curvature radius of the particle trajectories.dl is the segment length within time dt and dθ is the angle spanned by the adjacent velocity vectors along the particle trajectory.We compute the CR and ICS spectra for the primary particles along each trajectory.For the secondary pairs, the exact pair spectra need to be determined by performing a time-dependent pair cascade simulation.However, it is very difficult to perform a three-dimensional pair cascade simulation that includes time-dependent electromagnetic fields and particle production.Only the one-dimensional (1D) pair cascade simulation is available in a time-dependent manner (Timokhin & Arons 2013;Timokhin & Harding 2015).Therefore, we approximate the secondary pair spectra as a power-law exponential cutoff shape with where the γ is the energy of secondary pairs, s is the spectral index, and γ cut is the cutoff energy of the secondary pairs.The 1D pair cascade simulation found that the pair cascade operates only in the supercurrent, J/J GJ > 1, and anticurrent,J/J GJ < 0, regions.Therefore, in addition to injecting the secondary pairs in the range of r ovc , we further limit the pair injection to the supercurrent (J/J GJ > 1) and anticurrent (J/J GJ < 0) regions.For the Crab pulsar, we assume the pair Lorentz factors ranging from γ min = 20 to γ max = 5 × 10 6 .The spectral index s 1 = 2.2 and the cutoff lorentz factor γ cut1 = 7 × 10 4 for γ min ≤ γ < 2 × 10 5 ; s 2 = 3.6 and γ cut2 = 8 × 10 5 for 2 × 10 5 < γ ≤ γ max .For the Vela pulsar, the Lorentz factors range from γ min = 20 to γ max = 6 × 10 5 , s = 2.3 and γ cut = 1 × 10 4 .For the Geminga pulsar, the Lorentz factors range from γ min = 20 to γ max = 6 × 10 5 , s = 2.0 and γ cut = 4 × 10 3 .Their injected pair spectra adopted in the following calculations for the SR and SSC spectra of the secondary pairs are plotted in Figure 1.

The CR spectrum
For the particle with a Lorentz factor γ, the individual CR spectrum at each position along the trajectory will be given by the following expression: where x = E γ /E cur , E γ is the energy of emitting photon, RCR is the characteristic energy of the CR photon, R CR is the local curvature radius at each point along the trajectory, and the function F (x) is defined as where K 5/3 is the modified Bessel function of order 5/3.The total CR spectrum from all the radiation points will be the superposition of the individual particle spectra, which will be weighted by the surface primary particle flux as (Harding et al. 2008;Harding & Kalapotharakos 2015) where n GJ is Goldreich-Julian particle density at the pulsar surface, R ⋆ the pulsar radius, θ pc the pulsar PC half-angle, and ∆A the reduced surface area of the particle injection region relative to the PC area, which can be approximated by ∆A = (r max ovc ) 2 − (r min ovc ) 2 .

The SR spectrum
The individual particle SR spectrum emitted by a particle with a Lorentz factor γ at each emission point along the trajectory is given by the following formula: where B is the magnetic field at each point along the trajectory, α is the pitch angle, x = E γ /E SR , and E SR = 3ℏγ 2 q e Bsinα/2m e c is the critical energy of an SR photon.The SR spectrum from all the secondary pairs will be the integration of the individual particle SR spectra with the normalized secondary particle dis-tribution dN/N 0 dγ, given by the following expression: The total SR spectrum of the secondary pairs should be weighted by the primary Goldreich-Julian particle density of the PC surface with the regulation of the pair multiplicity M pair (the number of pairs produced by each Yang & Cao primary particle), which is used to adjust the pair SR spectra to match the optical to hard-X-ray data and reads as Ṅpair = M pair Ṅp .

The ICS/SSC spectrum
In the paper, we compute the isotropic ICS/SSC of the primary particles and the secondary pairs scattering with the isotropic SR photons of the secondary pairs.The scattering spectrum of an individual particle with a Lorentz factor γ at each emission point along the trajectory is given by Here The total primary particle ICS should be weighted by the primary particle flux Ṅp .The total SSC spectrum from all the secondary pairs will be the integration of the individual particle SSC spectra with the normalized secondary particle distribution and weighted by the same particle flux as described in Section 2.4.

The multiwavelength spectra and light curves
The location and geometry of the acceleration and radiation regions will be strongly imprinted on the pulsar multiwavelength radiation properties of the pulsar.
Therefore, the pulsar multiwavelength light curve and spectra modelings of the pulsar provide a very powerful tool to diagnose the magnetospheric structure, the location of particle acceleration, and the emission mechanisms.In our model, all the particles are injected from the neutron surface -the primary particles from the region of r ovc = 0.9 − 1 and for the secondary pairs from the adjacent region of r ovc = 0.8 − 0.9.Besides, the trajectories of both the primary particles and secondary pairs are separately tracked from the neutron surface up to 2.5R LC .The secondary pairs cannot emit SR unless they acquire some pitch angles.While stringently constrained by the strong magnetic field, the particles will travel almost along it.Kelner et al. (2015) discussed the evolution of the pitch angle in the strong magnetic field from analytical aspects.Lyubarskii & Petrova (1998) and Harding et al. (2008) proposed that the particles can acquire a small pitch angle at some height (r SR ) due to the cyclotron resonance absorption.Therefore, we have choose the SR emission heights to be r SR = 0.7R LC , 0.48R LC , and 0.4R LC for Crab, Vela, and Geminga pulsars, respectively.We first compute the pair SR and collect the pair SR density along the particle trajectories from r SR to 2.5R LC .The pair SSC emissions are then computed by using the collected pair SR density from r SR to 2.5R LC .The primary particles cannot emit the CR and ICS unless their trajectories enter into the equatorial current sheet.Therefore, we track the trajectories of the primary particles from the neutron surface to up to 2.5R LC and only compute their CR and ICS when their trajectories enter into the equatorial current sheet, where the collected pair SR density is used for the target photons for ICS of the primary particles.The similar calculation methods are also used by Harding & Kalapotharakos (2015) and Harding et al. (2021).
Assuming that the directions of the photon emission, η em , are in the direction of the particle motion, β = v/c, for the primary and secondary pairs, we can determine the direction of the photon emission from each particle species using where ζ = acos(µ em ) is the viewing angle.The observed phase is obtained by adding the rotation and time-delay corrections: where the first term, ϕ rot = Ω dt, is the rotation phase, the second term is the phase of the emitting photon, and the third term is the time-delay phase.The multiwavelength sky maps can be produced by collecting all the emitting photons from each particle species in the (ζ-ϕ) plane.The corresponding light curves can be constructed by cutting the sky maps with a given viewing angle ζ.More details about the constructions of the sky maps and the light curves can be found in those previous works (e.g., Kalapotharakos et al. 2014;Harding & Kalapotharakos 2015;Cao & Yang 2019;Yang & Cao 2021).The multiwavelength spectra can also be obtained by collecting all the radiation from both the primary particles and the secondary pairs at a given viewing angle along their trajectories using where d is the distance of the pulsar from the observer.

THE RESULTS
In the section, the dissipative model is used to compute the multiwavelength phase-averaged spectral energy distribution (SED), the sky maps, and the corresponding light curves for Crab, Vela, and Geminga pulsars.Two values of the pitch angle (α = 0.01 and 0.001) are chosen to compare their effects on the SEDs in the following results.The pair cascade multiplicity M pair is taken to be a free parameter to adjust the pair SR spectrum to match the optical to hard-X-ray data.Other parameters that affect the SED are magnetic inclination angle χ, the viewing angle ζ, and the secondary pair spectra.

Yang & Cao
tain the modeling results.In Figure 2, we plot the broadband phase-averaged SED of the Crab pulsar for the magnetic inclination χ = 75 • and viewing angle ζ = 57 • , with pitch angle α = 0.001 and pair multiplicity M pair = 3 × 10 5 , shown as the solid lines in the figure.The viewing angle is chosen by the better match with the Fermi γ-ray data by CR of the primary particles.We also see that the the SR of the secondary pair can match well with the observed spectrum from the optical to the hard-X-ray band if the assumed pair spectra in Figure 1 are used.The combination of the CR spectrum of the primary particles with the SSC spectrum of the pairs can well account for the transition from the GeV to the very-hight-energy spectrum up to around 1 TeV, observed by Fermi and MAGIC.Besides, the primary ICS is sensitive to the low-energy part of the pair SR, and scatters the SR photons in the Thompson limit.Therefore, the predicted ICS spectra from the primary particles produce a peak around 10 TeV that is close to the detection sensitivity of the High Altitude Water Cherenkov Observatory (HAWC).For comparison, the predicted SEDs for a larger pitch angle α = 0.01 are also shown as the dashed lines.When a larger α is used, a lower pair γ cut is required to reproduce the optical to the hard-X-ray spectra.Therefore, the predicted pair SSC spectra for a lower γ cut does not reach TeV energies to well match the MAGIC TeV spectra.Our results suggest that a low α ∼ 10 −3 is needed to reproduce the multiwavelength SEDs of the Crab pulsar.A similar low α is also used to explain the multiwavelength SEDs of the Crab pulsar by Harding et al. (2021) In fact, the inclination angles are not well constrained by observations and multiwavelength radiation modeling, but the viewing angle has been constrained from the measurement of the torus of the Crab wind nebula by Ng & Romani (2008) with ζ ∼ 61 • .It is found that our modeling inclination angle and viewing angle are similar to those of the previous multiwavelength radiation modeling, and our modeling viewing angle of ζ = 57 • is also close to that of the torus of the pulsar wind nebula with ζ ∼ 61 • obtained by Ng & Romani (2008).
In Figure 3, we show the sky maps and the corresponding light curves from the optical to γ-ray bands.We observe that the intensities in the sky maps decrease with increasing energy, and that the bright patterns in the sky maps do not change significantly.This is the caustic effect of the particles whose trajectories can approach the equatorial current sheet and become nearly radial outside the LC.Besides, comparing the modeled light curves with the observed ones, it can be seen that we can systematically well reproduce the observed trends whereby the light curves display two narrow peaks that are approximately aligned in phase for different energy ranges.This arises from the result that the radiation generated by the primary and secondary particles comes from different but adjoint regions near the equatorial current sheet, which leads to the radiation arriving at almost the same phase.Harding et al. (2021) also modeled multiwavelength light curves for the Crab pulsar by using the constant electric distribution in the current sheet of the FF magnetosphere.They found that a constant accelerating electric field distribution can not well reproduce the observed phases of multiwavelength light curves of the Crab pulsar.Our results show that a radially and azimuthally dependent accelerating electric field distributions from the simulations can provide an improved match to the multiwavelength light curves of the Crab pulsar.

Vela pulsar
The real observed parameters, surface magnetic field B = 4 × 10 12 G, period P = 0.089 s, and distance  Cao & Yang (2024).Similarly, the magnetic inclination is not well constrained by the current observation data and multiwavelength radiation modeling, but the viewing angle has been constrained from the measurement of the torus of the Vela wind nebula by (Ng & Romani 2008).We see that the modeled magnetic inclination and viewing angle are similar to those of the previous multiwavelength radiation modeling, and our modeled viewing angle ζ = 63 • is also close to the observed one ζ = 64 • obtained by Ng & Romani (2008).
In Figure 5, we show the sky maps and corresponding light curves from the optical to γ-ray bands.We can see that the bright caustics in the sky maps decrease with increasing energy, and that caustics produced by the secondary pairs distributed in more broader regions than those produced by the primary particles accelerated by the equatorial current sheet.This can be attributed to a significant inner emission from the secondary pair SR within the LC.Besides, the predicted optical light curves Yang & Cao have two narrower peaks than the γ-ray one, which can well match with that observed.The predicted GeV light curves can well match the observed peak separation and relative ratio of the two peaks of the Fermi γ-ray light curves.The observed sub-TeV γ-ray light curves only visible in the second peak can also be well reproduced by our model.It is seen that our results can well reproduce the observed trends of the multiwavelength light curves for the Vela pulsar.Moreover, our results can well reproduce the decreasing ratio of the first γ peak to the second one and the narrowing of the peak width toward the higher energies from the GeV to TeV bands.It is noted that we can give better match to the multiwavelength light curves of the Vela pulsar by including a radially and azimuthally dependent accelerating electric field distribution from the simulations than those obtained by Harding et al. (2021), where they used a uniform accelerating electric field independent of the radius and azimuth.

Geminga pulsar
The real Geminga observed parameters, surface magnetic field B = 1.8×10 12 G, period P = 0.237 s, and distance d = 0.25 kpc for, are used to produce the predicted multiwavelength light curves and spectra.In Figure 6, we plot the broadband SED for the Geminga pulsar with magnetic inclination χ = 70 • and viewing angle ζ = 66 • with α = 0.001 and M pair = 1 × 10 4 in the solid lines.It can be seen that our results well reproduce multiwavelength SEDs of the Geminga pulsar from the optical to γ-ray.We find that the optical data are produced by SR from the pairs, the Fermi data are dominated by the primary CR, and the MAGIC data are an extension of the primary CR.However, the predicted primary ICS spectra are below the current detection sensitivity.For comparison, the predicted SEDs for a larger pitch angle α = 0.01 are also shown as the dashed lines.We find that the pitch angle is not well constrained by observations due to the scarcity of the optical to hard-X-ray data.Therefore, a similar SED is also obtained by using a larger α value for the Geminga pulsar.We note that the magnetic inclination χ and viewing angle ζ were obtained by some previous multiwavelength radiation modeling, with χ = 50 Figure 7 shows the sky maps and the corresponding light curves in the GeV to sub-TeV range.We see that the γ-ray sky maps show the characteristics of the caustic emission from the equatorial current sheet.The predicted GeV light curve is consistent with the one observed by Fermi, and the observed phases and the relative strength can be well reproduced by using the radially and azimuthally dependent accelerating electric field distributions from the simulations.The predicted sub-TeV light curves are only visible in the second peak, and are in agreement with those observed by MAGIC.Moreover, the peaks of the sub-TeV light curves are nearly aligned in phase with those of the GeV light curves, because they are both produced by the primary CR accelerated in the equatorial current sheet.We also find that the ratio of the first γ-ray peak to the second one decreases and the peak width also becomes narrow with the increasing energy, which can be reflected by the decreasing caustics in sky maps.It is seen that our results can well reproduce the observed trends of the multiwavelength light curves for the Geminga pulsar

CONCLUSIONS AND DISCUSSIONS
We have explored the properties of the pulsar multiwavelength light curves and spectra in the combined FF and AE magnetospheres.The combined FF and AE magnetospheres are computed by the pseudo-spectral method in the comoving frame, which can produce the self-consistent accelerating electric field distribution in the equatorial current sheet outside the LC.We use the AE velocity to define the trajectories of the primary particles and secondary pairs in the combined FF and AE magnetospheres.The pulsar multiwavelength light curves and spectra are computed by using multiple radiation mechanisms from both the primary particles accelerated by the accelerating electric field distribution in the equatorial current sheet and an assumed pair spectrum injected from the PC surface.Then, we perform a direct comparison between the predicted multiwavelength light curves and spectra and those of the Crab, Vela, and Geminga pulsars.Our results can well reproduce the observed trends of the multiwavelength light curves and spectra for the these three pulsars.Our results indicate that the high-energy tails of the primary CR spectra for the Vela and Geminga pulsars have been detected by H.E.S.S.-II and MAGIC.The primary ICS component above 1 TeV is detectable for the Crab and Vela pulsars by the current HAWC and future CTA.Moreover, our results can also systematically well re-produce the decreasing ratio of the first γ-ray peak to the second one toward higher energies, especially for the Vela and Geminga pulsars.
The pulsar multiwavelength light curves and spectra are also modeled by using a constant accelerating electric field distribution in the FF magnetosphere by Harding et al. (2021).However, they cannot well reproduce the peak phases of the light curves for those three pulsars.Our model can provide an improved match to the multiwavelength light curves of these three pulsars by using a radially and azimuthally dependent accelerating electric field distribution from the simulations themselves.However, our model cannot explain the bridge emission of the Fermi light curves and the ∼ MeV spectra.It is suggested that these emission components may originate in the inner magnetosphere within the LC (Harding et al. 2021;Barnard et al. 2022).We will explore the observed bridge emission and ∼ MeV spectra by introducing the accelerating electric field dis-tribution not only in the equatorial current sheet outside the LC but also in the inner region within the LC in the future.Moreover, the polarization information can provide another powerful constraint on the locations of the emission regions and on the emission mechanisms.We will use the combined FF and AE magnetospheres to present the study of the pulsar polarization in the near future.

Figure 2 .
Figure 2. The broadband SEDs of the Crab pulsar for the magnetic inclination χ = 75 • , viewing angle ζ = 57 • , and the pair multiple number κ = 3 in the current density of equation (5).The pair multiplicity takes the value Mpair = 3 × 10 5 and the pitch angle α = 0.001 is shown as the solid lines.For comparison, the predicted SEDs for a larger pitch angle α = 0.01 are also shown as the dashed lines.We can see that the lower α produces the well matched SEDs with the multiwavelength SED.The primary CR spectra for these two α are identical because they are independent of the pitch angle.See the context in the results for more details about the discussions of the effect of the pitch angle on the SEDs.Observed data are taken from Kuiper et al. (2001), Abdo et al. (2013), Ansoldi et al. (2016), and Sollerman et al.(2019).The 50 hr HAWC observational sensitivity curve(Abeysekara et al. 2017) is also plotted.

Figure 3 .
Figure3.The broadband sky maps and light curves for the Crab pulsar with the same parameters as in Figure2for α = 0.001.Four energy bands, i.e., (1,10)eV, (0.1,100)GeV, (100,400)GeV, and (0.4,1.5)TeV, accommodating with the observed bands, are used to produce the modeling results.The brightness in sky maps is normalized to the same range, thus they share the same color bar.The horizontal white lines in the sky maps are the viewing angle of ζ = 57 • for generating the modeling light curves in the middle panels.The observed optical light curve is adopted fromSlowikowska et al. (2009).The GeV and sub-TeV light curves are adopted from Fermi(Abdo et al. 2013), and MAGIC (Ansoldi et al. 2016).

Figure 5 .
Figure 5.The broadband sky maps and light curves for the Vela pulsar with the same parameters as in Figure 4. Three energy bands accommodated with the ones observed, i.e., (0.1,10)eV, (0.1,20)GeV, and (20,100)GeV, are used to produce the sky maps and the light curves.The horizontal white line in the left panels is the viewing angle of ζ = 63 • for generating the counterpart light curves in the right panels.The observed light curves at those energies are taken from Harding et al. (2000), Abdo et al. (2013) and Abdalla et al. (2018).

Figure 7 .
Figure 7.The broadband sky maps and light curves for the Geminga pulsar with the same parameters as in Figure 6.Two energy bands accommodated with the corresponding observed ones, i.e., (0.1,15)GeV and (15,75)GeV are used to generate the sky maps and the light curves.The horizontal white line in the sky maps is the viewing angle of ζ = 66 • for generating the counterpart light curves.The observed light curves are taken from Fermi (Abdo et al. 2013) and MAGIC (Acciari et al. 2020).d = 0.29 kpc, are used to produce the multiwavelength light curves and spectra for the Vela pulsar.In Figure 4, we plot the broadband SED for the Vela pulsar with magnetic inclination χ = 60 • and viewing angle ζ = 63• with α = 0.001 and M pair = 2.4 × 10 4 as the solid lines.We see that the flux and the spectral index of the optical data can be well matched by the SR of the secondary pairs, and Fermi data are dominated by the CR of the primary particles accelerated in the equatorial current sheet.The high-energy tail of the primary CR spectra will approach or exceed the H.E.S.S.-II detection sensitivity, indicating that the H.E.S.S.-II has detected the high-energy tail of the primary CR component.However, our model cannot explain the observed emission in the MeV band, which may originate from primary SR radiation in the inner magnetosphere within the LC(Harding & Kalapotharakos 2015;Harding et al. 2021).Moreover, the primary ICS flux is very close to the H.E.S.S.-II detection sensitivity, therefore it is possible that the primary ICS component may partly explain observed TeV emission with the accumulated H.E.S.S.-II data in the future.For comparison, the predicted SEDs for a larger pitch angle α = 0.01 are also shown as the dashed lines.We find that the pitch angle is not well constrained by observations due to the scarcity of the optical to hard-X-ray data.Therefore, a similar SED is also obtained by using a larger α value for the Vela pulsar.The magnetic inclination angle χ and viewing • and ζ = 86 • obtained by Zhang & Cheng (2001), χ = 60 • and ζ = 90 • by Pétri (2009), χ = 45 • and ζ = 87 • by Brambilla et al. (2015), χ = 75 • and ζ = 55 • byHarding et al. (2021).It is found that the magnetic inclination angle and viewing angle were not well constrained; only the viewing angle is estimated as ζ > 60 • from the measurement of the Geminga X-ray image(Caraveo et al. 2003;Pierbattista et al. 2015).Although the magnetic inclination angles vary among different radiation models, our modeled viewing angle of ζ = 66 • lies well in the range of ζ > 60 • .