RELAXATION OF BLAZAR-INDUCED PAIR BEAMS IN COSMIC VOIDS

Pairs are most efficiently created just above the energy threshold for production, i.e., where

$$\begin{equation} s\equiv E_\gamma E_{\rm EBL} (1-\cos \phi)/2m_e^2c^4\ge 1, \end{equation} \tag{ 1 }$$

with ϕ the angle between the interacting photons, E_γ and E_EBL the energy of the incident and target EBL photon, respectively, and the relativistic invariant, s, the center of mass energy square in units $m_e^2c^4$. The mean free path for the process depends on the details of the EBL model (Kneiske et al. 2004; Franceschini et al. 2008) but is approximately

$$\begin{equation} \ell _{\gamma \gamma }\simeq 0.8\left(\frac{E_\gamma }{\rm TeV}\right)^{-1} (1+z)^{-\zeta }\,{\rm Gpc}, \end{equation} \tag{ 2 }$$

with ζ = 4.5 for z ⩽ 1, ζ = 0 otherwise (Neronov & Semikoz 2009; Broderick et al. 2012). The particle number density of the beam, n_b, is set by the balance of pair production rate, 2F_γ/ℓ_γγ, evaluated close to production threshold, and energy loss rate. If inverse Compton losses dominate then, at a distance D from the blazar such that F_γ = L_γ/4πD²,

$$\begin{eqnarray} \nonumber n_b&\simeq &2 \frac{F_\gamma }{c} \frac{\ell _{{\rm IC}}}{\ell _{\gamma \gamma }} \simeq 3\times 10^{-25}\,{\rm cm}^{-3} \left(\frac{E_\gamma L_\gamma }{10^{45}\,\rm erg\, s^{-1}}\right) \\ & & \times \left(\frac{D}{\rm Gpc}\right)^{-2} \left(\frac{E_\gamma }{\rm TeV}\right) (1+z)^{\zeta -4}, \end{eqnarray} \tag{ 3 }$$

where E_γL_γ is an estimate of the blazar's equivalent isotropic gamma-ray luminosity for a source at distance D (see Broderick et al. 2012). Each pair particle carries about half the energy of the incident gamma ray, so the beam Lorentz factor is

$$\begin{equation} \Gamma =\frac{E_\gamma }{2m_ec^2}\sim 10^6 \left(\frac{E_\gamma }{\rm TeV}\right). \end{equation} \tag{ 4 }$$

Another important characteristic quantity of the beam is its angular spread, Δθ, determined by the distribution of angles θ between the pair produced particles and the parent photon's direction. This can be found to be related to Γ and the relativistic invariant as

$$\begin{equation} \Delta \theta \le \frac{s}{\Gamma }\left(1-\frac{1}{s}\right)^\frac{1}{2}. \end{equation} \tag{ 5 }$$

In this section, we compute the characteristic quantities of a blazar-induced pair beam, using a Monte Carlo model of the electromagnetic cascade, fully described in Elyiv et al. (2009) and also applied in Neronov et al. (2010). For this purpose, the blazar's spectral emission is a typical power-law distribution of primary gamma-ray photons, $dn_\gamma /dE_\gamma \propto E_\gamma ^{-q}$, in the range 10³ ⩽ E_γ/m_ec² ⩽ 10⁸ and with q = 1.8 (Abdo et al. 2010). In addition, we use the nominal model of Aharonian (2001) for the EBL in the range 0.1–1000 μm. For the cosmic background radiations beyond 0.1 μm we used data from Hauser & Dwek (2001). Contrary to Elyiv et al. (2009), we did not consider the deflection of e⁺e⁻ pairs in the extragalactic magnetic field as well as the inverse Compton interactions with cosmic microwave background (CMB) photons. Here we took into account pair distributions resulting from just the first double photon collisions. The energy distribution and cross section of the relevant reactions were taken from Aharonian (2003).

Given the primary gamma-ray photon spectrum, we generated random gamma–photon interaction events. These are characterized by the random distance from the blazar at which the double photon collision occurs based on the EBL dependent mean free path of the high-energy photons. Next, we randomly generated the energies E_± of the produced e⁺e⁻ pairs and evaluated the proper angle θ between each pair and the direction of the incident gamma-ray photon. For these quantities we used the analytical expressions for μ_e = cos (θ) in Schlickeiser et al. (2012).

