© ESO, 2016

1. Introduction

Active galactic nuclei (AGN) at cosmological distances may emit photons with energies well above 1 TeV, but the mean free path of these photons is limited by the interaction with the extragalactic background light (EBL). Even though there are some uncertainties with regard to the spectrum of the EBL, the mean free path for photons with energies E_γ above 1 TeV can be estimated to be D ≈ 80(E_γ/ 10 TeV)^-1 Mpc (Neronov & Semikoz 2009). The interaction of high-energy photons with the EBL produces electron-positron pairs, where each of the leptons carries approximately half of the energy of the incident gamma-ray photon.

These very-high-energy pair beams are subject to several physical processes: deflection by the intergalactic magnetic field (Neronov & Vovk 2010), initiation of plasma instabilities in the thermal background medium (Broderick et al. 2012; Schlickeiser et al. 2012a,b; Miniati & Elyiv 2013), or inverse Compton scattering off the cosmic microwave background (CMB; Kneiske et al. 2004). These processes affect observable quantities. While the third process should yield photons upscattered to GeV energies, the other processes will add to the heating of the intergalactic medium instead.

Observations by the Fermi satellite (Abdo et al. 2009) have opened up the window for observations in the GeV photon range. This enables observers to look for a possible GeV signal from secondary radiation in distant AGN. Since the spectrum of the EBL is known to some degree and AGN are observed at TeV energies, the tentative spectrum of GeV photons can be calculated. The question at hand is which process dominates the evolution of the electron-positron pair beam. If magnetic field deflection is the dominating process, the beam will dissolve and the inverse Compton signal will be reduced substantially. On the other hand, the evolution of a plasma instability might be the fastest process. In that case, the question becomes whether the change of the particle distribution function is sufficient to prevent any of the other processes.

The absence of a GeV signal in certain blazars has been attributed solely to magnetic field deflection in order to determine the strength of the intergalactic magnetic field (Neronov & Vovk 2010). However, this approach does not take into account any effects arising from plasma instabilities. Owing to the highly non-thermal distribution function of electrons in the intergalactic medium (IGM), there is a lot of potential for plasma instabilities to develop.

While the growth rates of plasma instabilities caused by ultra-relativistic electrons can be calculated at the onset of the instability, later nonlinear phases require several assumptions and approximations in an analytic treatment. The ultimate outcome is determined, in particular, by the interaction of turbulence with the initial distribution function and other instabilities. Unfortunately, simulations cannot reflect the extreme parameter space of physical reality. Therefore, a simplified system needs to be constructed carefully in order to capture all the essential physical processes of the phenomena to be examined. In previous particle-in-cell (PiC) simulations Sironi & Giannios (2014) studied a wide range of density ratios α and beam Lorentz factor γ_b to allow for a reliable extrapolation. In this article we argue that the behavior of the system is governed principally by the energy density of the beam compared to the background energy density. Consequently, quantities like the high beam gamma factor and high density ratios need not be as problematic for PiC simulations as long as the full range of ϵ is studied, covering values both larger and smaller than unity. The latter case is especially interesting, because it has not been realized for any of the parameters discussed by Sironi & Giannios (2014).

2. Physical parameters

The physical setup of the system consists of hot, low-density ionized gas in thermal equilibrium and a high-intensity photon population with energies up to 10 TeV and a second photon population in the IR-UV range (the EBL). The parameters of the IGM are, to some degree, uncertain, but there is some consensus that the temperatures are in the range of 10⁴ K in cosmic voids (Hui & Gnedin 1997) to 10⁷ K in the intra-cluster medium (Mernier et al. 2015). Following the approach by Chang et al. (2012) we adopt a value of 10⁶ K and a density of n_background = 10^-7 cm^-3, which seems to be a good approximation when AGN heating plays a role. Since we are dealing with the interaction of pair beams emanating from AGN, it is a sensible approach to also take into account the heating effects. The magnetic field is hard to determine because it is too weak to be observable through effects such as Zeeman splitting or Faraday rotation, which are used to measure the strength of stronger, solar magnetic fields. This leaves the wide range between primordial magnetic fields that have not been amplified by dynamo processes and are as weak as 10^-24 G and an upper limit of about 10^-9 G (Kronberg 1994) based on quasar observations. More recently the spectrum (Essey et al. 2011) and surrounding halo (Chen et al. 2015) of distance AGN have been used to narrow this down to the range 10^-17 G...3 × 10^-14 G. This is consistent with the value of 10^-15 G favored by (Ando & Kusenko 2010).

