© The Authors 2025

1. Introduction

Very sensitive, modern radio observations have allowed for the detection of extended and diffuse emission from the extreme periphery of clusters of galaxies (Cuciti et al. 2022; Botteon et al. 2022) and from filaments (Govoni et al. 2019; Botteon et al. 2020; Vernstrom et al. 2021, 2023; Venturi et al. 2022). These observations provide evidence for substantial magnetic fields (B ∼ 10 − 100 nG), as well as for relativistic electrons (E ≥ 0.1 − 1 GeV) existing on these scales. By and large, such radio-emitting electrons cannot be accelerated in situ or directly from the thermal pool; instead, they are expected to be re-accelerated by large-scale plasma perturbations, either in shocks (e.g. Kang et al. 2012; Pinzke et al. 2013; van Weeren et al. 2017) or turbulence (e.g. Jaffe 1977; Brunetti & Lazarian 2010; Brunetti & Jones 2014; Miniati 2015). The existence of pockets of fossil radio plasma, whose synchrotron emissivity gets boosted by compression following accretion phenomena in clusters, is also the leading model used to explain the formation of more classical ‘radio phoenices’ (Brüggen & Kaiser 2002; de Gasperin et al. 2015). It has been suggested in the context of other recent observations of faint and diffuse remnant radio emission from the repeated activity of central radio galaxies (e.g. Brienza et al. 2021; Biava et al. 2021).

There are several mechanisms that can naturally lead to the injection of large amounts of cosmic ray electrons (CRe, hereafter) in large-scale structures. In this work, we focus on the three that are believed to be the most significant. Firstly, outflows and jets from active galactic nuclei (AGNs) can store a large fraction of their internal energy in the form of non-thermal components (cosmic rays and magnetic fields), although the balance between thermal and non-thermal components may depend on the radio galaxy type, as well as on the surrounding environment (e.g. Völk & Atoyan 2000; Croston et al. 2008, 2018; Vazza & Botteon 2024).

Secondly, the ubiquitous process of star formation can also enrich galactic halos with cosmic rays, through the collective inflation of expanding shocked shells produced by main sequence stars, red supergiant stars, Wolf-Rayet stars, and supernova remnants (Dorfi 2004; Pfrommer et al. 2017; Seo et al. 2018; Modak et al. 2023). Thirdly, the growth of cosmic structures is accompanied by the formation of shock waves, which dissipate a large fraction of infalling kinetic energy into gas heating (e.g. Sunyaev & Zeldovich 1972; Bykov et al. 2008). Strong quasi-stationary accretion shocks (ℳ ≥ 10 − 10² where ℳ is the sonic Mach number) are expected beyond the virial radius of clusters of galaxies (e.g. Molnar et al. 2009; Vallés-Pérez et al. 2024), as well as shocks connected to the interface between filamentary accretions and the intra-cluster or intra-group medium (e.g. Brown 2011), or internal merger shocks violently crossing the innermost regions of clusters (e.g. Markevitch & Vikhlinin 2007). Basically, all modern cosmological simulations agree that the bulk of kinetic energy in the cosmic volume gets dissipated by 2 ≤ ℳ ≤ 4 shocks (see e.g. Vazza et al. 2011, for a comparison of cosmological simulations) and that the ensemble of shocks developed in the cosmic web can steadily refill cosmic structures with cosmic rays (e.g. Miniati et al. 2000; Ryu et al. 2003; Pfrommer et al. 2006; Planelles & Quilis 2013; Vazza et al. 2012; Schaal et al. 2016) and also produce levels of radio emission that is potentially detectable with the next generation of radio interferometers (e.g. Vazza et al. 2015; Weltman et al. 2020; Böss et al. 2023a). However, the ℳ ≤ 5 regime is a one where making predictions of CR acceleration gets difficult, owing to a number of physical and numerical uncertainties (e.g. Bykov et al. 2019; Böss et al. 2023b, for recent discussions).

Modern numerical simulations have attempted, with various degrees of success, to cover the daunting physical complexities and vast range of scales needed to model the interplay between multi-phase accretion onto SMBHs, and their ejected relativistic jets (e.g. Perucho & López-Miralles 2023; Bourne & Yang 2023, for a few recent reviews). Cosmological simulations have been successfully used to simulate the ‘radio’-mode feedback, namely, mediated by kinetic jets dissipating their energy via heating of the gas reservoir undergoing cooling (e.g. Puchwein et al. 2008; Fabjan et al. 2010; Dubois et al. 2012; Tremmel et al. 2017; Weinberger et al. 2017; Dolag et al. 2017; Vogelsberger et al. 2020; Bourne & Yang 2023). However, to the best of our knowledge, no cosmological simulation has so far included an even simplistic model to co-evolve magnetic fields and cosmic ray electrons injected by radio galaxies, star formation and shocks at the same time.

Here, we introduce a new suite of cosmological simulations with cosmic-ray electrons, injected at run time by shock waves, star feedback, and AGNs, which are also sources of magnetic fields and allow us to carry out a run-time monitoring of how these non-thermal energy components are advected in the cosmic volume. Our paper is structured as follows: in the following section (Sect. 2), we give an overview of the new numerical implementations adopted in our runs. In Sect. 3, we give our main results, both on the simulation of galaxy properties and on the properties of the cosmic ray electron fluid and magnetic fields on various scales of the cosmic web. A discussion of the main limitations of our approach is given in Sect. 4, while our main conclusions are summarised in Sect. 5.

2. Methods

2.1. ENZO cosmological simulations

We produced new ideal magneto-hydrodynamical (MHD) cosmological simulations with the Eulerian code ENZO v2.6¹ (Bryan et al. 2014). These simulations include a number of ad-hoc modifications to inject a passive cosmic-ray electron fluid and to implement source feedback mechanisms for reproducing the population of active radio galaxies.

This version of ENZO allows us to take full advantage of the GPU architecture on the LEONARDO supercluster at CINECA, which ensures that our cosmological simulations are four times faster overall than their corresponding CPU-only version when working on the problem being analysed. For this project, we used the GPU-accelerated version of ENZO (Wang et al. 2010), which allowed us to run our simulations employing 4 GPU per node and 1024 MPI tasks in total. Larger simulations, using 2048³ cells and 1024³, are presently running and require 32 nodes instead, employing 4 GPU per node and 1024 MPI tasks in total. For a review of the basic properties of the code, we refer to the original paper on its method (Bryan et al. 2014), while for more specific details about the MHD version used here, we refer to our previous works (Vazza et al. 2017, 2021a).

We used the following cosmological parameters: h = 0.678, Ω_Λ = 0.692, Ω_M = 0.308 and Ω_b = 0.0478, based on the results from the study by Planck Collaboration XIII (2016). As a compromise between the requirements for a reasonable resolution of halo mass growth, feedback between galaxies and their surrounding medium and shock dynamics within halos, as well as their interface with filaments and cosmic voids, we resorted to a uniform resolution in our simulations. We used 1024³ cells and dark matter particles to cover a 42³ Mpc³ volume (additionally, to increase the statistics for an object, we complemented our suite also with a 2048³ simulation of a 85³ Mpc³ volume). The resulting spatial resolution (41.5 kpc comoving) is coarser than state-of-the-art cosmological simulations of galaxy evolution, in which typically adaptive-mesh refinement, smoothed particle hydrodynamics or moving-mesh techniques are employed (e.g. Vogelsberger et al. 2020). However, we calibrated sub-grid models of star formation and feedback to reproduce the most important global statistics of the stellar mass assembly in galaxies, as well as the energy input of AGNs from supermassive black holes (SMBHs) during cosmic evolution.

This approach has allowed us to study, for the first time, how reservoirs of cosmic-ray electrons build up in the cosmic web as a function of time. Moreover, we have been able to study how AGNs produce different enrichment patterns of fossil electrons depending on the AGN feedback mechanism, while broadly reproducing the bulk properties of radio galaxies across redshift and environments (e.g. cosmic star formation history, stellar distribution function of galaxies, etc.).

It is impossible to fully explore the large number of sub-grid parameters used to parametrise star formation and feedback, AGN feedback, and the injection of magnetic energy and CRe associated to these processes, which are non-linearly coupled to galaxy evolution. Here, we focus on the comparison of nine resimulations of the same cosmic volume, featuring plausible variations in sub-grid prescriptions for star formation or AGN feedback, and we compare them to key observable properties of galaxies (e.g. Sects. 3.1–3.3). As we show in the reminder of the analysis, there is one model that best reproduces all investigated properties (B4 in the following), while the others end up failing to reproduce some key properties of the galaxy population. Nevertheless, even failed models are useful for understanding the role of specific mechanisms in the seeding of CRe and magnetic fields on the scale of the cosmic web.

In detail, we tested: (a) three model variations for AGN feedback, where we kept the star formation recipe fixed (runs A1, A2, A3); (b) four model variations, where we varied the star formation recipe, and kept the AGN modelling fixed (runs B1, B2, B3, B4); and (c) two model variations of one of the previous models, where we added primordial magnetic fields (run C1) and also tested the obliquity-dependent injection of CRe by shocks (run C2). In all runs, for gas cooling, we assumed for simplicity a gas equilibrium cooling model, with a primordial gas composition and with no metal enrichment modelled or followed at run time. Table 1 gives the most important parameters describing the differences in our suite of simulations, while in the next sections, we introduce in detail the adopted numerical prescriptions for magnetic fields, shock acceleration, AGNs, and star formation, as well as for the advection and ageing of our CRe.

Figures 1 and 2 show simulated magnetic field and temperature distributions at three different epochs (z ≈ 2, ≈1, and ≈0), for our B4 model. The distribution of thermal and non-thermal energy components in the simulated cosmic web is, at all epochs, a combination of the large scale dynamics of the cosmic web (halos, filaments and voids components) and of intermittent and an ‘inside-out’ impulsive release of feedback energy related to the evolution of simulated galaxies. In the following, we explain in more detail how the different physical implementations affects the evolution of the quantities above (and of their associated CRe populations). Finally, we discuss our results and how variations in these procedure may affect the measured global trends.

Fig. 1. Projected mean magnetic field intensity (top panels, in units of comoving Gauss) and projected mean mass-weighted gas temperature (bottom panels) across the full simulated 42 Mpc3 volume in our B4 run at z ≈ 2 (left), for the z ≈ 1 (centre) and z ≈ 0 (right) epochs. Movies of some of our runs can be found online and at https://www.youtube.com/playlist?list=PL8ecsjnxOKP7LXQICrNBLHZkPgxh0jCzj

Fig. 2. Same as Fig. 1, but for a zoomed portion with side 15 × 15 Mpc2 (and 42 Mpc along the LoS) for our B4 run at three different epochs.

2.2. Primordial magnetic fields and sub-grid dynamo amplification

In most cases (runs A1, A2, A3, B1, B2, B3, and B4), the magnetic field in our simulated volume was initialised with a simple uniform volume-filling magnetic field with comoving amplitude of B₀ = 10⁻¹⁶ G at the beginning of the simulation (z = 40); namely, at the level of the lower limits from the non-detection of the inverse Compton emission from blazars (Neronov & Vovk 2010; Alves Batista & Saveliev 2021; Tjemsland et al. 2024). We initialised the magnetic fields with low amplitudes to make it easier to detect the implemented magnetic feedback from astrophysical sources in the simulation, since ‘magnetisation bubbles’ are expected to have fields with strengths of ≥10⁻¹² G (Arámburo-García et al. 2021; Bondarenko et al. 2022). In the C1 and C2 simulations, we used a tangled primordial magnetic field instead, with the fields’ scale dependence described by a power law spectrum: P_B(k) = P_B0k^n_B (n_B = −1.0), whose root-mean-square (rms) amplitudes after smoothing the fields falls within a scale λ = 1 Mpc (e.g. Paoletti & Finelli 2019) is set to B_1 Mpc = 0.37 nG. We notice that this value is about five times lower than the upper limits inferred by Paoletti & Finelli (2019) using priors from the analysis of the cosmic microwave background. This renormalisation was motivated by the latest results on extragalactic Faraday rotation using LOFAR (Carretti et al. 2023, 2024). Since our goal in this paper is to couple different CRe injection mechanisms with a plausible primordial magnetic field configuration, we focus here on the n_B = −1.0 spectrum, instead of using much flatter (i.e. the scale-invariant) or steeper (i.e. low-scale spectra from causal mechanisms) power spectra, which appear to be disfavoured by the comparisons of simulated and observed Faraday rotation data (e.g. Vazza et al. 2021b; Carretti et al. 2024; Mtchedlidze et al. 2024).

These simulations have an spatial resolution that is insufficient for resolving flows in halos with a high Reynolds number and can hardly capture any small-scale amplification within halos (e.g. Donnert et al. 2018, and references therein). However, having a realistic level of magnetisation in our halos is necessary to produce realistic estimates of synchrotron radio emission from our simulated cosmic web. Following previous works (e.g. Vazza et al. 2017, 2021b), we used a simplistic model to incorporate the expected effect of small-scale dynamo amplification, based on the run-time analysis of the solenoidal velocity field, which is the one mostly responsible for the amplification in halos (e.g. Ryu et al. 2008). In particular, the vorticity squared, (∇×v)² ≡ ϵ_ω (where v is the local gas velocity), is easy to measure at any timestep and it offers a convenient way to estimate the turbulent dissipation rate of solenoidal turbulence on the fly. This is expressed as F_t ≃ η_tρϵ_ω³/L, where ρ is the gas density in the cell, η_t = 0.014 in Kolomogorov turbulence and L is the stencil to compute the vorticity, as in Jones et al. (2011). We assumed dynamo amplification operates below the resolved scales in the simulation; thus, we generated new magnetic energy by converting, during a run time, a fraction of the kinetic energy power into the magnetic field energy: E_B, dyn = ϵ_dyn(ℳ_t)F_tΔt. The conversion factor must be calibrated based on dedicated simulations of dynamo amplification, and as tested in previous work (Vazza et al. 2017, 2021a,b), we used the analytical prescriptions derived by Federrath et al. (2014) for ϵ_dyn, which are given on the basis of the local turbulent Mach number.

