The Modeling of Pulsar Magnetosphere and Radiation

Abstract

We review the recent advances in the pulsar high-energy $\gamma$-ray observation and the electrodynamics of the pulsar magnetospheres from the early vacuum model to the recent plasma-filled models by numerical simulations. The numerical simulations have made significant progress toward the self-consistent modeling of the plasma-filled magnetosphere by including the particle acceleration and radiation. The current numerical simulations confirm a near force-free magnetosphere with the particle acceleration in the separatrix near the light cylinder and the current sheet outside the light cylinder, which can provide a good match to the recent high-energy $\gamma$-ray observations. The modeling of the combined multi-wavelength light curves, spectra, and polarization are expected to provide a stronger constrain on the geometry of the magnetic field lines, the location of the particle acceleration and the emission region, and the emission mechanism in the pulsar magnetospheres.

1. Introduction

Neutron stars are very fascinating astrophysical objects and excellent laboratories for studying the fundamental physics of the strong electromagnetic and gravitational fields. Pulsars are identified as rapidly rotating neutron stars with very strong surface magnetic fields (${10}^8$–${10}^{12}$ G) and rotation periods from milliseconds to seconds. It is believed that a pulsar will lose a substantial fraction of its rotational kinetic energy and convert it into particle acceleration and radiation. The periodic electromagnetic signals from these objects are visible across the entire electromagnetic spectrum from the radio to $\gamma$-ray bands. The multi-wavelength spectra from these objects are thought to be produced by the synchrotron, curvature, and inverse-Compton radiation from the high-energy particle accelerated by unscreened electric fields in the magnetosphere. However, it is still unclear where the particles are accelerated in the magnetosphere and how the subsequent radiation is produced. Although the pulsars were first discovered in the radio band [1], it is difficult to understand the origin of the radio emission due to the coherent nature in the radio band. Since the emission at the higher energies originates from the incoherent process, the origin of the pulsar emission and the particle acceleration in the magnetosphere can better be understood at the higher energies (for example, in the $\gamma$-ray band). The large area telescope (LAT) of the Fermi Gamma-Ray Space Telescope and ground-based Cherenkov telescopes have opened a new era in the study of the pulsar $\gamma$-ray physics. Fermi-LAT observations have discovered a large number of the $\gamma$-ray pulsars with the measurement of good-quality light curves and spectra in the GeV $\gamma$-ray bands [2,3,4]. Ground-based Cherenkov telescopes have also detected the sub-TeV $\gamma$-ray emission for a growing number of the Fermi pulsars [5,6,7,8]. The location and geometry of the particle acceleration and emission regions will strongly imprint on the multi-wavelength light curves and the phase-resolved spectra. The observed multi-wavelength light curves and spectra can provide the powerful probes of the emission region geometry and emission mechanics. The properties of the pulsar multi-wavelength emission are directly associated with the structure of the pulsar magnetospheres, which require a better and deeper understanding of the realistic pulsar magnetosphere.

In the pulsar magnetosphere, the electromagnetic field, the particle dynamics and the subsequent radiation mechanics are all coupled mutually. This requires us to perform a self-consistent calculation of the Maxwell equations, the particle dynamics, and the subsequent radiation to model the realistic pulsar magnetospheres. Significant progresses have been made toward the self-consistent modeling of the pulsar magnetosphere over the last few decades. The first solution of the pulsar magnetosphere is the vacuum fields, which has the advantage of the analytical expression derived by Deutsch (1955) [9]. However, the vacuum solution is not a real pulsar model, since it does not take into account the modification of the current and charge on the magnetosphere structure. It is well established that the pulsar magnetosphere would be filled with plasmas created by the pair cascades. A zeroth-order approximation of the plasma-filled magnetosphere is referred to as the force-free electrodynamics, where enough particles are produced to screen all of the accelerating electric field so that force-free condition $\mathit{E}·\mathit{B}=0$ holds everywhere in the magnetosphere. The force solutions have recently became available by the numerical simulation for the aligned rotator [10,11,12,13,14,15] and for the oblique rotator [16,17,18,19]. Although the force-free solutions are thought to be closest to realistic pulsar magnetospheres, they do not allow any particle acceleration and the production of radiation in the magnetospheres. The resistive pulsar magnetosphere is then developed by introducing the self-consistent accelerating electric field with a prescribed conductivity [20,21,22]. The resistive magnetospheres are still not self-consistent, since they cannot provide microscopic physical information that produces the magnetospheric conductivity. The kinetic particle-in-cell (PIC) models are later developed to self-consistently compute the feedback between the particle motions, the radiating photons and the electromagnetic fields [23,24,25,26,27,28,29]. The PIC models are not yet fully self-consistent, since they cannot use the real pulsar parameters to model the pulsar magnetosphere and predict the pulsar emission due to the large-scale separation between the plasma frequency and the pulsar rotation frequency. A clever method with the combined force-free and Aristotelian electrodynamics (AE) was recently developed to model the pulsar magnetosphere [30,31,32,33], which can include the back-reaction of the emitting photons onto particle dynamics and introduce the local accelerating electric field to produce the observed pulsar emission with the real pulsar parameters. All these magnetospheric simulations obtain a near force-free magnetosphere with the particle acceleration near the current sheet outside the light cylinder (LC). The good matches between the magnetospheric simulations and the $\gamma$-ray observations indicate that the particle acceleration and the $\gamma$-ray emission mainly originate from the regions near the current sheet outside the LC. However, there are some disagreements over whether the synchrotron or curvature radiation from the current sheet is the main Fermi $\gamma$-ray emission mechanism. The light curve and spectral information are not sufficient to distinguish between different emission region geometry and emission mechanisms in the pulsar magnetosphere. The pulsar polarization information can provide another independent constraint on the emission region geometry and emission mechanisms. The modeling of the combined multi-wavelength emission and polarization is expected to provide a stronger constrain on the geometry of the magnetic field lines, the location of the particle acceleration and the emission region, and the emission mechanism in the pulsar magnetospheres.

In this review, we summarize the recent advances in the pulsar high-energy $\gamma$-ray observations and the electrodynamics of the pulsar magnetospheres from the early vacuum model to recent plasma-filled models by the numerical simulations. We focus on the numerical simulations of recent plasma-filled magnetospheres and the fruitful feedback between the $\gamma$-ray observations and the magnetospheric simulations from the force-free, resistive, PIC and combined force-free and AE models. A similar excellent review was recently presented by Philippov and Kramer (2022) [34], where they focus on the PIC simulations of the pulsar magnetospheres and the meaningful connections between the radio and $\gamma$-ray observations and the PIC simulations. More recent reviews can be found in [35,36,37,38,39].

2. Pulsar Gamma-Ray Observations

The Fermi Gamma-Ray Space Telescope launched in 2008 provides an unprecedented opportunity to study the physics of the pulsar $\gamma$-ray emission. Three Fermi Pulsar Catalogs have been published in [2,3,4]. In fact, more than 40 $\gamma$-ray pulsars were detected in the first year of the Fermi large area telescope (LAT) operations [2]. In total, 117 $\gamma$-ray pulsars are listed in the Second Fermi Pulsar Catalog [3]. To date, Fermi-LAT have discovered 294 confirmed $\gamma$-ray pulsars, including 150 young gamma-ray pulsars and 144 millisecond pulsars (MSPs) [4]. The detected $\gamma$-ray pulsars can be divided into young radio-loud, young radio-quiet, and MSPs. Fermi observations have demonstrated that the pulsars are the dominant source of the GeV $\gamma$-rays in the Milky Way. A wealth of high-quality light curves, phase-averaged and phase-resolved spectra are now available. In the third Fermi Pulsar Catalog, the light curves show a wide variety of the pulsed profiles: single-peak, double-peak, and multi-peak profiles, where their ratios are ∼0.26, ∼0.53, and ∼0.21, separately. Moreover, in 120 Fermi pulsars with the double-peak profiles, the distribution of the phase separations ($\mathsf{\Delta }$) is not uniform, varying from ∼0.082 to ∼0.73, where 55% Fermi pulsars have $\mathsf{\Delta }\ge 0.4$. On the other hand, the radio peaks of young Fermi pulsars typically lead the first $\gamma$-ray peak by a small fraction of rotation period, but it is not the case for Fermi MSPs. The double-peaked $\gamma$-ray light curves usually show an energy evolution where the relative ratio of the first to the second $\gamma$-ray peak decreases toward the higher energies (see Figure 1). The observed $\gamma$-ray spectra are usually characterized by a power law with an exponential cutoff, and the cutoff energy typically lies in the range of 1–5 GeV. The large sample of the Fermi $\gamma$-ray pulsars also find some statistical trends and correlations among their observational quantities, typically including the radio lag and peak separation correlation, as well as the dependence of the $\gamma$-ray luminosity on spin-down power.

