Revisit of non-linear Landau damping for electrostatic instability driven by blazar-induced pair beams

TL;DR

This paper investigates the impact of non-linear Landau damping on the electrostatic instability of blazar-induced pair beams, revealing how wave energy transfer and IGM temperature influence beam relaxation times.

Contribution

It provides a realistic 2D model for the spectral evolution of electrostatic waves, demonstrating the significance of non-linear Landau damping in pair-beam relaxation.

Findings

Non-linear Landau damping transfers wave energy to small wavenumbers.

Relaxation time depends strongly on IGM temperature, exceeding inverse Compton time at low T_IGM.

Collisions and other non-linear processes can significantly alter relaxation times.

Abstract

We revisit the effect of non-linear Landau (NL) damping on the electrostatic instability of blazar-induced pair beams, using a realistic pair-beam distribution. We employ a simplified 2D model in ${\bf k}$-space to study the evolution of the electric-field spectrum and to calculate the relaxation time of the beam. We demonstrate that the 2D model is an adequate representation of the 3D physics. We find that non-linear Landau damping, once it operates efficiently, transports essentially the entire wave energy to small wavenumbers where wave driving is weak or absent. The relaxation time also strongly depends on the IGM temperature, $T_\mathrm{IGM}$, and for $T_\mathrm{IGM}\ll10$ eV, and in the absence of any other damping mechanism, the relaxation time of the pair beam is longer than the inverse Compton (IC) scattering time. The weak late-time beam energy losses arise from the accumulation of wave energy at small $k$, that non-linearly drains the wave energy at the resonant $\mathbf{k}$ of the pair-beam instability. Any other dissipation process operating at small $k$ would reduce that wave-energy drain and hence lead to stronger pair-beam energy losses. As an example, collisions reduce the relaxation time by an order of magnitude, although their rate is very small. Other non-linear processes, such as the modulation instability, could provide additional damping of the non-resonant waves and dramatically reduce the relaxation time of the pair beam. An accurate description of the spectral evolution of the electrostatic waves is crucial for calculating the relaxation time of the pair beam.

Figures

Tables

Table 1. Table 1: Simulation parameters and relaxation time

Case	$n_e$, cm-3	$T_e$, eV	${\tau }_{rel}$, s
1	${10}^{-7}$	1	$7\cdot {10}^{15}$
2	${10}^{-7}$	5	$7\cdot {10}^{13}$
3	${10}^{-7}$	8	$1.6\cdot {10}^{13}$

Table 2. Table 2: Simulations including particle collisions

Case	$n_e$, cm-3	$T_e$, eV	${\tau }_{rel}$, s
4	${10}^{-7}$	1	$2.7\cdot {10}^{14}$
5	${10}^{-7}$	0.8	$2.2\cdot {10}^{14}$
6	${10}^{-7}$	0.3	$8.3\cdot {10}^{13}$

Equations

$$\begin{gathered} {\omega }_i(k_{\perp },k_{\parallel })=\pi \frac{n_b}{n_e}{\omega }_{p,e}{\left(\frac{{\omega }_{p,e}}{kc}\right)}^3 \\ \times {\int }_{{\theta }_1}^{{\theta }_2}d\theta \frac{-2g\sin \theta +\left(\cos \theta -(kc/{\omega }_{p,e})\cos {\theta }^{\prime }\right)\frac{\partial g}{\partial \theta }}{{\left[(\cos \theta -\cos {\theta }_1)(\cos {\theta }_2-\cos \theta )\right]}^{1/2}}, \end{gathered}$$

$$g(\theta )=m_{e}c{\int }_0^{\infty }pf(p,\theta )dp,$$

$$\cos {\theta }_{1,2}=\frac{{\omega }_{p,e}}{kc}\left[\cos {\theta }^{\prime }\pm \sin {\theta }^{\prime }{\left({\left(\frac{kc}{{\omega }_{p,e}}\right)}^2-1\right)}^{1/2}\right].$$

$$\frac{dW(\mathbf{k})}{dt}=2({\omega }_i(\mathbf{k})+{\omega }_{LL}(\mathbf{k})+{\omega }_{NL}(\mathbf{k})+{\omega }_c)W(\mathbf{k}),$$

$${\omega }_{LL}(\mathbf{k})=-{\omega }_{p,e}\sqrt{\frac{\pi }{8}}{\left(\frac{{\omega }_{p,e}}{k{u}_e}\right)}^3\exp \left[-\frac{1}{2}{\left(\frac{{\omega }_{p,e}}{k{u}_e}\right)}^2\right],$$

$$\begin{gathered} {\omega }_{NL}=\frac{3(2\pi {)}^{1/2}}{64n_e{m}_e{u}_i}\int d^3{k}^{\prime }W({\mathbf{k}}^{\prime })\frac{(\mathbf{k}{\mathbf{k}}^{\prime }{)}^2}{(k^{\prime }k{)}^2}\frac{k^{\prime 2}-k^2}{|{\mathbf{k}}^{\prime }-\mathbf{k}|} \\ \times \exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k^{\prime 2}-k^2}{|{\mathbf{k}}^{\prime }-\mathbf{k}|}\right)}^2\right], \end{gathered}$$

$${\omega }_c=-1.45\cdot 10^{-6}{n}_e\lambda T_e^{-3/2}\text{ [s−1],}$$

$$W_{tot}=2\pi \int W(\mathbf{k})k_{\perp }d{k}_{\perp }d{k}_{\parallel }.$$

$${\omega }_i(k_{\perp })=\frac{{\omega }_i(k_{\perp }c/{\omega }_{p,e}=1)}{1+(k_{\perp }c/{\omega }_{p,e}{)}^{1.7}}$$

$$W(\mathbf{k},t=0)=W_0\approx k_B{T}_{\mathrm{IGM}}\approx 1.6\cdot 10^{-12}{T}_4,\text{ erg},$$

$$\begin{gathered} {\omega }_i(k_{\perp })=\frac{{\omega }_i(k_{\perp }c/{\omega }_{p,e}=1)\theta (k_{\perp }c/{\omega }_{p,e}-1)}{1+(k_{\perp }c/{\omega }_{p,e}{)}^{1.7}}\approx \\ 5.6\cdot 10^{-6}\frac{n_{b20}}{\sqrt{n_{e7}}}\frac{\theta (k_{\perp }c/{\omega }_{p,e}-1)}{1+(k_{\perp }c/{\omega }_{p,e}{)}^{1.7}},\text{ s−1}, \end{gathered}$$

$${\omega }_{LL}=-{\omega }_{p,e}\sqrt{\frac{\pi }{8}}{\left(\frac{{\omega }_{p,e}}{k_{\perp }{u}_e}\right)}^3\exp \left[-\frac{1}{2}{\left(\frac{{\omega }_{p,e}}{k_{\perp }{u}_e}\right)}^2\right],$$

