On the Maximum Pair Multiplicity of Pulsar Cascades

Dense electron–positron pair plasmas are an integral part of the standard model for rotation-powered pulsars, which was initially proposed by Goldreich & Julian (1969) and Sturrock (1971). According to the standard model, a pulsar magnetosphere is filled with dense pair plasma, which screens the electric field along magnetic field lines everywhere, except in some small zones responsible for particle acceleration and emission. The sharpness of peaks in pulsar light curves is a strong argument in favor of thin acceleration zones and screened electric field in most parts of the pulsar magnetosphere. There is also direct observational evidence for plasma creation in pulsars: the most energetic ones are surrounded by "cocoons" of dense relativistic plasma—pulsar wind nebulae (PWNe)—powered by plasma outflow from their pulsars. Understanding pair creation is important for unraveling the mystery of the pulsar emission mechanism(s) and understanding pulsar surroundings on both small and large scales. Pair plasma flows out of the magnetosphere, providing the radiating particles for PWNe, and could make a make significant contribution to the lepton component of cosmic rays.

The regions responsible for production of most of the pair plasma are believed to be pulsar polar caps (PCs)—small regions near the magnetic poles (Sturrock 1971; Ruderman & Sutherland 1975; Arons & Scharlemann 1979). Without dense plasma produced in the PCs, at the base of open magnetic field lines, the magnetosphere would have large volumes with unscreened electric field, as pair creation in, e.g., outer gaps (Cheng et al. 1976) cannot screen the electric field over the rest of the magnetosphere. The physical process responsible for pair production in the PCs is the conversion of high-energy γ-rays into electron–positron pairs in strong (≳10¹¹ G) magnetic fields. According to the recent self-consistent PC models, specific regions of the PC intermittently become charge starved, when the number density of charged particles is not enough to support both the change and the current density required by the global structure of the pulsar magnetosphere (Timokhin 2010; Timokhin & Arons 2013). This gives rise to a strong accelerating electric field and formation of (intermittent) accelerating zone(s). Some charged particles enter these zones, are accelerated to very high energies, and emit γ-rays, creating electron–positron pairs. The pairs can also emit pair-producing photons and so the avalanche develops until photons emitted by the last generation of pairs can no longer produce pairs and escape the magnetosphere.

The cascade process in pulsar PCs has been the subject of extensive studies (e.g., Daugherty & Harding 1982; Gurevich & Istomin 1985; Zhang & Harding 2000; Hibschman & Arons 2001; Medin & Lai 2010). Those works considered pair creation together with particle acceleration, and so the results were dependent on the acceleration model used. The most popular acceleration model assumed steady, time-independent acceleration of the primary particles in a flow with a relatively weak accelerating electric field (Arons & Scharlemann 1979; Muslimov & Tsygan 1992), which recently was shown to be incorrect by means of direct self-consistent numerical simulations (Timokhin & Arons 2013). The necessity of bringing PC pair-creation models up to date with the self-consistent description of pair acceleration motivated us to develop a simple semianalytical model for pair cascades in pulsar PCs that can be easily decoupled from the details of the particle acceleration model and that allows easy exploration of the parameter space (Timokhin & Harding 2015, hereafter Paper I). In Paper I, we considered cascades at PCs of young pulsars with moderate magnetic fields, when the dominant processes of high-energy photon emission are the curvature radiation (CR) of primary particles and the synchrotron radiation (SR) of secondary particles. This model agrees very well with the results of elaborate numerical simulations of pair cascade. We have shown that the maximum pair multiplicity achievable in pulsars does not exceed $\sim \mathrm{few}\times {10}^{5}$, which sets stringent limits on PWN models.

In this paper, we present a significant improvement to the semianalytical model from Paper I. The new model allows inclusion of additional emission processes, provides detailed information about the spatial distribution of pair creation, and is applicable for pulsars with strong, up to ∼10¹³ G, magnetic fields. More specifically, the new model (i) can include an arbitrary number of emission processes and collect detailed information about where and in what cascade branch pairs are created, (ii) incorporates strong field corrections to the expression for the attenuation coefficient for one-photon pair creation, and (iii) takes into account the effect of photons splitting on the cascade multiplicity.

The main question we try to answer in this paper is: what is the maximum number density of pair plasma a pulsar can generate? It has been shown in numerous studies that the highest cascade multiplicity is expected in young energetic pulsars with high accelerating electric fields, where cascades are initiated by CR of primary particles. In such pulsars, primary particles are accelerated to higher energies over a short distance; emit photons via CR, which is the most efficient radiation mechanism in such physical conditions; and those photons propagate a short distance before being absorbed in a still strong magnetic field. In Paper I, we limited ourselves to CR–synchrotron cascades, without considering resonant inverse Compton scattering (RICS) of thermal photons from the neutron star (NS) surface by particles in the cascade. As we argued in that paper, it was an adequate approximation for most young energetic pulsars. However, for B ≳ 10¹² G, right where CR–synchrotron cascades reach their highest multiplicity, RICS becomes an important emission mechanism, while inverse Compton scattering in the non-resonant regime remains irrelevant for PC cascades. In order to get an accurate limit on the maximum cascade multiplicity, RICS must be taken into account.

In this paper, we apply our new semianalytical model to cascades where pairs are created by photons emitted via CR of primary particles and by photons emitted by secondary particles via SR as well as RICS of soft X-rays from the NS surface.³ Similar to Paper I, we consider the physical processes in pair cascades and particle acceleration models separately to clearly set apart different factors influencing the efficiency of pair cascades. Results presented in this paper supersede results of Paper I for high-field (around ∼10¹³ G) pulsars and improve multiplicity estimates for pulsars with medium magnetic fields (∼10¹² G) covering all ranges of parameters for pulsars capable of generating high-multiplicity pair plasma.

We do not attempt to model how the entire magnetosphere is filled with plasma (like e.g., Philippov et al. 2015; Brambilla et al. 2018). We adopt the standard pulsar model and concentrate on the microphysics of the PC cascade zone—the most important supplier of pair plasma in the magnetosphere—to determine the upper limit on pair plasma density that a pulsar can generate.

The plan of the paper is as follows. In Section 2, we first discuss general properties of electron–positron cascades and then give an overview of the microphysical processes in PC cascades. In Section 3, we consider photon absorption in strong magnetic fields: single photon pair creation in Section 3.1, photon splitting in Section 3.2, and the energy of photons escaping from the cascade in Section 3.3. We discuss particle acceleration in Section 4. In Section 5, we give a simple estimate for the upper limit of the cascade multiplicity from first principles. Section 6 gives an overview of our semianalytical cascade model (with more technical details described in Appendix C). The main results are described in Section 7. We summarize our findings and discuss the limitations of our model in Section 8.

In pulsar magnetospheres, electron–positron pairs can be created by single-photon absorption in a strong magnetic field ($\gamma B$), which can happen only close to the NS where the magnetic field is strong enough, and in two-photon collisions (γγ), which are relevant mostly in the outer magnetosphere. In an electron–positron cascade, primary particles lose energy by some emission mechanism, creating high-energy photons that are absorbed in a pair-creation process and produce electron–positron pairs. Pairs can also emit high-energy photons, which then create the next generation of pairs. As the cascade develops, it "alternates" between electron/positron and photon states. At each step in the cascade, the energy of the parent particle is divided between secondary particles. Each subsequent generation of particles has smaller energies than the previous one. At some cascade generation, the energy of the photons drops below the pair formation threshold and the cascade terminates. If the energy of the parent particle is divided roughly equally between its secondaries, i.e., the photon's energy is roughly equally divided between electron and positron, and each pair member emits several hard photons of approximately the same energy, then at the last cascade step, the available energy will be approximately equally split between photons with energies somewhat above the pair formation threshold. These photons will create the last generation of pairs. The number of pairs in such a cascade grows as a geometric progression at each generation, and most of the pairs will be created at the last cascade step. In an ideal case, when both primary and secondary particles radiate all their energies as pair-producing photons, the multiplicity κ of such a cascade (the number of particles produced by each primary particle) would be

$$\begin{eqnarray}&&{\kappa }_{\max }\simeq 2\,\displaystyle \frac{{\epsilon }_{p}}{{\epsilon }_{\gamma ,\mathrm{esc}}},\end{eqnarray} \tag{ 1 }$$

where ${\epsilon }_{\gamma ,\mathrm{esc}}$ is the maximum energy of the photons that escape the cascade (or the minimum energy of pair-producing photons), and ${\epsilon }_{p}$ is the energy of primary particles. For convenience from here on, all particle and photon energies will be quoted in terms of m_ec². In a real cascade, both primary and secondary particles do not radiate all their kinetic energy as pair-producing photons, and ${\kappa }_{\max }$ can be considered as an upper limit on the multiplicity.

In the above ideal limit, the energy of the primary particle is divided into chunks of the size ${\epsilon }_{\gamma ,\mathrm{esc}}$. Usually, ${\epsilon }_{\gamma ,\mathrm{esc}}\gg 2$ and even in the ideal case, when energy is not lost at intermediate steps, the cascade multiplicity is much smaller than ${\epsilon }_{p}$ (in terms of m_ec²), which would be the case if all the energy of the primary particles went into the rest energy of pairs.⁴

In the de facto standard pulsar model, particles can be accelerated to very high energies and produce dense pair plasma in the PCs (Sturrock 1971), in thin regions along last closed magnetic field lines (the "slot gap" model of Arons 1983), and in the "outer gaps," regions in the outer magnetosphere along magnetic field lines crossing the surface where the Goldreich & Julian (1969) charge density changes sign (Cheng et al. 1976). The outer and slot gaps occupy only a relatively small volume of the magnetosphere so that most of the open magnetic field lines do not pass though them. All open field lines originate in the PCs, and a significant fraction of them pass through PC particle acceleration zones. The total number of primary particles in the PC cascades is much larger than that in cascades in the outer pulsar magnetosphere. Moreover, simulation of the outer gap cascades predicts multiplicities not higher than 10⁴ (e.g., Hirotani 2006). So, at least in the standard pulsar model, most of the pairs are produced in the PC cascades.

It was demonstrated in Timokhin & Arons (2013) and Timokhin (2010) that pair formation in pulsars is an intermittent process. Time periods of efficient particle acceleration and intense pair production alternate with periods of quiet plasma flow when dense plasma screens the accelerating electric field and no pairs are formed (more on this in Section 4). As in Paper I, here we consider cascades at the peak of the pair formation cycle, when their multiplicity is the highest, postponing discussion on the effects of intermittency to Section 8. Such cascades are generated by the primary particles accelerated in well-developed gaps. Screening of the accelerating electric field in the gap happens very quickly, well before the multiplicity reaches its maximum values. Once primary particles have produced the first generation of pairs that screen the accelerating electric field, they keep moving in the regions of screened electric field, radiating their energy away and giving rise to extensive pair cascades. So, the PC cascades can be considered as initiated by primary particles with given energies freely moving along magnetic field lines.

Figure 1 gives a schematic overview of processes involved in pair plasma generation in PC cascades;⁵ shown are the first two generations in a cascade initiated by a primary electron. Primary electrons emit CR photons (${\gamma }_{{\rm{cr}}}$) almost tangent to the magnetic field lines; primary electrons and CR photons are generation 0 particles in our notation.⁶ Magnetic field lines are curved, and the angle between the photon momentum and the magnetic field grows as the photon propagates farther from the emission point. When this angle becomes large enough, photons are absorbed and each photon creates an electron–positron pair—a generation 1 electron (${e}^{-}$) and positron (${e}^{+}$). The pair momentum is directed along the momentum of the parent photon, and at the moment of creation, the particles have non-zero momentum perpendicular to the magnetic field. They radiate this perpendicular momentum almost instantaneously via SR and move along magnetic field lines. The secondary particles can scatter thermal X-ray photons ${\gamma }_{X}$ coming from the NS surface and lose their momenta parallel to the magnetic field. Inverse Compton scattering of the thermal photons in the non-resonant regime is very ineffective in converting the energy of the parallel motion of pairs into pair-producing photons and can be ignored; see Appendix A. On the other hand, RICS, as it was first pointed out by Dermer (1990), can become a very efficient emission process in high-field pulsars. Although the secondary particles are relativistic, their energy is much lower than that of the primaries, and their curvature photons cannot create pairs. Generation 1 photons—synchrotron (${\gamma }_{\mathrm{syn}}$) and RICS (${\gamma }_{{\rm{rics}}}$) photons produced by generation 1 particles—are also emitted (almost) tangent to the magnetic field line, as the secondary particles are relativistic, and propagate some distance before acquiring the necessary angle to the magnetic field and creating generation 2 pairs. These pairs in their turn radiate their perpendicular momentum via SR and parallel momenta via RICS, emitting generation 2 photons. The cascade initiated by a single CR photon stops at a generation where the energy of synchrotron photons falls below ${\epsilon }_{\gamma ,\mathrm{esc}}$.

Zoom In Zoom Out Reset image size

Primary particles emit pair-producing CR photons throughout the entire cascade zone as they move along the field lines. Secondary particles emit all their pair-producing synchrotron photons after their creation almost instantaneously. RICS photons are emitted by secondary particles over some distance, which is usually much smaller than the size of the cascade zone.

In the following sections, we analyze the individual factors regulating the yield of electron–positron cascades and develop a semianalytical technique which models cascade development by following the general picture outline above.

3.1. Pair Creation

For the opacity for single photon pair production in a strong magnetic field $\gamma \to {e}^{+}{e}^{-}$, we use the prescription suggested by Daugherty & Harding (1983) which can be written as

$$\begin{eqnarray}&&{\alpha }_{\gamma \to \pm }({\epsilon }_{\gamma },\psi )=0.23\,\displaystyle \frac{{\alpha }_{f}}{{\unicode{x0019B}}_{{\rm{c}}}}\ b\,\sin \psi \,\exp \left(-\displaystyle \frac{4}{3\chi }\right){f}_{\alpha ,1},\end{eqnarray} \tag{ 2 }$$

