Simulating the Enrichment of Fossil Radio Electrons by Multiple Radio Galaxies

Abstract

We simulate the evolution of relativistic electrons injected into the intracluster medium by five radio galaxies. We study the spatial transport and the emission properties of the injected radio plasma over a ∼5 Gyr period, and the sequence of cooling and re-acceleration events experienced by electrons, using a Lagrangian approach joined with a numerical method to model the evolution of momentum spectra of relativistic electrons. When compared with electrons injected by shock waves, electrons injected by radio galaxies (here limited to a single injection event) in our tests are unable to fuel large ∼Mpc-sized radio relics with fossil electrons, as required by current theoretical models, while electrons previously seeded by other shocks can do this. On the other hand, the combination of seeding from radio galaxies and of re-acceleration events from plasma perturbation can produce detectable, small-scale, and filamentary emissions in the proximity (≤100–200 kpc) of radio galaxies.

1. Introduction

Radio galaxies are a spectacular example of how accretion discs surrounding supermassive black holes (SMBH) couple radically different astrophysical temporal and spatial scales, i.e., from the innermost accretion orbits of SMBH (≤${10}^{-3}\mathrm{pc}$) to the expanding jets, recently observed out to ∼2 Mpct from their host galaxy, e.g., [1]. In galaxy clusters and groups, radio galaxies represent a prominent and visible reservoir of magnetic fields and non-thermal particles, e.g., [2,3]. The recent decade or so suggested that a volume-filling distribution of fossil relativistic electrons is often required in order to explain the observed radio power of radio relics, i.e., elongated and polarised radio sources that are typically associated with merger shock waves, e.g., [4] for a review. Their low Mach numbers render direct Fermi I acceleration from the thermal pool challenging, e.g., [5,6,7,8]. The presence of a volume filling distribution of “fossil” ($\gamma \le {10}^3$, where $\gamma$ is the Lorentz factor) has been suggested to explain the formation of radio halos via Fermi II turbulent re-acceleration, e.g., [9] for a recent review. Even more recently, it has been invoked to explain the large-scale emission detected at the extreme periphery of a few interacting clusters of galaxies, e.g., [10,11] or on much larger scales than the classic radio halos, e.g., [12,13].

Tails of radio-emitting plasma from radio galaxies mix with diffuse radio emission from the ICM, e.g., [14,15]. Recent low-frequency observations are detecting complex morphologies of remnant plasma injected by radio galaxies, and in different stages of mixing their non-thermal content with their surrounding environment, e.g., [16,17,18,19,20].

The new generation of radio telescopes (e.g., LOFAR, MWA, ASKAP, MEERKAT) is capable of producing an unprecedented view of this phenomenon. At the same time, there presently are only a few simulations designed to connect diffuse re-acceleration mechanisms in the ICM to the ageing of radio lobes.

Simulating the multi-scale and long-term interplay between radio jets nearby and cluster weather still is a challenge, even for modern numerical simulations, which typically have to focus on specific ranges of spatial/temporal scales and physical aspects of the problem. This is particularly demanding for cosmological simulations, in which properly resolving the accretion regions of supermassive black holes is impossible. To give a few examples, Xu et al. [21,22] included magnetised outflows from radio galaxies in cosmological ENZO simulations, and studied the build-up of cluster magnetic fields from the injection of individual AGN jets, but there was no physical coupling between accretion and feedback powers. Mendygral et al. [23] used an Eulerian Magneto Hydrodynamical approach to follow in detail the effect of “cluster weather” on radio lobes, by injecting magnetised jets into a galaxy cluster extracted from a different (Smoothed Particle Hydrodynamics, SPH) cosmological simulation. Bourne and Sijacki [24] used the moving mesh AREPO code to realise high-resolution simulations of jets and predict X-ray and radio properties but had to assume concerning the post-processing distributions of magnetic fields and relativistic electrons.

More recently, the complex semi-analytic modelling of the radio emission from jets leading to Fanaroff Riley I and II systems have been presented by Turner et al. [25], and this model has been incorporated in simulations of jets released in the atmospheres of clusters of galaxies, extracted from a large set of SPH cosmological simulations by Yates-Jones et al. [26]. As for similar previous work, e.g., [23], the usage of a composite resimulation allows us to greatly increase the physical realism that can be modelled in the interaction of jets with their environment; however, the same strategy prevents following the evolution of lobe materials for timescales longer than ∼100–200 Myr since the injection, as it neglects the long term effect of gas self-gravity, as well as further cluster perturbation from outside the resimulation region.

Our series of recent works on this topic [27,28] aims at producing a complementary view on the above studies by implementing the macroscopic effects of radio jets at a coarser spatial resolution, but also following the spatial and spectral energy evolution of the relativistic particles over billions of years after their injection, in fully cosmological MHD simulations. First, applications of our simulations to real radio sources have recently been presented in [20,29,30].

In this paper, we apply the same framework to study the spatial transport of electrons injected by five radio sources, all with different feedback powers (in our case, initially assigned as thermal and magnetic power). These radio galaxies are activated at the same time, within an already formed ∼${10}^{14}{M}_{\odot }$ cluster of galaxies. With this new analysis, we mainly want to address the following questions:

• Under which conditions can radio galaxies fill the ICM with relativistic electrons to a significant degree?
• Is the predicted distribution of re-accelerated fossil relativistic electrons suitable to produce diffuse radio emission, and is this compatible with radio observations?

Our paper is structured as follows: In Section 2, we present the cosmological simulations and all numerical methods employed in this paper. In Section 3, we give our main results, first focusing on the evolution of electrons seeded by radio jets (Section 3.1), and then by comparing with the evolution of electrons injected by merger shocks (Section 3.2). Our conclusions are given in Section 4.

