Inverse Compton Cascades in Pair-Producing Gaps: Effects of Triplet Pair Production

TL;DR

This paper studies how triplet pair production influences inverse Compton cascades in accelerating gaps near supermassive black holes, providing analytical models and implications for gamma-ray observations of active galactic nuclei.

Contribution

It introduces a detailed analysis of triplet pair production effects on pair cascades in accelerating gaps, with analytical expressions for growth rates and electric fields.

Findings

Triplet pair production can regulate pair cascades in certain astrophysical environments.

Analytical formulas for cascade growth rate and electric field are derived.

Implications for gamma-ray emissions from active galactic nuclei are discussed.

Abstract

Inverse Compton-pair cascades are initiated when gamma-rays are absorbed on an ambient soft photon field to produce relativistic pairs, which in turn up-scatter the same soft photons to produce more gamma-rays. If the Compton scatterings take place in the deep Klein-Nishina regime, then triplet pair production ($e\gamma_b \rightarrow ee^{+}e^{-}$) becomes relevant and may even regulate the development of the cascade. We investigate the properties of pair-Compton cascades with triplet pair production in accelerating gaps, i.e., regions with an unscreened electric field. Using the method of transport equations for the particle evolution, we compute the growth rate of the pair cascade as a function of the accelerating electric field in the presence of black-body and power-law ambient photon fields. Informed by the numerical results, we derive simple analytical expressions for the peak growth rate and the corresponding electric field. We show that for certain parameters, which can be realized in the vicinity of accreting supermassive black holes at the centers of active galactic nuclei, the pair cascade may well be regulated by inverse Compton scattering in the deep Klein-Nishina regime and triplet pair production. We present indicative examples of the escaping gamma-ray radiation from the gap, and discuss our results in application to the TeV observations of radio galaxy M87.

Equations

$${\gamma }_{eq}={\left(\frac{3q_e{E}_0}{4{\sigma }_{\mathrm{T}}{U}_{ph}}\right)}^{1/2}.$$

$$E_0^{(\mathrm{KN})}=\frac{4{\sigma }_{\mathrm{T}}{U}_{ph}}{3q_e}{\left(\frac{m_e{c}^2}{⟨ϵ⟩}\right)}^2,$$

$$\frac{\partial n_e^{(1)}}{\partial t}+q_{e}c{E}_0\frac{\partial n_e^{(1)}}{\partial \gamma }+{\mathcal{L}}_e^{ics}+{\mathcal{L}}_e^{tpp}$$

$$\frac{\partial n_e^{(2)}}{\partial t}+q_{e}c{E}_0\frac{\partial n_e^{(2)}}{\partial \gamma }+{\mathcal{L}}_e^{ics}$$

$$\frac{\partial n_{\gamma }}{\partial t}+{\mathcal{L}}_{\gamma }^{\gamma \gamma }$$