$$\begin{gathered} {\omega }_{NL}=\frac{3(2\pi {)}^{1/2}}{64n_e{m}_e{u}_i}\int d^2{k}_{\perp }^{\prime }W(k_{\perp }^{\prime },t)\frac{({\mathbf{k}}_{\perp }{{\mathbf{k}}_{\perp }}^{\prime }{)}^2}{(k_{\perp }^{\prime }{k}_{\perp }{)}^2}\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{|{\mathbf{k}}_{\perp }^{\prime }-{\mathbf{k}}_{\perp }|} \\ \times \exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{|{\mathbf{k}}_{\perp }^{\prime }-{\mathbf{k}}_{\perp }|}\right)}^2\right], \end{gathered}$$

$$W_{tot}(t)=2\pi \int W(k_{\perp },t)k_{\perp }d{k}_{\perp }.$$

$$W(k_{\perp },t=0)=W_{0,\perp }\approx k_B{T}_{\mathrm{IGM}}\Delta k_{\parallel },$$

$${\int }_0^{\tau (\alpha )}P(t)dt=\alpha W_b(t=0)$$

$$P(t)=8\pi \int {\omega }_{i}W(k_{\perp },t)k_{\perp }d{k}_{\perp }$$

$${\int }_0^{\tau (\alpha )}P(t)dt\approx \frac{\tau }{T}{\int }_{t_{NL}}^{t_{NL}+T}P(t)dt,$$

$${\tau }_{rel}=\frac{T{W}_b(t=0)}{2}{\left({\int }_{t_{NL}}^{t_{NL}+T}P(t)dt\right)}^{-1}.$$

$${\omega }_{NL}(k_{\parallel },k_{\perp })=\frac{3(2\pi {)}^{1/2}}{64n_e{m}_e{u}_i}\int k_{\perp }^{\prime }d{k}_{\perp }^{\prime }d{\varphi }^{\prime }d{k}_{\parallel }^{\prime }IW(k_{\parallel }^{\prime },k_{\perp }^{\prime })\frac{k^{\prime 2}-k^2}{(k^{\prime }k{)}^2},$$

$$I=\frac{(k_{\parallel }^{\prime }{k}_{\parallel }+k_{\perp }^{\prime }{k}_{\perp }\cos {\varphi }^{\prime }{)}^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+k_{\perp }^{\prime 2}+k_{\perp }^2-2k_{\perp }^{\prime }{k}_{\perp }\cos {\varphi }^{\prime }\right]}^{1/2}}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k^{\prime 2}-k^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+k_{\perp }^{\prime 2}+k_{\perp }^2-2k_{\perp }^{\prime }{k}_{\perp }\cos {\varphi }^{\prime }\right]}^{1/2}}\right)}^2\right].$$

$$k_{\perp }^{\prime }\gg k_{\perp }\vee k_{\perp }^{\prime }\ll k_{\perp }\vee |k_{\parallel }^{\prime }-k_{\parallel }|\gg \sqrt{2k_{\perp }^{\prime }{k}_{\perp }}$$

$${\int }_0^{2\pi }Id{\varphi }^{\prime }=\pi \frac{2(k_{\parallel }^{\prime }{k}_{\parallel }{)}^2+(k_{\perp }^{\prime }{k}_{\perp }{)}^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+k_{\perp }^{\prime 2}+k_{\perp }^2\right]}^{1/2}}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k^{\prime 2}-k^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+k_{\perp }^{\prime 2}+k_{\perp }^2\right]}^{1/2}}\right)}^2\right],$$

$${\int }_0^{2\pi }Id{\varphi }^{\prime }\approx {\int }_0^{2\pi }d{\varphi }^{\prime }\frac{(k_{\parallel }^{\prime }{k}_{\parallel }+k_{\perp }^2\cos {\varphi }^{\prime }{)}^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+2k_{\perp }^2(1-\cos {\varphi }^{\prime })\right]}^{1/2}}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k_{\parallel }^{\prime 2}-k_{\parallel }^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+2k_{\perp }^2(1-\cos {\varphi }^{\prime })\right]}^{1/2}}\right)}^2\right].$$

$$\begin{gathered} {\int }_{0.3\pi }^{1.7\pi }Id{\varphi }^{\prime }\sim {\int }_0^{2\pi }d{\varphi }^{\prime }\frac{(k_{\parallel }^{\prime }{k}_{\parallel }+k_{\perp }^2\cos {\varphi }^{\prime }{)}^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+2k_{\perp }^2\right]}^{1/2}}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k_{\parallel }^{\prime 2}-k_{\parallel }^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+2k_{\perp }^2\right]}^{1/2}}\right)}^2\right]= \\ \pi \frac{2(k_{\parallel }^{\prime }{k}_{\parallel }{)}^2+(k_{\perp }{)}^4}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+2k_{\perp }^2\right]}^{1/2}}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k_{\parallel }^{\prime 2}-k_{\parallel }^2}{{\left[(k_{\parallel }^{\prime }-k_{\parallel }{)}^2+2k_{\perp }^2\right]}^{1/2}}\right)}^2\right]. \end{gathered}$$

$${\int }_{2\pi -\Delta }^{0+\Delta }Id{\varphi }^{\prime }\approx \Delta \frac{(k_{\parallel }^{\prime }{k}_{\parallel }+k_{\perp }^2{)}^2}{|k_{\parallel }^{\prime }-k_{\parallel }|}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k_{\parallel }^{\prime 2}-k_{\parallel }^2}{|k_{\parallel }^{\prime }-k_{\parallel }|}\right)}^2\right]$$

$${\omega }_{NL}(k_{\perp })=\frac{3(2\pi {)}^{1/2}}{64n_e{m}_e{u}_i}\int k_{\perp }^{\prime }d{k}_{\perp }^{\prime }d{\varphi }^{\prime }GW(k_{\perp }^{\prime })\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{(k_{\perp }^{\prime }{k}_{\perp }{)}^2},$$

$$G=\frac{(k_{\perp }^{\prime }{k}_{\perp }\cos {\varphi }^{\prime }{)}^2}{{\left[k_{\perp }^{\prime 2}+k_{\perp }^2-2k_{\perp }^{\prime }{k}_{\perp }\cos {\varphi }^{\prime }\right]}^{1/2}}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{{\left[k_{\perp }^{\prime 2}+k_{\perp }^2-2k_{\perp }^{\prime }{k}_{\perp }\cos {\varphi }^{\prime }\right]}^{1/2}}\right)}^2\right].$$

$${\int }_0^{2\pi }Gd{\varphi }^{\prime }\approx \pi \frac{(k_{\perp }^{\prime }{k}_{\perp }{)}^2}{{\left(k_{\perp }^{\prime 2}+k_{\perp }^2\right)}^{1/2}}\exp \left[-a{\left(\frac{c}{{\omega }_{p,e}}\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{{\left(k_{\perp }^{\prime 2}+k_{\perp }^2\right)}^{1/2}}\right)}^2\right].$$