Throughout this paper, we use the following cosmological parameters: $h=0.678$, ${\Omega }_{\mathsf{\Lambda }}=0.692$, ${\Omega }_{\mathrm{M}}=0.308$, and ${\Omega }_{\mathrm{b}}=0.0478$, based on the results from [31].

2. Methods & Material

In order to make this section as compact as possible, we highlight only the key technical aspects of our simulations, while referring the reader to our previous publications [27,28] for the technical details.

2.1. Cosmological Simulations & Feedback from Active Galactic Nuclei

We produced cosmological, adaptive mesh refinement ENZO-MHD [32] simulations using the MHD solver with a Lax–Friedrichs Riemann solver to compute the fluxes in the Piece-wise Linear Method, combined with the Dedner cleaning method [33] implemented by [34]. We start from a simple uniform magnetic field at $z=50$ with a value of $B_0=0.1\mathrm{nG}$ in each direction. This simplistic model has been shown to produce radio signatures which are in tension with a few radio observations of the low-density cosmic web, making tangled initial magnetic fields preferred [35], yet the differences expected within clusters of galaxies are negligible, owing to the effect of small-scale dynamo amplification [36,37,38].

The total simulated volume is of (50 Mpc$/h$)${}^3$ and is sampled with a root grid of ${128}^3$ cells and dark matter particles. We further use four additional nested regions with increasing spatial resolution (up to $\mathsf{\Delta }x=24.4$ kpc$/h$) within a cubic (5.2 Mpc)${}^3$ sub-volume, where a $M_{100}\approx 1.5·{10}^{14}{M}_{\odot }$ cluster forms. During run-time, two additional levels of mesh refinement are added using a local gas/DM overdensity criterion ($\mathsf{\Delta }\rho /\rho \ge 3$), allowing the simulation to reach a maximum resolution of ≈$6\mathrm{kpc}/\mathrm{h}=8.86$ kpc. As a result of our nested grid approach, the mass resolution for dark matter in our cluster formation region is of $m_{\mathrm{DM}}=2.82·{10}^6{M}_{\odot }$ per dark matter particle, for the highest resolution particles that are used to fill the innermost AMR level since the start of the simulation. No radiative gas cooling or other galaxy-formation-related physics (e.g., star formation or feedback from supernovae) was included, besides the single jet feedback episode from eight different radio galaxies, explained in the next Subsection.

The thermal and magnetic properties of this cluster have been described in detail in Vazza et al. [27] and Vazza et al. [28], and are largely independent of the specific kind of AGN feedback. The galaxy cluster has a virial mass $M_{100}\approx 1.1·{10}^{14}{M}_{\odot }$ at $z=0.5$, which grows to $M_{100}\approx 1.5·{10}^{14}{M}_{\odot }$ by $z=0$ through the accretion of subsequent minor mergers, the most prominent ones being at $z\sim 0.3$ and $z\sim 0.1$ and mostly directed along one of the diagonals of the computational box.

At $z=0.5$, we identify the five most massive halos in the innermost high resolution $5.2^3{\mathrm{Mpc}}^3$ volume, in the ∼${10}^9-{10}^{11}{M}_{\odot }$ range, and place an SMBH particle at the centre of each of them. We assume a black hole mass of $M_{\mathrm{BH}}=0.01M_{200}$, where $M_{200}$ is the total (gas+DM) mass of each halo.1

SMBH particles have initial masses in the $8.1·{10}^6$–$8.7·{10}^8{M}_{\odot }$ range, and after injection, we estimate their gas accretion rate, based on the usual Bondi–Hoyle formalism:

$${\dot{M}}_{\mathrm{BH}}=\frac{4\pi {\alpha }_B{G}^2{M}_{\mathrm{BH}}^2\rho }{{(v_{\mathrm{rel}}^2+c_s^2)}^{3/2}},$$

(1)

in which $\rho$ is the local gas density, $v_{\mathrm{rel}}$ is the SMBH velocity relative to the gas, $c_s$ is the sound speed at Bondi radius (which we assume to be relative to a fixed ${10}^6\mathrm{K}$ temperature, considering that our simulation lacks the resolution and physical ingredients to predict the temperature in this region, also considering the lack of radiative cooling), and $\rho$ is the local gas density. ${\alpha }_B$ is a free normalisation factor, often employed in simulations incapable of resolving the Bondi radius, e.g., [40,41,42], which we here simply leave to ${\alpha }_B=1.0$ for all simulated SMBHs.

As shown in

Table 1. Main parameters of SMBH and jet/radio power associated with the jets they emit (in which case the power is the sum of the power in the two jets).

ID	${\mathit{M}}_{\mathbf{BH}}$	${\dot{\mathit{M}}}_{\mathbf{BH}}/{\dot{\mathit{M}}}_{\mathbf{Edd}}$	${\mathit{P}}_{\mathit{j}}$	${\mathit{L}}_{140}$
	[${\mathit{M}}_{\odot }$]		[erg/s]	[erg/s/Hz ]
A	$8.7·{10}^8$	0.0026	$1.4·{10}^{43}$	$3.4·{10}^{32}$
B	$3.5·{10}^7$	0.0004	$1.0·{10}^{41}$	$1.6·{10}^{31}$
C	$1.4·{10}^7$	0.00004	$1.1·{10}^{40}$	$5.3·{10}^{30}$
D	$1.3·{10}^7$	0.00007	$1.3·{10}^{39}$	$2.9·{10}^{31}$
E	$8.1·{10}^6$	0.00004	$7.1·{10}^{38}$	$1.9·{10}^{31}$

