ELECTRIC FIELD SCREENING WITH BACKFLOW AT PULSAR POLAR CAP

In a pulsar magnetosphere, particles are significantly accelerated at the given regions, and emit electromagnetic radiation from radio to γ-ray wavelengths. Observations by the Fermi Gamma-Ray Space Telescope have shown that the differential spectra above 200 MeV are described well by the power-law functions with an exponential cutoff, and that cutoff shapes that are sharper than the simple exponential cutoff are rejected with high significance (e.g., Abdo et al. 2009). This rules out the near-surface emission proposed in the polar cap cascade model (Daugherty & Harding 1996), which would exhibit a much sharper spectral cutoff due to the attenuation of the magnetic pair creation. Hence, the detected γ-ray pulse emission should originate from the outer region of the magnetosphere, as considered in the outer gap model (e.g., Cheng et al. 1986; Romani 1996; Hirotani 2006, 2015; Takata et al. 2006, 2016), as well as the current sheet model (e.g., Kirk et al. 2002; Bai & Spitkovsky 2010; Kalapotharakos et al. 2014; Brambilla et al. 2015; Cerutti et al. 2016).

On the other hand, the region just above the neutron star surface (NSS) has been considered as the site of the radio pulsed emission (e.g., Noutsos et al. 2015). The mechanism of pulsar radio emission is established as a coherent process, so that the plasma dynamics near the NSS would be strongly related to the emission mechanism (e.g., Sturrock 1971). The two-stream instability is a promising process for creating the plasma bunches (e.g., Ruderman & Sutherland 1975). The curvature radiation from the bunches is usually discussed as the mechanism of the coherent radio emission (e.g., Saggion 1975). In order to investigate the possible instabilities near the NSS, we should take into account the non-stationary effects in the plasma flows.

The dynamics of the plasma and the electromagnetic field near the NSS highly depend on the ratio of the current density parameter along the magnetic field, ${j}_{{\rm{m}}}$, to the Goldreich–Julian (GJ) value, ${j}_{{\rm{GJ}}}={\rho }_{{\rm{GJ}}}c$ (Mestel et al. 1985; Shibata 1997; Timokhin & Arons 2013), which is characterized by the GJ charge density ${\rho }_{{\rm{GJ}}}=-{\boldsymbol{\Omega }}\cdot {\boldsymbol{B}}/2\pi c$ (Goldreich & Julian 1969), where ${\boldsymbol{\Omega }}$ is the stellar angular velocity vector and ${\boldsymbol{B}}$ is the local magnetic field vector. The parameter ${j}_{{\rm{m}}}$ is regulated by the twist of the magnetic field (${\rm{\nabla }}\times {\boldsymbol{B}}$) imposed by the global stress balance of the pulsar magnetosphere (e.g., Shibata 1991). In the polar cap region near the NSS, an accelerating electric field spontaneously develops to adjust the current and charge densities to the current density parameter ${j}_{{\rm{m}}}$ and the GJ charge density ${\rho }_{{\rm{GJ}}}$. In the cases ${j}_{{\rm{m}}}/{j}_{{\rm{GJ}}}\leqslant 0$ and ${j}_{{\rm{m}}}/{j}_{{\rm{GJ}}}\geqslant 1$, outflowing particles from the NSS alone cannot adjust the current and charge densities to ${j}_{{\rm{m}}}$ and ${\rho }_{{\rm{GJ}}}$ simultaneously (e.g., Mestel et al. 1985). In such cases, a significant accelerating electric field develops and causes copious pair creation. The newly created pairs then screens the accelerating electric field for a temporary period of time (e.g., Sturrock 1971; Levinson et al. 2005; Timokhin & Arons 2013).

The backflowing particles from the outer acceleration region (OAR) modify the above description of the dynamics near the NSS. As a result of discharge at the OAR in the magnetosphere (e.g., outer gap or current sheet), some fraction of charged particles comes back to the NSS. Such backflowing particles are actually seen in numerical studies (e.g., Hirotani 2006; Wada & Shibata 2007; Chen & Beloborodov 2014; Cerutti et al. 2015). The existence of the backflow is favorable for explaining the observed pulse profiles in the non-thermal soft γ-ray, X-ray, and optical wavelengths, whose peaks are not aligned with the GeV γ-ray peak (Takata et al. 2008; Kisaka & Kojima 2011; Wang et al. 2013). The backflowing particles have also been considered to heat the NSS around the magnetic pole, whose signature is observed as the thermal pulsed emission in the soft X-ray band (e.g., Zavlin & Pavlov 2004). The threshold of the occurrence of pair cascades near the NSS depends on the contribution of the backflowing particles to the current and the charge densities (Beloborodov 2008). The outflow from the NSS would also affect the dynamics in the OAR. The outgoing particles from the NSS contribute to the particle injection rate from the inner boundary of the OAR (e.g., Takata et al. 2006).

If an outflow from the NSS affects the accelerating electric field in the OAR, the current and number densities of the backflowing particles change, and the resulting particle outflow from the NSS may be also modified. Through such non-linear interplay between the NSS and OAR, the global magnetosphere is expected to reach the steady or quasi-steady state (e.g., periodic behavior). Leung et al. (2014) and Takata et al. (2016) suggest that a non-stationary outer gap model is favored to reproduce the sub-exponential cutoff feature in the GeV γ-ray spectrum observed with Fermi. In order to understand the global behavior of the magnetosphere, we need to link the dynamics of the region above the NSS and the OAR.

In the first step, we focus on the dynamics of only the restricted region just above the NSS for given backflowing particles. If a steady electric field just above the NSS exists, copious electron–positron pairs are produced via electromagnetic cascades. Such pairs flow into the OAR, and may screen out the electric field in the OAR. Therefore, the electric field just above the NSS should almost be screened out to activate the OAR.

Local simulations have been performed to investigate the particle acceleration and the pair creation processes near the NSS (Beloborodov & Thompson 2007; Lyubarsky 2009; Timokhin 2010; Barzilay 2011; Chen & Beloborodov 2013; Timokhin & Arons 2013; Timokhin & Harding 2015). Since the realistic pair creation process with the actual mass-to-charge ratio is difficult to include in the present global simulations (Spitkovsky & Arons 2002; Wada & Shibata 2007, 2011; Yuki & Shibata 2012; Chen & Beloborodov 2014; Philippov & Spitkovsky 2014; Belyaev 2015; Cerutti et al. 2015; Philippov et al. 2015a, 2015b), the local simulations are complementary. In order to link the local region above the NSS to the global magnetospheric structures, it is useful to model the properties of the outflowing particles from the NSS for an arbitrary ratio ${j}_{{\rm{m}}}/{j}_{{\rm{GJ}}}$ and the backflow from the OAR.

Timokhin (2010), Timokhin & Arons (2013), and Timokhin & Harding (2015) performed local particle simulations to investigate the pair cascade near the NSS. In their regimes, a large number density of pairs is supplied to the OAR so that the electric field in the OAR is screened by the copious pairs. Then, the backflow from the OAR is suppressed. In this context, the effect of the backflowing particle has not been investigated in the local particle simulations so far. However, the OAR serving as a source of backflowing particles should occur if the pair cascade near the NSS fails to supply enough pairs.

In this paper, we study the screening of the accelerating electric field above the NSS, taking into account the effect of the backflowing particles from the OAR. As we have mentioned above, we consider the screening of the electric field near the NSS as a necessary condition for activating the OAR, because too much pair supply from the inner magnetosphere via a strong electric field would choke the OAR. The local condition of electric field screening near the NSS in the absence of the pair cascade is investigated for a given ratio, ${j}_{{\rm{m}}}/{j}_{{\rm{GJ}}}$, which is imposed in the global magnetospheric structure. In Section 2, we introduce our model with a particle outflow from the NSS, where the number density and momentum distribution of the backflowing particles are given as model parameters. In Section 3, we analytically show the screening condition for the velocity of the plasma flow from the NSS in the case where the backflowing particles are ultra-relativistic ($\gamma \gg 1$). We see that the development of the accelerating electric field cannot be avoided for some combinations of the total current density and the contribution of the backflowing particles. In Section 4, we introduce additional components of backflowing particles, electron–positron pairs with a mildly relativistic temperature. Particle-in-cell simulations are performed to investigate the screening conditions near the NSS. The implications of our results for the pulsar radio emission are discussed in Section 5. We summarize our work in Section 6.

We consider a local problem of how the accelerating electric field is screened near the NSS for given parameters. This study is motivated by the idea that the field screening near the NSS may be an essential condition for activating the OAR, from which high-energy photons are emitted. Our model is one-dimensional (1D), with a spatial axis along magnetic field lines. We assume that charged particles move along straight magnetic field lines, which are perpendicular to the NSS. This assumption is justified for the following reasons. One is that charged particles in the strong magnetic field are in the first Landau level and move strictly along magnetic field lines. The other is that the length scale L of our calculation domain is much smaller than the radius of the field line curvature ${R}_{{\rm{cur}}}$ and the radius of the polar cap ${r}_{{\rm{pc}}}$.⁵ We neglect the induced variations of the magnetic field that accompany variable electric fields parallel to the magnetic field, E_∥.