where $b\equiv B/{B}_{q}$ is the local magnetic field strength B normalized to the critical quantum magnetic field ${B}_{q}=e/{\alpha }_{f}{\unicode{x0019B}}_{{\rm{c}}}^{2}\,=4.41\times {10}^{13}$ G, ψ is the angle between the photon momentum and the local magnetic field, ${\alpha }_{f}={e}^{2}/{\hslash }c\approx 1/137$ is the fine structure constant, and ${\unicode{x0019B}}_{{\rm{c}}}={\hslash }/{mc}=3.86\times {10}^{-11}$ cm is the reduced Compton wavelength. The parameter χ is defined as

$$\begin{eqnarray}&&\chi \equiv \displaystyle \frac{1}{2}\,{\epsilon }_{\gamma }b\,\sin \psi ,\end{eqnarray} \tag{ 3 }$$

where ${\epsilon }_{\gamma }$ is the photon energy in units of m_ec². Expression (2) differs from the usual Erber (1966) formula through the term

$$\begin{eqnarray}{f}_{\alpha ,1}=\left\{\begin{array}{ll}\exp \left(-0.56\,\displaystyle \frac{{b}^{2.6962}}{{\chi }^{3.7}}\right), & \mathrm{if}\ {\epsilon }_{\gamma }\sin \psi \geqslant 2\\ 0, & \mathrm{if}\ {\epsilon }_{\gamma }\sin \psi \lt 2\end{array}\right..\end{eqnarray} \tag{ 4 }$$

The function ${f}_{\alpha ,1}$ ensures that the attenuation coefficient for pair production becomes zero below the threshold ${\epsilon }_{\gamma }\sin \psi =2$ and corrects ${\alpha }_{\gamma \to \pm }$ for the case when absorption happens near the threshold. The threshold condition for $\gamma B$ pair production can be expressed in terms of χ as

$$\begin{eqnarray}&&\chi \geqslant b.\end{eqnarray} \tag{ 5 }$$

Expression (2) works for high, $B\gt 3\times {10}^{12}$ G, magnetic fields; for weaker fields, it reduces to the well-known Erber's formula (see Appendix B).

The optical depths for pair creation by a high-energy photon in a strong magnetic field after propagating a distance l is

$$\begin{eqnarray}&&\tau ({\epsilon }_{\gamma },l)={\int }_{0}^{l}{\alpha }_{\gamma \to \pm }({\epsilon }_{\gamma },\psi (x)){dx},\end{eqnarray} \tag{ 6 }$$

where integration is along the photon's trajectory. For photons emitted tangent to the magnetic field line, ${dx}={\rho }_{{\rm{c}}}d\psi$, where ${\rho }_{{\rm{c}}}$ is the radius of curvature of the magnetic field lines. From Equation (3), we have $\psi =2\chi /{\epsilon }_{\gamma }b$, and substituting it into Equation (6), we can express the optical depth τ to pair production as an integral over χ as

$$\begin{eqnarray}&&\tau ({\epsilon }_{\gamma },\tilde{\chi })=\displaystyle \frac{{A}_{\tau }}{{\epsilon }_{\gamma }^{2}}\,{\int }_{0}^{\tilde{\chi }}\displaystyle \frac{{\rho }_{{\rm{c}}}}{b}\chi \exp \left(-\displaystyle \frac{4}{3\chi }\right){f}_{\alpha ,1}\,d\chi ,\end{eqnarray} \tag{ 7 }$$

where ${A}_{\tau }\equiv 0.92{\alpha }_{f}/{\unicode{x0019B}}_{{\rm{c}}}\approx 1.74\times {10}^{8}\,{{\rm{cm}}}^{-1}$. The optical depth depends exponentially on χ, and the main contribution to the integral comes from the values of χ close to the upper boundary $\tilde{\chi }$. For a wide range of photon energies and field strengths, the value of χ at the point where the photon is absorbed, ${\chi }_{a}$, changes slowly. The mean free path (mfp) of photons can be estimated from Equation (3) as

$$\begin{eqnarray}&&{\lambda }_{\gamma \to \pm }=2{\rho }_{{\rm{c}}}{\chi }_{a}\displaystyle \frac{1}{b{\epsilon }_{\gamma }}.\end{eqnarray} \tag{ 8 }$$

Both ${\chi }_{a}$ and ${\rho }_{{\rm{c}}}$ change slower than b and ${\epsilon }_{\gamma }$ as the cascade develops. In each cascade generation, the energy of particles and photons is smaller than that in the preceding generation. The photon mfp ${\lambda }_{\gamma \to \pm }$ increases because of this. If ${\lambda }_{\gamma \to \pm }$ becomes comparable to the characteristic scale of the magnetic field variation L_B, then the increase of ${\lambda }_{\gamma \to \pm }$ for the next-generation photons will be compounded by an additional decrease of the magnetic field b as well, by at least ∼an order of magnitude (for a dipolar field). In most cases, the magnetic field at the anticipated absorption point for the next generation of photons will drop below the pair formation threshold (5). Hence, the cascade generation for which ${\lambda }_{\gamma \to \pm }\sim {L}_{B}$ should be the final one.

We consider strong cascades with large multiplicities; such cascades fully develop before ${\lambda }_{\gamma \to \pm }\sim {L}_{B}$. For such cascades in the region where most of the pairs are produced, the magnetic field b and the radius of curvature of the magnetic field lines ${\rho }_{{\rm{c}}}$ are approximately constant. In approximation of the constant b and ${\rho }_{{\rm{c}}}$, Equation (7) can be written as

$$\begin{eqnarray}&&\tau (\chi )={A}_{\tau }\displaystyle \frac{{\rho }_{{\rm{c}}}}{{\epsilon }_{\gamma }^{2}\,b}\,{\int }_{0}^{{\chi }_{a}}\chi \exp \left(-\displaystyle \frac{4}{3\chi }\right){f}_{\alpha ,1}\,d\chi \,\end{eqnarray} \tag{ 9 }$$

and integrated analytically. The resulting expression is quite cumbersome; it is derived in Appendix B and given by Equation (50).

Because the opacity to pair production depends exponentially on χ, it is a reasonable approximation that all photons are absorbed when they have traveled the distance where $\tau =1$. We define ${\chi }_{a}$ as the value of χ, where the optical depth reaches 1 through

$$\begin{eqnarray}&&{\chi }_{a}:\tau ({\chi }_{a})=1.\end{eqnarray} \tag{ 10 }$$

${\chi }_{a}$ is a solution of the nonlinear Equation (10), with τ given by Equation (50). We solved Equation (10) numerically for different values of ${\epsilon }_{\gamma }$, b, and ${\rho }_{{\rm{c}}}$. In Figure 2, we plot contours of $1/{\chi }_{a}$ as functions of $\mathrm{log}({\epsilon }_{\gamma })$ and $\mathrm{log}(B)$ for three different values of the radius of curvature of the magnetic field lines, ${\rho }_{{\rm{c}}}={10}^{6},{10}^{7},{\rm{and}}\,{10}^{8}\,{\rm{cm}}$. The smallest value of ${\rho }_{{\rm{c}}}$ corresponds to a strongly multipolar PC magnetic field, when the radius of curvature is comparable to the NS radius. The largest value corresponds to the radius of curvature of dipolar magnetic field lines, which at the NS surface is given by

$$\begin{eqnarray}&&{\rho }_{{\rm{c}},\mathrm{dip}}=\displaystyle \frac{4}{3}\displaystyle \frac{{R}_{{\rm{ns}}}}{\theta }=9.2\times {10}^{7}{\left(\displaystyle \frac{\theta }{{\theta }_{\mathrm{pc}}}\right)}^{-1}{P}^{1/2}\,{\rm{cm}},\end{eqnarray} \tag{ 11 }$$

where θ is the colatitude of the footpoint of the magnetic field line, ${\theta }_{\mathrm{pc}}=1.45\times {10}^{-2}{P}^{-1/2}$—the colatitude of the PC boundary, P—the pulsar period in seconds, and ${R}_{{\rm{ns}}}={10}^{6}$ cm—the NS radius.

Zoom In Zoom Out Reset image size

$1/{\chi }_{a}$ is a smooth function of $\mathrm{log}({\epsilon }_{\gamma })$, $\mathrm{log}(B)$, and $\mathrm{log}{\rho }_{{\rm{c}}}$—as it is to be expected from ${\alpha }_{\gamma \to \pm }\propto \exp (-1/\chi )$—and can be accurately approximated using a modest-size numerical table. The change in the behavior of $1/{\chi }_{a}$ for large magnetic fields, when the contour lines become horizontal, is due to photon being absorbed close to the pair production threshold Equation (5).

The values of ${\chi }_{a}$ we obtained here using a more accurate expression for the $\gamma B$ opacity, Equation (50), are up to 40% higher than that from Paper I where we used Erber's formula and made a simple correction for the pair formation threshold by setting the upper limit on χ according to Equation (5). This difference is larger than 10% only for magnetic fields with $B\,\gtrsim 1.5\times {10}^{12}{\rm{G}}$ (compare Figure 2 with Figure 3 from Paper I). Also note that the values of ${\chi }_{a}$ differ significantly from the often-used value ${\chi }_{a}=1/15$ first suggested by Ruderman & Sutherland (1975), especially for higher energy photons.

The contour plots of the mfp of the photons ${\lambda }_{\gamma \to \pm }$ emitted tangentially to the magnetic field lines are shown in Figure 3; ${\lambda }_{\gamma \to \pm }$ was calculated according to Equation (8). As expected, it scales linearly with $1/{\epsilon }_{\gamma }$ and $1/B$; the deviation from the linear behavior is seen only for the combination of ${\epsilon }_{\gamma }$ and B when pair formation happens near the threshold, at $B\gtrsim 1.5\times {10}^{12}$ G.

Zoom In Zoom Out Reset image size

3.2. Photon Splitting

Single photon pair creation is not the only process responsible for photon attenuation in a strong magnetic field, albeit it is the most significant one. The most important competing process to pair creation is magnetic photon splitting $\gamma \to \gamma \gamma$. Aside from the end product, the major differences between pair creation and photon splitting are (i) pair creation is a first-order QED process and splitting is a third-order one; therefore, splitting is weaker than pair creation on the order of ${\alpha }_{f}^{2};$ (ii) in contrast to pair creation, splitting has no threshold for photon's energy, (iii) while in moderately strong magnetic fields $B\lesssim {B}_{\mathrm{cr}}$ for pair creation, both modes of photon polarization (∥—when the photon's electric field is parallel to the plane containing ${\boldsymbol{B}}$ and photon's momentum, ⊥—the photon's electric field is perpendicular to that plane) have similar cross-sections and threshold conditions, below the pair formation threshold, photon splitting is allowed only for the process $\perp \to \parallel \parallel$ (Adler 1971; Usov 2002).

Radiation processes relevant for secondary particles in PC cascades of energetic pulsars (SR, RICS) produce predominantly ⊥ polarized photons. Despite the inherently smaller cross-section of magnetic splitting, the absence of an energy threshold could allow photons to split before acquiring large-enough angles to the magnetic field to produce pairs, thus reducing cascade multiplicity. As we are interested in the most efficient cascades, a regime where magnetic splitting becomes important is beyond the scope of this paper. More details on cascade kinetics in the presence of photon splitting can be found in Harding et al. (1997) and Baring & Harding (1997, 2001). Here we want to establish the boundary in the parameter space where photon splitting starts affecting cascade multiplicity. To do so, we consider a case of photon splitting $\perp \to \parallel \parallel$ for photons below the pair formation threshold.

The attenuation coefficient for photon splitting $\perp \to \parallel \parallel$ is (Baring & Harding 2001; Baring 2008)

$$\begin{eqnarray}&&{\alpha }_{\gamma \to \gamma \gamma }({\epsilon }_{\gamma },\psi )=\displaystyle \frac{{\alpha }_{f}^{3}c}{60{\pi }^{2}{\unicode{x0019B}}_{{\rm{c}}}}{\epsilon }_{\gamma }^{5}{b}^{6}\sin {\psi }^{6}{{ \mathcal M }}_{1}^{2}\end{eqnarray} \tag{ 12 }$$

at low values of the magnetic field perpendicular to the photon's trajectory—for photons below the pair formation threshold—the scattering amplitude ${{ \mathcal M }}_{1}$ is a constant independent of b, ${{ \mathcal M }}_{1}\approx 26/315$. Integration of ${\alpha }_{\gamma \to \gamma \gamma }$ over the distance gives the optical depth for photon splitting (see Equation (6)); the mfp for splitting can then be estimated as

$$\begin{eqnarray}&&{\lambda }_{\gamma \to \gamma \gamma }=1.8\,{\epsilon }_{\gamma }^{-5/7}{b}^{-6/7}{\rho }_{{\rm{c}}}^{6/7}\,{\rm{cm}}.\end{eqnarray} \tag{ 13 }$$

When the mfp for splitting becomes smaller than the mfp for pair formation, ${\lambda }_{\gamma \to \gamma \gamma }\lt {\lambda }_{\gamma \to \pm }$, the photon splits before producing a pair. In Figure 4, we plot the ratio ${\lambda }_{\gamma \to \gamma \gamma }$/${\lambda }_{\gamma \to \pm }$ as a function of magnetic field strength and photon energy, for three values of the radius of curvature of magnetic field lines. It is evident from these plots that splitting is an important attenuation mechanism only for strong magnetic fields and low-energy photons. With the increase of the magnetic field, photon splitting starts affecting first the last cascade generation where the energy of the photons becomes low. If these photons split, the resulting photons will be below the pair formation threshold. In the PC cascades, most of the pairs are produced in the last cascade generation, and when photon splitting becomes important, the cascade multiplicity can drop significantly. The exact fraction of perpendicularly polarized photons in the cascade—which are subject to splitting—in general depends on particle energy distributions, but it is more than 50% (e.g., Baring & Harding 2001). Hence, the multiplicity of the pair cascade will drop by at least a factor of 2 when photon splitting becomes important.

