# The role of resonant plasma instabilities in the evolution of blazar induced
pair beams

Roy Perry1, Yuri Lyubarsky1
1Physics Department, Ben-Gurion University, P.O.B. 653, Beer-Sheva 84105,
Israel

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The fate of relativistic pair beams produced in the intergalactic medium by
very-high energy emission from blazars remains controversial in the
literature. The possible role of resonance beam plasma instability has been
studied both analytically and numerically but no consensus has been reached.
In this paper, we thoroughly analyze the development of this type of
instability. This analysis takes into account that a highly relativistic beam
loses energy only due to interactions with the plasma waves propagating within
the opening angle of the beam (we call them parallel waves), whereas
excitation of oblique waves results merely in an angular spreading of the
beam, which reduces the instability growth rate. For parallel waves, the
growth rate is a few times larger than for oblique ones, so they grow faster
than oblique waves and drain energy from the beam before it expands. However,
the specific property of extragalactic beams is that they are extraordinarily
narrow; the opening angle is only $\Delta\theta\sim 10^{-6}-10^{-5}$. In this
case, the width of the resonance for parallel waves,
$\propto\Delta\theta^{2}$, is too small for them to grow in realistic
conditions. We perform both analytical estimates and numerical simulations in
the quasilinear regime. These show that for extragalactic beams, the growth of
the waves is incapable of taking a significant portion of the beam’s energy.
This type of instability could at best lead to an expansion of the beam by
some factor but the beam’s energy remains nearly intact.

###### keywords:

instabilities, plasmas, intergalactic medium, jets

††pubyear: 2021††pagerange: The role of resonant plasma instabilities in the
evolution of blazar induced pair beams–13

##  1 Introduction

The high-energy relativistic jets from blazars emit gamma ray radiation, which
travels along in the same direction as the jets. During their course through
the intergalactic medium (IGM), the very high energy (VHE) gamma photons
interact with the extragalactic background light (EBL) and the cosmic
microwave background (CMB). This interaction has a mean free path ranging from
several Mpc to 1 Gpc, which leads to the production of ultra-relativistic
electron-positron pairs with Lorentz factors of about $\gamma\sim
10^{4}-10^{7}$. These relativistic beams propagate in the IGM and are expected
to lose energy by inverse-Compton scattering upon the CMB. However, the
observed gamma-ray bump associated with the product photons of inverse-Compton
cascades (ICC) is lower than predicted (e.g. Aharonian et al. 2006; Neronov &
Vovk 2010; Broderick et al. 2012). One possible explanation is pair beam
deflection by the inter-galactic magnetic field. Measurement of the GeV flux
on telescopes in conjunction with the estimation of blazar-originated gamma
rays flux, a lower bound on the intergalactic magnetic field is predicted. The
inconsistencies between the predicted and measured flux are attributed to the
deflection of particles by the magnetic field (Archambault et al., 2017; Tiede
et al., 2017; Ackermann et al., 2018). An alternate possibility (Broderick et
al., 2012; Schlickeiser et al., 2012) is that before the beam pairs up-scatter
photons, they lose energy through excitation of plasma waves via the resonance
beam-plasma instability. The beam-plasma instability results in excitation of
plasma oscillations (Langmuir waves) with a frequency close to the plasma
frequency (e.g., Krall & Trivelpiece 1973). The IGM density is of the order of
$10^{-7}\text{cm}^{-3}$, so the plasma frequency is of the order of a few Hz,
while the beam density ranges over many orders of magnitude, though estimated
to be no more than $10^{-21}\text{cm}^{-3}$ (e.g., Broderick et al. 2012). The
instability could develop in two different regimes: reactive (hydrodynamic)
and kinetic (e.g., Boyd et al. 2003). The first applies to an essentially cold
beam, whose thermal spread is small enough such that all the particles are in
resonance with the wave. The second corresponds to the case of a wide spread
in the distribution function of the beam. In that case, only a fraction of the
particles contribute to the growth of the wave. If the beam is initially
narrow enough such that the instability develops in the hydrodynamic regime,
then mostly oblique waves are excited. Therefore the beam expands while losing
almost no energy (Fainberg et al., 1970). Eventually, the momentum spread of
the beam becomes large enough for the kinetic regime to come into play. During
the kinetic regime, the beam could lose most of its energy by exciting plasma
waves which eventually heats the background plasma (Fainberg et al., 1970;
Rudakov, 1971; Breǐzman & Ryutov, 1971). The literature on the collective
plasma processes in extragalactic pair beams is extensive but its conclusions
still remain controversial. Miniati & Elyiv (2013) studied different
mechanisms that affect the development of resonance instability and concluded
that for typical beam parameters, either non-linear wave interactions or
inhomogeneities of the background plasma suppress this instability. However,
no consensus has been reached about the role of the collective plasma
processes in the evolution of extragalactic beams (e.g., Broderick et al.
2012; Broderick et al. 2014; Schlickeiser et al. 2012, 2013; Miniati & Elyiv
2013; Chang et al. 2014; Krakau & Schlickeiser 2014; Supsar & Schlickeiser
2014; Chang et al. 2016; Sironi & Giannios 2014; Kempf et al. 2016; Rafighi et
al. 2017; Broderick et al. 2018; Shalaby et al. 2018; Vafin et al. 2018, 2019;
Yan et al. 2019; Alves Batista et al. 2019; Shalaby et al. 2020).

The results of PIC simulations could not be interpreted unambiguously because
the width of the wave-particle resonance is so small in extragalactic beams
that it could not be resolved in simulations. Moreover, the numerical noise in
simulations is high so the beam parameters had to be chosen such that the
level of the excited plasma waves were high enough. In this case different
non-linear processes play a significant role, so the results could hardly be
rescaled to the realistic range of parameters.

Up until now, the discussion was focused on the conditions for the development
of beam instability. It has been taken for granted that if the instability
develops in the kinetic regime, it is capable of decelerating the beam. In
this paper, we take into account a subtle fact that only waves propagating
within the open angle of the beam are responsible for energy loss in the beam
whereas excitation of oblique waves merely results in an increase in its
angular spread, which reduces the instability rate. Our main goal is to study
the pair beam’s evolution and behaviour in time, paying special attention to
the relevant range of parameters for intergalactic beams. Namely, an extremely
small growth rate and an extremely small angular opening of the beam. On the
one hand, the last is sufficient for a transition to the kinetic instability
regime (due to the very small growth rate), but on the other hand, it makes
the width of the wave particle resonance extremely small, which is crucial in
the evolution of the instability.

The instability of relativistic beams was comprehensively studied in the early
plasma literature (Bludman et al., 1960; Watson et al., 1960; Fainberg et al.,
1970; Breǐzman & Ryutov, 1971; Rudakov, 1971). However, most of the attention
was given to relatively powerful beams, for which the kinetic regime is
established only when the beam angular width is not too small. In these cases
the beam efficiently loses energy in the kinetic regime. For weak
intergalactic beams, the instability develops in the kinetic regime when the
beams are extremely narrow, in which case the conditions for energy transfer
from the beam to the plasma waves become unrealistically restrictive. We show
that even if the instability develops, the beam’s energy remains practically
unchanged, and its angular speard expands just a few times over.

This paper is organized as follows: In section 2 we introduce the beam-plasma
instability mechanism. In section 3 we describe the hydrodynamic and kinetic
regimes of the instability and give explicit formulas for the growth rates. In
section 4 we estimate the conditions for the instability and for the effective
energy loss of the beam for realistic IGM parameters. In section 5 we present
numerical simulations of the beam evolution by solving the kinetic equations
in the quasi-linear regime. In section 6 we conclude and summarize our work.

##  2 Preliminary Considerations

###  2.1 Beam-wave resonance

The beam-plasma instability (e.g., Krall & Trivelpiece 1973) results in
excitation of plasma oscillations (Langmuir waves) satisfying the resonance
(Cerenkov) condition:

$\omega={\bf{k}}\cdot{\bf{v}},$ (1)

where $\omega$ and ${\bf{k}}$ are the frequency and wave vector of the wave.
For dilute beams, $n_{b}\ll n_{p}$, the frequency of the excited waves is
close to the plasma frequency,

$\omega_{p}=\sqrt{\frac{4\pi n_{p}e^{2}}{m_{e}}},$ (2)

where $m_{e}$, $e$ are the electron mass and charge, respectively, and $n_{p}$
is the plasma density.

The extragalactic beams are highly energetic; the typical particle Lorentz
factor is $\gamma\sim 10^{4}-10^{7}$. Even though the beam energy spread is
large, $\Delta\gamma\sim\gamma$, the velocity spread is negligibly small,
$\Delta v\sim\Delta\gamma/\gamma^{3}$. Therefore we can describe each beam
particle velocity in spherical coordinates (on the azimuthal plane $\phi=0$)
as

${\bf{v}}=c\left(1-\theta^{2}/2,\theta\right),$ (3)

where the first and the second components correspond to the particle velocity
components along and perpendicular to the beam direction, respectively. The
wave-vector, ${\bf{k}}$, can be described in a three-dimensional spherical
coordinates space as:

${\bf{k}}=\left(k,\theta_{0},\phi_{0}\right).$ (4)

The resonance condition (1) can then be written as:

$\cos\phi_{0}=\frac{\omega_{p}/kc-\cos\theta\cos\theta_{0}}{\sin\theta\sin\theta_{0}}.$
(5)

One can therefore describe the resonance wave-vector $\bf{k}$ with two
independent parameters, $k$ and $\theta_{0}$.

The geometry of the particle-wave resonance in velocity space is shown in
Figure 1. The arrows represent the directions of propagation of the Langmuir
waves, where $\theta_{0},\theta^{\prime}_{0}$ are the angles of these vectors
with respect to the beam axis, chosen as the horizontal axis. The arc
represents the distribution of the pair-beam particles, ranging over an
angular spread of $\Delta\theta\ll 1$. The range of resonant waves is
signified with the bordered arrow between the two dashed lines. For any
$\theta_{0}$, only waves in this range of $\bf{k}$ can exchange energy and
momentum with the beam particles. The figure shows that the range of resonant
waves, $\Delta k$, decreases with decreasing wave vector angle. Indeed, to
first order in $\Delta\theta$, the range of resonance $\Delta k$ behaves as
$\sim\tan\theta_{0}\Delta\theta$. Therefore $\Delta k/k\sim\Delta\theta$ at
$\theta_{0}\sim 1$, but when the wave propagates within the opening angle of
the beam, $\theta_{0}\sim\Delta\theta$, the resonance becomes very narrow,
$\Delta k/k\sim\Delta\theta^{2}$.

Figure 1: (Color online) The beam-wave resonance condition in velocity space.
Shown are the particle velocity distribution (thick blue arc with an opening
angle $\Delta\theta\ll 1$) and two plasma waves at two angles
$\theta_{0}\gg\Delta\theta$ and $\theta^{\prime}_{0}\sim\Delta\theta$ (red and
green arrows). The range of the phase velocities which are in resonance with
the beam is larger for the first wave than for the second one.

The angular spread of intergalactic beams is very small. After pair-
production, the angular distribution of the particles is quite narrow because
of their direct emission from a point-like source, far away from the
intergalactic plasma. Initially, the spread is about
$\Delta\theta\sim\gamma^{-1}$. The beam could expand as a result of the
instability but not by much. Therefore the width of the resonance is very
small even for oblique waves, whereas for waves nearly aligned with the beam,
it is extraordinarily small. We show below that this has a profound effect on
the evolution of the beam because on the one hand, only nearly aligned waves
may be responsible for any energy loss of the beam, but on the other hand,
they easily violate the resonance conditions. Therefore, their growth may be
easily suppressed by weak perturbations, such as inhomogeneity of the
background plasma or/and non-linear wave interactions.

###  2.2 Energy loss vs expansion; qualitative consideration

The particle distribution of the beam evolves because of the induced emission
and absorption of plasma waves. The energy and momentum conservation
principles in a single emission/absorption event are described via the
following equations:

$E=E^{\prime}\pm\hbar\omega_{p};\qquad{\bf{p}}={\bf{p}}^{\prime}\pm\hbar{\bf{k}}.$
(6)

Here we denote the parameters of the beam particle before and after the wave
excitation by unprimed and primed quantities, respectively. For concreteness,
consider an emission event and assume for simplicity that the beam axis, the
particle velocity and the wave vector are all in the same plane. This imples:

$\delta\gamma=\gamma^{\prime}-\gamma=-\frac{\hbar\omega_{p}}{mc^{2}},$ (7)
${}\begin{array}[]{l}m\gamma^{\prime}v^{\prime}\sin\theta^{\prime}-m\gamma
v\sin\theta=-\hbar k\sin{\theta_{0}},\\\
m\gamma^{\prime}v^{\prime}\cos\theta^{\prime}-m\gamma v\cos\theta=-\hbar
k\cos{\theta_{0}},\\\ \end{array}$ (8)

where $\theta$, $\theta_{0}$ are the particle and wave vector angles with
respect to the beam axis, respectively. We consider a very narrow,
$\theta,\theta^{\prime}\ll 1$, highly relativistic, $\gamma\gg 1$, beam and
expand up to first order. The pair of momentum equations now become:

$\displaystyle\cos\left({\theta-{\theta_{0}}}\right)\left({\gamma^{\prime}\theta^{\prime}-\gamma\theta}\right)=-\frac{{\hbar{\omega_{p}}}}{{m{c^{2}}}}\sin{\theta_{0}},$
(9)
$\displaystyle\cos\left({\theta-{\theta_{0}}}\right)\left({\gamma^{\prime}-\gamma}\right)=-\frac{{\hbar{\omega_{p}}}}{{m{c^{2}}}}\cos{\theta_{0}},$

Upon dividing the first equation in (9) by the second one and rearranging
terms, we get the ratio between relative angular and energetic change,

${}\left|{\frac{{\delta\theta}}{\theta}\frac{\gamma}{{\delta\gamma}}}\right|=\frac{{\tan{\theta_{0}}}}{\theta}.$
(10)

Eq. (10) shows that the outcome of the emission/absorption process strongly
depends on the relative angles between the particle and the wave. Namely, the
relative change in energy is comparable to the relative change in angle only
when the wave is directed within the opening angle of the beam,
$\theta_{0}\sim\Delta\theta$. In the case $\theta_{0}\gg\Delta\theta$, the
particle’s angle varies significantly more than its energy.

This conclusion may be considered a preliminary one because it is based only
on analysis of a single emission/absorption event. The particle distribution
function evolves as a result of many interactions. The evolution may be
described as diffusion in the momentum space (e.g., Brejzman 1990). In section
5.1, we show rigorously, by analysing the momentum diffusion coefficients,
that only quasi-parallel waves ($\theta_{0}\sim\Delta\theta$) take away any
energy efficiently whereas the interaction with oblique waves
($\theta_{0}\gg\Delta\theta$) leads only to beam expansion.

The evolution of the beam is governed by the total spectrum of resonant plasma
oscillations. The dominant part of the spectrum will determine whether the
beam loses energy or expands in angle. The instability mechanism will enhance
the spectrum of waves in each direction according to their respective growth
rate. Those which grow the fastest are the ones who primarily control the
evolution of the beam. Since the growth is exponential, even a moderate excess
in growth rate will make the corresponding range of the wave spectrum
dominant. It follows from eq. (10), that the beam will lose its energy, if
small angle waves ($\theta_{0}\sim\Delta\theta$) grow faster. Indeed,
according to kinetic linear theory, the growth rate is larger for such waves
(Breǐzman & Ryutov 1971; Rudakov 1971, see also section 3.2). However, the
intergalactic beams are very narrow. Therefore the resonance width is
especially narrow, $\sim\Delta\theta^{2}$, for these waves, as we have seen in
the previous subsection. In this case, even a weak inhomogeneity of the IGM
or/and weak nonlinear processes result in loss of the resonance between these
waves and the beam, which suppresses their growth rate.

###  2.3 Effect of plasma inhomogeneity

When a wave propagates through a steady but weakly inhomogeneous plasma, the
wave vector varies according to the geometrical optics prescription

$\frac{d\mathbf{k}}{dt}=-\nabla\omega_{p}.$ (11)

Taking into account that $k\sim\omega_{p}/c$ one sees that if the
characteristic scale of the density variation is $L$, such that
$|\nabla\omega_{p}|=\omega_{p}/L$, the relative change in the wave vector
after time $t$ would be

$\frac{\delta k}{k}\sim\frac{ct}{L}.$ (12)

Let a wave be excited initially in the resonance range, bounded by the two
dashed lines in Figure 1. As it travels through the plasma, its wave number
moves along that line, towards one of the resonance boundaries, until it
crosses that boundary and exits the resonance region. If the exit time is
smaller than the instability time, the instability does not develop (Breizman
& Ryutov, 1970). Miniati & Elyiv (2013) have shown that inhomogeneity and
other effects (non-linearity, collisional damping) impose strong constraints
on the development of the instability of intergalactic beams. Since then,
these effects were studied by many authors quoted in our introduction. In this
paper, we show that the constraints on the energy exchange between the beam
and the background plasma are much more severe because of a special role
played by quasi-parallel waves. We have seen that the resonance width
decreases as the angle between the wave and beam directions decreases. Hence,
the inhomogeneity causes a loss of resonance much faster for quasi-parallel
waves than for oblique waves. We concentrate on the effect of inhomogeneity.
Non-linear effects could also suppress the instability by removing the waves
from the resonance. However, we show in section 4 that in the context of
intergalactic beams, inhomogeneity alone completely suppresses the generation
of quasiparallel waves, responsible for energy loss of the beam.

It was claimed recently (Shalaby et al., 2018, 2020) that the inhomogeneity
does not suppress the instability. The authors considered a monochromatic 1D
beam propagating through the density distribution with a local minimum and
applied periodic boundary conditions. They claim that within a "bowl", a
standing plasma wave is formed between two turning points, which remain in
resonance with the beam all the time. This artificial example could not be
considered a disproof of the general effect outlined above. First of all, the
instability develops only in a small region in the vicinity of the local
density minimum point. In the the 3D case, the effects of wave locking are
irrelevant for most of the system’s volume, even if there are many density
"bowls". Moreover, in 1D simulations with periodic boundary conditions, the
beam enters the bowl already appropriately modulated during the previous
passages, therefore the growth of perturbations proceeds further, which does
not happen in the real 3D world.

Additionally, it is worth stressing that the work by Shalaby et al. 2018, 2020
is based on an implicit assumption that standing waves are formed in the
region of interest. In 1D, standing waves are formed by back and forward
directed waves bouncing between the turning points. In 3D, the wave changes
direction at the turning point and therefore could not remain in resonance
with the beam. The group velocity of such a non-resonant plasma wave is
extremely small so that the global standing wave, in which the phases and the
wave vectors of local oscillations are adjusted within the whole region, are
not formed in the system.In any case, even if the theory by Shalaby et al.
(2018, 2020) is universally correct, their instability growth rate is still
smaller than the collisional decay rate of the plasma waves in the IGM. In
sect. 4 below, where parameters of the beam are discussed, we present the
corresponding estimate.

##  3 Beam-plasma instability

The growth rate of the instability $\Gamma$ can be calculated by making use of
Maxwell’s equations and the linearized Vlasov equation (e.g., Krall &
Trivelpiece 1973). It is common in the literature to describe this instability
mechanism in two different regimes, i.e. the hydro (reactive) regime and the
kinetic regime.

###  3.1 Hydrodynamic regime

The hydrodynamic regime applies to waves which resonate with the beam as a
whole. In other words, the beam is narrow enough such that the parallel
velocity component of all the particles matches the phase velocity of the
wave. For every wave vector whose parallel component satisfies
$v_{b}k_{\parallel}=\omega$, where $v_{b}$ is the typical velocity of the beam
particles, the interaction is unstable but waves with different perpendicular
components grow at different rates.

The dispersion relation for unstable Langmuir waves in the electrostatic
approximation looks like (Bludman et al., 1960):

$1-\frac{{\omega_{p}^{2}}}{{{\omega^{2}}}}-\frac{{\omega_{b}^{2}}}{{{{\left({\omega-{k_{\parallel}}{{\rm{v}}_{b}}}\right)}^{2}}}}\left({\frac{{{\gamma^{2}}{{\sin}^{2}}{\theta_{0}}+{{\cos}^{2}}{\theta_{0}}}}{{{\gamma^{3}}}}}\right)=0,$
(13)

where $\omega_{b}=\sqrt{4\pi e^{2}n_{b}/m_{e}}$ is the beam plasma frequency.
There is a wide range of unstable solutions of this equation but in the case
of interest, $n_{b}/n_{p}\sim 10^{-15}-10^{-18}$, $\gamma\sim 10^{6}$, only
the plasma waves in resonance with the beam grow fast enough (see Appendix A).
The maximal growth rate of the resonant plasma wave, is (Bludman et al., 1960;
Fainberg et al., 1970)

${}{\Gamma_{{\rm{hydro}}}}=\frac{{\sqrt{3}}}{{{2^{4/3}}}}\frac{{{\omega_{p}}}}{\gamma}{\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)^{1/3}}{\left({\cos{\theta_{0}}}\right)^{2/3}}{\left({1+{\gamma^{2}}{{\tan}^{2}}{\theta_{0}}}\right)^{1/3}}.$
(14)

This implies that in the hydrodynamic regime, the fastest growing waves are
oblique ones, i.e. waves at angle $\theta_{0}\sim 1$. The growth rate for such
waves is roughly

${\Gamma_{{\text{hydro}}}}\left(\theta_{0}\sim
1\right)\sim{\omega_{p}}{\left({\frac{{{n_{b}}}}{{\gamma{n_{p}}}}}\right)^{1/3}},$
(15)

while for waves which are almost parallel to the beam axis, i.e.
$\theta_{0}\sim\Delta\theta\sim\gamma^{-1}$, the growth rate reduces to

${}{\Gamma_{{\rm{hydro}}}}\sim\frac{{{\omega_{p}}}}{\gamma}{\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)^{1/3}}.$
(16)

Therefore the beam evolution is dominated by oblique waves. In this case, the
beam angular spread will expand, without significantly losing energy (Fainberg
et al., 1970).

As the beam broadens, the transverse velocity dispersion, $\Delta
v_{\perp}\sim c\Delta\theta$ increases accordingly and eventually the beam
could not be considered ‘cold’. The hydrodynamic regime is valid as long as
all the particles interact with the same wave, which implies:

$\Gamma\gg{{\bf{k}}\cdot\Delta{\bf{v}}},$ (17)

where $\Delta{\bf{v}}$ is the width of the particle velocity distribution. For
a beam with angular spread $\Delta\theta$, the parallel and transverse
velocity spreads are $\Delta v_{\perp}=c\Delta\theta$ and $\Delta
v_{\parallel}=c{\Delta\theta^{2}}/2$, respectively, so that the hydrodynamic
regime condition reduces to

${\omega_{p}}\left({\frac{{\Delta{\theta^{2}}}}{2}+\tan{\theta_{0}}\Delta\theta}\right)\ll{\Gamma_{hydro}},$
(18)

where we have set the resonance condition, $k_{\parallel}=\omega_{p}/c$. Now
one can write

$\displaystyle\Delta\theta\ll{\left({\frac{{{n_{b}}}}{{\gamma{n_{p}}}}}\right)^{1/3}}$
$\displaystyle{;\quad{\theta_{0}}\sim 1},$ (19)
$\displaystyle{\gamma\Delta\theta^{2}\ll{\left(\frac{n_{b}}{n_{p}}\right)}^{1/3}}$
$\displaystyle{;\quad{\theta_{0}}\sim\Delta\theta}.$

We show in section 4 that in agreement with Miniati & Elyiv (2013), the
condition for the hydrodynamic regime is mostly violated for characteristic
IGM and beam parameters so that the only applicable regime for the instability
mechanism is the kinetic one.

###  3.2 Kinetic regime

The growth rate of a Langmuir wave, in the kinetic regime is given by (e.g.,
Brejzman 1990)

$\Gamma={\omega_{p}}\frac{{2{\pi^{2}}{e^{2}}}}{{{k^{2}}}}\int{\left({\bf{k}}\cdot\frac{{\partial
f}}{{\partial{\bf{p}}}}\right)\delta\left({{\omega_{p}}-{\bf{k}}\cdot{\bf{v}}}\right){d^{3}}p}.$
(20)

Here the beam distribution function is normalized such that

$\int{f{d^{3}}p}=n_{b}.$ (21)

The beam is highly relativistic, which means that the magnitude of ${\bf{v}}$
is $c$ up to second order in $\gamma^{-1}$. Hence the delta function argument
contains only the angles of the beam momentum vector. Therefore, integration
of the delta function is performed over the azimuthal angle, $\phi$, and the
growth rate function could be written as an integral solely over the polar
angle. If one introduces the angular distribution function:

${}g\left(\theta\right)=\int_{0}^{2\pi}{d\phi}\int_{0}^{\infty}{f\left({\mathbf{p}}\right)pdp},$
(22)

the growth rate can be written as

${}\Gamma=\frac{{{\omega_{p}}}}{2}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right){\left({\frac{{{\omega_{p}}}}{{ck}}}\right)^{3}}\int_{{\theta_{1}}}^{{\theta_{2}}}{\frac{{\left[{\frac{1}{\theta}\frac{{\partial
g}}{{\partial\theta}}\left({\cos\theta-\frac{{ck}}{{{\omega_{p}}}}\cos{\theta_{0}}}\right)-2g}\right]}}{{\sqrt{\left({\cos{\theta_{2}}-\cos\theta}\right)\left({\cos\theta-\cos{\theta_{1}}}\right)}}}\theta
d\theta},$ (23)

where the boundaries $\theta_{1,2}$ are given by:

${}\cos{\theta_{1,2}}=\frac{{{\omega_{p}}}}{{ck}}\cos{\theta_{0}}\pm\sin{\theta_{0}}\sqrt{1-{{\left({\frac{{{\omega_{p}}}}{{ck}}}\right)}^{2}}}.$
(24)

There is a simple estimate (Breǐzman & Ryutov 1971) for the maximal growth
rate for a general wave vector angle, $\theta_{0}$:

${}\Gamma\sim\omega_{p}\left(\frac{n_{b}}{n_{p}}\right)\frac{\cos^{2}\theta_{0}}{\gamma\Delta\theta^{2}}.$
(25)

One sees that the growth rate decreases with increasing beam angular width,
$\Delta\theta$.

The expression (25) is a good estimate but the exact dependence of the growth
rate on $\theta_{0}$ is rather important since the growth of the oscillations
is exponential: a factor of two in the growth rate value will eventually lead
to difference by orders of magnitude between two waves’ energies. The exact
growth rate depends on the shape of the distribution function.

To demonstrate the dependence of the growth rate on the wave vector, we
consider a Gaussian distribution,

$g(\theta)=\frac{{{n_{b}}}}{{\gamma
mc}}\frac{2}{{\Delta{\theta^{2}}}}\exp\left({-\frac{{{\theta^{2}}}}{{\Delta{\theta^{2}}}}}\right).$
(26)

In this case, the growth rate may be found in the form (see Appendix B)

$\begin{gathered}\Gamma={\omega_{p}}\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\frac{{2\pi}}{{\gamma\Delta{\theta^{2}}}}{\cos^{3}}{\theta_{0}}\exp\left({-\frac{{\theta_{1}^{2}+\theta_{2}^{2}}}{{2\Delta{\theta^{2}}}}}\right)\times\hfill\\\
\left[{\left({{\xi^{2}}{{\left({\frac{{\Delta\theta}}{{2{\theta_{0}}}}}\right)}^{2}}-1}\right){I_{0}}\left(\xi\right)-\frac{\xi}{2}{I_{1}}\left(\xi\right)}\right],\hfill\\\
\end{gathered},$ (27)

where
$\xi\equiv\left({\theta_{2}^{2}-\theta_{1}^{2}}\right)/2\Delta{\theta^{2}}$
and $I_{0,1}\left(x\right)$ are the modified Bessel functions of the first
kind. The growth rate depends on the direction of the wave, $\theta_{0}$, and
on the absolute value of the wave vector, $k$. Figure 2 shows the maximum
growth rate with respect to $k$ as a function of $\theta_{0}$. For large
angles, $\theta_{0}\gg\Delta\theta$, the maximum growth rate varies as
$\cos^{2}\theta_{0}$, according to the general estimate (25). Because of that
dependence, the growth rate is constant in a wide range of angles,
$\Delta\theta\ll\theta_{0}\ll 1$. The shape of the maximum growth rate curve
for $\theta_{0}\sim\Delta\theta$ is sensitive to the shape of the distribution
function. Specifically for the Gaussian case, one gets a valley near
$\theta_{0}\sim\Delta\theta$ and the growth rate increases again to a plateau
at $\theta_{0}\ll\Delta\theta$. The plateau at very small angles is larger
than the intermediate angles plateau $\Delta\theta\ll\theta_{0}\ll 1$, for
this case by roughly 10%.

For every wave vector directed at the angle $\theta_{0}$, there exists a
narrow range of $k$ which is within resonance with the narrow particle beam.
One can conveniently use instead of $k$, the variable

${}x=\frac{\frac{\omega_{p}}{ck}-\cos\theta_{0}}{\sqrt{1-\frac{\omega_{p}}{ck}\cos\theta_{0}}}.$
(28)

For the remainder of this analysis, we shall use the wave vector resonance
parameter, $x$, as the wave parameter instead of $k$. The transition from $k$
to $x$ makes the various integration formulas more convenient. For large wave-
vector angles, $\theta_{0}\gg\theta$,

${}\left|x\right|\leq\theta,$ (29)

which is independent of $\theta_{0}$. For simple estimations, one can usually
assume that $\left|x\right|\sim\Delta\theta$, where $\Delta\theta$ is the
angular spread of the beam velocity distribution. Since $x$ is of the same
order as $\theta$ and $\Delta\theta$, one can safely assume that $x$ is a
small parameter, which considerably simplifies the derivation (see Appendix
C).

For $\Delta\theta\ll\theta_{0}$ and $\xi\gg 1$ equation (27) is reduced to

$\Gamma\left({{\theta_{0}}\gg\left|x\right|}\right)=\omega_{p}\sqrt{\pi}\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\left({\frac{{-x}}{{\Delta\theta}}}\right)\frac{{\exp\left({-\frac{{{x^{2}}}}{{\Delta{\theta^{2}}}}}\right)}}{{\gamma\Delta{\theta^{2}}}}{\cos^{2}}{\theta_{0}}.$
(30)

This growth rate is maximized for $x=-\frac{\Delta\theta}{\sqrt{2}}$ whose
value is:

${}{\Gamma^{{\rm{max}}}}\left({{\theta_{0}}\gg\Delta\theta}\right)={\omega_{p}}\sqrt{\frac{\pi}{{2e}}}\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\frac{{{{\cos}^{2}}{\theta_{0}}}}{{\gamma\Delta{\theta^{2}}}}.$
(31)

For the other extreme case where $\left|x\right|\gg\theta_{0}$ and $x<0$,
which corresponds to very small wave vector angles, the growth rate is

$\Gamma\left({\left|x\right|\gg{\theta_{0}}}\right)=2\pi{\omega_{p}}\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\left({\frac{{2{x^{2}}}}{{\Delta{\theta^{2}}}}-1}\right)\frac{{{e^{-\frac{{2{x^{2}}}}{{\Delta{\theta^{2}}}}}}}}{{\gamma\Delta{\theta^{2}}}},$
(32)

which is maximized for $x=-\Delta\theta$ and its value is:

${\Gamma^{{\rm{max}}}}\left({{\theta_{0}}\to
0}\right)=\frac{{2\pi}}{{{e^{2}}}}{\omega_{p}}\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\frac{1}{{\gamma\Delta{\theta^{2}}}}.$
(33)

One can show that the modes $x>0$ are stable. We therefore do not consider
positive values for $x$.

Figure 2: The maximum growth rate, $\Gamma^{\text{max}}$, as a function of
wave direction, $\theta_{0}$, for a Gaussian distribution with
$\Delta\theta=10^{-6}$. For any $\theta_{0}$, the maximum was found with
respect to the resonance parameter $x$. The growth rate is normalized with
respect to
$\Gamma_{0}=\omega_{p}\left(n_{b}/n_{p}\right)\left(\gamma\Delta\theta^{2}\right)^{-1}$.

The behavior of the maximum curve at $\theta_{0}\sim\Delta\theta$
significantly depends on the shape of the distribution function (Breǐzman &
Ryutov 1971; Rudakov 1971). In the general case, the growth rate could only be
found numerically. Taking into account a rather non-trivial behavior of the
growth rate at $\theta_{0}\sim\Delta\theta$, we could not use the same general
formula (23) and the same grid in order to find the growth rate in the entire
range of possible resonance waves. Instead, we introduce three overlapping
regions in the $(\theta_{0},x)$ plane. Within each region we expand the
general formula to leading order in the small parameters. After a long
derivation (see Appendix D) we get:

  1. 1.

$\theta_{0}\gg\left|x\right|$

$\Gamma=\frac{{{\omega_{p}}}}{2}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right){\rm{sgn}}(x){\cos^{2}}{\theta_{0}}\int_{1}^{\infty}{\frac{{\frac{{\partial
g}}{{\partial z}}}}{{\sqrt{z-1}}}dz},$ (34)

where

$z=\theta^{2}/x^{2}.$ (35)

Note that this function is anti-symmetric in $x$.

  2. 2.

$\theta_{0}\sim\left|x\right|$
For brevity one can define new parameters:

$\mu\equiv\theta_{0}^{2}/x^{2};\quad\rho\equiv{\mathop{\rm
sgn}}(x)\sqrt{1+4\mu}-1.$ (36)

The boundaries of the integral, eq. (24), are rewritten in the variable $z$ as

$z_{1,2}=2\mu-\rho\pm 2\sqrt{\mu\left(\mu-\rho\right)},$ (37)

The resonance range of $x$ guarantees that the roots $z_{1,2}$ are real. The
growth rate in this case is given by

$\Gamma={\omega_{p}}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\int_{{z_{1}}}^{{z_{2}}}{\frac{{\frac{1}{2}\frac{{\partial
g}}{{\partial
z}}\left({\rho-z}\right)-g}}{{\sqrt{\left({z-{z_{1}}}\right)\left({{z_{2}}-z}\right)}}}dz}.$
(38)

  3. 3.

$\theta_{0}\ll\left|x\right|$, $\quad x<0$
For this limit, the integration boundaries approach the same value, so the
integration effectively becomes a delta-integration at $z=2$. The growth rate
reduces to a simpler form:

$\Gamma=\pi{\omega_{p}}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right){\left[{\left({-2\frac{{\partial
g}}{{\partial z}}}\right)-g}\right]_{z=2}}.$ (39)

The three cases show that the growth rate is always proportional to the same
basic value $\Gamma_{0}$:

$\Gamma_{0}={\frac{n_{b}}{n_{p}}}\frac{1}{\gamma\Delta\theta^{2}}{\omega_{p}},$
(40)

which is consistent with the estimate given by eq. (25).

As an example of a non-Gaussian angular distribution, we calculate the
instability growth rate for the function

$g=g_{0}\exp{\left[-\left(\frac{\theta}{\Delta\theta}\right)^{5}\right]},$
(41)

which describes the self-similar expansion of the beam (Rudakov, 1971). Figure
3 shows the maximal growth rate with respect to the variable $x$ as a function
of $\theta_{0}$ for the case $\Delta\theta=10^{-6}$. One sees that for
$\theta_{0}\gg\Delta\theta$, the maximum growth rate is proportional to
$\cos^{2}\theta_{0}$, according to the general estimate (25); this part is not
sensitive to the shape of the distribution function. The effect of the beam
distribution is seen clearly at $\theta_{0}\sim\Delta\theta$. In this range,
the maximal growth rate increases by a factor of three. This is an effect of a
steeper gradient in the distribution function at $\theta\sim\Delta\theta$. For
distributions with a steep front, the growth rate for the angles
$\theta_{0}\sim\Delta\theta$ turns out to be larger than for moderate angles
$\theta_{0}\sim 1$ (Rudakov, 1971). Hence the ratio between the two ‘plateau’
values which are obtained for the maximal growth rate highly depends on the
steepness of the distribution function.

Figure 3: The same as in Fig. 2 for the beam distribution function (41) with
$\Delta\theta=10^{-6}$.

The larger growth rate at $\theta_{0}\sim\Delta\theta$ implies that the
spectrum of excited waves is dominated by the waves propagating within the
beam opening angle. These are the same waves which are capable of efficient
energy exchange with the beam. This subtle fact explains why in the kinetic
regime, relativistic beams efficiently lose their energy (Fainberg et al.,
1970; Breǐzman & Ryutov, 1971; Rudakov, 1971). However, we have seen above
that in very narrow beams, the width of resonance for these waves is
especially low. Therefore in this crucial range of angles, the instability is
easily suppressed. Below we show that no such condition exists which allows
the excitation of these waves by relativistic intergalactic beams.