From the AGN photons and the EBL photon field the distribution function of electrons in the beam can be calculated. This will typically produce a power-law of the pairs resembling the power law of the gamma rays from the source (Schlickeiser et al. 2012b). Since the minimal energy needed to produce pairs is already present at TeV energies, the width of the production spectrum is typically very small.

Because of the limited sensitivity of Air Cerenkov Telescopes and of the EBL absorption, it is not possible to make definite statements about the upper end of the AGN spectrum. As soon as the EBL has attenuated the AGN spectrum enough to make it fall under the telescope’s sensitivity, only models can give a clue to that part of the spectrum.

In order to describe the system, we make use of the parameters α = n_jet/n_background and ϵ = e_jet/e_background, where e is the energy density (e_jet = n_jetγm_ec² and e_background = n_backgroundk_BT). Typical parameters for the pair beam as derived by Schlickeiser et al. (2012a) are γ = 10⁶ and n_jet = 10^-22 cm^-3. The actual parameters can differ by several orders of magnitude owing to the extremely different AGN, but also on the EBL spectrum at the redshift of the source.

3. Methods

The simulations presented in this paper are prepared using the PiC (Hockney & Eastwood 1988) code ACRONYM (Kilian et al. 2012). A simple 2D Cartesian grid topology (2D3V) with periodic boundary conditions in both spatial directions is used. Moreover, the electric and magnetic field information is stored according to the standard Yee (Yee 1966) scheme. In addition, the current density is calculated following the Esirkepov scheme (Esirkepov 2001). Current and charge density are deposited on the grid via a triangular shaped cloud (TSC) form factor, providing a good balance of computational speed and numerical accuracy (Kilian et al. 2013). In the electromagnetic case, Maxwell’s equations are solved with an explicit second-order leap-frog scheme. For the electrostatic case, the electric field is obtained as a solution to Poisson’s equation provided by a Fourier solver. Field interpolation to particle locations proceeds with the TSC form factor used above. The relativistic equation of motion is applied through the implicit method described by Vay (Vay 2008), although the standard Boris push (Boris 1970) proved to be sufficient during testing.

In order to suppress unphysical field fluctuations arising from the grid Cherenkov instability, a Friedmann filter with a filtering parameter Θ = 0.3 is introduced in the electromagnetic case (Greenwood et al. 2004). Furthermore, the current density is filtered with a spatial binomial filter including a compensation pass.

4. Numerical setup

All simulations are performed with the following setup. Four particle populations are distributed homogeneously throughout the simulation volume with equal macro-particle numbers per cell. The background consists of electrons and protons of natural mass ratio at thermal equilibrium at temperature T. A given number density n_e = n_p is achieved by scaling the macro-particles appropriately, as required for a PiC simulation.The pair beam consists of positrons and electrons with number density αn_e each. Its distribution is a Maxwellian of temperature T with per dimension, drifting along the x-axis with a relativistic speed given by γ_b. This particle setup effects no net charges or currents. Initially, electric and magnetic field components are set to zero. Moreover, a cell size Δx is chosen so that it is small enough to resolve the electron inertial length with several cells. Lastly, the time step size Δt is given by the CFL criterion. Parameters for a basic simulation run are given in Table 1.

Below, simulations of varying beam Lorentz factor γ, density ratio α, and energy density ratio ϵ are presented. Comparison simulations with larger box sizes and different simulation sizes are also provided.

In Chang et al. (2012) four relevant processes have been identified: growing electrostatic fluctuations and aperiodic fluctuations (which convert beam energy into kinetic energy of the background plasma) and the modulational instability and nonlinear Landau damping (which turn the plasma energy into heat).