The detections of the pulsed TeV emission from the ground-based Cherenkov telescopes can provide some new clues on the particle acceleration and emission mechanics in the pulsar magnetosphere. The ground-based Cherenkov telescopes have detected the pulsed $\gamma$-ray emission at the sub-TeV energies for a growing number of Fermi pulsars. These pulsars include the Crab, Vela, Geminga, and PSR B1706-44 with the first two having pulsed emission above 1 TeV. MAGIC have detected the pulsed emission to up to 1.5 TeV from the Crab pulsar [5]. H.E.S.S-II have reported the detection of the pulsed emission from 20 to 100 GeV and also above 3 TeV from the Vela pulsar [6], and in the energy range below 100 GeV from PSR B1706-44 [7]. MAGIC also have announced the detection of pulsed emission between 15 and 75 GeV with the significance high up to 6.3 $\sigma$ from the Geminga pulsar [8]. The observed light curves show only the second $\gamma$-ray peak, which confirms the energy evolution of the Fermi light curves with the decreasing ratio of the first to second $\gamma$-ray peaks toward the higher energies. The observed $\gamma$-ray spectrum is smoothly connected to the Fermi spectrum, which shows a power-law extension of the Fermi spectrum for the Vela, Geminga, and PSR B1706-44 pulsars and a broken power-law extension of the Fermi spectrum for the Crab pulsar (see Figure 2). The detected TeV photon energy directly measures the maximum energy of the accelerated particles, which can provide an excellent diagnostic of the particle acceleration and emission mechanics in the pulsar magnetospheres.

3. Time-Dependent Maxwell Equations

The pulsar magnetosphere can be obtained by numerically solving the time-dependent Maxwell equations

$$\begin{array}{c}\hfill \frac{\mathsf{\partial }\mathit{B}}{\mathsf{\partial }t}=-\nabla \times \mathit{E},\\ \hfill \frac{\mathsf{\partial }\mathit{E}}{\mathsf{\partial }t}=\nabla \times \mathit{B}-\mathit{J},\\ \hfill \nabla ·\mathit{B}=0,\\ \hfill \nabla ·\mathit{E}={\rho }_{\mathrm{e}},\end{array}$$

(1)

where $\mathit{J}$ is the current density and ${\rho }_{\mathrm{e}}$ is the charge density. The time-dependent Maxwell equations can be numerically solved by specifying the current density $\mathit{J}$, which can be defined as a function of the fields ($\mathit{E}$, $\mathit{B}$) in the MHD models or be obtained by solving the Vlasov equation in the PIC models.

The magnetospheric simulation is set up with an initialized dipole magnetic field in vacuum, whose magnetic moment $\mu$ is inclined at an angle $\alpha$ with respect to the rotation axis. The neutron star is well approximated by a perfect conductor; the inner boundary condition at the stellar surface is enforced with a rotating electric field

$$\begin{array}{c}\hfill \mathit{E}=-(\Omega \times \mathbf{r})\times \mathit{B}.\end{array}$$

(2)

A non-reflecting or absorbing outer boundary condition is imposed to prevent artificial reflection from the outer boundary. The pulsar magnetosphere can be obtained by relaxing the time-dependent Maxwell equations toward a stationary state.

4. The Vacuum Dipole Magnetospheres

4.1. Field Structure of the Vacuum Magnetospheres

The simplest model of the pulsar magnetosphere is the vacuum dipole magnetosphere. It is assumed that the neutron star is surrounded by vacuum. There are no plasma or particles outside the neutron star, which implies ${\rho }_{\mathrm{e}}=0$ and $\mathit{J}=0$ in the time-dependent Maxwell equations. Therefore, the time-dependent Maxwell Equation (1) is linear in the vacuum case, which can be analytically solved by using the inner boundary condition (2) at the stellar surface. Deutsch (1955) [9] first derived the exact analytical expression for a vacuum dipole magnetosphere with a finite stellar radius. The general-relativistic extension of the Deutsch vacuum solution including the space curvature and frame-dragging effects is also derived by [41]. The exact analytical expression for the Deutsch vacuum solution is given by [9,18,42].

$$\begin{array}{cc}\hfill B_r(\mathit{r},t)& =2B\left[\frac{R^3}{r^3}cos\alpha cos\theta +\frac{R}{r}\frac{h_1^{\left(1\right)}\left(kr\right)}{h_1^{\left(1\right)}\left(kR\right)}sin\alpha sin\theta {\mathrm{e}}^{i\psi }\right],\hfill \\ \hfill B_{\theta }(\mathit{r},t)& =B\left[\frac{R^3}{r^3}cos\alpha sin\theta +\left(\frac{R}{r}\frac{\frac{d}{dr}\left(r{h}_1^{\left(1\right)}\left(kr\right)\right)}{h_1^{\left(1\right)}\left(kR\right)}+\frac{R^2}{r_{\mathrm{L}}^2}\frac{h_2^{\left(1\right)}\left(kr\right)}{\frac{d}{dr}\left(r{h}_2^{\left(1\right)}\left(kr\right)\right){\mid }_R}\right)sin\alpha cos\theta {\mathrm{e}}^{i\psi }\right],\hfill \\ \hfill B_{\varphi }(\mathit{r},t)& =B\left[\frac{R}{r}\frac{\frac{d}{dr}\left(r{h}_1^{\left(1\right)}\left(kr\right)\right)}{h_1^{\left(1\right)}\left(kR\right)}+\frac{R^2}{r_{\mathrm{L}}^2}\frac{h_2^{\left(1\right)}\left(kr\right)}{\frac{d}{dr}\left(r{h}_2^{\left(1\right)}\left(kr\right)\right){\mid }_R}cos2\theta \right]isin\alpha {\mathrm{e}}^{i\psi }.\hfill \\ \hfill E_r(\mathit{r},t)& =\frac{\mathsf{\Omega }BR}{c}\left[-\frac{R^4}{r^4}(3{cos}^2\theta -1)cos\alpha +3sin\alpha sin2\theta {\mathrm{e}}^{i\psi }\frac{R}{r}\frac{h_2^{\left(1\right)}\left(kr\right)}{\frac{d}{dr}\left(r{h}_2^{\left(1\right)}\left(kr\right)\right){\mid }_R}\right],\hfill \\ \hfill E_{\theta }(\mathit{r},t)& =\frac{\mathsf{\Omega }BR}{c}\left[sin\alpha {\mathrm{e}}^{i\psi }\left(\frac{R}{r}\frac{\frac{d}{dr}\left(r{h}_2^{\left(1\right)}\left(kr\right)\right)}{\frac{d}{dr}\left(r{h}_2^{\left(1\right)}\left(kr\right)\right){\mid }_R}cos2\theta -\frac{h_1^{\left(1\right)}\left(kr\right)}{h_1^{\left(1\right)}\left(kR\right)}\right)-\frac{R^4}{r^4}sin2\theta cos\alpha \right],\hfill \\ \hfill E_{\varphi }(\mathit{r},t)& =\frac{\mathsf{\Omega }BR}{c}\left[\frac{R}{r}\frac{\frac{d}{dr}\left(r{h}_2^{\left(1\right)}\left(kr\right)\right)}{\frac{d}{dr}\left(r{h}_2^{\left(1\right)}\left(kr\right)\right){\mid }_R}-\frac{h_1^{\left(1\right)}\left(kr\right)}{h_1^{\left(1\right)}\left(kR\right)}\right]isin\alpha cos\theta {\mathrm{e}}^{i\psi },\hfill \end{array}$$

(3)

where B, $\mathsf{\Omega }$ and R are the surface magnetic field, the rotation speed, and the stellar radius, respectively, $h_l^{\left(1\right)}$ is the spherical Hankel functions of order l, $\psi =\varphi -\mathsf{\Omega }t$ is the instantaneous phase at time t, $k=1/r_{\mathrm{L}}$ is the wavenumber, and $r_{\mathrm{L}}$=c/$\mathsf{\Omega }$ is the light cylinder. The physical solution is obtained by taking the real parts of each component.

The Deutsch vacuum solution can be used as a benchmark to investigate a plasma-filled pulsar magnetosphere by including the plasma feedback onto the vacuum fields. It can also serve as a reference solution for checking the algorithm and accuracy of the numerical solution. A spectral method was developed by Pétri (2012a) [18] and Cao et al. (2016a) [15] to compute the structure of pulsar magnetosphere; they tested the spectral algorithm by comparing the numerical solution and the Deutsch vacuum solution. Figure 3 shows the magnetic field lines of the perpendicular vacuum rotator in the equatorial plane for the numerical solution and analytical solutions. It is shown that the spectral method can accurately compute vacuum electromagnetic fields. The Deutsch vacuum solution is very close to the static dipole near the stellar surface, but the sweepback of the field lines produces a significant toroidal magnetic field by the displacement currents outside the LC. The rotating magnetic field induces a very strong quadrupole electric field near the stellar surface, which falls off as $1/r^4$ much faster than the magnetic field $1/r^3$. The sweepback of magnetic field lines distorts the polar cap and shifts it toward the trailing side. It is well believed that the pulsars lose their rotational kinetic energy by particle acceleration and radiation, leading to the spin down of the pulsar. A rotating vacuum dipole radiates an outward electromagnetic wave by magnetic dipole radiation, which is proposed to explain the luminosity of the Crab nebula [43]. The Poynting flux for a rotating vacuum dipole is given by