The results of the calculation are shown in Figure 1 for a blazar with equivalent isotropic gamma-ray luminosity of 10⁴⁵ erg s⁻¹. The top panel shows the spectral energy distribution of the production rate of the pairs at several distances from the blazar, ranging from 3.5 Mpc (top, black line) to 1 Gpc (bottom, yellow line). The bottom panel shows the peak (dashed line) and mean (solid line) energy of the generated pairs, in units of m_ec², as a function of distance from the blazar. The Lorentz factor of the pairs is a monotonically decreasing function of distance in the range Γ = 10⁵–10⁶ up to a Gpc away from the blazar. This shows that pair production typically peaks at near-infrared (E_EBL ≃ 0.1 eV) EBL target photons interacting with gamma rays with energy E_γ ∼ (0.1–1) TeV. The vertical bars indicate the energy range encompassing 68% of the particles. In fact, there is considerable energy spread about the mean value, which decreases toward larger distances as the energy range of gamma-ray photons interacting with the EBL is also reduced. However, the beam particles always remain ultrarelativistic, a detail relevant in the analysis below.

Zoom In Zoom Out Reset image size

Figure 2 shows the angular distribution of beam pairs at distances 244 Mpc (black dash) from the blazar. The beam angular spread is of order Δθ ≃ 10⁻⁵, consistent with Equation (5) and the value of Γ estimated above. The colored (dashed and dotted) curves show the angular distribution of pairs in eight different energy bins, indicating that the beam angular spread is primarily determined by pairs with Γ ∼ 10⁴–10⁶.

Zoom In Zoom Out Reset image size

The distribution that we use for the analysis of the pair beam stability below is the steady state one obtained by balancing the production rate given in Figure 1 with inverse Compton losses, i.e.,

$$\begin{equation} f(\Gamma)=\int _0^\Gamma \frac{\tau _{{\rm IC}}(\varepsilon)}{\varepsilon }\left(\frac{dn(\varepsilon)}{d \varepsilon } \right)_{\gamma \gamma } d\varepsilon, \end{equation} \tag{ 6 }$$

where τ_IC(Γ) = ℓ_IC/c is the energy dependent timescale for inverse Compton losses.

In Table 1, we report as a function of distance, D, from the blazar, the main properties of the beam, which are relevant to the analysis below. These include the number density of the beam pairs, n_b; the mean value of the inverse of the pairs' Lorentz factor, 〈Γ⁻¹〉⁻¹, which enters the estimate of the maximum growth rate of the instability; the mean value of the pairs' Lorentz factor, 〈Γ〉, which determines the mean energy of the beam; the mean square value of the pairs' Lorentz factor divided by the mean value of the same, 〈Γ²〉/〈Γ〉, which enters the estimate of the inverse Compton timescale; and the rms of the opening angle of the pairs, Δθ, which also enters the growth rate of the instability. These Lorentz gamma factors differ considerably, and they all monotonically decrease with distance as in the bottom panel of Figure 1.

Although the results in this and the next sections assume a blazar equivalent isotropic gamma-ray luminosity of 10⁴⁵ erg s⁻¹, they can be generalized to other luminosities by rescaling the pair beam number density according to n_b → n_b(E_γL_γ/10⁴⁵ erg s⁻¹).

The analysis of the pair beam instability also depends on the thermodynamic properties of the plasma in cosmic voids, namely, the number density of free electrons and their temperature. The number density of free electrons can be expressed as n_v ≃ 2 × 10⁻⁷(1 + δ)(1 + z)³ cm⁻³. The typical overdensity δ is taken to be the value at which the cumulative distribution of the IGM gas is 0.5. Using the simulation results presented in Section 4.2, and in particular the insets in Figure 9, we estimate as representative value for the voids δ_v = −0.9(1 + z), where the redshift dependence is approximate but sufficient for our purposes. Note that this implies a redshift evolution of the bulk IGM density in voids, n_v∝(1 + z)⁴. As for the gas temperature, we assume T_v ≃ a few × 10³ K (1 + z)^1.5, which reproduces the IGM temperature at mean density of a few × 10⁴ at redshift 3. This redshift dependence, while again a rough approximation, is acceptable for our purposes.