Once the additional magnetic energy is computed, we generated a three-dimensional magnetic field, B_turb, to add the existing magnetic field. We simply imposed the condition that B_turb must be parallel to gas vorticity, so that the new generated field is also (with good approximation) solenodial by construction. The equivalent amount of kinetic energy is removed from the same cell and momentum is removed assuming an isotropic dissipation of the small-scale velocity vectors. This procedure is manifestly simpler than more sophisticated subgrid models available, which measure the sub-grid turbulent energy via Favre filtering, and incorporate the electromotive force in a self-consistent way (Grete et al. 2016). Although our method has an advantage of producing a realistic level of magnetic fields in the simulated halos. We did not allow this mechanism to operate at cosmic density levels lower than that of filaments, that is, with respect to gas densities lower than ten times the cosmic mean gas density. This was done on the basis that no dynamo amplification has ever been found in this environment even with high-resolution MHD simulations; as described in, for example, Vazza et al. (2014).

The second row of Fig. 3 shows the distribution of magnetic fields at z = 0.02 for all our runs. We see that which illustrates well the much more efficient filling of cosmic volume by primordial magnetic fields (runs C1 and C2), as well as slightly different scales and intensities of magnetisation by such fields than each different run (e.g. where we used combined feedback from stars and AGNs) could produce.

Fig. 3. Top: Projected mean (mass weighted) gas temperature (in units of [K]) across the entire simulated volume for our seven runs at z = 0.02. Bottom: Projected mean (mass weighted) magnetic field strength (in units of [G]) for the same epoch and volume selection. The selected volume in all panels is about 21.25 × 42.5 Mpc2 and the length along the LoS is 42.5 Mpc.

2.3. The spatial advection and the age determination of relativistic electrons

Unlike CR protons, CR electrons can be effectively treated like a passive fluid (owing to their small pressure compared to the thermal gas and kinetic pressure, as well as to the magnetic pressure everywhere in the cosmic structures) which is frozen to the gas fluid because of the extremely small gyro-radius (r_L, e ∼ 10⁻⁶ pc for a μG magnetic field and γ ∼ 10⁶ electrons, where γ is the Lorentz factor), which is about the highest energy which electron can reach before radiating their energy via synchrotron and inverse Compton, in less than 1 Myr, (e.g. Sarazin 1999). With our modification of the ENZO code, we separately track the advection of the three kinds of CRe: the ones injected by shocks, the ones injected by AGNs, and the ones injected by star formation (see sections below for details on the injection procedures). The field that is advected for each family is the referred to as the number density (similarly to what ENZO already does for any other chemical species). The choice of the number density as the advected field is crucial here; this quantity is conserved even in the presence of re-acceleration processes (not included in our model). The thermalisation of CRe can proceed in a short time only at very high particle densities (n ≫ 0.1/cm³) via ionisation or Coulomb losses. This means that these loss processes can affect the simulated budget of CRe in the interstellar medium within galaxies, while instead our predictions concerning the radio emission by shocks or by radio jets are unaffected by the thermalisation of CRe because they always propagate through a much lower-density medium.

Our model lacks a spectral ageing model for CRe (numerically very expensive), yet we did want to recover the typical evolutionary timescales of our injected CRe everywhere in the simulation. Thus, we chose to explore, for the first time, the application of a simple, yet, powerful approach applied by Beckmann et al. (2019) in AGN feedback simulations of single objects. Specifically, with our ENZO modification we evolve, together with each of the three families of CRe, along with a second ‘mirror’ fluid, which is injected and advected exactly in the same way, but additionally subject to an artificial exponential decay with a fixed, arbitrary timescale, τ, so that $n{\prime }_{\mathrm{CRe}}\propto e^{-t/\tau }$. Numerically, this means that all CRe decaying species are advanced in time as $n{\prime }_{\mathrm{CRe}}(t+\mathrm{\Delta }t)=n{\prime }_{\mathrm{CRe}}(t)(1-\mathrm{\Delta }t/\tau )$. This allows us, at any stage in the simulation, to retrieve the time elapsed since the last injection of CRe in all cells (separately for each CRe species) as

$$\begin{array}{c}\hfill t_{\mathrm{age}}=-\tau log\left(\frac{n_{\mathrm{CRe}}^{\prime }(t)}{n_{\mathrm{CRe}}(t)}\right),\end{array}$$(1)

where n_CRe is the number density for each of the primary CRs (i.e. not subjected to the artificial decay law). To the best of our knowledge, this approach has never been applied to the simulation of CRe dynamics; it has the advantage of providing an efficient way to estimate of the typical age of CRs in 100% of the simulated volume. Provided that the expected evolution of CRs is simple, this information can be used to produce realistic predictions for the full energy spectrum of particles tracked by the total CRe density field (see the section below).

It should be noted that as long as τ is a known number, it can have any value to allow us to derive the effective age of the CRe in the cells, at least for the case of a single injection epoch for CRe. However, when multiple injections are involved, as in our simulations, it is important to use a physically motivated value of τ, so that the weighting of older and younger populations is realistic as well. Therefore, in our models we fixed τ = 0.1 Gyr, which is the representative value of the cooling time of relativistic electrons in the γ ∼ 10⁴ range, which dominate synchrotron radio emission for magnetic field intensities in the B ∼ 0.1 − 10 μG.

In detail, the cooling time is given by (e.g. Sarazin 1999; Brunetti & Jones 2014):

$$\begin{array}{c}\hfill {\tau }_{\mathrm{cool}}=\frac{0.77\mathrm{Gyr}}{(\gamma /300)[{\left(\frac{B}{3.25\mathrm{\mu }\mathrm{G}}\right)}^2+{(1+z)}^4]}·\end{array}$$(2)

The range of z and γ for which the cooling time (considering synchrotron as well as inverse Compton losses) is of the order of τ is shown in Fig. 4, for B = 0.1 μG and B = 10 μG. Our choice of τ = 0.1 Gyr thus ensures that, with good approximation, whenever we compute the radio emission from our distribution of CRe (Sect. 2.7), we are correctly weighting the contribution from the CRe in the cell that will actually dominate the radio emission in the real case – especially for long evolutionary timescales (≥0.1 Gyr). On the other hand, in the case of multiple AGN bursts on very short timescales and in the presence of significant losses in dense environments, the presence of multiple complex component may make our approach too simplistic. However, this limitation is not too severe for the large-scale global analysis performed in this work. Additional tests in Appendix A.1 show that the exact choice of this parameter (at least in the 0.01 − 1 Gyr range) is not crucial in the age determination of all CRe species analysed in this work.

Fig. 4. Cooling time for relativistic electrons underoging synchrotron and inverse compton losses (as in Eq. (2)), as function of their Lorentz factor and cosmic redshift, for two B = 0.1 μG and B = 10 μG. The solid black lines mark the region corresponding to our choice of the timescale τ for the age determination of CRe injected at run-time.

The different visual patterns of the elapsed time since their last injection for the CRe associated to different mechanisms is shown in Figures 5–7, for our runs at z = 0.02. Here, we produce maps of the age of CRe within each cell, projected along the LoS by weighting for the CRe density. The quantitative information is partially lost due to the smearing of the age distribution along the LoS, yet the maps give a qualitative view of the differences in ages related to the various mechanisms.

Fig. 5. Projected mean (CRe-weighted, only considering CRe from shocks) age of CRe since their last injection by shocks (in units of [Gyr]) across the entire simulated volume for our seven runs at z = 0.02.

Fig. 6. Projected mean (CRe-weighted, only considering CRe from AGNs) age of CRe since their last injection by AGNs (in units of [Gyr]) across the entire simulated volume for our seven runs at z = 0.02.

Fig. 7. Projected mean (CRe-weighted, only considering CRe from star formation) age of CRe since their last injection by star formation (in units of [Gyr]) across the entire simulated volume for our seven runs at z = 0.02.

Compared to the other mechanisms, shocks inject (on average) the ‘freshest’ electrons at all epochs, given the ubiquitous process of shock formation. This yields to a relatively young (≤0.1 − 0.4 Gyr in the maps) population of CRe at the periphery of halos and filaments and occasionally associated to the launching of shocks within halos (Fig. 5). Most of the projected volume has been filled, at some time in the past, with the injection of small quantity of CRe, since the process of shock formation begins early and can sweep a large fraction of the cosmic web. It should be noticed that the distribution of shock age and morphologies is not the same in all runs, due to galaxy feedback processes. Indeed, beside directly injecting new CRe in jets or winds (see the next sections), feedback also triggers new shocks, which in turn inject additional CRe as they expand into the intergalactic medium, an effect that has been already identified by previous cosmological simulations (e.g. Kang et al. 2007; Schaal et al. 2016).

Figure 6 illustrates the typical bubble structure of the regions enriched with CRe by AGN feedback instead. In all our models, AGN feedback becomes significant only for z ≤ 2, which explains why there are projected regions with no enrichment of CRe at all, even in this very overdense selection of the simulation. As it will be discussed later in Sect. 3.5, the three-dimensional filling factor of AGNs and stellar feedback bubbles is typically ≤10 − 35% in our runs. While the first expanding bubbles of CRe have long injection timescales (≥2 − 3 Gyr in the maps), the recurrent repeated activity of AGNs is marked by the smaller and darker blobs at the centre of bubbles, which are associated with very recent (≤0.1 Gyr in the maps) feedback episodes.

In the case of the age distribution of CRe injected by stellar feedback (Fig. 7) the bubbles are confined to a smaller volume, and they visually show a larger correlation with the underlying cosmic web distribution. Their average age is also much larger compared to the other two mechanisms, following from the fact that the peak of star formation was at z ∼ 2 in all runs. More recent episodes of injection of CRe by stellar feedback tend to be hard to be seen in projection, as they are covered by the information of older populations along the LoS.

This approach to monitor the injection and circulation of relativistic electrons has a strong predictive power, despite its great simplicity, due to aforementioned fact that the number of CRe does not evolve even in the presence of processes which are not included (yet) in our model, like Fermi reacceleration or loss processes of CRe. Of course, the actual radio emission from CRe crucially depends on their spectral energy distribution, and in this case we can only approximately guess it, as it will be discussed in Sect. 2.7.

This fluid approach to study CRe is also complementary to previous work by our group (Wittor et al. 2017; Vazza et al. 2021c, 2023; Smolinski et al. 2023), in which we employed Lagrangian particles to track relativistic electrons, injected in post-processing in our snapshots. Here, the application of a Eulerian fluid approach allows us to monitor the evolution of the fluid in 100% of the simulated domain, unlike what can be done through the application of a finite number of tracers – for which only the central regions of single galaxy clusters could be well sampled. Moreover, this approach allows us to evolve electrons at run-time and use the parallelization as well as the GPU acceleration available in the ENZO code, and apply this approach to very large simulations. This means that we can study the distribution of electrons in a O(1) fraction of the ∼1 − 8 ⋅ 10⁹ cells of our simulations at all timesteps, namely, four to five orders of magnitude more data that we can do with Lagrangian tracers.

2.4. The injection of relativistic electrons by shocks

To inject at run-time a fluid of fresh CRe accelerated by shocks based on Diffusive Shock Acceleration (DSA), we use an on-the-fly shock finder, which we coded in ENZO and used already in previous works (Vazza et al. 2012), which modifies the original algorithm by Ryu et al. (2003). In essence: (a) we select candidate shocked cells based on their negative velocity divergence; (b) we perform one-dimensional scans along the three spatial axes and check whether the local changes in gas temperature and entropy are equiverse (e.g. Ryu et al. 2003; c) we use the local gradient of temperature sets the shock propagation verse, and (d) we compute the local Mach number from the pressure jump using ideal Rankin- Hugoniot jump conditions.

For CR electrons to be injected via DSA, their Larmor radius must be of the same order of that of CR protons, implying a dependence on the local plasma conditions (i.e. downstream gas temperature, magnetic field strength and topology). Determining the exact injection momentum there (and, if there is a clean cut injection momentum at all) is still a theoretical challenge, which particle in cell (PIC) simulations have only since recently started to tackled with dedicated simulations (e.g. Amano & Hoshino 2007; Arbutina & Zeković 2021; Grassi et al. 2023; Gupta et al. 2024). A general finding is that, while the injection of CR protons from the thermal distribution is favoured in quasi-parallel shocks (e.g. Ha et al. 2018; Ryu et al. 2019), the situation for CR electrons is more complex, as their injection can proceed both in quasi-parallel or quasi-perpendicular shock geometries, through a combination of different acceleration mechanisms (e.g. Xu et al. 2020; Shalaby et al. 2022; Gupta et al. 2024). However, the lack of fully self-consistent three-dimensional PIC simulations capable of scanning the entire range of parameters relevant for collisionless non-relativistic structure formation shocks makes it difficult to generalise from the above work.

Given the present uncertainties, here we rely on a simpler semi-analytical recipe, in which the injection of CR electrons is inferred based on the one of CR protons (e.g. Kang 2024). In this view both, protons and electrons need to be pre-accelerated to suprathermal momenta greater than the so-called injection momentum, p_inj = x_injp_th, p, in order to diffuse across the shock transition layer and fully participate in the diffusive shock acceleration process, where $p_{\mathrm{th},\mathrm{p}}=\sqrt{2m_{\mathrm{p}}{k}_{\mathrm{B}}{T}_{\mathrm{p}}}$ is the thermal momentum of protons in the shock downstream and x_inj the injection parameter, which is not yet constrained by the DSA theory. In a shock, a fraction of the post-shock thermal particle density (n_p) becomes a number density of accelerated CR electrons (n_CRe), as (e.g. Kang 2018, 2020, 2021, 2024),

$$\begin{array}{c}\hfill n_{\mathrm{CRe}}={\xi }_e(\mathcal{M})n_{\mathrm{p}}=n_{\mathrm{p}}\frac{4x_{\mathrm{inj}}^3}{{\pi }^{3/2}({\alpha }_{\mathrm{inj}}-3)}exp(-x_{\mathrm{inj}}^2).\end{array}$$(3)