Zoom In Zoom Out Reset image size

The critical magnetic field strength ${B}_{\gamma \to \gamma \gamma }$ above which cascade multiplicity becomes affected by photon splitting is the field strength when ${\lambda }_{\gamma \to \gamma \gamma }\lt {\lambda }_{\gamma \to \pm }$ for the last-generation photons, i.e., the photons with the escaping energy ${\epsilon }_{\gamma ,\mathrm{esc}}$ (which we calculate in the next section),

$$\begin{eqnarray}&&{B}_{\gamma \to \gamma \gamma }:{\lambda }_{\gamma \to \gamma \gamma }({\epsilon }_{\gamma ,\mathrm{esc}})={\lambda }_{\gamma \to \pm }({\epsilon }_{\gamma ,\mathrm{esc}}).\end{eqnarray} \tag{ 14 }$$

3.3. Energy of Escaping Photons

As we discussed above, photons escaping the cascade are those with mfp larger than the characteristic scale of magnetic field variation ${\lambda }_{\gamma \to \pm }\gt {L}_{B}$. The formal criterion we use to calculate the energy of escaping photons ${\epsilon }_{\gamma ,\mathrm{esc}}$ is ${\lambda }_{\gamma \to \pm }({\epsilon }_{\gamma ,\mathrm{esc}})={s}_{\mathrm{esc}}{R}_{{\rm{ns}}};$ ${s}_{\mathrm{esc}}$ is a dimensionless parameter quantifying the escaping distance in units of ${R}_{{\rm{ns}}}$. From the expression for mfp, Equation (8), we get a (nonlinear)⁷ equation for ${\epsilon }_{\gamma ,\mathrm{esc}}$,

$$\begin{eqnarray}&&{\epsilon }_{\gamma ,\mathrm{esc}}=2\,\displaystyle \frac{{\rho }_{{\rm{c}}}}{{s}_{\mathrm{esc}}{R}_{{\rm{ns}}}}\displaystyle \frac{{\chi }_{a}}{b}.\end{eqnarray} \tag{ 15 }$$

Any global NS magnetic field near the surface decays with distance as ${(r/{R}_{{\rm{ns}}})}^{-\delta }$, $\delta \geqslant 3;$ a dipole field, δ = 3, is often considered as a reasonable assumption. A pure dipole, however, seems to be too idealized an approximation, as the NS magnetic field is slightly disturbed by the currents flowing in the magnetosphere. PC cascade models should consider at least near-dipole magnetic fields with different curvatures of magnetic field lines. Hence, a reasonable estimate for L_B would be the distance of the order of the NS radius ${R}_{{\rm{ns}}}$. For our approximation of constant B and ${\rho }_{{\rm{c}}}$, we found that the value ${L}_{B}={R}_{{\rm{ns}}}/2$—at that distance from the NS the magnetic field decays by at least a factor of 3—provides a good fit to the results of the numerical simulations described in Paper I.⁸ Results described in this paper are obtained assuming ${s}_{\mathrm{esc}}=0.5$.

In Figure 5, we plot the energy of escaping photons, $\mathrm{log}{\epsilon }_{\gamma ,\mathrm{esc}}$, as a function of the radius of curvature of magnetic field lines ${\rho }_{{\rm{c}}}$ and magnetic field strength B for ${s}_{\mathrm{esc}}=1$. This figure shows a (obvious) trend wherein for higher magnetic fields and smaller radii of curvature, the energy of escaping photons is lower. The deviation of contours from straight lines for $B\gtrsim 3\times {10}^{12}{\rm{G}}$ is due to the change in ${\chi }_{a}$ behavior near the pair formation threshold (see Figure 2 and the subsequent paragraph). For different values of ${s}_{\mathrm{esc}}$, the escape energy ${\epsilon }_{\gamma ,\mathrm{esc}}$ could be estimated from Figure 3.

Zoom In Zoom Out Reset image size

The critical magnetic field ${B}_{\gamma \to \gamma \gamma }$ above which photon splitting starts affecting cascade multiplicity is shown in Figure 5 by the dashed line; ${B}_{\gamma \to \gamma \gamma }$ is calculated from Equation (14). Cascade multiplicity will drop due to photon splitting at $B\gtrsim 4.4\times {10}^{12}{\rm{G}}$ for ${\rho }_{{\rm{c}}}={10}^{6}$cm and at $B\gtrsim 1.3\times {10}^{13}{\rm{G}}$ for ${\rho }_{{\rm{c}}}={10}^{8}$ cm. The increase of ${B}_{\gamma \to \gamma \gamma }$ for larger ${\rho }_{{\rm{c}}}$ is due to the increase of the energy of escaping photons—${\lambda }_{\gamma \to \pm }$ has a stronger dependence on ${\epsilon }_{\gamma }$ than ${\lambda }_{\gamma \to \gamma \gamma }$, which leads to the increase of the value of ${B}_{\gamma \to \gamma \gamma }$ according to Equations (14), (13), and (8).

It would be useful to have an analytical expression for the energy of the escaping photons; to obtain it, we construct an approximation for ${\chi }_{a,\mathrm{esc}}$. We solved Equation (15) to find ${\epsilon }_{\gamma ,\mathrm{esc}};$ now, using the interpolation formula for $1/{\chi }_{a}$, we can find $1/{\chi }_{a,\mathrm{esc}}\equiv 1/{\chi }_{a}({\epsilon }_{\gamma ,\mathrm{esc}})$. In Figure 6, we show contours of $1/{\chi }_{a,\mathrm{esc}}$ as a function of B and ${\rho }_{{\rm{c}}}$. There are two distinct regions on this plot: for $B\lesssim 3\times {10}^{12}{\rm{G}}$, ${\chi }_{a,\mathrm{esc}}$ changes very slowly, while for larger values of B, it changes significantly but does not depend on ${\rho }_{{\rm{c}}}$. For $B\gtrsim 3\times {10}^{12}{\rm{G}}$, the absorption of the last-generation photons happens near the pair formation threshold, when ${\epsilon }_{\perp }={\epsilon }_{\gamma ,\mathrm{esc}}({R}_{{\rm{ns}}}/{\rho }_{{\rm{c}}})\simeq 2$, and so ${\chi }_{a,\mathrm{esc}}\simeq b$. For weaker magnetic fields, the opacity for near-threshold photons is too small for them to be absorbed after traveling a distance ${R}_{{\rm{ns}}}$, so the last-generation photons have energies larger than the pair formation threshold and are absorbed at very similar values of ${\chi }_{a}$. We find that the following approximation works quite well:

$$\begin{eqnarray}{\chi }_{a,\mathrm{esc}}=\left\{\begin{array}{ll}b, & \mathrm{if}\ b\gt 1/15\\ 1/15, & \mathrm{if}\ b\leqslant 1/15\end{array}\right..\end{eqnarray} \tag{ 16 }$$

The energy of the escaping photons can be expressed following Equation (15) as

$$\begin{eqnarray}&&{\epsilon }_{\gamma ,\mathrm{esc}}\approx 1.8\times {10}^{3}\,\displaystyle \frac{{\rho }_{{\rm{c}},7}}{{B}_{12}}{\left(\displaystyle \frac{{s}_{\mathrm{esc}}}{0.5}\right)}^{-1}{\chi }_{a,\mathrm{esc}},\end{eqnarray} \tag{ 17 }$$

where ${\chi }_{a,\mathrm{esc}}$ is given by Equation (16), and ${B}_{12}\equiv B/{10}^{12}{\rm{G}}$ and ${\rho }_{{\rm{c}},7}\equiv {\rho }_{{\rm{c}}}/{10}^{7}$ cm. This simple prescription for ${\epsilon }_{\gamma ,\mathrm{esc}}$ deviates from the numerical values shown in Figure 5 by no more than 20% for $B\lesssim 2\times {10}^{12}{\rm{G}}$ and $B\gtrsim 8\times {10}^{12}$ G; the largest deviation is ∼60% at $B\simeq 3\times {10}^{12}$ G.

Zoom In Zoom Out Reset image size

For very small values of ${\chi }_{a,\mathrm{esc}}$, it is possible to obtain a more accurate analytical expression for ${\chi }_{a,\mathrm{esc}}$ (and ${\epsilon }_{\gamma ,\mathrm{esc}}$). Using the asymptotic expression for the optical depth in the limit of $\chi \ll 1$ derived in Appendix B, Equation (53), and substituting for the photon energy the energy of escaping photons, Equation (15), we get

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{{\chi }_{a,\mathrm{esc}}} & \approx & \displaystyle \frac{3}{4}\mathrm{ln}\left(\displaystyle \frac{3{A}_{\tau }{R}_{{\rm{ns}}}^{2}}{16}\displaystyle \frac{b}{{\rho }_{{\rm{c}}}}{\chi }_{a,\mathrm{esc}}\,{s}_{\mathrm{esc}}^{2}\right)\\ & \approx & 15.7+1.7\mathrm{log}\left(\displaystyle \frac{{B}_{12}}{{\rho }_{7}}\right)+3.5\mathrm{log}\left(\displaystyle \frac{{s}_{\mathrm{esc}}}{0.5}\right).\end{array}\end{eqnarray} \tag{ 18 }$$

Self-consistent modeling of accelerating zones in pulsar PCs (Timokhin 2010; Timokhin & Arons 2013) demonstrated that particle acceleration and pair formation are always non-stationary. Each period of intense particle acceleration and pair formation is followed by a period of quiet plasma flow when the accelerating electric field is screened, and no pairs are formed. At the end of the quiet phase, an accelerating gap begins to form—a region where plasma density is significantly smaller than the local Goldreich–Julian (GJ) number density ${\eta }_{{\rm{gj}}}/e$. Accelerating electric field in the gap increases linearly with distance as the gap grows. Charged particles entering the gap are accelerated to very high energies and emit gamma-rays, which give rise to electron–positron pair cascades. Dense pair plasma created in the cascades screens the electric field, stopping the growth of the gap. The gap does not stay at the same place but moves along magnetic field lines, roughly preserving its size (and the potential drop) for a while. Most of the pair plasma is created at or behind the trailing edge of the gap, where the high-energy particles are located.⁹ These particles move into the magnetosphere, emitting gamma-rays which convert into electron–positron pairs. Both primary and secondary particles are relativistic and so they move together, forming a blob of pair plasma whose density increases as pair formation continues. Some low-energy particles, however, "leak" from the blob, creating a tail of mildly relativistic plasma which screens the electric field behind the blob. When the blob with primary particles moves away from the PC and pair formation stops, the dense pair plasma from the tail keeps the electric field screened for a while until most of it has left the PC zones and a new cycle of pair formation begins.¹⁰

Whether the pair formation along given magnetic field lines occurs and how and it is depend on the ratio ${j}_{{\rm{m}}}/{j}_{{\rm{gj}}}$ of the current density required to support the twist of magnetic field lines in the pulsar magnetosphere (e.g., Timokhin 2006; Bai & Spitkovsky 2010), ${j}_{{\rm{m}}}\equiv (c/4\pi )| {\rm{\nabla }}\times {\boldsymbol{B}}|$, to the local GJ current density, ${j}_{{\rm{gj}}}\equiv {\eta }_{{\rm{gj}}}c$, where ${\eta }_{{\rm{gj}}}=B/{Pc}$ is the GJ charge density. Regardless of the ability of the NS surface to supply charged particles, i.e., in both the space charge-limited flow model of Arons & Scharlemann (1979) and the no-particle extraction model of Ruderman & Sutherland (1975), particle acceleration happens in essentially the same way. For the Ruderman & Sutherland (1975) regime, effective particle acceleration and pair formation are possible for almost all values of ${j}_{{\rm{m}}}/{j}_{{\rm{gj}}}$. In the space charge-limited flow regime, pair formation is not possible if $0\lt {j}_{{\rm{m}}}/{j}_{{\rm{gj}}}\lt 1$, but it is possible for all other values of ${j}_{{\rm{m}}}/{j}_{{\rm{gj}}}$. A detailed description of particle acceleration in pulsar PCs is given in Timokhin (2010) for the no-particle extraction regime and in Timokhin & Arons (2013) for the space charge-limited flow regime.

Although the character of plasma flow inferred from the self-consistent simulation of Timokhin (2010) and Timokhin & Arons (2013) qualitatively differs from that assumed in both the Ruderman & Sutherland (1975) and Arons & Scharlemann (1979) type models, the physics of particle acceleration in the gap is similar to that of the accelerating gap in the Ruderman & Sutherland (1975) model. Namely, due to the significant deviation of the charge and current densities from GJ values, the electric field in the gap is comparable to the vacuum electric field ($\sim h{\rm{\Omega }}B/c$, h is the size of the gap), and particles are accelerated over a short gap by the strong electric field, which increases linearly with the distance. In Section 6.2 of Paper I, we analyzed the physics of particle acceleration and derived an analytical expression for the energy of the primary particles accelerated in non-stationary cascades. According to Equation (41) in Paper I, the final energy of particles accelerated in the gap

$$\begin{eqnarray}\begin{array}{rcl}{\epsilon }_{\pm ,\mathrm{acc}} & = & \displaystyle \frac{49}{18}{\left(\displaystyle \frac{\pi {B}_{q}}{{\unicode{x0019B}}_{{\rm{c}}}^{3}c}\right)}^{1/7}{\chi }_{a,\mathrm{acc}}^{2/7}\,{\xi }_{j}^{1/7}{P}^{-1/7}{B}^{-1/7}{\rho }_{{\rm{c}}}^{4/7}\\ & \simeq & 5\times {10}^{7}\ {\chi }_{a,\mathrm{acc}}^{2/7}\,{\xi }_{j}^{1/7}{P}^{-1/7}{B}_{12}^{-1/7}{\rho }_{{\rm{c}},7}^{4/7}.\end{array}\end{eqnarray} \tag{ 19 }$$

${\chi }_{a,\mathrm{acc}}$ is the value of the parameter χ for photons that create pairs terminating the gap. ${\xi }_{j}$ is a factor that shows how the electric field in the gap is stronger/weaker compared to the field in a static vacuum gap of an aligned rotator,

$$\begin{eqnarray}&&{\xi }_{j}\equiv \displaystyle \frac{| j-{j}_{{\rm{m}}}| }{{j}_{{\rm{gj}}}^{0}}\left(1+\displaystyle \frac{c}{v}\right)\approx 2\,\displaystyle \frac{{j}_{{\rm{m}}}}{{j}_{{\rm{gj}}}^{0}}.\end{eqnarray} \tag{ 20 }$$

where j is the current density in the gap. In most cases, $| j-{j}_{{\rm{m}}}| \simeq {j}_{{\rm{m}}};$ ${j}_{{\rm{gj}}}^{0}$ is the GJ current density in an aligned rotator,

$$\begin{eqnarray}&&{j}_{{\rm{gj}}}^{0}\equiv {\eta }_{{\rm{gj}}}^{0}c=\displaystyle \frac{B}{P}.\end{eqnarray} \tag{ 21 }$$

v is the velocity of the gap; in most cases, the gap moves with relativistic velocities, so $v\simeq c$. Taking into account these approximations, we get the second expression for ${\xi }_{j}$ in Equation (20). ${j}_{{\rm{m}}}$, and thus ${\xi }_{j}$, depends on the pulsar inclination angle and the position of the given magnetic field line inside the PC (see, e.g., Figure 1 in Timokhin & Arons 2013). In cascades along magnetic field lines where ${j}_{{\rm{m}}}$ is close to the local value of ${j}_{{\rm{gj}}}$ in an aligned rotator ${\xi }_{j}\sim 2$, for the same situation in a pulsar with an inclination angle of 60^◦, ${\xi }_{j}\sim 1$. The energy of primary particles, Equation (19), has the same dependence on ${\rho }_{{\rm{c}}}$, P, and B as the expression for the potential drop in the gap derived by Ruderman & Sutherland (1975), their Equation (23). This is to be expected as in both cases, particles are accelerated by the electric field, which grows linearly with the distance, and the size of the gap is regulated by absorption on curvature photons in the magnetic field. The difference is in the presence of the factor ${\xi }_{j}$ and a different numerical factor.

Here we slightly improve the accuracy of this expression by calculating self-consistently the value of the parameter ${\chi }_{a,\mathrm{acc}}$. Using numerical interpolation for χ (see Section 3.1), it is easy to obtain self-consistently the energy of primary particles accelerated in the gap by solving numerically the equation

$$\begin{eqnarray}&&{\epsilon }_{\pm ,\mathrm{acc}}({\chi }_{a,\mathrm{acc}}(\epsilon ))=\epsilon ,\end{eqnarray} \tag{ 22 }$$

i.e., to find the energy of particles that emit photons terminating the gap growth, taking into account the dependence of ${\chi }_{a,\mathrm{acc}}$ on particle energy (via the energy of emitted photons). In Figure 7, we show the energy of accelerated particles for a pulsar with P = 33 ms. The contours of constant ${\chi }_{a,\mathrm{acc}}$ deviate only slightly from straight lines corresponding to $\propto {\rho }_{{\rm{c}}}^{4/7}{B}^{-1/7}$ for higher values of B. This deviation is due to the variation of ${\chi }_{a,\mathrm{acc}}$ near the pair formation threshold, and so expression (19) can be safely used in many cases with a constant value of ${\chi }_{a,\mathrm{acc}}$. ${\chi }_{a,\mathrm{acc}}$ itself varies slowly with gap parameters. In Figure 8, we show $1/{\chi }_{a,\mathrm{acc}}$ for a pulsar with P = 33 ms; it was obtained together with values from Equation (22). It is evident that the value of ${\chi }_{a,\mathrm{acc}}=1/5.5$ can be used in Equation (19) for pulsars with P = 33 ms. The dependence of ${\chi }_{a,\mathrm{acc}}$ on the pulsar period P and parameter ${\xi }_{j}$ is also quite weak; see Table 1 where we show the variation of ${\chi }_{a,\mathrm{acc}}$ for $B={10}^{12}$ G, ${\rho }_{{\rm{c}}}={10}^{7}$ cm with P, and ${\xi }_{j}$. For estimates of the primary particle energies, one can use ${\chi }_{a,\mathrm{acc}}=1/7$. The dependence of the particle energy is weak, and along most of the magnetic field lines in the PC, ${\xi }_{j}$ is no more than ∼an order of magnitude lower than 2. Assuming ${\chi }_{a,\mathrm{acc}}=1/7$ and ${\xi }_{j}=2$, Equation(19) can be written as

$$\begin{eqnarray}&&{\epsilon }_{\pm ,\mathrm{acc}}\simeq 3.2\times {10}^{7}\ {P}^{-1/7}{B}_{12}^{-1/7}{\rho }_{{\rm{c}},7}^{4/7}.\end{eqnarray} \tag{ 23 }$$

This estimate for the primary particle energy can be used for a wide range of parameters of young energetic pulsars. It is ≃4 times higher than that given by Equation (23) in Ruderman & Sutherland (1975).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  $1/{\chi }_{a,\mathrm{acc}}$ for $B={10}^{12}$ G, ${\rho }_{{\rm{c}}}={10}^{7}$ cm, and Different Values of P (in sec.) and ${\xi }_{j}$

P	${\xi }_{j}=0.25$	${\xi }_{j}=2$
0.033	6.6	5.5
0.33	7.9	6.6

Download table as:  ASCIITypeset image

The primary particle energy has a weak dependence on pulsar period, inclination angle (via ${\xi }_{j}$), and magnetic field strength; its strongest dependence is on the radius of curvature of magnetic field lines. These trends are consequences of the fact that the potential drop across the acceleration gap is regulated by pair formation. The gap terminates when particles reach energies high enough to emit pair-producing photons. A gap with a weak accelerating electric field due to, e.g., a weaker magnetic field and/or longer period and/or smaller current density ${j}_{{\rm{m}}}$ will have a larger height than a gap with a strong accelerating field. A larger height of the gap also results in longer distances traveled by photons; this largely alleviates the dependence of the energy of the pair-producing photons on the magnetic field strength, leaving the curvature of the magnetic field lines as the strongest factor in determining the energy of primary particles.

Now we can make a simple estimate of the maximum cascade multiplicity during a burst of pair creation. As discussed in Section 2, in a hypothetical ideal cascade, all of the kinetic energy of the primary particle is divided into the energies of the pairs produced by the photons with energies just above the escape energy; in such a cascade, the multiplicity is given by Equation (1)—twice the energy of the primary particle divided by the energy of the escaping photons. In Figure 9, we show the estimates for the multiplicity of an ideal cascade $\mathrm{log}{\kappa }_{\max }$ as a function of the magnetic field strength and the radius of curvature of the magnetic field lines for three sets of gap parameters $(P[{\rm{s}}],{\xi }_{j})$: (0.033, 2), (0.033, 0.25), and (0.33, 2). The energies of the primary particles and escaping photons are calculated according to Sections 3.3 and 4. The maximum cascade multiplicity is not very sensitive to pulsar period and inclination angle (via ${\xi }_{j}$); the strongest dependence is on the magnetic field strength. The maximum value of ${\kappa }_{\max }$ is about $3\times {10}^{6}$, which is the absolute upper limit on the PC cascade multiplicity.

Zoom In Zoom Out Reset image size

The analytical expression for the maximum cascade multiplicity can be obtained using expressions for ${\epsilon }_{\pm ,\mathrm{acc}}$ and ${\epsilon }_{\gamma ,\mathrm{esc}}$ from Sections 4 and 3.3. Substituting Equations (19) and (17) into Equation (1), we get an estimate of the upper limit of the cascade multiplicity,

$$\begin{eqnarray}\begin{array}{rcl}{\kappa }_{\max } & = & 5.7\times {10}^{4}\,{P}^{-1/7}{\rho }_{{\rm{c}},7}^{-3/7}{B}_{12}^{6/7}\\ & & \times \,{\chi }_{a,\mathrm{acc}}^{2/7}\,{\xi }_{j}^{1/7}{\chi }_{a,\mathrm{esc}}^{-1}\left(\displaystyle \frac{{s}_{\mathrm{esc}}}{0.5}\right).\end{array}\end{eqnarray} \tag{ 24 }$$

The weak dependence of ${\kappa }_{\max }$ on the pulsar period and inclination angle (via ${\xi }_{j}$) is evident from this formula; this is a consequence of the weak dependence of ${\epsilon }_{\pm ,\mathrm{acc}}$ on these parameters. The strong dependence of ${\kappa }_{\max }$ on the magnetic field strength is due to the strong dependence of ${\epsilon }_{\gamma ,\mathrm{esc}}$ on B. Using values for ${\chi }_{a,\mathrm{acc}}$, ${\xi }_{j}$ assumed in Section 4, and the approximation for ${\chi }_{a,\mathrm{esc}}$ given by Equation (16), we get two final expressions for ${\kappa }_{\max }$, valid for $B\lesssim 3\times {10}^{12}$ G,

$$\begin{eqnarray}&&{\kappa }_{\max }=5.4\times {10}^{5}{\rho }_{{\rm{c}},7}^{-3/7}{P}^{-1/7}{B}_{12}^{6/7}\end{eqnarray} \tag{ 25 }$$

and for $B\gtrsim 3\times {10}^{12}$ G,

$$\begin{eqnarray}&&{\kappa }_{\max }=1.6\times {10}^{6}{\rho }_{{\rm{c}},7}^{-3/7}{P}^{-1/7}{B}_{12}^{-1/7}.\end{eqnarray} \tag{ 26 }$$

For higher magnetic field strengths and smaller radii of curvature of magnetic field lines, the energy of the primary particles is larger and the energy of escaping photons smaller. The energy available for the cascade, and hence, the maximum cascade multiplicity, increases toward higher B and lower ${\rho }_{{\rm{c}}}$ values. This dependence on B saturates at $B\sim 3\times {10}^{12}{\rm{G}}$ because photon absorption at higher field strengths will happen near the pair formation threshold, limiting the decrease of the escaping photons' energy.

In real pulsar cascades, the multiplicity will be (substantially) smaller than ${\kappa }_{\max }$ mainly because (i) not all of the kinetic energy of primary and secondary particles is transferred to pair-producing photons, (ii) the last-generation photons have energies above the pair formation threshold, and (iii) pair production is intermittent; no pairs are produced during the quiet cascade phase. The first two issues are related to the physics of the cascade, and we will address them in the next section. The last issue is directly related to the physics of the screening of the electric field and plasma physics in the blob of freshly formed pair plasma; it can be addressed only by means of self-consistent high-resolution simulations such as those in Timokhin (2010) and Timokhin & Arons (2013) and will be the subject of future research. We can provide only very rough estimates on the effect of pair formation intermittency on the effective PC cascade multiplicity.

We can see from the results of this section that the absolute upper limit on the cascade multiplicity in a single burst of pair formation is ${\kappa }_{\max }\lesssim 3\times {10}^{6}$, with the real effective multiplicity being significantly smaller. This already excludes the possibility of extremely high cascade multiplicities $\sim {10}^{6}\mbox{--}{10}^{7}$ assumed in some theories of PWNe and pulsar high-energy emission (e.g., Bucciantini et al. 2011; Lyutikov 2013).

In Paper I, we developed a simple semianalytical cascade model that allowed us to simulate non-branching cascades—when only a single emission process is involved—and used it to explore CR–synchrotron cascades. As we argued in Paper I, such cascades develop in PCs of moderately magnetized pulsars (B ≲ 10¹² G) where SR of secondary particles is the only source of photons creating pairs in the next cascade generation. In this paper, we are interested in extending the range of applicability of our model for higher field pulsars as well as improving its accuracy. Our new model differs from the one in Paper I in two aspects: (i) it applies to cascades with arbitrary emission/absorption processes—cascade branches can be arbitrarily complex—and (ii) it can account for the fact that emission mechanisms can be broadband, and not all emitted photons are able to create pairs.

6.1. General Algorithm

The spectral energy distribution of synchrotron and curvature radiation is broadband, with a significant amount of energy emitted well below the peak energy ${\epsilon }_{\mathrm{peak}}$,

$$\begin{eqnarray}&&F(\epsilon )\propto \displaystyle \frac{\epsilon }{{\epsilon }_{\mathrm{peak}}}\ {\int }_{0}^{\epsilon /{\epsilon }_{\mathrm{peak}}}{K}_{5/3}(\zeta )d\zeta ,\end{eqnarray} \tag{ 27 }$$

where K is the modified Bessel function of the order 5/3. In Paper I, we used a monoenergetic approximation for these processes—all of the energy is emitted as photons with energies ${\epsilon }_{\mathrm{peak}}$. In our current model, we divide the spectrum into three spectral bins, $\{[0,0.3{\epsilon }_{\mathrm{peak}}],[0.3{\epsilon }_{\mathrm{peak}},1.5{\epsilon }_{\mathrm{peak}}],[1.5{\epsilon }_{\mathrm{peak}},\infty ]\}$.¹¹ CR and synchrotron emission of particles are modeled as emission of photons in each of the three spectral bins $j=1,2,3$ with energies

$$\begin{eqnarray}&&{\epsilon }^{j}\equiv {f}_{\epsilon }^{j}\ {\epsilon }_{\mathrm{peak}},\end{eqnarray} \tag{ 28 }$$

where the number of photons emitted in each spectral bin is equal to the energy emitted by the particle in that bin, ${W}^{i}={f}_{w}^{i}\,W$ (W is the total emission rate), divided by the characteristic energy of the photons,

$$\begin{eqnarray}&&{n}^{j}=\displaystyle \frac{{W}^{j}}{{\epsilon }^{j}}=\displaystyle \frac{{f}_{w}^{j}}{{f}_{\epsilon }^{j}}\displaystyle \frac{W}{{\epsilon }_{\mathrm{peak}}}.\end{eqnarray} \tag{ 29 }$$

The coefficients $({f}_{\epsilon }^{j},{f}_{w}^{j})$ for the energy bins we used are $\{(0.3,0.152)$, $(1,0.518)$, and $(1.5,0.33)\};$ they are calculated by integrating the spectral energy distribution $F(\epsilon )$, Equation (27), over the spectral bins.

In our algorithm, both leptons and photons are macroparticles; the statistical weight of each particle is the number of real particles it represents. We start by calculating the energy of the primary particle accelerated in the gap according to Section 4 and follow this particle as it moves along magnetic field lines, losing energy-emitting CR photons. Each CR photon initiates an electron–positron cascade with secondary particles emitting the next generation of pair-producing photons via SR and RICS of soft X-ray photons from the NS surface. We follow every generation of photons until their energy falls below the escaping energy according to Section 3.3 and compute the number of pairs created by each cascade generation. The diagram in Figure 10 shows the chain of physical processes initiated by a single CR photon. In Figure 10, "rows" represent different cascade generations (particles with the same number of iterations of particles/photons before their creation, starting with primary particles), while "columns" correspond to branches (particles of the same generations produced by different emission processes). In each generation of the cascade, pairs ${e}^{\pm }$ produce photons which are then turned into pairs of the next generation: ${\gamma }_{{\rm{cr}}}$ via CR (shown by solid lines with a double arrowhead, labeled "1"), ${\gamma }_{\mathrm{syn}}$ via SR (shown by solid lines with a single arrowhead, labeled "2"), and ${\gamma }_{\mathrm{RICS}}$ via RICS (shown by dashed lines with an arrowhead, labeled "3"). Numbers in parentheses show the origin of each particle; for example, (1, 3, 2) means that this pair was produced by a synchrotron photon (second generation) emitted by a pair produced by an RICS photon (first generation) created by a CR photon (zeroth generation).

Zoom In Zoom Out Reset image size

Our algorithm is described in detail in Appendix C; here we give a brief overview. The central part of our algorithm is the recursive function PairCreation(), Algorithm 2 (Figure 26). For each photon, PairCreation() calculates whether and where it will be absorbed to create a pair. The photon is counted as absorbed if its mfp is less than the escaping distance, ${\lambda }_{\gamma \to \pm }\leqslant {s}_{\mathrm{esc}}{R}_{{\rm{ns}}}$, and its absorption point x is still inside the cascade zone, $x\leqslant {s}_{\mathrm{cascade}}{R}_{{\rm{ns}}}$. Then, for each emission process, it calculates the energy of the next-generation photons emitted by the pair by calling function emissionFun() for each emission process. Finally, it recursively calls itself for each of the next-generation photons. We follow the primary particle as it moves along magnetic field lines, losing energy and emitting photons via CR—Algorithm 1 in Appendix C (Figure 25). For each CR photon, PairCreation() is called and through its successive recursive calls, the algorithm follows every branch of the cascade.¹² The total cascade multiplicity is calculated by integrating the number of particles produced in cascades generated by each CR photon (computed by recursive calls of PairCreation()) over the distance within the cascade zone. We assume that the size of the cascade zone is equal to the NS radius ${R}_{{\rm{ns}}}$, so ${s}_{\mathrm{cascade}}=1$.

6.2. Microscopic Processes

We analyzed the microphysics of PC cascades of young energetic pulsars in Paper I in great detail. Here we give a brief overview of how we treat the cascade microphysics.

At the distance s after exiting the acceleration zone (hereafter all distances are normalized to ${R}_{{\rm{ns}}}$), the energy of the primary particle is (Equation (19) in Paper I)

$$\begin{eqnarray}&&{\epsilon }_{\pm }(s)={\epsilon }_{\pm }^{0}{\left[1+3H\displaystyle \frac{{\left({\epsilon }_{\pm }^{0}\right)}^{3}}{{\rho }_{{\rm{c}}}^{2}}s\right]}^{-1/3},\end{eqnarray} \tag{ 30 }$$

where ${\epsilon }_{\pm }^{0}$ is the initial particle energy, $H=(2/3){R}_{{\rm{ns}}}{r}_{e}\,\approx 1.88\times {10}^{-7}\,{{\rm{cm}}}^{2}$, and ${r}_{e}={e}^{2}/{m}_{e}{c}^{2}$ is the classical electron radius. While traveling along a segment of the length ds, the energy emitted by the primary particle via CR (all energy quantities are normalized to m_ec²) is

$$\begin{eqnarray}&&{W}_{{\rm{cr}}}=\displaystyle \frac{3}{2}\,{\alpha }_{f}\displaystyle \frac{{\unicode{x0019B}}_{{\rm{c}}}{R}_{{\rm{ns}}}}{{\rho }_{{\rm{c}}}^{2}}\,{\epsilon }_{\pm }^{4}(s){ds}.\end{eqnarray} \tag{ 31 }$$

The peak energy of CR radiation photons is

$$\begin{eqnarray}&&{\epsilon }_{{\rm{CR}},\mathrm{peak}}=\displaystyle \frac{3}{2}\displaystyle \frac{{\unicode{x0019B}}_{{\rm{c}}}}{{\rho }_{{\rm{c}}}}{\epsilon }_{\pm }{}^{3}.\end{eqnarray} \tag{ 32 }$$

The energies and statistical weights of macroparticles representing CR photons are calculated from Equations (28) and (29) using Equations (31) and (32).

We follow the evolution of cascades initiated by CR photons by tracing the pair-producing photons as their energy degrades with each successive cascade generation. The energy of each pair-creating photon ${\epsilon }_{\gamma }$ is transferred to an electron–positron pair, which is always created with a non-zero momentum perpendicular to the magnetic field. The perpendicular energy is emitted by the pair via SR shortly after pair creation. The energy emitted as synchrotron photons is (see Equation (13) in Paper I)

$$\begin{eqnarray}&&{W}_{\mathrm{syn}}={\epsilon }_{\gamma }\left\{1-{\left[1+{\left(\displaystyle \frac{{\chi }_{a}}{b}\right)}^{2}\right]}^{-1/2}\right\},\end{eqnarray} \tag{ 33 }$$

and the peak energy of the synchrotron radiation is

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{syn},\mathrm{peak}}=\displaystyle \frac{3}{4}{\chi }_{a}{\epsilon }_{\gamma },\end{eqnarray} \tag{ 34 }$$