For the boundary condition on the NSS, we assume that electrons or ions with non-relativistic velocities are freely extracted from the NSS by the electric field E_∥ just above the NSS. This assumption is reasonable because the work function is small compared with thermal energy for most pulsars (Jones 1980; Müller 1984; Neuhauser et al. 1986).

One of the most important parameters is the magnetospheric current density parameter ${j}_{{\rm{m}}}=(c/4\pi ){\rm{\nabla }}\times {\boldsymbol{B}}$. The parameter ${j}_{{\rm{m}}}$ is induced to balance the global stress, since open magnetic field lines that pass through the light cylinder are twisted. When the actual current density j coincides with ${j}_{{\rm{m}}}$, the electric field becomes stationary. However, the condition $j={j}_{{\rm{m}}}$ is not always assured. In general, ${j}_{{\rm{m}}}$ can take any value, depending on some global conditions, because the local accelerator and the pulsar wind interplay with each other through the current. Thus, we cannot deduce the value of ${j}_{{\rm{m}}}$ from any local model. In our local model, we regard the value of ${j}_{{\rm{m}}}$ as a model parameter.

The local current density, j, tends to be quickly adjusted to ${j}_{{\rm{m}}}$ via the generation of the accelerating electric field E_∥. As a result, the character of the accelerator near the NSS strongly depends on ${j}_{{\rm{m}}}$ (Mestel et al. 1985; Shibata 1997; Beloborodov 2008; Timokhin & Arons 2013). Since the parameter ${j}_{{\rm{m}}}$ should have the magnetospheric timescale, it can be regarded as stationary on the dynamical timescale that is typical for the surface region, which is much longer than the local timescales L/c or ${r}_{{\rm{pc}}}/c$. We assume that the imposed parameter ${j}_{{\rm{m}}}$ is constant, while the local electric field and plasma flows develop with time (Timokhin 2010; Timokhin & Arons 2013).

The main targets in this paper are the effects of the backflowing particles from the OAR on the plasma dynamics near the NSS (Figure 1). Thus we consider the energetic pulsars whose OARs are sustained by the pair cascade. If electron–positron pairs whose number densities are much larger than the GJ value, ${n}_{{\rm{GJ}}}=| {\rho }_{{}_{{\rm{GJ}}}}| /e$, are flowing back to the NSS, the electric field should be easily screened (e.g., Lyubarsky 2009). Here, we focus on the cases where the number density of the backflowing particle is an order of the GJ value. For a pulsar with active OAR, this is naturally expected as follows. The OAR resides on the last-open field line, which is the boundary between the open and closed field lines (shaded region in Figure 1; e.g., Cheng et al. 1986). The returning current ($j/{\rho }_{{\rm{GJ}}}c\lt 0$) should flow along the field lines in the OAR to close the current circuit (Figure 1).⁶ The particles accelerated in the OAR should then dominantly contribute to the return current, so that the expected value of the returning current density is an order of the GJ value. The typical current density in the results of the global simulations is an order of GJ value (e.g., Kalapotharakos et al. 2012).

Zoom In Zoom Out Reset image size

In addition to the accelerated particles, the electron–positron pairs created between the OAR and NSS are also expected to come back to the NSS. Recent time-dependent particle simulations show the existence of two components in the momentum distribution, which are stationary and screen the electric field E_∥ outside the OAR (Timokhin 2010; Timokhin & Arons 2013). A beam component with a high Lorentz factor accounts for the current density as $j={j}_{{\rm{m}}}$, while another component with a quasithermal momentum distribution and an average velocity $\bar{\beta }\sim 0$ adjusts the charge density as $\rho ={\rho }_{{\rm{GJ}}}$ (Timokhin 2010; Chen & Beloborodov 2013; Timokhin & Arons 2013). In the outer gap model, the electric field E_∥ outside the OAR is screened via photon-photon pair creation near the null-charge surface (Cheng et al. 1986; Hirotani 2006; Takata et al. 2006). Most pairs are created above the OAR, because the curvature photons are emitted to the tangential direction of the field line. On the other hand, at the inner boundary of the OAR, the number density of created pairs may not be much larger than the GJ value. When the screening process with the two plasma components works near the null-charge surface, the non-relativistic particles are also expected to flow toward the NSS from the OAR. In what follows, we consider two cases for the momentum distribution of the backflowing particles with a comparable density to the GJ number density.

In case 1, we assume that only particles with a large Lorentz factor, $\gamma \gg 1$, are flowing back to the NSS. The schematic pictures of the case with a backflow are shown in Figure 2(a). The velocity of the backflowing particles is almost light speed until reaching the NSS. Hereafter, we call these particles a beam component. The beam component consists of electrons or positrons, and is accelerated in the OAR. The number density of the beam component would be an order of GJ value, ${n}_{{\rm{GJ}}}$. The Lorentz factor γ is regulated by the acceleration process in the OAR and the energy loss due to the curvature and synchrotron radiations during the travel to the NSS. As seen in the numerical results by Hirotani (2006), particles may be slightly decelerated before reaching the NSS ($\gamma \sim {10}^{5}\mbox{--}{10}^{6}$). We neglect the particle creation by the curvature photons emitted by the beam component. We analytically investigate the screening condition for case 1 in Section 3.

Zoom In Zoom Out Reset image size

In case 2, we consider the possibility that there is another component of the backflowing particles with a quasithermal momentum distribution and an average velocity $\bar{\beta }\sim 0$ (see Figure 2(b)). Hereafter we call this second component a thermal component. The number density of the thermal component is assumed to be an order of ${n}_{{\rm{GJ}}}$, as we mentioned. For case 2, we study the screening of E_∥ with kinetic time-dependent particle simulations in Section 4.

Here, we analytically describe the screening conditions in the 1D model. The momentum dispersion of the flows is neglected in this section. Note that the results do not depend on the sign of the GJ charge density ${\rho }_{{\rm{GJ}}}$. Hereafter, we denote the quantity Q due to the outflow from the NSS as ${Q}_{{\rm{ns}}}$, and that due to the backflow as ${Q}_{{\rm{bk}}}$.

Let us introduce dimensionless current densities for each component ${j}_{{\rm{k}}}$ (e.g., ${j}_{{\rm{ns}}}$ and ${j}_{{\rm{bk}}}$ in Figure 2) and the current density parameter ${j}_{{\rm{m}}}$, as

$$\begin{eqnarray}&&{\alpha }_{{\rm{k}}}=\displaystyle \frac{{j}_{{\rm{k}}}}{c{\rho }_{{\rm{GJ}}}},\quad {\alpha }_{{\rm{m}}}=\displaystyle \frac{{j}_{{\rm{m}}}}{c{\rho }_{{\rm{GJ}}}}.\end{eqnarray} \tag{ 1 }$$

We define the charge density for each component, ${\rho }_{{\rm{k}}}$. The dimensionless average velocity of a current component, ${\beta }_{{\rm{k}}}={j}_{{\rm{k}}}/{\rho }_{{\rm{k}}}c$, relates to ${\alpha }_{{\rm{k}}}$ as

$$\begin{eqnarray}&&{\beta }_{{\rm{k}}}=\displaystyle \frac{{\alpha }_{{\rm{k}}}}{({\rho }_{{\rm{k}}}/{\rho }_{{\rm{GJ}}})}.\end{eqnarray} \tag{ 2 }$$

The positive (negative) direction of ${\beta }_{{\rm{k}}}$ corresponds to the outward direction to (inward direction from) the magnetosphere.

In the 1D case, the electric field satisfies Maxwell's equations, which lead to (e.g., Levinson et al. 2005)

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {E}_{\parallel }=4\pi (\rho -{\rho }_{{\rm{GJ}}}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\parallel }}{\partial t}=4\pi c{\rho }_{{\rm{GJ}}}({\alpha }_{{\rm{m}}}-\alpha ),\end{eqnarray} \tag{ 4 }$$

in the co-rotating frame of the star, where $\rho \equiv {{\rm{\Sigma }}}_{{\rm{k}}}{\rho }_{{\rm{k}}}$ and $\alpha \equiv {{\rm{\Sigma }}}_{{\rm{k}}}{\alpha }_{{\rm{k}}}$ are the total charge and dimensionless current densities, respectively.

If the condition,

$$\begin{eqnarray}&&\rho ={\rho }_{{\rm{GJ}}},\,\,\alpha ={\alpha }_{{\rm{m}}},\,\,\end{eqnarray} \tag{ 5 }$$

is satisfied, we have a null solution, namely a steady and an uniform solution with ${E}_{\parallel }=0$, where no particle acceleration occurs. In the case that particles are supplied only from the NSS, Equations (2) and (5) require