For particles of an ultrarelativistic beam with modest angular spread, Δθ ≪ 1, we can assume v_∥ ≃ c and v_⊥ ≃ cΔθ. If the energy spread ΔE/E ≲ 1, then the longitudinal velocity spread of the beam is

$$\begin{equation} \Delta v_\parallel \simeq c\frac{\Delta E}{\langle \Gamma \rangle ^2 E}+c\Delta \theta ^2. \end{equation} \tag{ 16 }$$

Thus, for the angular spread and Lorentz factor characteristic of the blazar-induced beam, the longitudinal velocity spread is negligible with respect to the perpendicular velocity spread. It turns out that the first and second terms in Equation (16) are comparable to within a factor of a few, so for the sake of simplicity in the following we retain the second term only, neglecting any fudge factor. If we then use Equations (9) and (16), with the estimates for the beam angular spread and bulk Lorentz factor from the previous section, we find that virtually all modes require a kinetic description unless

$$\begin{equation} \frac{k_\perp }{k}\lesssim \times 10^{-5} \left(\frac{n_b/n_{\rm v}}{10^{-15}}\right) \left(\frac{\langle \Gamma \rangle }{10^5}\right)^{-1} \left(\frac{\Delta \theta }{10^{-5}}\right)^{-3}. \end{equation} \tag{ 17 }$$

The maximum growth rate occurs at k_⊥ provided by the above estimate. For smaller values, we enter the reactive regimes and relativistic inertia increases. For larger values, we are in the kinetic regime where the growth rate decreases due to the increasing velocity spread along k although the decrease becomes significant only for k_⊥ ⩾ ω_p/c. The fastest growth rate for modes with k_⊥ ⩽ ω_p/c is therefore given by

$$\begin{eqnarray} \gamma _{{\rm max}}&\simeq & \omega _p\langle \Gamma ^{-1}\rangle \frac{n_b}{n_{\rm v}} \frac{1}{\Delta \theta ^2} = 4\times 10^{-12}\,{\rm s^{-1}} \left(\!\frac{n_{\rm v}}{\rm 2\times 10^{-8}\,cm^{-3}} \!\right)^{-\frac{1}{2}} \nonumber\\ && \times \left(\frac{\langle \Gamma ^{-1}\rangle }{10^{-4}}\right)^{-1} \left(\frac{\Delta \theta }{10^{-4}}\right)^{-2} \left(\frac{D}{\rm Gpc}\right)^{-2}, \end{eqnarray} \tag{ 18 }$$

where we have taken n_b ≃ 10⁻²⁴ cm⁻³ at a Gpc from a blazar of luminosity E_γL_γ = 10⁴⁵ erg s⁻¹. A basic condition for the growth of an instability is that its growth rate exceeds the collisional damping rate, i.e., γ_max ≫ ν_c, where (Huba 2009)

$$\begin{equation} \nu _c\simeq 10^{-11}\,{\rm s^{-1}}\left(\frac{n_{\rm v}}{\rm 2\times 10^{-8}\, cm^{-3}}\right)\left(\frac{T_{\rm v}}{\rm 3\times 10^3 \,K}\right)^{-\frac{3}{2}}, \end{equation} \tag{ 19 }$$

and for the Coulomb logarithm we have used Λ_c = 27.4 (Huba 2009). In Figure 6, we plot the ratio γ_max/ν_c as a function of distance from the blazar, using the values reported in Table 1, which again applies for a blazar of equivalent isotropic gamma-ray luminosity of 10⁴⁵ erg s⁻¹. The solid, dashed, and long-dashed curves correspond to redshifts z = 0, z = 1, and z = 3, respectively. The redshift dependence is obtained by using the void average density and temperature redshift dependences discussed in Section 2.3, together with the redshift dependence of n_b given in Section 2.1. The shaded area corresponds to the region where the instability is inhibited by collisions. The plot shows that the instability can only develop at distances of less than a 50 Mpc at redshift z = 0 and about 20 physical Mpc z = 3.