where ${\chi }_{a}$ and b are the values of the parameter ${\chi }_{a}$ and the normalized magnetic field strength at the absorption point of the parent photon, respectively. As for CR, photon energies and statistical weights of macroparticles representing synchrotron photons are calculated from Equations (28) and (29) using Equations (33) and (34).

Pair particles can also scatter thermal photons from the NS surface. If it happens in the cascade zone, the kinetic energy associated with the motion of the particle parallel to the magnetic field is transferred to the next generation of pair-producing photons. The maximum energy that can be emitted as RICS photons is the pair's kinetic energy left after the emission of synchrotron photons,

$$\begin{eqnarray}&&{W}_{{\rm{rics}}}^{0}={\epsilon }_{\gamma }{\left[1+{\left(\displaystyle \frac{{\chi }_{a}}{b}\right)}^{2}\right]}^{-1/2}.\end{eqnarray} \tag{ 35 }$$

The mfp for RICS is given by (Dermer 1990; Sturner 1995; Zhang & Harding 2000)

$$\begin{eqnarray}\begin{array}{rcl}{\lambda }_{{\rm{rics}}} & = & -0.061\,{\epsilon }_{\pm }^{2}\,{T}_{6}^{-1}{B}_{12}^{-2}\times \\ & & {\mathrm{ln}}^{-1}\left[1-\exp \left(-\displaystyle \frac{134{B}_{12}}{{\epsilon }_{\pm }{T}_{6}(1-{\mu }_{s})}\right)\right]\,{\rm{cm}},\end{array}\end{eqnarray} \tag{ 36 }$$

where T₆ is the temperature of the NS surface in units of 10⁶K, B₁₂ is the magnetic field strength in units of 10¹² G, and ${\mu }_{s}=\cos {\theta }_{s}$, where ${\theta }_{s}$ is the angle between the momenta of the scattering photon and particle in the laboratory frame. Equation (36) implicitly takes into account the condition that soft photons must be in cyclotron resonance to be scattered; it is obtained by integrating the resonant cross-section with a blackbody spectrum of target photons (Dermer 1990). If the mfp for RICS is larger than the size of the cascade zone, ${\lambda }_{\mathrm{RICS}}\gt {R}_{{\rm{ns}}}$, we assume that no RICS pair-producing photons are emitted. As particles move away from the NS, the probability of RICS decreases relative to that at the injection point due to the decrease in the magnetic field strength and the number density of soft photons. To account for this effect, we assume that if ${\lambda }_{\mathrm{RICS}}\lt 0.1\,{R}_{{\rm{ns}}}$, all of the pair's kinetic energy left after the emission of synchrotron photons is transferred to RICS photons; this fraction linearly decreases as ${\lambda }_{\mathrm{RICS}}$ is getting bigger $\propto 0.1{R}_{{\rm{ns}}}/{\lambda }_{\mathrm{RICS}}$, so that the energy emitted as RICS photons is