$$\begin{gathered} {\omega }_{NL}^{3D}(k\approx {\omega }_{p,e}/c)\approx \frac{(2\pi {)}^{1/2}}{64n_e{m}_e{u}_i}\int d^2{k}_{\perp }^{\prime }{W}_{3D}({\mathbf{k}}^{\prime })\Delta k_{\parallel }\left(\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{|{\mathbf{k}}_{\perp }^{\prime }-{\mathbf{k}}_{\perp }|}+\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{|{\mathbf{k}}_{\perp }^{\prime }-{\mathbf{k}}_{\perp }-2{\mathbf{k}}_{\parallel }|}\right)= \\ \frac{{\pi }^{1/2}}{64n_e{m}_e{u}_i}\int d^2{k}_{\perp }^{\prime }{W}_{3D}({\mathbf{k}}^{\prime })\Delta k_{\parallel }\frac{k_{\perp }^{\prime 2}-k_{\perp }^2}{|{\mathbf{k}}_{\perp }^{\prime }-{\mathbf{k}}_{\perp }|}\left(1+\frac{|{\mathbf{k}}_{\perp }^{\prime }-{\mathbf{k}}_{\perp }|}{|{\mathbf{k}}_{\perp }^{\prime }-{\mathbf{k}}_{\perp }-2{\mathbf{k}}_{\parallel }|}\right), \end{gathered}$$

Full text

gamma rays, magnetic fields, instabilities, waves, relativistic processes

1 Introduction

The propagation of very high energy gamma-radiation ($E_{\gamma}>100$ GeV) and its reaction on the inter-galactic medium (IGM) have been actively studied in recent years both observationally (Neronov & Semikoz, 2009; Neronov & Vovk, 2010; Ackermann et al., 2012; Frank et al., 2013; de Naurois, 2015; Funk, 2015) and theoretically (Elyiv et al., 2009; Taylor et al., 2011; Broderick et al., 2012; Puchwein et al., 2012; Schlickeiser et al., 2012a, b; Miniati & Elyiv, 2013; Sironi & Giannios, 2014; Broderick et al., 2016; Rafighi et al., 2017; Vafin et al., 2018). An important feature of high-energy photons is the ability to interact with the extra-galactic background light (EBL) producing ultra-relativistic electron-positron beams (Gould & Schréder, 1966). Due to the inverse Compton (IC) scattering, these pairs will emit secondary photons (Aharonian et al., 1994) that can be also energetic enough to efficiently interact with the EBL. The observational analysis indicates indicates, however, that the measured gamma-ray signal in the GeV energy band is smaller than the predicted cascade emission assuming that the pairs lose their energy only due to the IC scattering (Neronov & Vovk, 2010; Tavecchio et al., 2010, 2011). Thus, some other dissipation processes must be in play.

One possible explanation relies on the existence of a large-scale hypothetical magnetic field in cosmic voids (Elyiv et al., 2009; Neronov & Vovk, 2010; Taylor et al., 2011) which remains questionable, as no direct magnetic field observations are yet available, but would, if true, havle profound consequences for magnetogenesis in the universe. The predicted magnetic field should result in GeV halos (about TeV sources) that are, however, absent as was shown recently (Broderick et al., 2018).

An alternative model uses only the fact that an electrositron beam propagating through the IGM plasma is subject to the electrostatic (two-stream) instability resulting in beam energy dissipation into the plasma waves (Broderick et al., 2012; Schlickeiser et al., 2012b; Vafin et al., 2018; Shalaby et al., 2018). Whether or not that is faster than IC cooling depends on the ratio of the beam-relaxation time and the IC scattering time. The former is determined by the instability growth rate as well as non-linear processes such as the non-linear Landau (NL) damping (Breizman et al., 1972; Breizman, 1990) and the modulation instability (Zakharov, 1972; Papadopoulos, 1975; Baikov, 1977). Using a simplified model, Miniati & Elyiv (2013) claimed that the NL damping can stabilize the realistic blazar-induced pair beam on the IC timescale making the effect of plasma instabilities negligible. Later, Chang et al. (2014) rigorously solving the kinetic wave equation came to the opposite conclusion. Their result relies on the fact that there is a region at large wave numbers where the electrostatic instability dominates over the NL as well as LL (linear Landau) damping. In that region, the electrostatic instability can efficiently reduce the beam energy. However, the work by Chang et al. (2014) seems to overestimate the growth rate assuming it to be constant at wave numbers $|{\bf k}|\geq\omega_{p,e}/c$ ($\omega_{p,e}=\sqrt{4\pi n_{e}e^{2}/m_{e}}$ is the electron plasma frequency, $n_{e}$ the electron number density, and $c$ the speed of light). This fact contradicts the realistic growth rate which rapidly decreases with the wave number (Miniati & Elyiv, 2013; Vafin et al., 2018). Therefore, the primary goal of this paper is to investigate the effect of the NL damping taking into account the realistic behavior of the electrostatic growth rate at large wave numbers for blazar-induced pair beams.

Other dissipation processes, such as the modulation instability and collisional damping, can also have a dramatic effect on the relaxation time of the pair beams. In the current work, we focus on the collisional damping in addition to the NL damping, whereas the modulation instability requires a more careful study going beyond the scope of the current manuscript.

The rest of the paper is organized as follows. Section 2 explains our theoretical model. Simulation results without and with the collisional damping are presented in section 3 and 4, respectively. The final summary of the work is given in section 5.

2 Theoretical model

2.1 Electrostatic growth rate

The growth rate of the electrostatic instability for an ultra-relativistic beam with the homogeneous number density $n_{b}$ and the axis-symmetric momentum distribution $f_{b}(p,\theta)$ in a spherical coordinate system reads

[TABLE]

where

[TABLE]

[TABLE]

Here, $\theta^{\prime}$ is the angle between the wave vector ${\bf k}$ and the beam direction, $c$ the speed of light, $n_{b}$ the beam density, $n_{e}$ the electron density in the IGM. Using the approximation for the pair-beam distribution at a distance 50 Mpc from a fiducial blazar obtained in our previous work (Eqs. (24)-(25), and (56) from Vafin et al. (2018) ), we calculated the electrostatic growth rate Eq. (1) shown in Fig. 1. It is well seen that the growth rate decreases towards large $k_{\perp}$ in contrast to the growth rate used by Chang et al. (2014). Thus, our model provides a more realistic description of the instability growth.

2.2 Three-dimensional (3D) evolution of the wave spectrum

Then, the non-linear wave kinetic equation including the LL, NL and collisional damping reads

[TABLE]

where $W({\bf k})$ is the spectral energy density of the electric field and

[TABLE]

[TABLE]

are, respectively, the LL and NL damping rates (Breizman et al., 1972). Here, $a=(9/8)[u_{e}^{2}/(cu_{i})]^{2}$ where $u_{i,e}=\sqrt{T_{\mathrm{IGM}}/m_{i,e}}$ denotes the IGM ion and electron thermal speeds and $T_{\mathrm{IGM}}$ is the IGM temperature. For $T_{\mathrm{IGM}}=10^{4}T_{4}\ \mathrm{K}$: $a\approx 3.6\cdot 10^{-3}T_{4}$. The collisional damping rate reads (Alexandrov et al., 1984; Miniati & Elyiv, 2013; Huba, 2016):

[TABLE]

where $\lambda=23.5-\ln(n_{e}^{1/2}T_{\mathrm{IGM}}^{-5/4})-[10^{-5}+(\ln T_{\mathrm{IGM}}-2)^{2}/16]^{1/2}$. Here, $n_{e}$ and $T_{\mathrm{IGM}}$ are in units cm*-3* and eV, respectively. Moreover, we have implicitly assumed that the real part of the frequency $\omega$ is approximately equal to the the plasma frequency $\omega_{p,e}$ in Eqs. (1) and (5). It is valid for the wave vectors $kc/\omega_{p,e}\ll c/u_{e}\approx 0.7\cdot 10^{3}/\sqrt{T_{4}}$.