$$\begin{eqnarray}&&{\beta }_{{\rm{ns}}}={\alpha }_{{\rm{m}}},\end{eqnarray} \tag{ 6 }$$

to screen the electric field $(j={j}_{{\rm{ns}}})$. Then, the condition $0\lt {\alpha }_{{\rm{m}}}\lt 1$ (hereafter, sub-GJ condition; Shibata 1997) should be satisfied for the null solution to be possible.⁷

If the condition $0\lt {\alpha }_{{\rm{m}}}\lt 1$ is not satisfied for the case without backflows, particles from the NSS are accelerated and the pair cascade is ignited. The condition ${\alpha }_{{\rm{m}}}\gt 1$ (hereafter super-GJ condition; Shibata 1997) requires superluminal velocity for particles from the NSS (Equation (6)), so that the velocity ${\beta }_{{\rm{ns}}}$ cannot equal to ${\alpha }_{{\rm{m}}}$, and then the system develops high voltage drops, causing γ-ray emissions due to curvature radiation and intense bursts of pair creation (Mestel et al. 1985; Shibata 1997; Beloborodov 2008) near the NSS. The system with ${\alpha }_{{\rm{m}}}\leqslant 0$ (hereafter anti-GJ condition; Timokhin & Arons 2013) also develops high voltage. A recent numerical study (Timokhin & Arons 2013) supports these speculations.

We focus on only the anti-GJ condition. Since the current in the OAR should be the return current (${\alpha }_{{\rm{m}}}\lt 0;$ Figure 1), as we mentioned in Section 2, the current at the inner boundary of the OAR is also expected to be an anti-GJ value. In this case, when the effect of the backflowing particles is negligible, a significant number of pairs is created near the NSS because of the unscreened electric field (Timokhin & Arons 2013). Here, we discuss the possibility that the backflows can assist the screening of the electric field near the NSS without pair creation.

In case 1, particles with $\gamma \gg 1$ (beam component) flow back to the NSS. The current densities of the flowing particles from the NSS ${j}_{{\rm{ns}}}$ and the OAR ${j}_{{\rm{bk,bm}}}$ are described by dimensionless current densities ${\alpha }_{{\rm{ns}}}$ and ${\alpha }_{{\rm{bk,bm}}}$ as follows:

$$\begin{eqnarray}{j}_{{\rm{ns}}} & = & {\alpha }_{{\rm{ns}}}c{\rho }_{{\rm{GJ}}},\\ {j}_{{\rm{bk,bm}}} & = & {\alpha }_{{\rm{bk,bm}}}c{\rho }_{{\rm{GJ}}}.\end{eqnarray} \tag{ 7 }$$

Note that the current density of the beam component ${\alpha }_{{\rm{bk,bm}}}$ is a parameter and is determined independently of the global current density parameter ${\alpha }_{{\rm{m}}}$. We also define the average velocity of the particles from the NSS and OAR,

$$\begin{eqnarray}{\beta }_{{\rm{ns}}} & = & \displaystyle \frac{{\alpha }_{{\rm{ns}}}}{({\rho }_{{\rm{ns}}}/{\rho }_{{\rm{GJ}}})},\\ {\beta }_{{\rm{bk,bm}}} & = & \displaystyle \frac{{\alpha }_{{\rm{bk,bm}}}}{({\rho }_{{\rm{bk,bm}}}/{\rho }_{{\rm{GJ}}})},\end{eqnarray} \tag{ 8 }$$

where ${\rho }_{{\rm{ns}}}$ and ${\rho }_{{\rm{bk,bm}}}$ represent the charge densities of the above two components. In steady state, the sum of the two current densities should satisfy

$$\begin{eqnarray}&&{\alpha }_{{\rm{ns}}}+{\alpha }_{{\rm{bk,bm}}}={\alpha }_{{\rm{m}}}.\end{eqnarray} \tag{ 9 }$$

We only consider the case ${\alpha }_{{\rm{bk,bm}}}\lt 0$ ($-1\lt {\beta }_{{\rm{bk,bm}}}\lt 0$) and ${\alpha }_{{\rm{m}}}\lt 0$, assuming that the OAR exists on the same field line. To screen the electric field uniformly, Equation (3) requires $\rho ={\rho }_{{\rm{ns}}}+{\rho }_{{\rm{bk,bm}}}={\rho }_{{\rm{GJ}}}$, which can be described as,

$$\begin{eqnarray}&&\displaystyle \frac{{\alpha }_{{\rm{ns}}}}{{\beta }_{{\rm{ns}}}^{{\rm{req}}}}+\displaystyle \frac{{\alpha }_{{\rm{bk,bm}}}}{{\beta }_{{\rm{bk,bm}}}}=1.\end{eqnarray} \tag{ 10 }$$

From this equation, the required value of ${\beta }_{{\rm{ns}}}$ to screen the field is given by,

$$\begin{eqnarray}&&{\beta }_{{\rm{ns}}}^{{\rm{req}}}=\displaystyle \frac{{\alpha }_{{\rm{ns}}}}{1-({\alpha }_{{\rm{bk,bm}}}/{\beta }_{{\rm{bk,bm}}})},\end{eqnarray} \tag{ 11 }$$

In case 1, since the current density ${j}_{{\rm{bk,bm}}}$ is carried by ultra-relativistic particles (${\beta }_{{\rm{bk,bm}}}\to -1$), Equation (11) becomes

$$\begin{eqnarray}{\beta }_{{\rm{ns}}}^{{\rm{req}}} & = & \displaystyle \frac{{\alpha }_{{\rm{ns}}}}{1+{\alpha }_{{\rm{bk,bm}}}}\\ & = & \displaystyle \frac{{\alpha }_{{\rm{m}}}-{\alpha }_{{\rm{bk,bm}}}}{1+{\alpha }_{{\rm{bk,bm}}}}.\end{eqnarray} \tag{ 12 }$$

For the flow from the NSS, we use ${\beta }_{{\rm{ns}}}^{{\rm{req}}}$ to characterize the system, including the effects of the beam component. Although the total current density is anti-GJ (${\alpha }_{{\rm{m}}}\lt 0$), ${\beta }_{{\rm{ns}}}^{{\rm{req}}}$ can have arbitrary value. The cases ${\beta }_{{\rm{ns}}}^{{\rm{req}}}\lt 0$, $0\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\lt 1$, and ${\beta }_{{\rm{ns}}}^{{\rm{req}}}\gt 1$ are similar to the situations of anti-GJ, sub-GJ and super-GJ in the system without the backflows, respectively. Hereafter, those revised versions of the conditions on the current are called such as R-anti-GJ etc., respectively. When the R-sub-GJ condition is satisfied, the electric field can be screened out by the contribution of the backflow. In this case, the OAR is steadily maintained with ${\alpha }_{{\rm{m}}}\lt 0$.

Figure 3 shows the diagrams for ${\alpha }_{{\rm{m}}}$ (a) and ${\beta }_{{\rm{ns}}}^{{\rm{req}}}$ (b) as a function of ${\alpha }_{{\rm{bk}}}(={\alpha }_{{\rm{bk,bm}}}$ in case 1) and ${\alpha }_{{\rm{ns}}}$. We focus on the cases where the number and current densities are an order of the GJ value, so that the ranges shown in Figure 3 are around ${\alpha }_{{\rm{m}}}\sim -1$ and ${\alpha }_{{\rm{bk}}}\sim -1$. Since ${\alpha }_{{\rm{m}}}$ in the steady state is the sum of the two current densities (Equation (9)), the lines corresponding to ${\alpha }_{{\rm{m}}}=\mathrm{const}$ are expressed by the dashed diagonal lines from the upper left to the lower right in this diagram. From Equation (12), the lines ${\alpha }_{{\rm{ns}}}=0,-{\alpha }_{{\rm{bk}}}-1,{\alpha }_{{\rm{bk}}}+1$ imply ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=0,-1,1$, respectively. For a certain value of ${\alpha }_{{\rm{m}}}$, there are two regions for the R-anti-GJ condition as shown in Figure 3(b) (blue; ${\beta }_{{\rm{ns}}}^{{\rm{req}}}\leqslant -1$, and $-1\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\leqslant 0$). In the red and blue regions in Figure 3(b),  significant pair creation would occur.

Zoom In Zoom Out Reset image size

These two diagrams show that even if the total current density is anti-GJ, certain ranges of the current density ${\alpha }_{{\rm{bk}}}$ can achieve R-sub-GJ. However, the electric field can be screened in only some particular combinations of the current densities ${\alpha }_{{\rm{bk}}}$ and ${\alpha }_{{\rm{m}}}$. Such an accidental combination may be rarely satisfied. In order to maintain the OAR steadily, another particle component is required unless the combination of ${\alpha }_{{\rm{bk,bm}}}$ and ${\alpha }_{{\rm{m}}}$ is adjusted in the particular range by some kind of regulation mechanism.