$$\begin{eqnarray}{W}_{{\rm{rics}}}=\left\{\begin{array}{ll}{W}_{{\rm{rics}}}^{0}, & \mathrm{if}\ {\lambda }_{{\rm{rics}}}\leqslant 0.1\,{R}_{{\rm{ns}}}\\ \displaystyle \frac{0.1{R}_{{\rm{ns}}}}{{\lambda }_{\mathrm{RICS}}}\ {W}_{{\rm{rics}}}^{0}, & \mathrm{if}\ 0.1\,{R}_{{\rm{ns}}}\lt {\lambda }_{{\rm{rics}}}\leqslant {R}_{{\rm{ns}}}\\ 0, & \mathrm{if}\ {\lambda }_{{\rm{rics}}}\gt {R}_{{\rm{ns}}}.\end{array}\right.\end{eqnarray} \tag{ 37 }$$

The spectrum of RICS radiation is narrowband, and we approximate this process as emission of monochromatic photons with the energy (see Equation (49) in Paper I)

$$\begin{eqnarray}&&{\epsilon }_{\lambda ,\mathrm{RICS}}={\epsilon }_{\gamma }b{\left[1+{\left(\displaystyle \frac{{\chi }_{a}}{b}\right)}^{2}\right]}^{-1/2}.\end{eqnarray} \tag{ 38 }$$

This number of RICS photons emitted by each pair is

$$\begin{eqnarray}&&{n}_{{\rm{rics}}}=\displaystyle \frac{{W}_{{\rm{rics}}}}{{\epsilon }_{\lambda ,\mathrm{RICS}}}.\end{eqnarray} \tag{ 39 }$$

6.3. Model Applicability

The limitations of our model come from the assumptions used in the derivation of the energy of primary particles and the magnetic field strengths at which physical processes different from the ones we consider here become important. Equation (19) for the energy of primary particles ${\epsilon }_{\pm ,\mathrm{acc}}$ was derived under the assumptions that (i) particles are accelerated freely, i.e., radiation reaction can be ignored, and (ii) the length of the gap is much smaller than the PC radius, so that a one-dimensional approximation can be used. We also do not model cascades for which photon splitting is important, so our model is applicable (iii) for magnetic field strengths below ${B}_{\gamma \to \gamma \gamma }$, calculated according to Equation (14). Constraints (i), (ii), and (iii) together define the range of pulsar parameters where our cascade model is applicable. Constraints (i) and (ii) remain the same as in Paper I; they are derived in Appendices B and C of Paper I correspondingly. Constraint (iii)—the magnetic field strength above which photon splitting becomes important for photons near the threshold of pair formation—depends on the radius of curvature of magnetic field lines ${\rho }_{{\rm{c}}}$. In Figure 11, we show the range of pulsar parameters for which our model is applicable superimposed on the $P\dot{P}$ diagram. The one-dimensional approximation (ii) is valid to the left of the solid line, the approximation (i) of free acceleration above the dotted–dashed line (given by Equation (54) in Paper I). Pulsars with $B\lt {B}_{\gamma \to \gamma \gamma }$ are below the dashed lines. The line with short dashes corresponds to ${\rho }_{{\rm{c}}}={10}^{6}$ cm and ${B}_{\gamma \to \gamma \gamma }=4.4\,\times {10}^{12}$ G, the line with medium dashes to ${\rho }_{{\rm{c}}}={10}^{7}$ cm and ${B}_{\gamma \to \gamma \gamma }=7.7\times {10}^{12}$ G, and the line with long dashes to ${\rho }_{{\rm{c}}}={10}^{8}$ cm and ${B}_{\gamma \to \gamma \gamma }=1.3\times {10}^{13}$ G. The yellow region shows the range of pulsar periods and period derivatives where these three assumptions are valid. We see that most young normal pulsars, including gamma-ray pulsars from the Fermi second pulsar catalog, fall in this range. Technically, the range of pulsar parameters for which our current model is applicable is only slightly different from that of Paper I (the limits on B are less restrictive now; see Figure 13 in Paper I), but our current model offers a considerably better treatment of cascades for B ≳ 10¹² G.

Zoom In Zoom Out Reset image size

For a wide range of pulsar parameters, we computed maximum multiplicities of PC cascades with $\gamma B$ as the pair-creation mechanism, and the CR, SR, and RICS of soft thermal photons from the NS surface as emission mechanisms according to the algorithm described in previous sections.

In Figure 12, we show contours of the cascade multiplicity $\mathrm{log}(\kappa )$ as a function of the magnetic field strength $\mathrm{log}(B)$ in Gauss and the radius of curvature of magnetic field lines $\mathrm{log}({\rho }_{{\rm{c}}})$ in centimeters for an NS with a uniform surface temperature of $T={10}^{6}$ K. The dashed line indicates the parameter space above which photon splitting starts (negatively) to affect the cascade multiplicity. The three plots in Figure 12 are for different pulsar periods $P=33,330$ ms and two different filling factors $\xi =2,0.25$. The electric field in the gap for the case (P = 33 ms, ξ = 2) is an order of magnitude larger than in the other two cases, but the cascade multiplicity is only moderately higher—it is greater by less than ∼2 times. Compared to the simple estimates from Section 5, the total multiplicity for the same pulsar parameters is smaller by the factor of ∼4–5; compare Figure 12 with Figure 9, where ${\kappa }_{\max }$ is plotted for the same combination of parameters. The maximum value for the multiplicity reaches $6\times {10}^{5}$ for smaller radii of curvature of magnetic field lines ${\rho }_{{\rm{c}}}\sim {10}^{6}$ cm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case of a pure CR–synchrotron cascade discussed in Paper I (when the contribution of RICS is ignored), the cascade multiplicity is the highest for magnetic fields around B ≳ 10¹² G and drops for both higher and lower magnetic field strengths (Figure 14 of Paper I). For cascades with RICS, considered in this paper, the multiplicity decrease for higher magnetic fields, B ≳ 10¹² G, is much smaller, because the energy of pair parallel motion is returned back to the cascade by RICS. For lower magnetic fields $B\lt {10}^{12}$ G, the differences in multiplicities between the Figure 12 here and Figure 14 of Paper I vary from a few percent for the case (P = 33 ms, ξ = 2) to ∼70% for the case (P = 330 ms, ξ = 2). These differences are because of the more accurate treatment of emission processes in this paper compared to Paper I.

In Figure 13, we show plots illustrating the dependence of the cascade multiplicity on the NS surface temperature. Contours of $\mathrm{log}(\kappa )$ are plotted for three different temperatures of the NS surface $T=5\times {10}^{5},{10}^{6},{\rm{and}}\,3\times {10}^{6}$ K. The number of soft X-ray photons available for RICS changes dramatically with the temperature and so does the contribution of RICS to cascade multiplicity. For $T=5\times {10}^{5}$ K, the multiplicity profiles are very similar to those of CR–synchrotron cascades. For magnetic fields $B\gtrsim 3\times {10}^{12}$ G, photons are absorbed after propagating a short distance, and pairs are created with small momenta perpendicular to B, which leaves a large fraction of the parent photon's energy in the pairs' parallel motion (this was discussed in detail in Paper I). If there are not enough soft X-ray photons to be upscattered via RICS, this kinetic energy is lost from the cascade, and the total multiplicity diminishes. For higher surface temperatures, the increasing number of soft photons makes RICS increasingly efficient, which leads to the decrease of energy "leaks" from the cascade, and the multiplicity κ becomes similar to the maximum multiplicity ${\kappa }_{\max }$. Indeed, for low temperature $T=5\times {10}^{5}{\rm{K}}$ (left panel of Figure 13), the plot for κ is similar to that of the κ of the CR–synchrotron cascade, and for high temperature $T=3\times {10}^{6}{\rm{K}}$ (left panel of Figure 13), the shape of the κ contours are similar to the ones of ${\kappa }_{\max }$ shown on Figure 9 (left panel).

The cascade efficiency $\kappa /{\kappa }_{\max }$, which characterizes how the energy available to the cascade is converted into pairs, is shown in Figure 14, for the same parameters as the κ in Figure 13. For magnetic fields below ∼10¹² G, the efficiency for all three surface temperatures is similar and ≳20%; for these magnetic fields, the RICS contribution is negligible and so the cascade behavior does not depend on the temperature. For stronger magnetic fields, the efficiency increases with the temperature. For high NS temperatures, the cascade efficiency can be as large as 30%. We should note that because of photon splitting (which mostly affects RICS photons), the real cascade efficiency above the dashed lines is smaller than the values shown in Figure 14.

Zoom In Zoom Out Reset image size

To get a better understanding of PC cascades near their peak multiplicities, we consider here several particular cascades and analyze their properties in detail. We consider cascades in PCs of pulsar with P = 33 ms, assume ξ = 2, and use two values of the radius of curvature of magnetic field lines, ${\rho }_{{\rm{c}}}={10}^{7}\,\mathrm{cm}$ and ${\rho }_{{\rm{c}}}={10}^{7.9}\approx 7.94\times {10}^{7}$ cm. For each value of ${\rho }_{{\rm{c}}}$, we analyze cascade properties for three values of the magnetic field strength to illustrate how the RICS contribution changes with B. These examples represent cascades near their highest multiplicities in PCs of pulsars when (i) there is a significant non-dipolar component of the magnetic field (${\rho }_{{\rm{c}}}={10}^{7}$ cm) as well as (ii) pulsars with a nearly unadulterated dipolar field (${\rho }_{{\rm{c}}}={10}^{7.9}$ cm). As we discussed earlier, the dependence of the cascade multiplicity on pulsar parameters is weak, so these examples should be representative for cascades in pulsars with a broad range of periods and filling factors ξ.

We visualize cascade development with three types of plots: (i) the cascade graphs (Figures 15 and 18), the cumulative pair injection rates decomposed according to (ii) emission mechanisms (Figures 16 and 19) and (iii) cascade generations (Figures 17 and 20). These plots show properties of the entire cascade, i.e., all pairs created by a single primary particle as it moves along magnetic field lines and emits CR photons, which initiate multiple individual cascades. The cascade graphs in Figures 15 and 18 are the quantitative representations of the cascade diagram shown in Figure 10. In these graphs, the vertices represent pairs with the same origin—the chain of emission processes which led to the emission of these pairs' parent photons is the same. The area of the circles at the graphs' vertices is proportional to the number of pairs; however, their sizes are consistent only within the same graph, and the sizes of the vertices of different graphs are not related. The cumulative pair injection rates $N(\lt x)$ for a given distance x, shown in Figures 16, 17, 19 and 20, are the total number of pairs created at distances less than x.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

PC cascades have five to eight generations at most. At each cascade generation, pairs emit several photons which split the energy of the particle. This leads to rapid energy degradation through cascade generations, and the cascade dies after several generations. On average, the most contribution comes from generation 4 (see Figures 17, 20).¹³ The total multiplicity does not necessarily increase with the number of cascade branches—despite there being more branches and generations in case (d), the total multiplicity of these cascades is ∼2.5 times lower than that of cascades for case (a).

Cases (a)–(c) and (e)–(f) represent cascades with the same values of ${\rho }_{{\rm{c}}}$ (10⁷ cm and ${10}^{7.9}$ cm, respectively) and decreasing magnetic field strengths. The role of RICS in cascades decreases with decreasing B. For case (a), RICS is responsible for a comparable, and for case (d), an even slightly larger, number of pairs than SR; see Figures 16(a) and (d). As a result, cascades (a) and (d) have more branches than their counterparts with lower magnetic fields ((b), (c) and (e), (f), respectively); see Figures 15 and 18. RICS never becomes a dominating process, but it can contribute an amount of pairs to the cascade comparable to synchrotron radiation for high magnetic fields.

Another illustration of the fact that the number of cascade branches is not directly related to the total multiplicity is provided by the comparison of cascades (d)–(f). The complexity of the cascades changes significantly (due to the diminishing role of RICS), but the total multiplicities of cascades (e) and (f) are only ∼30% and 55% smaller than that of cascade (d). What matters most is the amount of energy available for the cascade, i.e., the relation of ${\epsilon }_{\pm ,\mathrm{acc}}$ and ${\epsilon }_{\gamma ,\mathrm{esc}}$; see Figure 9.

We used different values for the photon escape distance ${s}_{\mathrm{esc}}=0.5$ and the size of the cascade zone ${s}_{\mathrm{cascade}}=1$. The value of ${s}_{\mathrm{esc}}$ has a direct impact on the energy of escaping photons according to Equation (15), and the multiplicity dependence is close to $\propto 1/{s}_{\mathrm{esc}}$. On the other hand, the exact value of ${s}_{\mathrm{cascade}}$, which limits the range where photons can be emitted and absorbed, has little impact on the total multiplicity as long as ${s}_{\mathrm{cascade}}\geqslant {s}_{\mathrm{esc}}$. As is evident from Figures 17 and 20, most pairs are produced at distances $\lt {s}_{\mathrm{esc}}$ from the NS which implies that their parent photons are emitted at distances a few time smaller than ${s}_{\mathrm{cascade}}$. The contribution from photons emitted at distances approaching ${s}_{\mathrm{cascade}}$ is rather small—primary particles have lost a large fraction of their energy, and resulting cascades have only one to two generations.

Synchrotron photons have mostly ⊥ polarization and RICS ∥ one. The photons susceptible to splitting are the ones with ⊥ polarization. At the field strength where photon splitting becomes important, RICS photons provide a comparable or slightly larger number of pairs than SR photons, hence the splitting should significantly affect cascade multiplicity. The highest cascade multiplicity is reached for B and ${\rho }_{{\rm{c}}}$ values near the dashed lines in Figures 12 and 13. So, the pulsars with the highest multiplicities should have $B\sim 4.4\times {10}^{12}-1.4\,\times {10}^{13}$ G, depending on ${\rho }_{{\rm{c}}}$, and surface temperature T ≳ 10⁶ K. On the $P\dot{P}$ diagram, Figure 11, the pulsars with the highest maximum multiplicity are those near the dashed lines.

The spatial distribution of the pair injection rate shows that most of the pairs produced by RICS photons are created at distances comparable to ${R}_{{\rm{ns}}}$, which are much larger than the PC size ${r}_{\mathrm{pc}}\approx 1.45\times {10}^{4}{P}^{-1/2}$ cm. If soft X-ray photons are emitted from the entire surface of the NS, as assumed in our model, their number density does not decrease dramatically throughout the cascade zone. However, if soft photons are emitted from a hot PC, and the rest of the NS surface is cold, $T\lt {10}^{6}$ K, the number density of soft photons at distances $\gg {r}_{\mathrm{pc}}$ will be much smaller than assumed here. In that case, the role of RICS is reduced, and the PC cascade will operate in the CR–synchrotron regime. In the latter case, the multiplicity will reach its peak at $B\sim {10}^{12}{\rm{G}}$ (for a detailed analysis of CR–synchrotron cascades, see Paper I).

In our previous paper, Paper I, we limited ourselves to CR–synchrotron cascades, which was an adequate approximation for most young energetic pulsars. However, right where CR–synchrotron cascades reach their highest multiplicity, RICS becomes an important emission mechanism, and in order to get an accurate limit on the maximum cascade multiplicity, it must be taken into account. In this study, we included in our model all three processes leading to the emission of pair-producing photons in PCs of energetic pulsars—RICS, SR, and CR—and considered the effect on photon splitting on the cascade multiplicity. The treatment of the radiation is improved, dividing the spectrum into three energy bands instead of the delta-function approximation used in Paper I. We used a more accurate prescription for the single photon pair production which takes into account the decrease of the attenuation coefficient near the pair formation threshold—an important correction for pair formation in magnetic fields with B ≳ 10¹² G. We also used a more consistent treatment of the particle acceleration by finding both the energy and the primary particles ${\epsilon }_{\pm ,\mathrm{acc}}$ and parameter ${\chi }_{a,\mathrm{acc}}$, which regulate pair injection and termination of the acceleration zone, simultaneously. We developed a new semianalytical algorithm, which can incorporate an arbitrary number of microphysical processes, model spatial evolution of cascades and allow fast exploration of cascade parameter space. The improvements upon the model presented in Paper I allowed us to conduct a reasonably accurate study of PC cascades in the regime where they reach their highest multiplicity. Our current model includes the most important microphysical processes relevant for PC cascades in young energetic pulsars with the highest pair yield.

The goal of our study was to find the upper limit on the multiplicity of electron–positron cascades in pulsars. We have performed a systematic study of pair cascades above pulsar PCs for a variety of input parameters including surface magnetic field, pulsar rotation period, primary particle energy, magnetic field radius of curvature, and the temperature of the NS surface. We used the modern description for particle acceleration derived from self-consistent models of the PC acceleration zone, i.e., those that are capable of generating currents consistent with global models of the pulsar magnetosphere. In our model, we do not address directly the non-stationary nature of pulsar cascades. We considered pair cascades generated by primary particles accelerated at the peak of the pair formation burst. The intermittency of the pair formation process reduces the total pair yield, and by studying cascades generated by the most energetic primary particles, we achieved our goal of finding the limit on pair cascade multiplicity.

We find that pair multiplicity is maximized for pulsars with hot $T\geqslant {10}^{6}$ K surfaces. These must be very young pulsars that have not yet cooled down. For such pulsars, cascade multiplicity (almost) monotonically grows with increasing B and decreasing ${\rho }_{{\rm{c}}}$ until photon splitting becomes more important than pair production, which happens first in the last cascade generation. In young hot pulsars, pair cascades reach their highest multiplicities near magnetic field strengths where photon splitting becomes more efficient than single photon pair production for the last-generation pair-producing photons. The maximum multiplicity is in the range $\kappa \sim {10}^{6}-3\times {10}^{5}$ for magnetic field strengths $3\times {10}^{12}\,{\rm{G}}\lesssim \,B\,\lesssim {10}^{13}\,{\rm{G}}$ depending on the radius of curvature of magnetic field lines—the lower and upper limits on B are for $\rho ={10}^{6}$ cm and 10⁸ cm, respectively. For older pulsars, whose surfaces have cooled down below 10⁶ K, the maximum cascade multiplicity is in the range $\kappa \sim 5\times {10}^{5}-{10}^{5}$, and it is achieved at $B\sim {10}^{12}$ G. Even if old pulsars have hot PCs, the density of soft photons at distances comparable to the NS radius will be too small to sustain efficient RICS, and so the cascade operates in the CR–synchrotron regime even for high magnetic field strengths.

Here we ignored geometrical effects caused by the curvature of magnetic field lines. For the smallest values of ${\rho }_{{\rm{c}}}$, the magnetic field lines at large distances from the NS within the cascade zone can bend rather significantly. Such bending causes the displacement of particles in the lateral direction, which, however, would have a negligible effect on the cascade multiplicity in pulsars with hot surfaces because it does not affect CR and SR, and the variation of the incident angle of incoming thermal photons is washed out by the large solid angle these photons are coming from. It would only affect the lateral spreading of the cascade. In long-period pulsars, where the NS surface is cold and the only source of the soft X-rays is the hot PC, the effects of field line bending might potentially increase the multiplicity of the cascade. The pairs' momenta in that case would have a larger angle with soft photons that could increase the cross-section to photon scattering. However, the efficiency of RICS as an emission process is based in part on the wide range of angles thermal photons are coming from. Due to the large range of photon incident angles, pairs within a wide range of energies can still scatter thermal photons (with a relatively narrow energy distribution) in the resonant regime—pairs of different energies scatter photons coming from different directions (Dermer 1990). In the case of hot PCs, the range of photon incident angles will be small (e.g., compared to the case of the hot NS surface), and the potential increase of the scattering cross-section could benefit pairs of a single generation at best, thus making this effect of little importance.

The multiplicity at the peak of the cascade cycle is not very sensitive to the pulsar period, magnetic field, and radius of curvature of the magnetic field lines. The multiplicity varies by less than ∼an order of magnitude for the range of pulsar parameters spanning two or more orders of magnitude. The reason for this is the self-regulation of the accelerator by pair creation: for pulsar parameters resulting in more efficient pair production, the size of the acceleration gap is smaller, and the primary particle energy is lower and vice versa.

Even the most efficient cascades typically have only several generations. High multiplicity is achieved because particles of each generations emit multiple pair-producing photons. The biggest contribution comes from individual cascades with three to four generations. RICS can play an important role in PC cascades; it can provide an even larger number of pairs than SR, although SR never becomes a negligible process. When RICS is an important process, the cascade can have many branches, but the multiplicity does not directly depend on cascade complexity (number of branches and generations). A more important factor is the energy available for the cascade process; the exact way of how this energy is distributed among the final pair population, i.e., via synchrotron or RICS branches, plays a secondary role.

The main factor determining the total yield of PC cascades is, however, the variation of the flux of primary particles due to the intermittency of the particle acceleration in time-dependent cascades. Simple estimates for the "duty cycle" of the particle acceleration, presented in Section 10 of Paper I, predict that the total pair yield of the PC cascade will be lower than the multiplicity values of cascades at the peak of the pair formation burst by a factor ranging from a ∼few for the case of space charge-limited flow (free particle extraction from the NS) with the super GJ current density ($j/{j}_{{\rm{gj}}}\gt 1$) and up to ∼several hundred for the case of Ruderman–Sutherland gaps (no-particle extraction from the PC) and space charge-limited flow with anti-GJ current density ($j/{j}_{{\rm{gj}}}\lt 0$). The effective pair multiplicity then cannot exceed a $\sim {\rm{few}}\times {10}^{5}$ even under the most favorable conditions: in young hot pulsars with a high magnetic field $B\sim 3\times {10}^{12}\,{\rm{G}}$ and a significant non-dipolar component of the magnetic field in PCs so that ${\rho }_{{\rm{c}}}\lesssim {10}^{7}$ cm.

An accurate estimate of the duty cycle requires self-consistent modeling of particle acceleration with a much higher numerical resolution than was done in Timokhin & Arons (2013) and Timokhin (2010), which will be a subject of a future paper. The results of this paper can then be easily adopted into a consistent model of pair supply in pulsars by scaling the multiplicity values obtained here by the duty cycle of particle acceleration. The semianalytical model presented here can also guide more accurate (and time-consuming) numerical simulations of cascades.

In terms of the direct astrophysical implication of this work, our main message is simple—under no circumstances can the pair yield of pulsars be greater than $\sim {\rm{several}}\times {10}^{5}{n}_{{\rm{gj}}}$. This should be taken into account, for example, in modeling PWNe and lepton components of cosmic rays.

This work was supported by NSF grant 1616632 as well as the Chandra and Fermi Guest Investigator programs.

Pairs can scatter soft X-ray photons emitted by the NS surface in the resonant as well as non-resonant regime. However, as we show below, the non-resonant Inverse Compton scattering (NRICS) is a very inefficient emission mechanism and can be ignored in comparison with scattering in the resonant regime.

An emission process could play a role in PC cascades if the distance over which a particle loses a substantial part of its energy to that emission process is less than the size of the cascade zone, which in this work is assumed to be ${R}_{{\rm{ns}}}$. In Figure 21, we plot contours of the mfp ${\lambda }_{\mathrm{NRICS}}$ (in centimeters) of an electron/positron to NRICS of soft photons emitted by the NS surface as a function of the particle's energy and the temperature of the NS. The mfp was calculated by integrating the full ICS cross-section over the non-isotropic distribution of photons emitted by the NS surface, and photons are coming in the solid angle that is centered around the particle's momentum and limited by $(0\leqslant \theta \lt {\theta }_{\max };0\leqslant \phi \leqslant 2\pi )$, with $\cos {\theta }_{\max }=0.5$ (in this paper, we use the same solid angle for modeling RICS). The spectral energy distribution of thermal photons was modeled as the Rayleigh–Jeans power law with the high-energy cutoff chosen in such a way that the total emitted energy is consistent with the Stefan–Boltzmann law.¹⁴ It is easy to see from that plot that the mfp to NRICS becomes less than the NS radius only for very high temperatures of the NS surface, $T\gtrsim 2\times {10}^{6}$ K, and even for $T\gtrsim 3\times {10}^{6}\,{\rm{K}}$ the mfp is only ${\lambda }_{\mathrm{NRICS}}\simeq 0.2\,{R}_{{\rm{ns}}}$. According to common models of NS cooling, even the youngest pulsars should have surface temperatures less than $3\times {10}^{6}\,{\rm{K}}$ (e.g., Haensel et al. 2007).¹⁵ This implies that it is highly unlikely that in PCs of pulsars, particles could lose any significant fraction of their energy to NRICS.

Zoom In Zoom Out Reset image size

But even if the NS surface is very hot, at the upper limit of the predicted temperature range, $T\simeq 3\times {10}^{6}$ K, NRICS would still be of very limited relevance for cascade physics. The reason is as follows. From Figure 21, it is clear that NRICS might be relevant for particles with energies ${10}^{3}\lesssim {\epsilon }_{\pm }\lesssim {10}^{4}$. Let us compare now the efficiency of resonant and non-resonant ICS. In Figure 22, we plot contours of the logarithm of the ratio of the mfp of a particle to non-resonant and resonant ICS ${\lambda }_{\mathrm{NRICS}}/{\lambda }_{\mathrm{RICS}}$ as a function of particle energy and the strength of the magnetic field for the surface temperature $T=3.5\times {10}^{6}$ K. In the energy range ${10}^{3}\lesssim {\epsilon }_{\pm }\lesssim {10}^{4}$, where a significant part of the particle's momentum could be radiated via NRICS within the cascade zone, ${\lambda }_{\mathrm{NRICS}}$ is less than ${\lambda }_{\mathrm{RICS}}$ only for magnetic fields $B\lesssim 3\times {10}^{12}$ G. For magnetic field strengths $\lesssim {10}^{12}\,{\rm{G}}$, the fraction of the parent photon's energy going into the parallel momentum of the created pair is smaller than $\sim 30 \%$—most of the energy is emitted as synchrotron photons; see Paper I, Section 4, Figure 5. In the narrow range of magnetic field strengths ${10}^{12}\,{\rm{G}}\lesssim B\lesssim 3\times {10}^{12}\,{\rm{G}}$, NRICS might become an important process, but only for particles of a single cascade generation. Indeed, because of a very steep dependence of ${\lambda }_{\mathrm{NRICS}}/{\lambda }_{\mathrm{RICS}}$ on particle energy for a given value of B, even if pairs of some generation with energy in the range ${10}^{3}\lesssim {\epsilon }_{\pm }\lesssim {10}^{4}$ do create photons more efficiently in the non-resonant regime, the next generation of pairs will scatter photons more efficiently in the resonant regime.

Zoom In Zoom Out Reset image size

In Paper I, we adopted the widely used Erber (1966) formula for the opacity for single photon pair production in a strong magnetic field,

$$\begin{eqnarray}&&{\alpha }_{\gamma \to \pm }({\epsilon }_{\gamma },\psi )=0.23\,\displaystyle \frac{{\alpha }_{f}}{{\unicode{x0019B}}_{{\rm{c}}}}\ b\,\sin \psi \,\exp \left(-\displaystyle \frac{4}{3\chi }\right),\end{eqnarray} \tag{ 40 }$$

where $b\equiv B/{B}_{q}$ is the local magnetic field strength B normalized to the critical quantum magnetic field ${B}_{q}=e/{\alpha }_{f}{\unicode{x0019B}}_{{\rm{c}}}^{2}\,=4.41\times {10}^{13}$ G, ψ is the angle between the photon momentum and the local magnetic field, ${\alpha }_{f}={e}^{2}/{\hslash }c\approx 1/137$ is the fine structure constant, and ${\unicode{x0019B}}_{{\rm{c}}}={\hslash }/{mc}=3.86\times {10}^{-11}$ cm is the reduced Compton wavelength. The parameter χ is defined as

$$\begin{eqnarray}&&\chi \equiv \displaystyle \frac{1}{2}\,{\epsilon }_{\gamma }b\,\sin \psi ,\end{eqnarray} \tag{ 41 }$$

where ${\epsilon }_{\gamma }$ is the photon energy in units of m_ec². Expression (40) has been obtained in the asymptotic limit of $\chi \ll 1$ and $b\ll 1$. While $\chi \ll 1$ is a good approximation for gamma-rays absorbed in PCs in all pulsars (see Section 3.1), the approximation $b\ll 1$ can become too restrictive for pulsars with higher magnetic fields. For pairs created near the kinematic threshold,

$$\begin{eqnarray}&&{\epsilon }_{\gamma }\sin \psi =2.\end{eqnarray} \tag{ 42 }$$

Equation (40) can overestimate the opacity by a factor of a few for pulsars with magnetic fields $\gtrsim 3\times {10}^{12}\,{\rm{G}}$ (Daugherty & Harding 1983). The discrepancy becomes larger for stronger fields. An accurate treatment of the pair-creation cross-section for non-small b and/or near-threshold pair creation requires summation over a finite number of cyclotron energy levels of created pairs, which results in unwieldy expressions (like Equation (6) in Daugherty & Harding 1983). For our semianalytical model, such treatment would be an overkill, resulting in unnecessary complication of the model. Instead we use the numerical fit to the exact expression for the opacity suggested by Daugherty & Harding (1983; their Equation (24)),

$$\begin{eqnarray}&&\begin{array}{l}{\alpha }_{\gamma \to \pm }({\epsilon }_{\gamma },\psi )\\ =\,\left\{\begin{array}{ll}0.23\,\displaystyle \frac{{\alpha }_{f}}{{\unicode{x0019B}}_{{\rm{c}}}}\ b\,\sin \psi \,\exp \left(-\displaystyle \frac{4{f}_{\alpha }({\epsilon }_{\gamma },b)}{3\chi }\right), & \mathrm{if}\ {\epsilon }_{\gamma }\sin \psi \geqslant 2\\ 0, & \mathrm{if}\ {\epsilon }_{\gamma }\sin \psi \lt 2\end{array}\right.,\end{array}\end{eqnarray} \tag{ 43 }$$

where

$$\begin{eqnarray}&&{f}_{\alpha }=1+0.42{\left(\displaystyle \frac{1}{2}{\epsilon }_{\gamma }\sin \psi \right)}^{-2.7}{b}^{-0.0038}.\end{eqnarray} \tag{ 44 }$$

The second term in f_α is significant only for pair creation close to the threshold (42). The non-zero part of expression (43) for ${\alpha }_{\gamma \to \pm }$ can be written as

$$\begin{eqnarray}&&\begin{array}{l}{\alpha }_{\gamma \to \pm }({\epsilon }_{\gamma },\psi )=0.23\,\displaystyle \frac{{\alpha }_{f}}{{\unicode{x0019B}}_{{\rm{c}}}}\ b\,\sin \psi \\ \times \,\exp \left(-\displaystyle \frac{4}{3\chi }\right)\exp \left(-0.56\,\displaystyle \frac{{b}^{2.6962}}{{\chi }^{3.7}}\right),\end{array}\end{eqnarray} \tag{ 45 }$$

i.e., it differs from the usual Erber's formula (40) by the exponential term

$$\begin{eqnarray}&&{f}_{\alpha ,1}(b,\chi )\equiv \exp \left(-0.56\,\displaystyle \frac{{b}^{2.6962}}{{\chi }^{3.7}}\right).\end{eqnarray} \tag{ 46 }$$

This term significantly differs from 1 when pair formation occurs close to the threshold.

The optical depths for pair creation by a photon in a strong magnetic field after propagating a distance l is

$$\begin{eqnarray}&&\tau ({\epsilon }_{\gamma },l)={\int }_{0}^{l}{\alpha }_{\gamma \to \pm }({\epsilon }_{\gamma },\psi (x)){dx},\end{eqnarray} \tag{ 47 }$$

where integration is along the photon's trajectory. For photons emitted tangent to the magnetic field line, ${dx}={\rho }_{{\rm{c}}}d\psi$, where ${\rho }_{{\rm{c}}}$ is the radius of curvature of magnetic field lines; the angle ψ is always small so the approximation $\sin \psi \approx \psi$ is very accurate. In our approximation of constant magnetic field, both b and ${\rho }_{{\rm{c}}}$ are constants. From Equation (41), we have $\psi =2\chi /{\epsilon }_{\gamma }b$, and substituting it into Equation (47) we can express the optical depth τ to pair production as an integral over χ,

$$\begin{eqnarray}&&\tau ({\epsilon }_{\gamma },l)={A}_{\tau }\,\displaystyle \frac{{\rho }_{{\rm{c}}}}{{\epsilon }_{\gamma }^{2}b}\,{\int }_{0}^{\chi ({\epsilon }_{\gamma },\psi (l))}\chi \exp \left(-\displaystyle \frac{4}{3\chi }\right){f}_{\alpha ,1}(b,\chi )\ d\chi ,\end{eqnarray} \tag{ 48 }$$

where ${A}_{\tau }\equiv 0.92{\alpha }_{f}/{\unicode{x0019B}}_{{\rm{c}}}\approx 1.74\times {10}^{8}\,{{\rm{cm}}}^{-1}$. If we expand the term ${f}_{\alpha ,1}$ in a Taylor series, the integral in Equation (48) can be represented as a sum of integrals $\int {x}^{\xi }\exp (-x){dx};$ each of those integrals can be integrated analytically.

We expand ${f}_{\alpha ,1}$ in a Taylor series as ${e}^{-x}\approx 1-x\,+{x}^{2}/2-...$ and keep up to the fifth term

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\alpha ,1}(b,\chi )\approx 1-0.56\ \displaystyle \frac{{b}^{2.6962}}{{\chi }^{3.7}}\\ +\,0.1568\ \displaystyle \frac{{b}^{5.3924}}{{\chi }^{7.4}}-0.02927\ \displaystyle \frac{{b}^{8.0886}}{{\chi }^{11.1}}\\ +\,4.0977\times {10}^{-3}\ \displaystyle \frac{{b}^{10.7848}}{{\chi }^{14.8}}-4.5894\times {10}^{-4}\ \displaystyle \frac{{b}^{13.481}}{{\chi }^{18.5}}.\end{array}\end{eqnarray} \tag{ 49 }$$

Only odd-order term expansions are monotonic, and we find that a third-order series is not a good-enough approximation near the threshold. Equation (49) provides a good fit with the minimum number of terms—see Figure 23.

Zoom In Zoom Out Reset image size

Substituting expansion (49) into Equation (48), we get an analytical expression for the optical depth as a series of special functions:

$$\begin{eqnarray}&&\begin{array}{l}\tau (\chi )={A}_{\tau }\,\displaystyle \frac{{\rho }_{{\rm{c}}}}{{\epsilon }_{\gamma }^{2}b}\left[\displaystyle \frac{16}{9}\ {\rm{\Gamma }}\left(-2,\displaystyle \frac{4}{3\chi }\right)\right.\\ -\,0.3434\ {b}^{2.6962}\,{\rm{\Gamma }}\left(1.7,\displaystyle \frac{4}{3\chi }\right)\\ +\,3.3165\times {10}^{-2}\ {b}^{5.3924}\,{\rm{\Gamma }}\left(5.4,\displaystyle \frac{4}{3\chi }\right)\\ -\,2.1354\times {10}^{-3}\ {b}^{8.0886}\,{\rm{\Gamma }}\left(9.1,\displaystyle \frac{4}{3\chi }\right)\\ +\,1.0312\times {10}^{-4}\ {b}^{10.7848}\,{\rm{\Gamma }}\left(12.8,\displaystyle \frac{4}{3\chi }\right)\\ \left.-\,3.9835\times {10}^{-6}\ {b}^{13.481}\,{\rm{\Gamma }}\left(16.5,\displaystyle \frac{4}{3\chi }\right)\right].\end{array}\end{eqnarray} \tag{ 50 }$$

The special function ${\rm{\Gamma }}(a,x)$ is the so-called upper incomplete gamma function, defined as ${\rm{\Gamma }}(a,x)={\int }_{x}^{\infty }{t}^{a-1}\exp (-t){dt}$. There are efficient numerical algorithms for its calculation implemented in numerical libraries and scientific software tools; using Equation (50) for the calculation of the optical depths results in more efficient numerical codes than direct integration of Equation (48).

For $b\ll 1$, when the opacity is given by Erber's formula (40), the optical depths to pair creation is given by the first term in Equation (50),

$$\begin{eqnarray}&&\tau (\chi )=\displaystyle \frac{16}{9}{A}_{\tau }\,\displaystyle \frac{{\rho }_{{\rm{c}}}}{{\epsilon }_{\gamma }^{2}b}\,{\rm{\Gamma }}\left(-2,\displaystyle \frac{4}{3\chi }\right).\end{eqnarray} \tag{ 51 }$$

This expression is a more compact from of the expression for $\tau (\chi )$ from Paper I (Equation (6) in that paper) as

$$\begin{eqnarray}&&\displaystyle \frac{16}{9}\ {\rm{\Gamma }}\left(-2,\displaystyle \frac{4}{3\chi }\right)=\left[\displaystyle \frac{{\chi }^{2}}{2}\left(1-\displaystyle \frac{4}{3\chi }\right){e}^{-\tfrac{4}{3\chi }}-\displaystyle \frac{8}{9}\,{\rm{Ei}}\left(-\displaystyle \frac{4}{3\chi }\right)\right].\end{eqnarray} \tag{ 52 }$$

For $\chi \ll 1$, Equation (51) can be simplified further by expanding it into a Taylor series around $\chi =0$ and retaining the first term,

$$\begin{eqnarray}&&\tau (\chi )=\displaystyle \frac{3}{4}{A}_{\tau }\,\displaystyle \frac{{\rho }_{{\rm{c}}}}{{\epsilon }_{\gamma }^{2}b}\,{\chi }^{3}{e}^{-\tfrac{4}{3\chi }}.\end{eqnarray} \tag{ 53 }$$

Here we show pseudo-codes of the algorithms used to compute cascade multiplicity. For the calculation of ${\chi }_{a}({\epsilon }_{\gamma },B,{\rho }_{{\rm{c}}})$, we computed and stored a table of $1/{\chi }_{a}$ values for a uniformly divided grid $77\times 30\times 20$ in $\mathrm{log}{\epsilon }_{\gamma }\,\times \mathrm{log}B\times \mathrm{log}{\rho }_{{\rm{c}}}$ space, and then used cubic piece-polynomial interpolation to get ${\chi }_{a}$ for parameter values required by expressions used in the algorithms.

Our seminumerical algorithm is built around the data structure we call the "cascade matrix." The cascade matrix contains information about the spatial distribution of the pair injection rate in the cascade ordered by cascade generations and emission processes which lead to the creation of pairs. The spatial distribution of pair injection rates is stored as arrays $\left\langle {\kappa }_{x}\right\rangle$ of (fixed) length nx; in our simulations, we usually use nx = 10. We divide our calculation domain $[0,{s}_{\mathrm{esc}}]$ into nx intervals; each value of ${\kappa }_{i}$ is the number of pairs injected in the interval $[{s}_{i},{s}_{i+1}]$.

The structure of a cascade matrix for a cascade initiated by CR when both synchrotron and RICS photons create secondary pairs is shown on Figure 24. Row ${i}_{\mathrm{gen}}$ contains distributions for all pairs created in cascade branches of generation ${i}_{\mathrm{gen}}$. Column ${i}_{\mathrm{proc}}$ contains the distributions for pairs created by photons emitted via the same emission mechanism with id = ${i}_{\mathrm{proc}}$, i.e., pairs created in cascade branches ending in id = ${i}_{\mathrm{proc}}$. Element $[{i}_{\mathrm{gen}},{i}_{\mathrm{proc}}]$ of the cascade matrix is an associative array consisting of entries $({\mathtt{orgn}}):\left\langle {\kappa }_{x}\right\rangle$, where the tuple $({\mathtt{orgn}})$ is the pair's "origin" (the label of the cascade branch which leads to the injection of the pairs), and array $\left\langle {\kappa }_{x}\right\rangle$ is the spatial distribution of pairs created by a given cascade branch. Although the tuple $({\mathtt{orgn}})$ is enough to identify the position of the pair distribution regarding cascade generation and the emission process—the tuple's length is the number of the cascade generation and the last entry is the id of the process producing the pair-creating photon—keeping entries sorted according to ${i}_{\mathrm{gen}}$ and ${i}_{\mathrm{proc}}$ makes interpretation and visualization of the simulation results much easier.

Zoom In Zoom Out Reset image size

Mathematical operations on the cascade matrix¹⁶ are defined as element-wise operations on the pairs spatial distribution arrays $\left\langle {\kappa }_{x}\right\rangle$ with the same $({\mathtt{orgn}})$, e.g., for the addition of two arrays ${\left\langle {\kappa }_{x}\right\rangle }^{1}$ and ${\left\langle {\kappa }_{x}\right\rangle }^{2}$, the resulting array is defined as

$$\begin{eqnarray}&&{\left\langle {\kappa }_{x}\right\rangle }^{\mathrm{res}}={\left\langle {\kappa }_{x}\right\rangle }^{1}+{\left\langle {\kappa }_{x}\right\rangle }^{2}:\,{\kappa }_{i}^{\mathrm{res}}={\kappa }_{i}^{1}+{\kappa }_{i}^{2}\,\mathrm{for}\,i=1...{nx}.\end{eqnarray} \tag{ 54 }$$

Neither the position of the element $({\mathtt{orgn}}):\left\langle {\kappa }_{x}\right\rangle$ nor the value of $({\mathtt{orgn}})$ changes. If an expression involving the cascade matrix entry $[i,j][({\mathtt{orgn}})]$ is missing in one of the matrices but not in the other(s), a zero-filled array of length nx is inserted at the position $[i,j][({{\mathtt{orgn}}}_{i})]$ of the missing element and then the operation is performed as in Equation (54), e.g., for the addition of two cascade matrices M¹ and M², each element of the resulting matrix ${M}^{\mathrm{res}}$ is given by

$$\begin{eqnarray}&&\begin{array}{l}{M}^{\mathrm{res}}[i,j][({\mathtt{orgn}})]\\ =\,\left\{\begin{array}{l}{M}^{1}[i,j][({\mathtt{orgn}})]\\ \,+\,{M}^{2}[i,j][({\mathtt{orgn}})] & \mathrm{for}\ \mathrm{all}\ ({\mathtt{orgn}})\mathrm{in\; both}\,{M}^{1}\,\mathrm{and}\,{M}^{2}\\ {M}^{1}[i,j][({\mathtt{orgn}})] & \mathrm{for}\ \mathrm{all}\ ({\mathtt{orgn}})\mathrm{only}\ \mathrm{in}\,{M}^{1}\\ {M}^{2}[i,j][({\mathtt{orgn}})] & \mathrm{for}\ \mathrm{all}\ ({\mathtt{orgn}})\mathrm{only}\ \mathrm{in}\ {M}^{2}\end{array}\right..\end{array}\end{eqnarray} \tag{ 55 }$$

Such data structure is relatively straightforward to implement in modern scripting languages (e.g., as a list of dictionaries in Python or a list of associations in MATHEMATICA).

The algorithm has two structural parts. The main function shown in Figure 25 as Algorithm 1 follows a primary particle as it propagates through the calculation domain and emits CR photons, integrating contributions of cascades initiated by CR photons. This function fills an (initially empty) cascade matrix M with the data from all cascades initiated by CR photons emitted by the primary particle. The calculation domain is divided in N logarithmic segments¹⁷ and the energy of the primary particle ${\epsilon }_{\pm }$ is calculated in each segment according to Equation (30), which takes into account particle energy losses due to CR. ${\epsilon }_{\gamma }$ is used to calculate the spectrum of CR, which is divided into ${n}_{\mathrm{CR}}$ spectral bins. The recursive function PairCreation is called for each spectral bin to calculate cascades initiated by CR photons. Each CR photon belongs to generation 0 of a new cascade, so PairCreation is called with generation number ${i}_{\mathrm{gen}}=0$, photon's origin ${\mathtt{orgn}}=(0)$, and an empty (auxiliary) cascade matrix M_1. Matrix M_1 filled in recursive call(s) of PairCreation is added to the matrix M according to the trapezoidal integration rule.

Zoom In Zoom Out Reset image size

Function PairCreation, shown in Figure 26 as Algorithm 2, fills a cascade matrix M, which is passed to it as an argument, with data about a photon-initiated cascade. Among its arguments, there is a list of functions emissionFun_list; each function in that list calculates the spectrum of the next-generation cascade photons for one emission mechanism. Emission processes relevant for the full cascade are set by the content of this list. In this paper, the list consists of two functions, the first one calculates spectrum of SR, and the second one the spectrum of RICS photons,¹⁸ but it can be extended to include additional processes, such as, e.g., NRICS. Other arguments given to PairCreation contain information about photon(s) initiating the next cascade generation (their energy ${\epsilon }_{\gamma }$, number of such photons n_γ, their origin ${\mathtt{orgn}}$, number of the current cascade generation ${i}_{\mathrm{gen}}$, coordinate of photon emission point x_e), and parameters of the cascade zone.

Zoom In Zoom Out Reset image size

PairCreation calculates the coordinate where the photon will be absorbed. If this point is outside of the calculation domain, control is returned to the calling program; all recursive calls of this function are terminated in this way. If absorption happens inside the domain, a cascade branch is created—the cascade generation counter ${i}_{\mathrm{gen}}$ is increased and a label for the new branch is created by adding the ID of the emission process to the tuple ${\mathtt{orgn}}$ received as a parameter. The position of injected pairs in array $\left\langle {\kappa }_{x}\right\rangle$ is computed, and an array $\left\langle {\kappa }_{x}\right\rangle$ is created for this branch of the cascade. This array is added to the cascade matrix either to an existing array for the same branch, or, if the branch is not yet present in the matrix, a new entry is created. Then, the function proceeds to the next cascade generation. For each emission process, the spectrum of the next-generation photons is computed by interaction over the function list emissionFun_list. A new instance of PairCreation is called for each of the spectral bins in the spectra of the next-generation photons. At this place, the algorithm gives control to the next cascade generation.

The general algorithms can handle a cascade process of any complexity, which depends on the number of emission processes given to it in the list emissionFun_list and is straightforward to parallelize. The resulting cascade matrix is easy to analyze, i.e., to get the total multiplicity of the cascade all entries of the matrix should be summed; to plot the cascade tree as shown in Figures 15 and 18, all arrays $\left\langle {\kappa }_{x}\right\rangle$ are summed and the matrix is traversed row-wise to create the tree, etc.