$$\begin{array}{c}\hfill L_{\mathrm{vac}}=\frac{2}{3}\frac{{\mu }^2{\mathsf{\Omega }}^4}{c^3}{\sin }^2\alpha .\end{array}$$

(4)

It is found that the vacuum-aligned rotator is not an active pulsar magnetosphere, $L_{\mathrm{alinged},\mathrm{vac}}=0$, which does not dissipate the Poynting flux and hence does not spin down.

4.2. The Radiation Models from the Vacuum Magnetospheres

The vacuum dipole field is the first solution of the global pulsar magnetosphere, which has the advantage of the analytical expression given by Equation (3). Therefore, it is widely used as the background magnetic field of the pulsar radiation modeling in the early study of the pulsar emission, even though Goldreich and Julian (1969) [44] pointed out that the realistic pulsar magnetosphere cannot be surrounded by vacuum. Some acceleration gap models are developed to explain the observed pulsar emission based on the vacuum dipole field. These gap models usually use the vacuum dipole magnetosphere as an approximation to the force-free magnetosphere. A local deviation with the force-free condition forms an acceleration gap with an unscreened accelerating electric field, which accelerates the particles and produces the observed pulsed emission. The accelerating electric field, ${\mathit{E}}_{\parallel }$, is found by solving Poisson’s equation

$$\begin{array}{c}\hfill \nabla ·{\mathit{E}}_{\parallel }=4\pi (\rho -{\rho }_{\mathrm{GJ}}),\end{array}$$

(5)

where $\rho$ is the charge density and ${\rho }_{\mathrm{GJ}}$ is the Goldreich–Julian charge density, which is the minimum density required to locally screen the accelerating electric field. However, the accelerating electric fields in these gap models do not self-consistently come from a global magnetosphere model.

The acceleration and radiation of the particles in the gap models are generally assumed to occur inside the LC. Such gap models have been developed to model the pulsar emission phenomena, including the polar gap (PC) [45,46], slot gap (SG, or the two-pole caustic model, TPC) [47,48] and outer gap (OG) [49,50,51] models. The PC model locates the particle acceleration at about several stellar radius above the PC surface. It is assumed that the particles are freely extracted from the neutron star surface, which leaves a charge deficit to produce an ${\mathit{E}}_{\parallel }$ above the PC surface. These particles are accelerated by the ${\mathit{E}}_{\parallel }$ to radiate the $\gamma$-ray photons by the curvature radiation. The produced $\gamma$-ray photons induce a pair cascade by one-photon magnetic ($\gamma$-B) pair production in the strong magnetic field. The pair cascades produce a low-altitude pair formation front, above which all ${\mathit{E}}_{\parallel }$ are screened by the pairs. If the pair formation front is extended to the high latitude, an SG gap is formed between the last closed lines and the pair formation front. These particles are not accelerated to enough high energy to produce the pairs at the low altitude, because the ${\mathit{E}}_{\parallel }$ near the last open field line is weak. However, they can be continuously accelerated to high altitude and produce the $\gamma$-ray photons by the curvature or inverse-Compton radiation. The OG model is a vacuum gap between the null-charge surfaces and the last closed field lines, where ${\mathit{E}}_{\parallel }$ are developed because of charge depletion. This gap accelerates the particles to produce the $\gamma$-ray photons by the curvature radiation. These $\gamma$-ray photons can interact with the surface X-ray photons to initiate a pair cascade by the photon–photon ($\gamma$-$\gamma$) interaction. The size of the gap is limited by the screening of the ${\mathit{E}}_{\parallel }$ by the pair cascades. The OG electrodynamical simulations also show that the inner boundary of the outer gap can shift toward the stellar surface from the null-charge surface when the current through the gap is taken into account [52,53,54]. The current sheet or the stripped wind (SW) model is also proposed as a potential site of the particle acceleration in addition to the traditional gap regions [55,56,57,58,59,60], where the particle acceleration and high-energy radiation are located in the current sheet near or outside the LC. An annular gap model is also developed to model the pulsar emission [61,62], which locates the emission region between the critical and last open field lines. The schematic diagram of different radiation models is depicted in Figure 4.

4.3. Light Curve and Spectra Modeling from the Vacuum Magnetospheres

The Fermi $\gamma$-ray observations have accumulated a wealth of the high-quality pulsar light curves over the last few decades. A geometric light curve model is usually used to fit the pulsar Fermi light curves based on the PC, SG and OG models with the vacuum dipole field, which can help us test different emission models and constrain the geometry of different emission regions. The geometric light curve can be produced by assuming the uniform emissivity along the field lines in the comoving frame. The direction of the photon emission can be obtained by a revised aberration formula given by [63]. The bright caustic is a common feature of the 3D emission pattern, which is formed when the phase delays from the field line sweepback cancel those due to the aberration and time-delay effect. The emission from different field lines arrives to an observer in the same phase. One or two sharp peaks in the light curves are seen when the viewing angle of the observer crosses the caustic. Figure 5 shows the sky maps and light curves from the PC, SG, and OG model at different inclination angles and viewing angles. The PC model predicts a hollow-cone emission pattern with the peak at phases 0 and 0.5, which produces a near-phase alignment of the radio and gamma-ray profiles. These models have difficulty producing the widely separated $\gamma$-ray light curves seen by Fermi-LAT. The original OG model can produce the double-peak light curve profiles with little off-peak emission from only one magnetic pole, because the gap forms above the null-charge surface where $\zeta ={90}^{\circ }$ (but the revised version of the OG model can solve the problem (e.g., [64])). The SG model also can produce the double-peak light curve profiles but from the opposite poles. Compared with the original OG model, the SG model has the emission below and above the null-charge surface, which can produce a significant off-peak component seen in some young radio-faint pulsars and MSPs.

The geometric light curve models based on different emission models were used to fit a sample of the $\gamma$-ray pulsar light curves in the 2PC [65,66]. It is found that the OG and SG models can generally fit the observed $\gamma$-ray light curves. These studies can constrain the magnetosphere geometry of each pulsar described by magnetic inclination angles $\alpha$ and viewing angles $\zeta$. More constraints can be obtained by jointly fitting the pulsar radio and $\gamma$-ray light curves [67,68,69,70,71]. The geometric light curve model cannot give the information of the pulsar energetics and spectra. The modeling of pulsar energy-dependent light curves and spectra can provide a better constraint on the modeling parameters and radiation mechanisms. The PC model predicts a super-exponential spectral shape due to magnetic pair production near the neutron star surface. The measurement of the exponential cutoff spectra by Fermi observation directly ruled out the super-exponential spectral shape predicted by the PC models. Therefore, the SG and OG models are more favored, since these models predict an exponential cutoff curvature spectrum at several GeV energies in agreement with the Fermi observation. The energy-dependent $\gamma$-ray light curves and spectra of the Vela pulsar were modeled by the curvature radiation for the OG [72] and annular gap model [73]. The multi-wavelength radiation of the Crab pulsar was also modeled by the OG [74], SG [75] and annular gap [76] models. Moreover, the pulsar polarization properties were explored by the OG [77] and SG models [78]. These models generally predict a fast swing of the position angle (PA) and the depolarization dips during the pulse peaks, which are similar to the measured polarization of the Crab optical emission. It is noted that the pulsar emission cannot be modeled beyond the LC in these gap models, since the particle trajectory does not follow the magnetic line beyond the LC where the particle velocity exceeds the speed of light.

5. The Force-Free Magnetospheres

5.1. Field Structure of the Force-Free Magnetospheres

After the discovery of the Deutsch vacuum solution, Goldreich and Julian (1969) [44] realized that the Deutsch vacuum solution produces an ${\mathit{E}}_{\parallel }$ accelerating electric field along the magnetic line at the stellar surface. This accelerating electric field lifts the particles off the star surface to fill the magnetosphere. On the other hand, these lifted particles are accelerated by the ${\mathit{E}}_{\parallel }$ to radiate the $\gamma$-ray photons by the curvature radiation and inverse-Compton scattering of the surface thermal X-ray photons, which initiates a pair cascade by the $\gamma$-B interaction near the stellar surface and $\gamma$-$\gamma$ interaction near the LC. The pair cascades produce more secondary pairs to fill the magnetosphere. If the plasma density is much higher than the GJ density, all accelerating electric fields are shorted out so that the force-free condition $\mathit{E}·\mathit{B}=0$ holds in the magnetosphere everywhere. This corresponds to the zeroth-order approximation of the plasma-filled magnetosphere, which is referred to as force-free electrodynamics. Observations of pulsar wind nebulae also indicate that a large amount of plasma is presented in the magnetosphere [79]. Therefore, a realistic pulsar magnetosphere should be closer to the force-free field than the vacuum dipole field.