The total electric field energy density is calculated as

[TABLE]

It is to be noted that the NL damping does not affect the total energy density of the electric field. It can be easily seen from the integration of Eq. (4) over all ${\bf k}$ that the term with the NL damping $\int\omega_{NL}({\bf k})W({\bf k})dk=0$ actually disappears. Thus, the process, traditionally called NL damping, does not actually damp waves, but only scatters them from larger to smaller wave numbers.

The spectrum evolution qualitatively consists of several stages. First, the wave energy grows exponentially during the linear instability phase and NL damping can be neglected. Afterward, NL damping starts efficiently scattering the waves towards smaller wave numbers. Then a rapid transition happens when most of the wave energy is transferred out of the instability region in ${\bf k}$-space, and as a consequence, the wave energy approaches a quasi-stationary level. However, the electric field energy continues to increase (but now much slower as compared to the linear phase), since there exists such a region in ${\bf k}$-space where both NL and LL damping are suppressed (Chang et al., 2014). In this region, the electrostatic waves are still unstable and can further reduce the beam energy.

Due to the cylindrical symmetry of the problem, the solution $W({\bf k}^{\prime})$ depends only on $k_{\parallel}^{\prime}$ and $k_{\perp}^{\prime}$. Then the integral with respect to the angle between ${\bf k}^{\prime}$ and ${\bf k}$ in Eq. (6) can be performed before solving Eq. (4). This integral is approximately evaluated in Appendix A.

As it is seen from Fig. 1, the electrostatic instability at parallel wave vectors develops in a very narrow band around $k_{\parallel}c/\omega_{p,e}\approx 1$. Therefore, in the parallel direction, we may limit the problem only to the region $k_{\parallel}c/\omega_{p,e}\lesssim 1$.

In the perpendicular direction to the beam, the growth rate extends to much larger wave numbers, and $k_{\perp}$ should be considered up to the value where $\omega_{i}\approx\omega_{LL}$. Fig. 2 demonstrates the growth rate as a function of $k_{\perp}$, while for each $k_{\perp}$ we took the maximum growth rate in the parallel direction. The growth rate maximized in the parallel direction can be well approximated by:

[TABLE]

which is shown by the thick line in Fig. 2. For the pair-beam number density $n_{b}=3\cdot 10^{-22}$ cm*-3* found in our previous work and corresponding to the growth rate presented in Fig. 1 and for the typical IGM parameters $n_{e}=10^{-7}$ cm*-3*, $T_{\mathrm{IGM}}=10^{4}$ K, Fig. 3 compares the growth rate Eq. (9) with the LL and collisional damping rate, Eq. (5) and Eq. 7. Since $k_{\parallel}c/\omega_{p,e}\lesssim 1$ and the LL rate is significant at $kc/\omega_{p,e}\gg 1$, we may set $k\approx k_{\perp}$ in Eq. (5). Fig. 3 shows that the LL damping suppresses the electrostatic instability at $k_{\perp}c/\omega_{p,e}\gtrsim 90$, so we may take into account only $k_{\perp}c/\omega_{p,e}\lesssim 100$. For higher temperatures, say $T_{\mathrm{IGM}}\approx$ $6\cdot 10^{4}$ K and $9\cdot 10^{4}$ K, the LL damping suppresses the electrostatic instability at $k_{\perp}c/\omega_{p,e}\gtrsim 45$ and $k_{\perp}c/\omega_{p,e}\gtrsim 35$, respectively. The range of $k_{\perp}c/\omega_{p,e}$, for our simulations, is chosen accordingly.

As the initial condition for Eq. (4), the constant value given by the discrete particle fluctuations can be used (Klimontovich, 1982; Schlickeiser et al., 2012b):

[TABLE]

where $T_{\mathrm{IGM}}=10^{4}T_{4}$ K. Eq. (10) provides the noise level of electric field turbulence during the time evolution. As there is no term responsible for the noise generation in Eq. (4), the solution $W({\bf k})$ must be forced to stay above $W_{0}$ by the condition $W({\bf k})\geq W_{0}$ in couple with Eq. (4). However, this provides an artificial energy source, but it contributes a much smaller energy fraction compared to the effect of electrostatic instability.

2.3 Reduction to two-dimensional (2D) model

To investigate Eq. (4), we utilize a 2D model ($W=W(k_{\perp})$) for analyzing the problem in the plane $k_{\parallel}c/\omega_{p,e}\approx 1$ of the ${\bf k}$-space. To preserve the fact that the waves are scattered out to the non-resonant wave vectors, i.e., out of the instability region, we modify the growth rate to

[TABLE]

where $\theta$ is the Heaviside step function. The growth rate is truncated at $k_{\perp}c/\omega_{p,e}=1$, for waves at small $k_{\perp}$ to be stable. At the same time, Eqs. (5)-(6) and Eq. (8) can be reduced to

[TABLE]

[TABLE]

[TABLE]

The angular integral in Eq. (13) can be evaluated similar to Appendix A with the result presented in Appendix B.

The initial condition for Eq. (4) in 2D can be taken as

[TABLE]

where $\Delta k_{\parallel}\approx 10^{-6}\omega_{p,e}/c$ is approximately the width of the instability region in the parallel direction (see Fig. 1). Eq. (15) can be a rather rough estimate, but the exact estimation of the noise energy density is unimportant for the later evolution as pointed out by Chang et al. (2014). Appendix C demonstrates analytically that the 2D model reflects the main features of the original 3D one and can be used to study the spectral equation (4). The results of a 3D test simulation are also discussed in Appendix D in comparison to those of our 2D model.

3 Numerical results without collisions

We first analyze the system of Eqs. (4) and (11)-(13) numerically using the finite-difference method and neglecting particle collisions. The pair beam density and average gamma-factor are, respectively, $n_{b}=3\cdot 10^{-22}$ cm*-3* and $<\gamma_{b}>=4\cdot 10^{6}$ that correspond to the pairs produced within the IC mean free path at distance $50$ Mpc from a fiducial source of high-energy gamma-rays (Vafin et al., 2018). For this beam setup, we then consider the effect of the IGM temperature on the energy loss time of the beam $\tau(\alpha)$ defined as

[TABLE]

where

[TABLE]

is the power loss of the beam, $W_{b}(t=0)=<\gamma_{b}>n_{b}m_{e}c^{2}$ the initial beam energy density, and $\alpha$ the fraction of initial beam energy lost at time $\tau(\alpha)$. Here, we account for the equipartition of beam energy into kinetic energy and electrostatic energy of the plasma waves. In the case of a cold beam, the fraction $\alpha$ can be as high as 0.5 (Schlickeiser et al., 2002) which we will also use for our estimation of the relaxation time of the beam $\tau_{rel}=\tau(0.5)$.