Zoom In Zoom Out Reset image size

In this section, we consider in some detail the main process that we believe compensates the growth of Langmuir waves, i.e., induced scattering off plasma ions, also known as nonlinear Landau damping (Tsytovich & Shapiro 1965; Breizman et al. 1972; Lesch & Schlickeiser 1987). In this process, a thermal ion, with characteristic velocity v_ti, interacts with the beam wave produced by two Langmuir oscillations, $\omega (\boldsymbol {k}), \omega (\boldsymbol {k}^\prime)$, under the condition for Cherenkov interaction, i.e.,

$$\begin{equation} \omega ({\bf k})-\omega ({\bf k}^\prime)=({\bf k}-{\bf k}^\prime)\cdot {\bf v}_{ti}. \end{equation} \tag{ 20 }$$

The rate of induced scattering of Langmuir waves off thermal ions in a Maxwellian plasma with number density n and ion/electron temperature T_i/T_e, respectively, is (e.g., Melrose 1989)

$$\begin{eqnarray} \gamma _{\rm nl}({\bf k}) &=& \frac{3(2\pi)^\frac{1}{2}}{2}\frac{T_iT_e}{(T_i+T_e)^2} \int \frac{d^3{\bf k}^\prime }{(2\pi)^3} \frac{\hat{W}(k^\prime)}{nm_ev_{ti}} \nonumber\\ && \times \left(\!\frac{{\bf k}\cdot {\bf k}^\prime }{kk^\prime } \!\right)^2 \frac{k^{\prime 2}-k^2}{|{\bf k}^\prime -{\bf k}|} \exp {\left[\!-\frac{1}{2}\left(\!\frac{3}{2}\frac{v_{te}^2}{\omega _pv_{ti}} \frac{k^{\prime 2}-k^2}{|{\bf k}^\prime -{\bf k}|} \!\right)^2 \!\right]}, \nonumber\\ \end{eqnarray} \tag{ 21 }$$