The force-free approximations assume the negligible plasma inertia and pressure; the Lorentz force acting on the plasma fluid element will vanish,

$$\begin{array}{c}\hfill {\mathit{\rho }}_{\mathit{e}}\mathit{E}+\mathit{J}\times \mathit{B}=0.\end{array}$$

(6)

By taking the cross-product to Equation (6) with $\mathit{B}$, we have

$$\begin{array}{c}\hfill \mathit{J}={\mathit{\rho }}_{\mathit{e}}\frac{\mathit{E}\times \mathit{B}}{B^2}+\frac{(\mathit{J}·\mathit{B})\mathit{B}}{B^2}=0.\end{array}$$

(7)

where $B^2\ne 0$ has been assumed. By using the force-free condition $\mathit{E}·\mathit{B}=0$, we have

$$\begin{array}{c}\hfill \frac{\mathsf{\partial }}{\mathsf{\partial }t}(\mathit{E}·\mathit{B})=\frac{\mathsf{\partial }\mathit{E}}{\mathsf{\partial }t}·\mathit{B}+\mathit{E}·\frac{\mathsf{\partial }\mathit{B}}{\mathsf{\partial }t}=0.\end{array}$$

(8)

By substituting the time-derivative terms of $\mathit{E}$ and $\mathit{B}$ with the time-dependent Maxwell Equation (1), Equation (8) becomes

$$\begin{array}{c}\hfill (\nabla \times \mathbf{B})·\mathit{B}-\mathit{J}·\mathit{B}-\mathit{E}·(\nabla \times \mathbf{B})=0.\end{array}$$

(9)

By substituting the $\mathit{J}·\mathit{B}$ term in Equation (7) into Equation (9), the force-free current density is given by [80,81,82].

$$\begin{array}{c}\hfill \mathit{J}={\mathit{\rho }}_{\mathit{E}}\frac{\mathit{E}\times \mathit{B}}{B^2}+\frac{(\mathit{B}·\nabla \times \mathbf{B}+\mathit{E}·\nabla \times \mathbf{B})\mathit{B}}{B^2}.\end{array}$$

(10)

The time-dependent force-free Equations (1) and (10) can be reduced to a one-dimensional pulsar equation for an axisymmetric aligned rotator [83]. An exact analytical solution to the pulsar equation was found for a rotating magnetic monopole [84], which gives the far-field solution with the asymptotically radial field lines. However, no analytical solution is found for a rotating dipole. It was realized that the force-free pulsar magnetospheres cannot be analytically solved due to the problem of nonlinearity and singularity. This renews new interest in the numerical simulation of the pulsar magnetosphere thanks to the development of numerical techniques and increasing computer power. The numerical solution to the pulsar equation was first obtained by Contopoulos et al. (1999, hereafter CKF) [10], who used a clever iterative algorithm to determine the poloidal current distribution that allows the magnetic field lines to smoothly cross the LC. The CKF solution consists of the closed and field line region separated by a Y-shaped current sheet outside the the light cylinder. It was later found that steady-state solutions can be constructed with the Y-point inside the light cylinder [85,86]. Therefore, it is compulsory for performing the time-dependent force-free simulation to solve the problem of the uniqueness and stability of the CKF solution. The time-dependent force-free equations for the oblique rotator were first solved by Spitkovsky (2006) [16] through the finite-difference time-domain (FDTD) method in the Cartesian coordinate and followed by Kalapotharakos and Contopoulos (2009) [17] with the incorporation of the absorbing outer boundary in the FDTD method. A 3D spectral method in the spherical coordinate was also developed by Pétri (2012) [18] and Cao et al. (2016b) [22] to solve the time-dependent force-free equations, which can impose the exact boundary condition by the use of the spherical coordinate. Moreover, the force-free electrodynamic was also extended to the full magnetohydrodynamic regime that takes plasma inertial and pressures into account [11,87] and the general-relativistic force-free regime that takes the frame-dragging and time–space curvature effects into account [88,89,90]. All these time-dependent force-free simulations confirmed the closed–open CKF solution separated by a separatrix near the LC and a current sheet outside the the LC.

Figure 6 shows the distributions of the force-free magnetic field line and the current density parallel to the magnetic field for different magnetic inclination angles. The presence of charges and currents modifies the field structure of the pulsar magnetosphere. This effect produces a toroidal magnetic field, which increases a fraction of the open field lines across the LC. The magnetic field is close to the vacuum dipole with a subdominant toroidal magnetic field component inside the LC. The magnetic field becomes gradually dominated by the toroidal component and becomes asymptotically monopolar outside the LC. The extra toroidal component from the current increases the sweepback of the field lines compared to the vacuum dipole. A generic feature of the force-free magnetosphere is the formation of the equatorial current sheet outside the LC, where the magnetic field changes polarity and decreases to zero. For the aligned rotator, the current sheet is axisymmetric along the rotational equatorial plane outside the LC. For the oblique rotator, the current sheet becomes asymmetric with the undulating shape of the striped wind, which oscillates about the rotational equatorial plane with the angular amplitude $2\alpha$ and a wavelength of $2\pi r_{\mathrm{L}}$.

Figure 7 shows the distributions of the current density parallel to the magnetic field on the force-free polar caps (the solid lines) for different magnetic inclination angles. The current flows out of one region of the polar cap and returns to another region of the polar cap. For small $\alpha$, the current flows out along the open field lines and returns to the star surface by the Y-shaped current sheet. As $\alpha$ increases, a gradually smaller fraction of the return current reaches the stellar surface by the Y-shaped current sheet. For $\alpha ={90}^{\circ }$, all the current reaching the stellar surface is distributed over half of the polar caps, implying that the equatorial current sheet outside the LC is not directly connected to the star surface. It is also shown that the force-free field has a larger PC size than that of the vacuum fields because of the larger open field volume. Moreover, the force-free PC rim is more shifted toward the trailing side due to the higher field line sweepback relative to the vacuum field. The Poynting flux for a force-free rotator can be well approximated by

$$\begin{array}{c}\hfill L_{\mathrm{FF}}=\frac{{\mu }^2{\mathsf{\Omega }}^4}{c^3}(1+{\sin }^2\alpha ).\end{array}$$

(11)

It is found that the force-free aligned rotator radiates with the Poynting flux of $L_{\mathrm{aligned},\mathrm{FF}}=\frac{{\mu }^2{\mathsf{\Omega }}^4}{c^3}$. This contrasts with the vacuum aligned rotator that does not radiate.

5.2. Light Curve and Spectra Modeling from the Force-Free Magnetospheres

The force-free solution provides more realistic field structures of the pulsar magnetosphere than that of the vacuum solution, which can be used to model the pulsar light curves and spectra based on different radiation models with the assumed ${\mathit{E}}_{\parallel }$ distributions. However, these models are not self-consistent since the force-free solutions do not allow any particle acceleration along magnetic field lines with ${\mathit{E}}_{\parallel }=0$. The force-free field structures were used to model the geometric $\gamma$-ray light curves by assuming the uniform emissivity along the particle trajectory in the observer frame [91]. The particle trajectory is defined as the combination of a drift velocity and one along the field line, which can help us compute the pulsar emission from the stellar surface to beyond the light cylinder. The projections of particle trajectories on the sky map were computed and the geometric $\gamma$-ray light curves were produced based on the ‘separatrix’ emission model extending from the stellar surface to beyond the LC [91]. The “separatrix” model forms the sky map stagnation (SMS) caustics, because the backward drift of the particle trajectory along the asymptotically monopole field near or outside the LC compensates for the aberration and time-delay effects. The emission from the same field line near the current sheet will arrive in phase to produce the caustics. It is shown that the “separatrix” model can naturally produce the double-peak $\gamma$-ray light curve in agreement with the Fermi observed ones. The “separatrix” model with the opposite hemisphere emission is similar to the SG model and is different from the one-hemisphere OG model. The geometric $\gamma$-ray light curves of the SG and OG model were produced in the force-free and vacuum magnetosphere [92]. It is found that the peak phases of the force-free light curves are shifted to the later phases with the magnetic poles compared to the vacuum ones, which is caused by a larger offset toward the trailing side of the PC due to the larger PC size and higher field line sweepback (see Figure 7). The force-free light curves have larger phase lags than those observed by Fermi-LAT, indicating that either the field structure is not close to the force-free one or the emission along the field lines is non-uniform. The force-free field was used to produce the radio and $\gamma$-ray light curves; the magnetic inclination angles and viewing angles are then constrained through simultaneously fitting time-aligned radio and $\gamma$-ray light curves [93,94]. It is also shown that additional constraints can be obtained by fitting the radio polarization data with a rotating vector model.