, our simulated SMBHs have a range of accretion powers from ∼$0.26\%$ to ∼$0.0004\%$ of the Eddington limited accretion rate. In this sense, they all represent a relatively quiet population of SMBHs in the epoch when the jets are launched in the ICM. The quietness of our population can be partially ascribed to the absence of large mass accretion rates, which in turn follow from the absence of radiative cooling (which tends to increase density and reduce the sound speed/temperature in Equation (1)); however, this ensures that our results can represent rather average sources in clusters of galaxies.

The bolometric luminosity of each SMBH is the standard $L_{\mathrm{BH}}={ϵ}_r{\dot{M}}_{\mathrm{BH}}{c}^2$, where we assumed the standard values in the literature of ${ϵ}_r=0.1$ for the radiative efficiency of the SMBH. This is used to compute the jet power released in the feedback stage (which only last, by construction, one root grid timestep, i.e., ∼$32\mathrm{Myr}$) as $P_j={ϵ}_r{ϵ}_{\mathrm{BH}}{\dot{M}}_{\mathrm{BH}}{c}^2$, where ${ϵ}_{\mathrm{BH}}$ is the factor that converts the bolometric luminosity to the thermal feedback energy (fixed to 0.05, also in line with what is done in most cosmological simulations).

Jets are introduced at the same time ($z=0.495$) for all SMBH particles, following the same recipes presented in [28].

SMBH particles inject $10\%$ of the total feedback energy as magnetic field energy (with the remaining 90% injected as thermal energy), by introducing two magnetic loops ($2\times 2$ cells at the highest resolution level), located at $\pm 1$ cell along the z-direction from each SMBH, i.e., at $\pm 8.86\mathrm{kpc}$ separation). Since the thermal feedback energy gets isotropically spread over a larger volume (i.e., 27 cells), the typical magnetic energy within the active jets is a factor of a few larger than the kinetic and thermal energy within the same region, hence our jets are initially effectively magnetic energy dominated, at least in the first ≤$100\mathrm{Myr}$ since their creation.

Even though imposing a fixed jet orientation along the z axis for all our SMBH jets is artificial, in the long term, the initial direction of jets is observed to play no role. We also tested in Vazza et al. [28] that different choices in the initial direction of jets produce no significant differences in the statistics of the thermodynamical properties of the ICM and of the related radio emission.

Our approach is similar to early work by [21,43], who initialised toroidal magnetic fields on opposite sides of the SMBHs, showing that this procedure naturally led to self-collimating outflows and a supersonic expansion in the cluster core. It shall be noticed that our chosen $2\times 2$ cell configuration for the initial magnetic field loops is not enough to properly resolve the field loops in the model of [43]; however, this is not a source of concern because the focus of our work is more on the energetics and distribution of tracer particles, and not on the details of jet morphology.

The distribution of the simulated radio power at 140 MHz, measured ∼$32\mathrm{Myr}$ since the start of the jet active stage and computed after the injection of relativistic electrons and with the synchrotron approach explained in the next Section, is plotted in Figure 1 as a function of the X-ray luminosity within ≤$500\mathrm{kpc}$ of each source, to compare with the recent LOFAR-HBA observations by Pasini et al. [44]. This suggests that the range of radio powers spanned by our sample is well within the one of observed real sources, even if our sample also includes non-central ones. Recent work by [45] used radio and X-ray observations of six rich clusters of galaxies to study the role of AGN feedback by non-central extended radio galaxies, and concluded that while cluster-central sources are the dominant factor to balance cooling over the long-term, non-central sources also have a significant impact.

2.2. Injection, Tracking and Energy Evolution of Relativistic Electrons from Radio Galaxies

We track the relativistic electrons with a post-processing Lagrangian scheme, running the CRATER code [46] over all snapshots of our simulation (221 timesteps, with a time resolution of about $\mathsf{\Delta }t\sim 32\mathrm{Myr}$).

We injected and propagated ∼$5·{10}^4$ particles, assigning them according to the local gas density distribution at ≤$8\mathrm{kpc}$ from each active jet after we discretized the tracer mass resolution to $m_{\mathrm{trac}}=8·{10}^5{M}_{\odot }$. The various physical quantities of interest (gas density, velocity, divergence, vorticity, temperature, and magnetic field intensity) are assigned to the tracers using a cloud-in-cell (CIC) interpolation method, while shocks were detected based on temporal temperature discontinuities.

Electron spectra for each tracer were initialised so as to reproduce the observed radio spectra of radio galaxies, through the Continuous Injection-Off model, e.g., [47], with initial normalisation following from assuming a $0.5\%$ fixed fraction between the injected relativistic electrons and the thermal proton density in the jet region, e.g., [23]. For tests of this approach in our model, we refer the reader to Vazza et al. [28].

After initialisation, we solve for the time-dependent diffusion-loss equation of radio-emitting particles with the parallel ROGER solver2, using momentum bins equally spaced in $\log \left(\mathrm{p}\right)$, in the $p_{\min }\le p\le p_{\max }$ momentum range (where $p=\gamma m_{e}v/\left(m_{e}c\right)$ is the normalised momentum of CR electrons). In all production runs, we used $p_{\min }=5$ and $p_{\max }={10}^6$ and $dlog\left(p\right)=0.05$.