In all our simulations we assumed x_inj = 3.5, which is also in line with the latest indications from one-dimensional kinetic PIC simulations (e.g. Gupta et al. 2024). This means that for typical ∼10⁶ − 10⁸ K post-shock gas temperatures, CR electrons are injected with a Lorentz factor γ ∼ 3 − 30. Therefore it is worth stressing that, in the reminder of the paper we shall loosely use the term “Cosmic Ray electrons” (CRe) to refer to all electrons accelerated by DSA in this way, from electrons with γ ∼ a few, up to the super-relativistic regime of electrons emitting synchrotron radiation in the radio band (γ ∼ 10⁴), or beyond. The injection spectral index is computed from the DSA theory: α_inj = 4ℳ²/(ℳ² − 1), which is appropriate for plane shocks in the test particle regime (e.g. Kang 2011, and references therein). The trend of the ξ_e(ℳ) used in this paper, and the comparison with small different variations of x_inj are shown in Fig. 8.

Fig. 8. Trend of injection efficiency of CRe (given in terms of the number of injected CRe with respect to the number of thermal particles) as a function of shock Mach number, for three different choices of the injection momentum, xinj (the baseline model assumed in our simulations is in red).

In line with Kang (2024) we also considered a low Mach number limiter for the injection, ℳ_thr = 2.3 because, since several PIC simulations haver suggested a suppression of CR proton acceleration for low Mach number shocks in high β plasmas below ℳ ∼ 2.25 (e.g. Ha et al. 2018; Ryu et al. 2019). In choosing this low Mach number limiter, we are also motivated by the fact that very weak shocks are typically difficult to measure in cosmological simulations, owing to the amount of transonic and supersonic velocity fluctuations typically found in the same environment (e.g. Ryu et al. 2003; Vazza et al. 2011; Schaal et al. 2016).

The shock obliquity (i.e. the angle between the upstream magnetic field and the shock velocity vector) is a key parameter ruling the acceleration of cosmic rays. Several recent works, mostly using particle-in-cell simulations, have explored the relations between shock obliquity, Mach number and the efficiency in the injection of both relativistic protons (e.g. Caprioli & Spitkovsky 2014; Park et al. 2015; Ha et al. 2018; Ryu et al. 2019; Bohdan et al. 2020) and electrons (Guo et al. 2014a,b; Matsumoto et al. 2017; Bohdan et al. 2017, 2019; Xu et al. 2020; Arbutina & Zeković 2021). In the case of electron, a robust finding is that for quasi-perpendicular shocks (θ_B ≥ 45°) and strong enough Mach numbers, electron are pre-accelerated by the motional electric field by drifting along the shock front, via ‘shock drift acceleration’, and a significant fraction of these electrons eventually takes part to diffusive shock acceleration, and it becomes highly relativistic.

Therefore, in one of our runs (C2) we tested the role of shock obliquity by allowing the injection of the CRe fluid only for quasi-perpendicular shocks (θ_B ≥ 45°), by measuring at run-time the angle between the pre-shock magnetic field vector and the shock velocity vector, inferred at run-time by our shock finder. The result of this model are compared with the more standard implementation, in which the role of obliquity is not considered and all shocks equally contribute solely as a function of their Mach number (C1). Both the C1 and C2 runs were run assuming the same primordial stochastic model with P_B(k) = A_Bk⁻¹, in order to work with a more realistic distribution of angles between shocks and magnetic fields, compared to the more simple setup with a uniform primordial magnetic field used in all other runs.

The first column in Fig. 9 gives the evolving spatial distribution of CRe injected by shocks, for the example of our B4 model, at three different epochs. As already discussed above, shocks inject overall CRe in a very large fraction of the cosmic volume. DSA converts a nearly constant fraction of the gas matter into CRe (∼10⁻⁴ − 10⁻³, see Fig. 8), and it starts almost together with the structure formation processes, hence the large-scale distribution of CRe looks very similar to the standard gas cosmic web, although on small-scales some differences appear, due to the injection of recent and powerful shocks, triggered either by mergers or by AGNs.

Fig. 9. Projected number of CRe across the full simulated 42.5 Mpc3 volume in our B4 run (considering only 5 Mpc along the LoS to better isolate single structures), for the z ≈ 2, z ≈ 1 and z ≈ 0 epochs (from top to bottom). The first column shows CRe injected by shocks, the second by AGN feedback and the third by star formation. The full movie of these fields can be found online and at https://youtu.be/xMTSK_rxZZU?si=hj4Yckl3hBuhYF-v

2.5. The injection of relativistic electrons and magnetic fields by AGNS

In most cosmological simulations, SMBHs are modelled with Lagrangian ‘sink’ particles, which accrete mass from the surrounding cells and need to be injected in the simulation at arbitrary times and positions (e.g. Kim et al. 2011, for the case of ENZO). The feedback from SMBHs is added either by introducing isotropic thermal heating, or kinetic energy through directional jets along arbitrary directions (e.g. Li et al. 2015).

This approach comes with a number of numerical challenges: from the need of stopping and restarting simulations to inject SMBH seeds (which kills performances in large simulations), to the problems connected in SMBHs drifting away from the shallow potential of their host galaxy if simulated with insufficient resolution (see e.g. Damiano et al. 2024), to the ad-hoc prescriptions needed to produce realistic SMBH masses (e.g. Habouzit et al. 2022). Considering that our main purpose here is to produce realistic forecasts for the distribution of CRe injected by AGNs, without necessarily resolving the puzzle of how SMBHs co-evolve with their host environment, we explored a simplified approach which bypasses all aforementioned challenges (see Fig. 10).

Fig. 10. Example of six different CR electrons spectra computed by our post-processing model, for different (random) variations of the combinations between gas density, magnetic field strength, redshift, and time elapsed since the injection of the initial power-law distribution of CRs. Such spectra are extracted from our tabulated list of spectra generated with the ROGER code and applied to compute the synchrotron emission from every cell in our simulated volumes in post-processing. Finally, we note that our approach neglects the spatial diffusion of cosmic rays, which is mostly negligible on the scales of interest here (this is addressed in more detail in Sect. 4).