The phase-resolved spectra and energy-dependent light curves of the Vela pulsar were produced in an extended slot gap and current-sheet model for the force-free magnetosphere [95]. The decreasing flux of the first $\gamma$-ray peak relative to the second one with increasing energy is attributed to the systematically larger curvature radius of the second peak than the first peak one. However, an additional phase shift needs to be added to match the peaks phases of the observed $\gamma$-ray light curves. Therefore, an azimuthal-dependent ${\mathit{E}}_{\parallel }$ distribution is suggested to give an improved match over the observed $\gamma$-ray light curves [95]. The particle trajectory method was used to model the multi-wavelength emission from the young pulsars and MSPs in the force-free magnetosphere [96,97] (see Figure 8), where the synchrotron, curvature, and synchrotron-self Compton (SSC) emissions from both the primary particles accelerated with an assumed ${\mathit{E}}_{\parallel }$ value near the current sheet and electron–positron pairs produced in the PC cascades were computed.

The synchrotron emission from the pairs produces the optical to hard X-ray spectrum. The SSC emission from the pairs produces a broad spectrum with typically peak energy between 1 and 10 GeV for most young pulsars and around 100 GeV for MSPs, and it extends to around 1 TeV for Crab. The pair SSC emission from MSPs peaks at a higher energy than most young pulsars due to the higher pair energy, but the MSPs SSC component is suppressed to the low flux level by the Klein–Nishina effects. Therefore, the SSC emission from the pairs has a negligible contribution to the observed spectrum for most young pulsars and MSPs; only the Crab-like pulsars produce the detectable SSC components. The curvature emission from the primary particles contributes to the Fermi spectrum from several GeV up to 100 GeV, and the tail of this emission can explain the HESS-II and MAGIC spectrum from Vela, Geminga and B1706-44. The SSC emissions from the primary particles produce a very high-energy spectrum up to around 30 TeV [98]. This spectrum component is well below current detection thresholds for most pulsars, but it is still detectable for the Crab and other pulsars by the High-Altitude Water Cherenkov Observatory and Cherenkov Telescope Array. The detected maximum photon energy can provide a direct lower limit on the maximum energy of the accelerated particles. The maximum photon energies above 1 TeV from Crab and Vela require the very high particle energies of $\gamma >{10}^7$, which tends to favor the curvature radiation over the synchrotron radiation as the dominant Fermi $\gamma$-ray emission mechanism. The pulsar polarization measurements can provide an independent constraint on the location and geometry of the emission region and emission mechanism. The pulsar multi-wavelength polarization was also modeled in the force-free field structure [99], including synchrotron emission at the optical to hard X-ray band from the cascade pairs and the synchrotron or curvature emission in the $\gamma$-ray band from the primary particles. It is found that a fast position angle (PA) swing and depolarization is produced during the light-curve peaks in all energy bands. The PA swings and depolarization increase as the emission radius increases; there is a nearly ${180}^{\circ }$ PA swing for emission in the current sheet outside the light cylinder due to the change of the magnetic polarity. The curvature emission is predicted to produce a higher polarization degree than the synchrotron or inverse-Compton emission. Therefore, the $\gamma$-ray polarization measurements can distinguish whether the synchrotron or curvature emission is the dominant Fermi $\gamma$-ray emission mechanism. It is also suggested that the detection of a sudden change in PA and a sharp rise in polarization degree at 1–100 MeV may indicate that curvature emission is the Fermi $\gamma$-ray emission mechanism.

6. The Resistive Magnetospheres

6.1. Field Structure of the Resistive Magnetospheres

The force-free solutions provide a useful first step to understand a more realistic plasma-filled magnetosphere. However, the force-free solutions cannot allow any accelerating electric fields along magnetic field lines, and they cannot provide any information about the particle acceleration and the production of radiation. A more realistic pulsar magnetosphere should allow for a local dissipative region in the magnetosphere to produce the observed pulsed emission. A new degree of freedom is required to include the dissipation effect in the magnetosphere. The vacuum solutions have ${\mathit{E}}_{\parallel }$ but have no particles to be accelerated, while the force-free solutions have enough particles but no ${\mathit{E}}_{\parallel }$ to accelerate the particles. Thus, the realistic pulsar magnetosphere should lie between the vacuum limit and force-free limit. Therefore, the resistive magnetospheres are developed to include the dissipation effect by introducing a conductivity parameter, which can span the magnetospheric solutions from the vacuum to force-free limit, and the spatial distribution of the ${\mathit{E}}_{\parallel }$ value is self-consistently controlled by the $\sigma$ value. The current density can be defined by the conductivity parameter $\sigma$ in an Ohm’s law form. However, there is no unique prescription for the exact expression of the resistive current density. A formulation of the resistive current density can be defined by Ohm’s law in the fluid rest frame where the $\mathit{E}$ and $\mathit{B}$ fields are parallel, and the current density $\mathit{J}$ in the observer frame can be obtained by the Lorentz transform with the minimal velocity hypothesis [20].

$$\begin{array}{c}\hfill \mathit{J}=\frac{{\rho }_{\mathrm{e}}\mathit{E}\times \mathit{B}+\sqrt{\frac{B^2+E_0^2}{B_0^2+E_0^2}}\sigma E_0(B_0\mathit{B}+E_0\mathit{E})}{B^2+E_0^2},\end{array}$$

(12)

where $E_0$ and $B_0$ are defined by the Lorentz invariant relations

$$\begin{array}{c}\hfill B_0^2-E_0^2={\mathit{B}}^2-{\mathit{E}}^2,E_0{B}_0=\mathit{E}·\mathit{B},E_0\ge 0.\end{array}$$

(13)

A simpler current density is also defined in Ohm’s law in the observer frame by [21,22].

$$\begin{array}{c}\hfill \mathit{J}={\rho }_{\mathrm{e}}\frac{\mathit{E}\times \mathit{B}}{B^2+E_0^2}+\sigma {\mathit{E}}_{\parallel },\end{array}$$

(14)

It is expected that different current density descriptions should obtain similar field line structures and ${\mathit{E}}_{\parallel }$ distributions.

The time-dependent Maxwell Equation (1) with the resistive current density (12) and (14) were numerically solved by the FDTD method [20,21] and followed by the spectral method [22]. All these time-dependent simulations obtain a range of similar resistive solutions, which spans from the vacuum solution to the force-free solution with increasing $\sigma$ and shows the self-consistent ${\mathit{E}}_{\parallel }$ distributions in the magnetospheres. Figure 9 shows the distributions of the magnetic field line and accelerating electric field for a ${45}^{\circ }$ resistive magnetosphere with $\sigma =1.5\mathsf{\Omega }$ and $\sigma =24\mathsf{\Omega }$ at the x−z plane. It is found that the field line structure is close to the vacuum solution and a strong ${\mathit{E}}_{\parallel }$ region appears only inside the LC for low $\sigma$. As the conductivity $\sigma$ increases, the resistive solution tends to the force-free solution with increasing PC size and field line sweepback, which shift the PC toward the trailing side. A strong ${\mathit{E}}_{\parallel }$ region also appears near the Y-shaped current sheet outside the LC. The ${\mathit{E}}_{\parallel }$ decreases with increasing $\sigma$, and the ${\mathit{E}}_{\parallel }$ region only exists in the PC region inside the LC and near the Y-shaped current sheet outside the LC for the high $\sigma$. Moreover, the Poynting flux also increases with increasing conductivity $\sigma$ and inclination angle $\alpha$.

6.2. Light Curve and Spectra Modeling from the Resistive Magnetospheres

The resistive magnetospheres not only provide the field structure geometries but also the self-consistent ${\mathit{E}}_{\parallel }$ distributions from the vacuum to the force-free solutions, which can be used to model the pulsar $\gamma$-ray light curves and spectra by the direct comparison with observations. By assuming the uniform emissivity along the field lines in the comoving frame, the resistive magnetospheres were used to model the geometry $\gamma$-ray light curves [100]. The geometry $\gamma$-ray light curves of the OG and SG model were produced for a range of the resistive magnetospheres from the vacuum to the force-free solution. It is found that the peak phases of the $\gamma$-ray light curves progressively shift to the later phase with respect to the magnetic pole as the $\sigma$ increases, which is attributed to a larger offset toward the trailing side of the PC with increasing $\sigma$. The geometry $\gamma$-ray light curves were also produced from the spectral method simulation in the resistive magnetospheres [101], and a similar progression in the $\gamma$-ray light curve shape as $\sigma$ increases is also found (see Figure 10). Based on the self-consistent ${\mathit{E}}_{\parallel }$ distributions from the resistive magnetospheres, the $\gamma$-ray light curves were first modeled through computing the curvature radiation along the particle trajectory [100]. It is found that these light curves generally have larger phase lags than those in the geometry model. A further exploration of the $\gamma$-ray emission pattern and light curves from the curvature radiation in resistive pulsar magnetospheres was performed by [102]; it is found that all the emission comes from the inner magnetosphere for the low $\sigma$ values. As the $\sigma$ increases, the emission gradually moves outward to the high attitude, and a significant part of the emission is produced in regions near the Y-shape current sheet outside the LC. All the emission is produced in regions near the Y-shape current sheet outside the LC for the very high $\sigma$ values. The light curves are broad for the low $\sigma$ values and significantly become narrower as the $\sigma$ increases. By performing the comparison of the model radio-lag ($\delta$) versus peak-separation ($\mathsf{\Delta }$) distribution and the observed one, it was found that the low $\sigma$ models give the poorest match to the observed $\delta -\mathsf{\Delta }$ distribution [102]. A relatively good match over the observed $\delta -\mathsf{\Delta }$ distribution can be obtained for the high $\sigma$ models with the emission from the current sheet outside the LC , but there are still some model points on the $\delta -\mathsf{\Delta }$ diagram that do not lie near the observed ones. Motivated by the results of the high $\sigma$ models, the Force-free Inside Dissipative Outside (FIDO) model with the near-force regime inside the LC and the dissipative regime outside the LC is developed to produce the $\gamma$-ray light curves. It is found that the FIDO model can best explain the light curve characteristics of the Fermi $\gamma$-ray pulsars, particularly the observed $\delta -\mathsf{\Delta }$ distribution (see Figure 11). The FIDO model shows the emission from the 3D current sheet is non-uniform and non-axisymmetric. The self-consistent ${\mathit{E}}_{\parallel }$ distribution from the resistive magnetosphere was used to model the $\gamma$-ray light curves from the curvature radiation [101], and it is found that the FIDO model with the non-uniform emission in the current sheet provides a better match to the observed $\gamma$-ray light curves.