In principle, the relaxation time can be estimated from the simulation. However, if $\tau_{rel}$ is large, it can take a prohibitively long time to achieve it numerically, even in 2D. A more efficient way of calculating $\tau_{rel}$ is to extrapolate simulation results in the following manner. It will be demonstrated below that during the quasi-saturation stage the power loss function $P(t)$ changes almost periodically. Since the system turns to a quasi-saturation at a time $t_{NL}$ much shorter than $\tau_{rel}$ and the beam energy loss at $t=t_{NL}$ is only about 1%, we can evaluate the integral Eq. (16) approximately as

[TABLE]

where $T$ is the oscillation period of $P(t)$ at $t>t_{NL}$. Then

[TABLE]

The simulation parameters are summarized in Table 1. Case 1 treats the beam evolution in the IGM with the average parameters, while cases 2 and 3 consider higher temperatures of the IGM. The energy loss time becomes shorter for a higher IGM temperature, because of smaller NL damping. The latter has to do with the exponential cut-off in integrand of $\omega_{NL}$ where the factor $a\propto T_{e}$. Thus, the NL damping is weaker for higher plasma temperatures. Cases 2 and 3 with higher IGM temperatures refer to the IGM heating scenario by blazar-induced pair beams proposed by Puchwein et al. (2012). Therefore, it is crucial to understand qualitatively how the IGM temperature effects the relaxation time of the pair beam.

The highest value of $k_{\perp}$ in each simulation is determined by equality of the LL damping rate and the linear growth rate (see Fig. 2). The lowest value of $k_{\perp}$ was chosen small enough ($10^{-8}\omega_{p,e}/c$) to prevent peaking of the solution at the single left grid point within the simulated time interval. The number of grid points was adjusted automatically at each time step to maintain the accuracy of the solution within 10 %.

Fig. 4 illustrates the spectrum evolution for case 1. During early times $t\lesssim 10^{8}$ s, there is one pronounced maximum in the spectrum (black curve in Fig. 4) resulting from the modes with the highest growth rate. Afterward, NL damping becomes important and efficiently scatters waves to lower $k_{\perp}$. As a result, the peak at $k_{\perp}c/\omega_{p,e}\sim 1$ disappears, and the spectral density in the region of the highest growth rate $W(k_{\perp}c/\omega_{p,e}\sim 1)$ becomes equal to the noise level $W_{0,\perp}$. Any subsequent growth at $k_{\perp}c/\omega_{p,e}\sim 1$ is efficiently suppressed by NL damping, and W($k_{\perp}c/\omega_{p,e}\sim 1$) remains at the noise level during the remainder of the simulation time. In this case, the NL damping is determined by the non-resonant waves at $k_{\perp}c/\omega_{p,e}<1$.

Therefore, as mentioned above, the long-term beam evolution is determined by the spectral properties at large $k_{\perp}$, where both the LL and NL damping rates can be smaller than the growth rate. Fig. 5 demonstrates the spectrum at large $k_{\perp}$ for much longer simulation times compared to Fig. 4. The spectrum at large $k_{\perp}$ evolves somewhat periodically: the waves continuously grow and get scattered to smaller wave numbers. The right maximum in Fig. 5 corresponds to the region where the waves are amplified and the left one arises due to the NL damping of the waves in the right peak. NL damping transfers this energy to the low-$k$ region at later times. However, the saturation energy density at low-$k_{\perp}$ remains effectively constant since the energy density at large $k_{\perp}$ is much smaller compared to the energy of non-resonant waves at $k_{\perp}c/\omega_{p,e}<1$. The spectrum at $k_{\perp}c/\omega_{p,e}<1$ stays almost constant, and the spectrum exhibits a flat peak at $k_{\perp}\simeq 0$.

The evolution of the total energy density of the electric field and the power losses of the beam for case 1 are shown in Figs. 6 and 7, respectively. It is well seen from Fig. 6 that the beam loses only a small fraction of its energy ($\sim 10^{-2}W_{b}$) during the linear instability phase $t\lesssim 3\cdot 10^{8}$ s. At the same time, $P(t)$ achieves its highest value at $t\approx 3.3\cdot 10^{8}$ s. However, almost immediately afterward the power losses undergo a rapid decline due to the fact that most of the waves are scattered to the region $k_{\perp}c/\omega_{p,e}<1$ where the instability growth rate is zero, $\omega_{i}=0$, in our model. The wave growth at large $k_{\perp}$ results in the increase of $P(t)$ again at $t>10^{10}$ s. Eventually, $P(t)$ turns to a periodic solution.

Fig. 8 compares the power losses of the beam for cases 1-3. Here, it is well seen that the function $P(t)$ changes periodically during the quasi-saturation phase and, therefore, our extrapolation of the time integral given by Eq. (18) is justified. However, since the oscillation amplitude of $P(t)$ changes with time, we considered a time interval $T$ where $P(t)$ achieves the maximum multiple times to estimate Eq. (18). We also note that the power losses are strongly affected by the plasma temperature. The power losses grow with the increase in plasma temperature because the NL damping rate becomes suppressed due to a stronger exponential cutoff in the NL damping rate as evident from Fig. 8. The relaxation time Eq. (19) for cases 1-3 is summarized in Table 1. The relaxation time should be compared to the IC scattering time $\tau_{IC}\approx 10^{20}/\gamma_{b}=2.5\cdot 10^{13}$ s for $\gamma_{b}=4\cdot 10^{6}$. We conclude that only for a hot enough IGM the relaxation time can become slightly smaller than the IC time.

Thus, if the NL damping is the only non-linear stabilization process, the electrostatic instability has a negligible effect on the beam energy dissipation within the IC time scale, $\tau_{IC}$, and other effects should be invoked to explain the lack of the cascade GeV signal (Neronov & Vovk, 2010). The reason for such a strong suppression of the instability is that the NL damping rate in the region with the highest growth rate ($k_{\perp}c/\omega_{p,e}\simeq 1$) is proportional to the energy density of the long-wave non-resonant oscillations at $k_{\perp}c/\omega_{p,e}\ll 1$. Obviously, an effective dissipation of energy in the long-wave part of the spectrum will result in a smaller NL damping rate and, consequently, in a shorter energy loss time of the pair beam modifying the conclusion of our simulations above. Possible dissipation mechanisms are the modulation instability and collisional damping. While the former requires a thorough investigation going well beyond the scope of the current work, the collisional rate can be easily estimated for electrons with $n_{e}=10^{-7}$ cm*-3* and $T_{e}=1$ eV to be $\nu_{c}\approx 10^{-11}$ s*-1*. Since the IC time for the beam considered above is $\tau_{IC}\approx 2.5\cdot 10^{13}$ s, the collisional dissipation can significantly affect the value of $\tau_{loss}$ in case 1. It is to be noted that for case 3, where $\tau_{loss}<\tau_{IC}$ s, $\nu_{c}\approx 2.2\cdot 10^{-13}$ s and collisions do not play a considerable role at $t\lesssim\tau_{IC}$.

4 Effect of particle collisions

Presence of particle collisions results in an additional damping rate in Eq. (4). In view of the fact that the effect of particle collisions becomes important for low IGM temperatures, we carried out three simulations with parameters given by Table 2.