Their respective growth rates are given by electrostatic fluctuation (1)aperiodic (Weibel) fluctuations (2)modulation instability (3)and nonlinear Landau damping (4)The parameters N₇, n₂₂, T₄, and Γ₆ describe the physical system with typical numbers, where n_bg = N₇10^-7 cm^-3, n_jet = n₂₂10^-22 cm^-3, γ = 10⁶Γ₆, and T = 10⁴ K. In our scenario α and ϵ are the relevant quantities. Here N₇ = 1 and T₄ = 100 are constants. Rewriting the above equations using our set of variables yields For the typical values of α and ϵ used in this paper the ratio of (9)is larger than unity, as it is in the original paper. In addition, the ordering of modulational instability and nonlinear Landau damping are conserved.

5. Results

5.1. Basic simulation

The simulation setup was first tested with a basic simulation: ϵ = 1, α = 2.5 × 10^-5, γ = 10. This setup can be regarded as an extreme scenario since in typical physical environments the beam strength is expected to be smaller than the thermal energy of the background.

We have identified several observables as relevant for the development of the instability and the evolution of the distribution function. In Fig. 1 the energy of the electric field is plotted over time. The magnetic field energy follows basically the same curve. A sharp increase can be seen until , which then relaxes until the electric field energy saturates at around .

The growth rates are about γ = 3 × 10^-4ω_pe at the beginning and grow to γ = 7 × 10^-4ω_pe. The analytical values predict a maximum growth rate of γ_E = 1.9 × 10^-2ω_pe for the electrostatic instability and γ_W = 3.8 × 10^-3ω_pe for the aperiodic fluctuation, while the modulation instability should kick in at .

Fig. 1 Electric field energy ϵ = 1, α = 2.5 × 10-5, γ = 10.

On the one hand, it should be noted that the growth rates are really only maximum growth rates, so that in this respect the observed growth rates are in line with the analytical predictions; on the other hand, the average growth rates in the simulations are in a similar range to the analytical calculations. A detailed look at Fig. 2 highlights a sharp peak of the Fourier transformed electric field energy at 4 × 10^-10 cm^-1 and another bump at 2.3 × 10^-10 cm^-1. The 2D Fourier transform shows (Fig. 3) that the former peak is linked to a structure at fixed k_∥ over a wide range of k_⊥. This is not predicted by Chang et al. (2012), who claim that the maximum growth occurs at a fixed angle independent of the absolute of k. The position in k-space is compatible with k_∥u ≈ k_∥c = ω_pe.

Fig. 2 E-field 1D Fourier ϵ = 1, α = 2.5 × 10-5, γ = 10 at .

Fig. 3 E-field 2D Fourier ϵ = 1, α = 2.5 × 10-5, γ = 10 at .

When observing the change of the distribution function as outlined in Fig. 4, it can be seen that the peak of the beam decreases in amplitude, but even after more than 10⁴ plasma timescales it just moves to a plateau. The prediction of Chang et al. (2012), which invokes results by Grognard (1975), says that after around 5000 plasma timescales the peak should have vanished. We conclude that the plateauing of the distribution function slows down the instability, which brings the system to an almost stable situation after around 7000 plasma timescales.

Fig. 4 Particle histogram ϵ = 1, α = 2.5 × 10-5, γ = 10.

Fig. 5 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), E-field 1D Fourier, resolution 1024 × 1024 at .

Fig. 6 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), electric field energy, resolution 1024 × 1024.

5.2. Faster beams

Since the basic simulation is limited to a beam Lorentz factor of only 10, simulations with faster beams are necessary to make statements about the evolution of beams in the actual physical setting. Our setting for a mildly faster beam was ϵ = 1, α = 1.25 × 10^-5, γ = 20. Early attempts with a similar resolution (Figs. 5, 6) led to results showing a significantly different behavior from that of the slower beam. It turned out that faster beams require a higher resolution. The results analogous to the γ = 10 case are shown in Figs. 7–9.

When calculating the theoretically predicted values, the results are only slightly different to the basic case: γ_E = 1.2 × 10^-2ω_pe for the electrostatic instability and γ_W = 1.9 × 10^-3ω_pe for the aperiodic fluctuation. For the well-resolved simulation there is a decline in the electric field energy at around 6000 plasma timescales, but the distribution function changes basically in the same way as before.

Fig. 7 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), E-field 1D Fourier, resolution 2048 × 2048 at .

Fig. 8 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), electric field energy, resolution 2048 × 2048.

Fig. 9 Faster beam, E-field 2D Fourier, resolution 2048 × 2048 at .

