Abstract

Electron–positron pair production by collision of photons is investigated in view of application to pulsar physics. We compute the absorption rate of individual gamma-ray photons by an arbitrary anisotropic distribution of softer photons, and the energy and angular spectrum of the outgoing leptons. We work analytically within the approximation that 1 ≫ mc2/E > ε/E, with E and ε the gamma-ray and soft-photon maximum energy and mc2 the electron mass energy. We give results at leading order in these small parameters. For practical purposes, we provide expressions in the form of Laurent series which give correct reaction rates in the isotropic case within an average error of ∼ 7 per cent. We apply this formalism to gamma-rays flying downward or upward from a hot neutron star thermally radiating at a uniform temperature of 106 K. Other temperatures can be easily deduced using the relevant scaling laws. We find differences in absorption between these two extreme directions of almost two orders of magnitude, much larger than our error estimate. The magnetosphere appears completely opaque to downward gamma-rays while there are up to ∼ 10 per cent chances of absorbing an upward gamma-ray. We provide energy and angular spectra for both upward and downward gamma-rays. Energy spectra show a typical double peak, with larger separation at larger gamma-ray energies. Angular spectra are very narrow, with an opening angle ranging from 10−3 to 10−7 radians with increasing gamma-ray energies.

1 INTRODUCTION

Electron–positron pair creation by collision of two photons, also called Breit–Wheeler process, is important in a series of astrophysical questions (Ruffini, Vereshchagin & Xue 2010). Amongst them is the filling of recycled pulsar magnetospheres with plasmas.

The cross-section of two-photon-pair creation has been derived in Berestetskii, Lifshitz & Pitaevski (1982). This is a function of the four momentum of both electrons. In pulsar magnetospheres, there is generally a huge reservoir of low-energy photons and a small number of high-energy photons. In order to decrease computational cost compared to pairwise calculations, the cross-section is integrated over the distribution of the low-energy photons. The exact formula for the reaction rate on an isotropic soft-photon background was first derived in Nikishov (1962) to estimate the absorption of gamma-rays by the extragalactic background light. Numerical integration was needed to obtain practical results. In contexts such as active galactic nuclei or X-ray binaries, various formulations and approximations were developed. Approximated analytical expressions were given in Bonometto & Rees (1971) and Agaronyan, Atoyan & Nagapetyan (1983) in the case of an isotropic soft-photon background distribution and averaging over outgoing angles of the produced leptons. The expression of Agaronyan et al. (1983) also applies for a bi-isotropic photon distribution (both strong and weak photon distributions are isotropic) without angle averaging over leptons. In these papers, the authors provide the energy spectrum of the outgoing leptons. An exact expressions in the case of bi-isotropic photon distribution is derived in Boettcher & Schlickeiser (1997), as well as a comparison to the previous approximations that favors the approach in Agaronyan et al. (1983) for its better accuracy.

The standard picture of a pulsar magnetosphere assumes that its inner part is filled with plasma and corotates with the neutron star with angular velocity Ω*. The primary plasma is made of matter lifted from the neutron-star surface by electric fields (Goldreich & Julian 1969). These particles have highly relativistic energies; their motion in the neutron-star magnetic field generates synchrotron and curvature gamma-ray photons. In addition to primary particles, Sturrock (1971) has shown that electron–positron pairs are created in or near the acceleration regions of the magnetosphere. This provides plasma capable of screening the electric field component parallel to the magnetic field. There are two processes of pair creation: two-photon process and one-photon process in the presence of a strong magnetic field. The one-photon process is the most efficient with young and standard pulsars, of which magnetic field is in the range B ∼ 106–108 T (Burns & Harding 1984). The photon–photon pair-creation process can become more important with high-temperature polar caps, and when the magnetic field is below 106 T as in recycled pulsar magnetospheres. Anisotropy of the soft-photon sources is prone to be important as they are expected to come either from the star (hot spots) or from synchrotron radiation in magnetospheric gaps. That is the main reason of our present investigation.

Many detailed studies of pair-creation cascades in pulsar magnetospheres are based only on the one-photon process. This is for instance the case in the recent studies in Timokhin & Harding (2015). Others take the two reactions into account (Chen & Beloborodov 2014; Harding, Muslimov & Zhang 2002).

In numerical simulations of pulsars, the pair-creation rate is generally estimated with simple proxies. For instance, in Chen & Beloborodov (2014), a mean-free path l = 0.2R* is used for the one-photon process, and l = 2R* for the two-photon process. The rate of creation of pairs is not explicited as a function of the electron (or positron) momentum, neither of the local photon background. Instead, pair creations are supposed to be abundant enough to supply electric charges and current densities. The authors write that this approximation is somehow similar to the force-free approximation. In Harding et al. (2002), both one-photon and two-photon processes are taken into account, and the two-photon process is controlled by a mean-free path derived from Zhang & Qiao (1998), where anisotropy is partially taken into account : the energy integral has a lower limit that depends on the angular size of the hot cap providing the soft-photon background. Besides, these authors do not provide spectra for the created pairs although the energy distribution of the outgoing particles are important for the dynamics of pair cascades. A more complete model needs an integration over every local surface element with a threshold that depends on the location of each elementary emitter. This is what the results of the present paper allow us to do within some approximation, together with angular and energy spectra of the outgoing pairs.

Pair creation by two photons is also important in high energy gamma-ray astrophysics. Many papers about gamma-ray bursts and active galactic nuclei refer to Svensson (1987) and the integrated mean-free path in this paper is also based on Nikishov (1962). Actually, spectra of TeV radiation observed from distant (beyond 100 Mpc) extragalactic objects suffer essential deformation during the passage through the intergalactic medium, caused by energy-dependent absorption of gamma-rays interacting with the diffuse extragalactic background light (Nikishov 1962; Gould & Schréder 1966). This effect drastically limits the horizon of the gamma-ray universe, and this has been taken into account in the science case of high-energy gamma-ray observatories (Vassiliev 2000).

In this paper, we revisit the computation of the two-photon pair-creation rate with the aim of dealing with arbitrarily anisotropic soft-photon background distribution. In addition, we give formulas for angle and energy spectra in order to be able to determine in which state pairs are created. After an introduction to the two-photon pair creation equations in Section 2, the integral over the low-energy photons is defined in Section 3. Practical expressions for spectra are derived in Section 4, and applications to the cosmic microwave background and to a hot neutron star are developed in Section 5.2.

2 THE TWO-PHOTON-QUANTUM- ELECTRODYNAMICS REACTION

When not specified, we use a unit system where the speed of light c = 1.

2.1 General formalism

Any quantum-electrodynamics reaction from an initial quantum state |i〉 to an outgoing state |o〉 can be represented as the decomposition on a final states basis {|fk〉} of the evolved state |$\hat{S}{|i\rangle }$|⁠, |$\hat{S}$| being the evolution operator,
\begin{equation} {|o\rangle } = \sum _k {\langle f_k|}\hat{S}{|i\rangle } {|f_k\rangle } \end{equation}
(1)
From that starting point, if one is able to derive the appropriate evolution operator, one can then determine the probability of transition from a given state to any state of the final basis. We are interested in the reaction which yields an electron e and a positron e+ from the encounter of two photons. Common applications take place in a frame where one is ‘strong’, that is high energy, and the other is ‘weak’. Hence, we call them γs and γw, and
\begin{equation} \gamma _{\rm s} + \gamma _{\rm w} \rightarrow e_- + e_+. \end{equation}
(2)

The state of a free photon can be decomposed on a plane-wave basis parametrized by four momentum and polarization. The common assumption is that the effective state of a photon is very well approximated by one plane wave at the time of the encounter. Such a state is not physical in itself, because it cannot be normalized, i.e. it does not belong to the L2 space, or more physically because the Heinsenberg inequality imposes to the wavefunction to be entirely spread through space as a consequence of the perfect determination of momentum. Though, this assumption should be valid over a local four-volume of space–time Vδt, where the interaction through operator |$\hat{S}$| takes place.