Let us consider case 4 which treats the same IGM parameters as case 1 but now including the effect of particle collisions. In this case, collisions are negligible up to the time $\approx 10^{10}$ s, but afterward, they considerably dissipate the electric field energy by almost two orders of magnitude which is illustrated in Fig. 9. If we assume for a rough estimation that collisions dissipate the energy exponentially $\propto\exp(-2|\omega_{c}|t)$, then the energy will be reduced by two orders of magnitude after the time $\Delta t=\ln(10)/|\omega_{c}|\approx 2.3/|\omega_{c}|\approx 5.1\cdot 10^{11}$ s which perfectly agrees with our numerical result. Collisions reduce the energy density of the non-resonant waves at small $k_{\perp}c/\omega_{p,e}<1$, as the collision frequency can be neglected compared to the linear growth rate at $k_{\perp}c/\omega_{p,e}>1$ (see Fig. 3). This energy dissipation modifies the NL damping rate which is now also reduced by almost two orders of magnitude (see Fig. 10). As a result, the NL damping rate can become smaller than the maximum instability growth rate near $k_{\perp}c/\omega_{p,e}=1$ as it is illustrated in the zoomed view in Fig. 10. It naturally leads to an electric field growth at $k_{\perp}c/\omega_{p,e}=1$ at later times in addition to the electrostatic waves at large $k_{\perp}$. The corresponding spectrum of these waves is shown in Fig. 11 for several instants of time (we do not show the region $k_{\perp}c/\omega_{p,e}\gg 1$ because it remains similar to case 1, Fig. 5). It is to be emphasized that a similar peak at later times did not appear in case 1, as the NL damping rate at $k_{\perp}c/\omega_{p,e}=1$ for case 1 stays above the growth rate during the non-linear stage. The spectral peak in Fig. 11 periodically grows and drops with time which drastically modifies the power losses of the beam. Fig. 12 shown the time evolution of the power losses for case 4 at later times when the effect of collisions becomes noticeable. In Fig. 12, we can recognize two periodic structures: sharp large amplitude and broad small amplitude peaks. The latter is similar to case 1 and is caused by the electric field at $k_{\perp}c/\omega_{p,e}\gg 1$, while the former results is an effect of collisions and is related to the spectral peak at $k_{\perp}c/\omega_{p,e}=1$ (Fig. 11). Due to collisions, the power losses achieve a much higher value compared to case 1. This reduces the relaxation time by factor 20 and it becomes $2.7\cdot 10^{14}$ s. Thus, the beam losses its energy much faster (Fig. 13). Table 2 compares the relaxation time for smaller IGM temperatures. From what follows that the relaxation time is well above the IC time for small IGM temperatures even if collisions are present.

5 Summary and conclusions

We revisited the effect of NL damping on the relaxation time of the blazar-induced pair beam, $\tau_{rel}$. The relaxation time of the pair beams is given by the wave growth rate and the wave intensity in the resonant band of the spectrum, and so one needs to know the saturation level of the wave intensity and its stability, as the relaxation time is far longer than the time needed to reach saturation. NL damping is an important factor in the nonlinear evolution of the intensity of the electrostatic plasma waves that the pair beams drive. During the quasi-stationary phase, its efficiency primarily depends on the wave intensity at small $k$ (where it is much higher than anywhere else in the spectrum), which it provides itself by transferring wave energy from the resonant band of the spectrum, implying a high level of nonlinearity. We investigate this problem numerically using a simplified 2D model which we demonstrate to be an acceptable approximation of the real 3D physics. We also consider a more realistic growth rate of the electrostatic instability than is found in, e.g., Chang et al. (2014).

In an earlier study (Vafin et al., 2018) we used a simplified description of NL damping and found it weaker than the modulation instability and both permitting a saturation intensity, for which the relaxation time of the pair beam is shorter than the IC time ($\tau_{IC}=2.5\cdot 10^{13}$ s). This new study is based on a detailed calculation of the intensity spectrum and its evolution, which should provide an accurate estimate of the relaxation time of the beam.

Puchwein et al. (2012) argued that the dissipation of pair-beam energy may be an efficient means of IGM heating. Since the NL damping rate exponentially depends on the plasma temperature, we performed several simulations to explore the effect of the IGM temperature on $\tau_{rel}$. At the average IGM temperature of 1 eV, the relaxation time of the pair beam is much longer than at $T_{\mathrm{IGM}}=8$ eV.

We find that quasi-perpendicular waves ($k_{\parallel}c/\omega_{p,e}\simeq 1$ and $k_{\perp}c/\omega_{p,e}\gg 1$) display a complicated and time-variable spectrum whose intensity is typically too low to impose significant beam dissipation. In the spectral band where the growth rate is highest ($k_{\parallel}c/\omega_{p,e}\simeq 1$ and $k_{\perp}c/\omega_{p,e}\approx 1$), the wave intensity is completely drained by NL damping driven by very high wave intensity at small $k$. Taken at face value, these results suggest that the beam dissipates a negligible portion of its energy during the IC cooling time. However, if other dissipation mechanisms are present that can efficiently dissipate the non-resonant waves at small $k$ that are responsible for a high NL damping rate, then the relaxation time of the beam would be significantly reduced. Among the possible candidates are the modulation instability and the collisional damping. In the current work, we considered only the latter.

Since the effect of collisions becomes important at small plasma temperatures, we performed several calculations with $T_{\mathrm{IGM}}<1$ eV including collisional damping. We found that the relaxation time is indeed considerably shorter than without collisions. This is surprising in view of the collision rate being about four orders of magnitude smaller than the peak growth rate of the electrostatic waves, and it reflects the delayed build-up of wave intensity at small $k$ by NL damping. Formally, even for $T_{\mathrm{IGM}}=0.3$ eV, the relaxation time is still longer than the IC time. The large impact of a seemingly sub-dominant process indicates that the beam dissipation rate is extremely sensitive to the dissipation processes under consideration and to the accuracy of their treatment. The latter is best seen in the difference in estimated relaxation time between approximating NL damping by a simple decay term (as in Vafin et al. (2018)) and explicitly following the evolution of the wave spectrum, as done here. Thus, any other damping process (e.g., the modulation instability) will further modify the dissipation rate, and a full spectral treatment of all damping and cascading processes is needed to reliably calculate the beam dissipation rate, but that is a topic of future investigations.

The effect of collisional damping in the beam dissipation process is substantially small in the sense that the relaxation time is much larger than the IC scattering time even in their presence. Moreover, it is inversely dependent on the IGM temperature. In principle, for a sufficiently cold plasma, collisions could play a major role in beam energy dissipation, despite it being one of the slowest processes. Other processes like wave-wave scattering or wave-particle scattering produce higher rates thereby dominating the effect of particle collisions (Chang et al., 2014). Other commonly ignored processes may play a role as well, and so an accurate estimate of the dissipation rate, and consequently, the relaxation time requires one to consider all possible damping mechanisms. Since collisional damping is not a predominant factor in the beam dissipation process, a direct constraint on the IGM temperature from its effect appears to be out of reach.