5.3. Strong beam

From the scenario of the faster beam we developed a strong beam scenario which also has a beam with Lorentz factor γ = 20, but a different energy density ratio. We assumed ϵ = 10 with a constant density ratio of α = 1.25 × 10^-4 and again with a resolution of 2048 × 2048.

The results of this simulation run differ considerably from the previous simulation; the growth rates are much higher than in the basic simulation (7 × 10^-4ω_pe rising to 7 × 10^-3ω_pe). This is only partially expected from the theoretical calculations. The expected growth rate for the electrostatic fluctuations are only marginally higher (2.6 × 10^-2ω_pe), while the aperiodic fluctuation rate is even lower (6 × 10^-4ω_pe). A close look at the 2D Fourier transform in Fig. 11 shows the already well-known features at k_∥ = 3 × 10^-10 cm^-1, but additionally it shows a ringlike structure at | k | = 5 × 10^-10 cm^-1. The fact that this feature only appears for ϵ> 1, which can be considered an unphysical scenario, puts the results of Sironi & Giannios (2014) into question.

Fig. 10 Strong beam (ϵ = 10, α = 1.25 × 10-4, γ = 20), electric field energy.

Fig. 11 Strong beam (ϵ = 10, α = 1.25 × 10-4, γ = 20 at ), E-field 2D Fourier.

Fig. 12 Strong beam (ϵ = 10, α = 1.25 × 10-4, γ = 20), particle histogram.

An even more interesting result can be seen when inspecting the histogram in Fig. 12. The peak shows a sharp decrease in amplitude (which would only be a qualitative difference to the basic simulation), but it also shows a shift in the peak position. This is a quantitive difference and the only case in this series of simulations where the beam suffers an actual loss of energy.

5.4. Weak beam

Fig. 13 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, E-field 1D Fourier at .

Fig. 14 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, electric field energy.

Fig. 15 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, E-field 2D Fourier at .

The last setup is the weak beam with ϵ = 0.1, α = 2.5 × 10^-6, γ = 10. This is a mostly realistic scenario with regard to the energy density ratio.

The observed growth rate (Fig. 14) of the instability here is approximately 2 × 10^-4ω_pe with theoretical values of 8.8 × 10^-3ω_pe (electrostatic) and 1.2 × 10^-3ω_pe (aperiodic). A saturation stage is reached very early. The 1D Fourier transform (Fig. 13) and 2D Fourier transform (Fig. 15) are similar to the basic simulation, but the instability peak is much less pronounced.

An important feature seen here is the negligible change in the pair beam distribution function as shown in Fig. 16.

Fig. 16 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, particle histogram.

6. Discussion

In the simulations presented in this article we were able to recover the onset of an instability caused by electron-positron pair beams. The general existence of this instability is not affected by the change of the energy density ratio of beam and background. The actual development of the instability depends strongly on the energy density ratio ϵ.

Common to all simulations, regardless of the beam Lorentz factor γ and the energy density ratio ϵ is the excitation at k_∥ = 3 × 10^-10 cm^-1 over a whole range of perpendicular wave lengths k_⊥. The strength of this emission depends on the beam strength (stronger beams lead to stronger emission) and the beam speed (faster beams produce stronger emission). The two opposite cases of strong and weak beams produce interesting opposing results. Owing to the weakness of the beam, Fig. 15 is mostly dominated by noise, while the strong beam produces a clear signal along the described k_∥ axis as seen in Fig. 11. An additional feature seen there is a ring-like structure in the k-space with a radius of | k | = 5 × 10^-10 cm^-1.

The 1D energy spectra reflect this finding. Regardless of ϵ and γ they all show falling spectra with a peak around | k | = 3 × 10^-10 cm^-1. In the case of ϵ = 10 another peak related to the ring-like structure can be seen. The under-resolved case in Fig. 5 shows an additional hump at | k | = 2.2 × 10^-10 cm^-1 that vanishes when the velocity increases. This simulation has, therefore, been rejected as unphysical.

One important thing to note is that in all simulations a linear stage can be observed on timescales shorter than ; after that the rise of the instability slows down and eventually reaches a maximum. For the case of a strong and fast beam, a decline in electric field energy is also observable. The difference between these different runs can be explained by taking a look at the distribution functions.