In the frame of the FIDO model with a range of conductivity $\sigma$ outside the LC, the properties of the energy-dependent $\gamma$-ray light curves and spectra from Fermi pulsars were explored by [104]. It was found that the variation of the spectral index and cutoff energy as a function of pulse phase can be explained reasonably. The FIDO model can also reproduce the observed trend of the $\gamma$-ray light curves with the decreasing ratio of the first to second $\gamma$-ray peak toward higher energies. It is indicated that the conductivity $\sigma$ increases with the spin-down power $\dot{\mathcal{E}}$, and it is expected that the higher $\dot{\mathcal{E}}$ pulsars have more efficient pair cascades that can provide the magnetosphere conductivity [104]. However, this study only uses approximate ${\mathit{E}}_{\parallel }$ values from the corresponding force-free solutions instead of the resistive solutions. To restrict the ${\mathit{E}}_{\parallel }$ regions only near the current sheet outside the LC, the FIDO model was refined through applying magnetic-field-line-dependent conductivity $\sigma$ [103]. A relation between $\sigma$ and $\dot{\mathcal{E}}$ was found by reproducing the observed cutoff energy with the refined FIDO model. It is shown that the $\sigma$ increases with $\dot{\mathcal{E}}$ for the high $\dot{\mathcal{E}}$ and saturates for the low $\dot{\mathcal{E}}$. The refined FIDO model can reproduce the observed trends of the $\gamma$-ray luminosity $L_{\gamma }$ as a function of $\dot{\mathcal{E}}$. It is suggested that higher emitting particle multiplicities are needed to reproduce the observed $L_{\gamma }$ at the high $\dot{\mathcal{E}}$ (see Figure 11). Moreover, the FIDO model with the self-consistent ${\mathit{E}}_{\parallel }$ values from the resistive magnetosphere was used to explore the energy-dependent $\gamma$-ray light curves and spectra from the Fermi pulsars [105]. It is also found that the emission from the FIDO model is non-uniform and asymmetric over the PC, the high emission is concentrated on the leading side of the PC edges, and the FIDO model can reproduce the phase-averaged spectra and energy-dependent $\gamma$-ray light curves for the Vela and Crab pulsar, especially the observed trend of the decreasing ratio of the first to second $\gamma$-ray peak with increasing energies [105].

7. The PIC Magnetospheres

7.1. Field Structure of the PIC Magnetosphere

The resistive magnetospheres can produce the accelerating electric fields that are self-consistent with the magnetic field structures by introducing a macroscopic conductivity. However, the resistive magnetospheres are still not self-consistent, since they cannot model the source of particles that provide magnetosphere conductivity. The kinetic plasma simulation is required to self-consistently treat the feedback between the particle motions, radiating photons, and electromagnetic fields. The kinetic model provides a method to self-consistently compute the feedback between the particle motions, radiating photons and electromagnetic fields by solving the time-dependent Vlasov–Maxwell equations. The time-dependent Vlasov equations are given by

$$\begin{array}{c}\hfill \frac{\mathsf{\partial }f}{\mathsf{\partial }t}+\frac{\mathit{p}}{\gamma m}\frac{\mathsf{\partial }f}{\mathsf{\partial }\mathit{r}}+q\left(\mathit{E}+\frac{\mathit{v}\times \mathit{B}}{c}\right)\frac{\mathsf{\partial }f}{\mathsf{\partial }\mathit{r}}=0,\end{array}$$

(15)

where $\mathit{p}$ is the particle momentum and $f(\mathit{r},\mathit{p})$ is the particle distribution function. It is difficult to directly solve the time-dependent Vlasov–Maxwell equations in the six-dimensional phase space $f(\mathit{r},\mathit{p})$. A particle-in-cell (PIC) method can be used to indirectly solve the time-dependent Vlasov–Maxwell equations by sampling the distribution function in phase space with the particles. Each particle represents a large number of real particles that approximately follow the exact same path in phase space. The Vlasov equation can be solved along characteristics curves by a set of ordinary differential equations

$$\begin{array}{cc}\hfill \frac{d\mathit{r}}{dt}& =\mathit{v},\hfill \\ \hfill \frac{d\mathit{v}}{dt}& =q\left(\mathit{E}+\frac{\mathit{v}\times \mathit{B}}{c}\right).\hfill \end{array}$$

(16)

which corresponds to the particle motion equations. The particle motion in Equation (13) produces a set of characteristic curves representing a surface in phase space, which is a solution of the Vlasov equation. The currents and charges from the particle motions are deposited in the numerical grid by following the charge-conserving scheme

$$\begin{array}{cc}\hfill \rho & =\sum _j{q}_j\delta (\mathit{r}-{\mathit{r}}_j),\hfill \\ \hfill \mathit{J}& =\sum _j{q}_j{\mathit{v}}_j\delta (\mathit{r}-{\mathit{r}}_j).\hfill \end{array}$$

(17)

which can be used as the sources to solve the Maxwell equations.

The time-dependent Maxwell Equation (1) and the motion Equation (17) for an aligned rotator were first solved through the 2D Cartesian PIC Tristan-MP code by Philippov and Spitkovsky (2014) [23]; the structure of pulsar magnetospheres was explored by injecting the particles only from the stellar surface or everywhere in the magnetosphere. The PIC simulation of an aligned rotator was performed through the 2D spherical PIC code with the incorporation of the full pair cascades from the stellar surface to beyond the LC [24]. The PIC simulations of an aligned rotator by the 2D spherical PIC code were also performed with the particle injection from the ${\mathit{E}}_{\parallel }\ne 0$ region [24] or only above the stellar surface [26]. Furthermore, the PIC simulations of an oblique rotator were first performed through the 3D Cartesian PIC Tristan-MP code by Philippov and Spitkovsky (2015) [27]. In this case, approximate pair cascades are implemented by injecting the particles whose energy exceeds a given threshold. The general-relativistic corrections with the frame-dragging and time–space curvature were later included in the 3D Cartesian PIC Tristan-MP code [106,107]. The PIC simulations of an oblique rotator were also performed by the 3D spherical PIC Zeltron code with the particle injection from only above the stellar surface [108]. These PIC simulations show that magnetospheres transition from a charge-separated “electrosphere” solution with the disk–dome structure for the low particle injection to a near force-free solution with the particle acceleration near the current sheet for the high particle injection. It is also found that the pair cascades from the stellar surface cannot fill the magnetosphere to produce a near force-free solution; the pair cascades in the current sheet are needed to produce a near force-free solution. Recently, a 3D Cartesian PIC code was developed to perform the PIC simulations of an oblique rotator [28], and the magnetosphere transition from the near-vacuum solution to the near force-free solution was explored by applying the particle injection rates from the low to high values. Figure 12 shows the distributions of the magnetic field line and accelerating electric field in the PIC simulation of a ${45}^{\circ }$ inclined rotator with different particle injection rates at the x−z plane. It is shown that the ${\mathit{E}}_{\parallel }$ regions are confined only near the current sheet outside the LC as the particle injection rate increases, while the volume of the ${\mathit{E}}_{\parallel }$ region decreases with the increasing particle injection rate. A similar trend of the ${\mathit{E}}_{\parallel }$ distribution with increasing particle injection rate was also found by applying the magnetic-field-line-dependent particle injection in the PIC simulation [109]. It is also found that a near force-free magnetosphere can be produced at all inclination angles without the need for pair cascades near the current sheet [29]. All the PIC simulations show a near force-free magnetosphere with the particle acceleration near the current sheet, revealing that the particle acceleration and the high-energy emission is mainly produced near the current sheet outside the LC. Moreover, the reported Poynting flux from the PIC simulation is similar to the MHD expectation. All the PIC simulations show the ∼1–20% dissipation of the Poynting flux outside the LC, which is converted into particle acceleration and radiation in the current sheet.