where $\hat{W}(k)$ indicates the spectral energy density of Langmuir waves. The growth rate, γ_nl(k), bears the sign of (k' − k). This indicates that as a result of induced scattering, Langmuir waves cascade toward regions of phase space of lower wavevectors, i.e., lower energies, the energy difference being absorbed by the thermal ions. Eventually, the wave energy is transferred to modes with wavenumber, k, small enough that the wave phase-speed, ω/k > c, exceeds the speed of light, and resonance with the beam particles is lost. The wavenumbers allowed in the scattering process are constrained by the integral expression in Equation (21). In particular, the following condition must be fulfilled:

$$\begin{equation} \frac{|k^{\prime 2}-k^2|}{|{\bf k}^\prime -{\bf k}|} \le \omega _p \frac{v_{ti}}{v_{te}^2} \simeq 35\times \frac{\omega _p}{c}\left(\frac{T_{\rm v}}{3\times 10^3\,{\rm K}}\right)^{-\frac{1}{2}}. \end{equation} \tag{ 22 }$$

The above constrain is satisfied for the case of differential scattering, i.e., Δk/k ≪ 1, whereby k ∼ k' and k' ∼ −k. In this case, Langmuir waves, generated with wavenumber, k ∼ ω_p/c, parallel to the beam, are isotropized. This reduces the level of resonant energy density by a factor ∼Δθ²/4π, somewhat increasing the lifetime of the beam. However, given the low temperature of the IGM in voids, the more efficient integral scattering, with k' ≪ k, is also allowed. In this case, Langmuir waves are mostly kicked out of resonance in a single scattering event, dramatically reducing the level of resonant energy density and suppressing the instability.

The general solution for the evolution of the energy density in plasma waves is a non-trivial task as it requires solving for integro-differential equations that describe the detailed energy transfer of the wave energy across different modes. However, for a conservative estimate, we can neglect differential scattering and evaluate the rate of induced integral scattering from Equation (21) using the condition k' ≪ k. We thus obtain

$$\begin{equation} \gamma _{\rm nl}\simeq \omega _p \frac{W_{{\rm nr}}}{n_{\rm v}k_BT_{\rm v}}\frac{v_{te}^2}{v_{ti}\,c}, \end{equation} \tag{ 23 }$$

where W_nr is the total energy density in Langmuir waves at k ≪ ω_p/c, i.e., non-resonant with the beam. This energy density is excited by the nonlinear scattering process and for the most part is dissipated by Coulomb collisions at a rate ν_c (see further discussion below). Thus it evolves according to

$$\begin{equation} \frac{\partial W_{{\rm nr}}}{\partial t}=2\tilde{\gamma }_{\rm nl}W_{{\rm nr}} W_{r}-\nu _c W_{{\rm nr}}, \end{equation} \tag{ 24 }$$

where we have used $\tilde{\gamma }_{\rm nl}\equiv \gamma _{\rm nl}/W_{{\rm nr}}= \omega _p(1/n_{\rm v}k_BT_{\rm v})(v_{te}^2/v_{ti}c)$ and W_r is the total energy density in Langmuir waves at k ∼ ω_p/c, i.e., resonant with the beam. The latter obviously evolves according to

$$\begin{equation} \frac{\partial W_r}{\partial t}=2\gamma _{{\rm max}} W_r-2\tilde{\gamma }_{\rm nl}W_{{\rm nr}}W_r, \end{equation} \tag{ 25 }$$

where we have neglected the role of collisions (i.e., we assume γ_max ≫ ν_c, as required for the existence of the instability). Equations (24) and (25) form a well-known Lotka–Volterra system of coupled nonlinear differential equations, which has stable periodic solutions, with the following average values for the energy densities:

$$\begin{equation} \overline{W}_{{\rm nr}}=\frac{\gamma _{{\rm max}}}{\tilde{\gamma }_{{\rm nl}}}, \quad \overline{W}_{r}=\frac{\nu _{c}}{2\tilde{\gamma }_{{\rm nl}}}. \end{equation} \tag{ 26 }$$

In this regime, the transfer rate of Langmuir waves out of resonance by nonlinear Landau damping equals on average their production rate, i.e., γ_nl ≃ γ_max. Thus, the beam emission of Langmuir waves is only linear in time, with an average power $P(W_r)=2\overline{\gamma }_{{\rm max}}W_{r}$, and the beam relaxation timescale at redshift z = 0 is

$$\begin{eqnarray} \tau _{{\rm beam}} &\simeq & \frac{n_b\langle \Gamma \rangle m_e c^2}{2\gamma _{{\rm max}}\overline{W}_{r}}= 1.5\times 10^{9}\,{\rm yr} \left(\frac{n_{\rm v}}{\rm 2\times 10^{-8}\, cm^{-3}}\right)^{-1} \nonumber\\ && \times \left(\frac{\langle \Gamma ^{-1}\rangle }{10^{-4}}\right) \left(\frac{\langle \Gamma \rangle }{10^{5}}\right) \left(\frac{\Delta \theta }{10^{-4}}\right)^2 \left(\frac{T_{\rm v}}{3\times 10^{3}\,\rm K}\right). \quad \end{eqnarray} \tag{ 27 }$$

The above timescale should be compared with the pairs' cooling time on the CMB, τ_IC = ℓ_IC/c ≃ 3 × 10⁶(E_±/TeV)⁻¹(1  +  z)⁻⁴ yr. The ratio of these timescales is plotted in Figure 7 using the values reported in Table 1 as a function of distance from our reference blazar with isotropic gamma-ray luminosity of 10⁴⁵ erg s⁻¹. At redshift z = 0 (solid line), the beam appears to be stable on significantly longer timescales than the inverse Compton emission energy loss timescale, particularly within 100 Mpc from the blazar, where the average value of Γ of the pairs tends to be higher. This conclusion is reinforced at higher redshifts (dashed line for z = 3), where the redshift dependence is inferred as in Section 3.1.

Zoom In Zoom Out Reset image size

The above analysis works in the weak turbulence regimes, which require that the energy density of resonant Langmuir waves be a small fraction of the beam energy density. For the typical values of IGM gas and beam parameters used above, this requirement is readily fulfilled as $\overline{W}_{{\rm nr}}/n_b\Gamma m_ec^2\simeq 3\times 10^{-6}$, warranting our approach limited to second-order processes.

Another consistency check to be performed concerns the assumption of collisional dissipation of the long-wavelength Langmuir waves. In fact, accumulation of energy in these non-resonant waves can generate modulation instability of the background plasma if (Breizman 1990)

$$\begin{equation} \frac{\overline{W}_{{\rm nr}}}{n_{\rm v}k_BT_{\rm v}} \ge k^2\lambda _{D}^2\sim k_BT_{\rm v}/mc^2. \end{equation} \tag{ 28 }$$

In Figure 8, the ratio on the left-hand side of the above equation is plotted as a function of distance from our reference blazar using the values in Table 1. The solid line corresponds to redshift z = 0 and the horizontal line is the nominal threshold value for the onset of the modulation instability. For redshifts higher than z = 0, we plot for comparison with the same threshold line the same ratio on the left-hand side of Equation (28) but divided by a factor (1 + z)³ to account for the IGM temperature redshift dependence (see Section 2.3).

Zoom In Zoom Out Reset image size

The plot shows that the assumption of collisional dissipation of the long-wavelength Langmuir waves is always valid except at short distances from low-redshift blazars. The modulation instability deserves more attention than the scope of the current paper can afford. Here we notice that, while the modulation instability could help stabilize the beam (Nishikawa & Ryutov 1976), it could also provide an effective dissipation rate that is more efficient than collisions. In this case, the level of energy density of resonant waves, W_r, will increase with consequent reduction of the beam lifetime (see Equations (26) and (27)). However, because the threshold condition for triggering the modulation instability depends quadratically on the temperature (see Equation (28)), should the background plasma suffer even modest heating caused by the beam relaxation, the modulation instability will quickly stabilize, restoring the conditions for collisional dissipation of the non-resonant waves.

Therefore, in conclusion, from the above analysis together with the findings in Section 3.1, it appears that the beam is stable at basically all relevant distances from the blazar, and the beam instability plays only a secondary role on the electromagnetic shower, the beam dynamics, and the thermal history of the IGM.

In addition to the kinetic effects described above, the energy density of the plasma waves evolves in time due to spatial gradient effects according to (Kaplan & Tsytovich 1973)

$$\begin{equation} \frac{d}{dt}\hat{W}(x,t,k) = \frac{\partial \hat{W}}{\partial t} + {\bf v}_g\nabla _{\bf x}{\hat{W}} -\nabla _{\bf x}\omega \cdot \nabla _{\bf k}{\hat{W}}, \end{equation} \tag{ 29 }$$

where for the rate of change of the wavevector we have used the equation of geometric optics

$$\begin{equation} \frac{d\bf k}{dt}=-\nabla _{\bf x}\omega. \end{equation} \tag{ 30 }$$

The first and second terms on the right-hand side of Equation (29) describe, as usual, explicit time dependence and the effects of spatial gradients discussed at the beginning of Section 4. The last term describes the change in $\hat{W}$ associated with modifications of the waves' wavevector as a result of inhomogeneities. This term is important because, just like induced scattering by thermal ions, it transfers the excited Langmuir waves to wavemodes that are out of resonance with the beam particles, therefore suppressing the instability (Breizman & Ryutov 1971; Nishikawa & Ryutov 1976).

For Langmuir waves, the most important contribution to ∇_xω comes from density inhomogeneities. In addition, the beam stabilization mainly results from changes in the longitudinal component of the wavevector. Therefore, we restrict our analysis to this case only, and write

$$\begin{equation} \frac{d k_\parallel }{dt}\simeq \frac{1}{2}\frac{\omega _p}{\lambda _\parallel }, \end{equation} \tag{ 31 }$$

with $\lambda _\parallel =n_{\rm v}/(\vec{\nabla }n_{\rm v})_\parallel$, the length scale of the density gradient along the beam.

In order to estimate the scale lengths of IGM density gradients, λ, we have carried out a cosmological simulation of structure formation including hydrodynamics, dark matter, and self-gravity as described in Miniati & Colella (2007). For the cosmological model, we adopted a flat ΛCDM universe with the following parameters: total mass density, normalized to the critical value for closure, Ω_m = 0.2792; normalized baryonic mass density, Ω_b = 0.0462; normalized vacuum energy density, Ω_Λ = 1 − Ω_m = 0.7208; Hubble constant H₀ = 70.1 km s⁻¹ Mpc⁻¹; spectral index of primordial perturbation, n_s = 0.96; and rms linear density fluctuation within a sphere of comoving radius of 8 h⁻¹ Mpc, σ₈ = 0.817, where h ≡ H₀/100 (Komatsu et al. 2009). The computational box has a comoving size L = 50 h⁻¹ Mpc, which is discretized with 512³ comoving cells, corresponding to a nominal spatial resolution of 100 h⁻¹ comoving kpc. The collisionless dark matter component is represented with 512³ particles with mass 6 × 10⁵ h⁻¹ M_☉.

Figure 9 shows the range of scale lengths of IGM density gradients as a function of IGM gas overdensity for three different cosmological redshifts, z = 0 (top), z = 1 (middle), and z = 3 (bottom). Accordingly, the distance covered by the beam particles during the fastest growth time, ${\sim } c \gamma _{{\rm max}}^{-1}\simeq$ 1 kpc, is much shorter than the typical scale length of density gradients. This corresponds to the case of regular, as opposed to random, inhomogeneities.

Zoom In Zoom Out Reset image size

It is clear that in order for the excited waves to have an effect on the beam, the beam–wave interaction under resonant conditions must continue for a sufficiently long time. Therefore, the condition for wave excitation is expressed as (Breizman & Ryutov 1971)

$$\begin{equation} \gamma _{{\rm max}}\frac{\Delta k_\parallel }{|dk_\parallel /dt|}= \gamma _{{\rm max}}\frac{2\lambda _\parallel \Delta k_\parallel }{\omega _p}>\Lambda _c, \end{equation} \tag{ 32 }$$

where Λ_c is the Coulomb logarithm and Δk_∥ the change in longitudinal component of the wavevector allowed by the resonant condition (Equation (7)). Using Equation (7) (and neglecting the term ΔE/E Γ²) to obtain Δk_∥ ≲ (ω_p/c)Δθ² + k_⊥Δθ, Equation (32) can be solved to express the condition for wave excitation in terms of λ_∥, i.e.,

$$\begin{eqnarray} \lambda _\parallel &\ge & \frac{c}{2\omega _p}\langle \Gamma ^{-1}\rangle \frac{n_{\rm v}}{n_b}\Lambda _c \left(1+\frac{k_\perp }{k_\parallel \Delta \theta }\right)^{-1} > 10^{6}\,{\rm kpc} \nonumber\\ && \times \left(\frac{D}{\rm Gpc}\right)^2 \left(\frac{\langle \Gamma ^{-1}\rangle }{10^{-4}}\right) \left(\frac{\Delta \theta }{10^{-4}}\right) \left(\frac{\Lambda _c}{30}\right)(1+z)^2, \quad \end{eqnarray} \tag{ 33 }$$

where in the second inequality we have used κ_⊥ ≃ κ_∥, which corresponds to the most favorable case for the instability growth in the presence of inhomogeneities, and again the redshift dependence is derived as described in Section 3.1. We can again plot the minimal values of λ_∥ allowed for the growth of the beam instability using the parameter values for the beam from Table 1. This is shown by the oblique lines in Figure 10 for redshifts z = 0 (solid line), z = 1 (dashed line), and z = 3 (long dashed line). For the same redshifts (and with line-style matching that of the corresponding oblique lines), the three horizontal thin lines represent the mean scale length of density inhomogeneities at typical void overdensity (i.e., where the cumulative gas distribution function is 0.5), extracted from Figure 9. The figure shows that the growth of Langmuir waves is constrained by the presence of inhomogeneities, to distances D < 30, 6, and1 Mpc from the blazar, for z = 0, 1, and 3, respectively. Inhomogeneities provide another independent argument against the growth of Langmuir waves and the unstable behavior of the pair beam. While nonlinear Landau damping weakens with distance from the blazar (see Figure 6), the impact of inhomogeneities becomes stronger (see Figure 10), so that the stabilization effects of the two processes compensate each other at different distances.

Zoom In Zoom Out Reset image size