There is observational evidence for beam dissipation arising from the absence of GeV haloes around the jets of misaligned AGN (Broderick et al., 2018). This paper demonstrates that theoretically calculating the beam dissipation rate is very challenging on account of the inherent non-linearity, and simple estimates may be misleading. It is conceivable that beam dissipation is efficient only for beams of a certain density, which translates to a maximum distance, $L_{\mathrm{max}}$, from the AGN that depends on its multi-TeV gamma-ray flux. We are not yet in the position to reliably calculate that distance and hence to estimate which part of the cascade would be quenched (at $L<L_{\mathrm{max}}$) and what spectrum the observable cascade produced at $L>L_{\mathrm{max}}$ would have.

Appendix A Approximate expression for $\omega_{NL}({\bf{k}})$ in 3D model

Eq. (6) for $\omega_{NL}({\bf k})$ is a three dimensional integral that can be written as

[TABLE]

where

[TABLE]

It is easy to see that for

[TABLE]

one can neglect $\cos\phi^{\prime}$ under the square roots in Eq. (A2) yielding

[TABLE]

while in the region $k_{\perp}^{\prime}\sim k_{\perp}\wedge|k_{\parallel}^{\prime}-k_{\parallel}|\lesssim\sqrt{2}k_{\perp}$

[TABLE]

If $\phi^{\prime}$ is not very close to 0 so that $|1-\cos\phi^{\prime}|$ is order of 1 ($0.3\pi\lesssim\phi^{\prime}\lesssim 1.7\pi$, $\cos\phi\lesssim 0.5$) then the integral (A5) can be again estimated neglecting $\cos\phi^{\prime}$ under the square roots with the result which is order of Eq. (A4) at $k_{\perp}^{\prime}\sim k_{\perp}$:

[TABLE]

At the same time, the contribution to the integral (A5) from the angles satisfying $|1-\cos\phi^{\prime}|\ll(k_{\parallel}^{\prime}-k_{\parallel})^{2}/(2k_{\perp}^{2})$ ($\phi^{\prime}\lesssim\Delta\ll\sqrt{2}|k_{\parallel}^{\prime}-k_{\parallel}|/k_{\perp}$) is

[TABLE]

which is much smaller compared to Eq. (A6). Thus, Eq. (A4) can be used in the whole $(k_{\perp}^{\prime},k_{\parallel}^{\prime})$-plane to approximately calculate the $\phi^{\prime}$-integral in the expression for $\omega_{NL}({\bf k})$.

Appendix B Approximate expression for $\omega_{NL}({\bf{k}})$ in 2D model

In the 2D model

[TABLE]

where

[TABLE]

The same proceeding as in Appendix A yields

[TABLE]

Appendix C Applicability of 2D model

Here, we provide arguments for the applicability of the 2D model to study Eq. (4). The evolution of the spectrum in 3D as well as 2D exhibits two characteristic stages: linear and non-linear. During the linear stage, the non-linear Landau damping is negligible compared to the growth rate, because the spectral energy $W({\bf k})$ is too small. The system turns to the non-linear regime when the NL damping rate becomes comparable to the maximum growth rate of the instability. Let us denote the time when it takes place by $t=t_{NL}$. We want to demonstrate now that the 2D model reflects the main features of the original 3D one during the both linear and non-linear regimes.

C.1 Linear regime

During the linear regime $t\ll t_{NL}$, the wave energy growth exponentially $\propto\exp(2\max(\omega_{i})t)$ in both 3D and 2D model where the maximum growth rate $\max(\omega_{i})$ in 2D is the same as in 3D.

At the transition time $t\approx t_{NL}$, the NL damping starts to balance the instability growth. We can estimate the electric filed energy at this time. Let us consider a time $t$ close enough to $t_{NL}$ and such that the electric filed energy is still concentrated near the wave vector with the maximum growth rate (at $t\gg t_{NL}$, the most waves will be scattered to small wave vectors where the growth rate is zero). In the 3D model, the most unstable wave number is located near $k\sim\omega_{p,e}/c$ ($k_{\perp}\sim\omega_{p,e}/c$ in 2D). In the vicinity of this wave vector, we can neglect the exponential function in Eq. (6) (Eq. (13) in 2D) because the factor $a\ll 1$. Additionally, we replace $({\bf k}^{\prime}{\bf k})^{2}/(k^{\prime 2}k^{2})$ in Eq. (6) by its angle-averaged value $1/3$ ($({\bf k}_{\perp}^{\prime}{\bf k}_{\perp})^{2}/(k_{\perp}^{\prime 2}k_{\perp}^{2})$ in Eq. (13) is replaced by $1/2$). Then, since the width of growth rate in the parallel direction $\Delta k_{\parallel}$ in 3D is extremely narrow ($\Delta k_{\parallel}\ll\omega_{p,e}/c$), we also adopt ${\bf k}_{\parallel}^{\prime}\approx{\bf k}_{\parallel}$ and ${\bf k}_{\parallel}^{\prime}\approx-{\bf k}_{\parallel}$ for positive and negative $k_{\parallel}^{\prime}$, accordingly. Finally, the NL damping rate in 3D and 2D cases becomes, respectively,

[TABLE]

and

[TABLE]

Since the term $|{\bf k}_{\perp}^{\prime}-{\bf k}_{\perp}|/|{\bf k}_{\perp}^{\prime}-{\bf k}_{\perp}-2{\bf k}_{\parallel}|<1$ and $W_{2D}(k_{\perp})\approx W_{3D}({\bf k})\Delta k_{\parallel}$, Eq. (C2) differs from Eq. (C1) at most by factor 3/2. Let us continue to work with Eq. (C2) to derive the electric field energy at $t\approx t_{NL}$. Using the same technique as in Appendix A and B, Eq. (C2) yields

[TABLE]

where the integral is calculated in the vicinity the most unstable mode $k_{\perp}\approx\omega_{p,e}/c$. Then, $k_{\perp}^{\prime}$ is close to $k_{\perp}$ and Eq. (C3) becomes

[TABLE]

where $W_{E}^{tot}$ is the total electric field energy density and $\Delta k_{\perp}$ in the characteristic spectral width in the perpendicular direction at $t\approx t_{NL}$. At the transition point, $\omega_{NL}\approx\max(\omega_{i})$ and

[TABLE]

Since Eq. C1 differs from Eq. C2 only by factor $2/3$, $W_{E}^{tot}$ in the 3D model will be of the same order as Eq. (C5) derived for the 2D one. Furthermore, our numerical results (see section 3) show that $W_{E}^{tot}$ constitutes about $1$ % of the the initial beam energy density. Therefore, an additional factor $2/3$ in the 3D model won’t principally affect this result. Therefore, the the wave scattering in the ${\bf k}_{\parallel}$-direction absent in the 2D model is not crucial for the beam energy loses during the linear phase and can be, indeed, neglected.

C.2 Non-linear regime