In case 2, the thermal component with a relatively small Lorentz factor also comes back to the NSS from the OAR [Figure 2(b)]. The thermal component is generally required to screen the electric field at the inner boundary of the OAR, as seen in the simulations in Timokhin (2010). In order to treat the particles with finite momentum distribution (including non-relativistic cases), we perform the numerical simulations of the encountering plasma just above the NSS with the 1D electrostatic Particle-in-cell code (e.g., Birdsall & Langdon 1985; Timokhin 2010; Timokhin & Arons 2013). First, we confirm the analytical results for case 1 [Figure 2(a)] discussed in Section 3. Then, we consider case 2, where both the beam and thermal components are injected from the outer boundary under the anti-GJ condition, and investigate how many thermal particles are needed to be injected from the outer boundary in order to screen the electric field.

4.1. Numerical Setup

In this subsection, we describe the numerical setup. From here, we omit the subscript ∥ from all quantities. A coordinate axis x is directed along the field lines; its origin is at the NSS and the positive direction is toward the OAR.

In the 1D model the only changing component of electromagnetic fields is the electric field component $E(x,t)$ parallel to the x-axis. We solve the evolutionary equation for the electric field $E(x,t)$ with

$$\begin{eqnarray}&&\displaystyle \frac{\partial E(x,t)}{\partial t}=-4\pi (j(x,t)-{j}_{{\rm{m}}}).\end{eqnarray} \tag{ 13 }$$

Here, $j(x,t)$ is the current density at the point x and time t. For the calculation of the current density $j(x,t)$, we use a 1D version of the charge conservative algorithm proposed by Villasenor & Buneman (1992). In order to obtain the initial value of the electric field $E(x,0)$, we solve the Poisson equation

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dx}}(x,0)=4\pi (\rho (x,0)-{\rho }_{{\rm{GJ}}}).\end{eqnarray} \tag{ 14 }$$

For the GJ charge density ${\rho }_{{\rm{GJ}}}$, we assume a spatially and temporary constant value in the calculation domain.

To model the free emission of the particles from the NSS (x = 0), we adopt the same method proposed by Timokhin & Arons (2013). At the beginning of each time-step, electrons and ions with a certain equal number ${N}_{{\rm{ns}}}$ are injected into a cell just outside the numerical domain $x\lt 0$ (ghost cell) to carry out the electric field calculation. We adopt a significant number of ${N}_{{\rm{ns}}}$ to correctly simulate the space-charge limitation condition (for details, see Appendix C in Timokhin & Arons 2013). The momentum of each injected particle is sampled from an uniform distribution in the interval $[-{p}_{{\rm{ns}}},{p}_{{\rm{ns}}}]$ to model non-relativistic and finite temperature of the particles. We adopt the momentum ${p}_{{\rm{ns}}}={10}^{-2}{m}_{{\rm{e}}}c$, though the results do not basically depend on the specific value of ${p}_{{\rm{ns}}}(\lt {m}_{{\rm{e}}}c)$. We take the mass ratio of ion to electron as ${m}_{{\rm{i}}}/{m}_{{\rm{e}}}=1836$.

We introduce ${\alpha }_{{\rm{bk,th}}}$ as the current density of the thermal component. Then, the current density of backflowing particles is described by

$$\begin{eqnarray}&&{\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{bk,bm}}}+{\alpha }_{{\rm{bk,th}}}.\end{eqnarray} \tag{ 15 }$$

In the steady state, the current densities should satisfy the relation

$$\begin{eqnarray}&&{\alpha }_{{\rm{ns}}}+{\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{m}}}.\end{eqnarray} \tag{ 16 }$$

We inject the backflowing particles just outside the numerical domain $x\gt 1$ (the last cell) in our simulations, where x is normalized by the length of the calculation domain L. For the beam component, we inject electrons (for ${\rho }_{{\rm{GJ}}}\lt 0$) or positrons (for ${\rho }_{{\rm{GJ}}}\gt 0$) with Lorentz factor $\gamma ={10}^{6}$ into the calculation domain. The number density of the injected particles is determined by the parameter ${\alpha }_{{\rm{bk,bm}}}$. The particles are injected with a time-step $\delta t$, and the injection point for each particle is distributed randomly in space into the last cell. The thermal particles from the OAR consist of electrons and positrons with the same number density. At the injection point, the average velocity of the thermal component is ${\bar{\beta }}_{{\rm{bk,th}}}=0$ and the temperature is assumed as ${{kT}}_{{\rm{bk,th}}}={m}_{{\rm{e}}}{c}^{2}$. In this setup, the charge and current densities are ${\rho }_{{\rm{bk,th}}}=0$ and ${\alpha }_{{\rm{bk,th}}}=0$ at the injection point. The model parameter for the thermal component is only the number density ${n}_{{\rm{bk,th}}}$ in the last cell.

In the 1D model, charged particles are represented by a thin sheath with an infinite extend in the direction perpendicular to the x-axis. The equation of motion for a particle i is

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}_{i}}{{dt}}={v}_{i}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dp}}_{i}}{{dt}}=\displaystyle \frac{{eE}({x}_{i})}{{mc}},i=1,\ldots ,\,{N}_{{\rm{p}}},\end{eqnarray} \tag{ 18 }$$

where x_i and p_i are the position and momentum of the ith particle in units of ${\lambda }_{{\rm{pe}}}$ and mc, respectively. The length scale ${\lambda }_{{\rm{pe}}}$ is the electron skin depth for the GJ density,

$$\begin{eqnarray}&&{\lambda }_{{\rm{pe}}}\equiv \displaystyle \frac{c}{{\omega }_{{\rm{pe}}}}=c{\left(\displaystyle \frac{4\pi | {\rho }_{{\rm{GJ}}}| e}{{m}_{{\rm{e}}}}\right)}^{-1/2}.\end{eqnarray} \tag{ 19 }$$

Since we focus on the nearly screened state of the electric field, particles do not obtain high momentum (at most $\sim {10}^{6}{m}_{{\rm{e}}}c$ even for the beam particles) in our calculation domain within the reasonable range of ${\alpha }_{{\rm{m}}}$ (an order of unity). The radiation drag force starts to work at $\gamma \gt {10}^{6}$, so we can neglect the reaction force in the equation of motion (18) (Timokhin 2010). We have to choose the total number of the simulation particles ${N}_{{\rm{p}}}$ and the length of the calculation domain L in order to satisfy the conditions

$$\begin{eqnarray}&&L\gg {\lambda }_{{\rm{pe}}},{N}_{{\rm{p}}}\gg \displaystyle \frac{L}{{\lambda }_{{\rm{pe}}}}.\end{eqnarray} \tag{ 20 }$$

In this limit, the results are expected to be independent of the choice of ${N}_{{\rm{p}}}$ and L. In our simulations, typically $L=({10}^{2}$–${10}^{3}){\lambda }_{{\rm{pe}}}$ and ${N}_{{\rm{p}}}\sim {10}^{6}$. Another requirement is a significantly short time-step in the simulation, $\delta t\ll {\omega }_{{\rm{pe}}}$, to resolve plasma oscillations. In our simulations, non-relativistic particles play an important role in screening the electric field. Our choices for the time-step $\delta t$ and each grid size $\delta x$ are always $\delta t=0.02{\omega }_{{\rm{pe}}}$ and $\delta x=0.1{\lambda }_{{\rm{pe}}}$.

Initially, we set a condition that there is no particle in the calculation domain ($\rho (x,0)=0$). We have checked that the results basically do not depend on the initial particle distribution once the system reaches the quasi-steady state.

The extracted particles (electrons or ions) from the NSS in the quasi-steady state depend on the signs of the parameters, ${\alpha }_{{\rm{m}}}-{\alpha }_{{\rm{bk,bm}}}$ and ${\rho }_{{\rm{GJ}}}$.⁸ Equation (9) indicates that the signs of ${\alpha }_{{\rm{ns}}}$ and ${\alpha }_{{\rm{m}}}-{\alpha }_{{\rm{bk,bm}}}$ are the same in the quasi-steady state. Since ${\beta }_{{\rm{ns}}}$ is always positive, the sign of ${\rho }_{{\rm{ns}}}$ is the same as (opposite to) that of ${\rho }_{{\rm{GJ}}}$ for positive (negative) ${\alpha }_{{\rm{ns}}}$ (see Equation (8)) or equivalently positive (negative) ${\alpha }_{{\rm{m}}}-{\alpha }_{{\rm{bk,bm}}}$. The parameter sets we adopt are summarized in Table 1.

Table 1.  Simulation Model Parameters