We consider radiative losses, Coulomb losses and expansion (compression), as well as the acceleration from diffusive shock acceleration (DSA). The shock kinetic energy flux that we assumed to be converted into the acceleration of cosmic rays is: ${\mathsf{\Psi }}_{\mathrm{CR}}={\xi }_e\eta \left(\mathcal{M}\right){\rho }_u\left(V_s^{3}d{x}_t^2\right)/2$, in which ${\rho }_u$ is the pre-shock gas density, $V_s$ is the shock velocity and the combination ${\xi }_e\eta \left(\mathcal{M}\right)$ gives the CR electrons acceleration efficiency.

The latter is a combination of the energy going into cosmic rays, $\eta \left(\mathcal{M}\right)$ (for which we use the convenient polynomial approximation by [48]), and the electron-to-proton acceleration rate, ${\xi }_e$, which we fix by requiring an equal number density of supra-thermal cosmic-ray electrons and protons above the injection momentum, ${\xi }_e={(m_p/m_e)}^{(1-{\delta }_{\mathrm{inj}})/2}$, consistent with [6]. For simplicity, we neglect dependencies on the shock obliquity, e.g., [49], on the basis that our previous work showed that the majority of internal shocks in the ICM are quasi-perpendicular and hence suitable for electron acceleration, e.g., [50,51].

The injection momentum of electrons $P_{\mathrm{inj}}$ is linked to the thermal momentum of particles, i.e., $P_{\mathrm{inj}}=\xi P_{\mathrm{th}}$ ($P_{\mathrm{th}}=\sqrt{2k_b{T}_d{m}_p}$), measured locally for each Lagrangian tracer.

The working surface associated with each tracer, $d{x}_t^2$, is adjusted at run-time as $d{x}_t^3=d{x}^3/n_{\mathrm{tracers}}$, where $d{x}^3$ is the volume initially associated with every tracer and ($n_{\mathrm{tracer}}$ is the number of tracers in every cell).

Supra-thermal electrons are injected by shocks with a power-law momentum distribution, e.g., [52], $Q_{\mathrm{inj}}\left(p\right)=K_{\mathrm{inj}}{p}^{-{\delta }_{\mathrm{inj}}}{(1-p/p_{\mathrm{cut}})}^{{\delta }_{\mathrm{inj}}-2}$ where the initial slope of the input momentum spectrum, ${\delta }_{\mathrm{inj}}$, follows from the standard DSA prediction, ${\delta }_{\mathrm{inj}}=2({\mathcal{M}}^2+1)/({\mathcal{M}}^2-1)$. $p_{\mathrm{cut}}$ is the cut-off momentum, which we set for every shocked tracer as the maximum momentum beyond which the radiative cooling time scale gets shorter than the acceleration time scale, ${\tau }_{\mathrm{DSA}}$:

$${\tau }_{\mathrm{DSA}}=\frac{3D\left(E\right)}{V_s^2}·\frac{r(r+1)}{r-1},$$

(2)

in which r is the shock compression factor and $D\left(E\right)$ is the electron diffusion coefficient as a function of energy, which is poorly constrained as it depends on the microphysical conditions of the turbulent plasma. However, all plausible choices of $D\left(E\right)$ in Equation (2) give an acceleration timescale which is many orders of magnitude smaller than the typical cooling time of radio-emitting electrons, e.g., [5]; hence we can assume the momentum distribution at injection can be assumed to follow a power law within our momentum range of interest.

The normalisation factor for shock-injected electrons, $K_{\mathrm{inj}}$, follows from equating the cosmic ray energy flux crossing each tracer volume element, and the product between the total energy of cosmic rays ($E_{\mathrm{CR}}$) advected with a post-shock velocity ($v_d$): ${\mathsf{\Psi }}_{\mathrm{CR}}d{x}_t=v_d{E}_{\mathrm{CR}}$ ($v_d$ is the post-shock velocity), similar to [6].

Additionally, we model re-acceleration by shocks, e.g., [5,53,54]: $N\left(p\right)=({\delta }_{\mathrm{inj}}+2)·p^{-{\delta }_{\mathrm{inj}}}{\int }_{p_{min}}^p{N}_0\left(x\right)x^{{\delta }_{\mathrm{inj}}-1}dx$ where ${\delta }_{\mathrm{inj}}$ is the local slope within each energy bin and $N_0\left(p\right)$ is the input spectrum of electrons before the re-acceleration.

We also include Fermi II re-acceleration via stochastic interaction with diffusing magnetic field lines in super-Alfvenic turbulence. The re-acceleration is computed following the Adiabatic Stochastic Acceleration (ASA) model [55], and it depends on the amplitude of the local turbulent velocity, $\delta V_{\mathrm{turb}}$, measured from the gas vorticity, i.e., $\delta V_{\mathrm{turb}}=|\nabla \times \overrightarrow{v}|L$, for which we used the same fixed reference scale of $L\approx 27$ kpc to compute vorticity via finite differences (i.e., 3 cells on the high-resolution mesh). The acceleration term for electrons is $dp/dt\approx p/{\tau }_{\mathrm{ASA}}$, for which ${\tau }_{\mathrm{ASA}}\left(p\right)=p^2/4D_{pp}$, which is in the ${\tau }_{\mathrm{ASA}}\sim 0.1$–$1\mathrm{Gyr}$ range for the typical dynamics of these systems. The diffusion coefficient in momentum space in this scenario is (see [11] for the derivation) $D_{pp}\approx \left(48F_{\mathrm{turb}}{p}^2\right)/\left(c{\rho }_{ICM}{V}_A\right)$, which depends on the solenoidal turbulent energy flux (conserved in the Kolmogorov model of turbulence): $F_{\mathrm{turb}}=\frac{\rho \delta V_{\mathrm{turb}}^3}{2L}$. For the sake of brevity, we omitted here the full derivation of the parameters above, in connection with the quantities recorded by our Lagrangian tracers, which can be found in longer detail in [28].