Equivalently, free electrons and positrons live on a basis of plane-wave spinors parametrized by a four momentum and a spin. From now on, the leptons are characterized by their charges and their four-momenta |$P_+ = (P^0_+, \vec{p_+})$| and |$P_-= (P^0_-, \vec{p_-})$| and photons by their four-momenta |$K_{\rm s} = (K^0_{\rm s}, \vec{k_{\rm s}})$|⁠, and |$K_{\rm w}= (K^0_{\rm w}, \vec{k_{\rm w}})$|⁠. We consider that their distributions are averaged over spin and polarization, respectively. Following Berestetskii et al. (1982), the Lorentz-invariant cross-section equations are derived in terms of kinematic invariants (also called Mandelstam variables), defined as
\begin{eqnarray} s &=& (P_- - K_{\rm s})^{2} = (P_+ - K_{\rm w})^{2}, \nonumber \\ t &=& (K_{\rm w} + K_{\rm s})^{2} = (P_+ + P_-)^{2}, \nonumber \\ u &=& ((P_- - K_{\rm w})^{2} = (P_+ - K_{\rm s})^{2}. \end{eqnarray}
(3)
The conservation of four momentum writes
\begin{equation} s + t + u = 2m^2, \end{equation}
(4)
where m is the mass of the electron.
The probability dw per unit time of making a pair is
\begin{equation} \mathrm{d}w = \mathrm{d}^2{\sigma } \times j, \end{equation}
(5)
where dσ is the Lorentz-invariant cross-section corresponding to the Feynman diagrams of Fig. 1
\begin{eqnarray} \mathrm{d}^2{\sigma } &=& - \mathrm{d}s 8\pi r_{\rm e}^2 \frac{m^2}{t^2} \times \left[ \left( \frac{m^2}{s-m^2} + \frac{m^2}{u-m^2}\right)^2 \right.\nonumber\\ &&\left. +\left( \frac{m^2}{s-m^2} + \frac{m^2}{u-m^2}\right) \!-\! \frac{1}{4}\left(\frac{s-m^2}{u-m^2} + \frac{u-m^2}{s-m^2}\right) \right]\!, \end{eqnarray}
(6)
where re is the classical radius of the electron1, ds is the differential of the kinematic invariant s at |$P_-^0$| and |$K_{\rm s}^0$| fixed,
\begin{eqnarray} \mathrm{d}s=2\mathrm{d}(\vec{p_-}\cdot \vec{k_{\rm s}}) = 2\left\Vert \vec{p_-} \right\Vert \Vert \vec{k_{\rm s}} \Vert \sin (\vec{p_-},\vec{k_{\rm s}}) \mathrm{d}(\vec{p_-},\vec{k_{\rm s}}) \end{eqnarray}
(7)
and j is the elementary two-particle flux of the reaction,
\begin{equation} j = \frac{1}{V}\frac{K_{\rm s}\cdot K_{\rm w}}{K^0_{\rm s} K^0_{\rm w}}, \end{equation}
(8)
and V is the interaction volume previously defined2. Only the current j is frame-dependent. In particular it reads j = 2/V in the centre of mass (CM) of the reaction.
Figure 1.

Reaction of electron–positron pair creation from a pair of photons represented to first order by Feynman diagrams. Photons have four-momenta Ks and Kw while electron and positron have, respectively, P and P+.

Let us notice that a reaction is possible only if the energy of the two photons exceeds the mass energy of the electron and of the positron. One shows that the kinematic invariant t (3) is equal to the square of the energy in the CM. This allows us to define the frame-invariant criterion
\begin{equation} \sqrt{t} \ge 2m, \end{equation}
(9)
which turns into
\begin{equation} K^0_{\rm s} K^0_{\rm w}\left(1-\cos \xi \right) \ge 2m^2 \end{equation}
(10)
where ξ is the angle between the two photons.

A few important properties of dw can be evidenced by taking a look at cross-section (6) averaged over every possible direction of the outgoing lepton (Berestetskii et al. 1982). As a result, the averaged cross-section depends only on the kinematic invariant τ = t/(4m2) the ratio between the CM energy and the threshold energy. Without loss of generality in the present discussion, we can assume that the reaction takes place in the CM frame, such that the elementary current j = 2/V. The ultra-relativistic limit (Berestetskii et al. 1982) shows that the cross-section vanishes like log τ/τ. This kind of decrease with energy is a common feature of quantum mechanical cross-sections. Moreover, one can numerically estimate the CM energy corresponding to the maximum of the reaction rate to be |$\sqrt{t} \simeq 1.4(2m)$|⁠.

Concerning the angular dependency, leptons are created almost isotropically when the reaction is near threshold while their momenta become aligned with those of the progenitor photons when going to higher energies (see e.g. Berestetskii et al. 1982). In the observer's frame, this translates in a larger angular dispersion for reactions close to threshold.

Equations (38) fully describe the interaction for a given pair of photons; but in a pulsar's magnetosphere, there is a huge amount of photon pairs. In a simulation, it is not possible to compute dw for each pair; we need a statistical approach and a kind of ‘collective’ reaction rate dw. We define it in the next section.

2.2 The pair reactions that count in a pulsar magnetosphere

In a pulsar magnetosphere, the weak photons Kw are mostly caused by the blackbody radiation of the neutron star, or possibly by synchrotron from secondary pair cascades. Their energies range in the X-ray domain. The strong photons are caused by the synchro-curvature radiation of energetic particles (electrons, positrons and possibly ions). Their energies are in the gamma-ray domain. They are more scarce than weak photons. Let's follow a ‘rare’ high-energy, strong photon taken from a phase-space distribution fs. We assume that it flies through an abundant stream of low-energy, weak photons with a distribution fw.

Strong photons negligibly interact with other strong photons because they are not abundant, and because the reaction would likely be far above threshold in equation (10), and therefore inefficient.

Weak photons do not interact with other weak photons since their energies are under the reaction threshold.

Thus, only weak/strong interactions remain, but weak photons are so numerous that a reaction negligibly changes their distribution. Because strong photons are less abundant, pair creations can change their distribution. Hence, for the simulation of a pulsar's magnetosphere, we need to compute the probability of interaction of a strong photon on the background distribution fw of weak photons. Indeed, it does not matter which weak photon is annihilated but we want to update fs as well as the lepton distributions with the outcome of the reactions. With our representation of the involved particles, this amounts to compute the probability dw of creating a lepton of four-momentum P from a photon Ks,
\begin{equation} \mathrm{d}W = \mathrm{d}W_{K_{\rm s} \rightarrow P}. \end{equation}
(11)
For example, one could think of high-energy synchrotron or curvature photon emitted above the polar cap of a pulsar and flying through a stream of thermal photons emitted by the crust. Let us notice that the probabilities of making a positron or an electron are the same, and that a four momentum has four components but only three are independent since ‖P2 = m2c4. These three free parameters can be parametrized by one direction (two parameters) and the energy of the particle.

3 PROBABILITY OF REACTION FOR A GIVEN PHOTON DISTRIBUTION

Quite naturally, the desired probability is the sum over all the possible reactions involving a photon Ks from the background, that would produce an electron at P (respectively, a positron at P+),
\begin{eqnarray} \mathrm{d}W_{f_{\rm w}}(K_{\rm s}, P_-) \!=\!\!\!\! \sum _{L_- \!=\! \lbrace (K_{\rm w}, P_+) / K_{\rm s} \rightarrow P_-\rbrace } \mathrm{d}w(K_{\rm w}, K_{\rm s}, P) \!\times N_{\rm w} \times N_{\rm s}, \end{eqnarray}
(12)
where Nw and Ns are the number of photons of four momentums Kw and Ks, respectively, within the interaction volume V. In spite of greater simplicity in the CM frame, we must use equation (6) in the laboratory frame, because the CM frame would be different for each of the summed pairs of photons. Since low energy photons are parametrized continuously, we must change our sum for an integral, which yields
\begin{eqnarray} N_{\rm w} & \rightarrow & f_{\rm w}(\vec{x_{\rm w}}, \vec{k_{\rm w}}) \mathrm{d}^{3}\vec{x}\mathrm{d}^{3}\vec{k_{\rm w}}, \nonumber \\ \mathrm{d}W_{f_{\rm w}}(K_{\rm s}, P_-) & = & N_{\rm s} \int _{L_-} \frac{c}{V}\frac{K_{\rm s}\cdot K_{\rm w}}{K^0_{\rm s} K^0_{\rm w}} f_{\rm w}(\vec{x}, \vec{k_{\rm w}}) \mathrm{d}^{6}{\Omega }, \nonumber \\ \mbox{ where }\mathrm{d}^{6}{\Omega }& = & \mathrm{d}^2 \sigma \mathrm{d}^{3}\vec{x}\mathrm{d}^{3}\vec{k_{\rm w}}. \end{eqnarray}
(13)
We assume that strong photons are spread out in space such that their density does not vary on the interaction volume V such that their local density is ns = Ns/V. Consequently, the differential probability of interaction per unit time reads
\begin{eqnarray} \mathrm{d}W_{f_{\rm w}}(K_{\rm s}, P_-) &=& n_{\rm s} W_{\rm k}, \nonumber\\ W_{\rm k}&=&\int _{L_-} \mathrm{d}^2 \sigma \frac{c K_{\rm s}\cdot K_{\rm w}}{K^0_{\rm s} K^0_{\rm w}} f_{\rm w}(\vec{x}, \vec{k_{\rm w}}) \mathrm{d}^{3}\vec{k_{\rm w}}, \end{eqnarray}
(14)
where the volume element |$\mathrm{d}V = \mathrm{d}^{3}\vec{x}$|⁠.

3.1 The domain of integration

Let us precisely define the domain L of integration. We note |$\Pi _\alpha = \lbrace P \in \mathbb {R}^4 : \left\Vert P \right\Vert = \alpha ^2c^4, P^0 \ge \left|\alpha \right|\rbrace$| such that Πm is the set of lepton four-momenta (m being the mass of the electron) and Π0 is the set of photon four-momenta. Then,
\begin{eqnarray} L_-(P_-,K_{\rm s}) \!=\! \lbrace (P_+, K_{\rm w}) \in \Pi _{\rm m}\times \Pi _0 : K_{\rm w} \!-\! P_+ = P_- \!-\! K_{\rm s} \rbrace . \end{eqnarray}
(15)
Equivalently, L is the subset of |$\mathbb {R}^4 \times \mathbb {R}^4$| parametrized by Kw with the following constraints:
\begin{equation} \left\lbrace \begin{array}{rclr}P_+ & = & K_{\rm w} - (P_- - K_{\rm s}) & \mathrm{(a)},\\ \left\Vert P_+ \right\Vert ^2 & = & m^2 c^4 & \mathrm{(b)}, \\ \left\Vert K_{\rm w} \right\Vert ^2 & = & 0 & \mathrm{(c)}, \\ K_{\rm w}^0 & \ge & 0 & \mathrm{(d)},\\ P_+^0 & \ge & m c^2& \mathrm{(e)}. \end{array} \right. \end{equation}
(16)
Condition (a) expresses the conservation of four momentum. Conditions (b) and (e) come from |$P_+^0 \in \Pi _{\rm m}$|⁠, and conditions (c) and (d) come from Kw ∈ Π0. We can compute the number of degrees of freedom in L. The set L is a subset of Πm × Π0 of dimension 8. The condition (a) on quadrivectors substracts 4 degrees of freedom. The conditions (b) and (c) both substract 1 degree of freedom. We are left with a set L of dimension 2.
Some of the conditions in equation (16) are already incorporated in the solution of our problem. Condition (c) is already implicitly met in equation (14). Condition (a) is also implicitly met by the set of variables used. Only (b) is not straightforward, since P+ is not directly part of the variables of integration. One can still convert it into a condition on the three other four-vectors by putting (a) into (b) and using (c), the three following equalities being equivalent:
\begin{eqnarray} \left\Vert P_+ \right\Vert ^2 = m^2 c^4 \nonumber \\ \left\Vert K_w - (P_- - K_{\rm s}) \right\Vert ^2 = m^2c^4 \end{eqnarray}
(17)
\begin{eqnarray} K_{\rm w}^0(K_{\rm s}^0 - P_-^0) - \Vert \vec{k_{\rm s}} - \vec{p_-}\Vert K_{\rm w}^0 \cos \xi = K_{\rm s}\cdot P_-, \end{eqnarray}
(18)
\begin{eqnarray} \xi = \mbox{angle}({\vec{k_{\rm s}} - \vec{p_-},\vec{k_{\rm w}}}), \end{eqnarray}
(19)
where |$\Vert \vec{k_{\rm w}}\Vert = K_{\rm w}^0$|⁠. The limit case where |$\vec{k_{\rm s}} = \vec{p_-}$|⁠, for which cos ξ is not defined, is physically impossible because equation (18) would imply |$K_{\rm w}^0<0$|⁠, in contradiction with condition (d). With some algebra, we can show that Ks · P ≥ 0 and that the condition |cos ξ| ≤ 1 imposes |$K_{\rm w}^0 > \epsilon _{{\rm min}}$|⁠, where
\begin{equation} \epsilon _{{\rm min}}= \frac{K_{\rm s}\cdot P_-}{\left\Vert \vec{k_{\rm s}} - \vec{p_-} \right\Vert + K_{\rm s}^0 - P_-^0} \end{equation}
(20)
More precisely |$K_{\rm w}^0([-1,\cos \xi _0[) = [\epsilon _{{\rm min}},+\infty [$| and |$K_{\rm w}^0(\cos \xi > \cos \xi _0) < 0$|⁠.
We can distinguish three regimes of approximation:
\begin{eqnarray} K_{\rm s}^0 >> p_- & \mathrm{ : } & \epsilon _{{\rm min}}\sim \sqrt{m^2 + p_-^2} - p_-\cos \theta , \end{eqnarray}
(21)
\begin{eqnarray} K_{\rm s}^0 \sim p_- & \mathrm{ : } & \epsilon _{{\rm min}}\sim K_{\rm s}^0\sqrt{\frac{1 - \cos \theta }{2}} \nonumber , \\ K_{\rm s}^0 << p_- & \mathrm{ : } & \epsilon _{{\rm min}}\sim p_-. \end{eqnarray}
(22)
For further approximations, we consider that the weak-photon distribution has a cut-off at ε = εmax < m/4.
\begin{eqnarray} k \gg m/4 > \epsilon _{\rm max}=128 \mbox{ {}keV}. \end{eqnarray}
(23)
Because the weak distribution function is in the most extreme case composed of thermal X-rays typically in the range 1–10 keV for a pulsar, this approximation is reasonable.

4 GENERAL SOLUTION

4.1 Energy spectrum

The probability of interaction depends on the integral Wk defined in equation (14). In this section, Wk is directly expressed as a multiple integral with explicit boundaries. The results exposed in this section can be used directly for applications. The path followed to compute them are described in Appendix A. A summary of the notations and useful relations is given in Appendix B. The new expression of Wk involves new variables that appear both in the integrand and in the boundaries of the integral. First, new notations are introduced, for shorter formulas,
\begin{eqnarray} K_{\rm s} & \equiv &(k, \vec{k} = k\vec{z}), \end{eqnarray}
(24)
\begin{eqnarray} P & \equiv & \left(P^0, \vec{p}\right), \end{eqnarray}
(25)
\begin{eqnarray} K_{\rm w} & \equiv & \left(\epsilon , \vec{x} =(x,y,z) \right), \end{eqnarray}
(26)
\begin{eqnarray} \theta & \equiv & \mbox{angle} (\vec{k},\vec{p}). \end{eqnarray}
(27)
With the new notations related to P and to Ks, the integration set L(P, Ks) can be rewritten L(p, cos θ, k). Let Ω be the set of angular components of the electron P, we rewrite Wk as
\begin{equation} W_{\vec{k}} = c \int _{\Omega } \mathrm{d}\Omega \int _{L_-} \frac{\mathrm{d}^2\sigma }{\mathrm{d}\Omega } \frac{K_{\rm s}\cdot K_{\rm w}}{K_{\rm s}^0 K_{\rm w}^0} f_{\rm w}(\vec{k_{\rm w}}) \mathrm{d}^3\vec{k_{\rm w}}. \end{equation}
(28)
We wish to compute the probability of making a pair of which the electron P is in a volume of phase space defined by
\begin{eqnarray*} k/2 < p_1 < & p & < p_2 < k \nonumber \\ (\cos \theta , \phi ) \in \Omega & = & [C_1, C_2 ] \times [0,{2\pi }]\, {\rm with }\,C_{\rm min} \le C_1 < C_2 < 1 \nonumber , \end{eqnarray*}
where p1, p2, C1 and C2 can be set arbitrarily as long as the above inequalities are correct. After the computations exposed in Section A, Wk is transformed into a multiple integral with explicit boundaries. Before showing it, a new set of variables is introduced. The parameter μ parametrizes cos θ,
\begin{equation} \cos \theta = 1 - \left( \frac{2(k-p)}{kp}\mu \epsilon _{\rm max}- \frac{m^2}{2p^2} \right). \end{equation}
(29)
It varies in an interval μ ∈ [μmin, 1], where μ = 1 corresponds to cos θ = cmin and μmin is such that cos θ = 1,
\begin{equation} \mu _\mathrm{min}= {\frac{1}{4}}\frac{km}{p(k-p)} \frac{m}{\epsilon _{\rm max}} \end{equation}
(30)
We define the dimensionless coefficients ai(p),
\begin{eqnarray} a_1(p) & = & -\frac{m^2 \left(k^2-2 k p+2 p^2\right)}{8 k p^2 (k-p)}, \end{eqnarray}
(31)
\begin{eqnarray} a_2(p) & = & -\frac{m^4 \left(k^2-4 k p+2 p^2\right)}{16 \epsilon _{\rm max}k p^3 (p-k)}, \end{eqnarray}
(32)
\begin{eqnarray} a_3(p) & = & \frac{m^6 (3 k-2 p)}{32 \epsilon _{\rm max}^2 p^3 (k-p)^2}, \end{eqnarray}
(33)
\begin{eqnarray} a_4(p) & = & -\frac{k m^8}{64 \epsilon _{\rm max}^3 p^4 (k-p)^2}, \end{eqnarray}
(34)
and
\begin{equation} R = 2\epsilon _{\rm max}\sqrt{\mu (1 - \mu )}. \end{equation}
(35)
In the following, we do not write the p dependance of the ai coefficients except when otherwise stated. It is convenient to express the weak-photon three momentum in cylindrical coordinates |$\vec{x} =(r,\phi _w, z)$|⁠. Then, only the distribution function fw depends on the angle ϕw, which allows a direct integration defining the function (see also A50)
\begin{equation} F_{\rm w}(r, \mu ) = \int _{\phi _{\rm w} = 0}^{2\pi }f_{\rm w}\left(r, \phi _{\rm w}, z(r^2,\mu )\right) \mathrm{d}\phi _{\rm w}, \end{equation}
(36)
where z(r2, μ) is defined in equation (A47) by
\begin{equation} z(r^2,\mu ) = \frac{k}{4}\left(\frac{r^2}{\mu k\epsilon _{\rm max}} - 4 \mu \frac{\epsilon _{\rm max}}{k}\right). \end{equation}
(37)
The integral Wk in (28) is approximated by
\begin{eqnarray} W_{\vec{k}} = c{2\pi }\int _{p_1}^{p_2} \mathrm{d}p \int _{\mu _1}^{\mu _2} \mathrm{d}\mu \sum _{i=1}^{4} \frac{a_i}{\mu ^i} \int _{r = 0}^{R} 2 F_{\rm w}(r, \mu )r \mathrm{d}r. \end{eqnarray}
(38)
Here, the boundaries of the integration domain are left arbitrarily. The reaction probability integrated other every outgoing momenta can be computed as well. In this case, the μ integral is taken from μmin to 1 and p ranges between k/2 and a maximum pmax defined such as μmin(p = pmax) = 1. We find
\begin{equation} p_{\mathrm{max}}= \frac{k}{2}\left(1 + \sqrt{1 - \frac{m^2}{k\epsilon _{\rm max}}}\right) \end{equation}
(39)
The spectrum of outgoing lepton energy is readily obtained as
\begin{eqnarray} &&{\frac{\mathrm{d}W_{\vec{k}} }{\mathrm{d}p}\left(\max \left(p, k-p \right)\right) = c{2\pi }}\nonumber\\ &&\times\int _{\mu _\mathrm{min}}^{1} \mathrm{d}\mu \sum _{i=1}^{4} \frac{a_i}{\mu ^i} \int _{r = 0}^{R} 2 F_{\rm w}(r, p, \mu )r \mathrm{d}r. \end{eqnarray}
(40)

4.2 Angular spectrum

It is also possible to compute the angular spectrum of the outgoing leptons. The problem has to be split in two, whether one consider the higher energy particle (p > k/2) or the lower-energy particle (p < k/2).

For the higher energy particle, one takes equation (38) and changes variable μ to cθ = 1 − cos θ using equation (29). One then obtains
\begin{equation} \frac{\mathrm{d}W_{\vec{k}}}{\mathrm{d}{c_\theta }} = {2\pi }c \int _{p_1}^{p_2} \mathrm{d}p\sum _{i=1}^{4} \frac{a_i^{\prime }}{c^i}\int _{r=0}^{R} 2F_{\rm w}(r,p,{c_\theta }) r\mathrm{d}r, \end{equation}
(41)
where
\begin{equation} a_i^{\prime } = a_i \frac{\mathrm{d}\mu }{\mathrm{d}{c_\theta }} = a_i \frac{kp}{2(k-p)\epsilon _{\rm max}}, \end{equation}
(42)
and the domain of integration has the following limits
\begin{equation} p_{\min }= k/2 \le p_1 \le p_2 \le p_{\mathrm{max}} \end{equation}
(43)
with
\begin{equation} p_{\mathrm{max}}= \frac{\epsilon _{\rm max}k}{m} \frac{1 + \sqrt{1 - \frac{m^2}{\epsilon _{\rm max}k} - \frac{{c_\theta }m^2}{2 \epsilon _{\rm max}^2} } }{2\epsilon _{\rm max}/m + {c_\theta }k/m}, \end{equation}
(44)
obtained by inverting equation (30). The limits for cθ are given by
\begin{equation} 0 \le {c_\theta }\le {c_\theta }_{\max } \end{equation}
(45)
with
\begin{equation} {c_\theta }_{\max } = 2 \left(\frac{\epsilon _{\rm max}}{k} - \frac{m^2}{ k^2}\right), \end{equation}
(46)
for which pmax = k/2 + ◯(1).

In virtue of 37, Fw now depends explicitly on p, hence the dependence in (41).

For the lower energy lepton, we need first to establish the kinematic relation between its outgoing angle defined by |$\cos \theta ^{\prime } = 1 - {c_\theta ^{\prime }}$| and the higher energy lepton variables. All primed quantities refer to the lower energy lepton. Taking the strong-photon direction along the z-axis we have the relation
\begin{equation} 1-{c_\theta ^{\prime }}= \frac{p_z^{\prime }}{p^{\prime }}. \end{equation}
(47)
Using the conservation of momentum (16 a) one can express |${c_\theta ^{\prime }}$| to leading order
\begin{equation} 1 -{c_\theta ^{\prime }}= \frac{1 - \frac{\tilde{p}}{1-\tilde{p}}\frac{\tilde{m}^2}{\tilde{p}^2} - \tilde{p}}{ \left(\left(1 - \tilde{p}- \frac{\tilde{m}^2}{2\tilde{p}}\right)^2 - \tilde{m}^2\right)^{1/2}} +\bigcirc (1), \end{equation}
(48)
where every quantities-tilded quantity is expressed in unit of k, |$\tilde{a} =a/k$|⁠. Since there is no dependence on cθ, one can directly deduce that
\begin{equation} \frac{\mathrm{d}W_{\vec{k}}}{\mathrm{d}{c_\theta ^{\prime }}} =\frac{\mathrm{d}W_{\vec{k}}}{\mathrm{d}p} \left|\frac{\mathrm{d}p}{\mathrm{d}{c_\theta ^{\prime }}} \right|, \end{equation}
(49)
which can be expressed using
\begin{eqnarray} \frac{\mathrm{d}{c_\theta ^{\prime }}}{\mathrm{d}p} \!=\! \frac{2\tilde{m}^2\left(\tilde{m}^4 \!+\! 2 (1-\tilde{p})^2\tilde{p}(3\tilde{p}-1) \!+\! \tilde{m}^2 (1 \!-\! \tilde{p}(\tilde{p}(5-2\tilde{p}) \!+\! 2)\right)}{k(1-\tilde{p})^2\left(\tilde{m}^4-4\tilde{m}^2\tilde{p}\!+\! 4(1-\tilde{p})^2\tilde{p}^2\right)^{3/2}} \nonumber\\ \end{eqnarray}
(50)
after numerical inversion of (48). One finds that |${c_\theta ^{\prime }}$| is a monotonously increasing function of p and that
\begin{eqnarray} {c_\theta ^{\prime }}(k/2) & = & 4\frac{m^2}{k^2}, \end{eqnarray}
(51)
\begin{eqnarray} {c_\theta ^{\prime }}(p_{\mathrm{max}}) & = & \frac{\epsilon _{\rm max}^2}{m^2}\left(1 + \sqrt{1 - \frac{m^2}{k\epsilon _{\rm max}}}\right). \end{eqnarray}
(52)

4.3 With conic boundary conditions

We consider the case where the soft photon distribution is defined everywhere between two cones of axis |$\vec{k}$| and of half-apertures 0 ≤ ξ1 < ξ2 ≤ π, and the distribution Fw (36) is given as a function of the coordinates |$(\bar{\epsilon }=\epsilon /\epsilon _{\rm max}, C_{\xi }=\cos \xi )$|⁠, where we deduce from (A48) and (A45)
\begin{eqnarray} \bar{\epsilon }& = & \frac{\bar{r}^2}{4\mu } + \mu , \end{eqnarray}
(53)
\begin{eqnarray} C_{\xi }& = & 1 - \frac{2\mu }{\bar{\epsilon }}, \end{eqnarray}
(54)
where |$\bar{r}=r/\epsilon _{\rm max}$|⁠.
We are now looking for the appropriate boundary conditions to apply to integral (38). Using the fact that
\begin{equation} \tan \left(\frac{\pi }{2} - \xi \right) = \frac{z(r^2,\mu )}{r}, \end{equation}
(55)
where z(r2, μ) is defined in equation (37), one finds the new boundaries in r by inverting this relation. The resulting r boundaries are given by
\begin{equation} r_{\xi _{1,2}} = 2\epsilon _{\rm max}\mu \left(\frac{1}{\tan \xi _{1,2}} + \frac{1}{\sin \xi _{1,2}}\right), \end{equation}
(56)
where one checks that |$r_{\xi _{2}} < r_{\xi _{1}}$|⁠. We need the intersection |$[r_{\xi _{2}}, r_{\xi _{1}}]\cap [0,R_{\max }]$| which implies solving for |$R_{\max }(\mu _{\xi _{1,2}}) = r_{\xi _{1,2}}$|⁠, which gives us
\begin{equation} \mu _{\xi _{1,2}} = \frac{1}{2} \frac{\sin ^2\xi _{1,2}}{1 + \cos \xi _{1,2}}, \end{equation}
(57)
where one checks that |$\mu _{\xi _{2}} > \mu _{\xi _{1}}$|⁠. We can rewrite the energy spectrum (40) as
\begin{eqnarray} &&{\frac{\mathrm{d}W_{\vec{k}} }{\mathrm{d}p}\left(\max \left(p, k-p \right)\right) = }\nonumber\\ && {2\pi }cr_{\rm e}^2\left( \int _{\mu _\mathrm{min}}^{\min \left(\max \left(\mu _\mathrm{min}, \mu _{\xi _1}\right),1\right)} \mathrm{d}\mu \int _{\max \left(0, r_{\xi _2}\right)}^{r_{\xi _1}}\mathrm{d}r \, \right. \nonumber \\ && \left.+ \int _{\min \left(\max \left(\mu _\mathrm{min}, \mu _{\xi _1}\right),1\right)}^{\min \left(\mu _{\xi _2},1\right)} \mathrm{d}\mu \int _{\max \left(0, r_{\xi _2}\right)}^{R}\mathrm{d}r \right) \sum _{i=1}^{4} \frac{a_i}{\mu ^i} 2 F_{\rm w}(r, \mu )r. \end{eqnarray}
(58)
Concerning the angular spectrum, nothing more needs to be done for lower energy leptons, and for higher energy leptons we proceed similarly as for the energy spectrum above. Starting from (41), one needs to replace μ by its expression as a function of cθ and p in R. This allows us to define the p analogoues of |$\mu _{\xi _{1,2}}$| by
\begin{equation} p_{\xi _{1,2}} = k\frac{1 + \bar{r}_{\xi _{1,2}}^2\frac{\epsilon _{\rm max}}{{c_\theta }k} + \sqrt{1-\bar{r}_{\xi _{1,2}}^2}}{2 + {c_\theta }\frac{k}{\epsilon _{\rm max}} + \bar{r}_{\xi _{1,2}}^2\frac{\epsilon _{\rm max}}{{c_\theta }k}}, \end{equation}
(59)
where one shows that |$p_{\xi _{1}} <p_{\xi _{2}}$|⁠. The spectrum for higher energy leptons is then obtained from (41)
\begin{eqnarray} &&{\frac{\mathrm{d}W_{\vec{k}} }{\mathrm{d}c} = }\nonumber\\ && c{2\pi }\left( \int _{p_{\min }}^{\min \left(\max \left(p_{\min }, p_{\xi _1}\right),p_{\mathrm{max}}\right)} \mathrm{d}p \int _{\max \left(0, r_{\xi _2}\right)}^{r_{\xi _1}}\mathrm{d}r \, \right. \nonumber \\ && \left.+ \int _{\min \left(\max \left(p_{\min }, p_{\xi _1}\right),p_{\mathrm{max}}\right)}^{\min \left(p_{\xi _2},p_{\mathrm{max}}\right)} \mathrm{d}p \int _{\max \left(0, r_{\xi _2}\right)}^{R}\mathrm{d}r \right) \sum _{i=1}^{4} \frac{a_i^{\prime }}{\mu ^i} 2 F_{\rm w}(r, \mu )r. \nonumber\\ \end{eqnarray}
(60)
The weak-photon distribution (36) is defined by
\begin{equation} F_{\rm w}(\bar{\epsilon }, C_{\xi }) = \sum _{n,m} F_{\rm w}^{(n,m)} \bar{\epsilon }^nC_{\xi }^m. \end{equation}
(61)
Integrations over r in (58) and (60) yield expressions of the type
\begin{eqnarray} \begin{array}{l}\sum _{i=1}^{4}\frac{a_i}{\mu ^i}\int _{r_1}^{r_2} 2F_{\rm w}(r,\mu ) r\mathrm{d}r \!=\! \\ 2\epsilon _{\rm max}^2\sum _{i=}^{4} a_i\sum _{n,m}\sum _{l=0}^{m}\left({\begin{array}{cc}m \\ l \end{array}} \right) (-2)^{l} \mu ^{l-i + 1} \left[\bar{r}\bar{\epsilon }(\bar{r},\mu )^{n-l}\right]_{r_1}^{r_2}, \end{array}\nonumber\\ \end{eqnarray}
(62)
where
\begin{eqnarray} \left[\bar{r}\bar{\epsilon }\left(\bar{r},\mu \right)^{n-l}\right]_{r_1}^{r_2} \!=\! \left\lbrace \begin{array}{ll}2\mu \frac{\bar{\epsilon }\left(\bar{r}_2,\mu \right)^{n-l+1} \!-\! \bar{\epsilon }\left(\bar{r}_1,\mu \right)^{n-l+1}}{n-l+1} & \mathrm{if}\, n-l \ne -1\\ 2\mu \log \left( \frac{\bar{\epsilon }\left(\bar{r}_2,\mu \right)}{\bar{\epsilon }\left(\bar{r}_1,\mu \right)}\right) & \mathrm{if}\, n-l=1 \end{array}\right.. \end{eqnarray}
(63)

To obtain the final spectrum (58) (resp. 60), integration over μ (resp. over p) is possible analytically: the first line of (63) is a rational fraction that can be integrated through partial fraction decomposition and the second line yields expressions of the type ∫xklog (polynomial(x))dx (where k is integral) which values are given in most relevant textbooks such as (Gradshteyn et al. 2000). However, the resulting expressions may be lengthy and a direct numerical integration might sometimes be more efficient.

5 APPLICATIONS

5.1 Isotropic blackbody background distribution

Here, we propose to check our approximation equation (38) against the exact isotropic case described in Nikishov (1962), Agaronyan et al. (1983) and Boettcher & Schlickeiser (1997). We assume a high-energy photon hitting on a thermal soft photon background given by
\begin{equation} f_\mathrm{bb}(\epsilon ) = \frac{2}{(\hbar c {2\pi })^3}\frac{1}{e^{ \epsilon / k_{\rm B} T} - 1}, \end{equation}
(64)
where T is the temperature of the body and kB the Boltzmann constant. We choose a cut-off εmax = 20T (see equation 23) such that the neglected part of the blackbody spectrum (64) represents less than ∼e−20 ∼ 10−9, the total amount of background photons. We perform a Chebyshev interpolation (see e.g. Grandclément & Novak 2009) of ε2fbb(ε) on the 25 first Chebyshev polynomials achieving a relative accuracy better than 1000th everywhere and better than 10−6 for 0 ≤ ε/T ≤ 10, energies between which most of the photons are. This then allows us to derive the coefficients of the Laurent serie describing fbb with poles of order 1 and 2. Then, we produce the spectra of Figs 3 and 2.
Figure 2.

Spectra of outgoing leptons (electron or positron) for different strong-photon momenta k on a blackbody background at temperature Tbb (top panel). m is the mass of the electron, the speed of light and the Boltzmann constant are taken to be unity. The amplitudes are normalized to the amplitude of the peaks of each spectrum, and these amplitudes are reported on the lower panel. As in Fig. 2, these amplitudes are normalized to correspond to the cosmic-microwave background. The most intense peaks arise around a momentum k such that its reaction with a background photon at Tbb is at threshold, i.e. kTbb ∼ m2. The more above threshold, the more separated, narrow and low the peaks are. The separation of the peaks can be understood as a mere relativistic-frame effect, by analogy with a two-photon collision. [A color version of this figure is available online.]

Figure 3.

Comparison of absorption of high-energy photons on a blackbody background with Nikishov's formula (dashed line) and with our's (VMB, plain line). The scaling (66) is that of a blackbody at Tbb = 2.7 K (see formula 66) to give an estimate of the effect of the cosmic-microwave background. In this case, the energy of the strong photons ranges between k ∼ 100 TeV and k ∼ 108 TeV. The bottom panel shows the ratio between the two theories. The ratio of probabilities averaged over k is about 1.07. The peak of our curve occurs around 2.6m2/T while Nikishov's is around 1.9m2/T. The ratio between the two curves at the position of our peak is approximately 1.01. [A color version of this figure is available online.]

On the top panel of Fig. 2, we plot the total probability of absorbing a strong-photon of energy k as a function of
\begin{equation} \zeta = \frac{k k_{\rm b} T}{(mc^2)^2}. \end{equation}
(65)
This parametrization by ζ makes the temperature dependency simple
\begin{equation} W_{\vec{k} } \propto \left(\frac{k_{\rm B} T}{mc^2}\right)^3. \end{equation}
(66)
Here, we choose to take T = 2.7 K, which allows us to reproduce the result of Gould & Schréder (1966) (dashed line) concerning absorption on the cosmic microwave background. The lower panel of Fig. 2 shows the ratio between our formula and the exact formula of Nikishov (1962). It shows that our result is 50 percent off at ζ < 1 and asymptotically tends to the correct value for large ζ, the difference between the two curves is ∼ 10 per cent around the maximum of the curve located at ζ ∼ 2. On average on the range plotted on Fig. 2, our formula overestimates the reaction rate by 7 per cent.
A toy model can help us understand the shape of this curve. The peak of a blackbody spectrum is roughly at εbb ≃ 5kbT. The cross-section peaks when the CM energy is 1.4(2m), so if one approximates the blackbody spectrum to its peak one gets
\begin{equation} \epsilon _\mathrm{{\rm bb}} k(1-\cos \xi ) \simeq 3.9 m^2. \end{equation}
(67)
For an isotropic distribution of soft photons, collisions take place at every angle ξ ∈ [0, π]. Taking the intermediate value ξ = π/2, one obtains from (67) an estimate of ζ for the peak of reactions
\begin{equation} \zeta = \frac{k k_{\rm b} T}{m^2} \simeq 0.8 \end{equation}
(68)
which is the right order of magnitude. One could argue that at such energies reactions would occur more face-on, meaning ξ < π/2 which is consistent with the higher peak position found on Fig. 2. We now proceed to the computation of pair energy spectra.

Fig. 3 shows the pair-creation spectra for different values of ζ. Those spectra are directly computed using equation (38) and expressed as a function of p/k which allows the same scaling law as in equation (66), with p the momentum of one of the created leptons, and normalize each spectrum to unity such that the obtained spectral shape are universal, i.e. do not depend on the temperature of the blackbody or on the absolute value of k, but only on ζ. The shape and evolution of the spectra with the strong-photon energy is consistent with Agaronyan et al. (1983). In this paper, the authors consider spectra resulting from the reaction of two isotropic monoenergetic photon distributions with energies ε and k that are symmetrical with respect to (k + ε)/2. Here, every spectrum is symmetrical with respect to p/k = 0.5 as a result of neglecting ◯(ε/k) terms. Besides the shape of these spectra is very reminiscing of pair-creation in the photon-plus-magnetic-field process that is well known in the field of neutron-star magnetospheres (Daugherty & Harding 1983). The analogy is not fortuitous since the latter process can in principle be seen as the interaction of a strong photon with an assembly of magnetic-field photons. We see in Fig. 3 that each spectrum is made of two peaks that move apart and become narrower and weaker as the reaction occurs farther above threshold. Note that the narrowing is relative to the momentum span and not absolute.

The separation of the peaks at higher energies results from the fact that the cross-section favors alignment of ingoing and outgoing particles in the CM frame if the energy is much larger than the threshold energy. It follows that a Lorentz boost to the observer's frame along this axis results in a low-energy and a high-energy particle. The intensity of the peaks of course depends on the background distribution, but also on the cross-section which decays as log (τ)/τ (see Section 2.1). The latter dependency explains the above-threshold decrease of the peak intensity and the former explains the below-threshold decrease, as shown on the lower panel of Fig. 3. One notices that spectra are not smooth in their centre, which is naturally explained by our approximations that ensure continuity at the centre but not continuity of derivatives.

5.2 Above a hot neutron star

In this section, we consider a homogeneously hot neutron star at temperature T and two kinds of photons: the down photons and the up photons. Down photons are going radially towards the centre of the star while up photons are going in the opposite direction, away from the star. This configuration aims at approximating a pulsar magnetic pole. Indeed, in a pulsar magnetosphere high-energy photons are expected to be mostly created by curvature radiation of electrons and positrons flowing along magnetic field lines that can be considered radial at low altitudes above the poles. Note that we do not consider only a hot cap here but the full star, as means of geometrical simplification.

The case of pair production from photon–photon collisions in pulsar magnetospheres was studied by authors such as Zhang & Qiao (1998) and Harding et al. (2002). In these papers, the authors generalize the formula of Nikishov (1962) with a minimum energy threshold for the background distribution corrected by a factor (1 − cos θc)−1, where θc is the maximum viewing angle on the hot polar cap of the star. In other words, they consider an isotropic blackbody distribution where only photons within the viewing angle of the cap contribute, however with a threshold energy that corresponds to the largest incidence angle only since the threshold does not depend on the location of the emitter on the cap. Therefore, this approximation overestimates the threshold which generally translates in underestimating the reaction rate. This has little consequences when the viewing angle is wide, which is the case very close to the cap. However, one expects a faster decrease as one goes away from the cap and the factor (1 − cos θc)−1 grows larger. As an example, the authors of Zhang & Qiao (1998) compute a maximum reaction probability of 5.7 × 10−5 m−1 at a viewing angle of 90° when we get 6.7 × 10−5m−1 (see peak of the down-photon h = 10−3 curve on Fig. 5 for an estimate), but they obtain only 6.3 × 10−6 m−1 at 45° when we still have a probability of 4.3 × 10−5 m−1 (see their equation 9, for T = 106 K).

Besides, an interest of our formalism is that it can in principle deal with any other orientation of the strong photon with respect to the star, and in particular the up photons.

In this configuration, the distribution of soft photons is still given by equation (64) except that it is now zero when the angle ξ between the soft and the strong photon is beyond the horizon of the star as seen from the strong photon (see Fig. 4). For a photon going upward, the horizon is defined by
\begin{equation} \sin \xi < \frac{R_*}{R_* + h} = \sin \xi _{\mathrm{horizon}}, \end{equation}
(69)
where R* is the radius of the star (typically 10km) and h is the height above its surface. Consequently, we use equations (58) and (60) with angles ξ1 = 0, ξ2 = ξhorizon for an up photon and ξ1 = π − ξhorizon, ξ2 = π for a down photon.
Figure 4.

A neutron star of centre O and radius R* above which an up photon and a down photon are represented by radial arrows of opposite directions. Both photons are represented at a height h above the surface of the star. From this height, they can interact with soft photons coming from the surface of the star within a cone of aperture ξhorizon, equation (69), represented by dashed lines. The incidence angle between the strong photon and soft photons therefore lies between 0 and ξhorizon for the up photon (purple upward arrow), and between |$\pi - \xi _\mathrm{horizon}$| and π for the down photon (blue downward arrow). [A color version of this figure is available online.]

Fig. 5 shows the probability of reaction per unit length (we will sometimes say ‘reaction rate’) as a function of ζ at various heights h above the cap (left-hand panel), and as a function of h at various ζ (right-hand panel). As in the previous subsection, the temperature dependance is T3 for a given value of ζ. All the figures in this section are made with a fiducial temperature of 106 K. With this value, the conversion from ζ to k is: k ≃ 5.9 × 103ζmc2. At the lowest altitude we computed,  h = 10−3R*, the peak of the reaction rate is around ζ = 1.6 for down photons and about an order of magnitude higher for up photons ζ ≃ 16. This is a direct consequence of the threshold equation (67) given the less favourable incidence angles of up photons. Another point is that the position of the maximum shifts to lower ζ as height increases for down photons, but to larger ζs for up photons. As can be seen on the right-hand-side panel, the reaction probability per unit length is fairly stable (within a factor of two) until ∼1R*, after which it decays very sharply. The decay is sharper as ζ increases for down photons and smoother for up photons, which explains the crossing between some curves on the right-hand-side panel.

Figure 5.

Probability of reaction per metre of a strong photon of momentum k as a function of ζ = (kkBT)/mc2)2 and height h above a star of radius R*. Up-triangle markers represent photons going radially up from the star. Down-triangle markers represent photons going down to the star. The probability scales like T3, according to equation (66), and is here represented using a fiducial T = 106 K. Left-hand-side panel: probability as a function of ζ at various heights. Right-hand-side panel: probability as a function of height h at various ζ. [A color version of this figure is available online.]

A qualitative reasoning explains these behaviours. For a low-energy down photon (i.e. ζ ≲ 1), most of the soft photons most likely to react are in a narrow cone almost face-on with the strong photon. The aperture of the cone defines the limit beyond which the reaction is below threshold. When the strong photon is higher, the almost face-on soft photons are the last to disappear because of the shrinking of the viewing angle. As energy rises, this cone becomes wider since soft photons provoking a near-threshold reaction are located at a wider angle according to formula (67). Inside the first cone also appears a co-axial cone with a narrower aperture inside which photons are not contributing significantly anymore, since reactions are too far above threshold (and therefore the cross-section is too small) because of small incidence angles. At large strong-photon energies (ζ ≫ 1), the soft photons close to the outer cone are the first to disappear when the viewing angle shrinks because of a larger height. This explains the faster decay of the reaction probability with h for larger ζs of down-photon curves on Fig. 5. The same kind of reasoning applies for up photons. Because soft photons are arriving ‘from behind’, there is always an inner cone inside which the reactions are below threshold, and an outer cone limited by the angle beyond which the cross-section is too small if ζ ≫ 1 or the viewing angle of the strong-photon energy is small enough. The lower the energy of the strong photon the wider the outer cone and the most sensitive to viewing angle the reaction rate is. That explains why, contrary to down photons, the reaction rate decays slower with altitude when ζ is larger in Fig. 5. With this reasoning, one also understands why the energy of the reaction-rate peak (left-hand panel) is quite stable at low altitudes and becomes smaller for down photons at high altitudes (≳1R*) or larger for up photons.

Fig. 6 shows the optical depth of strong photons as a function of ζ through 10R* from the surface. The optical depth is defined by
\begin{equation} \tau _\zeta (10 R_*) = \int _0^{10 R_*} W_\zeta (h) \mathrm{d}h. \end{equation}
(70)
Because of the effects mentioned above the peak for down photons is slightly shifted downward at ζ ≃ 1.4 while upward for up photons at ζ ≃ 30. The corresponding typical Lorentz factors of the created particles are 8 × 103 and 2 × 105, respectively. The peak optical depths are, respectively, τ ≃ 8.8 and τ ≃ 0.23 at a temperature of 106 Kelvins. One concludes that at this temperature more than three out of four up photons at the peak energy escape the magnetosphere if no other reaction or source of soft photons opacifies it. The magnetosphere may become opaque if the star is hotter than ∼1.6 × 106 K, temperature for which the maximum optical depth reaches 1 owing to the T3 dependence of the reaction rate. Down photons with ζ between ∼0.25 and ∼64 have optical depths larger than unity and therefore are absorbed before they hit the star except if they are emitted at very low altitudes hR*. The maximum optical depth of down photons is below 1, namely the magnetosphere is transparent, for a temperature below 0.5 × 106 K. Let's notice that our approximation of a uniformly hot star obviously leads to overestimating the optical depth on distances larger than the size of an actual hotspot.
Figure 6.

Reaction rate in equation (28) integrated from 0 to 10R* for up and down photons as a function ζ. The amplitudes are valid for a homogeneously hot star of temperature T = 106 K, and can be converted to other temperatures using the T3 scaling law (66). For photons going down towards the star, the peak is at ζ ≃ 1.4 with an amplitude of ≃8.8. For photons going up away from the star, the peak is at ζ ≃ 30 with an amplitude of ≃0.23. [A color version of this figure is available online.]

Fig. 7 shows the energy spectra of the created leptons (left-hand panel) and the evolution of the position and widths of the peaks as a function of ζ at various heights (right-hand panel). The spectra have the same double-peaked structure as in the isotropic case (Fig. 3) but evolve differently depending on the orientation of the strong photon. The general principle is the same: the more above threshold the more separated peaks, with the consequence that they narrow when they get close to the limits of p/k ∈ [0, 1]. For down photons, the width of the peaks wp (in unit of k) has very little dependence on altitudes which is due to the fact that for the range of ζ ≲ 20 visible on this plot (right-hand panel), the efficient soft photons are mostly face-on and suffer no effect of viewing angle. The same thing applies for the position of the most energetic peak pp (and the least energetic at k − pp). Down-photon peaks are wide wp ∼ 0.45 for ζ ≲ 2 and then sharply narrow while their position smoothly goes from pp/k ∼ 0.8 to pp/k ≲ 1 at large ζs. On the contrary, up photons are very sensitive to altitude, which is explained by the fact that the higher above the star, the narrower the viewing angle and therefore the incidence angle, and the more energetic up photons need to be for the reaction to be at or above threshold. As a consequence up-photons peaks are very centred at low values of ζs, with pp/k ∼ 0.6, and are even more centred at higher altitudes. With ζ rising, the energy distribution becomes increasingly asymmetric as pp/k → 1, although it takes a larger ζ at higher altitude. Similarly, peak widths are growing with ζ until a maximum wp ∼ 0.45 at a ζ all the more large that altitude is high, after which wp drops sharply. This sharp change of slope happens because the two peaks separate (see comment of Fig. 7).

Figure 7.

Left-hand-side panel: example of two normalized spectra of energy of created particles. These spectra are normalized to the amplitude of the largest peak, and the energy in abscissa is normalized to the energy of the incident strong photon k. In both cases, k corresponds to ζ = 10 at an altitude h = 0.5R*, and the only difference resides in the up or down orientation of the strong photon. As in the isotropic case 3, spectra are generally made of two peaks more or less thin and separated. The width at half-maximum of peaks wp is defined in the two possible cases: if one side of a peak never reaches its half before rising again to another peak in which case the width is taken to be half of the double-peak width, or if the peak is well defined on both sides in which case the definition is straightforward. The position of the most energetic peak pp/k is defined as well. Right-hand-side panel: evolution with ζ of positions pp/k of the higher energy peak (curves on the higher part of the plot), and widths at half-maximum wp (curves on the lower part of the plot) for up and down photons at various heights h (in units of R*). Positions are ranging from 0.5 at low ζs which corresponds to a perfectly centred peak or to a null spectrum when a reaction is below threshold (lowest energies of up photons), to ≃0.98 at large ζs. The horizontal dotted line shows the positions at which the ratios between the two peaks is 10. Widths at half-maximum are rising to ∼0.45 until the two peaks separate and drop sharply to ≃0.029. [A color version of this figure is available online.]

Fig. 8 shows the normalized angular spectra for both up and down photons, and both higher energy (p > k/2) and lower energy (p < k/2) outgoing leptons at various values of ζ. It is remarkable that apart from their amplitudes (not visible on this normalized plot), these spectra do not change much with height apart at large and very unlikely angles, and therefore we limit ourselves to only one height. These spectra are monotonously decreasing as the angle becomes larger, and the larger ζ the larger the outgoing angle. Lower energy leptons have larger outgoing angles than their higher energy counterpart and are not created below a minimum angle defined in equation (51). For a given ζ, leptons created from down photons are always going out at larger angles, and the difference is growing at larger angles of the spectrum. In a pulsar magnetosphere, the outgoing may be important because the pairs will radiate more or less synchrotron radiation depending on their momentum perpendicular to the local magnetic field. We see here that the angles with respect to the progenitor strong photon are overall very small, which is expected from relativistic collimation. If one assumes that strong photons are produced though curvature radiation along the magnetic-field lines, then the angle distributions presented on Fig. 8 matter only if the mean-free path is much shorter than the radius of curvature of the field line. This is not the case with the parameters presented in this section, and would probably require an extra source of photons.

Figure 8.

Angular spectra for various values of |$\zeta = \frac{k k_{\rm b} T}{(mc^2)^2}$| and both up and down orientations (respectively, up-triangle and down-triangle markers) normalized to their maximum values. Solid lines correspond to angular spectra of the higher energy leptons (p > k/2) and dashed-line spectra to their lower-energy counterparts (p < k/2). The legend shows only solid lines as the corresponding dashed line always directly follows at larger angle. Spectra are ordered from smaller to larger angles by decreasing value of ζ. All these spectra correspond to a height above the star h = 0.001R* which is representative for all other heights. Indeed, their amplitudes significantly change with height but their shapes (and therefore the shown normalized spectra) barely change. [A color version of this figure is available online.]

6 DISCUSSION

Recent simulations of aligned millisecond-pulsar magnetospheres indicate that significant pair production may occur near the so-called separatrix gap and y point (see Cerutti & Beloborodov 2016 and references therein) near the light cylinder. This implies that the source of pairs be photon–photon collisions. However, in the most detailed modelling of pair creation realized by Chen & Beloborodov (2014), photon–photon pairs are created with a constant and uniform mean-free path of 2R*. If one assumes that the source of soft photons is only provided by the star, this assumption seems reasonable close to the star, h < 2R*, but greatly underestimated beyond owing to the exponential cut-off of the reaction rate with altitude (Fig. 5). This issue can be overcome if another source of soft photons can be found, resulting for example from synchrotron radiation near the light cylinder. Moreover, in these simulations, the direction of strong photons relative to the soft-photon sources is not taken into account, which can have an effect of several orders of magnitude on reaction rates with a strong dependence on strong-photon energies (see Fig. 6). The energy separation of the two outgoing leptons (Fig. 7) may also have an important impact on the subsequent synchrotron radiated by the pair. Indeed, the synchrotron peak frequency scales like γ2, where γ is the Lorentz factor of the particle around the magnetic field. Therefore, a typical situation in which the higher energy lepton takes 10 times more energy than the other (dotted line on Fig. 5) results in two synchrotron peaks radiated two orders of magnitude apart. This situation is reached at values of ζ for which the optical depth on Fig. 6 is still high, i.e. more than half the peak value. Note that we implicitly assume here that both particles share the same angle with respect to the local magnetic field, which is justified by small outgoing angles shown in Fig. 8.

The photon–photon mechanism competes with the photon-magnetic-field mechanism |$\gamma + \vec{B} \rightarrow e_+ e_-$| for the creation of pairs in pulsar magnetospheres. For our present discussion, we focus on magnetic fields smaller than the critical magnetic field Bc = 4.4 × 109 Teslas. Additionally the photon energies produced by curvature or synchrocurvature radiations in pulsar magnetospheres cannot exceed ∼100 GeV owing to radiation reaction (see e.g. Viganò et al. 2014). With these two limits, the photon-magnetic-field reaction rate can be computed with the asymptotic expression (Tsai & Erber 1975; Daugherty & Harding 1983)
\begin{eqnarray} W_{\gamma \vec{B}} \!\mathop{\simeq }\limits_ {\chi \ll 1} 4.3\!\times\! 10^{9} \frac{B\sin \theta }{B_c} \exp \left(\!-\!\frac{4}{3\chi }\!\right) \mathrm{m}^{-1} \, \mathrm{with} \, \chi \!=\! \frac{\hbar \omega }{2mc^2}\frac{B\sin \theta }{B_c},\nonumber\\ \end{eqnarray}
(71)
where ℏω is the energy of the gamma photon, B the intensity of the local magnetic field and θ the angle between the direction of the magnetic field and the direction of the photon. The photon-magnetic-field optical depth heavily depends on the magnetic-field geometry: typically, a photon in the pulsar magnetosphere is emitted parallel to the local magnetic field due to relativistic beaming, thus starting with a reaction rate |$W_{\gamma \vec{B}} = 0$| which increases along the propagation as θ increases. The upper limit of the reaction rate of the photon can be estimated by considering sin θ = 1 everywhere, although this is bound to largely overestimate the probability of creating a pair.

Fig. 9 shows a comparison between the reaction rate of equation (71) for a range of values of Bsin θ versus the photon–photon reaction rates computed at various altitudes of Fig. 5. Considering that the range of probable photon energies lies below 100 GeV, one sees that the photon–photon mechanism can dominate near the surface of millisecond pulsars where Bsin θ ≲ 105 Teslas, and clearly dominates for down photons at altitudes h ≥ 10R* where, assuming a dipolar magnetic field where B ∝ (R* + h)−3, one has Bsin θ ≲ 100 Teslas for millisecond pulsars. For up photons, it is less clear which mechanism dominates without taking into account a particular magnetic geometry. It should also be noted that, if a comparable soft-photon density can be achieved in the outer magnetosphere as is necessary in some recent simulations (e.g. Chen & Beloborodov 2014), the photon–photon mechanism would be largely dominating the photon-magnetic-field mechanism in this region of the magnetosphere.

Figure 9.

Comparison of reaction rates of the γγ → e+e process (solid lines) versus the |$\gamma + \vec{B} \rightarrow e_+e_-$| process (dashed lines). The photon–photon reaction rates are identical to those on the left-hand-side panel of Fig. 5, but a temperature of T = 106 K is taken giving the photon energies in GeV reported on the upper horizontal axis corresponding to the ζ values on the lower axis. The photon-magnetic-field reaction rates are given in the χ ≪ 1 approximation (see the text) within which they depend only on the the photon energy and the magnetic intensity perpendicular to the photon direction Bsin θ (in Teslas), where θ is the angle between the local magnetic field and the photon direction. [A color version of this figure is available online.]

We have neglected two effects possibly important in our application to a hot neutron star (Section 5.2): general relativity and the effect of the strong magnetic field on the pair-production cross-section. The former effect, general relativity, redshifts the spectrum of soft photons, and enlarges the visible horizon of the star. Indeed, light bending due to the gravitational field of the star curves the trajectories of soft photons in such a way that part of the surface beyond the geometrical horizon of the star becomes ‘visible’ by the strong photon. Simple analytical formulas for the effective soft-photon distribution with general-relativistic effects have been given by Beloborodov (2002) and Turolla & Nobili (2013). However, these expression are valid for an observer (in our case the strong photon) located at infinity and are therefore inapplicable in the present case where we consider reactions between 10−3R* and 100R*. Instead, one would need to compute numerically the geodesics followed by the soft photons between the surface and the strong photon.

The second neglected effect is the effect of a strong magnetic field on the cross-section for photon–photon pair creation. It has been worked out by Kozlenkov & Mitrofanov (1986). This cross-section shows a sawtooth behaviour at energies corresponding to the quantified Landau levels of the outgoing leptons. Unfortunately, this cross-section is very unwieldy for practical calculations. However, its effect is low or moderate in magnetic fields much lower than the critical field Bc, which fortunately corresponds to the range of parameters where photon–photon pair creation can dominate over photon-magnetic-field pair production (see above and Fig. 9), and therefore the range of interest of our formalism. As mentioned before, the other important domain of application of our formalism is the outer magnetosphere where the magnetic field is also much smaller than Bc, including in the case of young pulsars such as the Crab.

7 CONCLUSION

We propose a formalism to analyse photon–photon pair creations with an arbitrarily anisotropic soft-photon background. This formalism allows us to calculate energy and angle spectra of outgoing pairs, as given by formulas (58) and (60), respectively.

Calculations are carried using two approximations: the first being that the strong photon is much more energetic than the soft-photon cut-off energy εmax, and the second that the outgoing higher energy lepton of momentum |$\vec{p}$| be very aligned with the progenitor strong photon of momentum |$\vec{k}$| in the sense that |$(\vec{k} - \vec{p})_{\perp } / (\vec{k} - \vec{p})_{\parallel } \ll 1$|⁠, (A22), where perpendicular and parallel components are taken with respect to |$\vec{k}$|⁠. This latter approximation is the most stringent one. Indeed, one can show that the inequality itself (<1) is always true within the frame of our first approximation, but its large validity (≪ 1) comes if the reaction is far above threshold, i.e. kεmax/m2 ≫ 1.

In Section 5.1, we compare our formalism with the exact formula that can be found in the literature (Nikishov 1962, or Agaronyan et al. 1983 equations 4 and 5 for a more detailed formulation), and show that our approximated formulation gives results accurate at ∼ 7 per cent on average, with ∼ 10 per cent near the peak and asymptotically tend to the exact value at high energies. However, the difference can be as large as ∼ 50 per cent at low energies. We show pair spectra that are consistent with those of Agaronyan et al. (1983) in the isotropic case.

In Section 5.2, we show that the differences created by the strong anisotropy of radiation near a hot neutron star are much more important than a few percent, potentially reaching several orders of magnitude depending on energy, direction of the strong photon, and altitude above the star. We consider two directions for strong photons: radially towards the star (down photons) and away from the star (up photons). In both cases, reaction rates are stable until 1R*, before undergoing an exponential cut-off. However, the peak of strong-photon absorption occurs at an energy ∼10 times larger for up photons. Energy pair spectra show two peaks symmetric with respect to k/2, similarly to the isotropic case. These peaks separate as the energy of the reaction rises. We show that such a difference in energy between the two outgoing leptons can importantly affect the synchrotron emission of the pairs for a large range of strong-photon energy compared to a simple model in which both components of a pair take away the same energy.

These findings are meant to contribute to a better modelling of pair creation from photon–photon collisions in pulsar magnetospheres. Recent millisecond-pulsar-magnetosphere simulations gave an important role to this pair-production mechanism (Cerutti & Beloborodov 2016). However, the current state of modelling leaves an important uncertainty on the amount of soft photons needed to sustain such pair discharges. The results of this work provide means to estimate the mean-free path on a soft-photon background resulting from a homogeneously hot neutron star. Moreover, it provides formulas to obtain results with virtually any soft photon distribution, in particular resulting from secondary synchrotron close to the light cylinder. The possibility to generate energy spectra allows us to differentiate between the two components of a pair and therefore to differentiate their synchrotron emissions.

1

In international units |$r_{\rm e} = \frac{e^2}{4\pi \epsilon _0 m c^2} \simeq 2.8179 \times 10^{-15}$| metres, with ε0 the electric permittivity of vacuum.

2

Include a factor c in the definition of j when it is not assumed that c = 1.

3

The energy of one photon in the CM of two photons, one at energy k and another at energy ε, is |$\sqrt{k\epsilon (1 - cos\omega )}$|⁠, where ω is the angle between the two photons three-momenta

4

Then one shows easily from (A30) that x2 + y2 has a minimum for |$z_m = -z_0\frac{L/z_0}{1 + L^2/z_0^2} = -z_0 L/z_0 + \bigcirc \left(L^2/z_0^2\right)$| (⁠|$\bigcirc \left(L^2/z_0^2\right) = \bigcirc \left(1\right)$|⁠). Hence, |$z_m \simeq - \frac{A}{\alpha }\frac{\beta _\Vert }{\alpha }$| while |$z_\delta = \frac{A}{\beta _\Vert }$|⁠. Using the fact that |$\alpha \lesssim \beta _\Vert$| one gets that zmzδ, which means zm is slightly under the bottom of the hyperboloid. Since x2 + z2 can be easily shown to be a growing function of z, its smallest value can only be zδ.

REFERENCES

Agaronyan
F. A.
,
Atoyan
A. M.
,
Nagapetyan
A. M.
,
1983
,
Astrophysics
,
19
,
187

Beloborodov
A. M.
,
2002
,
Astrophys. J. Lett.
,
566
,
L85

Berestetskii
V. B.
Lifshitz
E. M.
Pitaevski
L. P.
,
1982
,
Quantum Electrodynamics
.
Oxford Univ. Press
,
Oxford, Royaume-Uni

Boettcher
M.
,
Schlickeiser
R.
,
1997
,
A&A
,
325
,
866

Bonometto
S.
,
Rees
M. J.
,
1971
,
MNRAS
,
152
,
21

Burns
M. L.
,
Harding
A. K.
,
1984
,
AJ
,
285
,
747

Cerutti
B.
,
Beloborodov
A. M.
,
2016
,
Space Sci. Rev.
, pp
1
26

Chen
A. Y.
,
Beloborodov
A. M.
,
2014
,
Astrophys. J. Lett.
,
795
,
L22

Daugherty
J. K.
,
Harding
A. K.
,
1983
,
AJ
,
273
,
761

Goldreich
P.
,
Julian
W. H.
,
1969
,
AJ
,
157
,
869

Gould
R. J.
,
Schréder
G.
,
1966
,
Phys. Rev. Lett.
,
16
,
252

Gradshteyn
I. S.
Ryzhik
I. M.
Jeffrey
A.
Zwillinger
D.
,
2000
,
Table of Integrals, Series, and Products
.
Academic Press
,
New York

Grandclément
P.
,
Novak
J.
,
2009
,
Living Rev. Relativ.
,
12
,
1

Harding
A. K.
,
Muslimov
A. G.
,
Zhang
B.
,
2002
,
AJ
,
576
,
366

Kozlenkov
A. A.
,
Mitrofanov
I. G.
,
1986
,
Zh. Eksp. Teor. Fiz.
,
91
,
1978

Nikishov
A. I.
,
1962
,
Sov. Phys. JETP
,
14

Ruffini
R.
,
Vereshchagin
G.
,
Xue
S.-S.
,
2010
,
Phys. Rep.
,
487
,
1

Sturrock
P. A.
,
1971
,
AJ
,
164
,
529

Svensson
R.
,
1987
,
MNRAS
,
227
,
403

Timokhin
A. N.
,
Harding
A. K.
,
2015
,
AJ
,
810
,
144

Tsai
W.-y.
,
Erber
T.
,
1975
,
Phys. Rev. D
,
12
,
1132

Turolla
R.
,
Nobili
L.
,
2013
,
AJ
,
768
,
147

Vassiliev
V. V.
,
2000
,
Astropart. Phys.
,
12
,
217

Viganò
D.
,
Torres
D. F.
,
Hirotani
K.
,
Pessah
M. E.
,
2014
,
MNRAS
,
447
,
1164

Zhang
B.
,
Qiao
G. J.
,
1998
,
A&A
,
338
,
62

APPENDIX A: Derivation of the general result

We show that domain of integration can be approximated by a hyperboloid of revolution. Then, we compute the integral Wk over this surface assuming that the distribution function of weak photons is given by a polynomial (e.g. Taylor expansion). A variety of notations and relations are used and we summarize them in Appendix B.

A1 Parametrization of L by the three momentum of the weak photons

The spectrum of pair creation is the density of probability of making a pair as a function of the energy of one of the particles. By definition, it is symmetric with respect to half of the total energy k + ε ≃ k: if one of the particles has an energy p then the other has k − p as a result of energy conservation. Therefore, we consider only the upper half of the spectrum, for p > k/2 and
\begin{equation} \frac{\mathrm{d}W}{\mathrm{d}p}(p) = \frac{\mathrm{d}W}{\mathrm{d}p}(k-p) \end{equation}
(A1)
Therefore, we are left with the very helpful ordering
\begin{equation} k \gtrsim p \gg m, \epsilon _{\rm max}, \end{equation}
(A2)
which allows us to write:
\begin{equation} P^0 = p + \frac{m^2}{2p} + \bigcirc \left(\frac{m^2}{p^2}\right). \end{equation}
(A3)
Further, we learn from the angle-averaged cross-section (Berestetskii et al. 1982) that when the reaction is way above threshold, one of the particles of the pair takes most of the energy while the other takes almost nothing (Section 2), which reinforces our assumption. In the following calculation, we note ◯(n) a development up to a bounded function of
\begin{equation} \left(\frac{m}{k}\right)^n \sim \left(\frac{m}{p}\right)^n \sim \left(\frac{\epsilon _{\rm max}}{k}\right)^n \sim \left(\frac{\epsilon _{\rm max}}{p}\right)^n. \end{equation}
(A4)
This leads to the conclusion that |$\vec{p}$| is almost aligned with |$\vec{k}$|⁠. Indeed, we can show that εmin < εmax implies that momentum can be conserved only if cos θ > Cmin, where
\begin{eqnarray} C_{\rm min} &=& \frac{k+\epsilon _{\rm max}}{kp}(P^0 - \epsilon _{\rm max}) \nonumber \\ &\quad-& \frac{\epsilon _{\rm max}}{kp} \sqrt{k^2 + p^2 - 2kP^0 + 2k\epsilon _{\rm max}- 2P^0 \epsilon _{\rm max}+ \epsilon _{\rm max}^2}. \nonumber\\ \end{eqnarray}
(A5)
This is approximated as
\begin{equation} C_{\rm min} =1 + \frac{m^2}{2p^2} - 2 \frac{\epsilon _{\rm max}(k - p)}{kp} + \bigcirc \left(3\right). \end{equation}
(A6)
Therefore, we set |$\vec{k}$| as the main axis of our coordinate system (Fig. A1), parallel to the unit vector |$\vec{z}$| of the direct triad |$(\vec{x}, \vec{y}, \vec{z})$|⁠. For a weak photon of energy ε:
\begin{equation} 1 - \cos \theta \le 2 \frac{\epsilon (k - p)}{kp} - \frac{m^2}{2p^2} + \bigcirc \left(3\right) \end{equation}
(A7)
By squaring relevantly the mass-shell constrain (16 b), one obtains the following quadratric constrain:
\begin{equation} \left(\epsilon \left(P^0 - k\right)\right)^2 = ( \vec{x} \cdot (\vec{k} - \vec{p}) + K\cdot P )^2 \end{equation}
(A8)
which can be rewritten as
\begin{equation} \vec{x}(\alpha ^2 1 - \overline{\beta }) \vec{x} - 2A\vec{\beta }\cdot \vec{x} = A^2, \end{equation}
(A9)
where
\begin{eqnarray} A & = & K\cdot ,P \end{eqnarray}
(A10)
\begin{eqnarray} \vec{\beta } & = & \vec{k} - \vec{p}, \end{eqnarray}
(A11)
\begin{eqnarray} \alpha & = & k - P^0, \end{eqnarray}
(A12)
and |$\overline{\beta }$| is defined by
\begin{equation} \vec{x}\overline{\beta }\vec{x} = (\vec{\beta }\cdot \vec{x})^2. \end{equation}
(A13)
We find
\begin{eqnarray} \overline{\beta } & = & \left(\begin{array}{ccc}\beta _x^2 & 0 & 0 \\ 2\beta _x\beta _y & \beta _y^2 & 0 \\ 2\beta _x\beta _z & 2\beta _y\beta _z & \beta _z^2 \end{array}\right). \end{eqnarray}
(A14)
Let's rewrite equation (A9) in a dimensionless form,
\begin{equation} \vec{x}\left(\frac{\alpha ^2}{A^2} 1 - \frac{\overline{\beta }}{A^2}\right) \vec{x} - 2\frac{\vec{\beta }}{A}\cdot \vec{x} = 1. \end{equation}
(A15)
The three proper values of this quadratic form are
\begin{equation} \left(\frac{\alpha ^2}{A^2} - \frac{\beta _x^2}{A^2}, \frac{\alpha ^2}{A^2} - \frac{\beta _y^2}{A^2}, \frac{\alpha ^2}{A^2} - \frac{\beta _z^2}{A^2}\right). \end{equation}
(A16)
The geometrical type of this quadratic form is determined by the signs of its proper values. For this, we express the different quantities using the approximation defined in equation (A2),
\begin{eqnarray} A & = & m^2\left(\frac{1}{2}\frac{k}{p} + \frac{kp}{m^2}(1 - \cos \theta ) \right)+ \bigcirc \left(1\right) \end{eqnarray}
(A17)
\begin{eqnarray} & = & 2m^2 \frac{k-p}{m}\frac{\epsilon _{\rm max}}{m}\mu \nonumber \\ \alpha & = & m\left( \frac{k}{m}-\frac{p}{m} - \frac{m}{2p} \right) + \bigcirc \left(2\right) \end{eqnarray}
(A18)
\begin{eqnarray} \sqrt{\beta _x^2 + \beta _y^2} = \beta _\perp & = & p\sqrt{2(1 - \cos \theta )} + \bigcirc \left(3\right) \end{eqnarray}
(A19)
\begin{eqnarray} \beta _z = \beta _\Vert & = & p\left(\frac{k}{p} - \cos \theta \right) \end{eqnarray}
(A20)
It can be shown that, provided that |$\epsilon _{\rm max}< \frac{3}{8}m$| and for any relevant θ or p,
\begin{equation} \alpha > \beta _\perp . \end{equation}
(A21)
Similarly, provided that |$\epsilon _{\rm max}< \frac{1}{4}m$| (see equation 23),
\begin{equation} \frac{\beta _\perp }{\beta _\Vert } < 1 \end{equation}
(A22)
The smaller εmax with respect to m the more effective these constraints will be. (Note that the functions are monotonous on the appropriate range.) Moreover, the maximum value of 1 − cos θ is the limiting factor for εmax, and therefore these limits are less stringent if one considers creation of particles at smaller angles. Besides, |$\sqrt{k\epsilon _{\rm max}}$|3 is the higher bound of the energy of the two photons in the CM frame, and given our condition km, εmax close to m leads to an energy way above threshold in equation (10), and therefore very unlikely to happen (Section 2), although it depends on the angle of incidence of the weak photon on the strong one as well. For these reasons, we should consider that the higher limit for εmax is a ‘smooth’ one meaning that most photons of the weak distribution should actually not be close to εmax, even when εmax is close to the limit m/4, except if one has a very peculiar photon distribution. This discussion a posteriori justifies condition 23. The proper vectors associated with the proper values in equation (A16) are
\begin{eqnarray*} \vec{v}_1&=&(0,0,1), \nonumber \\ \vec{v}_2&=&(0,1,v_{2z}), \nonumber \\ \vec{v}_3&=&(1,\frac{2 \beta _x \beta _y}{\beta _y^2-\beta _x^2)},v_{3z}), \nonumber \end{eqnarray*}
with
\begin{eqnarray} v_{2z}&=&\frac{2 \beta _y \beta _z}{\beta _z^2-\beta _y^2} << 1, \end{eqnarray}
(A23)
\begin{eqnarray} v_{3z}&=&\left[ 2 \beta _x \beta _z -4 \frac{\beta _x \beta _z \beta _y}{\beta _z^2-\beta _y^2}\right] \frac{1}{\beta _z^2-\beta _x^2} <<1. \end{eqnarray}
(A24)
The above components are negligible in virtue of equations (A21) and (A22). Therefore, any vector parallel to the z-axis has its image parallel to the z-axis, and any vector perpendicular to the z-axis has its image roughly perpendicular to the z-axis.
Figure A1.

Coordinate system. In our approximation, |$\vec{k}$| is a quasi-symmetry axis.

We can simplify the orthogonal proper values in equation (A16), which are now both equal to:
\begin{equation} \frac{\alpha ^2}{A^2}. \end{equation}
(A25)
It can be shown that
\begin{eqnarray} \alpha & < & \beta _\Vert , \end{eqnarray}
(A26)
\begin{eqnarray} \frac{\beta _\Vert - \alpha }{m} & =& \bigcirc \left(1\right). \end{eqnarray}
(A27)
This implies that the parallel proper value (the third one in equation A16) is negative. Because the parallel proper value is negative while the two orthogonal values are positive, the quadratic form in equation (A9) describes a paraboloid of revolution.
The above remarks and equations (A21) and (A22) allow us to simplify the quadratic form in equation (A9). We are left with
\begin{eqnarray} \!-\!\left(\frac{\beta _z^2}{A^2} \!-\! \frac{\alpha ^2}{A^2} \right)(z + z_0)^2 \!+\! \frac{\alpha ^2}{A^2}(y^2 + z^2) \!=\! 1 \!-\! z_0^2\left(\frac{\beta _z^2}{A^2} \!-\! \frac{\alpha ^2}{A^2} \right),\nonumber\\ \end{eqnarray}
(A28)
where
\begin{equation} z_0 = \frac{A\beta _\Vert }{\beta _\Vert ^2 - \alpha ^2}. \end{equation}
(A29)
Dividing everything by |$\left(\frac{\beta _z^2}{A^2} - \frac{\alpha ^2}{A^2} \right)$| we get our final, although not fully standard, form of equation (A9):
\begin{equation} \frac{(z + z_0)^2}{z_0^2} - \frac{x^2 + y^2 }{L^2} = 1 - \delta ^2, \end{equation}
(A30)
where the characteristic orthogonal radius L and the displacement δ are
\begin{eqnarray} L^2 & = & \frac{A^2}{\alpha ^2}\frac{\beta _\Vert ^2}{\beta _\Vert ^2 - \alpha ^2}, \end{eqnarray}
(A31)
\begin{eqnarray} \delta ^2 & = & \frac{\beta _\Vert ^2 - \alpha ^2}{\beta _\Vert ^2}. \end{eqnarray}
(A32)
Equation (A28) describes a paraboloid of revolution of axis |$\vec{z}$|⁠, i.e. parallel to the strong photon three-momentum |$\vec{k}$|⁠. We notice that |$L^2/z_0^2 = \bigcirc \left(1\right)$| meaning that the hyperboloid is very steep around the parallel axis. Besides, the small displacement δ, δ2 = ◯(1), is responsible for a shift of the bottom of the paraboloid under the plane of zero parallel momentum. This corresponds to reactions with head-on weak photons that are in general of smaller energies, as shown on Fig. A2.
Figure A2.

Representation of the mass shell (16b) within approximation (23).

We must remember that the inclusion of L into the paraboloid defined in equation (A30) is derived from the condition in equation (16b). It must be completed with the condition in equation (16e) that reads
\begin{equation} \vec{x} \cdot \vec{\beta } + A \ge 0. \end{equation}
(A33)
From a geometrical point of view, this means that the relevant photons are those with momenta above the plane of normal vector |$\vec{\beta }$| of equation |$\vec{x} \cdot \vec{\beta } =- A$|⁠. Using relation (A22), this approximates to the plane orthogonal to the parallel direction, |$\vec{z}$|⁠, at position:
\begin{equation} z_{\rm min}= - \frac{A}{\beta _\Vert } \simeq -2 \epsilon _{\rm min}. \end{equation}
(A34)
As a consequence, only the upper sheet of the hyperboloid defined in equation (A30) corresponds to the physical mass shell. Indeed, this hyperboloid crosses the parallel axis |$\vec{z}$| at abscissa |$z_\delta ^\pm$| such that:
\begin{equation} z_\delta ^\pm = -z_0\left(1 \pm \sqrt{1 - \delta ^2}\right). \end{equation}
(A35)
This relation takes into account the fact that z0 ≫ |zmin|. Because δ< < 1,
\begin{equation} z_\delta = -z_0(1 - \sqrt{1 - \delta ^2}) \simeq -\frac{1}{2}z_0\delta ^2. \end{equation}
(A36)
Further, one may show that |$z_\delta = \min {\left\Vert \vec{x} \right\Vert }$| 4 which corresponds to the physical idea that the smallest weak photon that can produce a pair is the one that hits the strong photon head-on. How does it compare to εmin determined in the previous section? With the notations used in this section,
\begin{equation} \epsilon _{\rm min}= \frac{A}{\left\Vert \vec{\beta } \right\Vert + \alpha }. \end{equation}
(A37)
Using the approximation in equation (A22), |$\left\Vert \vec{\beta } \right\Vert = \beta _\Vert + \bigcirc \left(\beta _\perp / \beta _\Vert \right)$|⁠, we get
\begin{equation} \epsilon _{\rm min}= z_\delta + \bigcirc \left(\beta _\perp / \beta _\Vert \right). \end{equation}
(A38)
This is consistent with the definition of εmin in equation (20) as the minimum energy allowed in L.

A2 Integration

We now compute the integral Wk defined in equation (14). Within the frame of our approximations, equation (A30) shows that the probability of making a pair is symmetric around |$\vec{k}$|⁠. This leads to a first angular integration of ϕ that yields a 2π factor. The differential element |$\mathrm{d}^3\vec{k_w} =\mathrm{d}x\mathrm{d}y\mathrm{d}z$| is constrained by equation (A30) that defines L. Thus, we write z = z(x, y, p) through the constrain L(p, cos θ, k) and
\begin{equation} \mathrm{d}z = \left|\frac{\mathrm{\partial} z}{\mathrm{\partial} p} \right| \mathrm{d}p. \end{equation}
(A39)
We can now write (28) as follows
\begin{eqnarray} W_{\vec{k}} = c &{2\pi }& \int _{C_1}^{C_2} \mathrm{d}\cos \theta \int _{p_1}^{p_2} \mathrm{d}p \nonumber\\ & &\times\int _{L_-(x,y)} \frac{\mathrm{d}^2\sigma }{\mathrm{d}\Omega } \frac{K_{\rm s}\cdot K_{\rm w}}{K_{\rm s}^0 K_{\rm w}^0} f_{\rm w}(\vec{k_{\rm w}})\left|\frac{\mathrm{\partial} z}{\mathrm{\partial} p} \right| \mathrm{d}x\mathrm{d}y, \end{eqnarray}
(A40)
where L(x, y) is the projection of L on to the (x, y) plane. We need to expand the different quantities appearing in (A40). Let us start with an explicit projection of the upper hyperboloid on the plane |$(\vec{x}, \vec{y})$|⁠. This projection is a disc of radius
\begin{equation} R = 2\epsilon _{\rm max}\sqrt{\mu (1 - \mu )}, \end{equation}
(A41)
where μ is defined in equation (29). For the change of variable θ → μ, we need to switch the integration on p with the integration on cos θ in (A40), with
\begin{equation} \mathrm{d}\cos \theta = 2\frac{(k-p)}{p}\frac{\epsilon _{\rm max}}{k}\mathrm{d}\mu \end{equation}
(A42)
Moreover, the shape of the domain naturally suggests to use polar coordinates in the plane |$(\vec{x}, \vec{y})$|⁠, with radius |$r = \sqrt{x^2 + y^2}$|⁠, angle ϕw and dxdy = rdrw. The differential cross-section is
\begin{eqnarray} \frac{\mathrm{d}^2\sigma }{\mathrm{d}\Omega } & = & -\frac{r_e^2}{4} \frac{p m^2}{k\epsilon _{\rm max}^2\mu ^2} \left[ \left( \frac{m^2}{4\epsilon _{\rm max}p \mu } + \frac{m^2}{4\epsilon _{\rm max}(k-p) \mu } \right)^2 \right. \nonumber \\ && \left.- \frac{m^2}{4\epsilon _{\rm max}p \mu } - \frac{m^2}{4\epsilon _{\rm max}(k-p) \mu } - {\frac{1}{4}}\frac{p}{k-p} - {\frac{1}{4}}\frac{k-p}{p} \right].\nonumber\\ \end{eqnarray}
(A43)
The elementary current in equation (8) is
\begin{equation} \frac{K_s\cdot K_w}{K_s^0 K_w^0} = 1 - \cos \xi + \bigcirc \left(\beta _\perp /\beta _\Vert \right), \end{equation}
(A44)
where
\begin{equation} \cos \xi = 1 - 2\mu \frac{\epsilon _{\rm max}}{\epsilon } + \bigcirc \left(1\right). \end{equation}
(A45)
The expression of |$\epsilon = \sqrt{x^2 + y^2 + z^2}$| needs as well to be developed as a function of x, y and p, which implies to write a clear expression for z(p). From equation (A30),
\begin{equation} z = - z_0 + z_0 \sqrt{ 1 + \frac{x^2 + y^2}{L^2} - \delta ^2 }. \end{equation}
(A46)
One can show that under approximations in equation (23), |$\frac{x^2 + y^2}{L^2} < 1/16$| and should be in practice much smaller. Therefore,
\begin{eqnarray} z \simeq \frac{z_0}{2} \left( \frac{x^2 + y^2}{L^2} - \delta ^2 \right) \end{eqnarray}
(A47)
This allows us to make ε explicit,
\begin{equation} \epsilon = \frac{1}{4\mu \epsilon _{\rm max}}\left(x^2 + y^2 + 4\mu ^2 \epsilon _{\rm max}^2 \right), \end{equation}
(A48)
as well as the derivative of z:
\begin{equation} \frac{\mathrm{\partial} z}{\mathrm{\partial} p} = {\frac{k^2 }{2(k-p)p}\left(\frac{m^2}{kp} - \frac{2\epsilon _{\rm max}\mu }{k}\right) }{\left(1 + \frac{r^2}{4\epsilon _{\rm max}^2\mu ^2}\right) }. \end{equation}
(A49)
We note
\begin{eqnarray*} \left(\frac{\mathrm{\partial} z}{\mathrm{\partial} p}\right)_{{\rm left}} &=& {\frac{k^2 }{2(k-p)p}\left(\frac{m^2}{kp} - \frac{2\epsilon _{\rm max}\mu }{k}\right) },\nonumber \\ \left(\frac{\mathrm{\partial} z}{\mathrm{\partial} p}\right)_{{\rm right}} &=& {\left(1 + \frac{r^2}{4\epsilon _{\rm max}^2\mu ^2}\right) }.\nonumber \end{eqnarray*}
We need the absolute value of ∂z/∂p, and one can show from (A49) that it is always negative, so that we shall always take the opposite of equation (A49) and remove the absolute value in the following developments. The distribution function fw is the only element that depends on ϕw. Moreover, the integration over the hyperboloid L leads to get rid of z through equation (A47). For further developments, we explicitly keep track of the fact that L = L(p, cos θ) = L(p, μ)
\begin{equation} F_{\rm w}(r, p, \mu ) = \int _{\phi _{\rm w} = 0}^{2\pi }f_{\rm w}\left(r, \phi _{\rm w}, z(r^2,\mu )\right) \mathrm{d}\phi _{\rm w}. \end{equation}
(A50)
We can separate the integration of (28) in several parts. The parts with a dependance on r are to be found in the Jacobian |∂z/∂p| (see equation A49) of which we take only the rightmost factor, the current in equation (A44), the distribution function, and the differential element rdr. Parts that depend only on μ or p are the differential cross-section in equation (A43), the two first factors in the Jacobian |$\left|\frac{\mathrm{\partial} z}{\mathrm{\partial} p} \right|$| in equation (A49), and the Jacobian in equation (A42). The dependance of the integrated distribution function Fw is not known a priori. We obtain
\begin{eqnarray} W_{\vec{k}} &=& c{2\pi }\int _{p_1}^{p_2} \mathrm{d}p \int _{\mu _1}^{\mu _2} \mathrm{d}\mu \frac{\mathrm{\partial} \cos \theta }{\mathrm{\partial} \mu } \frac{\mathrm{d}^2\sigma }{\mathrm{d}\Omega } \left|\left.\frac{\mathrm{\partial} z}{\mathrm{\partial} p}\right|_{{\rm left}} \right| \nonumber\\ &&\times \int _{r = 0}^{R} \left.\frac{\mathrm{\partial} z}{\mathrm{\partial} p}\right|_{{\rm right}} \frac{K_{\rm s}\cdot K_{\rm w}}{K_{\rm s}^0 K_{\rm w}^0} F_w(r, p, \mu )r \mathrm{d}r. \end{eqnarray}
(A51)
The boundary conditions on the p integral, (pi)i = 1, 2 must be understood as pi = max (pi, k − pi) in virtue of symmetry (A1). At lowest order, one can show that the r part of the integrand is merely equal to 2Fw(r, p, μ)rdr and that the μ part can be reduced after a partial fraction decomposition to
\begin{equation} \sum _{i=1}^{4} \frac{a_i(p)}{\mu ^i}, \end{equation}
(A52)
where the dimensionless coefficients ai(p) are given by equation (31). Then equation (A51) can be formally reduced to equation (38).

APPENDIX B: FORMULA COMPENDIUM

The cosine of the angle θ between the strong photon |$\vec{k}$| and the outgoing lepton |$\vec{p}$| is parametrized below by
\begin{equation} \cos \theta = 1 - c \end{equation}
(B1)
and the following parametrization by μ can lead to significant simplifications
\begin{equation} c = 2\frac{k-p}{p} \frac{\epsilon _{\rm max}}{k}\mu - \frac{m^2}{2p^2}. \end{equation}
(B2)
The following quantities are used as intermediates in the derivation of the hyperboloid of integration,
\begin{eqnarray} A & = & K\cdot P \end{eqnarray}
(B3)
\begin{eqnarray} \vec{\beta } & = & \vec{k} - \vec{p} \end{eqnarray}
(B4)
\begin{eqnarray} \alpha & = & k - P^0 \end{eqnarray}
(B5)
\begin{eqnarray} \overline{\beta } & = & \left(\begin{array}{ccc}\beta _x^2 & 0 & 0 \\ 2\beta _x\beta _y & \beta _y^2 & 0 \\ 2\beta _x\beta _z & 2\beta _y\beta _z & \beta _z^2 \end{array}\right) \end{eqnarray}
(B6)
and can be explicited to relevant order (see A4) as
\begin{equation} \begin{array}{l}\begin{array}{lcl}A & = & m^2\left(\frac{1}{2}\frac{k}{p} + \frac{kp}{m^2}c\right)+ \bigcirc \left(1\right) = 2m^2 \frac{k-p}{m}\frac{\epsilon _{\rm max}}{m}\mu \nonumber\\ \alpha & = & m\left( \frac{k}{m}-\frac{p}{m} - \frac{m}{2p} \right) + \bigcirc \left(2\right) \nonumber\\ \end{array}\\ \begin{array}{rcl}\left|\beta _x \right|\sim \left|\beta _y \right| \sim \sqrt{\beta _x^2 + \beta _y^2} = \beta _\perp & = & p\sqrt{2c} + \bigcirc \left(3\right) \nonumber\\ \beta _z = \beta _\Vert & = & p\left(\frac{k}{p} - 1 + c\right). \end{array} \end{array} \end{equation}
(B7)
The characteristics of the hyperboloid (A30), are then related to the previous quantities by
\begin{eqnarray} z_0 & = & \displaystyle\frac{A\beta _\Vert }{\beta _\Vert ^2 - \alpha ^2} = \displaystyle\frac{k}{2(k - p)} (k - p + pc) = \displaystyle\frac{k}{2} + \bigcirc \left(1\right) \nonumber\\ L^2 & = & \displaystyle\frac{A^2}{\alpha ^2}\displaystyle\frac{\beta _\Vert ^2}{\beta _\Vert ^2 - \alpha ^2} = \displaystyle\frac{pk^2}{4(k - p)} \left(2c + \displaystyle\frac{m^2}{p^2} \right) = \mu k\epsilon _{\rm max}\nonumber\\ \delta ^2 & = & \displaystyle\frac{\beta _\Vert ^2 - \alpha ^2}{\beta _\Vert ^2} = \displaystyle\frac{p}{k - p + 2cp} \left(2c + \displaystyle\frac{m^2}{p^2} \right) = 4\mu \displaystyle\frac{\epsilon _{\rm max}}{k} {+}\bigcirc \left(2\right),\nonumber\\ \end{eqnarray}
(B8)
where c = 1 − cosθ.
From this, one finds the derivative of z:
\begin{equation} \frac{\mathrm{\partial} z}{\mathrm{\partial} p} = \frac{k^2 }{2(k-p)p}\left(\frac{m^2}{kp} - \frac{2\epsilon _{\rm max}\mu }{k}\right) \left(1 + \frac{r^2}{4\epsilon _{\rm max}^2\mu ^2}\right) \end{equation}
(B9)