We avoided the use of sink particles entirely, assuming instead that each gas density peak in the simulation harbours a SMBH at every timestep, with a fiducial mass assigned based on observational scaling relations. In detail, we measure at run time the total baryonic mass around density peaks in the simulation (only considering cells at least with a total matter density 100 larger than the cosmic mean density at the same epoch), within a fixed comoving radius of ≤83 kpc (equal to 2 cells), and we use a scaling relation to compute the average SMBH mass that would be hosted by a real galaxy at that location. Among the many possible best-fit relations derived by Gaspari et al. (2019) for several observable proxies of SMBHs, we use the relation between the enclosed gas mass within the galactic potential, as in their Fig. 6; that is M_BH ≈ 2 ⋅ 10⁷ M_⊙(M_g/(10⁶ M_⊙)^0.57, which is easy to measure at run-time. To avoid overcrowding of SMBHs in the same galaxies, we implemented a local exclusion criterion, which prevents us from considering more than one single SMBH cell within the same cube of 7³ cells (i.e. closer than a distance of 3 cells, or 124.5 kpc, from the central gas density peak). This choice somewhat limits our capability of properly computing the SMBH masses of halos undergoing close mergers, since only one of the two gas density peaks (in case they are separated by less than three cells) is assigned with a SMBH mass. In Appendix A.2 we present a model variation of this scheme, showing that varying the minimum allowed separation of gas density peaks does not have a significant impact on our simulations.

Of course, given the coarse resolution of our simulation, a very detailed measurement of the gas mass may be prone to numerical uncertainties, and other quantities can be used as well to guess the SMBH mass (i.e. total stellar mass, central gas temperature, central gas density etc.). We defer a more detailed study of the potential use of different scaling relations to future work, where a higher spatial resolution will be employed. We impose M_BH ≥ 10⁶ M_⊙ as an additional requirement to trigger SMBH feedback (because such low mass halos are dominated anyway by stellar feedback).

For every identified SMBH location, during run-time we compute the instantaneous mass accretion rate with the standard Bondi-Hoyle formalism: ${\dot{M}}_{\mathrm{BH}}=4\pi {\alpha }_{\mathrm{B}}{G}^2{M}_{\mathrm{BH}}^2\rho /{(v_{\mathrm{rel}}^2+c_{\mathrm{s}}^2)}^{3/2}$, where c_s is the sound speed of the gas at the SMBH’s location, v_rel is the relative velocity between the accreted gas and the SMBH (which we assume = 0 here), ρ is the local gas density and α_B is an ad-hoc parameter meant to compensate for the lack of resolution around the Bondi radius (e.g. Booth & Schaye 2009; Gaspari et al. 2012; Tremblay et al. 2016). The bolometric luminosity of the black hole is defined as L_BH = ϵ_rṀ_BHc², where ϵ_r is the radiative efficiency of the SMBH, and P_j = ϵ_BH L_BH = ϵ_BHϵ_rṀ_BHc² is the feedback power, in which ϵ_BH is the factor that converts the bolometric luminosity of the SMBH into the thermal feedback energy (we use ϵ_BH = 0.05 in line with the literature), and ϵ_r = 0.1 is the radiative efficiency of the SMBH.

We followed the approach of other cosmological simulations and adopted two distinct feedback modalities depending on the temperature of the accreted matter (e.g. Steinborn et al. 2018):

• “cold gas accretion” when T ≤ 5 ⋅ 10⁵ K, where we use a low value for α_B and distribute most of the feedback energy in the form of dipolar dumps of thermal energy, with velocity vectors pointing at the opposite sides of the SMBH; roughly mimicking “quasar feedback” in halos accreting cold gas via mergers.

• “hot gas accretion” when T > 5 ⋅ 10⁵ K, where we use a large value orf α_B and assume that most of the feedback energy is in the form of kinetic jets; thus mimicking ‘radio feedback’ in halos mostly accreting hot gas.

The actual values of α_B, cold and α_B, hot used in all runs is given in Table 1. To estimate the energy ratios going to each energy form, we first consider in all cases that ϵ_B, AGN = 10% is the fraction of energy released in magnetic energy, through pairs of magnetised loops wrapped around the direction of kinetic jets. Every time jets are activated, their direction is randomly drawn along any coordinate axis of the simulation. In all cases we injected new CRe within the jets, with a number density equal to a fixed fraction of the number density of the thermal gas in jets, n_CRe, AGN = ξ_AGNn_jet, with ξ_AGN = 10⁻³ in our runs. This specific value of ξ_AGN has been calibrated a-posteriori, based on the comparison between the simulated radio luminosity function of our AGNs with the observed one (Sect. 3.3), yet it is close to typical choices in the literature (e.g. Mendygral et al. 2012; Nolting et al. 2019).

For the thermal and kinetic energies, we assumed that (a) in “cold gas accretion” feedback, the remaining 72% of energy is thermal and 18% is kinetic, while in (b) the “hot gas accretion” feedback, 81% of the remaining energy is kinetic and the rest is thermal. Similar numbers are quoted in the recent literature (e.g. Steinborn et al. 2018; Davé et al. 2019) and even though they are susceptible to changes when working at a different resolution, they still represent parameters giving the best results in our runs (see Sect. 3).

The second column in Fig. 9 gives the evolving spatial distribution of CRe injected by AGNs, for the example of our B4 model, at three different epochs. The distribution of CRe injected by AGNs is, as expected, much less volume filling than the one of CRe injected by shocks, but it reaches large values in the proximity of very active AGNs.

2.6. The injection of relativistic electrons and magnetic fields by stellar feedback

To couple the injection of CRe to star formation, we used the Kravtsov (2003) model, implemented in ENZO. Its key parameters are the minimum gas density threshold to form stars, n_*, the dynamical timescale of the star formation process, t_*, and the minimum mass, m_*, for newly formed stars. Whenever n ≥ n_*, the gas is contracting (∇ ⋅ v < 0), the local cooling time is smaller than the dynamical timescale (t_cool ≤ t_*) and the baryonic mass in the cell is larger than m_*, a new “star particle” (actually representing an entire stellar mass distribution) is formed with a mass m_* = m_bΔt/t_*. This recipe has been proposed to reproduce the observed Kennicutt’s law (Kennicutt 1998) and in our simulations can well match the observed cosmic star formation history (e.g. Vazza et al. 2017). Table 1 gives the parameters that overall best reproduce the cosmic star formation history in our runs, and are meant just to give an effective sub-grid model prescription of the star formation process, averaged over the large 41.5³ kpc³ volume resolution of our runs.

The feedback from star formation depends on the assumed fractions of energy, momentum, and mass ejected per each formed star particles, E_SN = ϵ_SFm_*c², with ϵ_SF = 10⁻⁸ − 5 ⋅ 10⁻⁸ as fiducial parameter, as in previous work (Vazza et al. 2017). 90% of the feedback energy is released in the thermal form (i.e. hot supernovae-driven winds), distributed among the 27 nearest cells around the star particle, while 10% in the form of magnetic energy, assigned to dipoles during each feedback episode.

CRe are injected into the same cells, using a fixed fraction of the local gas density, ξ_SF = n_CRe/n_g. This fraction depends on several processes: (a) the direct injection of cosmic rays by shocks driven supernova remnants and pulsar wind nebulae and (b) the continuous injection of secondary CRe from the hadronic collisions between thermal protons and CR protons; (c) the thermalisation of CRe due to collisional and ionisation losses in the dense interstellar medium (e.g. Pfrommer et al. 2017). With preliminary tests with smaller volumes at the same spatial resolution (e.g. as in the appendix), we calibrate a reasonable value for such “macroscopic” injection efficiency to ξ_SF = 10⁻⁵ as fiducial value, as this ensures that the synchrotron emission from star-forming galaxies in the high luminosity end of our distribution is in line with observations (e.g. Cochrane et al. 2023; Heesen et al. 2023). However, as we shall discuss in Sect. 3.3, our simplistic treatment of CRe from star formation cannot produce a match of the entire observed power distribution of star-forming radio galaxies at different redshifts (see also Sect. 4 for further discussion).

The last column in Fig. 9 gives the evolving spatial distribution of CRe injected by star formation for our B4 model, at three different epochs. The distribution of CRe injected by this mechanism has a low filling factor, and has overall a more patchy distribution compared to the CRe injected by the other two mechanisms. This is because the sources of the stellar feedback are more than sources of CRe from active AGNs, considering that stellar feedback is prominent in the majority of low-mass galaxies, which dominate the galaxy distribution.

2.7. Synchrotron emission

While the radio emission from cosmic shocks can be computed based on the prompt acceleration of relativistic particles (e.g. Hoeft & Brüggen 2007), here the presence of multiple injections from different processes require a more sophisticated treatment.

Ideally one should continuously update the CRe spectra under the effect of cooling and acceleration processes (e.g. Sarazin 1999); however, this approach is presently only doable for a few objects at a time, using Lagrangian tracers injected in Eulerian simulations (e.g. ZuHone et al. 2013; Wittor et al. 2017; Vazza et al. 2021c) or limited to the same number of gas resolution elements, like in Smoother Particle Hydrodynamics (e.g. Böss et al. 2023b,c).

Here, we explored a different approach, which has the advantage of allowing us to keep track of the time elapsed since the last injection of CRs in each simulated cell. In a nutshell, we applied pre-tabulated CRe spectra to predict the momentum distribution of CRe in all cells, based on their local value of density, magnetic fields, redshift and time since last injection of CRe, and use these spectra to predict the synchrotron radio emission. For the tabulated spectra, we used the ROGER Fokker-Planck solver² developed by Vazza et al. (2023) to produce a grid of template electron spectra and of their synchrotron emission spectra, in which we scanned a large range of discrete possible values of the magnetic field (10⁻⁹ G ≤ |B|≤10⁻⁵ G), gas density (10⁻⁵ ≤ n ≤ 10⁻² part/cm³), redshift (0 ≤ z ≤ 6), and age (0.01 Gyr ≤ t_age ≤ 10 Gyr) since the injection. A total of 276 480 possible combinations of the above parameters, using 10² momentum bins equally spaced in momentum bins and with up to 100 subcycling timestep, was used to produce a large library of momentum spectra for CRe, and for their synchrotron emission spectra (in the 50 MHz–5 GHz range). For every snapshot in which radio emission is needed, we selected for every cell in the simulation, the momentum and synchrotron spectra relative to the closest bin in density, magnetic field, age, and redshift. Then we normalised the CRe spectra and the synchrotron emission based on the actual CRe number density in the ENZO simulation. The same procedure was repeated separately for the three CRe fluid fields, which allowed us to compute, with a reasonable approximation, the radio emission produced by all CRe in the simulation. Examples of electron spectra generated with this approach, for a range of variations in the cosmic environment and epochs, are given in Fig. 10.

Some caveat on this procedure are worth noticing. First, it is reasonable to compute the emission from our CRe at a given epoch as long as the physical condition in the given cell did not change much since the first injection of CRe. In the presence of fast advection of CRe and/or of CRe injected since a long time, our solution for the synchrotron emission is likely to overestimate the true emission, as strong adiabatic losses, or the effect of stronger magnetic fields in the past, are underestimated. Our approach also neglects the role of shock and turbulent re-acceleration processes, which on the other hand may further stretch the high energy part of the CRe energy distribution on long timescales (e.g. Beduzzi et al. 2023).

In any case, CRe undergoing cooling losses converge onto the same evolving distribution for a large range of initial different parameters: after a first rapid (≤0.1 Gyr) cooling stage which mostly affects the highest energy tail of their distribution (usually γ ∼ 10⁴), the rest of the distribution settles to a slow evolution on long (∼1 − 10 Gyr) timescales, since collisional or adiabatic losses are negligible for the majority of the gas density values occupied by CRe (e.g. Sarazin 1999). Therefore, our approach gives an accurate answer only as long as the evolution of CRe is either dominated by recent injections (e.g. prompt shock emission or radio jets), or else for very long timescales during which loss processes dominate over re-acceleration terms.

3. Results

We divided our results in two main parts: first, a validation part where we show how our modelling compares with the basic known properties of galaxies; and second part where we highlight the important aspects of magnetic fields and cosmic ray electrons evolution that could be studied, for the first time, with our new simulations.

3.1. Properties of galaxies

The evolution of the cosmic star formation rate density (SFR), namely, the rate of star formation per unit time and volume, is given in Fig. 11. This is the first global statistics to assess our sub-grid implementation of star formation physics. While we show only the best results from our implementations and discarded all failed attempts in this work, we notice that even small parameter variations (i.e. the minimum mass of the formed star particles, the strength of the stellar feedback and the combination with different choices for the AGN feedback model) are found to significantly alter the total star formation across time. The simulated SFR is compared with the collection of observational data by Madau & Dickinson (2014).

Fig. 11. History of the simulated cosmic star formation density (i.e. normalised for Mpc3 comoving) in our suite of simulations, as a function of cosmic time. The grey points with error bars show the observed cosmic star formation derived in Madau & Dickinson (2014).

As a noticeable difference, the first three models (A1, A2, and A3) show a drop of the cosmic SFR at ∼5 Gyr since the start of the simulation (z ∼ 1.3), while the other four (B1, B2, B3, and B4) show a more prolonged star formation rate; in particular, models B3 and B4 roughly match also the star formation rate in the Local Universe. We deem model B4 to be the best performing here, because compared to B3 it can also best reproduce the peak of the cosmic star formation history at t ∼ 3 Gyr (z ∼ 2). The above differences can be understood because of the largest threshold mass used for star formation in the last four models, and the longer assumed dynamical time for star formation in the latest two in particular (B3, B4). As a result, less cold gas is consumed at high redshift, yielding a more realistic balance between star formation and the circulation of cold and hot baryons until the end of the run.

The last two models (C1 and C2) are not showed here for clarity, as in all galaxy properties they can just overlapped to model B4, as they employed exactly the same galaxy formation physics. It should be noted that our volume (42.5³ Mpc³) is still too small to be entirely free from cosmic variance effects, which can be dominated by the rare tail of high mass galaxies. In any case, the inefficient star formation in the low-z part of our simulated evolution is a deficiency in our models, which is likely to underestimate the amount of CRs injection from star formation in the local Universe. On the other hand, considering that the peak of the cosmic star formation rate is well captured in all these runs, and that the overall injection of CRs from galaxies is dominated by AGN feedback during most of the simulated evolution, we consider the latter an acceptable failure of our model, to be addressed in future work.

Next, we compute the distribution of the stellar mass and of the stellar mass fraction in all halos identified in the simulation, which is a more stringent test for the relative effects of star formation and stellar and AGN feedback across the full mass range of our galaxy distribution. The first three panels of Fig. 12 show the galaxy stellar mass distribution at three different epochs (z = 1.84, z = 1.065 and z = 0.02), which are compared to the corresponding best-fit relation obtained by analysing the Hubble Space Telescope observations in similar redshift bins by McLeod et al. (2021). The more prolonged star formation in B- runs implies an increased amount of halos with a large stellar mass compared to the first three models, and in general all models at z = 0.02 show a slight excess of massive galaxies compared to the best fit of observations (e.g. Behroozi et al. 2019). This is better shown by the last panel, which gives the final distribution of the stellar mass fraction as a function of the halo mass at z = 0.02 for our galaxies. All models have the tendency of overproducing the stellar mass fraction, although our ‘favoured’ model B4 has this problem only at the high mass end of the distribution. We notice that similar trends with halo masses are found also by recent simulations adopting more sophisticated prescriptions for galaxy physics, for instance, by Scharré et al. (2024) in the SIMBA simulation.

Fig. 12. Stellar mass luminosity function for our simulated galaxies in all models and for three different epochs, compared with the best-fit derived by McLeod et al. (2021) for observations at approximately equal redshift bins. The last panel shows the distribution of the stellar mass fraction as a function of the host halo mass at z = 0.02, compared to the best-fit relation derived from observations by Behroozi et al. (2019).

Far from being an attempt to self-consistently model the process of star formation in the Universe, our approach seeks an optimal combination of numerical parameters which can reasonably reproduce the observable properties of star formation on large scales, using the physical information available at the spatial and temporal resolution of our simulations. Based on these tests, we can conclude that with the right tuning of sub-grid parameters, despite its coarse spatial resolution our model can reasonably well capture the time evolution of the global SFR, as well as distribute the stellar component across a large range of halo masses, reasonably close to reality (with the exception of very large masses). Of all investigated models, the combination of parameters in run B4 can produce a realistic stellar mass distribution, a realistic history of cosmic star formation as well as a realistic correlation between black hole mass and halo baryonic mass. Therefore, this is the model in which our predictions for the total magnetic energy and the total amount of CRe injected by galaxy feedback can be considered as best anchored to real galaxies. Additional tests on the role played by our subgrid parameters for star formation and AGN physics on the simulated cosmic star formation history are discussed in Appendix A.2.

3.2. Properties of the simulated black hole population

Next, we computed the scaling relations between the properties of simulated SMBH, to test the result of our implementation of SMBH evolution and growth (Sect. 2.5). As we will do later also with other comparisons, we focus here on two representative examples of our models, namely, model A3 which is characterised by an intense AGN feedback activity in the past, and model B4 which presents a much more moderate one (runs C1 and C2 presenting only negligible differences in galaxy properties compared to run B4).

The first column in Fig. 13 gives the relation between the simulated black hole mass and the gas mass of the host galaxy for three sample epochs, in which the solid line gives the best fit of the observational trend, as reconstructed by Gaspari et al. (2019) for low-z objects, $M_{\mathrm{BH}}\propto M_{\mathrm{g}}^{0.57}$. As expected given our procedure, we can recover at all time steps the theoretical relation between the two quantities, meaning that the mass associated with our simulated SMBH is realistic at all epochs. It should be noticed that our points do not always stay exactly on top of the theoretical relation. The reason is that, in this plot, the gas mass is measured in post-processing, after the code has generated the output of one evolutionary timestep since the creation of SMBHs. Therefore, the gas mass is measured at that location after a small time delay (typically a few Myr) during which evolutionary processes of further growth and feedback have started affecting the surrounding gas mass. The second column of Fig. 13 gives the relation between the mass accretion rate onto SMBHs, and the third column shows the injected jet velocity from the same objects.

Fig. 13. Simulated scaling relations for the SMBHs modelled at run-time in our A3 and B4 runs, at three different epochs. The first column show the relation between the gas mass of the host galaxy and the SMBH mass, compared with the scaling relation inferred from observations by Gaspari et al. (2019). The second column shows the relation between the SMBH mass and the mass accretion rate onto the SMBH; the third column gives the relation between the SMBH mass and the jet velocity associated with feedback events.

The two simulated models clearly differ in their simulated AGN activity over time, with the A3 model being much more active for z ≥ 1, and the B4 producing more active SMBHs (albeit with a low matter accretion rate, ≤1 M_⊙/yr) down to z = 0. Such differences are entirely ascribed to the largest minimum SMBH mass assumed in run B4, which resulted into an overall lower power activity by SMBHs, but more prolonged in time, since a moderate amount of gas was allowed to remain closer to a typical SMBH until the end of the simulation. The impulsive vs more self-regulated AGN feedback of these two models is well captured by the distribution of Ṁ; it is systematically lower at all times in the B4 run, while in run A3 correspondingly higher jet velocities are reached, which clearly had the effect of expelling a larger fraction of gas away from halos already at high redshift.

Figure 14 gives the evolution of the mass accreted onto star-forming regions or onto AGNs, for the A3 and B4 cases. In the A3 model, there is an initial very efficient stage of stellar mass accretion, followed by a sharp drop due to the effect of stellar feedback, as well as by bursty cycles of AGN feedback. In the last Gyr of evolution, star formation is negligible while the AGN accretion is poorly self-regulated and characterized by episodic bursts. In the B4 model instead, the stellar mass accretion is more prolonged in time, and the balance between stellar and AGN feedback produces a more regular pattern of matter accretion onto AGN, as well as significant star formation also at later times.

Fig. 14. Evolution with cosmic time of the matter accretion rate onto star-forming particles (red lines) and of the matter accretion rate onto AGNs (blue lines) for our A3 and B4 models.

3.3. Radio emission from galaxies

We extracted the distribution of radio power emitted from our simulated radio sources, by generating their integrated emission at different epochs, based on the normalisation of the CRe density and of the local magnetic field strength, as well as on the average age since the injection of CRe in each cells (Sect. 2.7), and then by computing the emission within the projected area covered by the R₅₀₀³ of each galaxy identified in the simulation. Since CRe injected by AGNs or by star formation are differently tracked, we computed the emission produced by the two mechanisms separately, which further allows us to roughly compare with the distribution of emission from star-forming and radio galaxies.

Figure 15 gives the example of a crowded fields in our simulations at z ≈ 1, in runs A1, A2, A3, B1, B2 and B4 (B3 being extremely similar to B4). The contours mark the radio detectable emission from the CRe injected by AGNs at 150 MHz, assuming detection being larger than 3σ_noise where σ_noise = 1.05 μJy/arcsec² is of the order of the sensitivity of LOFAR for the LOTSS survey for a θ = 25″ beam. The colours give the mass-weighted gas temperature for the same volume instead. Several double lobe objects are visible in the image, often formed at the same location in all runs but with morphological differences (both in power and jets orientation) due to the stochasticity of the gas accretion events and launching jet directions from the simulated SMBH. By knowing the typical time elapsed since the injection by AGNs of the CRe emitting in each cell, as allowed by our method (see e.g. Fig. 5), we observe that all our detectable radio sources are ‘young’ (i.e. ≤0.1 − 0.2 Gyr). The latter thus represents the typical duty cycle of our simulated radio sources at low redshift, which is in good agreement with the typical duty cycle inferred from recent observations of radio galaxies in galaxy clusters and groups (e.g. Sabater et al. 2019; Biava et al. 2021; Brienza et al. 2022).

Fig. 15. Projected mass weighted gas temperature (colours, in units of log10[K]) and over imposed contours of detectable synchrotron radio emission at 150 MHz from CRe injected by AGNs in a subselection of our volume (about 10 Mpc across) in six runs at z ≈ 1.0.

The basic feature of a double lobed radio emission, loosely resembling Fanaroff-Riley type II radio galaxies in the cosmic web, is correctly captured by our implementation. However, the small-scale details and morphology are obviously poorly captured by the coarse spatial resolution of our runs.

Limited to the radio emission from CRe injected by star formation (only relevant in the faint luminosity end of the distribution) we introduced an ad-hoc fix to our procedure, to account for the fact that large fraction (i.e. from 60 to 80%) of the radio emission from galaxy halos is (likely) dominated by ‘fresh’ secondary electrons from the hadronic collision process (e.g. Werhahn et al. 2021). In reality, the continuous injection of secondary electrons well accounts for the observed flat distribution of radio spectral indices (α_R ≤ 1.0 where I(ν)∝ν^−α_R is the emitted radio power). In our case instead, the absence of the continuous injection of secondary CRe makes the population of CRe injected by star formation events dominated by z ≥ 1 events, which results into long average age for the CRe in the cells, and consequently too steep spectra associated with this population in particular. As a simple fix, we imposed an upper value of 10 Myr to their age (limited to this population) to compute their emission spectrum.

Figure 16 shows the evolving distribution function of the radio power of our simulated galaxies at 150 MHz, dividing the contribution from CRe injected by star formation (top panels) or by AGNs (bottom panels). In all cases, we compare with the recent results of LOFAR observations, derived by Cochrane et al. (2023) and by Kondapally et al. (2022) for star-forming galaxies or AGNs, respectively. Within the scatter related to the small simulated volume, the B4 model appears as the best in reproducing the real radio luminosity function of AGNs, across nearly four orders of magnitude in power both at z ≈ 0 and z ≈ 1, while all models underpredict the real AGN radio luminosity function at z ≈ 2. The comparison between the radio luminosity function of star-forming galaxies and of our star formation-related CRe is less satisfactory, with an under-prediction for low-power galaxies, as well as an over-prediction of high luminosity objects. The definition of radio luminosity functions based on physical processes is not straightforward, leading to biases in observed samples (e.g. Morabito et al. 2025). However, the mismatch of our simulated radio luminosity function of star-forming galaxies is larger than any realistic observational bias. Hence, it likely implies that our modelling of the injection of CRe by star formation processes is too crude (also considering the fact that we do not follow the production of secondary CRe, nor follow their diffusion and streaming through the multi-phase interstellar medium in halos and discs; see e.g. Armillotta et al. 2021; Hopkins et al. 2022).

Fig. 16. Top row: Evolving distribution function of radio power at 150 MHz from our simulated galaxies, in which only the radio emission from CRe injected by star formation is considered. Bottom row: Evolving distribution function of radio power at 150 MHz from our simulated galaxies, in which only the radio emission from CRe injected by AGNs is considered. The dotted grey lines in the top rows give the best-fit relation of observed data with LOFAR, in the same redshift bins, from Cochrane et al. (2023) (top) and Kondapally et al. (2022) (bottom).

One remarkable recent result from recent radio surveys (e.g. with LOFAR) is the hint of a relation between the radio galaxy morphological class (e.g. Fanaroff-Riley class I and II) and the specific star formation rate (i.e. star formation rate divided by the galaxy stellar mass) of host galaxies, while links between morphological class and stellar mass, environment, and so on, are more elusive to detect (Mingo et al. 2022). Inspired by this work, we measured (as shown in Fig. 17) the relation between the radio emission from our galaxies and their stellar mass (top panel) or their specific star formation rates (bottom panel) again in run A3 and B4. Here, we used the z = 1.065 epoch, which where most of the objects analysed by Mingo et al. (2022) come from. As in real observations, in both these models there is no clear relation between radio power and stellar mass (e.g. Fig. 16 left of Mingo et al. 2022). However, our simulated catalogues also do not show any clear correlated trend between radio power and specific star formation rate, which is in disagreement with the observations. While in principle the cosmic volume sampled by our simulations is rather small and the high power radio galaxies studied in Mingo et al. (2022) are not present here, this lack of correlation may suggest that additional variations of our jet feedback models, combined with a higher resolution and larger volumes, are necessary to reproduce reality.

Fig. 17. Relation between the emitted radio power at 150 MHz and the stellar mass (top panel) or the specific star formation rate (bottom panel) for our A3 and B4 run at z = 1.065.

While very valuable attempts have instead been produced using post-processing on the simulated distribution of galaxies (e.g. Thomas et al. 2021), our simulation (to the best of our knowledge) is the first in the real radio luminosity function of radio galaxies is reasonably well matched by the result of direct cosmological simulations, lunching magnetised jets and CRe at every feedback event.

3.4. Analysis of the network of galaxies

The schematisation of the continuous matter distribution of the cosmic web as a network, through its reduced representation by means of nodes connected by edges, has been recently proliferating in cosmological simulations. This has led to powerful new probes of cosmological parameters and galaxy clustering (e.g. de Regt et al. 2018; Tsizh et al. 2020, 2023; Vazza & Feletti 2020; Biagetti et al. 2021; Ouellette et al. 2023; Bermejo et al. 2024).

Here, we explore the application of network analysis to highlight potentially detectable differences between the distributions of radio galaxies and star-forming galaxies at z = 0.02, comparing model A3 and B4 as before. We consider whether a network analysis would be able to tell these two models apart, based on either the distribution of radio emission from star-forming galaxies or AGNs.

We used the three-dimensional (3D) distributions of identified matter halos for both runs, and with a standard procedure (see e.g. de Regt et al. 2018; Tsizh et al. 2020, for discussions) we built a connected network by assuming that all closer than a r = 1 Mpc distance are connected by an edge, or are disconnected otherwise. In order to test whether different radio emission components to the total radio luminosity functions are characterized by the same network, for each model we further doubled the catalogue of halos, and we alternatively removed the halos without detectable radio emission (using P ≥ 10²⁰ W/Hz at 150 MHz as a limiter for simplicity) from CRe injected by star formation or by AGNs.

We used the same algorithm by Tsizh et al. (2020) (to which we refer the curious reader for the full mathematical prescription of all parameters) to compute the average values of all most important network parameters, namely: (a) the average node degree, which is the average number of edges connected to a node; (b) the average clustering coefficient, which gives the fraction of connections between neighbours, compared to all possible connection between them, therefore quantifying the existence of structure in the local vicinity of nodes; (c) the average neighbour degree, that is: the average degree of the neighbours of a given node, normalized by number of neighbours; and (d) the average number of triangles, which gives the number of triangles, formed by edges and having nodes at the vertices, that include a given node. The list of the average parameters for the two models and their further subdivision into star-forming (SF) galaxies and AGN radio galaxies is given in Table 2.

A first striking difference to notice is that (regardless of the model) the networks traced by SF galaxies or by AGNs are very different in all parameters, with the second showing a smaller number of nodes and a higher degree of clustering and density of nodes. This fully reflects the approximate mass segregation of the two effects (with star formation, and its related radio emission, being more relevant in the low mass end of the galaxy mass distribution) and it echoes the topological mass bias (e.g. Bermejo et al. 2024), according to which halos with different mass ranges of halos samples give a biased view of the underlying cosmic network.

Furthermore, the analysis shows that the distribution of radio AGNs in model A3 is significantly denser than the distribution of radio AGNs in the B4 model (as marked by the much higher average node degree, neighbouring degree and triangles), despite the number of nodes is significantly larger (by 17%) in model B4. On the other hand, the network traced by the radio emission by the SF galaxies has ∼3 more nodes than the network of AGN galaxies in both models, but it is significantly less dense than the latter, with average parameters varying by less than ∼10% in all cases.

This first pilot exploration of the simulated network of radio emitting galaxies already shows the potential of our suite of simulations to, both, assess the biases associated with the use of specific samples of radio emitting galaxies (considering that this will be one of the main aims of the Square Kilometre Array, e.g. Prandoni & Seymour 2015), and it also highlights the impact of different prescriptions for galaxy evolution on the reconstructed network parameters.

3.5. Volume filling properties of magnetic fields and of cosmic ray electrons

The interpretation of several key observational probes of primordial magnetic fields crucially depends on the accurate removal of the possible magnetic field contamination by astrophysical processes (e.g. Vazza et al. 2017, 2021a; Alves Batista & Saveliev 2021, for recent discussions and tests). A simple way to assess the possible contribution of primordial and astrophysical mechanisms to the magnetisation of the Universe is to compute the volume filling fraction of magnetic fields produced by a given process. For example, Tjemsland et al. (2024) have recently estimated that a volume filling factor f ≥ 0.67 with a magnetic fields strength B ≥ 10⁻¹⁵ G is required to explain the observed lack of Inverse Compton Cascade emission from blazars at ≥GeV energies (e.g. Neronov & Vovk 2010).

The volume filling factors of magnetic fields for our runs at at the usual three epochs are shown in Fig. 18. The solid lines refer to all cells in the simulation, and show the very clear effect of primordial magnetic fields (runs C1 and C2) compared to all other models: strong enough primordial magnetic fields (3.7 ⋅ 10⁻¹⁰ G in this case) can clearly dominate the magnetisation of virtually ∼100% of the cosmic volume (with the exception of the tiny volume fraction occupied by galaxies). In runs with negligible primordial field (=10⁻¹⁶ G), on the other hand, astrophysical sources can dominate the magnetisation of a significant fraction of the cosmic volume. To quantify how much, and since there still is a minimal amount of primordial magnetisation in all our runs, we use our tracking of CRe fields to compute exactly which fraction of the cosmic volume has been magnetised by galaxies in all runs. Of course a fraction of the magnetic fields in these cells will also be of primordial origin, but by construction these regions were mostly contaminated by feedback processes.

Fig. 18. Volume-filling factors of comoving magnetic field strength in our runs at three different epochs. The solid lines give the cumulative distribution for all cells in the simulation, while the dashed lines show the contribution from magnetic fields only injected by ‘astrophysical’ processes (i.e. AGN feedback and star formation).

Therefore, in Fig. 18 we also mark the volume filling factor only of cells affected by feedback in our runs with dashed lines. At the reference value of 10⁻¹⁵ G, all models without a strong primordial magnetic field have a filling factor f < 20%, in line with Tjemsland et al. (2024), with the exception of the A3 model which has f ∼ 35% due to its enhanced AGN activity. These filling factors are reached only for z ≤ 1, as the magnetisation bubbles emerge from their host galaxies and start percolating only at evolved cosmic epochs At 10⁻¹² G, the filling factor is f ≤ 15% again with the exception of run A3 (f ∼ 37%), which is in line with similar statistics computed by Arámburo-García et al. (2021) for Illustris TNG simulations (f ∼ 11 − 15%). In all cases, such low values of f supports the idea that blazar observations can be a powerful probe of primordial magnetism (e.g. Neronov & Vovk 2010), because no plausible feedback scenario could magnetise a large enough fraction of the Universe’s volume to compete with primordial magnetic fields (e.g. Bondarenko et al. 2022; Tjemsland et al. 2024).

A similar volumetric analysis of the distribution of CRe with respect to thermal gas particles is useful to assess how many fossil electrons are available across the cosmic web, as a result of all simulated processes. The top panel of Fig. 19 gives average ratio between CRe density (of all species of injected CRe) and the thermal gas proton density as a function of the gas overdensity in the simulation by the end of our runs.

Fig. 19. Top panel: Average ratio between CRe density and thermal gas proton density as a function of the gas over-density in the simulation by the end of our runs. Middle panel: Cumulative absolute number of CRe normalised for a comoving Mpc3 volume and averaging over an increasing gas overdensity. Bottom panel: Average time elapsed since the last injection of CRe as a function of gas overdensity. The solid lines give the trend for CRe injected by shocks, the dotted lines for CRe injected by AGNs and the dashed lines for CRe injected by star formation.

A very general finding is that, regardless of model-to-model variations, at large overdensities the energy density of CRe is always dominated by the time-integrated injection by shocks (solid lines). The relative abundance of models is similar to the trend of the luminosity functions (Sect. 3.3), in which the smallest amount of CRe injected by AGNs was found for model A3, and the highest for model B4 Interestingly however, the abundance of CRe injected by AGNs in the A3 model is the highest of all, but only when gas overdensities ≤1 are included: this means that the stronger feedback activity found in the A3 model at z ≥ 1 has expelled a lot of gas enriched with CRe in the extreme periphery of halos, or beyond. This trend is also consistent with the similarly high trend of the volume filling factor of magnetic fields, found above.

Shocks are instead found to dominate the injection of CRe at all overdensities, with increasing importance with decreasing density. While all runs without strong primordial magnetic fields display a drop in the ratio of CRe density to thermal gas density, owing to the fact that shocks are almost never formed in the low density range typical of voids, in our C1 and C2 run with primordial stochastic magnetic fields there is a non negligible injection of CRe also at very low density, due to a population of shocks driven by the dynamics of magnetic field fluctuations in voids, as first reported in our previous work (Vazza et al. 2021b).

Star formation is never the dominant injection mechanism of CRe, even if its ratio compared to the AGN injection depends on the details of the implemented galaxy formation prescriptions, and can be more important than shocks for a large, ≥10² overdensity. Considering the limited resolution of our runs and the simplistic adopted star formation model, this result may need further confirmation with more sophisticated galaxy formation models in the future.

The middle panel of Fig. 19 gives the cumulative absolute number of CRe, injected by the different mechanisms, normalised for a comoving Mpc³ volume and averaging over an increasing gas overdensity. This shows that through most of cosmic volume the budget of CRe is dominated by shock acceleration. While a magnetic obliquity dependant acceleration efficiency (see also Sect. 3.7) does not change this picture, the total number of CRe injected by shocks, especially at low densities, is affected by the additional shock dynamics induced by galaxy formation-driven winds and outflows (e.g. Kang et al. 2007; Schaal et al. 2016).

The ratio between the number of CRe injected by various mechanisms within the R₅₀₀ of halos in the volume, and the number of thermal gas protons within the same volume is given in Fig. 20 for the three epochs considered already elsewhere. By the end of the simulation, the number density of CRe is dominated by the contribution from shocks. On the other hand, the CRe injected by star formation is small and only competes with the other two at z = 1.84, namely, close to the peak of the cosmic star formation history in all models. The trend with halo mass is approximately constant for shocks injected CRe, while we report a clear increase with mass of the abundance of CRe injected by AGNs, consistent with the dominance of this feedback mechanisms in high mass halos. In summary, in all simulated halos for z ≤ 2 the number density of CRe is ∼0.02 − 0.1% of the total number density of thermal protons, and it is dominated by the production at shocks. Once more it is worth stressing that total number of CRe includes very low energy electrons, with momenta a few times above the thermal momentum of protons (∼1 − 10 MeV), as well as super-relativistic electrons typically emitting synchrotron radiation in the radio band (∼1 − 10 GeV), or beyond, with energy distribution depending on the local spectral evolution.

Fig. 20. Distribution of the ratio between CRe density and baryon density within R500 of our halos at three different epochs and for all runs. The solid lines show the distribution of CRe injected by shocks, the dotted lines by AGNs and the dashed lines by star formation.

Finally, the bottom panel of Fig. 19 gives the average time elapsed since the last injection of CRe, as a function of environment, separately for the three families of CRe and for all runs. This shows in a more systematic way the trends we have highlighted already: even at z = 0.02 the injection by shocks is a very active process, and the average time since the last shock injection is ≤0.5 Gyr at all densities. CRe injected by AGNs have instead average ages in the ∼1 − 3 Gyr range, with variations depending on the assumed routine for triggering AGN activity. Instead in all our models, the average age of CRe injected by star formation is much higher; namely, ∼2 Gyr within halos. This increases to ∼10 Gyr in voids and filaments, which were contaminated only in the very first star formation episodes that could have driven some amount of CRe via outflows, before the full formation of cosmic structures.

3.6. Radio properties of the cosmic web

We studied some potentially detectable signatures of the diffuse magnetic fields and CRe injected by the different mechanisms in the radio band. First, the Faraday rotation measure (RM) of polarised extragalactic sources, integrated over long (≥Gpc) cosmological LoS, and removed by the contamination from the Galactic foreground, is a powerful probe of cosmic magnetism, thanks to the sensitivity of modern radio surveys (e.g. Vernstrom et al. 2019; O’Sullivan et al. 2020; Arámburo-García et al. 2021; Carretti et al. 2022; Pomakov et al. 2022; Carretti et al. 2023; Mtchedlidze et al. 2024). We thus produced mock RM statistics of cosmological LoS through the simulated volume of our runs, to constrain the observable effects of our implemented feedback physics on the magnetisation of large scale structures, and compare it to the contribution from purely primordial magnetic fields. Following a procedure already introduced elsewhere in the comparison with real LOFAR data (Carretti et al. 2023, 2024), we built 100 LoSs through the simulation, by replicating the simulated volume 153 times, from a set of 17 snapshots saved at nearly equally spaced redshifts, in order to sample the distribution of thermal electron density and of the magnetic field component parallel to the LoS, from z = 2.98 to z = 0. Each full LoS is ≈6.502 comoving Gpc (containing 156 672 cells in total) long and it is assembled by randomly varying the volume-to-volume crossing position to minimise the production of periodic artefacts, moving along the same Cartesian axis of the simulation. Additionally, all cells in the simulation with a gas matter over-density larger than 50 times the cosmic average gas density were flagged and excluded from the integration to minimise the RM contamination by cluster halos. This is meant to mimic the more complex procedure applied to real LOFAR observations, where the information of spectroscopic sample of more than 150 000 halos from Wen & Han (2015) is used to excise the LoS crossing too close to the possible contamination by galaxies (see Carretti et al. 2023, and discussion therein).

For all LoSs, we integrated

$$\begin{array}{c}\hfill {\mathrm{RRM}}_f[\mathrm{rad}/{\mathrm{m}}^2]=0.812{\displaystyle {\int }_z^0}\frac{n_{\mathrm{e}}{B}_{\Vert }}{{(1+z)}^2}\mathrm{d}l,\end{array}$$(4)

where n_e [cm⁻³] is the density of thermal electrons, dl [pc] is the differential path length and B_∥[μG] is the magnetic field component along the integration path.

The panels of Fig. 21 show the example of the full RRM_f integrated within the snapshot of a 42.5³ Mpc³ box at z = 0.02 for all our runs. With the exception of the C1 and C2 runs, which included a primordial stochastic magnetic field (with spectrum P_B ∝ k⁻¹ and normalisation B_1 Mpc = 0.37 nG) in the initial conditions, all models produce a similar amplitude of Faraday Rotation, in all cases strongly enhanced in the proximity of dense halos. More subtle threads of enhanced Faraday Rotation can be seen also in filaments in runs C1 and C2, highlighting the larger volume filling factor of magnetic fields in primordial models and the fact that only primordial magnetic field models are capable of significantly magnetise cosmic filaments. Such differences are best highlighted when integrated for the full long LoS, and after excising the contamination by halos, as in Fig. 22, in which we computed the full RRM_f integrated up to z = 3 limited to the interesting comparison of models A3, B4 and C1. This is motivated by the fact that, in all previous tests, model A3 stood out as the one producing the largest volume filling factor of astrophysical magnetic fields in voids. Hence it is the one in our models which can be most competitive with primordial magnetic fields. However, for reasons discussed above, the amount of star formation and feedback developed in A3 are unrealistically high, while they are more realistic in B4. The latter, featuring a more prolonged activity by star formation and AGNs at low redshifts, produces higher values of RRM_f(z) for z ≤ 0.5, while the more bursty A3 model dominates the rest of RRM_f(z) for 0.5 ≤ z ≤ 2.5. In Fig. 22 we compare our simulations to the recent measurements of residual Rotation Measure obtained using LOFAR by Carretti et al. (2024). The latter include the subtraction, at all redshifts, of a local intervening astrophysical contribution at the location of polarised sources (which cannot be modelled in our simulations), with an assumed, ∝1.1/(1 + z) [rad/m²] redshift trend as motivated in Carretti et al. (2024). In all cases, both purely astrophysical scenarios for the seeding of extragalactic magnetic fields produce values of RRM_f(z) much lower, at all redshifts, than LOFAR observations. Model C1, which also features the stochastic primordial magnetic field, yields instead values in the ∼1 − 3 rad/m² range for z ≥ 0.5, and thus it provide the best reproduction of observed LOFAR among our models (also including the additional primordial models tested in Carretti et al. 2024).

Fig. 21. Projected RM maps (in units of [rad/m2]) across the entire simulated volume for our seven runs at z = 0.02. Each panel is approximately 21.25 × 42.5 Mpc2 large.

Fig. 22. Trend of simulated rms of the residual Faraday RM, observed at z = 0 and for integrating out to increasingly larger redshift (up to z = 2.5), for our A3, B4, and C1 runs. The lines give the mean dispersion of the RRMf for 100 randomly selected lines of sight through the simulated volume, while the black line with larger error bars shows the observed RRM trend by LOFAR, taken from Carretti et al. (2024).

Next, we compute the synchrotron radio emission from CRe accelerated by cosmic shocks, which gives another promising way to detect the cosmic web, especially relevant for the next generation of radio observations, culminating with the Square Kilometre Array (SKA e.g. Brown et al. 2011; Vazza et al. 2015; Weltman et al. 2020; Hodgson et al. 2021; Oei et al. 2022).

The panels in Fig. 23 show the projected synchrotron radio emission from the CRe accelerated by shocks in all our runs at z = 0.02, with additional red contours showing the detectable radio emission, assuming as before ≥3σ_noise detection where σ_noise = 1.05 μJy/arcsec² is of the order of the sensitivity of LOFAR for the LOTSS survey for a θ = 25″ beam. Once more, the much more volume filling distribution of magnetic fields predicted for the primordial models in run C1 and C2 greatly increases the average synchrotron emissivity of filaments, which produce a nearly fully connected radio cosmic web across scales of tens of Megaparsecs. However, the bulk of the emission from CRe accelerated by shocks around filaments remains a factor ∼10 or more below the LOFAR HBA detection threshold, with the occasional exception of the peripheral regions of some halos, and of a few very few short-range intra-group bridges, which are barely detectable with LOFAR HBA here. This is consistent with present constraints from recent radio observations: the stacking detection of (likely) cosmic filaments by Vernstrom et al. (2021, 2023) have been obtained using up to ∼6 ⋅ 10⁵ pairs of physically connected halos of luminous red galaxies. Based on this, the imaging detection of filaments with LOFAR HBA should be possible only with a ∼25 better sensitivity, that is, $\sigma {\prime }_{\mathrm{noise}}=0.014\mathrm{\mu }$Jy/arcsec² at 150 MHz for a θ = 25″ beam. This estimate is also in line with the latest limits from the non-detection of filaments in LOFAR HBA using a smaller sample of massive filaments, by Hoang et al. (2023).

Fig. 23. Projected mass weighted gas temperature (colours, in units of log10[K]) and over imposed contours of detectable synchrotron radio emission at 150 MHz from CRe accelerated by shocks in a subselection of our volume (about 20 Mpc across) in our nine runs at z = 0.02.

3.7. Role of shock obliquity

Finally, we also tested the impact of an obliquity-dependent acceleration of CRe in shocks, which can be assessed by comparing our runs C1 (no obliquity dependence) and C2 (injection of CRe by shocks only for θ_B ≥ 45°).

Both the projected maps of radio emission (not shown) and the properties of CRe distribution within halos (Fig. 24) show that in general the differences of the two models are extremely small, with the C2 model having only a tiny (≤10%) reduction in the amount of injected CRs within halos of all masses, compared to the C1 model. This is consistent with previous results from the few numerical simulations which investigated the effect of shock obliquity in cosmological simulations. They reported that most of the energy dissipation in the cosmic web happens through the oblique shocks (Banfi et al. 2020; Böss et al. 2023b), due to the fact that, within halos, merger shocks sweep a quasi-random distribution of magnetic field orientations (Wittor et al. 2017), while the magnetic field around accretion shocks in filaments and outside clusters is predominantly aligned with the shock surface, as an effect of the local gas dynamics (Banfi et al. 2021). Finally, very recent kinetic plasma simulations by Gupta et al. (2024) have reported large levels of acceleration of relativistic electrons even by quasi-parallel shocks, due to the onset of resonant and non-resonant streaming instabilities induced by the acceleration of cosmic ray protons. Even though they limited to one dimension, the implications of these models is that the acceleration of electrons by cosmic shocks might be larger than expected, solely based on diffusive shock acceleration. In this case, our approach is flexible enough to allow us to test different CRe injection efficiencies than the one tested here (Sect. 2.4).

Fig. 24. Distribution of the ratio between shock-injected CRe density and baryon density (bottom row) within R500 of our halos for the C1 and C2 runs at the end of the simulation.

4. Discussion

We comment here on a few important physical processes that are not included in our models. First, our simulations neglect the effects of cosmic-ray diffusion or cosmic ray streaming in our simulation. Predicting the diffusion of CRs in a tangled magnetic field is a non trivial problem in astrophysics (e.g. Bhattacharjee 2000; Lazarian et al. 2023; Harari et al. 2014, for reviews). However, even the most energetic CR electrons contributing to the energy density of relativistic particles in the cosmic web (∼10 GeV) can only diffuse over a spatial scale that is several orders of magnitude smaller than the length scales typically covered by advection during the same time. A diffusion coefficient of cosmic rays in the range of D ≤ 10³⁰ − 10³¹ cm²/s has been derived for the intracluster medium (e.g. Brunetti & Jones 2014, and discussion therein). this implies that the diffusion timescale over a L_D ≈ 1 Mpc scale is much longer than the Hubble time, and longer than any other relevant energy loss process for these cosmic rays: τ_diff ∼ L_D²/(4D) = 7.5 − 75 Gyr. Moreover, considering that CRs sourced by galactic evolution processes are injected in highly magnetised regions (AGNs or star-forming regions), their diffusion length scale is very likely to be even smaller than this (because L_D ∝ 1/|B|) indicating that the effect of CR diffusion can be reasonably neglected compared to the one of advection, for the sake of our analysis. However, on the scale of galaxies the effects CR diffusion and streaming can be relevant to correctly reproduce their radio emission. Our simulations are not designed to match observations on those spatial scales (≤40 kpc). There the modelling of CR transport, diffusion and streaming in the multi-phase interstellar medium, as well as the interplay between primary and secondary cosmic rays (e.g. Evoli et al. 2018; Armillotta et al. 2021; Hopkins et al. 2022; Girichidis et al. 2022) adds complexity to our simplistic treatment and will deserve more investigation in the future.

Second, the recipes of CRe injection adopted in this work are rather simplistic, and can made more sophisticated in several ways and based on the (ongoing) exploration of shocks in collisionless plasmas which particle-in-cell simulations make possible. While we already explored the role of shock obliquity in the budget of accelerated CRe (Sect. 3.7), different efficiencies and injection spectra from strong shocks have been proposed in the literature (e.g. Caprioli 2012; Xu et al. 2020). Moreover, more complex and realistic theoretical recipes to model the injection of electrons from low-Mach number and oblique shocks, in which shock-drift acceleration can operate to pre-accelerate electron into DSA,can be considered (e.g. Amano & Hoshino 2022). While this is not expected to change our estimates of radio emission from filaments and cluster outskirts, where the direct injection of shocks should be dominant and the magnetic obliquity should be close to perpendicular (e.g. Banfi et al. 2020), it might increase the level of emission from weak internal shocks driven my mergers in the intracluster medium, where the expected DSA efficiency is low (Guo et al. 2014a).

Our simulations also neglected (see discussion in Sect. 3.3) the continuous injection of secondary CRe by hadronic collisions, which are surely important within galaxies (e.g. Werhahn et al. 2021). Although this mechanism is crudely taken care of in our estimate of radio emission (Sect. 3.3) and it is deemed to be generally subdominant for the seeding of CRe on large scales (compared to AGNs and shocks), future works may include also this additional seeding mechanism for CRe at run-time.

Finally, another caveat is the role of spatial resolution in our predictions for CRe and magnetic field-related statistics. While in previous works, we already tested whether the statistics of cosmic rays by accretion-driven shocks is well converged already for ≤100 kpc resolution (e.g. Vazza et al. 2011, 2012), it remains to be assessed with further studies whether the volume-filling factor of magnetic fields and CRe injected by galaxies are equally well converged or otherwise.

5. Conclusions

We introduce new ENZO cosmological simulations to study the injection and the evolution of relativistic electrons in the cosmic web. We focused on their number density and on their average age and we simulated their injection by three main processes: structure formation shocks, AGN feedback, and star formation. Our procedure tracks the evolution of relativistic electrons with typical minimum energies of ∼1 − 10 MeV and accelerated up to ∼1 − 10 GeV or beyond. The spectral energy distribution (SED) can be reconstructed in post-processing, thanks to the combination of our spectral ageing model (Sect. 2.7) and of the additional tracking of the average time elapsed since the injection of cosmic rays into any particular cell (Sect. 2.3).

Given the spatial and mass resolution limits of our simulations, as well as of the approximations used to mimic the growth and feedback from a SMBH, we calibrated sub-grid physics to reproduce known observable relations; rather than developing a full physical model to solve galaxy-scale physics or the formation and growth of SMBHs. We have shown that it is possible to calibrate this procedure to aptly reproduce several key observations (e.g. the stellar mass distribution in galaxies, the correlation between SMBHs and host galaxy masses and the radio power distribution function of galaxies) and to use this information to robustly predict the amount of CRe injected in large-scale structures by cosmic evolution. In particular, we report the first and very good reproduction of the luminosity distribution of radio emission from galaxies in the simulation and up to z ∼ 2, using a distribution of magnetic fields and CRe evolved at run-time (focusing on the number density of CRe and tracking their average age with our simple exponential-decay age tracing model) and using ideal MHD simulations. Our main results concerning the evolving distribution of these two important non-thermal components can be summarised as follows:

• Structure formation shocks dominate the injection of CRe in the majority of the Universe: Even for conservative assumptions on the CRe acceleration efficiency from diffusive shock acceleration (calibrated to reproduce known radio emission on the scale of clusters of galaxies), shocks dominate the volume-filling population of CRe in the cosmic web, as well as the total amount of injected CRe in all investigated models. This finding is not altered by including the obliquity-dependent injection efficiency by shocks.

• AGN feedback dominates the injection of CRe with respect to stellar feedback both within halos and also when volume averaged statistics are considered: While stellar feedback starts to enrich the Universe earlier with CRe and magnetic fields, the more powerful activity by AGNs until late cosmic times end up dominating the overall production of CRe, both in high- and low-mass halos.

• Known astrophysical sources can only magnetise a tiny fraction of the Universe: By computing the filling factor of magnetic fields in the different runs, we estimate f ≤ 15% of the Universe at all redshift can be realistically magnetised above 10⁻¹⁵ G by galaxies (both including AGN feedback and star formation-driven outflows); thus, any primordial magnetic field larger than this would dominate observations. Even for our most extreme feedback model (A3), the resulting filling factor (∼37%) is not enough to account for the non detection of inverse Compton cascade emission from blazars or to explain the trends of Faraday RM obtained with LOFAR.

Our new numerical framework (which is scalable to larger, or more resolved simulations) represents a first important step towards the self-consistent coupling of galaxy formation processes to their observable or invisible non-thermal output. While this approach can still be subject to improvements on the side of the implemented galaxy evolution physics, this is key to investigate the important role played by several sources in seeding of large-scale emissions with a ‘fossil’ population of mildly relativistic electrons, which Fermi I and Fermi II processes can further re-energise (e.g. Brunetti & Jones 2014, for a review).

Based on the results of this work, we are now able to quantify the total budget of CRe injected within halos by all mechanisms considered in our simulations by z = 0. The budget of CRe injected by each mechanism scales with the assumed acceleration efficiencies: ξ_{e, ℳ ≥ 5} is the CRe injection efficiency for strong shocks, ξ_AGN is the CRe injection efficiency assumed for AGN feedback, and ξ_SF is our assumed CRe injection efficiency from star formation (averaged over our 41.5³ kpc³ resolution). There are, of course, several other dependencies on the assumed physical implementations of AGN physics and SMBH growth, star formation, and feedback, which cannot be easily prescribed in this way due to the non-linearity of their couplings.

In order to give an order of magnitude estimate of the various contribution to CRe from all considered mechanisms within matter halos, we use our best model B4 (and based on Fig. 20), we can parametrise the most important model dependencies as follows:

$$\begin{array}{c}\hfill \frac{N_{\mathrm{CRe}}}{N_{\mathrm{th}}}=5.51·{10}^{-4}\frac{{\xi }_{\mathrm{e},\mathcal{M}\ge 5}}{4.6·{10}^{-4}}+4.9·{10}^{-5}\frac{{\xi }_{\mathrm{AGN}}}{2·{10}^{-4}}+6·{10}^{-6}\frac{{\xi }_{\mathrm{SF}}}{{10}^{-5}},\end{array}$$(5)

for our largest, ∼5 ⋅ 10¹³ M_⊙ halos at z = 0, and

$$\begin{array}{c}\hfill \frac{N_{\mathrm{CRe}}}{N_{\mathrm{th}}}=6.26·{10}^{-4}\frac{{\xi }_{\mathrm{e},\mathcal{M}\ge 5}}{4.6·{10}^{-4}}+2.1·{10}^{-5}\frac{{\xi }_{\mathrm{AGN}}}{2·{10}^{-4}}+2.1·{10}^{-6}\frac{{\xi }_{\mathrm{SF}}}{{10}^{-5}},\end{array}$$(6)

for the smallest ∼5 ⋅ 10⁹ M_⊙ halos at z = 0. While the details of their star formation rate density or galaxy luminosity functions, may differ significantly at specific redshifts, their final integrated budget of CRe only differs by ∼60%. The total amount of injected CRe in large-scale structures appears as a robust prediction of the model, despite the several assumptions and open theoretical problems related to galaxy physics and to the evolution and formation of SMBHs.

It is now important to ask whether these amount of CRe are sufficient to account for the budget of fossil radio emitting electrons required by the modelling of diffuse radio emission in clusters of galaxies. In a recent review (Vazza & Botteon 2024), we computed the number of CRe that must be contained within the volume of observed radio halos at the centre of clusters of galaxies, as implied by simplistic assumptions of energy equipartition between cosmic rays and magnetic field, and for a variety of observed radio spectra and (unknown) ratio between the energy density of CR protons and electrons. For a M = 10¹⁵ M_⊙ cluster of galaxies, between N_CRe ∼ 5 ⋅ 10⁶⁶ and N_CRe ∼ 10⁶⁸ electrons within a Mpc³ volume are required to explain the observed radio emission of a medium power radio halo; this depends on the spectra energy distribution of the electron population and on their minimum energy, as well as on the amount of CR protons (which enters the derivation of the equipartition magnetic field).

Our new simulations show that the amount of CRe in massive halos at low redshift is generally in line with the above estimates and, thus, that the combination of the CRe seeding mechanisms explored in this work is more than enough to fuel Fermi processes in their production of diffuse radio emissions in clusters. Moreover, the same approach will also allow us to better understand the acceleration physics in more peripheral and less dense regions of the cosmic web, which were recently probed by recent radio observations (Govoni et al. 2019; Botteon et al. 2020, 2022; Vernstrom et al. 2021, 2023; Cuciti et al. 2022), considering that it is presently unclear how much fossil CRe, or magnetic fields, are present at such large scales. This appears timely and crucial, considering that a systematic survey of these regions will be one of the future main goals of the SKA (e.g. Weltman et al. 2020).

Data availability

Movies associated to Figs. 1 and 9 are available at https://www.aanda.org

Acknowledgments

We thank our anonymous reviewer for the valuable help in improving the quality of this work. F.V. acknowledges the CINECA award “IsB28_RADGALEO” and “IsCc4_FINRADG” under the ISCRA initiative, for the availability of high-performance computing resources and support. F.V. also gratefully acknowledges the usage of computing at the CSCS-ETH supercomputing centre (Lugano, Switzerland) on Piz-Daint under project ID s1096. F.V. and S.M.’s research has been supported by Fondazione Cariplo and Fondazione CDP, through grant n° Rif: 2022-2088 CUP J33C22004310003 for “BREAKTHRU” project. F.V. also acknowledges financial support from the European Union’s Horizon 2020 program under the ERC Starting Grant “MAGCOW”, no. 714196. MB acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306 and the DFG Research Group “Relativistic Jets”. M.T. was supported by the National Research Foundation of Ukraine under Project No. 2020.02/0073. M.T. would like to thank the Armed Forces of Ukraine for providing security to perform this work. We acknowledge fruitful scientific discussions with Y. Dubois, G. Brunetti, V. Heesen, A. Botteon, M. S. S. L. Oei and H. Rottgering.

References

Appendix A: Test re-simulations with physical variations

Here, we present several important numerical tests in which we resimulated the same cosmic volume (smaller than the one used in the main paper) to compare the outcome of important physical or numerical variations of our algorithms. This is meant to assess the robustness of the conclusions of the main paper and their dependence on the many sub-grid physical prescriptions for galaxy formation. In detail, we used a 21.25³Mpc³ volume using 512³ cells and 512³ dark matter particles, hence these test simulations cover the same spatial and mass resolution of the runs employed in the main paper, but have a 8 times smaller volume. All explored variations of parameters are given in Table A.1, and in the next subsections we comment on each of these variations in more detail.

A.1. Tests on the relativistic electrons modelling

First, we tested the robustness of our new proposed method to track the advection and ageing of relativistic electrons in ENZO cosmological simulations (Sect. 2). While the specific choice of the timescale τ used for the age reconstruction in Eq.1 is arbitrary and, a single injection of CRe would always allow for perfect recovery with this technique, multiple injections complicate the process as older or younger population of CRe may in principle differently weight the resulting age, depending on the choice of τ. We thus tested the use of τ = 10 Myr (model test_time001 in Table A.1), τ = 100 Myr (as in the main paper, model testA) and τ = 1 Gyr (model test_time1), in all cases using the same physical parameters of run B3. Figure A.1 gives the projected CRe-weighted maps for the time elapsed since the last injection of CRe, separately for the three different families of CRe and for the two extreme 10 Myr and 1 Gyr cases. Besides small random differences (due to our implementation of random axis orientations for feedback events), the distributions of ages look extremely similar in both cases, for all considered CR species. We thus conclude that, for the distribution of ages concerned with the physical evolution of our simulations, the choice of this parameter is not critical for a correct a-posteriori reconstruction of the time elapsed since the injection of CRs.

Fig. A.1. Projected mean (CRe-weighted) age of CRe since their last injection (in units of [Gyr]) by shocks (top row), AGNs (middle row), or star formation (bottom row) for two control resimulations at z = 0.02, which employed τ = 0.01 Gyr (model test_001, left panels) or τ = 1 Gyr (model test_time1, right panels) to solve Eq.1. Each panel is 21.25 Mpc across. No significant differences in the reconstructed age of the different CRe population are found.

A.2. Tests on variations of the AGN modelling

Our run-time procedure for AGN feedback, without relying on actual SMBH Lagrangian sink particles, is based on a few assumptions to generate realistic injection sites of feedback events. As detailed in Sect. 2.5, we identify at run time local high baryonic density peaks (100 times larger than the baryonic mean density), and there we measure the total mass contained within a fixed comoving radius (2 cells, i.e. 83 kpc). We also adopted a local exclusion criterion to avoid generating multiple SMBHs within a radius of three cells (124.5 kpc) from the central gas density peak.

With a test run (test_AGN_sep in Table A.1), we resimulated the same 21.5³Mpc³ volume using the same local exclusion criterion of 3 cells as in the baseline simulation of the main paper, or using four cells instead (so 166 kpc from the central gas density peak). The top panel of Fig. A.2 shows the integrated star formation history in this volume, which does not show significant differences in the two cases. Other quantities related to this (i.e. injection of CRe by star formation or AGNs, halo scaling relations, etc) equally show negligible differences, suggesting the adopted minimum separation between SMBH host cells in our run-time procedure is not critical.

The same figure also shows the effect of other variations related to our SMBH implementation and AGN feedback: we varied the values of the α_B "boost" parameter used in the Bondi formalism to compute the matter accretion from the hot and the cold gas phase (Sect. 2.5) in models testB and testC, or we increased the minimum density (ρ_AGN) to form SMBHs at our run-times (testD and testE). The effects of all variations in the star formation history is minimal, while our fiducial choice of parameters ρ_AGN = 20, α_B, cold = 1 and α_B, cold = 10 (e.g. in the B4 model of the main paper) overall gives the best results when coupled to the prediction of luminosity functions of radio galaxies (see Sect. 3.3).

A.3. Tests on variations of the stellar formation modelling

With other tests, we varied the dynamical timescale associated with the formation of stellar particles, which is a fundamental parameter in our model (Sect. 2.6), while keeping instead the AGN feedback prescription fixed. In this case the effects are more evident: a short t_* = 1 Myr timescale allows the fast formation of a large quantity of stars, reaching a peak already ∼1 Gyr after the start of the simulation, followed by a steady decline. Conversely, a much longer t_* = 100 Myr timescale delays the formation of the first stars to ∼3 Gyr after the star of the simulation. As shown by the bottom panel of Fig. A.2 both trends of cosmic star formation are quite far from the observed ones, and the best combination of parameters in our resolution regime is given by the testH model (t_* = 20 Myr, as in the B4 model of the main paper).

A.4. The interplay between AGN, star formation, and shocks

With other three variations, we resimulated the same 21.5³Mpc³ volume alternatively by switching off AGN activity, star formation, or both – to assess the contribution by galaxy evolution processes to CRe and magnetic fields, as well as the interplay between feedback events and the additional generation of shocks. The reference galaxy formation model is again testA here ("baseline" in the figure). Figure A.6 shows the projected maps of the number of CRe only injected by shocks in the four relevant models (testA, test_noAGN, test_noSF and test_noSFnoAGN) at z = 0. The general distribution of CRe is very similar in all cases, with a few notable exceptions related to the inclusion or to the absence of AGN feedback, marked by the red circles: AGNs release powerful shocks around their host halos, which both inject new CRs via diffusive shock acceleration, as well as spread already existing CRe from previous shocks on a larger volume. At low redshift, the stellar feedback is in general too weak in our simulations to produce similar effects. Figure A.3 give the line profile for the projected CRe number distribution in four cases (for the green lines in the previous figure), showing the subtle modulation of the by AGN, which typically smoothes the distribution of CRe outside from the centre of halos.

We conclude that shocks driven by structure formation are, by far, the main contributor of CRe from the DSA mechanism. However, the feedback from AGNs can, beside directly injecting new CRe through their jets/outflows, also drive powerful shocks in the intergalactic medium, which further add CRe in coccoon-like structures surrounding halos.

A.5. The contribution of different astrophysical sources to magnetic fields

By comparing the same above four models (testA, test_noAGN, test_noSF and test_noSFnoAGN) we can also better visualise the role of each mechanism to the magnetisation of the cosmic web. Figure A.5 shows the mean (mass-weighted) magnetic field strength along the LoS for these models at three different epochs (z = 3, = 1 and = 0). Above z = 3, the magnetisation in and around halos is dominated by star formation-driven outflows. If both AGN activity and star formation are off, only the compression and (sub-grid) amplification of magnetic fields from structure formation are active. In a second stage (z ∼ 1) AGN feedback takes over as it generates large outflows of baryons, which carry magnetic fields with them. Finally, in the present-day Universe (z = 0), the cosmic web is magnetised at all scales, but regardless of the magnetisation from AGN-driven outflows, most of the cosmic volume remains dominated by its initial primordial magnetic field. The typical magnetic field strength at the virial radius of our most massive halos here (∼5 ⋅ 10¹³ M_⊙) is ∼0.1μG if AGN feedback is activated (with or without stellar feedback), ∼0.01μG if only stellar feedback is active, and ∼0.001μG if no feedback sources are active.

In filaments, on the other hand, the magnetic field amplitude by the end of our runs is always similar to the test_noSFnoAGN run, meaning that its value is mostly anchored to the evolution of primordial magnetic fields and that no contamination from astrophysical fields can outshine it.

Fig. A.2. History of the simulated cosmic star formation density (i.e. normalised for a Mpc3 comoving volume) for several resimulations of a 21.53 Mpc3 volume, in which we varied the parameters of our star formation feedback implementation (top panel), or of our prescriptions for AGN feedback (bottom panel). See Appendix A.2 for explanations.

Appendix B: Tests on the spectral modelling of CRe

We tested the validity of our procedure to model the ageing of the spectra of injected CRe, which combines our treatment of the time elapsed since the injection of CRe in cells (see Sect. 2.3) with the use of template CRe momentum spectra, depending on the local thermo-dynamical values in cells (see Sect. 2.7), against simpler approaches with a fixed spectrum for CRe. The top panel of Fig. A.6 gives the radio spectral index of the synchrotron radio emission from shock accelerated CRe, measured between 1400 and 150 MHz for our test volume at z = 0.

The bottom panel shows instead the ratio between the radio power at 150 MHz obtained with a simpler constant spectral index for all cells (i.e. α = 2 in momentum spectrum, yielding a I(ν)∝ν⁻¹ radio emission spectrum in DSA), and the power computed with our procedure to evolve CRe spectra as in the main paper.

Fig. B.1. Maps of the total number of CRe (injected by shocks only, in units of [log10(NCRe)]) along the LoS of a 10 × 10Mpc2 zoomed region in our test runs, where we included all sources of feedback in the simulation (’baseline’ model, as in the main paper) or we switched off AGN feedback (’no AGN’) or stellar feedback (’no SF’), or both (’no SF no AGN’). The red circles identify regions in which obvious differences are visible between models, while the green lines show the direction for which we generated radial profiles of the number of CRe in Fig. A.4.

If we adjusted a fixed spectrum to α = 2, the simplified calculation can reasonably reproduces the radio emission from recent merger shocks internal to halos and filaments, but it would dramatically exceed (by orders of magnitude) the radio emission in less dense environments, dominated by older populations of CRe. As is obvious from the distribution of spectral indices, no single value α can reproduce the realistic range of spectral indices of CRe across the cosmic web. We remark that assuming a fixed spectrum is not the standard procedure used in many works in the literature (e.g. Skillman et al. 2011; Vazza et al. 2015; Nuza et al. 2017). However, even in such work only the prompt emission from shock injected CRe can be estimated. Our simulations instead show that different populations of can dominate the budget of CRe in different regions of the cosmic web and, hence, that a realistic spectral modelling of the emission (such as the one proposed here) is necessary.

Fig. B.2. Line profile of the total number of CRe (injected by shocks only) for the four variations and the same line in Fig. A.3.

Fig. B.3. Projected mass-weighted magnetic field strength for runs test_noSFnoAGN, test_noAGN, test_noSF and testA, at three different epochs:z = 3 (top), z = 1 (middle) and z = 0 (bottom).

Fig. B.4. Top panel: Distribution of radio spectral index between 1400 and 150 MHz for CRe injected by shocks, for our fiducial post-processing treatment of CRe spectra, for a 21.253 Mpc3 volume at z = 0. Bottom panel: Map of the ratio between the synchrotron radio power at 150 MHz For the same volume and epoch, computed assuming a fixed α = −2 for the momentum spectrum of CRe injected by shocks, and the power from CRe assuming their realistic spectrum as in Sect. 2.7.

All Tables

All Figures

Fig. 1. Projected mean magnetic field intensity (top panels, in units of comoving Gauss) and projected mean mass-weighted gas temperature (bottom panels) across the full simulated 42 Mpc3 volume in our B4 run at z ≈ 2 (left), for the z ≈ 1 (centre) and z ≈ 0 (right) epochs. Movies of some of our runs can be found online and at https://www.youtube.com/playlist?list=PL8ecsjnxOKP7LXQICrNBLHZkPgxh0jCzj

In the text

	Fig. 2. Same as Fig. 1, but for a zoomed portion with side 15 × 15 Mpc2 (and 42 Mpc along the LoS) for our B4 run at three different epochs.
In the text

	Fig. 3. Top: Projected mean (mass weighted) gas temperature (in units of [K]) across the entire simulated volume for our seven runs at z = 0.02. Bottom: Projected mean (mass weighted) magnetic field strength (in units of [G]) for the same epoch and volume selection. The selected volume in all panels is about 21.25 × 42.5 Mpc2 and the length along the LoS is 42.5 Mpc.
In the text

	Fig. 4. Cooling time for relativistic electrons underoging synchrotron and inverse compton losses (as in Eq. (2)), as function of their Lorentz factor and cosmic redshift, for two B = 0.1 μG and B = 10 μG. The solid black lines mark the region corresponding to our choice of the timescale τ for the age determination of CRe injected at run-time.
In the text

	Fig. 5. Projected mean (CRe-weighted, only considering CRe from shocks) age of CRe since their last injection by shocks (in units of [Gyr]) across the entire simulated volume for our seven runs at z = 0.02.
In the text

	Fig. 6. Projected mean (CRe-weighted, only considering CRe from AGNs) age of CRe since their last injection by AGNs (in units of [Gyr]) across the entire simulated volume for our seven runs at z = 0.02.
In the text

	Fig. 7. Projected mean (CRe-weighted, only considering CRe from star formation) age of CRe since their last injection by star formation (in units of [Gyr]) across the entire simulated volume for our seven runs at z = 0.02.
In the text

	Fig. 8. Trend of injection efficiency of CRe (given in terms of the number of injected CRe with respect to the number of thermal particles) as a function of shock Mach number, for three different choices of the injection momentum, xinj (the baseline model assumed in our simulations is in red).
In the text

Fig. 9. Projected number of CRe across the full simulated 42.5 Mpc3 volume in our B4 run (considering only 5 Mpc along the LoS to better isolate single structures), for the z ≈ 2, z ≈ 1 and z ≈ 0 epochs (from top to bottom). The first column shows CRe injected by shocks, the second by AGN feedback and the third by star formation. The full movie of these fields can be found online and at https://youtu.be/xMTSK_rxZZU?si=hj4Yckl3hBuhYF-v

In the text

Fig. 10. Example of six different CR electrons spectra computed by our post-processing model, for different (random) variations of the combinations between gas density, magnetic field strength, redshift, and time elapsed since the injection of the initial power-law distribution of CRs. Such spectra are extracted from our tabulated list of spectra generated with the ROGER code and applied to compute the synchrotron emission from every cell in our simulated volumes in post-processing. Finally, we note that our approach neglects the spatial diffusion of cosmic rays, which is mostly negligible on the scales of interest here (this is addressed in more detail in Sect. 4).

In the text

	Fig. 11. History of the simulated cosmic star formation density (i.e. normalised for Mpc3 comoving) in our suite of simulations, as a function of cosmic time. The grey points with error bars show the observed cosmic star formation derived in Madau & Dickinson (2014).
In the text

Fig. 12. Stellar mass luminosity function for our simulated galaxies in all models and for three different epochs, compared with the best-fit derived by McLeod et al. (2021) for observations at approximately equal redshift bins. The last panel shows the distribution of the stellar mass fraction as a function of the host halo mass at z = 0.02, compared to the best-fit relation derived from observations by Behroozi et al. (2019).

In the text

Fig. 13. Simulated scaling relations for the SMBHs modelled at run-time in our A3 and B4 runs, at three different epochs. The first column show the relation between the gas mass of the host galaxy and the SMBH mass, compared with the scaling relation inferred from observations by Gaspari et al. (2019). The second column shows the relation between the SMBH mass and the mass accretion rate onto the SMBH; the third column gives the relation between the SMBH mass and the jet velocity associated with feedback events.

In the text

	Fig. 14. Evolution with cosmic time of the matter accretion rate onto star-forming particles (red lines) and of the matter accretion rate onto AGNs (blue lines) for our A3 and B4 models.
In the text

	Fig. 15. Projected mass weighted gas temperature (colours, in units of log10[K]) and over imposed contours of detectable synchrotron radio emission at 150 MHz from CRe injected by AGNs in a subselection of our volume (about 10 Mpc across) in six runs at z ≈ 1.0.
In the text

Fig. 16. Top row: Evolving distribution function of radio power at 150 MHz from our simulated galaxies, in which only the radio emission from CRe injected by star formation is considered. Bottom row: Evolving distribution function of radio power at 150 MHz from our simulated galaxies, in which only the radio emission from CRe injected by AGNs is considered. The dotted grey lines in the top rows give the best-fit relation of observed data with LOFAR, in the same redshift bins, from Cochrane et al. (2023) (top) and Kondapally et al. (2022) (bottom).

In the text

	Fig. 17. Relation between the emitted radio power at 150 MHz and the stellar mass (top panel) or the specific star formation rate (bottom panel) for our A3 and B4 run at z = 1.065.
In the text

	Fig. 18. Volume-filling factors of comoving magnetic field strength in our runs at three different epochs. The solid lines give the cumulative distribution for all cells in the simulation, while the dashed lines show the contribution from magnetic fields only injected by ‘astrophysical’ processes (i.e. AGN feedback and star formation).
In the text

Fig. 19. Top panel: Average ratio between CRe density and thermal gas proton density as a function of the gas over-density in the simulation by the end of our runs. Middle panel: Cumulative absolute number of CRe normalised for a comoving Mpc3 volume and averaging over an increasing gas overdensity. Bottom panel: Average time elapsed since the last injection of CRe as a function of gas overdensity. The solid lines give the trend for CRe injected by shocks, the dotted lines for CRe injected by AGNs and the dashed lines for CRe injected by star formation.

In the text

	Fig. 20. Distribution of the ratio between CRe density and baryon density within R500 of our halos at three different epochs and for all runs. The solid lines show the distribution of CRe injected by shocks, the dotted lines by AGNs and the dashed lines by star formation.
In the text

	Fig. 21. Projected RM maps (in units of [rad/m2]) across the entire simulated volume for our seven runs at z = 0.02. Each panel is approximately 21.25 × 42.5 Mpc2 large.
In the text

Fig. 22. Trend of simulated rms of the residual Faraday RM, observed at z = 0 and for integrating out to increasingly larger redshift (up to z = 2.5), for our A3, B4, and C1 runs. The lines give the mean dispersion of the RRMf for 100 randomly selected lines of sight through the simulated volume, while the black line with larger error bars shows the observed RRM trend by LOFAR, taken from Carretti et al. (2024).

In the text

	Fig. 23. Projected mass weighted gas temperature (colours, in units of log10[K]) and over imposed contours of detectable synchrotron radio emission at 150 MHz from CRe accelerated by shocks in a subselection of our volume (about 20 Mpc across) in our nine runs at z = 0.02.
In the text

	Fig. 24. Distribution of the ratio between shock-injected CRe density and baryon density (bottom row) within R500 of our halos for the C1 and C2 runs at the end of the simulation.
In the text

Fig. A.1. Projected mean (CRe-weighted) age of CRe since their last injection (in units of [Gyr]) by shocks (top row), AGNs (middle row), or star formation (bottom row) for two control resimulations at z = 0.02, which employed τ = 0.01 Gyr (model test_001, left panels) or τ = 1 Gyr (model test_time1, right panels) to solve Eq.1. Each panel is 21.25 Mpc across. No significant differences in the reconstructed age of the different CRe population are found.

In the text

	Fig. A.2. History of the simulated cosmic star formation density (i.e. normalised for a Mpc3 comoving volume) for several resimulations of a 21.53 Mpc3 volume, in which we varied the parameters of our star formation feedback implementation (top panel), or of our prescriptions for AGN feedback (bottom panel). See Appendix A.2 for explanations.
In the text

Fig. B.1. Maps of the total number of CRe (injected by shocks only, in units of [log10(NCRe)]) along the LoS of a 10 × 10Mpc2 zoomed region in our test runs, where we included all sources of feedback in the simulation (’baseline’ model, as in the main paper) or we switched off AGN feedback (’no AGN’) or stellar feedback (’no SF’), or both (’no SF no AGN’). The red circles identify regions in which obvious differences are visible between models, while the green lines show the direction for which we generated radial profiles of the number of CRe in Fig. A.4.

In the text

	Fig. B.2. Line profile of the total number of CRe (injected by shocks only) for the four variations and the same line in Fig. A.3.
In the text

	Fig. B.3. Projected mass-weighted magnetic field strength for runs test_noSFnoAGN, test_noAGN, test_noSF and testA, at three different epochs:z = 3 (top), z = 1 (middle) and z = 0 (bottom).
In the text

Fig. B.4. Top panel: Distribution of radio spectral index between 1400 and 150 MHz for CRe injected by shocks, for our fiducial post-processing treatment of CRe spectra, for a 21.253 Mpc3 volume at z = 0. Bottom panel: Map of the ratio between the synchrotron radio power at 150 MHz For the same volume and epoch, computed assuming a fixed α = −2 for the momentum spectrum of CRe injected by shocks, and the power from CRe assuming their realistic spectrum as in Sect. 2.7.

In the text