##  4 Estimates for intergalactic beams

Now let us apply these general results to extragalactic beams in particular.
The IGM plasma density is (Broderick et al., 2012)

${}n_{p}=2.2\times 10^{-7}\left(1+z\right)^{3}\text{cm}^{-3},$ (42)

where $z$ is the redshift of the TeV blazar. The beam density $n_{b}$ is
determined by balance between cooling and creation processes. The upper limit
for the beam’s density is found when the cooling is determined by inverse
Compton (IC) scattering. The IC scattering rate is given by (Broderick et al.,
2012):

${\Gamma_{{\text{IC}}}}=1.4\times{10^{-20}}{\left({1+z}\right)^{4}}\gamma\quad{{\text{s}}^{-1}},$
(43)

and the corresponding upper limit for the beam density derived from the
cooling rate is:

${}n_{b}^{\text{max}}\sim 5.1\times
10^{-25}\left({1+z}\right)^{9.5}\left(\frac{L_{E}}{10^{45}\text{erg}\;\text{s}^{-1}}\right)\left(\frac{E}{\rm
TeV}\right)\text{cm}^{-3}.$ (44)

Here $L_{E}$ is the gamma-ray luminosity of the blazar at the photon energy
$E$. The Lorenz factors of the produced pairs are found in the range
$\gamma\sim 10^{4}-10^{7}$ with an average of $\gamma\sim 10^{5}$ (Miniati &
Elyiv, 2013). The pairs are produced at an angular spread
$\Delta\theta\sim\gamma^{-1}$.

Substituting the above plasma and beam densities into the condition for the
hydrodynamic instability regime, Eq. (19), and taking into account that
$\Delta\theta\sim\gamma^{-1}$, one finds that the hydrodynamic regime is valid
only for large Lorentz factor values:

$\displaystyle\gamma\gg
7.5\times{10^{7}}{\left(1+z\right)^{-13/6}}{\left({\frac{L_{E}}{{{{10}^{45}}{\text{erg}}\;{{\text{s}}^{-1}}}}}\right)^{-1/3}}$
$\displaystyle{;\quad{\theta_{0}}\sim 1},$ (45) $\displaystyle\gamma\gg
8\times{10^{5}}{\left(1+z\right)^{-13/8}}{\left({\frac{L_{E}}{{{{10}^{45}}{\text{erg}}\;{{\text{s}}^{-1}}}}}\right)^{-1/4}}$
$\displaystyle{;\quad{\theta_{0}}\sim\Delta\theta}.$

One sees that the condition for the hydrodynamic regime is not satisfied for
oblique waves and could only be marginally satisfied for quasi-parallel waves.
This means that the treatment of the beam instability in the kinetic regime is
justified for typical beams in the IGM. Now let us estimate the role of plasma
inhomogeneity on the development of the instability. It was discussed in
section 2.3 that the effects of plasma inhomogeneities on the instability
result in loss of resonance between the particles and the Langmuir wave caused
by changes in the wave vector at a rate given by (12). The resonant wave grows
exponentially with a rate $\Gamma$ from the initial thermal noise level. In
order for the wave to grow significantly, the exponent has to be rather large.
Therefore the time of growth must exceed $t=\Lambda/\Gamma$, where $\Lambda$
is the Coulomb logarithm.

It was shown in section 2.1 that the resonance width of wave numbers strongly
depends on the angle between the wave and the beam axis, $\theta_{0}$, and may
be estimated as