Model	${\rho }_{{\rm{GJ}}}$	${\alpha }_{{\rm{m}}}$	${\alpha }_{{\rm{bk,bm}}}$	${\beta }_{{\rm{ns}}}^{{\rm{req}}}$	Condition	${n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}$
A	negative	−0.4	−0.6	0.5	R-sub-GJ	0
B	negative	−0.4	−0.8	2.0	R-super-GJ	0
C1	negative	−0.4	−0.8	2.0	R-super-GJ	0.1
C2	negative	−0.4	−0.8	2.0	R-super-GJ	0.2
C3	negative	−0.4	−0.8	2.0	R-super-GJ	0.4
C4	negative	−0.4	−0.8	2.0	R-super-GJ	0.6
C5	negative	−0.4	−0.8	2.0	R-super-GJ	0.8
C6	negative	−0.4	−0.8	2.0	R-super-GJ	1.0
D1	negative	−0.4	−1.2	−4.0	R-anti-GJ	0.2
D2	negative	−0.4	−1.2	−4.0	R-anti-GJ	0.6
D3	negative	−0.4	−1.2	−4.0	R-anti-GJ	1.0
D4	negative	−0.4	−1.2	−4.0	R-anti-GJ	1.4
D5	negative	−0.4	−1.2	−4.0	R-anti-GJ	1.8
D6	negative	−0.4	−1.2	−4.0	R-anti-GJ	2.2
E1	positive	−0.4	−0.2	−0.25	R-anti-GJ	0.5
E2	positive	−0.4	−0.2	−0.25	R-anti-GJ	0.6
E3	positive	−0.4	−0.2	−0.25	R-anti-GJ	0.7
E4	positive	−0.4	−0.2	−0.25	R-anti-GJ	0.8
E5	positive	−0.4	−0.2	−0.25	R-anti-GJ	0.9
E6	positive	−0.4	−0.2	−0.25	R-anti-GJ	1.0
F1	positive	−0.4	−0.8	2.0	R-super-GJ	0.2
F2(a–e)	positive	−0.4	−0.8	2.0	R-super-GJ	0.4
F3	positive	−0.4	−0.8	2.0	R-super-GJ	0.6
F4	positive	−0.4	−0.8	2.0	R-super-GJ	1.1
F5	positive	−0.4	−0.8	2.0	R-super-GJ	1.6

Download table as:  ASCIITypeset image

Although we only consider ${\alpha }_{{\rm{m}}}=-0.4$ in Models A–F, the following results are the same as the case ${\alpha }_{{\rm{m}}}=-1.6$. The diagram in Figure 3(b) is symmetrical to the point $({\alpha }_{{\rm{bk}}},{\alpha }_{{\rm{ns}}})=(-1,0)$. For example, electrons with the average density ${\bar{n}}_{{\rm{ns}}}/{n}_{{\rm{GJ}}}=0.4$ and velocity ${\bar{\beta }}_{{\rm{ns}}}=0.5$ are extracted from the NSS in Model A. In the anti-pulsar case (${\rho }_{{\rm{GJ}}}\gt 0$) with the current densities ${\alpha }_{{\rm{m}}}=-1.6$ and ${\alpha }_{{\rm{bk,bm}}}=-1.4$, the electrons with average density ${\bar{n}}_{{\rm{ns}}}/{n}_{{\rm{GJ}}}=0.4$ and velocity ${\bar{\beta }}_{{\rm{ns}}}=0.5$ are also extracted from the NSS to screen the electric field from Equation (12). Then, the two cases, $({\alpha }_{{\rm{m}}},{\alpha }_{{\rm{bk,bm}}},{\rho }_{{\rm{GJ}}})$ and $({\alpha }_{{\rm{m}}}^{\prime },{\alpha }_{{\rm{bk,bm}}}^{\prime },{\rho }_{{\rm{GJ}}}^{\prime })$ = $(-{\alpha }_{{\rm{m}}}-2,-{\alpha }_{{\rm{bk,bm}}}-2,-{\rho }_{{\rm{GJ}}})$, result in the same solution for the particle flow from the NSS.

4.2. Single Beam Backflow

In the left panels of Figure 4, we show a snapshot for different physical quantities in the quasi-steady state of Model A (${\alpha }_{{\rm{m}}}=-0.4,{\alpha }_{{\rm{bk,bm}}}=-0.6$), for which ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=0.5$ satisfies the R-sub-GJ condition. In the case ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=0.5$, the extracted electrons may behave similarly to the sub-GJ case of ${\alpha }_{{\rm{m}}}=0.5$ without the beam component (Timokhin & Arons 2013; Chen & Beloborodov 2013). The figure shows that an electric field is screened in most of the calculation domain as we expected. In the momentum distribution of the extracted electrons, two components, the cold outgoing beam and the trapped thermal beams, are seen. The cold beam particles from the NSS adjust the total current density to ${\alpha }_{{\rm{m}}}$. Because particles from the NSS initially have a non-relativistic velocity, the absolute value of the normalized charge density $| \rho /{\rho }_{{\rm{GJ}}}|$ just above the surface is larger than unity in the steady state. Then, the non-relativistic particles are accelerated by the generating electric field and the normalized charge density $| \rho /{\rho }_{{\rm{GJ}}}|$ decreases toward the OAR. The accelerated particles are mainly composed of the outgoing beam component. The trapped component in the calculation domain adjusts the total charge density to ${\rho }_{{\rm{GJ}}}$, except for the region just above the NSS. The formation of this trapped thermal component is attributed to the instability as discussed by previous authors (Chen & Beloborodov 2013; Timokhin & Arons 2013), who show that the spatially oscillated solution in cold limits on the particle momentum (e.g., Shibata 1991) is unstable for the particles with non-negligible temperatures at the NSS. The trapped particles are reflected by the electric field near the boundaries, so most of them cannot take part in the outgoing beam component. The electric field structure within $0\lt x\lt 0.1$ can be regarded as the relativistic double layer, which is defined as a large potential drop ($\gg {m}_{{\rm{e}}}{c}^{2}$) maintained by an anode and a cathode neglecting gravity and thermal plasma motion at the boundaries (e.g., Carlqvist 1982). The electric field also appears near the outer boundary ($0.9\lt x\lt 1.0$) to trap the thermal particles. The results for the extracted particles from the NSS are close to the sub-GJ case without backflowing component (Chen & Beloborodov 2013; Timokhin & Arons 2013), so even if ${\alpha }_{{\rm{m}}}\lt 0$, the electric field is screened and no significant particle acceleration occurs when $0\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\lt 1$. The tiny momentum fluctuation of the beam component from the outer boundary (seen in the left second column of Figure 4) is due to the fluctuation in the charge density of the trapped particles. For the beam component from the OAR, the velocity (∼−c) and the current density are almost constant, as shown in the left panels of Figure 4.

Zoom In Zoom Out Reset image size

In the right panels of Figure 4, we also show the quasi-steady state for the R-super-GJ case (Model B). We set the sign of ${\rho }_{{\rm{GJ}}}$ as positive and inject positrons as beam particles from the outer boundary. Then, the extracted particles from the NSS are mainly electrons in Model B. Obviously, particles from both boundaries are continuously accelerated. Because our calculation domain is not so large, the change of Lorentz factor is $\sim {10}^{3}$. In reality, the value of the Lorentz factor reaches the acceleration-reaction limit. As a result, electrons and positrons emit γ-ray photons and significant pair creation should occur. Therefore, we confirm that the beam component from the OAR does not assist with the screening for the R-super-GJ case, as mentioned in Section 3. In order to screen the field without pair creation in the anti-GJ case, the R-sub-GJ condition is required, as long as only an ultra-relativistic component is considered as the backflow.

4.3. Beam and Thermal Backflows

In this subsection, we assume that both the beam and thermal components come back to the NSS under the anti-GJ condition (${\alpha }_{{\rm{m}}}=-0.4$). The thermal particles may coherently interact with the particles extracted from the NSS. Then, the results depend on the relative mass of the particles. We separately consider the electron-extracted case (Section 4.3.1) and the ion-extracted case (Section 4.3.2) from the NSS. Hereafter, we classify cases as R-anti-GJ, etc., based on Equation (12). As we will see later, if the number density of thermal particles is large enough, their contributions finally change the status of the whole calculation domain to R-sub-GJ.

4.3.1. Electron-extracted Case

For a given parameter ${\alpha }_{{\rm{m}}}$, a dotted line in Figure 3(b) shows that the system can be the R-super-GJ, R-sub-GJ, and two R-anti-GJ regions (${\alpha }_{{\rm{bk}}}\lt -1,{\alpha }_{{\rm{bk}}}\gt -1$) for different ${\alpha }_{{\rm{bk}}}$. We investigate the screening conditions for the R-super-GJ condition (Models C1-6; ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=2.0$) and for the two kinds of R-anti-GJ conditions (Models D1-6; ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=-4.0$, and Models E1-6; ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=-0.25$). The length of the calculation domain is $L=100{\lambda }_{{\rm{pe}}}$. Since the beam backflow is not significantly decelerated in our results, the charge and current densities of the beam component are constant in the calculation domain as seen in Section 4.2. We omit showing the momentum portrait for the beam component from the OAR in Section 4.3.