Finally, to account for the missing efficient amplification of magnetic fields (due to the under-resolved role of a small-scale dynamo) and avoid the unphysical amount of turbulent reacceleration in the ASA scheme (since ${\tau }_{\mathrm{ASA}}\propto B$) we assign each tracer a re-normalised estimate of the local magnetic field, in which we take the maximum between the magnetic field value directly produced by the MHD calculation, and the one assuming turbulent dynamo amplification, for a constant ${\eta }_B=2\%$ conversion efficiency between the local turbulent kinetic energy flux of the solenoidal component and the energy density of magnetic fields (see [28] for detailed tests on the outcomes of this approach).

The synchrotron emission from the electrons was computed as explained in Vazza et al. [28], i.e., using the (fast) fitting procedures introduced by [56], for five radio frequencies (50, 120, 610, 1400, and 5000 MHz). Our ROGER electron solver is designed to solve multiple populations of electrons under different physical scenarios. Hence we could simulate the outcome of different possible scenarios for the origin and the evolution of radio-emitting electrons, in particular:

• Electrons seeded by radio jets: we followed the evolution of electrons injected by our radio galaxies at $z=0.5$, solely under the influence of loss processes (e.g., radiative processes, Coulomb collision, ionisation losses, and adiabatic changes, model “C”), by including the additional injection of new electrons by diffusive shock acceleration, as well the re-acceleration by DSA on the previously injected population of electrons (model “CS”), or by additionally including also the Fermi II re-acceleration by solenoidal turbulence (model “CST”). This set of models followed $\approx 3·{10}^4$ tracers, all initially placed within the jet launching regions of our five radio galaxies.
• Electrons seeded by merger shocks: we initialised pools of tracers in the simulation at $z=0.4$ (i.e., well after the jet activity has ended for all considered radio sources), assigning them according to the gas density profile and giving them an initially negligible content of relativistic electrons. We let them evolve according to all loss processes, shock injection, and shock/turbulent re-acceleration. This is similar to the previous “CST” physical model, with the important difference that the electrons are only seeded by merger/accretion shock waves after the short active stage of jets. This model is meant to check the effectiveness of the multiple-shock scenario explored in Inchingolo et al. [57], in which a large fraction of the radio power from relics comes from shock re-acceleration of electrons previously injected by older structure formation shocks. In other words, this second model is meant to quantify the relative importance of fossil electrons injected by shocks, compared to fossil electrons injected by radio galaxies on longer timescales. This second set of models followed $\approx 2·{10}^5$ tracers, placed in the cluster at $z=0.4$ following the cluster density profile and only limited to $T\ge {10}^7\mathrm{K}$ regions.
• Prompt injection of electrons by shocks: we included a simplistic scenario in which we only compute the prompt injection of electrons, and their radio emission, at a single time-step and based on the DSA model outlined above, i.e., with injection efficiency scaling with the Mach number and normalisation depending on the shock kinetic energy flux. Otherwise, the electron spectra are deleted from one-time step to the next. This last model is just meant to compare with the standard approach to predict radio emission from shocks in the simulated ICM, in which single snapshots of simulations are used, assuming quasi-stationary shock conditions and neglecting the presence of fossil electrons, e.g., [58,59,60]. This last scenario was run on top of the same distribution of tracers used for the previous model.

3. Results

3.1. The Evolution of Electrons Seeded by Radio Jets

We begin by analysing the evolution of the three-dimensional location of the electrons seeded by radio jets. The time sequence shown in Figure 2 shows the evolving gas density (green colours), mass-weighted gas temperature (blue) and the projected location of all electrons injected by radio jets (red) starting from $z=0.5$ and, for a few interesting evolutionary steps, down to $z=0.032$.

Electrons rapidly expand into the cluster atmosphere, and get progressively dispersed and mixed on larger scales by being entrained by turbulent motions or large-scale advection flows induced by the accretion of cluster satellites.

With the exception of the central, dominant SMBH (which corresponds to the BCG in this system), the initial jet power from the other galaxies is not sufficient to push electrons to large distances from their host SMBH, and electrons first settle into bubble-like distributions, to be later dispersed into the main cluster while they follow the accretion of their host galaxy. The relative motions between the injected electrons and the ICM are large enough to thoroughly disperse the CRe throughout the cluster volume. By the end of our simulation, a large fraction of the cluster area (i.e., the cluster volume seen in the projection) appears to be filled by the electrons seeded by a single burst of our radio galaxies (we will quantify the volume-filling factor of these electrons in the next Section).

The fuelling of electrons in the innermost cluster region is thus largely dominated, at least for the first ∼$0.8\mathrm{Gyr}$ ($z\sim 0.5$–$0.4$) by the contribution of the central radio galaxy only. Figure 3 shows the projected radio emission at 50 MHz in the first stage of radio lobe expansion, i.e., from $z\approx 0.5$ to $z\approx 0.4$, shortly after which the remnant emission even from the most powerful radio galaxy becomes invisible (the rest of the evolution of the radio detectable emission is shown in Figure 8, in which we also add the emission from shock injected electrons). To estimate which emission would be detectable, we considered a fixed luminosity distance of $d_L=132\mathrm{Mpc}$ for all snapshots, in which case our simulated pixel size corresponds to the resolution beam of a LOFAR Low Band Antenna (LBA) at 50 MHz, of ≈$8\mathrm{kpc}$, considering a beam of $\theta ={12.5}^{″}$, and a sensitivity of $\sigma =5.7·{10}^{-4}\mathrm{Jy}/\mathrm{beam}$, as in [13].