7.2. Light Curve and Spectra Modeling from the PIC Magnetospheres

The PIC magnetospheres have the near force-free field structures with the ${\mathit{E}}_{\parallel }$ regions near the current sheet outside the LC, which can be connected to the pulsar observational signatures by simultaneously extracting the pulsar emission information from the current sheet in the PIC simulation. The PIC simulation cannot use the realistic physical parameters to model the pulsar magnetosphere because of the large-scale separation between the plasma frequency and the pulsar rotation frequency. The artificially low magnetic field of B∼${10}^5\mathrm{G}$ is usually used in the PIC simulations to resolve the plasma frequency, so that the particles are only accelerated to the low energy of $\gamma$∼${10}^3$. This is clearly not enough to produce the observed $\gamma$-ray photons by the synchrotron or curvature radiation. Therefore, a scale-up technique for the magnetic field and energy is required to model the pulsar $\gamma$-ray emission. The pulsar $\gamma$-ray emission was studied by modeling the synchrotron radiation from the PIC particles accelerated by magnetic reconnection in the current sheet [107,108], which can produce the $\gamma$-ray light curves and spectra by scaling up the synchrotron photons of the PIC particles to those of the real pulsar. It is found that the sky maps from the PIC model show the two bright caustics that are characteristic of radiation from the current sheet. The double-peak $\gamma$-ray light curves can generally be produced when the observer angle crosses the two bright caustics, but the phase lags from the magnetic pole are seemly larger that those seen in Fermi $\gamma$-ray light curves. The predicted $\gamma$-ray spectra show the power-law exponential cutoff shapes with the cutoff energies around several GeV, which are similar to those of the Fermi $\gamma$-ray spectra. The pulsar $\gamma$-ray emission was also explored by modeling the curvature radiation from the scaled-up PIC particles [28,109], where the surface magnetic field, the accelerating electric field, and the particle energy are scaling up to those of the real pulsar. It is found that the particles can be accelerated to the realistic emitting $\gamma$-ray energies of $\gamma$∼${10}^7$ in the radiation-reaction limited regime where the particle acceleration is balanced by the curvature losses. The predicted $\gamma$-ray spectra with the cutoff energies of several GeV are also similar to those of the Fermi $\gamma$-ray spectra. Moreover, the predicted $\gamma$-ray light curves have smaller phase lags from the magnetic pole than those of the PIC synchrotron model, which has better agreement with those of the Fermi $\gamma$-ray light curves. It is also found that there is a positive relation between the the particle injection rate and the spin-down power by reproducing the observed $\gamma$-ray cutoff energies and luminosities, which provides a physical link between the microscopic injection rate and macroscopic conductivities. The pulsar $\gamma$-ray emission was further explored by the refined PIC model with the magnetic-field-line dependent particle injection. It is shown that the PIC magnetosphere with the particle injection near the separatrix zone can better reproduce the broad properties of the Fermi $\gamma$-ray pulsar phenomenology, including the $\gamma$-ray light curve profiles, the spectra shapes, and the $\delta -\mathsf{\Delta }$ relation [109].

The pulsar $\gamma$-ray polarization was also predicted by modeling the synchrotron radiation from the current sheet in the PIC simulation [110]. Figure 13 shows the predicted $\gamma$-ray polarization properties from the PIC synchrotron model for the Crab-like pulsar. It is found that each pulse peak is accompanied by a nearly ${180}^{\circ }$ swing of the polarization angle. This is the natural result of the emission pattern from the current sheet where the the field line polarity changes as the line of sight crosses the current sheet. It is also shown that the synchrotron radiation predicts a lower polarization degree (≲$30\%$) than that of the curvature curvature (≳$40\%$). The future $\gamma$-ray polarization measurements can help us distinguish the curvature origin from the synchrotron one for the observed Fermi $\gamma$-ray emission.

8. The Combined Force-Free and AE Magnetospheres

8.1. Field Structure of the Combined Force-Free and AE Magnetospheres

The PIC simulation cannot use the real pulsar parameters to model the pulsar magnetosphere and predict the subsequent radiation because of the large ratio between the macroscopic scales and the microscopic scales. Therefore, the scale-up method by extrapolating the simulation parameters to real pulsar parameters is required to model the pulsar high-energy emission. However, it is dangerous for any extrapolation to the realistic values in the nonlinear emission scenarios. Therefore, a clever method is required to perform realistic simulation of the pulsar magnetosphere and predict the subsequent pulsed emission.

It is expected that the pulsar magnetosphere is filled with the $e^{\pm }$ pairs created by the pair cascades. These particles can be accelerated by the unscreened ${\mathit{E}}_{\parallel }$ to high energies and produce the multi-wavelength radiation photons by the synchrotron, curvature, and inverse-Compton scattering processes. These photons have a back-reaction onto the particle motion and produce a radiative friction against the motion. The radiative friction brakes the particles in a direction opposite to their motion. We can expect that the particle acceleration and radiation can reach a stationary balance in the magnetosphere, which is called Aristotelian electrodynamics or radiation reaction limit. The balance condition between the Lorentz force and radiative friction gives the following [37,111,112].

$$\begin{array}{c}\hfill q(\mathit{E}+\mathit{v}\times \mathit{B})=K\mathit{v}.\end{array}$$

(18)

where K is a positive constant parameter. By taking the cross-product to Equation (18) with $\mathit{B}$ and using a vector calculus, we can obtain the particle velocity

$$\begin{array}{c}\hfill \left(B^2+\frac{K^2}{q^2}\right)\mathit{v}=\frac{K}{q}\mathit{E}+\mathit{E}\times \mathit{B}+\frac{q}{K}(\mathit{E}·\mathit{B})\mathit{B}.\end{array}$$

(19)

By taking the dot product to Equation (19) with $\mathit{v}$ and using a vector calculus, Equation (19) becomes

$$\begin{array}{c}\hfill K^4{v}^2-q^2(E^2-v^2{B}^2)K^2-q^4{(\mathit{E}·\mathit{B})}^2=0.\end{array}$$

(20)

Assuming that the particle velocity is equal to the light velocity ($v=1$), we have

$$\begin{array}{c}\hfill K^2=\frac{q^2}{2}\left[E^2-B^2\pm \sqrt{{(E^2-B^2)}^2+4{(\mathit{E}·\mathit{B})}^2}\right].\end{array}$$

(21)

By using the Lorentz invariants (13), we find that the constant K is the solution of the following equations

$$\begin{array}{cc}\hfill E^2-B^2& =K^2/q^2-B_0^2,\hfill \\ \hfill \mathit{E}·\mathit{B}& =K{B}_0/\left|q\right|.\hfill \end{array}$$

(22)

Substituting Equation (22) into Equation (19), we can retrieve the AE velocity expression [113,114,115].

$$\begin{array}{c}\hfill {\mathit{v}}_{\pm }=\frac{\mathit{E}\times \mathit{B}\pm (B_0\mathit{B}+E_0\mathit{E})}{B^2+E_0^2},\end{array}$$

(23)

where the two signs correspond to positrons and electrons, they follow a different trajectory that only depends on the local electromagnetic fields in the magnetosphere. The AE current density can be derived by Equation (23).

$$\begin{array}{c}\hfill \mathit{J}=\frac{{\rho }_e\mathit{E}\times \mathit{B}+{\rho }_0(B_0\mathit{B}+E_0\mathit{E})}{B^2+E_0^2},\end{array}$$

(24)

where ${\rho }_0={\rho }_{+}+{\rho }_{-}$ is the total charge density. The AE current density can by defined by giving a description for the total charge density ${\rho }_0$. However, there is no unique prescription for the current density, since we do not know the exact expression of the total charge density. A formulation of the AE current density can be defined by introducing a pair multiplicity $\kappa$ [32,116].

$$\begin{array}{c}\hfill \mathit{J}={\rho }_e\frac{\mathit{E}\times \mathit{B}}{B^2+E_0^2}+(1+\kappa )\left|{\rho }_e\right|\frac{(B_0\mathit{B}+E_0\mathit{E})}{B^2+E_0^2},\end{array}$$

(25)

where $|{\rho }_e|$ can be understood as the charge density of the primary electrons, and $\kappa |{\rho }_e|$ can be understood as the charge density of the secondary pairs from the pair cascades. Therefore, the pair multiplicity $\kappa$ is physically related with the pair cascade processes in the magnetosphere. When $\kappa =0$, the AE current density is exactly consistent with the one given by [115].

