Physics Background

This page explains what grplinst computes and why. It is meant to be read before you run anything, so that the numbers you get out carry a physical meaning. None of it is needed to make the code compile, but it is needed to interpret the results honestly.

1. Blazar-induced electromagnetic cascades

TeV blazars emit very-high-energy (VHE) gamma rays. On their way to us, these photons do not travel through empty space: the Universe is filled, among other things, by the extragalactic background light (EBL) and the cosmic microwave background (CMB). A VHE photon can collide with a low-energy background photon and produce an electron–positron pair,
\(\qquad \gamma + \gamma_\mathrm{bg} \;\rightarrow\; e^{+} + e^{-} .\)
The pairs are ultra-relativistic (Lorentz factors \(\gamma \sim 10^{6}\)) and inherit the direction of the parent gamma ray, forming a narrow pair beam. Each pair then up-scatters CMB photons through inverse-Compton (IC) scattering,
\(\qquad e^{\pm} + \gamma_\mathrm{bg} \;\rightarrow\; e^{\pm} + \gamma\),
producing secondary gamma rays at GeV energies. Those can pair-produce again, and so on: the result is an electromagnetic cascade that reprocesses TeV power down to GeV energies.

The observable consequences are: a “pair echo”, a spatially extended “pair halo”, and the reprocessed spectrum. They are all sensitive to the intergalactic magnetic field (IGMF), due to the action of the field on the pairs. This makes blazar cascades one of the few probes of magnetic fields in cosmic voids. CRPropa’s EMPairProduction and EMInverseComptonScattering modules implement exactly this chain.

2. The plasma-instability effect

The picture above assumes the pairs lose energy only through inverse-Compton scattering. But a relativistic pair beam streaming through the ionised IGM is a two-stream configuration, and such configurations are generically unstable: the beam excites electrostatic (and other) plasma waves in the background medium, feeding its kinetic energy into those waves instead of into IC photons.

If this beam–plasma instability grows fast enough, it becomes a competitive (or even dominant) energy-loss channel for the pairs. The consequences are considerable:

Whether instabilities actually win over IC cooling is genuinely debated. The growth rate depends on the beam density, its angular and energy spread, the temperature and density of the background plasma, and non-linear saturation effects that linear theory does not capture. The literature spans “instabilities dominate” to “instabilities are irrelevant”, which is precisely why a tool that lets you swap between prescriptions is useful.

3. How grplinst models it: effective cooling

grplinst does not simulate the instability itself. Instead, it represents each literature prescription as an effective cooling term acting on the pairs, exactly like any other continuous energy-loss process in CRPropa. Only electrons and positrons are affected (\(|\mathrm{id}| = 11\)); photons are ignored by the code.

Each model provides a single quantity, a characteristic cooling time \(\tau(E, n_b, n, T, z)\), which the timescale on which a pair of energy \(E\) loses its energy to the instability. The base module turns this into an energy-loss length and applies it during propagation,
\(\qquad \dfrac{\mathrm{d}E}{\mathrm{d}x} = \eta \, \dfrac{E}{c\,\tau}\),
where \(c\) is the speed of light and \(\eta \in [0, 1]\) is a user-defined efficiency factor that lets you scale the effect (for sensitivity tests) or switch it off. For details, see PlasmaInstability::process.

Units convention. The module interprets \(\tau\) as a time in seconds. In fact, internally the code use S.I. units throughout.

4. The three physical inputs

Every prescription is a function of some subset of three environmental inputs, which you supply as separate, composable objects.

4.1 The pair-beam density — Flow

The single most important (and most uncertain) input is the beam number density (\(n_b\)). A Flow object stores properties related to the beam, including the beam luminosity (for normalisation), and the density, through getDensity(energy, position, redshift).

The default estimate, FlowHomogeneous, follows Broderick et al. (2012). It assumes the source injects all of its luminosity \(L\) into the beam and that IC cooling sets the pair lifetime, giving an upper limit of sorts on the beam density (their eq. 7):
\(\qquad n_b(E, z) = \dfrac{L}{2\pi \, \lambda_{\gamma\gamma}^{3} \, \Gamma_\mathrm{IC}} \, \dfrac{1}{E}\),
wherein \(\lambda_{\gamma\gamma}\) is the pair-production mean free path of the parent photon. If no CRPropa pair-production module is attached, the code uses the analytic approximation
\(\qquad \lambda_{\gamma\gamma} \approx 35\ \mathrm{Mpc} \, \dfrac{0.5\,\mathrm{TeV}}{E} \, \left(\dfrac{1+z}{2}\right)^{-4.5}\),
the factor \(0.5\) reflects that a pair member carries about half the parent photon energy. Here \(\Gamma_\mathrm{IC}\) is the inverse-Compton loss rate, which in the Thomson regime is given by,
\(\qquad \Gamma_\mathrm{IC} = \dfrac{4\sigma_\text{T} \, c}{3} \, \dfrac{E}{m_e c^2} \, \dfrac{u_\mathrm{CMB}}{m_e c^2} \, (1+z)^4\).

If you instead attach CRPropa’s EMPairProduction and EMInverseComptonScattering modules to the Flow (via setPairProduction / setInverseCompton), the corresponding actual rates are used in place of the analytic approximations, making the beam-density estimate consistent with the cross-sections and photon fields used in the rest of the simulation. Because the beam density enters the cooling time strongly (see the exponents in Models), the choice of Flow model is often the dominant systematic.

FlowJet1D lets you impose a tabulated \(n_b(r)\) profile instead, and you can sub-class Flow in Python for anything else.

4.2 The ambient density (MediumDensity)

\(n\) is the number density of the background IGM plasma. It sets the local plasma frequency
\(\qquad \omega_\mathrm{pl} = \sqrt{\dfrac{n e^2}{m_e \varepsilon_0}}\,\),
which controls how the medium responds to the beam. MediumDensityHomogeneous returns a constant comoving value scaled to the proper density,
\(\qquad n(z) = n_0 (1+z)^3\).

4.3 The ambient temperature (MediumTemperature)

\(T\) is the temperature of the background plasma. It matters because the thermal spread of the background electrons can quench the instability: a hot medium is harder to destabilise. Models such as Schlickeiser et al. (2012) and Vafin (2018) depend on \(T\) explicitly. MediumTemperatureHomogeneous returns
\(\qquad T(z) = T_0 (1+z)\),
and the base class can convert temperature into a thermal velocity \(v = \sqrt{k_B T / m}\) through getVelocity.

5. Redshift conventions

grplinst follows CRPropa’s convention that a candidate stores its observed (present-day) energy and that the environment is described in comoving terms. Internally each model evaluates the physics in the local frame at the candidate’s redshift \(z\):

Quantity Scaling used in the code
Particle energy \(E_\mathrm{local} = E_\mathrm{obs}\,(1+z)\)
Medium density \(n(z) = n_0\,(1+z)^3\)
Beam density \(n_b(z) = n_b^\mathrm{com}\,(1+z)^3\)
Medium temperature \(T(z) = T_0\,(1+z)\)

The step length is likewise de-redshifted,
\(\qquad \mathrm{d}x = \dfrac{\Delta x}{1 + z}\),
before the loss is applied. You normally do not apply these factors yourself; instead, you specify present-day comoving values and the module handles the scaling.

6. What this model is not

Please keep the following limitations in mind:

Used with these caveats in mind, grplinst is a controlled way to ask “how much could plasma instabilities change my cascade result, and how much does the answer depend on which prescription I believe?”