In Figure 5, we present snapshots of the quasi-steady state, starting from the R-super-GJ condition (Models C1-6; ${\alpha }_{{\rm{bk,bm}}}=-0.8$ and ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=2.0$) for different number densities of the thermal electrons and positrons at the outer boundary. The parameters ${\alpha }_{{\rm{m}}}$ and ${\alpha }_{{\rm{bk,bm}}}$ are the same as in Model B. The sign of ${\rho }_{{\rm{GJ}}}$ is negative, so the beam particles from the OAR are electrons. In Model C1 $({n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.1;$ the first row in Figure 5), the electric field is not screened. The value of the current density $-0.7\lt {\alpha }_{{\rm{bk}}}\lt -0.4$ is required to achieve the R-sub-GJ [Figure 3(b)]. Even if all the thermal positrons initially have a velocity ${\beta }_{{\rm{bk,th}}}\sim -1$ and act like the beam particles, the total current density of the flow from the OAR becomes ${\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{bk,bm}}}+{{en}}_{{\rm{bk,th}}}\times (-c)/({\rho }_{{\rm{GJ}}}c)=-0.7$ (${\beta }_{{\rm{ns}}}^{{\rm{req}}}=1$), which is just at the boundary of the required condition. Considering the slow average velocity $| {\bar{\beta }}_{{\rm{bk,th}}}| \ll 1$ at the outer boundary, the current density becomes ${\alpha }_{{\rm{bk}}}\lt -0.7$ due to the continuity equation. Furthermore, if the electric field is screened at the outer boundary, only half the number density of the thermal positrons $0.5{n}_{{\rm{bk,th}}}$ with an initial momentum ${p}_{{\rm{tp}}}\lt 0$ can enter the calculation domain. Thus, the injected number density of the thermal particles ${n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.1$ is too small to steadily screen the electric field. In order to clarify this result, we plot the value of the current densities ${\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{bk,bm}}}+{\alpha }_{{\rm{bk,th}}}$ and ${\alpha }_{{\rm{ns}}}$ in Figure 6 for each position x in the calculation domain on the same diagram as Figure 3(b). The dashed line denotes the condition ${\alpha }_{{\rm{ns}}}+{\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{m}}}=-0.4$. In Model C1 (red dots), the dots indicate the R-super-GJ region (red region). When the currents achieve the quasi-steady state, the dots should be along the dashed line. Actually the results show ${\alpha }_{{\rm{bk}}}+{\alpha }_{{\rm{ns}}}={\alpha }_{{\rm{m}}}$ on average and the system can be considered as the quasi-steady state ($\partial E/\partial t\sim 0$ from Equation (13)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Model C2 $({n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.2;$ the second row in Figure 5), the accelerating electric field still exists over the calculation domain. Because of the electric field near the outer boundary, the backflowing thermal positrons are accelerated toward the NSS (the third column in Figure 5). Most electrons of the thermal component do not contribute to the current density in the calculation domain (the second column in Figure 5). As shown in Figure 6, the center of the distribution of the current densities (green dots) indicates ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=1$ (solid diagonal line). Since particles cannot reach $| v| =c$ strictly, this solution still requires the electrons from the NSS to be accelerated continuously (e.g., Michel 1974).

When ${n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}\gt 0.2$, the electric field is almost screened over the calculation domain except for the regions near the inner and outer boundaries. In Figure 5, no significant particle acceleration is seen in Models C3-6 (${n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.4,0.6,0.8$ and 1.0, respectively). The momentum distributions of the thermal component from the OAR are mainly determined by the initial temperature ${{kT}}_{{\rm{bk,th}}}={m}_{{\rm{e}}}{c}^{2}$. Note that in contrast to the particles from the NSS, both thermal electrons and positrons from the OAR enter the calculation domain. In Model C3 (the third row in Figure 5), although the electric field near the outer boundary accelerates positrons up to ∼−c, the momentum of the accelerated particle is an order of $\sim -{m}_{{\rm{e}}}c$ in the steady state. Since the initial temperature at the outer boundary is comparable to the electron rest-mass energy, parts of them have a mildly relativistic velocity $\gamma \sim 3\mbox{--}5$. Such mildly relativistic electrons are not reflected by the electric field near the outer boundary and penetrate the calculation domain. On the other hand, at the inner boundary, injected particles have non-relativistic velocity, so only electrons accelerated by the electric field near the inner boundary can enter the calculation domain.

In Models C5 and C6, the thermal particles from the OAR largely contribute to the total charge and current densities (the fifth and sixth columns of Figure 5). These contributions adjust the the total current density to the stationary condition $\alpha ={\alpha }_{{\rm{m}}}$. As shown in Figure 6, almost all dots for ${n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}\geqslant 0.4$ indicate the R-sub-GJ region, except for a small number of dots in the R-super-GJ region. The electrons from the NSS contribute to the total current density (${\alpha }_{{\rm{ns}}}\ne 0;$ see the blue lines in the fifth column of Figure 5). Such electrons have to be accelerated to significantly contribute to the total current density. The electric field structure near the NSS can be regarded as the relativistic double layer, which is similar to the result in the R-sub-GJ case without the thermal components from the OAR (Figure 4). The dispersion of the dots in the R-sub-GJ region in Figure 6 is caused by the plasma oscillation of trapped particles in the calculation domain, which adjusts the charge density to the GJ value. The contribution of ${\alpha }_{{\rm{ns}}}$ to the total current density becomes small as ${n}_{{\rm{bk,th}}}$ increases. Especially for Model C6, ${\alpha }_{{\rm{ns}}}\sim 0$ in the calculation domain. In this case, only the thermal component from the OAR adjusts both the current and charge densities and the structure of the double layer near the NSS is no longer sustained. We also perform the simulations with a larger injection number density ${n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}\gt 1.0$ and confirm that only the thermal component adjusts the current and charge densities to $\alpha ={\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{m}}}$ and $\rho ={\rho }_{{\rm{GJ}}}$ in the R-sub-GJ region. This is because the thermal electrons from the OAR screen the electric field instead of the electrons from the NSS. This practically indicates elimination of the polar cap in the anti-GJ region; the entire region except for the OAR is steadily filled with a high-density plasma that screens the electric field.

Let us move on to the case starting from the R-anti-GJ condition. First, we investigate the case ${\beta }_{{\rm{ns}}}^{{\rm{req}}}\leqslant -1$ (Models D1-6; ${\alpha }_{{\rm{m}}}=-0.4$, ${\alpha }_{{\rm{bk,bm}}}=-1.2$, so that ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=-4.0$). In Figure 7, we present the current densities (${\alpha }_{{\rm{bk}}},{\alpha }_{{\rm{ns}}}$) at each position in the quasi-steady state for Models D1-6. Similar to the results in the case starting from the R-super-GJ condition (Models C1-6), the distribution of dots shifts toward the right side along the line of ${\alpha }_{{\rm{ns}}}+{\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{m}}}$ as ${n}_{{\rm{bk,th}}}$ increases. In Model D1 $({n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.2;$ red dots in Figure 7), the current density is ${\alpha }_{{\rm{bk}}}\sim -1.0$. This means that almost all injected thermal positrons are accelerated up to ∼−c due to the large electric field near the outer boundary. Because the average velocity of particles that enter the calculation domain is $| {\bar{\beta }}_{{\rm{bk,th}}}| \ll 1$ at the outer boundary and the continuity equation is satisfied, the current density ${\alpha }_{{\rm{bk}}}$ is a bit smaller than −1.0. As the density ${n}_{{\rm{bk,th}}}$ increases, the contribution of ${\alpha }_{{\rm{ns}}}$ decreases. The electric field is almost screened out when the dots are on the green region in Figure 7. Those results are almost the same as the case starting from the R-super-GJ condition (Models C1-6; Figure 5).

Zoom In Zoom Out Reset image size

Next, we consider another R-anti-GJ condition (Models E1-6; ${\alpha }_{{\rm{bk,bm}}}=-0.2$, ${\alpha }_{{\rm{m}}}=-0.4$ and ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=-0.25$), in which the sign of ${\alpha }_{{\rm{m}}}-{\alpha }_{{\rm{bk,bm}}}$ is negative. To investigate the electron-extracted case, we assume that the sign of ${\rho }_{{\rm{GJ}}}$ is positive and that the beam particles from the OAR are positrons. The thermal component from the OAR succeeds in decreasing the current density ${\alpha }_{{\rm{bk}}}$ down to −0.4. We show the results in the current density diagram of Figure 8. All the cases for different values of ${n}_{{\rm{bk,th}}}$ enter the R-sub-GJ region in the diagram, and the electric field is screened in the calculation domain. In the cases with ${\beta }_{{\rm{ns}}}^{{\rm{req}}}\lt -1$, there are steady-state solutions with ${\alpha }_{{\rm{ns}}}\ne 0$ in the R-sub-GJ region (Models D4 and D5). However, in the case $-1\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\lt 0$ (Models E1-6), the current density ${\alpha }_{{\rm{ns}}}$ is always ∼0 in the screened solutions. This tendency is irrespective of whether the extracted particles are ions or electrons.

Zoom In Zoom Out Reset image size

Our numerical results shown above provide the condition for screening the electric field. In the electron-extracted case, we divide non-R-sub-GJ regions into $| {\beta }_{{\rm{ns}}}^{{\rm{req}}}| \geqslant 1$ and $-1\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\lt 0$. In the former case, the injection of the thermal component tends to reduce $| {\alpha }_{{\rm{ns}}}|$, keeping ${\alpha }_{{\rm{ns}}}+{\alpha }_{{\rm{bk}}}={\alpha }_{{\rm{m}}}$. If we inject the thermal component enough to achieve ${\beta }_{{\rm{ns}}}^{{\rm{req}}}\lt 1$, the electric field is screened over the calculation domain. When the system enter the R-sub-GJ region with ${\alpha }_{{\rm{ns}}}\ne 0$, particles are extracted from the NSS. In order to derive the screening condition, we introduce a parameter, n₀, which is the density of particles entering the calculation domain without the electric field at the outer boundary (i.e., ${n}_{0}=0.5{n}_{{\rm{bk,th}}}$ for ${\bar{\beta }}_{{\rm{bk,th}}}=0$). If electrons or positrons of the thermal component are accelerated up to ∼−c in the calculation domain, the density n₀ is described as ${n}_{0}/{n}_{{\rm{GJ}}}=| {\alpha }_{{\rm{bk}}}-{\alpha }_{{\rm{bk,bm}}}|$. From Equation (12), the screening condition is

$$\begin{eqnarray}&&\displaystyle \frac{{\alpha }_{{\rm{m}}}-{\alpha }_{{\rm{bk}}}}{1+{\alpha }_{{\rm{bk}}}}\lt 1.\end{eqnarray} \tag{ 21 }$$

Then, we derive the required density of the thermal component n₀ as,

$$\begin{eqnarray}\displaystyle \frac{{n}_{0}}{{n}_{{\rm{GJ}}}}\gt \left\{\begin{array}{ll}\tfrac{1}{2}| ({\alpha }_{{\rm{bk,bm}}}+1)({\beta }_{{\rm{ns}}}^{{\rm{req}}}-1)| & \,{\rm{for}}\,\,{\alpha }_{{\rm{bk,bm}}}\ne -1,\\ \tfrac{1}{2}| {\alpha }_{{\rm{m}}}+1| & \,{\rm{for}}\,\,{\alpha }_{{\rm{bk,bm}}}=-1.\end{array}\right.\end{eqnarray} \tag{ 22 }$$

Note that ${\beta }_{{\rm{ns}}}^{{\rm{req}}}$ is defined in Equation (12), neglecting ${\alpha }_{{\rm{bk,th}}}$. On the other hand, for the case $-1\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\leqslant 0$, the screened state always implies ${\alpha }_{{\rm{ns}}}=0$ (${\alpha }_{{\rm{m}}}={\alpha }_{{\rm{bk}}}$). From Equation (12), the screening condition is

$$\begin{eqnarray}&&{\alpha }_{{\rm{m}}}\gt {\alpha }_{{\rm{bk}}},\end{eqnarray} \tag{ 23 }$$

so we derive the required density n₀ to screen the electric field as,

$$\begin{eqnarray}&&{n}_{0}\gt | {\beta }_{{\rm{ns}}}^{{\rm{req}}}({\alpha }_{{\rm{bk,bm}}}+1)| .\end{eqnarray} \tag{ 24 }$$

Note that this equation also has a singular point at ${\alpha }_{{\rm{bk,bm}}}=-1$. However, ${\alpha }_{{\rm{bk,bm}}}=-1$ corresponds to ${\beta }_{{\rm{ns}}}^{{\rm{req}}}=\pm \infty$, which is not included in the case $-1\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\leqslant 0$. The results do not depend on L and ${N}_{{\rm{p}}}$ if the conditions (20) are satisfied.

4.3.2. Ion-extracted Case

Here, we investigate the case in which extracted particles from the NSS are ions. We start from the R-super-GJ condition, in which the parameters of the current density are the same as  those for Models B and C1-6 (Models F1-5; ${\alpha }_{{\rm{bk,bm}}}=-0.8$ and ${\alpha }_{{\rm{m}}}=-0.4$). To consider the ion-extracted case, we assume that the GJ charge density ${\rho }_{{\rm{GJ}}}$ is positive. Then, the beam particles from the OAR are positrons. In the electron-extracted cases (Section 4.3.1), we set the length of the calculation domain $L=100{\lambda }_{{\rm{pe}}}$. Because the ion skin depth is large, ${\lambda }_{{\rm{pp}}}\sim 43{\lambda }_{{\rm{pe}}}$, we have to take a larger calculation domain to correctly simulate the dynamics including electrons and ions. According to the limited calculation time, we set $L={10}^{3}{\lambda }_{{\rm{pe}}}$ and the grid size is the same as the previous calculations.

Although we confirm that the screening condition of the electric field (Equation (22)) also works in the ion-extracted case, in some cases we obtain the quasi-periodic solutions, which are not seen in the electron-extracted cases. The typical cycles of the quasi-periodic behavior for Model F2 $({n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.4)$ are shown in Figure 9. Initially (see the first row) the thermal component from the outer boundary already reaches the NSS and screens the electric field over the entire calculation domain. During the electric field being screened out, newly extracted ions are not significantly accelerated ($\ll c$) and the current density ${\alpha }_{{\rm{ns}}}$ is suppressed (see the second row). In order to maintain the condition $\alpha ={\alpha }_{{\rm{m}}}$, the electric field near the NSS starts to increase and ions are extracted from the NSS (see the third row). However, since the ion inertia is larger compared to the electron inertia, ions are not immediately accelerated and the total current density α does not reach ${\alpha }_{{\rm{m}}}$ at that time. Then, the electric field is continuously grown, due to the difference $\alpha -{\alpha }_{{\rm{m}}}$. Although the thermal electrons are accelerated toward the NSS and positrons are reflected by the electric field, the electron/positron number density at the NSS is not enough to screen the electric field. Subsequently, the electric field significantly develops at the NSS (see the fourth row). When the electric field develops enough to accelerate ions up to relativistic velocity, accelerated ions catch up to slowly moving ions ejected earlier and a bunch of ions (i.e., over-density region) is generated (see the fifth row). This ion bunch reflects positrons and traps electrons in order to adjust the charge density $\rho \to {\rho }_{{\rm{GJ}}}$. This process increases the number density of the particles in the bunch. We show the evolution of the number density distribution for the thermally injected electrons/positrons in Figure 10. Finally, the bunch escapes from the calculation domain (see the sixth row in Figure 9) and again the thermal component starts to screen the electric field in the calculation domain. This sequence of the instability occurs quasi-periodically, as shown in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The quasi-periodic solutions are obtained only in limited parameter ranges. Figure 11 shows the temporal evolutions of the density of particles except for the beam component from the OAR at the outer boundary in Models F1-5. The quasi-periodic behavior is only seen for Models F2 $({n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.4;$ blue) and F3 $({n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}=0.6;$ magenta). Since Model F1 $({n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}\leqslant 0.2;$ red) is not screened, the electric field over the calculation domain, the thermal component develops to the beam particles. This solution is almost the same one in Model C2. If number densities of the thermal particles that are too large are injected at the outer boundary (${n}_{{\rm{bk,th}}}/{n}_{{\rm{GJ}}}\gtrsim 1.0$, Models F4 and F5), ${\alpha }_{{\rm{ns}}}\to 0$ and ions are not extracted from the NSS. Then, the solution becomes quasi-steady. Therefore, the quasi-periodic solution is obtained in the cases with ${\alpha }_{{\rm{ns}}}\ne 0$. This means that the cases with the R-sub-GJ and R-anti-GJ conditions possibly have the quasi-periodic solutions in some ranges of the density of the thermal component. On the other hand, the R-anti-GJ case with $-1\lt {\beta }_{{\rm{ns}}}^{{\rm{req}}}\leqslant 0$ is not expected to have the quasi-periodic solution, because the screened state of the electric field is always ${\alpha }_{{\rm{ns}}}=0$, as shown in the electron-extracted case (Figure 8).

Zoom In Zoom Out Reset image size

Although the qualitative results in the quasi-steady state (e.g., Models F4 and F5) do not depend on the length of the calculation domain L, the period of the sequence and the growth of the bunch may depend on L. We perform simulations with different domain lengths L. In all the cases, we fix the current densities and the injected number density as Model F2. We label the models F2(a) ($L=100{\lambda }_{{\rm{pe}}}$), F2(b) ($L=500{\lambda }_{{\rm{pe}}}$), F2(c) ($L=1000{\lambda }_{{\rm{pe}}}$), F2(d) ($L=2000{\lambda }_{{\rm{pe}}}$), and F2(e) ($L=3000{\lambda }_{{\rm{pe}}}$). Figure 12 shows the temporal evolution of the number density of particles, except for the beam component from the OAR at the outer boundary for various $L/{\lambda }_{{\rm{pe}}}$. Except for Model F2(a), the peak number density is about an order of magnitude larger than the number density of the thermal component at the injection. The period of the sequence ${P}_{{\rm{s}}}$ clearly depends on L. In Figure 13 we plot the relation between ${P}_{{\rm{s}}}$ and L. For Models F2(a–e), we obtain the relation between ${P}_{{\rm{s}}}$ and L as

$$\begin{eqnarray}&&{P}_{{\rm{s}}}\sim 3.9L/c+6.0\times {10}^{2}{\omega }_{{\rm{pe}}}^{-1}.\end{eqnarray} \tag{ 25 }$$

The timescale for the propagation of a bunch of particles is almost a few times the light crossing time. The intercept in Equation (25) is almost consistent with the acceleration time of ions, $t{\omega }_{{\rm{pe}}}\sim ({m}_{{\rm{i}}}/{m}_{{\rm{e}}})\sim {10}^{3}$ for $E/({m}_{{\rm{e}}}c{\omega }_{{\rm{pe}}}/e)\sim 2$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Coherent curvature radiation from a bunch of particles has been considered to be a possible mechanism of the observed pulsar radio emission (e.g., Ruderman & Sutherland 1975). We estimate the luminosity of the radio emission based on our results in the ion-extracted case, where the bunch of particles is formed quasi-periodically. From the simulation results, the size of a bunch is an order of $\sim 10{\lambda }_{{\rm{pe}}}\sim 20\,{\rm{cm}}$ and the number density of the bunching particles is several times the GJ number density. The period of the quasi-periodic solution is an order of the light crossing time of the calculation domain.

In our results in the previous section, the Lorentz factor of the bunch particles is up to $\gamma \sim 10$. However, this value should depend on the number density of extracted ions. In some parameter sets we performed, the bunching particles are heated up to $\gamma \sim {10}^{2}$. The generating electric field near the NSS accelerates ions up to relativistic velocity in order to adjust the current density. Then, the maximum Lorentz factor of the electron/positron, which is determined by the potential energy near the NSS, can be an order of the mass ratio $\sim {m}_{{\rm{i}}}/{m}_{{\rm{e}}}\sim {10}^{3}$.

The number of particles in a bunch is

$$\begin{eqnarray}&&N\sim \kappa {n}_{{\rm{GJ}}}\times {(10{\lambda }_{{\rm{e}}})}^{3},\end{eqnarray} \tag{ 26 }$$

where κ is the average density in the bunch normalized by ${n}_{{\rm{GJ}}}$. We assume that multiple bunches compose a thin shell whose surface area is comparable to the polar cap surface. In this case, the number of bunches in a shell is ${N}_{{\rm{bunch}}}\,\sim 10{\lambda }_{{\rm{e}}}\times \pi {r}_{{\rm{pc}}}^{2}/{(10{\lambda }_{{\rm{e}}})}^{3}$. We also assume ${P}_{{\rm{s}}}\sim {r}_{{\rm{pc}}}/c$ so that the number of shells we can observe in one rotational period is ${N}_{{\rm{shell}}}\sim {\rm{\Delta }}P/{P}_{{\rm{s}}}$, where ${\rm{\Delta }}P$ is the observed pulse width and we use ${\rm{\Delta }}P\sim {10}^{-2}{P}^{1/2}$ based on the observational results in Maciesiak & Gil (2011). Then, the coherent radio luminosity is estimated as

$$\begin{eqnarray}{L}_{{\rm{radio}}} & \sim & {\dot{E}}_{{\rm{cur}}}{N}^{2}{N}_{{\rm{bunch}}}{N}_{{\rm{shell}}}\\ & \sim & {10}^{28}{\left(\displaystyle \frac{\kappa }{10}\right)}^{2}{\left(\displaystyle \frac{\gamma }{{10}^{2}}\right)}^{4}{\left(\displaystyle \frac{{R}_{{\rm{cur}}}}{{10}^{7}{\rm{cm}}}\right)}^{-2}\ \mathrm{erg}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 27 }$$

where ${\dot{E}}_{{\rm{cur}}}$ is the power of curvature radiation for single particles and ${R}_{{\rm{cur}}}$ is the radius of the field line curvature. Note that the radio luminosity in Equation (27) does not explicitly depend on the spin-down luminosity. The observations also show that the radio luminosity is $\sim {10}^{27}\mbox{--}{10}^{31}$ erg s⁻¹, regardless of the position in the $P\mbox{--}\dot{P}$ diagram (Szary et al. 2014). The characteristic frequency is

$$\begin{eqnarray}&&{\nu }_{{\rm{c}}}\sim {10}^{9}\ {\rm{Hz}}{\left(\displaystyle \frac{\gamma }{{10}^{2}}\right)}^{3}{\left(\displaystyle \frac{{R}_{{\rm{cur}}}}{{10}^{7}{\rm{cm}}}\right)}^{-1}.\end{eqnarray} \tag{ 28 }$$

The obtained values are roughly consistent with the observations.

The important point is that the extraction of ions from the NSS occurs for both signs of ${\rho }_{{\rm{GJ}}}$, if the backflowing particles exist. Most studies including the ion extraction from the NSS focus on the anti-pulsar (${\rho }_{{\rm{GJ}}}\gt 0$ at the polar cap; e.g., Chen & Ruderman 1980; Jones 2010, 2016). In this paper, we consider the screened state with ${\alpha }_{{\rm{ns}}}\ne 0$. The current density of the beam component is either ${\beta }_{{\rm{bk,bm}}}^{{\rm{req}}}\lt -1$ or ${\beta }_{{\rm{bk,bm}}}\gt 0$. For ${\rho }_{{\rm{GJ}}}\lt 0$ and ${\alpha }_{{\rm{m}}}\lt -1$, the ions are extracted from the NSS. In this case, the quasi-periodic behavior in Section 4.3.2 is expected to occur. For ${\rho }_{{\rm{GJ}}}\gt 0$, because of the symmetry to the point $({\alpha }_{{\rm{bk}}},{\alpha }_{{\rm{ns}}})=(-1,0)$ in the ${\alpha }_{{\rm{bk}}}\mbox{--}{\alpha }_{{\rm{ns}}}$ diagram (Figure 3), ions are extracted from the NSS for the case ${\alpha }_{{\rm{m}}}\gt -1$ as discussed in Section 4.1. Therefore, the quasi-periodic behavior in Section 4.3.2 is possible for both signs of ${\rho }_{{\rm{GJ}}}$.

In order to activate the OAR, the electric field along the magnetic field just above the NSS should be screened out. In this paper, we investigate the condition on the electric field screening just above the NSS with backflowing particles. We have focused on the case of the so-called anti-GJ condition, which is established along the magnetic fields that connect to the OAR. Without the backflows, the electric field cannot be screened out by the particles extracted from the NSS alone.

First, we consider the case in which particles accelerated in the OAR with a Lorentz factor $\gamma \gg 1$ (the beam component) are flowing back to the NSS (Figure 2(a)). The current and charge densities of the beam component and particles from the NSS contribute to the total ones. Then, even in the anti-GJ case, there is some screened solutions in a certain region for the combination of the current densities ${\alpha }_{{\rm{m}}}$ and ${\alpha }_{{\rm{bk,bm}}}$ (Figure 3). We analytically introduce a parameter ${\beta }_{{\rm{ns}}}^{{\rm{req}}}$ to characterize the required flow speed of the particles extracted from the NSS. The electric field is not screened in the conditions that are so-called R-super-GJ and R-anti-GJ, which are defined by ${\beta }_{{\rm{ns}}}^{{\rm{req}}}$.

However, we can expect a quasithermal plasma with a mildly relativistic temperature to screen out the electric field between the OAR and the outer boundary of the polar cap region (Figure 2(b)). Some fraction of such thermal components can also flow back to the NSS. Using numerical simulations, we show that the thermal component can adjust both the current and charge densities to the required values. We obtain the minimum number density of the thermal component to screen the electric field in Equations (22) and (24).

We also investigate the ion-extracted case from the NSS. The difference of masses between electrons and ions causes bunches of particles, which are formed quasi-periodically. The period is linearly proportional to the length of the calculation domain. This may be an important process for coherent radio emission.

We are grateful to the anonymous referee for helpful comments. We would like to thank Kunihito Ioka, Tsunehiko N. Kato, Yutaka Ohira, Shinpei Shibata, and Shuta J. Tanaka for fruitful discussions. This work is supported by KAKENHI 24103006, 16J06773 (S.K.), 24540258, and 15K05069 (K.A., T.T.).