Now, we turn to the non-linear spectral evolution at times $t\gg t_{NL}$. After the NL damping becomes significant and scatters the most wave energy to lower wave numbers, the further energy growth of the electric field is regulated by the spectrum at large wave numbers ($k_{\perp}c/\omega_{p,e}\gg 1$) where both LL and NL damping rates are smaller than the instability growth rate (Chang et al., 2014). However, the spectrum at $k_{\perp}c/\omega_{p,e}\approx 1$ also contributes to the electric field growth at $t\gg t_{NL}$ if particle collisions are present as we found in our 2D model (see section 4). There should exist a similar effect in the 3D model near the region of the maximum growth rate $kc/\omega_{p,e}\approx 1$. Thus, there can be two resonant groups of waves driving the instability at $t\gg t_{NL}$. At the same time, there is also a group of non-resonant waves at small wave vectors $k\sim 0$ ($k_{\perp}\sim 0$ in 2D). We sketched these three different parts of the spectrum schematically in Fig. 14. We note that the position of the red circle (corresponding to resonant waves growing when the collisional damping is present) is slightly different for the 3D and 2D models. The reason simply because the maximum growth rate in our 2D model occurs at $k_{\perp}c/\omega_{p,e}=1$, while in the 3D one it is $k_{\perp}c/\omega_{p,e}<1$. Another difference between 2D and 3D is that the spectrum in 3D is actually symmetric for negative $k_{\parallel}$, the 2D model treats the spectrum evolution only in the plane $k_{\parallel}c/\omega_{p,e}=1$. We will show now that the 2D and 3D models provide close values of the NL damping rate in the resonant regions (blue and red circles in Fig. 14) which are decisive for the beam energy losses, because only the resonant waves can interact with the beam.

C.2.1 Resonant waves at large wave vectors

Let us first consider the resonant region marked by the blue circle in Fig. 14. It is rather straightforward that the regions located at much smaller wave vectors (green and and red circles) in 2D and 3D provide close contributions to the NL damping, since the difference between the positions of green and red circles in the 2D and 3D models is negligible for the NL damping rate at large wave vectors. To demonstrate this explicitly, we can neglect $|{\bf k}^{\prime}|\ll|{\bf k}|$ and approximately replace $({\bf k}^{\prime}{\bf k})^{2}/(k^{\prime 2}k^{2})$ in Eq. (6) by its angle-averaged value $1/3$. Then, Eq. (6) becomes

[TABLE]

A similar derivation in 2D leads to

[TABLE]

The total energies in Eqs. (C6) and (C7) are expected to be comparable, since we found above that the electric field gains close energies in the 2D and 3D model during the linear stage. Then the difference between Eq. (C6) and (C7) is only factor $2/3$. Thus, our 2D model provides a reasonable form of spectrum at small wave vectors that does not considerable modify the NL damping in the resonant region at large wave numbers. To finish the analysis of the NL damping in the resonant region at large wave vectors, we need to consider the contribution into NL damping rate from another resonant region at large wave vectors but with negative $k_{\parallel}<0$. We note that this region is absent in the 2D model, since we study the problem in the plane $k_{\parallel}c/\omega_{p,e}=1$. We will show below that this contribution in 3D is negligible and our 2D model is justified in this respect. Setting ${\bf k}_{\perp}^{\prime}\approx{\bf k}_{\perp}$ and $k_{\parallel}^{\prime}\approx-k_{\parallel}$, we arrive at

[TABLE]

Except from the region where the exponential factor in Eq. (C6) is small enough, Eq. (C6) gives a much larger contribution into the NL damping rate than Eq. (C8), as $\Delta k_{\parallel}\ll k$ and the most electric field energy energy is concentrated at small wave vectors. Summarizing the above derivation, the NL damping rate at large $kc/\omega_{p,e}\gg 1$ in 3D is close to that in 2D.

C.2.2 Resonant waves at small wave vectors due to collisional effect

Now, let us consider another resonant region that can appear at $t\gg t_{NL}$ due to collisional damping of the electric field energy (see section 4) and shown by the red circle in Fig. 14. First, we note that the position of the red circle that we expect in 3D is not far from that one in 2D, since the maximum growth rate in 3D occurs at $k_{\perp}c/\omega_{p,e}\approx 0.6$, while in 2D it is $k_{\perp}c/\omega_{p,e}\approx 1$. Second, similar to section C.2.1 it can be demonstrated that the integral contribution into the NL damping rate from resonant waves at large wave vectors (blue circle) differs slightly for the 2D and 3D models, because the difference in the wave vectors between those groups of waves is large. Finally, it is left to compare the NL damping rate due to the interaction with the non-resonant waves (green circle) (it can be proved in a similar manner as in section C.2.1 that the contribution from the symmetric region at $k_{\parallel}c/\omega_{p,e}\approx-1$ in 3D is smaller than from the non-resonant waves at $k\approx 0$ in 3D). To do this, we can utilize similar simplifications as we used above to derive Eq. (C1). Again, we neglect the exponential factor and replace $({\bf k}^{\prime}{\bf k})^{2}/(k^{\prime 2}k^{2})$ in Eq. (6) by its angle-averaged value $1/3$

[TABLE]

A similar derivation in 2D leads to

[TABLE]

Thus, the NL damping rates in 2D and 3D are close to each other also for resonant waves at small wave vectors. This completes the analysis of applicability of the 2D model.

Appendix D Numerical comparison of 3D and 2D models

To compare the 3D and 2D models, we chose the parameters of case 3 in Table 1, since the simulation region is the smallest ($k_{\perp}c/\omega_{p,e}<35$, $k_{\parallel}c/\omega_{p,e}\lesssim 1$) among all cases. Although the spatial resolution in 3D has been compromised owing to computational costs, the calculations should provide a reasonable estimate of the comparability of the 2D and the 3D models in the linear and early nonlinear stages of evolution. The region of instability ($\delta k_{||}c/\omega_{p,e}=10^{-6}$ at $k_{||}c/\omega_{p,e}=1$) is resolved with 7 linearly spaced grid points, and we employ 60 logarithmically spaced grid points between $k_{||}c/\omega_{p,e}=10^{-7}$ and $k_{||}c/\omega_{p,e}=1$ to cover the region to which NL damping transfers wave energy. In the transverse direction, we use 40 logarithmically spaced grid points to resolve the spectrum at $k_{\perp}c/\omega_{p,e}\leq 1$, where the growth rate remains essentially independent of $k_{\perp}$. In the band $k_{\perp}c/\omega_{p,e}\geq 1$, we have 110 linearly spaced grid points extending up to $k_{\perp}c/\omega_{p,e}=35$, beyond which linear Landau damping prohibits wave growth (cf. Fig. 3).

Fig. 15 (a) illustrates the total electric field energy density normalized to the beam energy density for the 2D and 3D scenarios. It is clearly seen that the energy evolution in 2D case replicates the 3D case except for the noise level. In 2D, the noise is given by $k_{B}T_{\mathrm{IGM}}\Delta k_{\parallel}$ to account for the energy concentrated in the actual region of instability ($\Delta k_{\parallel}c/\omega_{p,e}=10^{-6}$), whereas for 3D it is simply given by $k_{B}T_{\mathrm{IGM}}$. Nevertheless, the electric field energy turns to quasi-saturation at the value much smaller than the initial beam energy density. A similar comparison of the power losses is shown in Fig. 15 (b). The power losses in 3D were calculated similar to Eq. (17):

[TABLE]

Although the nonlinear effects set in slightly earlier in the 3D scenario, the power losses in the 2D model are in good agreement with the 3D one. Finally, we conclude that our 2D model can be reasonably used to study Eq. (4).