The resistive model can produce the local dissipations by relaxing the force-free condition, but they cannot self-consistently include the back-reaction of emission onto particle dynamics. The PIC simulation can model the pulsar magnetospheres by including the self-consistent feedback between the particle motions, the radiating photons, and the electromagnetic fields, but they cannot use the real pulsar parameters to model the pulsar magnetospheres due to the large ratio between the macroscopic scales and the microscopic scales. The combined force-free and AE model is a good approximation between the resistive and PIC models, which can include the back-reaction of the emitting photons onto particle dynamics and allows for the local dissipations to produce the observed pulsar emission with the real pulsar parameters. A clever method with the combined force-free and AE was then developed to model the pulsar magnetosphere by using the force-free description where $E\le B$ and the AE description where $E>B$. The combined force-free and AE magnetospheres for the oblique rotator in the limit of $\kappa =0$ were first computed by the finite-difference method [30]. This study was extended by introducing a new pair multiplicity, and the combined force-free and AE solutions for the aligned rotator are simulated by the spectral method [31]. The combined force-free and AE solutions of the oblique rotator with a prescribed pair multiplicity were first obtained by Cao and Yang (2020) [32] and then by Pétri (2022) [33], who both used the 3D spectra method to numerically solve the time-dependent Maxwell Equation (1) with the AE current density (31). The high-resolution simulations of the combined force-free and AE magnetosphere for the oblique rotator were also explored to accurately resolve the ${\mathit{E}}_{\parallel }$ regions in the current sheet by the 3D spectral method [117]. All these time-dependent simulations obtain a family of the similar field structures with the ${\mathit{E}}_{\parallel }$ regions near the current sheet outside the LC for a range of the $\kappa$ values. As the pair multiplicity $\kappa$ increases, the combined force-free and AE magnetosphere tends to the force-free solution with the current sheet outside the LC, and the ${\mathit{E}}_{\parallel }$ region is more restricted to the current sheet outside the LC. The ${\mathit{E}}_{\parallel }$ regions decrease with increasing $\kappa$ values, and the ${\mathit{E}}_{\parallel }$ region is only confined to the area near the current sheet outside the LC for the high $\kappa$ values. Figure 14 shows the distributions of the magnetic field lines and accelerating electric fields in the combined force-free and AE magnetospheres for different inclined angles with the pair multiplicity $\kappa =3$ at the x−z plane. It is seen that all these solutions have a near force-free magnetosphere with the ${\mathit{E}}_{\parallel }$ region only near the current sheet for all the inclination angles at the high $\kappa$ value, which indicates that the particle acceleration and the high-energy emission mainly occur at the current sheet near the LC.

8.2. Light Curve and Spectra Modeling from the Combined Force-Free and AE Magnetospheres

The combined force-free and AE magnetosphere provide a near force-free field structure with the ${\mathit{E}}_{\parallel }$ distribution only near the current sheet outside the LC. We can use these field structures and ${\mathit{E}}_{\parallel }$ distributions to predict the pulsar $\gamma$-ray emission by integrating the particle trajectory with the AE velocity. The AE velocity expressions are extended to new AE velocity expressions by including the Landau–Lifshitz terms. It is found that the new AE velocity expressions give the similar results in accuracy to the standard AE velocity expressions [112,118,119]. It is also demonstrated that the AE velocity can correctly follow the particle trajectory in the strong electromagnetic field, which can be implemented to compute realistic pulsar light curves and spectra in the arbitrary pulsar magnetosphere [116,120]. The pulsar $\gamma$-ray light curves and spectra can be computed by the synchro-curvature radiation from each radiating particle [50,121,122].

$$\begin{array}{c}\hfill \frac{d{N}_{sc}}{dE}=\frac{\sqrt{3}{e}^3\gamma y}{4\pi \hslash r_{\mathrm{eff}}}\left[(1+z)F\left(y\right)+(1-z)k_{2/3}\left(y\right)\right],\end{array}$$

(26)

with

$$\begin{array}{c}\hfill y=\frac{E}{E_c},E_c=\frac{3}{2}\hslash c{Q}_2{\gamma }^3,z={\left(Q_2{r}_{\mathrm{eff}}\right)}^{-2},\end{array}$$

(27)

$$\begin{array}{c}\hfill F\left(y\right)=y{\int }_y^{\infty }{K}_{5/3}\left(y\right)dy,\end{array}$$

(28)

$$\begin{array}{c}\hfill Q_2=\frac{{cos}^2\psi }{\rho }{\left(1+3\xi +{\xi }^2+\frac{r_g}{\rho }\right)}^{1/2},\end{array}$$

(29)

$$\begin{array}{c}\hfill r_{\mathrm{eff}}=\frac{\rho }{cos\psi }{\left(1+\xi +\frac{r_g}{\rho }\right)}^{-1},\end{array}$$

(30)

$$\begin{array}{c}\hfill r_g=\frac{m{c}^2\gamma sin\psi }{eB},\xi =\frac{\rho {sin}^2\psi }{r_g{cos}^2\psi },\end{array}$$

(31)

where E is the photon energy, $E_c$ is the characteristic energy of the synchro-curvature photon, $\gamma$ is the particle Lorentz factor, and $\rho$, $\psi$ and $r_g$ are the particle curvature radius, the pitch angle and the Larmor radius, respectively. The parameter $\xi$ determines the dominant synchro-curvature regime with $\xi$$>>$1 at the synchrotron regime and $\xi$$<<$1 at the curvature regime.

The properties of the pulsar $\gamma$-ray light curves and spectra from the curvature radiation in the combined force-free and AE magnetosphere were first explored by Cao and Yang (2022) [117], where the AE velocity was used to determine the particle trajectory, and the Lorentz factor along each trajectory was computed by including the self-consistent ${\mathit{E}}_{\parallel }$ values and the curvature radiation losses. When the direction of the photon emission along the direction of particle motion is assumed, the pulsar energy-dependent $\gamma$-ray light curves and spectra are produced by collecting the curvature photons from all of the radiating particles in sky maps. Figure 15 shows the sky maps and the corresponding light curves for different inclination angles and view angles in the combined force-free and AE magnetospheres with the pair multiplicity $\kappa =3$. The sky maps from the combined force-free and AE magnetosphere show the characteristics of the caustic emission from the current sheet, which are similar to those of the PIC model and the other current sheet emission model. The double-peak $\gamma$-ray light curves can generally be produced for a broad range of inclination angles and view angles, and the phase lags from the magnetic pole are in good agreement with the Fermi observation for the high pair multiplicity. The predicted $\gamma$-ray spectra can reach the $\gamma$-ray energies with the cutoff energies around several GeV in agreement with those of the Fermi $\gamma$-ray spectra, which do not need the scale-up method as in the PIC simulations, since they can be computed for realistic pulsar parameters. The pulsar $\gamma$-ray light curves and spectra are further explored by performing a direct comparison with the Fermi observation for the Vela pulsar [123]. It is found that the predict energy-dependent $\gamma$-ray light curves can well reproduce the observed trend of the $\gamma$-ray light curves with the energy evolution for the Vela pulsar (Figure 16); the decreasing ratio of the first peak to the second one with increasing energy is explained as the systematically larger curvature radius and the higher cutoff energy in the second peak. The predicted $\gamma$-ray spectra can also well explain the phase-averaged spectra and phase-resolved spectra of the first and second peaks for the Vela pulsar. The observed higher spectral flux and cutoff energy in the second peak are also well reproduced, because the second peak of the $\gamma$-ray light curves survives at higher energy than the first peak one.

9. Conclusions

The last few decades have witnessed significant advances in the pulsar $\gamma$-ray observations made by the Fermi Gamma-Ray Space Telescope and ground-based Cherenkov telescopes. These $\gamma$-ray observations have obtained valuable observational data about the pulsar $\gamma$-ray light curves and spectra in the GeV−TeV range. The modeling of the Fermi $\gamma$-ray light curves and spectra locates the $\gamma$-ray emission region at the outer magnetosphere near the LC or beyond the LC. The TeV $\gamma$-ray data from ground-based Cherenkov telescopes provide other important constraints on the particle acceleration and emission mechanics in the pulsar magnetospheres. In parallel to the advances in the pulsar $\gamma$-ray observations, the self-consistent modeling of the pulsar magnetosphere has also made significant breakthroughs from the early vacuum model to the recent plasma-filled magnetospheres by numerical simulations. The current numerical simulations confirmed the general picture that a typical pulsar magnetosphere has a near force-free field structure with the particle acceleration in the separatrix near the LC and the current sheet outside the LC, which can provide a good match to the recent high-energy $\gamma$-ray observations. However, there is disagreement over whether the $\gamma$-rays emission are produced by the synchrotron radiation or the curvature radiation from the current sheet near the LC. The light curve and spectral information is not sufficient to distinguish between different emission models and emission mechanisms in the magnetospheres. This degeneracy can be broken by using the additional constraints from the future polarization measurements. The modeling of the combined multi-wavelength light curves, spectra, and polarization would put a stronger constraint on the the location and geometry of the particle acceleration and emission region as well as the emission mechanism. Realistic PIC simulations are prohibited by the large ratio between the plasma frequency and the pulsar rotation frequency. A more clever idea is required to perform the realistic meaningful simulations of the pulsar magnetosphere and predict the realistic pulsar emission.