After the initial bright stage of jet emission, only a tiny fraction of the emission from the injected electrons remains visible at radio frequencies (red contours). Already after 0.5 Gyr since the injection (z = 0.436 in the images), the only visible emission is from the stretched remnant lobe structures, which have begun to mix in the innermost cluster atmosphere. This mixing occurs mostly in the vertical direction in the innermost cluster atmosphere, thus carrying the memory of the initial jet orientation. The detectable structures are mostly filamentary and patchy, with ≤$500\mathrm{kpc}$ length, until $z=0.38$, after which there is no detectable radio emission for a long while. However, the progressive mixing and dispersal of the different families of electrons—now become a “fossil” population—continue; later on (e.g., $z\approx 0.37$, $0.26$, $0.21$ and in Figure 8) patches of electrons that are crossed by shock waves, or are intersected by turbulence, get re-accelerated and become visible again. All these events are triggered by accretion events of substructures, which mostly happens along a diagonal in the image, and often release pairs of weak ($\mathcal{M}\le 3$) shock waves, and stir new subsonic turbulence in the main cluster centre.

The electron momentum spectra and the radio spectra (of all tracers located within a $r\le 1\mathrm{Mpc}$ radius from the moving cluster centre) at different epochs, shown in Figure 4, show that only when Fermi I and Fermi II terms are included (models CS and CST) do the electron spectra remain energetic enough to produce significant radio emission, typically with very steep ($\alpha \le 1.5$–2) radio spectra. Only at late times, turbulent re-acceleration is important to increase the budget of low-energy electrons and their low-frequency emission (CST). The net effect of turbulent re-acceleration is to systematically increase the budget of low energy electrons ($p\le {10}^3$), as well as the low-frequency radio emission ($\nu \le 610$ MHz), beyond the effect on shock re-acceleration alone, as shown by the small excess of the solid lines compared to the dotted lines (otherwise always superimposed) in Figure 4. However, it is difficult to state in a quantitative way the relative effects of turbulent re-acceleration in comparison with shock re-acceleration: we studied this in detail in Vazza et al. [27] and determined that this depended on the epoch of observation and on the specific dynamical history experienced by tracers, with a range of extra energy below $p\le {10}^3$, which can range from ∼2 to ∼10 depending on the time since the last shock event experienced by tracers.

We notice that, especially in the initial stage of lobe expansions, where the further injection of CRe from shocks is small, the fraction of the CRe population in the radio band may slightly change for a higher assumed initial injection fraction of CRe electrons in jets. However, our choice ($0.5\%$ of the number density of thermal particles) already is at the high end of the range typically assumed for these sources, e.g., [23,61], and moreover, a ∼10 higher normalisation of the initial CRe density would make our central BCG extremely powerful, e.g., nearly at the top of the distribution of radio power recently measured by Pasini et al. [44], which would be odd for a source in such a small mass cluster of galaxies.

Although the electrons injected by the central BCG clearly dominate the radio emission at the start of the simulation, its relative contribution when re-acceleration processes become dominant changes over time. Figure 5 shows the trajectories of particles initially belonging to different radio galaxies (with different colours) and their relative contribution to the total radio emission from the cluster (here without observational cuts for the CST model). While the emission is initially dominated by the powerful “A” source, the relative contribution of the electrons in its lobe mixes with time with the electrons coming from other radio sources. At later times also the pool of electrons released by sources “D” and “B” become the bright emitters in the field since their particles get crossed by shock waves. This shows that the seeding of electrons from sources other than the central powerful radio galaxy likely has a key role in enriching the ICM with additional fossil electrons, and it thus complements our previous results, which were limited to the impact of a single radio source in the cluster [27,28]. Moreover, detectable radio emission features can arise from the blending of the initially distinct population of radio electrons, similar to what has been proposed to explain a few puzzling steep spectrum radio sources [29]. However, the size of the detectable regions is very small, i.e., ≤100 kpc at all epochs for $z\le 0.4$ (see Figure 8). This means that in all our tests so far, it appears implausible that the much larger extent of giant radio relics (∼$2\mathrm{Mpc}$ in the largest cases, e.g., [62,63]) can be efficiently filled by a uniform pool of fossil electrons, and thus can alone provide the background of seed electrons for Fermi I re-acceleration in these objects. On the other hand, the seeding of electrons from radio galaxies appears a more plausible channel for refilling the central regions of classic radio halos, as already partially supported by direct observations, e.g., [17,64]. Of course, this result has to be considered only as tentative so far, given the single jet events we could simulate so far.

3.2. Shock Injection of Electrons vs. the Injection by Radio Jets

Next, we compare with the result for electrons seeded by merger shocks, in which we track the additional contribution from electrons injected by merger shock waves starting from $z=0.4$, to have a more conservative view of what shocks driven by accretion alone (i.e., without considering shocks triggered by AGN events) can do.

Even if this includes more CR electrons in our modelling, it still represents an underestimate, because we only track here the contribution from merger shocks in the innermost cluster regions, and neglect the injection of electrons from shocks in more peripheral regions. Indeed, given the large number of tracers to follow with our electron solver, we initialised ≈$2·{10}^5$ particles limited to $T\ge {10}^7\mathrm{K}$ cells (approximately within ≤$R_{500}$ for this system at that epoch), which is the environment in which merger shock develops and that is also expected to dominate the injection of CR in clusters, e.g., [65]. Therefore, if even more peripheral shocks would be accounted for in the seeding of electrons, the dominance of shock-seeded electrons over radio jet-seeded electrons would only increase. On the other hand, it is still possible that the time-integrated effects of seeding by radio jets, considering their entire duty cycle even from high redshift, can produce volume-filling factors approaching the ones by electrons seeded by merger shocks.