$\frac{\delta k}{k}\sim\left\\{\begin{array}[]{ll}\Delta\theta;&\theta_{0}\sim
1\\\ (\Delta\theta)^{2};&\theta_{0}\sim\Delta\theta.\end{array}\right.$ (46)

Now the condition for the development of the instability can be written as:

$\frac{\Lambda c}{L\Gamma}\leq\frac{\delta
k}{k}\sim\left\\{\begin{array}[]{ll}\Delta\theta;&\theta_{0}\sim 1\\\
(\Delta\theta)^{2};&\theta_{0}\sim\Delta\theta.\end{array}\right.$ (47)

Making use of the estimate (31) for the kinetic growth rate and introducing
the following short-hand notations:

$n_{b}=n_{b,-22}\times 10^{-22}\text{cm}^{-3}\>,\>n_{p}=n_{p,-7}\times
10^{-7}\text{cm}^{-3},$

we can write the condition in the form:

$\displaystyle{}0.1\left({\gamma\Delta\theta}\right)\left({\frac{\Lambda}{{20}}}\right)\left({\frac{{n_{p,-7}^{1/2}}}{{{n_{b,-22}}}}}\right)\left({\frac{{100Mpc}}{L}}\right)<1;\qquad\theta_{0}\sim
1$ (48) $\displaystyle
0.1\gamma\left({\frac{\Lambda}{{20}}}\right)\left({\frac{{n_{p,-7}^{1/2}}}{{{n_{b,-22}}}}}\right)\left({\frac{{100Mpc}}{L}}\right)<1;\qquad\theta_{0}\sim\Delta\theta.$
(49)

We normalize $L$ by the maximal inhomogeneity scale, of the order of the
cosmological void side. One can see that even in this case, only oblique waves
could grow. Interaction of the beam with these waves leads only to expansion
of the beam without noticeable loss of energy. The instability ceases when the
beam expands about an order of magnitude, $\gamma\Delta\theta\sim 10$. Nearly
parallel modes, $\theta_{0}\sim\Delta\theta$, which could take the energy off
the beam, would not grow at the relevant Lorentz factors of the beam,
$\gamma\sim 10^{4}-10^{7}$. These are the only waves which could take a
significant fraction of the beam energy, but one can see that their growth is
completely suppressed by the inhomogeneity condition. Therefore even the
weakest possible inhomogeneity makes any energy loss of the beam inefficient.

Now let us come back to the theory proposed by Shalaby et al. (2018, 2020),
which is discussed in sect. 4. Even if their theory were universlly correct,
their instability growth rate (see eq. 4.16 in Shalaby et al. 2020) would only
be

$\Gamma\sim\frac{1}{2}\left(\frac{2n_{b}}{n_{p}\gamma^{3}}\right)^{2/5}\omega_{p}=7.5\times
10^{-13}\frac{n^{0.4}_{b,-22}n^{0.1}_{p,-7}}{\gamma^{1.2}_{6}}\,\rm s^{-1},$
(50)

which is less rapid than the collisional decay time of the plasma waves (e.g.,
Miniati & Elyiv 2013)

${}\nu_{\rm coll}\approx 10^{-11}n_{p,-7}T_{4}^{-3/2}\,\rm s^{-1},$ (51)

where $T=T_{4}\times 10^{4}$ K. We stress out that their result was obtained
in the hydrodynamic approximation of the instability, even though it should be
described by kinetic theory at these parameters. Then the kinetic instability
is weaker than the hydrodynamic one at the same beam density because in the
kinetic regime, only a fraction of electrons resonantly excite the wave. This
means that the waves propagating along the beam are not excited in any case.
As to oblique waves, it was explained above and will be confirmed by the
rigorous theory in the next section that the excitation of these waves results
merely in the expansion of the beam, but doesn’t lead to energy loss. In the
next section, we demonstrate straightforwardly that even if the oblique
instability develops, the beam just expands, while energy loss remains
negligibly small.

##  5 Evolution of the beam; quantitative theory

###  5.1 Particle diffusion in momentum space

When excited waves’ energy is small in comparsion with the plasma’s thermal
energy, the evolution of the beam distribution function is governed by the
diffusion equation in momentum space (e.g. Breǐzman & Ryutov 1971; Brejzman
1990)

$\frac{\partial f}{\partial t}=\frac{\partial}{\partial
p_{\alpha}}\left(D_{\alpha\beta}\frac{\partial f}{\partial p_{\beta}}\right),$
(52)

where $D_{\alpha\beta}$ is the resonant momentum-diffusion tensor defined by

${D_{\alpha\beta}}={\pi
e^{2}}\int{W\left({{\bf{k}},t}\right)\frac{{{k_{\alpha}}{k_{\beta}}}}{{{k^{2}}}}\delta\left({{\bf{k}}\cdot{\bf{v}}-\omega}\right){d^{3}}k}.$
(53)

For the spherical, azimuthal-symmetric case, the diffusion equation is written
as (under the previous assumption that the beam is narrow, $\theta\ll 1$)

${}\begin{array}[]{l}\frac{{\partial f}}{{\partial
t}}=\frac{1}{{{p^{2}}\theta}}\frac{\partial}{{\partial\theta}}\left({\theta{D_{\theta\theta}}\frac{{\partial
f}}{{\partial\theta}}}\right)+\frac{1}{{p\theta}}\frac{\partial}{{\partial\theta}}\left({\theta{D_{\theta
p}}\frac{{\partial f}}{{\partial p}}}\right)\\\
+\frac{1}{{{p^{2}}}}\frac{\partial}{{\partial p}}\left({p{D_{\theta
p}}\frac{{\partial
f}}{{\partial\theta}}}\right)+\frac{1}{{{p^{2}}}}\frac{\partial}{{\partial
p}}\left({{p^{2}}{D_{pp}}\frac{{\partial f}}{{\partial
p}}}\right).\end{array}$ (54)

The diffusion coefficients become

${D_{\mu}}={\pi
e^{2}}\int{W\left({{\bf{k}},t}\right)\left(\frac{k_{\theta}}{k}\right)^{\mu}\delta\left({{\bf{k}}\cdot{\bf{v}}-\omega}\right){d^{3}}k},$
(55)

where $\mu=0,1,2$ for $D_{pp},D_{\theta p},D_{\theta\theta}$, respectively.

It follows immediately from eq. (55), that if the spectral energy of
oscillations, $W\left({\bf{k}}\right)$, is dominated by waves propagating
within the angle $\theta_{0}$, the diffusion coefficients are related as

${D_{\theta\theta}}:{D_{\theta p}}:{D_{pp}}\sim\theta_{0}^{2}:\theta_{0}:1.$
(56)

The intergalactic beams are extremely narrow, $\Delta\theta\sim 10^{-5}$,
whereas their energy spread is large $\Delta p\sim p$. In this case, the
different terms in the R.H.S. of eq. (54) are related as

$\begin{gathered}\frac{1}{\theta}\frac{\partial}{{\partial\theta}}\left({\theta{D_{\theta\theta}}\frac{{\partial
f}}{{\partial\theta}}}\right):\frac{p}{\theta}\frac{\partial}{{\partial\theta}}\left({\theta{D_{\theta
p}}\frac{{\partial f}}{{\partial p}}}\right)\hfill\\\
:\frac{\partial}{{\partial p}}\left({p{D_{\theta p}}\frac{{\partial
f}}{{\partial\theta}}}\right):\frac{\partial}{{\partial
p}}\left({{p^{2}}{D_{pp}}\frac{{\partial f}}{{\partial p}}}\right)\hfill\\\
\sim{\left({\frac{{{\theta_{0}}}}{{\Delta\theta}}}\right)^{2}}:\left({\frac{{{\theta_{0}}}}{{\Delta\theta}}}\right):\left({\frac{{{\theta_{0}}}}{{\Delta\theta}}}\right):1\hfill\\\
\end{gathered}$ (57)

One sees that if the energy of the excited waves is dominated by oblique
waves, $\theta_{0}\gg\Delta\theta$, the first term in the R.H.S. of eq. (54)
significantly exceeds the others, which means that the beam expands without
losing energy. Only if the waves are excited predominantly within the beam
opening angle, $\theta_{0}\sim\Delta\theta$, all the terms on the R.H.S. are
comparable, and then the beam’s energy is efficiently redistributed.

A special feature of the kinetic beam instability is that at
$\theta_{0}<\Delta\theta$, the growth rate is larger than at a larger
$\theta_{0}\sim\Delta\theta$ (Breǐzman & Ryutov 1971; Rudakov 1971, see also
section 3.2). Due to the exponential growth from a very low level of initial
thermal fluctuations, nearly parallel waves could reach a much larger energy
than oblique ones. This subtle fact leads to the widely known statement
(Fainberg et al., 1970; Breǐzman & Ryutov, 1971; Rudakov, 1971) that in the
kinetic regime, relativistic beams efficiently lose their energy. However, the
width of resonance for these waves is extremely small. Therefore the
instability of these waves may be easily suppressed by, e.g., weak
inhomogeneity (Breizman & Ryutov 1970, see also section 2.2). The estimates
show (see section 4) that for typical parameters of intergalactic beams, this
effect completely suppresses the excitation of nearly parallel waves, so the
beam could only expand without cooling. In section 5.3 we present numerical
simulations supporting the estimates above.

###  5.2 Evolution of the wave spectrum

The evolution equation for the oscillation spectrum density,
$W\left(\bf{k}\right)$, is written in the quasilinear approximation as

$\frac{\partial W}{{\partial t}}+{\bf v}_{g}\cdot\nabla
W-\nabla\omega\cdot\frac{{\partial
W}}{{\partial{\bf{k}}}}=2\left(\Gamma-\nu_{\rm coll}\right)W,$ (58)

where $\nu_{\rm coll}$ is the collisional decay rate of the plasma waves (51).
The l.h.s. of this equation is the full time derivative of the plasma wave
spectral energy density along the rays which obey the geometrical optics law:

$\frac{d\bf
r}{dt}=\frac{{\partial\omega}}{{\partial{\bf{k}}}};\qquad\frac{d\bf
k}{dt}=-\nabla\omega.$ (59)

The first term in (58) represents the time evolution, the second is spatial
advection, and the third term describes a spectral shift in $k$-space due to
inhomogeneity. The r.h.s. of this equation describes the growth or damping of
waves’ energy due to resonance interaction with the beam and decay due to
Coulomb collisions (see Eq. (51)).

We do not take into account non-linear processes such as induced scattering
and modulation instability of the plasma waves. These processes redistribute
the energy of the waves over the spectrum, thus removing them from the
resonance region. Therefore they provide additional mechanisms for the
stabilization of the beam instability. We show that the inhomogeneity solely
restricts the efficiency of energy loss via the beam instability. Thus, there
is no special need to consider additional processes.

A formal solution to eq. (58) is

$W(t,{\bf r},{\bf
k})=W_{0}\exp\left\\{2\int_{0}^{t}\left[\Gamma(t^{\prime},{\bf
r}^{\prime},{\bf k}^{\prime})-\nu_{\rm coll}({\bf
r}^{\prime})\right]dt^{\prime}\right\\},$ (60)

where the integral is along the characteristics ${\bf r}^{\prime}={\bf r}+{\bf
v}_{g}(t^{\prime}-t)$, ${\bf k}^{\prime}={\bf k}-\nabla\omega(t^{\prime}-t)$.
The growth rate $\Gamma$ is given by eq. (23). The time-dependence of the
growth rate, $\Gamma$, is determined by the time variation of the distribution
function. The dependence on $\bf k$ is very strong because the resonance range
$\Delta k$ is very narrow. Both $\Gamma$ and $\nu_{\rm coll}$ depend on the
spatial coordinate only via the plasma density, which enters the expressions
for $\Gamma$ and $\nu_{\rm coll}$ as a multiplier and varies at a very large
scale $L$. Taking into account that the light travel time across the spatial
inhomogeneity scale significantly exceeds the instability time, $t\ll L/c$,
one can neglect this dependence in the integrand. Assuming that the initial
fluctuation rate, $W_{0}$, is spatially homogeneous, one finds that the
dependence on $\bf r$, and therefore on ${\bf v}_{g}$, drops out of the
solution. This means that one can drop the second term in the l.h.s. of eq.
(58): the evolution of the wave density is determined by the time evolution of
$\bf k$ according to the second equation in (59) and not by the spatial
transfer.

For numerical purposes, we will use eq. (58) but not eq. (60).
Now we can rewrite this equation in a simplified form:

$\frac{{\partial W}}{{\partial t}}+\frac{\omega_{p}}{L}\frac{{\partial
W}}{{\partial{k}}}=2\left(\Gamma-\nu_{\rm coll}\right)W.$ (61)

Taking into account that the waves are concentrated in a very narrow range of
$k$, one can conveniently rewrite, for numerical purposes, the equation in
terms of the resonant parameter (29), which is roughly bounded by
$|x|\leq\Delta\theta$. Then for the case $\theta_{0}\gg\Delta\theta$, we get

${}\frac{{\partial W}}{{\partial
t}}-\frac{c}{L}\frac{{{{\cos}^{2}}{\theta_{0}}}}{{\sin{\theta_{0}}}}\frac{{\partial
W}}{{\partial x}}=2\left(\Gamma-\nu_{\rm coll}\right)W;$ (62)

whereas for $\theta_{0}\sim\Delta\theta$, the equation may be written as

$\frac{{\partial W}}{{\partial
t}}-\frac{2}{{{\theta^{2}}}}\frac{c}{L}\frac{{\partial W}}{{\partial
y}}=2\left(\Gamma-\nu_{\rm coll}\right)W,$ (63)

where

$y=\frac{{x\sqrt{{x^{2}}+4\theta_{0}^{2}}-{x^{2}}}}{{{\theta^{2}}}}+1.$ (64)

One can show that the resonance width $\Delta y$ is roughly
$\theta_{0}/\theta$ so that $\Delta y\sim 1$ for $x\sim\theta$.

###  5.3 Numerical simulation

In this section we present the results of our numerical simulation of the
beam-plasma interaction in the quasi-linear regime.

The extragalactic beams are extremely narrow, $\Delta\theta\sim 10^{-5}$. Our
previous estimates show that only oblique waves, $\theta_{0}\sim 1$, are
excited, in the presence of inhomogeneity. This implies that the first term on
the R.H.S. of the momentum diffusion equation (54) exceeds other terms by 5-10
orders of magnitude. Therefore we initially only consider the evolution of the
angular distribution of particles according to the equation

$\frac{{\partial g}}{{\partial t}}=\frac{1}{{{{\left({\gamma
mc}\right)}^{2}}}}\frac{\partial}{{\partial{\theta^{2}}}}\left({4{\theta^{2}}{D_{\theta\theta}}\frac{{\partial
g}}{{\partial{\theta^{2}}}}}\right).$ (65)

where the angular distribution function, $g$, has been defined in (22). Having
found the angular evolution of the beam, together with the spectral evolution
of the wave, we evaluate the other diffusion coefficients in the momentum
diffusion equation. If they are comparable to $D_{\theta\theta}$, then
according to our estimates the other diffusion processes are negligible, which
is self-consistent with our assumption.

For each time step, we first solve the diffusion equation (65) given the
diffusion coefficients from the last time step. After the distribution
function is updated, we calculate the growth rates for each wave vector,
making use of the formulas presented in section 3.2. After the growth rates
are calculated in the current time step, we next calculate the new values of
the oscillations energy density, $W$, according to the evolution equation (eq.
62 or, for small $\theta_{0}$, 63). Using the updated values of $W$, we
calculate the diffusion coefficients (the procedure is described in Appendix
E). The new values of the diffusion coefficients are used in the next time
step for a new cycle as described above.

Each separate simulation is characterized by a set of
$\left(\gamma,n_{p},n_{b},L\right)$. The simulation time is measured in units
of the characteristic instability time

$t_{0}=\Gamma_{0}^{-1}=\left(n_{p}/n_{b}\right)\omega_{p}^{-1}\gamma^{-1}.$
(66)

The initial angular distribution function is chosen in the form

${}g\left(\theta\right)=\exp\left(-0.2\left(\theta\gamma\right)^{5}\right),$
(67)

so that the initial angular spread is always $\Delta\theta\approx\gamma^{-1}$,
and the function has a sharp gradient around $\theta\sim\gamma^{-1}$. The
oscillations energy density is initialized uniformly to a low level such that
the angular diffusion rate due to the initial background is negligibly slow
(see below). In all of the simulations we performed, we have fixed the value
of the plasma density to be $n_{p}=2.2\times 10^{-7}\text{cm}^{-3}$, which
corresponds to no redshift. We consider the maximal possible plasma
inhomogeneity length scale, $L=100Mpc$, providing the minimal influence of
inhomogeneity.

Figure 4 shows the angular spread, $\Delta\theta$, and the maximum value of
the angular distribution function, ${\rm max}\left[g\right]$, as a function of
the normalized time. The beam is spreading in angle while the distribution
flattens, which means that particles gain transverse momentum and they are
drifting away from the initial direction of the beam. The spreading begins
after about a dozen instability times, when the resonant wave grew enough from
the initial low level. After the beam spreads by about an order of magnitude,
the growth rate of the instability falls below the collision decay rate of the
plasma waves, and the process is saturated.

In figure 5 we present simulations of the beam with the same parameters but
with different initial plasma wave energy levels, $\tilde{W}_{0}$, normalized
by

${}\tilde{W}_{0}=\frac{mc^{2}}{\omega_{p}t_{0}}\left(\frac{e^{2}}{mc^{2}}\frac{\omega_{p}}{c}\right)^{-1}W_{0},$
(68)

where $W_{0}$ is dimensionless. One sees that the results are weakly dependent
on $W_{0}$, the initial level of oscillations, provided it is small enough.

Figure 4: The angular spread of the beam distribution function (solid) as a
function of normalized time and $\max\left(g\right)$ (dashed) as a function of
normalized time. $\gamma=10^{6}$, $n_{b,-22}=1$, $L=100Mpc$. Figure 5: The
normalized angular spread, $\gamma\Delta\theta$, as a function of normalized
time for different initial levels of oscillations, $W_{0}$, see eq. (68). The
beam parameters are the same as in the previous figure.  Figure 6: The
normalized angular spread, $\gamma\Delta\theta$ as a function of normalized
time for $\gamma=10^{6},n_{b,-22}=1$ (solid line), $\gamma=10^{5},n_{b,-22}=1$
(dashed line), $\gamma=10^{6},n_{b,-22}=0.1$ (dotted line),
$\gamma=10^{5},n_{b,-22}=0.1$ (dashed-dotted line). $L=100Mpc$.

Figure 6 shows the evolution of the beam width for different values of beam
density and Lorentz factors. One sees that in all cases, the behavior of the
beam is essentially the same: the expansion saturates after the beam expands a
few times over. This agrees with the result of our qualitative estimates
presented in section 4.

In our numerical simulation, we calculate the growth and decay of the resonant
Langmuir waves energy. Figures 7 and 8 show the maximum oscillation energy
density with respect to $k$ as a function of $\theta_{0}$ at different times.
We measure the oscillations energy density with respect to their initial
value, $W_{0}$ which is considered as the dimensionless noise level. The first
figure shows results for $n_{b}=10^{-22}\text{cm}^{-3}$ and the second for
$n_{b}=10^{-23}\text{cm}^{-3}$. One sees that each curve has a maximum at some
angle $\theta_{0}\sim 1$. This peak is a result of the competition between two
terms: the growth, which is proportional to $\cos^{2}\theta_{0}$, thus
stronger for smaller angles, and the inhomogeneity term, which suppresses the
instability more efficiently for smaller angles. Both cases show that small
angles waves do not grow, as expected. Note that whereas the densities of two
beams differ by a factor of 10, the level of the oscillations energy differs
by four orders of magnitude. The reason is that the instability growth rate
for the smaller density beam is ten times smaller than that of the denser
beam. Therefore for the more dilute beam, inhomogeneity also succeeds in
suppressing the growth of oblique waves, $\theta_{0}\sim 1$ (see below).

Figure 7: The maximum oscillation energy density $W_{\rm osc}$ with respect to
the resonance parameter $x$ as a function of $\theta_{0}$ for different times
(normalized by the instability time (66)). $n_{b,-22}=1$ ,$\gamma=10^{6}$,
$L=100Mpc$. Figure 8: Same as figure 7, but with $n_{b,-22}=10^{-1}$.

It is interesting as well to examine the most energetic mode as a function of
time and how it depends on the beam parameters. This is shown in figure 9. We
have already seen in previous plots that the maximum value of $W_{\rm osc}$
initially grows but eventually begins to decrease. This occurs when
inhomogeneity overcomes the linear growth, yet the competition between the two
processes is dynamic and eventually the instability is completely suppressed.
The maximum energy density with respect to the whole resonant spectrum
strongly depends on the beam density, where a change in $n_{b}$ by a factor of
10 changes the possible growth of the energy density by orders of magnitude.

Figure 9: The maximum oscillation energy density $W_{osc}$, normalized by
$W_{0}$ with respect to the entire resonant spectrum as a function of
normalized time (66) for $\gamma=10^{6},n_{b,-22}=1$ (solid line),
$\gamma=10^{5},n_{b,-22}=1$ (dashed line), $\gamma=10^{6},n_{b,-22}=0.1$
(dotted line), $\gamma=10^{5},n_{b,-22}=0.1$ (dashed-dotted line). $L=100Mpc$.

In the above simulations, the beam expands while the instability growth rate
remains larger than the collisional decay time, after which the collisionless
evolution ceases. We also examine an artificial case where collisional decay
is absent, i.e. $\nu_{\text{coll}}=0$. For this case, waves will decay only if
they exit the resonance range due to inhomogeneity. The waves propagating
nearly perpendicularly to the direction of the density gradient still grow so
that the beam expands further out. The instability growth rate, and therefore
the rate of evolution, decreases with the beam width. Eventually the beam
expansion time becomes smaller than the Compton scattering time, after which
the collisionless relaxation does not affect the beam any more: the beam will
cool down and deposit its energy to the background light (Broderick et al.
2012; Sironi & Giannios 2014). Figure 10 shows the expansion of the beam for
this artificial case, for which $\gamma=10^{6}$ and $n_{b,-22}=0.1$. The
simulations were stopped when the expansion time became equal to the
characteristic Compton time. One sees that the beams expands a few times more
than in the case shown in figure 6. In any case, the beam evolution is
determined by obliquely propagating waves, therefore the beam practically does
not loose energy but just expands until the Compton scattering comes into the
play.

Figure 10: The evolution of the beam width when the collisional decay of the
plasma waves is absent. Shown is the normalized angular spread,
$\gamma\Delta\theta$, as a function of normalized time for
$\gamma=10^{6},n_{b,-22}=0.1$, $L=100$ Mpc. The typical inverse-Compton
scattering time is shown as a vertical black line.

###  5.4 Justification of the neglect of the energy evolution of the beam

In the previous subsection, we studied the angular expansion of the beam. In
addition to angular diffusion, the beam particles could exchange energy with
the excited oscillations, which leads to a radial diffusion in momentum space.
The energy exchange is governed by the small angles part of the spectrum of
oscillations. In section 5.1 we have estimated the relative rate of the
angular and energy diffusion and have shown that it depends both on the ratio
between the diffusion coefficients and on the beam width, so that if the
diffusion coefficients are of the same order of magnitude, then energy
diffusion is negligible with respect to angular diffusion, see eq. (57).
According to our estimates, the diffusion coefficients are of the same order
unless the plasma wave spectrum is dominated by small angle waves. Taking into
account that growth of these waves is totally suppressed even by the slightest
possible inhomogeneity, we retained only the angular terms in the momentum
diffusion equation. Now we could justify our conjecture by directly
calculating all the diffusion coefficients making use of the spectrum as
obtained in simulations.

Figure 11 shows the angular dependence of the diffusion coefficients for
different times whereas figure 12 shows the time evolution of the diffusion
coefficients for different values of $\theta$. One sees that the diffusion
coefficients remain of the same order all the time for all the relevant
particle angles. This shows that the terms on the R.H.S. of equation (54) are
related as $1:\Delta\theta:\Delta\theta:(\Delta\theta)^{2}$, correspondingly.
Taking into account that $\Delta\theta\sim 10^{-6}-10^{-5}$, this justifies
the neglect of the energy evolution of the beam.

Figure 11: Angular dependence of the diffusion coefficients for different
times: $t=103.34t_{0}$ (squares); $t=147.4t_{0}$ (circles); $t=213.16t_{0}$
(dashes); $t=497.4t_{0}$ (triangles). $n_{b,-22}=1$ ,$\gamma=10^{6}$,
$L=100Mpc$. Figure 12: The time evolution of the diffusion coefficients at
$\gamma\theta=0.1$ (solid line), $\gamma\theta=1$ (dashed) and
$\gamma\theta=10$ (dotted).Beam parameters are the same as in the previous
figure.

##  6 Conclusions

In this paper, we have examined the beam-plasma instability of Langmuir waves
in the intergalactic medium. The source of the instability is the resonant
energy exchange between a narrow, dilute, energetic beam of ultra-relativistic
pairs traveling through the plasma and the Langmuir oscillations of the
background plasma. The beams of electrons and positrons which travel through
the IGM are created by pair production processes of highly energetic gamma ray
photons with the background radiation in the plasma. The evolution of the pair
beams and the deposition of their energy to the IGM and back to the gamma ray
flux is important to understand in order to explain the lack of any GeV halos
around blazars (Aharonian et al., 2006; Neronov & Vovk, 2010; Broderick et
al., 2012).

The crucial point is that the energetic evolution of the beam, whether it
loses energy or gains transverse momentum, depends on the direction of the
wave vector: if the interaction is dominated by waves propagating within the
opening angle of the beam, $\theta_{0}\sim\Delta\theta$, then the beam loses a
large fraction of its energy. However, when the interaction is dominated by
oblique waves, the momentum loss is small and the beam simply expands.

For realistic parameters of the IGM, the beam instability develops in the
kinetic regime. The kinetic growth rate dependence on the wave vector angle
shows that the fastest growing waves are ones which travel at small angles
($\theta_{0}\sim\Delta\theta$); their growth rate is a few times larger than
that of oblique waves ($\theta_{0}\sim 1$). This property of the kinetic
instability explains the claim in the literature that the beam efficiently
loses its energy (Fainberg et al., 1970; Breǐzman & Ryutov, 1971; Rudakov,
1971) and therefore the instability is proposed as the principal mechanism of
the beam evolution. However, the width of the wave-beam resonance for the
nearly parallel waves is small $(\sim(\Delta\theta)^{2})$, which means that in
a very narrow beam, the growth of these waves is easily suppressed. We have
shown that even for the slightest density gradient with a length scale as
large as 100 Mpc, the quasi-parallel waves lose resonance at any reasonable
parameters of intergalactic beams. Therefore the growth of the same waves
which contribute to the decelerating mechanism is forbidden by the
characteristics of the inhomogenous medium.

We have solved the quasi-linear equations for the momentum diffusion of pair
beams with a time evolution of the plasma oscillations, including the
instability mechanism and the inhomogeneity term. We have found that the beam
broadens in angle by about an order of magnitude but practically does not lose
energy. This leads to the conclusion that the lack of GeV bumps in the blazar
spectra cannot be attributed to the beam decay due to plasma instabilities.

## Acknowledgements

This research was supported by the grant 2067/19 from the Israeli Science
Foundation.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Ackermann et al. (2018) Ackermann M., et al., 2018, The Astrophysical Journal Supplement Series, 237, 32
  * Aharonian et al. (2006) Aharonian F., et al., 2006, Nature, 440, 1018
  * Alves Batista et al. (2019) Alves Batista R., Saveliev A., de Gouveia Dal Pino E. M., 2019, MNRAS, 489, 3836
  * Archambault et al. (2017) Archambault S., et al., 2017, The Astrophysical Journal, 835, 288
  * Bludman et al. (1960) Bludman S., Watson K., Rosenbluth M., 1960, The Physics of Fluids, 3, 747
  * Boyd et al. (2003) Boyd T., Boyd T., Sanderson J., 2003, The physics of plasmas. Cambridge University Press
  * Breizman & Ryutov (1970) Breizman B. N., Ryutov D. D., 1970, Soviet Journal of Experimental and Theoretical Physics Letters, 11, 421
  * Brejzman (1990) Brejzman B., 1990, in , Vol. 15, Reviews of Plasma Physics, ed.: B.Kadomtsev. Springer, p. 61
  * Bret et al. (2010) Bret A., Gremillet L., Dieckmann M. E., 2010, Physics of Plasmas, 17, 120501
  * Breǐzman & Ryutov (1971) Breǐzman B. N., Ryutov D. D., 1971, Soviet Journal of Experimental and Theoretical Physics, 33, 220
  * Broderick et al. (2012) Broderick A. E., Chang P., Pfrommer C., 2012, The Astrophysical Journal, 752, 22
  * Broderick et al. (2014) Broderick A. E., Pfrommer C., Puchwein E., Chang P., 2014, The Astrophysical Journal, 790, 137
  * Broderick et al. (2018) Broderick A. E., Tiede P., Chang P., Lamberts A., Pfrommer C., Puchwein E., Shalaby M., Werhahn M., 2018, ApJ, 868, 87
  * Chang et al. (2014) Chang P., Broderick A. E., Pfrommer C., Puchwein E., Lamberts A., Shalaby M., 2014, The Astrophysical Journal, 797, 110
  * Chang et al. (2016) Chang P., Broderick A. E., Pfrommer C., Puchwein E., Lamberts A., Shalaby M., 2016, The Astrophysical Journal, 833, 12pp
  * Fainberg et al. (1970) Fainberg Y. B., Shapiro V., Shevchenko V., 1970, Soviet Phys. JETP, 30, 528
  * Kempf et al. (2016) Kempf A., Kilian P., Spanier F., 2016, A&A, 585, A132
  * Krakau & Schlickeiser (2014) Krakau S., Schlickeiser R., 2014, ApJ, 789, 84
  * Krall & Trivelpiece (1973) Krall N., Trivelpiece A., 1973, Principles of plasma physics. No. 0-911351 in International series in pure and applied physics, McGraw-Hill, https://books.google.co.il/books?id=b0BRAAAAMAAJ
  * Lifshitz & Pitaevskii (1981) Lifshitz E. M., Pitaevskii L. P., 1981, Physical kinetics. Pergamon Press
  * Miniati & Elyiv (2013) Miniati F., Elyiv A., 2013, The Astrophysical Journal, 770, 54
  * Neronov & Vovk (2010) Neronov A., Vovk I., 2010, Science, 328, 73
  * Rafighi et al. (2017) Rafighi I., Vafin S., Pohl M., Niemiec J., 2017, Astronomy & Astrophysics, 607, A112
  * Rudakov (1971) Rudakov L. I., 1971, Soviet Journal of Experimental and Theoretical Physics, 32, 1134
  * Schlickeiser et al. (2012) Schlickeiser R., Ibscher D., Supsar M., 2012, The Astrophysical Journal, 758, 102
  * Schlickeiser et al. (2013) Schlickeiser R., Krakau S., Supsar M., 2013, The Astrophysical Journal, 777, 49
  * Shalaby et al. (2018) Shalaby M., Broderick A. E., Chang P., Pfrommer C., Lamberts A., Puchwein E., 2018, The Astrophysical Journal, 859, 45
  * Shalaby et al. (2020) Shalaby M., Broderick A. E., Chang P., Pfrommer C., Puchwein E., Lamberts A., 2020, Journal of Plasma Physics, 86, 535860201
  * Sironi & Giannios (2014) Sironi L., Giannios D., 2014, The Astrophysical Journal, 787, 49
  * Supsar & Schlickeiser (2014) Supsar M., Schlickeiser R., 2014, The Astrophysical Journal, 783, 96
  * Tiede et al. (2017) Tiede P., Broderick A. E., Shalaby M., Pfrommer C., Puchwein E., Chang P., Lamberts A., 2017, The Astrophysical Journal, 850, 157
  * Vafin et al. (2018) Vafin S., Rafighi I., Pohl M., Niemiec J., 2018, The Astrophysical Journal, 857, 43
  * Vafin et al. (2019) Vafin S., Deka P. J., Pohl M., Bohdan A., 2019, ApJ, 873, 10
  * Watson et al. (1960) Watson K., Bludman S., Rosenbluth M., 1960, The Physics of Fluids, 3, 741
  * Yan et al. (2019) Yan D., Zhou J., Zhang P., Zhu Q., Wang J., 2019, ApJ, 870, 17

##  Appendix A Hydrodynamic instability

The dispersion equation (3.1) for a mono-energetic beam is easily solved if
$\omega_{b}\ll\omega_{p}$ (e.g., §61 in Lifshitz & Pitaevskii 1981). One can
conveniently write the dispersion equation in the form

${}1-\frac{\omega_{p}^{2}}{\omega^{2}}=\frac{\eta}{(\omega-\mathbf{k\cdot
v}_{b})^{2}};$ (69)

where

$\eta=\omega_{b}^{2}\frac{\gamma^{2}\sin\theta_{0}+\cos\theta_{0}}{\gamma^{3}}.$
(70)

At $\omega_{b}\ll\omega_{p}$, the rhs is non-negligible only if the
denominator is close to zero, therefore one can present the solution in the
form

$\omega=\mathbf{k\cdot v}_{b}+\delta,$ (71)

where $\delta\ll\omega$. Substituting this into (69), one finds

$\delta=\pm\sqrt{\frac{\eta}{1-\frac{\omega^{2}_{p}}{k^{2}v^{2}_{b}}}}.$ (72)

This solution represents the well known beam mode. The wave is unstable at
$kv<\omega_{p}$; the growth rate is typically

$\Gamma_{1}\sim\eta^{1/2},$ (73)

but it increases when $kv_{b}$ approaches $\omega_{p}$. The maximal growth
rate is achieved at $kv\approx\omega\approx\omega_{p}$. In this case, we deal
with the plasma wave in resonance with the beam. The solution is now presented
as $\omega=\omega_{p}+\xi$, where $\xi\ll\omega_{p}$. The standard procedure
yields the growth rate of the plasma wave

$\Gamma_{2}\sim(\eta\omega_{p})^{1/3}.$ (74)

The ratio of the growth rate of the beam mode and of the plasma wave is

${}\frac{\Gamma_{1}}{\Gamma_{2}}\sim\frac{\eta^{1/6}}{\omega_{p}^{1/3}}=\left(\frac{n_{b}}{n_{p}}\frac{\gamma^{2}\sin\theta_{0}+\cos\theta_{0}}{\gamma_{b}^{3}}\right)^{1/6}.$
(75)

Due to the power of $1/6$, the growth rates of two modes could be roughly of
the same order for not very large Lorentz factors and moderate density ratios.
Therefore in simulations, one can observe the beam modes growing together with
the resonant plasma wave (e.g., Bret et al. 2010). However, in the
extragalactic beams, $n_{b}/n_{p}\sim 10^{-15}-10^{-18}$, $\gamma_{b}\sim
10^{6}$; therefore the above ratio (75) is as small as $10^{-4}-10^{-6}$.

The instability with such a small growth rate as that of the beam modes could
not develop in the IGM because of the collisional damping. Regardless f this
reason, one can readily see that the beam modes do not in fact exist in the
IGM. Namely, the condition for these modes to exist is that the beam may be
considered as monoenergetic, so that all the beam particles are in resonance
with the wave,

$\Gamma_{1}>k\Delta v.$ (76)

In our case, $\Delta v\sim c\Delta\theta$, $k\sim\omega_{p}/c$. Even for
oblique waves, for which the growth rate is larger than for parallel ones, we
get

$\frac{\Gamma_{1}}{\omega\Delta\theta}\sim\sqrt{\frac{n_{b}}{n_{p}\gamma_{b}}}\frac{1}{\Delta\theta}\sim
10^{-6}-10^{-4},$ (77)

so that the condition (76) is violated by large margin.

The above consideration shows that only resonant plasma waves could be excited
in the extragalactic beam-plasma systems therefore we consider only these
waves in the paper.

##  Appendix B Kinetic growth rate for a Gaussian distribution

Substituting into the general formula for the growth rate (23) the
distribution function (26) and taking into account that $\theta\ll 1$, one
gets:

$\begin{gathered}\Gamma=8{\omega_{p}}\left({\frac{1}{{\gamma\Delta{\theta^{4}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right){\left({\frac{{{\omega_{p}}}}{{ck}}}\right)^{3}}\times\hfill\\\
\int_{\theta_{1}^{2}}^{\theta_{2}^{2}}{d{\theta^{2}}\frac{{\frac{{{\theta^{2}}}}{2}+\left({\frac{{2{{\cos}^{2}}{\theta_{0}}}}{{\theta_{2}^{2}+\theta_{1}^{2}}}-1-\Delta{\theta^{2}}}\right)}}{{\sqrt{\left({{\theta^{2}}-\theta_{1}^{2}}\right)\left({\theta_{2}^{2}-{\theta^{2}}}\right)}}}{e^{-\frac{{{\theta^{2}}}}{{\Delta{\theta^{2}}}}}}}\hfill\\\
\end{gathered}.$ (78)

Using the integral formulas of the modified Bessel functions of the first kind
$I_{n}(x)$,

${I_{n}}\left(x\right)=\frac{1}{\pi}\int_{0}^{\pi}{{e^{x\cos\theta}}\cos\left({n\theta}\right)d\theta},$
(79)

it is possible to evaluate the integrals

$\displaystyle\begin{split}&\int_{{x_{1}}}^{{x_{2}}}{\frac{{{e^{\frac{{2x}}{{\Delta{\theta^{2}}}}}}}}{{\sqrt{\left({x-{x_{1}}}\right)\left({{x_{2}}-x}\right)}}}dx}=\pi{e^{\frac{{{x_{2}}+{x_{1}}}}{{\Delta{\theta^{2}}}}}}{I_{0}}\left({\frac{{{x_{2}}-{x_{1}}}}{{\Delta{\theta^{2}}}}}\right)\\\
\end{split},$ (80a)
$\displaystyle\begin{split}&\int_{{x_{1}}}^{{x_{2}}}{\frac{{x{e^{\frac{{2x}}{{\Delta{\theta^{2}}}}}}}}{{\sqrt{\left({x-{x_{1}}}\right)\left({{x_{2}}-x}\right)}}}dx}=\pi{e^{\frac{{{x_{2}}+{x_{1}}}}{{\Delta{\theta^{2}}}}}}\times\\\
&\left[{\frac{{{x_{2}}-{x_{1}}}}{2}{I_{1}}\left({\frac{{{x_{2}}-{x_{1}}}}{{\Delta{\theta^{2}}}}}\right)+\frac{{{x_{2}}+{x_{1}}}}{2}{I_{0}}\left({\frac{{{x_{2}}-{x_{1}}}}{{\Delta{\theta^{2}}}}}\right)}\right].\end{split}$
(80b)

Finally, we now get

$\begin{gathered}\Gamma={\omega_{p}}\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\frac{{2\pi}}{{\gamma\Delta{\theta^{2}}}}{\cos^{3}}{\theta_{0}}\exp\left({-\frac{{\theta_{1}^{2}+\theta_{2}^{2}}}{{2\Delta{\theta^{2}}}}}\right)\times\hfill\\\
\left[{\left({{\xi^{2}}{{\left({\frac{{\Delta\theta}}{{2{\theta_{0}}}}}\right)}^{2}}-1}\right){I_{0}}\left(\xi\right)-\frac{\xi}{2}{I_{1}}\left(\xi\right)}\right],\hfill\\\
\end{gathered},$ (81)

where
$\xi\equiv\left({\theta_{2}^{2}-\theta_{1}^{2}}\right)/2\Delta{\theta^{2}}$.

##  Appendix C Resonance parameter

The resonance parameter, $x$, is defined for each wave vector $(k,\theta_{0})$
as

${}x=\frac{\frac{\omega_{p}}{ck}-\cos\theta_{0}}{\sqrt{1-\frac{\omega_{p}}{ck}\cos\theta_{0}}}.$
(82)

First, we prove that $\left|x\right|\leq\theta$ for resonant waves with
particles at a polar angle $\theta$. We begin with the resonance condition

$-1\leq\frac{\omega_{p}/kc-\cos\theta\cos\theta_{0}}{\sin\theta\sin\theta_{0}}\leq
1.$ (83)

Given that the polar angles are in the range $[0,\pi]$, one can get

${}-\theta\sin{\theta_{0}}-\cos{\theta_{0}}\frac{{{\theta^{2}}}}{2}\leq\frac{{{\omega_{p}}}}{{kc}}-\cos{\theta_{0}}\leq\theta\sin{\theta_{0}}-\cos{\theta_{0}}\frac{{{\theta^{2}}}}{2}.$
(84)

After some arithmetics, we find the complete resonance range of $x$ as

$-\frac{{\theta\left({1+\frac{\theta}{{2\tan{\theta_{0}}}}}\right)}}{{\sqrt{1+\frac{\theta}{{\tan{\theta_{0}}}}+\frac{{{\theta^{2}}}}{{2{{\tan}^{2}}{\theta_{0}}}}}}}\leq
x\leq\left\\{{\begin{array}[]{*{20}{c}}{\frac{{\theta\left({1-\frac{\theta}{{2\tan{\theta_{0}}}}}\right)}}{{\sqrt{1-\frac{\theta}{{\tan{\theta_{0}}}}+\frac{{{\theta^{2}}}}{{2{{\tan}^{2}}{\theta_{0}}}}}}}}&{{\theta_{0}}\geq\theta}\\\
{{\theta_{0}}/\sqrt{2}}&{{\theta_{0}}\leq\theta}\\\ \end{array}}\right..$ (85)

In order to express $k$ via $x$ and $\theta_{0}$, we write eq. (82) as a
quadratic equation:

${\left({\frac{{{\omega_{p}}}}{{ck}}}\right)^{2}}+\left({{x^{2}}-2}\right)\cos{\theta_{0}}\frac{{{\omega_{p}}}}{{ck}}+{\cos^{2}}{\theta_{0}}-{x^{2}}=0.$
(86)

The solution is

$\frac{{{\omega_{p}}}}{{ck}}=\left({1-\frac{{{x^{2}}}}{2}}\right)\cos{\theta_{0}}+\frac{x}{2}\sqrt{{x^{2}}{{\cos}^{2}}{\theta_{0}}+4{{\sin}^{2}}{\theta_{0}}},$
(87)

where we kept only the solution which is consistent with (82). We divide the
form of this solution to three different cases, or ranges of $x/\theta_{0}$:

$\displaystyle\frac{{{\omega_{p}}}}{{ck}}=\cos{\theta_{0}}\left(1+x\tan{\theta_{0}}\right)\quad;\quad{\theta_{0}}\gg\left|x\right|,$
(88)
$\displaystyle\frac{{\omega_{p}}}{ck}=1-\frac{\theta_{0}^{2}}{2}+\frac{x^{2}}{2}\left(\text{sgn}\left({x}\right)\sqrt{1+\frac{4\theta_{0}^{2}}{x^{2}}}-1\right)\quad;\quad{\theta_{0}}\sim\left|x\right|,$
$\displaystyle\frac{{{\omega_{p}}}}{{ck}}=1-{x^{2}}\quad;\quad{\theta_{0}}\ll\left|x\right|,x<0,$

where in the second and third cases we assumed that $\theta_{0}\ll 1$.

##  Appendix D Kinetic growth rate in the general case

Proceeding to a general distribution function, we first need to simplify the
different terms in the growth rate formula. We begin by expressing the
denominator in the growth rate formula (23) via $x$ and $\theta_{0}$:

$\displaystyle\left({\cos{\theta_{2}}-\cos\theta}\right)\left({\cos\theta-\cos{\theta_{1}}}\right)=$
(89)
$\displaystyle{\sin^{4}}{\theta_{0}}\left[\begin{array}[]{l}1-\frac{{2x}}{{\tan{\theta_{0}}}}+\frac{{{x^{2}}}}{{{{\tan}^{2}}{\theta_{0}}}}-{x^{2}}{\left({1-\frac{x}{{2\tan{\theta_{0}}}}}\right)^{2}}-\\\
{\left\\{{1-\frac{x}{{\tan{\theta_{0}}}}+\frac{{{x^{2}}}}{{2{{\tan}^{2}}{\theta_{0}}}}-\frac{{{\theta^{2}}}}{{2{{\sin}^{2}}{\theta_{0}}}}}\right\\}^{2}}\end{array}\right].$

We next consider the formula in three different regimes of $x/\theta_{0}$.

  1. 1.

$\left|x\right|\ll\theta_{0}$

$\frac{\omega_{p}}{ck}=\cos\theta_{0}+x\sin\theta_{0},$ (90)
$\displaystyle\left({\cos\theta-\cos{\theta_{1}}}\right)\left({\cos{\theta_{2}}-\cos\theta}\right)=$
(91)
$\displaystyle{\sin^{2}}{\theta_{0}}-{\cos^{2}}\theta-{\left({\frac{{{\omega_{p}}}}{{ck}}}\right)^{2}}+$
$\displaystyle\frac{{2{\omega_{p}}}}{{ck}}\cos{\theta_{0}}\cos\theta.$

After some arithmetics, one gets:

$\left({\cos\theta-\cos{\theta_{1}}}\right)\left({\cos{\theta_{2}}-\cos\theta}\right)={\sin^{2}}{\theta_{0}}\left({{\theta^{2}}-{x^{2}}}\right).$
(92)

The nominator can be rewritten with a change of variables:

$\begin{array}[]{l}\left[{\frac{1}{\theta}\frac{{\partial
g}}{{\partial\theta}}\left({\cos\theta-\frac{{ck}}{{{\omega_{p}}}}\cos{\theta_{0}}}\right)-2g}\right]\theta
d\theta=\\\ \left[{\frac{{\partial
g}}{{\partial\left({{\theta^{2}}/{x^{2}}}\right)}}\left({x\tan{\theta_{0}}-{x^{2}}\left({{{\tan}^{2}}{\theta_{0}}+\frac{{{\theta^{2}}}}{{2{x^{2}}}}}\right)}\right)-g{x^{2}}}\right]d\left({\frac{{{\theta^{2}}}}{{{x^{2}}}}}\right).\end{array}$
(93)

Combining all these together, the growth rate is written as an integral over a
single variable $z=\theta^{2}/x^{2}$:

$\displaystyle\Gamma=\frac{{{\omega_{p}}}}{2}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right){\cos^{2}}{\theta_{0}}\left({1+3x\tan{\theta_{0}}}\right)\times$
(94) $\displaystyle\int_{1}^{\infty}{\frac{{\frac{{\partial g}}{{\partial
z}}\left[{{\mathop{\rm
sgn}}(x)-\left({\frac{{\left|x\right|}}{{\tan{\theta_{0}}}}}\right)\left({{{\tan}^{2}}{\theta_{0}}+\frac{z}{2}}\right)}\right]-g\left({\frac{{\left|x\right|}}{{\tan{\theta_{0}}}}}\right)}}{{\sqrt{z-1}}}dz}$
$\displaystyle\approx\frac{{{\omega_{p}}}}{2}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right){\rm{sgn}}(x){\cos^{2}}{\theta_{0}}\int_{1}^{\infty}{\frac{{\frac{{\partial
g}}{{\partial z}}}}{{\sqrt{z-1}}}dz}.$

  2. 2.

$\left|x\right|\sim\theta_{0}$
First we evaluate the denominator terms:

$\begin{gathered}\cos{\theta_{1,2}}-\cos\theta=\frac{{{\omega_{p}}}}{{ck}}\cos{\theta_{0}}\pm\sin{\theta_{0}}\sqrt{1-{{\left({\frac{{{\omega_{p}}}}{{ck}}}\right)}^{2}}}-\cos\theta=\hfill\\\
\frac{{{\theta^{2}}}}{2}-\theta_{0}^{2}+\frac{{{x^{2}}}}{2}\left[{{\text{sgn}}(x)\sqrt{1+\frac{{4\theta_{0}^{2}}}{{{x^{2}}}}}-1}\right]\pm\hfill\\\
{\theta_{0}}\left|x\right|\sqrt{\frac{{\theta_{0}^{2}}}{{{x^{2}}}}+1-{\text{sgn}}(x)\sqrt{1+\frac{{4\theta_{0}^{2}}}{{{x^{2}}}}}}=\hfill\\\
=\frac{{{x^{2}}}}{2}\left\\{\begin{gathered}\frac{{{\theta^{2}}}}{{{x^{2}}}}-\frac{{2\theta_{0}^{2}}}{{{x^{2}}}}+\left[{{\text{sgn}}(x)\sqrt{1+\frac{{4\theta_{0}^{2}}}{{{x^{2}}}}}-1}\right]\pm\hfill\\\
\frac{{2{\theta_{0}}}}{{\left|x\right|}}\sqrt{\frac{{\theta_{0}^{2}}}{{{x^{2}}}}+1-{\text{sgn}}(x)\sqrt{1+\frac{{4\theta_{0}^{2}}}{{{x^{2}}}}}}\hfill\\\
\end{gathered}\right\\}.\hfill\\\ \end{gathered}.$ (95)

Then the denominator term reduces to

$\begin{gathered}\left({\cos{\theta_{2}}-\cos\theta}\right)\left({\cos\theta-\cos{\theta_{1}}}\right)=\hfill\\\
\frac{{{x^{4}}}}{4}\left[\begin{gathered}4\mu\left[{\mu+1-{\text{sgn}}(x)\sqrt{1+4\mu}}\right]-\hfill\\\
{\left({z-2\mu+\left[{{\text{sgn}}(x)\sqrt{1+4\mu}-1}\right]}\right)^{2}}\hfill\\\
\end{gathered}\right]\hfill\\\
=\frac{{{x^{4}}}}{2}\left[{\left({2\mu-\rho}\right)\left({z-1}\right)-\frac{{{z^{2}}}}{2}}\right],\hfill\\\
\end{gathered},$ (96)

where we have defined $z=\theta^{2}/x^{2}$, $\mu=\theta_{0}^{2}/x^{2}$ and
$\rho={\mathop{\rm sgn}}(x)\sqrt{1+4\mu}-1$.

The roots of the denominator, which are also the integral boundaries, are
given by solving the proper quadratic equation for $z$:

${}z_{1,2}=2\mu-\rho\pm 2\sqrt{\mu\left(\mu-\rho\right)}.$ (97)

The nominator is evaluated using the definitions of $z,\mu$ and $\rho$ above.
Given the expression

$\frac{{{\omega_{p}}}}{{ck}}=1-\frac{{{x^{2}}}}{2}\left({\mu-\rho}\right),$
(98)

the nominator can be written as

$\left[{\frac{1}{\theta}\frac{{\partial
g}}{{\partial\theta}}\left({\cos\theta-\frac{{ck}}{{{\omega_{p}}}}\cos{\theta_{0}}}\right)-2g}\right]\theta
d\theta=\left[{\frac{1}{2}\frac{{\partial g}}{{\partial
z}}\left({\rho-z}\right)-g}\right]{x^{2}}dz.$ (99)

Combining all these together, the expression for the growth rate is

${}\Gamma={\omega_{p}}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right)\int_{{z_{1}}}^{{z_{2}}}{\frac{{\frac{1}{2}\frac{{\partial
g}}{{\partial
z}}\left({\rho-z}\right)-g}}{{\sqrt{\left({z-{z_{1}}}\right)\left({{z_{2}}-z}\right)}}}dz}.$
(100)

  3. 3.

$\left|x\right|\gg\theta_{0}$ , $x<0$

$\displaystyle\frac{{{\omega_{p}}}}{{ck}}\approx
1-{x^{2}}-\frac{{3\theta_{0}^{2}}}{2},$ (101)
$\displaystyle\sqrt{\left({\cos{\theta_{2}}-\cos\theta}\right)\left({\cos\theta-\cos{\theta_{1}}}\right)}=\frac{{{x^{2}}}}{2}\sqrt{\left({z-2}\right)\left({2+8\mu-z}\right)},$
$\displaystyle\left[{\frac{1}{\theta}\frac{{\partial
g}}{{\partial\theta}}\left({\cos\theta-\frac{{ck}}{{{\omega_{p}}}}\cos{\theta_{0}}}\right)-2g}\right]\theta
d\theta=-\left[{\frac{{\partial g}}{{\partial
z}}\left({\frac{z}{2}+1}\right)+g}\right]{x^{2}}dz.$

All together, we get

$\Gamma=\pi{\omega_{p}}\left({\frac{{mc}}{{{n_{b}}}}}\right)\left({\frac{{{n_{b}}}}{{{n_{p}}}}}\right){\left[{\left({-2\frac{{\partial
g}}{{\partial z}}}\right)-g}\right]_{z=2}}.$ (102)

The evaluation of the growth rates in section 3.2 and in the numerical
simulation is done by using these three cases, given that the formulas match
each other at intermediate ranges of $x/\theta_{0}$.

##  Appendix E Diffusion coefficients

The evaluation of diffusion coefficients (55) as a function of the beam angle
$\theta$ requires an integration over the entire range of waves resonant with
a pair-beam particle with polar angle $\theta$. After the $\phi_{0}$
integration, using the delta function, the diffusion coefficient takes the
form

${D_{\mu}}=\pi{e^{2}}\int{\frac{{W\left({{\mathbf{k}},t}\right)}}{{{\theta^{\mu}}}}\frac{{{{\left({\frac{{{\omega_{p}}}}{{kc}}-\cos{\theta_{0}}-\frac{{{\theta^{2}}\cos{\theta_{0}}}}{2}}\right)}^{\mu}}{k^{2}}dkd\cos{\theta_{0}}}}{{kc\sqrt{\left[{\cos{\theta_{0}}-{{\left({\cos{\theta_{0}}}\right)}_{1}}}\right]\left[{{{\left({\cos{\theta_{0}}}\right)}_{2}}-\cos{\theta_{0}}}\right]}}}}.$
(103)

This expression is too cumbersome to evaluate, even numerically. One has to
simplify the expressions in the integral by dividing the integration space to
two non-overlapping sub-spaces, based on the magnitude of $x/\theta_{0}$ and
$\theta_{0}/\theta$. The complete diffusion coefficient integral is then
presented as a sum of two expressions,

${}D_{\mu}=D_{1,\mu}+D_{2,\mu},$ (104)

where subscript ’1’ corresponds to large values of $\theta_{0}/x$ or
$\theta/\theta_{0}$ and ’2’ to the $(x,\theta_{0})$ space complementary. The
first subspace, ’1’, contains all points $\left(x,\theta_{0}\right)$ for which
$\theta_{0}$ is much larger than $\theta$ or than $\left|x\right|$. The second
subspace contains all points $\left(x,\theta_{0}\right)$ for which
$\theta_{0}$ is not much larger than $\left|x\right|$ and $\theta$. Subspace
’1’ is separated further into two smaller non-overlapping sub-spaces ’1a’ and
’1b’ as follows: Subspace ’1a’ ,$V_{1a}$, contains all points $(x,\theta_{0})$
such that $50\left|x\right|<\theta_{0}$ and $\theta_{0}\leq 50\theta$.
Subspace ’1b’ ,$V_{1b}$, contains all points $(x,\theta_{0})$ such that
$\theta_{0}>50\theta$. We have found that the arbitrary choice of the large
factor of 50 is suitable for the approximations made in the simplifications of
the algebraic expressions. Subspace ’2’ contains all the rest. The sub-
division of subspace ’1’ is mainly for convenience, while the integration
formula is the same. An example for the integration space
$\left(x,\theta_{0}\right)$ for the diffusion coefficients is described in
figure 13. The solid curves correspond to the boundaries of the resonance
parameter, $x$, as determined by formula (85). The complete integration space
is then the area between these two curves (ranging over all values of
$\theta_{0}$). The dashed lines correspond to the separation between the two
sub-spaces, where the numbers ’1’ and ’2’ denote the subspace.

We now proceed to algebraic development of the diffusion integration formulas
for each subspace separately. In the subspace 1,
$\theta,\left|x\right|\ll\theta_{0}$. Expanding in small parameters to the
first non-vanishing order, one gets, after some work,

${}D_{1,\mu}=\frac{{\pi{e^{2}}\omega_{p}^{2}}}{{2{c^{3}}}}\int_{{V_{1}}\left(\theta\right)}{{{\left({\frac{{{{\tan}^{2}}{\theta_{0}}}}{{1+{{\tan}^{2}}{\theta_{0}}}}}\right)}^{\mu/2}}\frac{{{z^{\mu}}}}{{\sqrt{1-{z^{2}}}}}W\left({{\mathbf{k}},t}\right)dzd{{\tan}^{2}}{\theta_{0}}},$
(105)

where $z=x/\theta$.

In the subspace 2, we define new variables

${}y=\frac{{x\sqrt{{x^{2}}+4\theta_{0}^{2}}-{x^{2}}}}{{{\theta^{2}}}}+1;\qquad
w=\frac{\theta_{0}}{\theta}.$ (106)

Then the diffusion coefficient can be written, to zeroth order in $\theta$, as

${}D_{2,\mu}^{0}=\frac{{2\pi{e^{2}}\omega_{p}^{2}}}{{{c^{3}}}}{\left({\frac{\theta}{2}}\right)^{\mu+2}}\int_{{V_{2}}\left(\theta\right)}{W\left({{\mathbf{k}},t}\right)\frac{{{{\left({y-2}\right)}^{\mu}}}}{{\sqrt{4{w^{2}}-{y^{2}}}}}d{w^{2}}dy}.$
(107)

Figure 13: Integration space for the diffusion coefficient at
$\theta=\Delta\theta=10^{-6}$. The solid lines are the curves for the upper
and lower bounds of $x$ as a function of $\theta_{0}$, respectively. The
dashed lines represent the boundaries of the two integration subspaces,
denoted by $1$ and $2$ in the figure, respectively. The subdivision of
subspace ’1’ is denoted by ’1a’ and ’1b’ and their boundaries are also denoted
by dashed lines. The vertical dashed line represents the line
$\theta_{0}=50\theta$.