$$\displaystyle\langle\epsilon\rangle=\left\{\begin{array}[]{ll}\epsilon_{\min}\ln\left(\frac{\epsilon_{\max}}{\epsilon_{\min}}\right)\left(1-\frac{\epsilon_{\min}}{\epsilon_{\max}}\right)^{-1},&\Gamma=2\\ \\ \frac{-\Gamma+1}{-\Gamma+2}\frac{\epsilon_{\max}^{-\Gamma+2}-\epsilon_{\min}^{-\Gamma+2}}{\epsilon_{\max}^{-\Gamma+1}-\epsilon_{\min}^{-\Gamma+1}},&\Gamma\neq 2\end{array}\right.$$

$${\ell }_{\mathrm{IC}}(ϵ)\approx \frac{1}{{\sigma }_T{n}_b(ϵ)ϵ}=\frac{{ϵ}^{\Gamma -1}}{{\sigma }_T{n}_0{ϵ}_{\min }^{\Gamma }}.$$

$${\ell }_{\gamma \gamma }({ϵ}_t)\sim \frac{1}{{\sigma }_{\gamma \gamma }{n}_b({ϵ}_t){ϵ}_t}=\frac{{ϵ}_t^{\Gamma -1}}{{\sigma }_{\gamma \gamma }{n}_0{ϵ}_{\min }^{\Gamma }}.$$

$${\gamma }_{eq}{ϵ}_{\min }\lesssim m_e{c}^2{\left(\frac{4/3}{\Gamma -2}\right)}^{1/(\Gamma -2)},\Gamma >2$$

$$\ln \left(\frac{{ϵ}_{\max }}{{ϵ}_{\min }}\right)\gtrsim 3/4,\Gamma =2.$$

$$\kappa \propto \frac{c}{{\ell }_{\mathrm{IC}}(ϵ=m_e{c}^2/{\gamma }_{eq})}\propto {\gamma }_{eq}^{\Gamma -1}\propto E_0^{\frac{\Gamma -1}{2}},$$

$$\kappa \propto \frac{c}{{\ell }_{tpp}(\gamma ={\gamma }_{\mathrm{coff}}^{(tpp)})}\propto \ln E_0,$$

$$\lambda \approx {\lambda }_{\max }{\left(\frac{E_0}{E_0^{(\mathrm{KN})}}\right)}^{\frac{\Gamma -1}{2}},E\le E_0^{(\mathrm{KN})}$$

$${\lambda }_{\max }\approx 481\frac{U_{ph,6}}{⟨ϵ{⟩}_0}[{\mathrm{s}}^{-1}],$$

$$t_E\approx \frac{E_0}{\dot{E}}\approx {\left(2\pi \frac{E_0}{B_0}\right)}^{1/2}\frac{r_g}{c},$$

$$⟨ϵ{⟩}_0^2{U}_{ph,0}^{-3/2}\gtrsim 0.02M_9{B}_{0,2}^{-1/2},$$

$$P_{\mathrm{gap}}\sim \frac{16E_0}{B_0}\frac{h}{r_g}{P}_{\mathrm{jet}}.$$

$${\gamma }_{\mathrm{cur}}\simeq 10^{10}{E}_{0,2}^{1/4}{M}_9^{1/4}{\left(\frac{{\rho }_c}{r_g}\right)}^{1/2}$$

$${\gamma }_{\mathrm{coff}}^{(tpp)}$$

$$⟨ϵ{⟩}_0{U}_{ph,0}^{-1}\lesssim 7\times 10^{-4}{E}_{0,2}^{-3/4}{M}_9^{1/4}{\left(\frac{{\rho }_c}{r_g}\right)}^{1/2}.$$

$$\frac{E_0}{B_0}>{\left(\frac{3}{8\pi }\right)}^{1/2}{\left(\frac{r_B}{r_g}\right)}^{1/2}{\left(\frac{{\ell }_{\mathrm{IC}}^{\mathrm{tot}}}{r_g}\right)}^{1/2}{\left(\frac{m_e{c}^2}{⟨ϵ⟩}\right)}^{1/2},$$

$$\frac{E_0}{B_0}\gtrsim 0.01M_9^{-1}{U}_{ph,0}^{-2}⟨ϵ{⟩}_0^2,$$

Full text

acceleration of particle – gamma-rays: galaxies – galaxies: active, individual (M87) – radiation mechanisms: non-thermal

††journal: ApJ

1 Introduction

The production of a relativistic electron-positron pair from the absorption of a gamma-ray photon is the starting point of an avalanche of pairs and gamma rays, known as an electromagnetic (EM) cascade. Inverse Compton-pair cascades, in particular, are initiated when gamma-rays are absorbed on a soft photon field to produce relativistic pairs, which in turn inverse Compton scatter off the same soft photons to produce more gamma-rays.

In the simplest case, where the pairs are not accelerated, only a few generations of them are expected to be produced before the energy of gamma-rays is degraded to a value below the threshold for $\gamma\gamma$ pair creation, thus ceasing the EM cascade. At this point, gamma-rays can escape unimpeded from the region of their production, with their initial energy being shared among the secondary pairs and photons. As long as the energy density of the secondary photons is lower than that of the background soft photon field, the cascade is considered to be linear, i.e., the interactions between secondary pairs and photons are negligible.

Linear inverse Compton-pair cascades in isotropic photon fields111The development of these cascades in anisotropic photon fields was studied by Protheroe et al. (1992). have been discussed in the context of compact X-ray sources (Akharonian et al., 1985; Mastichiadis, 1986; Zdziarski et al., 2009) and applied to the propagation of ultra-high-energy (UHE) gamma-rays in the intergalactic medium (IGM) (Gould & Rephaeli, 1978; Protheroe, 1986; Zdziarski, 1988; Coppi & Aharonian, 1997). In the latter case, the initial Compton scatterings of UHE pairs on low-energy background photons (e.g., cosmic microwave background) can take place deep in the Klein-Nishina regime (i.e., $s\gg(m_{e}c^{2})^{2}$, where $s$ is the squared center-of-mass energy of the interaction). In this regime, triplet pair production (TPP; $e\gamma_{b}\rightarrow ee^{+}e^{-}$) can be an additional source of electron cooling. Although TPP is a third-order quantum electrodynamical process, it can operate in direct competition to inverse Compton scattering (ICS) for $s\gtrsim 100$ MeV2 (Haug, 1981; Mastichiadis et al., 1986; Mastichiadis, 1991; Anguelov et al., 1999). Because of the very high interaction energies required, the role of TPP has been mainly investigated in the context of UHE-photon induced cascades in the IGM (Lee, 1998; Settimo & De Domenico, 2015; Wang et al., 2017).

The vicinity of supermassive black holes in active galactic nuclei (AGN), where strong electric fields may be present and copious soft photons are emitted by the accretion disk, is another astrophysical environment where linear Compton-pair cascades have been studied (e.g., Bednarek & Kirk, 1995; Levinson & Rieger, 2011; Broderick & Tchekhovskoy, 2015; Hirotani & Pu, 2016; Chen et al., 2018; Katsoulakos & Rieger, 2018; Levinson & Cerutti, 2018). In this case, the physics of the cascade are coupled to the acceleration of particles by the component of the electric field parallel to the magnetic field over the length of a charge-starved region, the so-called gap222There is an analogy to pair cascades in pulsar magnetospheres (Ruderman & Sutherland, 1975), although the physical processes related to the production of gamma-ray photons are different (see e.g., Timokhin & Arons, 2013; Timokhin & Harding, 2019, and references therein).. A seed electron that enters the gap can accelerate to high enough energies to up-scatter disk photons and produce gamma-rays, which in turn pair-produce on the disk photons, thus initiating a pair cascade in the gap. In contrast to the cascades induced in the IGM, the secondary electrons are continuously accelerated inside the gap. This leads to an exponentiation of their number density until they provide the necessary charge density to sort out the electric field and quench acceleration.

In this paper, we study the generic properties of inverse Compton cascades in pair-producing gaps in the Thomson and Klein-Nishina regimes with the inclusion of triplet pair production. More specifically, we compute the growth rate of secondary pairs as a function of the electric field strength for a background soft photon field (with a black-body or power-law energy distribution) and study the effects of triplet pair production on the development of these cascades (e.g., growth rate and maximum particle energy). This paper is structured as follows. In Section 2 we qualitatively discuss the general properties of Compton-pair cascades with triplet pair production. In Section 3 we outline the adopted numerical scheme for the calculation of the cascade and present our results on the pair growth rate in Section 4. We discuss their astrophysical implications in Section 5 and conclude with a discussion on the model’s caveats in Section 6.

2 General considerations

We consider a stationary rectilinear particle accelerator immersed in an isotropic time-independent ambient photon field with energy density $U_{ph}\equiv\int d\epsilon\,\epsilon n_{b}(\epsilon)$, where $n_{b}(\epsilon)$ is the differential photon number density. The accelerator is characterized by a constant electric field with strength $E_{0}$. The electric field is assumed to be almost parallel to the local magnetic field, which has no curvature. The size of the acceleration region is taken to be much larger than all relevant mean free paths of relativistic pairs, thus allowing us to follow the full development of the pair cascade.

Within our toy model, we treat electrons and positrons as identical particles with initial Lorentz factor $\gamma_{0}$. Upon entering the acceleration region, pairs will gain energy from the electric field, and at the same time, lose energy due to ICS (or TPP) on the background soft photon field. Energy losses due to synchrotron and curvature radiation can be neglected under our assumptions (i.e., the particle pitch angle and the curvature of magnetic field lines are both approximately zero).

As long as the scatterings occur in the Thomson regime, particles can accelerate from their initial Lorentz factor to a maximum value $\gamma_{eq}$ where the continuous energy losses balance the energy gains:

[TABLE]

In the above equation $c$ is the speed of light, $\sigma_{\rm T}$ is the Thomson cross section, and $q_{e},m_{e}$ are the electron charge and mass. The scattering of soft photons with energy $\langle\epsilon\rangle\equiv\int d\epsilon\,\epsilon n_{b}(\epsilon)/\int d\epsilon\,n_{b}(\epsilon)$ by pairs with Lorentz factor $\gamma_{eq}$ will occur in the Klein-Nishina regime, if $E_{0}$ exceeds a critical value that depends solely on the properties of the target photon field, namely:

[TABLE]

where we used the condition $\gamma_{eq}\langle\epsilon\rangle\gtrsim(3/4)m_{e}c^{2}$ and Equation (1). In particular, for a black-body photon field with temperature $T$ we find $E_{0}^{\rm(KN)}\propto T^{2}$. Using this critical value of the electric field, we can distinguish the following two regimes.

Low-$E_{0}$ regime: for $E_{0}<E_{0}^{\rm(KN)}$ the acceleration process is balanced by energy losses mostly in the Thomson regime and the particle Lorentz factor saturates at the value given by Equation (1). As a result, the up-scattered photons, of energy $\epsilon_{1}\sim\gamma_{eq}^{2}\langle\epsilon\rangle$, are below the threshold for $\gamma\gamma$ production with the average energy soft photons (i.e., $\epsilon_{1}\langle\epsilon\rangle<2(m_{e}c^{2})^{2}$). Still, pair production is expected to occur at some level due to $\gamma\gamma$ interactions with photons from the high-energy tail of the distribution. We note that TPP is not relevant in this regime due to the much lower interaction rates (see e.g., Fig. 4 in Mastichiadis, 1991). The number density of photons at $\epsilon\gg\langle\epsilon\rangle$ together with the rate of acceleration will essentially determine the growth rate of the produced pairs. For a black-body photon field, there is essentially no free parameter other than $E_{0}$, since the photon number density at $\epsilon\gg\langle\epsilon\rangle\simeq 2.7k_{B}T$ is just given by the Wien part of the spectrum. However, in the more general case of a power-law photon field, the exact growth rate of the number of pairs is expected to depend on several free parameters: the maximum and minimun cutoff energies ($\epsilon_{\max}$, $\epsilon_{\min}$), the photon index ($\Gamma$), and the energy density ($U_{ph}$).

High-$E_{0}$ regime: for $E_{0}>E_{0}^{\rm(KN)}$, pairs can achieve much higher Lorentz factors than $\gamma_{eq}$ (see Equation (1)); the highest energy pairs can up-scatter soft photons in the Klein-Nishina regime, where the energy losses are not strong enough to balance the energy gains by the acceleration. As a result, the maximum particle energy does not saturate at $\gamma_{eq}$, but continues to increase till another process (e.g., physical escape from the system) settles in. Quickly upon entering the acceleration region, the seed pairs can achieve high enough energies to produce secondary pairs off the photon field either directly through TPP or indirectly through the combination of ICS and $\gamma\gamma$ absorption. In both cases, the produced pairs will also be accelerated by the strong electric field, and their number will exponentiate.

In both regimes, a pair cascade will develop in the acceleration region. The cascade is, in good approximation, one-dimensional as long as the magnetic deviations of the particle trajectories over the relevant particle interaction lengths are smaller than the emission cone of the pairs. Our goal is to compute the growth rate of the ensued linear pair cascade as a function of $E_{0}$ for different ambient photon fields. In the following section, we describe the numerical approach we adopted for the study of the pair cascade.

3 Numerical approach

The physical problem outlined in Section 2 can be mathematically formulated using the concept of transport equations that describe the evolution of three particle populations, namely the seed pairs that enter the acceleration region, the secondary pairs produced by $\gamma\gamma$ and TPP processes, and the photons emitted by both pair populations via ICS:

[TABLE]

Here, $n^{(1)}_{e},n^{(2)}_{e}$, and $n_{\gamma}$ denote respectively the differential number densities of seed pairs, secondary pairs, and photons produced by both populations of pairs. The operators $\mathcal{Q}_{i}^{j}$ and $\mathcal{L}_{i}^{j}$ denote injection and loss terms of the particle species $i$ due to the process $j$.

To solve the particle transport equations, we employ the numerical code of Mastichiadis & Kirk (1995), but with two main modifications: the implementation of the triplet pair production process and the use of a more accurate scheme for the treatment of electron losses in the Klein-Nishina regime. For completeness, we summarize the different terms appearing in eqs. (3)-(5) briefly.

Pair injection. Seed pairs with some low Lorentz factor are injected into the acceleration region at a rate $Q_{e}^{inj}=Q_{0}\delta(\gamma-\gamma_{0})H(t)$. Here, we adopted (without loss of generality) $\gamma_{0}=1$. We checked that the derived exponential growth rates (see Section 4) do not depend on the choice of $\gamma_{0}$ or of the injection profile (e.g., continuous versus instantaneous injection).

Inverse Compton scattering. The term $\mathcal{L}^{ics}_{e}$ takes into account energy losses in both the Thomson and Klein-Nishina regimes according to Eqs. (5.7) and (5.17) of Blumenthal & Gould (1970). For the photon production rate $\mathcal{Q}^{ics}_{\gamma}$ we use Eq. (2.48) of Blumenthal & Gould (1970).

$\gamma\gamma$ pair production. The pair injection term $\mathcal{Q}^{\gamma\gamma}_{e}$ and photon attenuation rate $\mathcal{L}^{\gamma\gamma}_{\gamma}$ are given respectively by Eqs. (57) and (54) of Mastichiadis & Kirk (1995).

Triplet pair production. We consider the energy losses of the recoiling electron (or positron) to be continuous. For the single-particle loss rate $q_{tpp}$ we adopt the Eq. (29) of Mastichiadis (1991), which holds away from the interaction threshold ($\gamma\epsilon>10^{3}m_{e}c^{2}$). The exact functional form of $q_{tpp}$ near the threshold is not important, since ICS losses dominate in this energy regime (Mastichiadis, 1991). The loss term appearing in Equation (3) is modeled as $\mathcal{L}_{e}^{tpp}(\gamma)=\int d\epsilon\,n_{b}(\epsilon)q_{tpp}(\gamma\epsilon/m_{e}c^{2})$. Following Lee (1998), we model the TPP production rate as $\mathcal{Q}_{e}^{tpp}(\gamma^{\prime})=q_{0}\gamma^{{}^{\prime}-\delta}\int d\gamma\,n^{(1)}_{e}(\gamma)\int d\epsilon\,\sigma_{tpp}(\gamma\epsilon/m_{e}c^{2})n_{b}(\epsilon)$, where $\sigma_{tpp}$ is the total TPP cross section (Haug, 1981; Mastichiadis, 1991). We adopt the asymptotic value of $7/4$ for the exponent $\delta$ (see also Fig. 2 in Mastichiadis, 1991), while noting that our results are not sensitive on the exact value of $\delta$. The normalization constant $q_{0}$ is determined by the requirement of energy conservation, i.e., $\int d\gamma\,n^{(2)}_{e}(\gamma)\gamma\mathcal{Q}_{e}^{tpp}(\gamma)=\int d\gamma\,n^{(1)}_{e}(\gamma)\gamma\mathcal{L}_{e}^{tpp}(\gamma)$.

Finally, we note that the numerical code solves for the combined distribution $n_{e}=n_{e}^{(1)}+n_{e}^{(2)}$. Although the code does not differentiate between seed pairs and secondary pairs, any exponentially growing solutions of $n_{e}$ are purely determined by the secondaries. This can be easily understood by considering the extreme case of no pair production ($n_{e}^{(2)}=0$). In this case, one can show by solving Equation (3) with $\mathcal{L}^{tpp}_{e}=0$ that the lepton number $\int n_{e}(\gamma)d\gamma$ increases only linearly with time.

4 Results

Using the numerical code described in the previous section, we compute the growth rate ($\lambda$) of exponentially increasing solutions for a wide range of electric field strengths $E_{0}$ and different ambient photon fields. More specifically, we perform a linear fit to $y(t)\equiv\ln\left[\int d\gamma\,n_{e}(\gamma,t)\right]$, since its slope gives directly the rate of exponentially growing solutions. In what follows, we compute $\lambda$ for a black-body photon field of temperature $T$ (Section 4.1) and for a power-law photon field with photon index $\Gamma$ and energy density $U_{ph}$, extending from $\epsilon_{\min}$ to $\epsilon_{\max}$ (Section 4.2).

4.1 Black-body photons

Figure 1 presents the exponential growth rate $\lambda$ computed numerically for different electric field strengths $E_{0}$ when: (i) both $\gamma\gamma$ and triplet pair production processes are taken into account (filled symbols), (ii) only triplet pair production is taken into account (open triangles), and (iii) only $\gamma\gamma$ pair production is included (open circles). Note that the growth rate in the absence of pair production asymptotes to a constant value instead of approaching zero (see e.g., filled symbols at $E_{0}\lesssim 10^{-2}$ statV cm*-1* and open circles at $E_{0}>10^{3}$ statV cm*-1*). This is a result of our choice for the injection of seed particles (see also Section 3). With a continuous injection, even in the absence of any pair production processes ($n_{e}^{(2)}=0$), the number of seeds is expected to increase linearly with time. These non-exponential growing solutions are reflected to the non-zero values of the numerically derived growth rates333We verified that $\lambda\rightarrow 0$ for sufficiently low or high (in the absence of TPP) electric field strengths, if seed particles are instantaneously injected into the gap.

The actual shape of the curve $\lambda(E_{0})$ is determined by the convolution of the ambient photon field energy distribution with the energy-dependent cross sections of the relevant physical processes. In the low-$E_{0}$ regime (i.e., for $E_{0}<E_{0}^{\rm(KN)}$), the growth rate is determined essentially by the $\gamma\gamma$ pair production process. However, in the high-$E_{0}$ regime ($E_{0}>E_{0}^{\rm(KN)}$), where the $\gamma\gamma$ cross section decreases (see open circles), triplet pair production begins to take over and eventually dominates the total growth rate of pairs at $\gtrsim 100\,E_{0}^{\rm(KN)}$ (open triangles). In this regime, the slow increase of $\lambda$ with increasing $E_{0}$ reflects the logarithmic dependence of the TPP cross section (see also Figure 3 and text in Section 4.2, for more details).

As an indicative example, we show in Figure 2 (left panel) the temporal evolution of the pair energy distribution computed for $E_{0}=1$ staV cm${}^{-1}\simeq E_{0}^{\rm(KN)}$, when both pair production processes are taken into account. The dashed vertical line marks the saturation Lorentz factor assuming ICS losses in the Thomson regime (see Equation (1)). The energy spectrum can be modeled as $\gamma n_{e}(\gamma)\propto\gamma e^{-\gamma/\gamma_{\rm coff}}$, where $\gamma_{\rm coff}$ is generally time-dependent. Before the number of pairs exponentiates (i.e., $t\lesssim\lambda^{-1}$), we find that $\gamma_{\rm coff}\propto t$ and $\gamma_{\rm coff}\gtrsim\gamma_{eq}$. The extension of the spectrum beyond $\gamma_{eq}$ suggests that the up-scattering of black-body photons with energy $\langle\epsilon\rangle=2.7k_{B}T$ by pairs with $\gamma\gtrsim\gamma_{eq}$ takes place in the Klein-Nishina regime, already from $E_{0}\sim E_{0}^{\rm(KN)}$.

Soon after the onset of the cascade (i.e., $t\gtrsim\lambda^{-1}$), the cutoff Lorentz factor of the distribution freezes at a value that depends on $E_{0}$. This is illustrated in the middle panel of Figure 2, where we plot the time-independent $\gamma_{\rm coff}$ determined numerically for different electric field strengths. The dashed red line denotes the saturation Lorentz factor considering particle losses only in the Thomson regime (see Equation (1)). As soon as pair production becomes important, we find that $\gamma_{\rm coff}>\gamma_{eq}$. For the adopted parameters, this occurs at $E_{0}\gtrsim 10^{-2}$ statV cm*-1* (see also Fig. 1). In the high-$E_{0}$ regime, in particular, the dependence of $\gamma_{\rm coff}$ on $E_{0}$ can be derived from the following condition: any pairs injected via the TPP process should reach the Lorentz factor $\gamma$ of the parent electron within the mean free path of the process, namely $m_{e}\gamma c^{2}\approx eE_{0}\ell_{tpp}(\gamma)$ with $\ell^{-1}_{tpp}(\gamma)\approx\int d{\epsilon}\,n_{b}(\epsilon)\sigma_{tpp}(\gamma\epsilon/m_{e}c^{2})$. For a black-body photon distribution, we find $\ell_{tpp}(\gamma)\approx$ constant and $\gamma_{\rm coff}^{(tpp)}\propto E_{0}$, in agreement with the numerical results for $E_{0}\gtrsim 30$ statV cm*-1*, where the pair production is regulated via the TPP process (see also Figure 1).

The right panel of Figure 2 shows snapshots of the photon energy spectra produced by ICS off an ambient black-body photon field by the pairs in the accelerating gap (see left panel of the same figure). For reference, the typical photon energy of $4\gamma_{eq}^{2}k_{B}T$ is also indicated (dashed grey line). Each snapshot corresponds to a certain time (in units of $\lambda^{-1}$) or, equivalently, a certain location within the gap. If the physical size of the latter is taken to be $h=14\,c/\lambda=10^{10}$ cm, then the last snapshot shown in the figure would correspond to the spectrum of photons escaping from the gap. Note that this would not necessarily be the observed spectrum, because other processes operating outside the gap could still alter it (for further discussion, see Section 5.1).

4.2 Power-law photons

Figure 3 presents the exponential growth rate $\lambda$ as a function of $E_{0}$ for an ambient photon field with a power-law energy distribution. Here, both $\gamma\gamma$ and triplet pair production are considered. Different panels (from left to right) demonstrate the effects of $U_{ph},\epsilon_{\min}$, and $\Gamma$ on the growth rate (for details, see figure caption and inset legends). By varying the photon energy density, while keeping all other parameters fixed (left panel), we obtain a self-similar family of curves for $\lambda(E_{0})$. For example, if we were to shift the black curve to the right and upwards by a factor of $10$ in each direction, we would obtain the blue curve. More specifically, the peak growth rate $\lambda_{\max}$ and the electric field strength at the peak scale linearly with $U_{ph}$. The latter scaling agrees with the predicted one for $E_{0}^{(\rm KN)}$ (see Equation (2)).

The growth rate depends also on $\epsilon_{\min}$, as shown in the middle panel, for fixed $U_{ph}$ and $\Gamma$. In particular, we find larger growth rates for lower $\epsilon_{\min}$, as the latter translates to a higher number density of low energy photons, acting as targets for ICS and pair production. Were the black curve shifted to the right by a factor of $\sim 30^{2}$ and upwards by a factor of $\sim 30$, it would match the blue curve. Thus, for fixed $U_{ph}$ and $\Gamma$ we find $\lambda_{\max}\propto\epsilon_{\min}^{-1}$. This also translates to $\lambda_{\max}\propto\langle\epsilon\rangle^{-1}$, where the mean photon energy for a power-law photon field is:

[TABLE]

Finally, the effects of the photon index on the growth rate are presented in the right panel of Figure 3. We find that changes on $\Gamma$ are imprinted on both the peak growth rate and the slope of the power-law segment of $\lambda(E_{0})$ at $E_{0}<E_{0}^{\rm(KN)}$. By changing the photon index and keeping all other parameters fixed, we are effectively varying the mean photon energy of the power-law distribution. This explains why the peak electric field $\approx E_{0}^{\rm(KN)}$ scales as $\langle\epsilon\rangle^{-2}$ (see Equation (2)). We also verified that the peak growth rates of the black and blue curves match, if the black curve is shifted downwards by a factor equal to the ratio of the mean photon energies. Our numerical results also suggest that $\lambda\propto E_{0}^{(\Gamma-1)/2}$ in the low-$E_{0}$ regime. In what follows, we present analytical arguments to explain the power-law dependence of the growth rate on the photon index.

The growth rate is proportional to the number of generations of pairs produced per unit time. The latter, denoted as $\kappa$, can be estimated as $\kappa\approx c/(\ell_{\rm acc}+\ell_{\rm IC}+\ell_{\gamma\gamma})$, where $\ell_{\rm acc}$ is the acceleration length, $\ell_{\rm IC}$ is the IC scattering mean free path, and $\ell_{\gamma\gamma}$ is the $\gamma\gamma$ collision mean free path. When ICS is in Thomson regime, particles reach a terminal Lorentz factor as determined from Equation (1). As these particles IC scatter soft photons from different segments of the power-law energy spectrum, their mean free paths will also be different. Let $n_{b}(\epsilon)=n_{0}(\epsilon/\epsilon_{\min})^{-\Gamma}$ be the differential number density of the ambient photon field. Then, the mean free path for ICS on photons of energy $\epsilon$ is roughly:

[TABLE]

The ICS on these soft photons produces high energy photons of energy $\epsilon_{1}\approx\gamma_{eq}^{2}\epsilon$, which can then interact with soft photons to produce pairs. The most likely target soft photons for $\gamma\gamma$ pair production are those near the threshold of the interaction, namely their energy is $\epsilon_{t}\approx(m_{e}c^{2})^{2}/\epsilon_{1}$, and the $\gamma\gamma$ mean free path is:

[TABLE]

The peak cross section of the $\gamma\gamma$ pair production process (near the threshold) is similar to the Thomson cross section (\sigma_{\gamma\gamma}$$\approx\sigma_{T}/5, Gould & Schréder, 1967). So, if the target photons for ICS have low $\epsilon$, then $\ell_{\rm IC}(\epsilon)$ is small, but $\ell_{\gamma\gamma}(\epsilon_{t})$ of the most-likely-to-be attenuated photons is large, and vice versa. As a result, pair production is most efficient on target photons that have similar mean free paths for ICS and $\gamma\gamma$ pair production. This essentially means that $\epsilon\approx\epsilon_{t}\approx m_{e}c^{2}/\gamma_{eq}$.

Moreover, in the low-$E_{0}$ regime, the acceleration length $\ell_{\rm acc}$ is typically much smaller than the interaction mean free paths. The condition $\ell_{\rm acc}(\gamma_{eq})\lesssim\ell_{\rm IC}(\epsilon=m_{e}c^{2}/\gamma_{eq})$ can be written as:

[TABLE]

Noting that $\gamma_{eq}\epsilon_{\min}<\gamma_{eq}\langle\epsilon\rangle<m_{e}c^{2}$, both conditions are easily satisfied in the low-$E_{0}$ regime. Thus, the generation rate of pairs can be estimated by:

[TABLE]

which roughly agrees with the slope we obtain from the numerical calculations (see right panel in Figure 3).

In the high-$E_{0}$ regime, and in particular at $E_{0}\gg E_{0}^{\rm(KN)}$, pairs are injected into the acceleration region via the TPP process. The mean free path of a high-energy electron with Lorentz factor $\gamma<2m_{e}c^{2}/\epsilon_{\min}$ is $\propto\gamma^{-\Gamma+1}$, while it has only a logarithmic dependence on $\gamma$ otherwise, i.e., $\ell_{tpp}(\gamma)\propto\ln^{-1}\left(2\gamma\epsilon_{\min}/m_{e}c^{2}\right)$. The latter case is typically satisfied in the high-$E_{0}$ regime, where pairs are ultra-relativistic and their acceleration length is not anymore negligible. Pair production is most efficient when $\ell_{\rm acc}(\gamma)\approx\ell_{tpp}(\gamma)\Rightarrow\gamma\propto E_{0}$, where the latter relation is derived after neglecting the logarithmic dependence of the TPP mean free path on the Lorentz factor. This condition also gives an estimate of the cutoff Lorentz factor of the pair distribution444When TPP is the dominant source of secondary pairs, the cutoff Lorentz factor of the pair distribution can be derived by $\ell_{\rm acc}(\gamma_{\rm coff})\approx\ell_{tpp}(\gamma_{\rm coff})$. Thus, $\gamma_{\rm coff}^{(tpp)}\propto E_{0}^{1/\Gamma}<2m_{e}c^{2}/\epsilon_{\min}$ and $\gamma_{\rm coff}^{(tpp)}\propto E_{0}>2m_{e}c^{2}/\epsilon_{\min}$. Typically, only the latter branch is relevant, since at low Lorentz factors pair production is regulated by ICS and $\gamma\gamma$ processes., as discussed in Sec. 4.1 for a black-body photon field (see also Fig. 2). The number of pairs in the deep high-$E_{0}$ regime is therefore estimated as:

[TABLE]

which agrees with the logarithmic scaling of the growth rate determined numerically (see dash-dotted line in the left panel of Figure 3).

Based on the self-similarity of the numerical results shown in Figure 3 and the analytical arguments for the dependence of the growth rate on the photon index in the low-$E_{0}$ regime, we may write the following analytical expression for the growth rate:

[TABLE]

and its maximum value:

[TABLE]

where $U_{ph}\equiv 10^{6}\,U_{ph,6}$ erg cm*-1* and $\langle\epsilon\rangle\equiv 10^{0}\,\langle\epsilon\rangle_{0}$ eV.

5 Astrophysical applications

A funnel-like region forming at the interface of a jet with the accretion flow of a black hole (e.g., Bednarek & Kirk, 1995) or a magnetospheric gap (e.g., Blandford & Znajek, 1977; Beskin et al., 1992; Hirotani & Okamoto, 1998) could be considered astrophysical analogues of our simplified particle accelerator.

In the case of a gap in the black-hole magnetosphere, it has been shown that the development of a non-ideal electric field and its screening can be highly time variable (Chen et al., 2018). Let us consider a flux tube threading the black hole horizon. At small distances within the tube particles fall into the black hole, while at large distances particles are flung outward by magneto-centrifugal effects. As a result, charge carriers get gradually depleted in the flux tube, leading to the growth of a non-ideal electric field when particles are not enough to carry the necessary current. Before any significant pair production happens, the electric field increases almost quadratically with time, so the time needed for the electric field to grow up to $E_{0}$ can be approximated as (Chen et al., 2018):

[TABLE]

where $r_{g}=2GM_{\rm BH}/c^{2}$ is the gravitational radius of a black hole with mass $M_{\rm BH}$ and $B_{0}$ is the strength of the poloidal magnetic field near the black-hole horizon. If the acceleration time of the particles is much shorter than the timescale for electric field growth, then the assumption of a constant electric field adopted here is valid (see Section 6 for a discussion on the validity of this assumption). A crude criterion for the screening of the gap is that the pair exponentiation timescale becomes comparable to or shorter than the timescale for the electric field growth given by Equation (18).

In Figure 4 we show the growth timescale of the Compton-pair cascade (colored curves) for a power-law background photon field with different energy densities and minimum energies (for details, see plot legends). The cascade timescale is large at small $E_{0}$, but becomes smaller than $t_{E}$ (dashed grey line) for sufficiently large $E_{0}$. The equality of the two timescales gives an estimation of the maximum electric field that can develop in an electrostatic gap before it is screened, i.e., before the number density of pairs exceeds the Goldreich-Julian density $n_{\rm GJ}=\Omega B_{0}/2\pi ec\approx\alpha B_{0}/4\pi er_{g}\simeq 1.7\times 10^{-5}\,B_{0,2}\,r^{-1}_{g,15}$ cm*-3*, where $\alpha\sim 1$ is the dimensionless spin of the black hole and $r_{g,15}\equiv r_{g}/10^{15}\,{\rm cm}$. For certain combinations of $\langle\epsilon\rangle$ and $U_{\rm ph}$, the screening of the electric field can happen in the regime where the pair cascade is governed by the TPP process and Compton scattering in the deep Klein-Nishina regime (see e.g., black and blue curves in the left and right panels of Figure 4, respectively):

[TABLE]

where we used the condition $\lambda_{\max}^{-1}>t_{E}$, Equations (2),(17)-(18), and $M_{9}\equiv M_{\rm BH}/10^{9}M_{\odot}$.

The characteristic timescale $t_{*}\equiv t_{E}=\lambda^{-1}$ can also be used to estimate the typical size of the gap as $h=c\,t_{*}$, which lies in the range $(0.01-1)\,r_{g}$ for all the cases displayed in Figure 4. The maximum electric field that can be developed in the gap dictates also the amount of dissipation in the magnetosphere. The gap power can be estimated as $P_{\rm gap}\sim E_{0}JAh$, where $A=\pi r_{g}^{2}$ is the cross sectional area of the jet funnel, and $J\sim Bc/(4\pi r_{g})$ is the Goldreich-Julian current density required by the jet. As a comparison, the jet power is roughly $P_{\rm jet}\sim\alpha B_{0}^{2}Ac/(64\pi)$, and one can write:

[TABLE]

5.1 Application to M87

As an indicative example, we consider the radio galaxy M87 (NGC 4486) located at a distance $\simeq 16.7$ Mpc (Mei et al., 2007). Its soft photon spectrum peaks at 300 GHz, with a flux around 1 Jy (corresponding luminosity $L\sim 10^{42}\,\rm{erg\,s}^{-1}$), and an estimated photon index $\Gamma\sim 2.2$ above 300 GHz (see Broderick & Tchekhovskoy, 2015, and references therein). However, the actual amount of the soft radiation in the direct vicinity of the gap is quite uncertain. Following Chen et al. (2018), we adopt the 300 GHz frequency as the low-energy end of the power-law photon distribution, i.e., $\epsilon_{\min}=1.2\times 10^{-3}$ eV$=2.4\times 10^{-9}m_{e}c^{2}$, and use $U_{ph}=L/(4\pi r^{2}c)\sim 0.1\,\rm{erg}\,\rm{cm}^{-3}$ as an estimate of the soft photon energy density at a distance $r\approx 3r_{g}$ from the black hole. The pair cascade growth rate for these parameters is shown by the red curve in the right panel of Fig. 4. The screening of the gap occurs in the low-$E_{0}$ (Thomson) regime (see crossover of the red and dashed grey lines) and the peak electric field is predicted to be $E_{0}\sim 5\times 10^{-3}\,\rm{statV}\,\rm{cm}^{-1}\sim 5\times 10^{-5}\,B_{0}$ for $B_{0}=100$ G, consistent with the estimation by Chen et al. (2018). In reality, the properties of the radiation field at the gap (e.g., $U_{ph}$ and $\epsilon_{\mathrm{min}}$) are highly uncertain. For example, if $\epsilon_{\mathrm{min}}\gtrsim 10^{-2}\,\mathrm{eV}$ or $U_{ph}\lesssim 10^{-3}\,\rm{erg}\,\rm{cm}^{-3}$, then the pair cascade would be regulated by triplet pair production and Klein-Nishina scatterings. Even though the mean free path for ICS scatterings in the Klein-Nishina regime is expected to be large, it is still possible to screen the electric field in the gap due to an extra source of pairs provided by the TPP process, an effect which has not been considered in the past for M87.

We next discuss the predictions of our gap model about the gamma-ray emission. To do so, we compute the escaping photon spectra for two indicative cases where the screening of the gap happens either in the low-$E_{0}$ regime (case A) or in the high-$E_{0}$ regime (case B). More specifically, we use the parameters that correspond to the intersection of the grey dashed line with the red-colored and blue-colored curves in the right panel of Figure 4. As discussed previously, the intersection of the curves provides an estimate of the gap size and power, which for cases A and B read respectively $h\sim 0.02\,r_{g}$, $P_{\rm gap}\simeq 3\times 10^{38}$ erg s*-1* and $h\sim 0.2\,r_{g}$, $P_{\rm gap}\simeq 6\times 10^{41}$ erg s*-1*.

The photon spectra escaping from the gap for the two cases are presented in Figure 5 (solid grey lines). When the gap is screened in the low-$E_{0}$ regime (see left panel), the growth of the pair cascade takes place mostly in the Thomson regime, where pairs cool efficiently via ICS. As a result, the escaping gamma-ray luminosity is comparable to the gap power and the photon spectra typically extend from tens to hundreds of TeV depending on the parameters. For case A, in particular, we find $L_{\gamma}\approx P_{\rm gap}/2$ and a peak photon energy of $\sim 10$ TeV (see left panel). On the contrary, if the gap is screened in the high-$E_{0}$ regime, where Klein-Nishina scatterings and TPP are important, the photon spectrum can extend to much higher energies (e.g., $\gtrsim$ PeV), because electron cooling is suppressed (see right panel). This also means that a large fraction of the gap power that is being “pumped” into the pairs by the electric field is not radiated away. In this particular example, we find that $L_{\gamma}\sim 0.01\,P_{\rm gap}$. Thus, in both cases, the direct gamma-ray emission from the gap cannot explain the TeV observations of M87 (see symbols in Figure 5), but for different reasons: either the gap is screened in the Thomson regime and its power is low (left panel) or the gap is screened in the Klein-Nishina-TPP regime and its power is higher, but the radiative efficiency is low.

The gamma-ray spectrum that escapes from the gap is not necessarily the same as the observed one. Photons still have to propagate through the ambient radiation field, and the degree of the gamma-ray flux attenuation will depend on the typical size of the radiation $R_{ph}$, as illustrated in Figure 5 (see dashed colored lines in both panels). Gamma-ray photons with energy $\sim 3/\epsilon_{\min}$ will have the shortest mean free path, and a large dip in the spectrum is expected at this energy, as shown in the right panel of the figure. Photons with either much lower or higher energy will suffer less attenuation.

The absorbed gamma-rays will initiate a post-gap pair cascade, and the absorbed luminosity will re-emerge at lower energies with a spectrum that will depend on the interplay of various cooling processes (e.g., curvature, synchrotron, ICS in the Thomson or Klein-Nishina regimes). Although the attenuated luminosity, for both cases considered here, will be a small fraction of the observed TeV luminosity of M87, the outflowing pairs from the gap will still carry most of the gap power in case B. Thus, if these ultra-relativistic pairs are able to lose their energy in the post-gap region via channels other than ICS and TPP (e.g., curvature radiation), a radiative signal with $L_{\gamma}\sim P_{\rm gap}$ is expected. A detailed study of the post-gap cascade emission will be presented elsewhere.

6 Discussion

We have presented numerical calculations of pair-Compton cascades in the Thomson and Klein-Nishina regimes and investigated the effects of triplet pair production. By adopting a simplified model for the acceleration region and using the method of transport equations for the particle evolution, we computed the growth rate of the pair cascade as a function of the accelerating electric field in the presence of black-body and power-law ambient photon fields. Informed by our numerical results, we derived simple analytical scalings of the growth rate (see Eqs. 16 and 17) and the corresponding electric field (Eq. 2) on the properties of the background photon field. For certain parameters, which can be realized in the vicinity of AGN black holes, the pair cascade may well be regulated by triplet pair production and Compton scatterings in the deep Klein-Nishina regime.

Even if the physical conditions in the acceleration region (i.e., background photon field and electric field strength) are such as to ensure the development of the pair cascade in the deep Klein-Nishina regime, this may not be realized if pairs can lose energy via a process other than TPP and ICS. Because of the ultra-relativistic energies involved in the high-$E_{0}$ regime (see right panel in Figure 1), the most relevant energy loss process becomes curvature radiation. The limiting Lorentz factor in this case is:

[TABLE]

where $\rho_{c}$ is the typical radius of curvature for a particle trajectory, which is assumed to be independent of the particle energy. In the high-$E_{0}$ (TPP-dominated) regime (i.e., for $E_{0}\gtrsim E_{0}^{(\rm KN)}\simeq 5\times 10^{-4}\,U_{ph,0}\,\langle\epsilon\rangle_{0}^{-2}$ statV cm*-1*), the cutoff Lorentz factor of pairs can be approximated by555The logarithmic dependence of the TPP cross section is neglected when deriving this expression.:

[TABLE]

where the notation $\epsilon_{j,x}=\epsilon_{j,x}/(10^{x}\,{\rm eV})$ was introduced, a power-law photon field with $\Gamma=2$ and $\epsilon_{\max}\gg\epsilon_{\min}$ was assumed, and Equation (9) was used to derive the second expression. Curvature losses can be neglected as long as $\gamma_{\rm coff}^{(tpp)}\lesssim\gamma_{\rm cur}$ or equivalently:

[TABLE]

Figure 6 shows the two characteristic Lorentz factors as a function of $E_{0}$ for different values of the minimum energy $\epsilon_{\min}$ of a power-law photon field with index 2. If, for a given $E_{0}$, $\gamma_{\rm coff}^{(tpp)}>\gamma_{\rm cur}$, then the actual growth rate of the pair cascade will be lower than the one obtained in Section 4. For the indicative parameters of the photon field in M87 that we adopted in the previous section (i.e., $U_{ph}=0.1$ erg cm*-3* and $\epsilon_{\min}=10^{-3}$ eV) we find that $\gamma_{\rm coff}^{(tpp)}$ becomes comparable to $\gamma_{\rm cur}$ only for $E_{0}\sim 0.1-1\,\mathrm{statV}\,\mathrm{cm}^{-1}$ (see red line in Figure 6). However, our simplified model predicts that the electric field will not be able to grow to such values and that the gap will be screened by the pair growth in the Thomson regime, as shown in Figure 4 (see red curve in right panel). Nevertheless, for other photon field parameters the pair cascade may develop in the TPP-dominated regime before the particle energy is limited by curvature losses (see e.g., black and blue lines in Figure 6).

In this work, we studied the pair cascade under the assumption of a time-independent electric field. This is a good approximation, if the particle acceleration timescale is much shorter than the timescale for the electric field to evolve. In a dynamic gap, if the electric field grows quadratically with time, as in Equation (18), the condition $t_{\rm acc}(\gamma_{eq})<t_{E}$ in the low-$E_{0}$ regime can be written as:

[TABLE]

where we used Equation (1), $r_{B}=m_{e}c^{2}/(q_{e}B_{0})\approx 17\,B^{-1}_{0,2}$ cm and $\ell_{\rm IC}^{\rm tot}=\langle\epsilon\rangle/\sigma_{T}U_{\rm ph}$ is the mean free path for IC scattering on the overall soft photon energy distribution. For $U_{\rm ph}$ and $B_{0}$ values typical for the black-hole environment of AGN, Equation (24) becomes $E_{0}/B_{0}>10^{-5}M_{9}^{-1}B_{0,2}^{-1/2}U_{\rm ph,0}^{-1/2}$, and can be easily satisfied. In the high-$E_{0}$ (TPP-dominated) regime, the condition on the timescales reads $t_{\rm acc}(\gamma_{\rm coff}^{(tpp)})<t_{E}$ or equivalently:

[TABLE]

where we used Equation (22) and considered a power-law soft photon field with $\Gamma=2$. This condition is more restrictive than the one in Equation (24), since the pairs in the TPP-dominated regime reach ultra-relativistic energies and thereby have much longer acceleration timescales. It is therefore possible that the assumption of constant electric field breaks down in the high-$E_{0}$ regime for certain parameters. In the case where the acceleration time is no more negligible, the growth of the electric field and its response to pair creation should be taken into account. Moreover, our simplifying assumption of identical particles, which disregards the deceleration of positrons, is questionable. In order to simulate the cascade in this case, a proper particle-in-cell (PIC) treatment (Levinson & Cerutti, 2018; Chen et al., 2018) with inclusion of both ICS and TPP processes is necessary. This will be carried out in a dedicated future work.

Here, we have studied the development of pair cascades initiated by $\gamma\gamma$ absorption of $\gamma$-ray photons produced by accelerated electrons in a gap. It is also possible that a non-relativistic proton enters the charge-starved region and accelerates to high energies (Bednarek & Protheroe, 1999). In this case, one expects that protons will start producing pairs via the Bethe-Heitler process ($p\gamma\rightarrow pe^{-}e^{+}$) on ambient soft photons with energy $\epsilon_{0}$ when they reach a Lorentz factor $\gamma^{(pe)}_{p}\gtrsim 2m_{e}c^{2}/\epsilon_{0}$. The pairs will be injected into the gap with relativistic Lorentz factors and a broad energy spectrum (Mastichiadis et al., 2005; Kelner & Aharonian, 2008; Petropoulou & Mastichiadis, 2015). Similar to the pairs produced by $\gamma\gamma$ pair production, the Bethe-Heitler pairs can initiate electromagnetic cascades leading to high multiplicities (i.e., large number of secondary pairs per proton). It is therefore likely that the electric field will be screened before the protons reach the required energy to produce pions and, in turn, neutrinos, as the threshold for photo-meson production ($p\gamma\rightarrow n\pi^{+}(p\pi^{0}$)) is higher, namely $\gamma^{(p\pi)}_{\rm p}\gtrsim m_{\pi}c^{2}/\epsilon_{0}$, where $m_{\pi}$ is the pion mass. However, if the background photon field has a wide energy spectrum (e.g., power law), then even lower energy protons can interact with sufficiently high-energy photons to produce neutrinos. The computation of the escaping neutrino and $\gamma$-ray fluxes in this scenario, which becomes relevant in light of the recent neutrino detections from the $\gamma$-ray AGN TXS 0506+056 (IceCube Collaboration et al., 2018a, b), requires a dedicated study of the ensued pair cascade.

The authors thank the anonymous referee for useful suggestions that helped to improve the manuscript. This research was motivated by discussions at the 2018 MIAPP workshop “The High Energy Universe: Gamma Ray, Neutrino, and Cosmic Ray Astronomy”. It is supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”. M.P. and Y.Y. acknowledge support from the Lyman Spitzer, Jr. Postdoctoral Fellowship. MP acknowledges support also from the Fermi Guest Investigation grant 80NSSC18K1745. A.C. acknowledges support from NASA grant NNX15AM30G.