The volume spanned by merger shocks is a significant fraction of the cluster volume, and hence already after ∼1 Gyr, since this mechanism is allowed in our simulation, a large portion of the innermost ICM is filled with some amount of cosmic ray electrons, as shown by Figure 6 which displays a much larger area covering factor of electrons, compared to the corresponding maps for the seeding from radio jets (last two panels of Figure 2).

We computed the 3-dimensional and 2-dimensional filling factors for our electron tracers released from the different mechanisms as a function of time. Figure 7 shows the evolution of the fraction of volume, or projected area, of radial shells (centred in the moving cluster density peak) filled by at least one simulated tracer. The two trends are opposite and clear. Electrons seeded by radio jets (and in particular by the central one, corresponding to the BCG) occupy a large fraction of the central cluster volume right after their injection epoch, and get advected to a larger volume, and their volume filling factor drops to ∼0.1–1% in the innermost $1{\mathrm{Mpc}}^3$ region, at most epochs. Conversely, electrons seeded by shocks are initially deposited preferentially in outer radial shells in the cluster, where shocks are formed, but then steadily increase their volume filling factor over time, reaching a ∼3–30% filling in the innermost $1{\mathrm{Mpc}}^3$ already at $z\sim 0.23$, i.e., ∼$1.5$ Gyr since they are first injected in our model. More relevant for the potential radio observations which can target such low energy electrons, the filling factors for the projected surface are higher in both cases, but again the filling factor of electrons seeded by shocks is much higher at later times, i.e., ≥5–10 more projected surface is filled with fossil electrons seeded by shocks at all radii, for $z\le 0.23$.

Again, the latter number is probably underestimated by neglecting the entire injection process of CR electrons by structure formation shocks, e.g., [66].

Since seed cosmic-ray electrons are spread over a larger volume, they also have more frequent chances of being re-accelerated by subsequent shocks sweeping the same system and producing radio-detectable emission, as shown by the sequence of radio images in Figure 8, in which we applied the same observational cuts of Figure 3. While fossil electrons released early on by radio jets are sometime re-accelerated by shocks and turbulence at late redshift, the detectable emission covers a wider area if electrons injected by shocks are included. In the latter case, they more naturally give rise to a few ∼500 kpc large relics in our simulation, with a more uniform surface brightness and radio spectra. Considering that here we can only track the evolution of particles first injected at $z=0.4$ in the innermost cluster regions, it is quite plausible to expect that the inclusion of all past injection events of electrons should produce even larger scale—and possibly smoother—radio emission regions.

A final comparison of the long-term radio evolution of emission from electrons seeded by merger shocks and the one by radio jets is given in Figure 9. Following from the previous trends of filling factors, while the contribution from electrons seeded by radio jets overall decreases with time, even in the presence of re-acceleration events, the one from electrons seeded by shocks overall increases with time, and by the latest epochs in the simulation, the total radio emission from the cluster is ∼${10}^2$ larger at 50 MHz if electrons injected by shocks are considered. Given the higher frequency of re-acceleration events experienced by the volume-filling population of shock-injected electrons, the emission from the latter is, at most epochs, significantly flatter than the one coming from re-accelerated electrons from radio galaxies.

Finally, the right panel also shows the radio emission only for the prompt injection of electrons by shocks (thick dashed lines). While at the beginning of the sequence, the contribution from freshly injected electrons dominates the total budget, this becomes, over time, secondary with respect to the effect of shock re-acceleration. Limited to this particular system, and its history of weak internal shocks ($\mathcal{M}\le 3$), this is in line with the frequent finding that radio relic emission observed in many clusters of galaxies is better explained by some form of re-acceleration of fossil electrons, e.g., [5,66,67,68], which in this case appear to be most likely a multiple shock scenario, e.g., [57]. On the contrary, observed radio relic emission compatible with the “simple” prompt seeding of electrons by shocks, based on DSA, appears to be compatible only with a very few objects observed in the radio band, e.g., [69].

4. Discussion & Conclusions

In this work, which is a follow-up study after our recent Vazza et al. [27] and Vazza et al. [28] work, we presented a first assessment of the relative role of multiple radio jets, and of merger shock waves, in (a) filling the ICM with fossil relativistic electrons, and (b) giving rise to detectable radio emission, based on state-of-the-art models of particle (re)acceleration.

Based on the simulation of a single system, our most significant conclusions can be summarised as follows:

• The seeding of electrons from sources other than the central powerful radio galaxy has a significant role in enriching the ICM with an additional amount of fossil electrons. During the evolution, powerful enough radio emission features can be produced from populations of electrons seeded by peripheral galaxies, which also blend and mix over time.
• Right after the active stage of radio jets, the remnant radio plasma typically dominates the radio emission at a ≤$0.5$ Mpc distance from sources, even up to ∼0.5–1 Gyr since its first injection, but only if quasi-continuous re-acceleration events (e.g., frequent weak shocks or turbulence) are active on particles. The emission from re-accelerated electrons seeded by radio galaxies can produce detectable emission at any time, in presence of re-acceleration events, but only leading to small (≤100 kpc in our test) and often filamentary radio features with steep radio spectra ($\alpha \ge$ 1.5–2). This makes it hard for radio galaxies alone to fuel giant radio relics on scales of ∼1–2 Mpc, with the required uniform population of fossil electrons, while they appear a sufficiently viable channel to fuel the regions of central radio halos.
• If merger shocks (here limited to those forming in $T\ge {10}^7\mathrm{K}$ gas, for computational limitations) also contribute seeds of cosmic rays, then the volume filling factor of fossil electrons in the central $r\le 1\mathrm{Mpc}$ region increases from 0.1–1% (in the case of CR electrons only injected by radio galaxies) to 3–30% if also shock-injected CR electrons are included. The injection by shocks occurs on much larger scales than an injection by radio galaxies. It can also naturally lead to large and correlated populations of fossil electrons that subsequent re-acceleration events can illuminate on ∼$\mathrm{Mpc}$ scales and with flatter spectra than the one from radio galaxies, owing to the shorter time elapsed between re-acceleration events.