For the basic simulations we see a broadening of the peak distribution function after a rather long timescale of more than . This spread continues, but even after the peak has only broadened by 20% and the peak amplitude is still a fifth of the initial amplitude. When this is compared to the strong beam scenario, it becomes clear that the evolution of the distribution function happens much faster and changes the distribution function much more drastically. Not only has the distribution function changed after less than , but the position of the maximum has also shifted. This may be linked to the strong decline of the instability observed in the electric field.

The weak beam scenario shows a drastically different behavior. Even after more than the distribution function of beam pairs is almost completely unchanged. We want to emphasize that we resolve the timescales for the modulation instability and nonlinear Landau damping in this simulation (cf. Schlickeiser et al. 2012b).

We conclude that for the simulation setup presented here, the change in the distribution is negligible over longer times when the energy density ratio is below 1. For actual physical scenarios we would typically expect even lower values of ϵ than those shown here. The loss channels of inverse Compton scattering and magnetic field deflection are, therefore, still open for the pair beams. This is consitent with assumptions that are usually made in the determination of the intergalactic magnetic field based on observations of distante point sources.

Given the great numerical effort involved, our simulations are limited to low γ values, but as the comparison of the basic scenario and the fast beam show, ϵ is the far more decisive parameter in the evolution of the system as such.

7. Conclusion

We have shown simulations of systems containing a hot, thermal proton-electron plasma and mildly relativistic electron-positron beams. Our simulations suggest that for low-energy density ratios (i.e., less energetic beams) instabilities are created that do not lead to a strong change in the distribution function of the beam.

When taking into account the physical scenario, we would conclude that while the instabilities in question may broaden the beam distribution, they do not provide enough energy loss to explain missing GeV photons. This brings us back to the original paper by Neronov & Vovk (2010). We do not observe upscattered GeV photons from EBL generated pair beams and the instability does not successfully remove the beam electrons; therefore, magnetic deflection may still be the governing process.

Acknowledgments

We acknowledge the use of the ACRONYM code and would like to thank the developers (Verein zur Förderung der kinetischen Plasmasimulationen e.V.) for their support. The work of A.K. was supported by the Deutsche Forschungsgemeinschaft through grant Schl 201/31-1. F.S. acknowledges support from NRF through the MWL program. This work is based upon research supported by the National Research Foundation and Department of Science and Technology. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and therefore the NRF and DST do not accept any liability in regard thereto.

References

All Tables

All Figures

	Fig. 1 Electric field energy ϵ = 1, α = 2.5 × 10-5, γ = 10.
In the text

	Fig. 2 E-field 1D Fourier ϵ = 1, α = 2.5 × 10-5, γ = 10 at .
In the text

	Fig. 3 E-field 2D Fourier ϵ = 1, α = 2.5 × 10-5, γ = 10 at .
In the text

	Fig. 4 Particle histogram ϵ = 1, α = 2.5 × 10-5, γ = 10.
In the text

	Fig. 5 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), E-field 1D Fourier, resolution 1024 × 1024 at .
In the text

	Fig. 6 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), electric field energy, resolution 1024 × 1024.
In the text

	Fig. 7 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), E-field 1D Fourier, resolution 2048 × 2048 at .
In the text

	Fig. 8 Faster beam (ϵ = 1, α = 1.25 × 10-5, γ = 20), electric field energy, resolution 2048 × 2048.
In the text

	Fig. 9 Faster beam, E-field 2D Fourier, resolution 2048 × 2048 at .
In the text

	Fig. 10 Strong beam (ϵ = 10, α = 1.25 × 10-4, γ = 20), electric field energy.
In the text

	Fig. 11 Strong beam (ϵ = 10, α = 1.25 × 10-4, γ = 20 at ), E-field 2D Fourier.
In the text

	Fig. 12 Strong beam (ϵ = 10, α = 1.25 × 10-4, γ = 20), particle histogram.
In the text

	Fig. 13 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, E-field 1D Fourier at .
In the text

	Fig. 14 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, electric field energy.
In the text

	Fig. 15 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, E-field 2D Fourier at .
In the text

	Fig. 16 Weak beam ϵ = 0.1, α = 2.5 × 10-6, γ = 10, particle histogram.
In the text