A number of limitations prevent us from generalising our results, which are a first attempt at modelling diffuse radio emission in clusters, caused by the seeding by radio galaxies.

First, our prescription for radio galaxies is simplistic, in the sense that we only considered relatively low-power jests, activated all at the same time and only once in the late cluster evolution. Radio galaxies with higher powers can inject more fossil electrons, even if our tests in [27] and Vazza et al. [28] actually showed that too-powerful jets tend to spread electrons too widely. On the other hand, recent low-frequency surveys have established that the most massive radio galaxies in clusters and groups of galaxies are “on” most of their lifetime [70]. We also did not include self-regulated feedback and the effect of multiple AGN bursts, leading to repeated radio jet activity, as observed, e.g., [20]. While it is plausible that the collective effect of all the above missing ingredients can significantly increase the volume filling factor of fossil electrons from radio galaxies, it still remains to be proven that radio galaxies alone can fill with seed CR electrons the ∼1–2 Mpc wide radio relics, and result into the uniform spectral and brightness distribution observed in the most spectacular examples. Our results imply already that the single activity of just a few radio sources is unlikely to do that.

Second, discretisation effects prevent us from following the additional diffusion of electrons on small scales. However, this should not dramatically alter the mixing pattern captured by our tracer approach. Nevertheless, it will increase the filling factor that we can measure (hence a lower limit). Our tracer technique is suitable to follow most of the turbulent diffusion acting on cosmic rays, considering that the propagation of tracers is sensitive to the largest eddies (∼100–200 kpc) resolved in the turbulent cascade by the adopted MHD solver, e.g., [71,72]. However, our numerical method cannot capture turbulent diffusion on the smallest scales, which can happen on timescales of the same order as our cooling and re-acceleration timescales. The CR diffusion timescale can be roughly quantified in the hydro regime as $l_{\mathrm{D}}\sim 2\sqrt{D_{\mathrm{turb},\mathrm{small}}·t}$, where the small-scale turbulent diffusion coefficient is $D_{\mathrm{turb},\mathrm{small}}\sim (V_{\mathrm{turb},\mathrm{small}}·\mathsf{\Delta }x)/3$ (with $\mathsf{\Delta }x$ and $V_{\mathrm{turb},\mathrm{small}}$ being, respectively, the effective numerical resolution at the tracer location and the turbulent velocity fluctuation within that scale). Based on our previous work, in which the same cluster was simulated, we can estimate a typical rms turbulent velocity of $V_{\mathrm{turb},\mathrm{small}}\le 5\mathrm{km}/\mathrm{s}$ for $\mathsf{\Delta }x\le 8-16\mathrm{kpc}$, which gives $l_{\mathrm{D}}\le 0.01\mathrm{kpc}$ over our typical timestep $\delta t\approx 32\mathrm{Myr}$, or ∼$2\mathrm{kpc}$ for the entire simulated evolution. However, despite the fact that the effect of unresolved CR diffusion should be small and subdominant compared to that of turbulent diffusion, a conclusive view of the circulation of CR electrons in the simulated ICM requires a fully fluid treatment of CR, e.g., [15]. Moreover, a few other obvious injectors of fossil radio electrons are missing in our model, most noticeably the injection from star formation winds, e.g., [73,74,75] whose impact is observed at radio frequencies, e.g., [76].

Finally, we must remark that the applicability of our conclusions is limited to the fact that we only analysed one cluster of galaxies, which prevents us from monitoring how the transport and ageing of particles change as a function of the host cluster mass and, probably most importantly, on the host cluster dynamical state. The latter is important, because of the striking observed association between disturbed and merging clusters of galaxies and central radio halos or radio relics, and between more relaxed and cool-core clusters of galaxies with mini-halo emission, e.g., [4]. Our previous works, employing Lagrangian tracers, have shown that the turbulent mixing and dispersion of tracers tracking mass accretions of simulated clusters with different dynamical histories are nearly indistinguishable after a few ${10}^2\mathrm{Myr}$ [71], and so are the typical statistics of turbulence and shocks recorded along their propagation history [72]. However, we did find some significant differences in the mixing pattern of tracers initially located within the core of already formed relaxed or perturbed clusters of galaxies and meant to track the mixing of metals with the ICM. Our analysis indeed measured a smaller degree of mixing of tracers outside of the cores of the most relaxed clusters in the simulated sample, even after a few Gyr of evolution, e.g., [71]. However, no AGN feedback was modelled in those early works, and it is likely that this artificially reduced the amount of turbulent and mixing motions in the core of the most relaxed clusters. Only with future resimulations, will we be in the position of generalizing the trends suggested by the work presented here.

In summary, our results show that in order to properly model the long-term evolution of radio emission in clusters and groups it appears necessary that all key ingredients of cosmic ray enrichment are included. Their relative importance can vary with spatial scale, proximity to sources, distance from the cluster centre and redshift. This, in turn, poses new challenges to the numerical modelling, as well as to the increased computational costs to co-evolve the thermal and non-thermal components of the ICM—jointly with the pervasive magnetic fields which link the two.