© The Authors 2024

1. Introduction

From observations of galaxy clusters, we can assume an efficient acceleration of cosmic ray (CR) electrons (CRes), which provide powerful tracers of the intra-cluster medium’s magnetic field such as radio relics and radio halos (see van Weeren et al. 2019, for a review). The recent discovery of “radio mega-halos” (Cuciti et al. 2022) and “radio bridges” (e.g., Govoni et al. 2019; Botteon et al. 2020a; Bonafede et al. 2021; Venturi et al. 2022; Radiconi et al. 2022) also indicates that there must be a substantial population of relativistic electrons even further from the cluster center than previously expected. The goal of this work is to explore if this is also true in cosmic web filaments.

The existence of a volume-filling, relativistic electron population poses a number of theoretical problems for the acceleration of electrons and their stabilization against energy losses, as well as magnetic field strength in these regimes. The canonical process of CR acceleration in galaxy clusters is diffusive shock acceleration (DSA) (see Drury 1983, for a review), where particles are accelerated from the thermal pool by scattering off magnetohydrodynamic (MHD) turbulence up and downstream of shock fronts, gaining energy at every crossing from upstream to downstream. This naturally produces a power-law distribution of CRs in momentum space, consistent with the power-law slope of the synchrotron spectrum of radio relics. Hybrid- and particle-in-cell (PIC) simulations of CR acceleration find that the acceleration efficiency at the shock depends strongly on the sonic Mach number (see e.g., studies by Kang & Ryu 2013; Caprioli & Spitkovsky 2014; Ryu et al. 2019), which would indicate a highly efficient acceleration of CRs at the periphery of the cluster; for example, at accretion shocks around galaxy clusters and cosmic filaments where high Mach number shocks are found in simulations (e.g., Ryu et al. 2003; Pfrommer et al. 2006; Vazza et al. 2009, 2011; Schaal & Springel 2015; Banfi et al. 2020; Ha et al. 2023). However, detailed studies on the efficiency of CRe acceleration in high-β (with $\beta =\frac{P_{\mathrm{th}}}{P_{\mathrm{B}}}\sim 50-{10}^3$) plasmas are typically performed for high-temperature, low-Mach-number shocks (e.g., Guo et al. 2014; Kang et al. 2019; Kobzar et al. 2021; Ha et al. 2021) to understand the discrepancy between the predicted electron acceleration efficiency and radio brightness of radio relics (see Botteon et al. 2020b, for a detailed discussion). Ha et al. (2023) recently performed PIC simulations with a high-β, low-temperature (T = 10⁴ K), and high-sonic-Mach-number (ℳ_s = 25 − 100) setup to study the acceleration efficiency of CRe in accretion shocks. Applying their model to post-process outputs from a cosmological simulation, they find emission of the order S_ν ∼ 1 − 10 mJy beam⁻¹ for accretion shocks around a Coma-like cluster, which is in agreement with emission by an accretion shock reported by Bonafede et al. (2022).

The second ingredient to potential synchrotron emission in filaments is the magnetic field strength in these environments. These field strengths can be estimated from rotation measurements, where the intrinsic polarization angle of a strongly polarized source, like synchrotron emission by electrons accelerated in an active galactic nucleus (AGN) jet, is rotated by the magnetic field along the line of sight to the observer. These observations indicate magnetic field strengths of 10–100 nG (e.g., Vernstrom et al. 2019; Carretti et al. 2023; O’Sullivan et al. 2019; O’Sullivan et al. 2020; O’Sullivan et al. 2023). Alternatively, the stacking of cluster pairs and filaments can provide upper limits on the synchrotron emission and with that the superposition of CRe population and magnetic field strength (e.g., Brown et al. 2017; Vernstrom et al. 2017, 2021; Locatelli et al. 2021; Hoang et al. 2023). Direct observations of radio emission from cosmic web filaments could help to disentangle this superposition, and therefore give an insight into the CRe acceleration, as well as the magnetic field strength and structure in cosmic web filaments. Several numerical studies have made predictions about such observations by taking into account different models for magnetic fields and CRe components and applying them to large-scale simulations.

Vazza et al. (2015a) used MHD simulations of a (50 Mpc)³ box (Vazza et al. 2014) and applied shock detection and the CR acceleration model of Hoeft & Brüggen (2007) to study potential synchrotron emission by shock-accelerated electrons in accretion shocks and internal shocks in cosmic web filaments. They varied their magnetic field strength by rescaling their simulated magnetic field with a density model and found in the high amplification model a diffuse emission from filaments on the order of ∼μJy beam⁻¹.

Brown et al. (2017) used a different approach by modeling CR electrons as secondaries of shock-accelerated protons. These protons are expected to live for longer than one Hubble time in the low-density cosmic web filaments and can scatter with thermal protons into charged pions, which decay into electrons and positrons (e.g., Blasi et al. 2007). These secondary models require a substantial primary proton population to produce considerable secondary electrons, which is in conflict with the current non-detection of γ-ray photons that should also originate from this process (see e.g., Wittor et al. 2019, and references therein). Brown et al. (2017) took the model by Dolag et al. (2004), which constructs radio power based on secondaries from a volume-filling proton population with a proton to thermal pressure ratio of X_cr = 0.01, and applied it to a constrained simulation resembling the local cosmic web (Dolag et al. 2004, 2005a). This allowed them to use S-PASS data to constrain the diffuse flux from synchrotron emission in cosmic web filaments to I_1.4 GHz < 0.073 μJy arcsec⁻², which is just below the current detection limit of, for example, the EMU survey (Norris et al. 2011).

Oei et al. (2022) used a (100 Mpc)³ cosmological MHD simulation by Vazza et al. (2019) to obtain synchrotron emission by accretion and merger shocks, again following the Hoeft & Brüggen (2007) model. They constructed a specific intensity function of the synchrotron cosmic web, by ray-tracing through their simulation box, while taking into account the cosmological distance effects of the emitted radiation. With this, they find specific intensities of the order I_νξ⁻¹ ∼ 0.1 Jy deg⁻² at ν = 150 MHz, with ξ being the poorly constrained efficiency parameter of electron acceleration.

Here, we take a different approach to a prediction of synchrotron emission from the cosmic web, by injecting and evolving a population of CRe at shocks and studying the magnetic field strengths required to obtain synchrotron emission compatible with current observational limits. In this work, we study potential emission from CRs accelerated at shocks in merger and accretion shocks of galaxy clusters and accretion shocks around cosmic web filaments. For this, we employed the first realization of a cosmological MHD simulation of a (500 h⁻¹ c Mpc)³ box with an on-the-fly Fokker-Planck solver to model CR proton and electron evolution. Here, we only focus on the electron populations, leaving the study of protons to future work. This paper is structured as follows. In Sect. 2, we describe the simulation and the code configuration employed for this work. In Sect. 3, we study full sky projections of the injected CR electrons and their synchrotron emission under different magnetic field models. In Sect. 4, we focus on the analysis of nonthermal emission from our replica of the Coma cluster and its surrounding filament structure. Section 5 extends this analysis to the full simulation domain to study the emissivity of all cosmic web filaments. In Sect. 6, we compare our results to the post-processing model of Hoeft & Brüggen (2007) and discuss the impact of free parameters in our CR model on our results and their observational prospects. Finally, Sect. 7 contains our conclusions and a summary of our work.

2. Methods

2.1. Initial conditions

The initial conditions of the simulations have been described extensively in Sorce (2018), Dolag et al. (2023), so we limit the description here to a short overview.

We set up constrained initial conditions based on velocity field reconstructions in the local Universe (Carlesi et al. 2016; Sorce 2018, and references therein). In this method, a Wiener filter algorithm (Zaroubi et al. 1995, 1998) is applied to reconstruct the 3D peculiar velocities of galaxies from the COSMICFLOWS-2 distance modulus survey (Tully et al. 2013). This survey undergoes several treatments upstream and downstream of this reconstruction (Doumler et al. 2013a,b,c; Sorce et al. 2014, 2017; Sorce 2015; Sorce & Tempel 2017, 2018). Subsequently, using the constrained realization algorithm (Hoffman & Ribak 1991), it can be used to construct a density field of the local Universe at an initial redshift; in other words, the constrained initial conditions. The resolution of the density field is then increased to reach our study requirement using the GINNUNGAGAP software¹.

These initial conditions have been used, for instance, to study individual clusters such as Virgo (Sorce et al. 2016, 2019; Olchanski & Sorce 2018), as well as the large-scale environment around the Local Group (Carlesi et al. 2016, 2017) and other clusters (Sorce et al. 2024), in particular the Coma cluster (Malavasi et al. 2023). Previous work on ultrahigh-energy CRs has shown how these constrained initial conditions help to provide predictions for observations (Hackstein et al. 2018).

In this work, we study the SLOW-CR3072³ simulation, which consists of 2 × 3072³ particles in a volume of V = (500 h⁻¹ c Mpc)³. This leads to a mass resolution of M_gas ≈ 8.5 × 10⁷ M_⊙ and M_DM ≈ 4.6 × 10⁸ M_⊙ for gas and DM particles, respectively. The maximum resolution of gas particles in our simulation corresponds to h_sml, min ≈ 6.24 kpc at z = 0.

To test the impact of resolution on our results, we also performed zoom-in simulations of the Coma cluster at eight times the resolution of the cosmological box. To create such zoom-in initial conditions, a dark-matter-only simulation has been extended far into the future, for example up to an expansion factor of a ≈ 1000, to identify the gravitationally bound region of the supercluster that the target cluster is a part of. Then, the spatial region in the simulation that contains all particles bound to the Coma supercluster was identified. This region was then traced back to the corresponding Lagrangian volume in the initial conditions, which was then sampled with the particles taken from much higher-resolution versions of the box; in this case, the one sampled with 6144³ particles in total. In order to smooth the resolution difference at the boundaries while minimizing the total particle number, boundary layers were added to the high-resolution region, using the lower-resolution versions of the full box in descending order of particle number (3072³, 1536³, 768³, and 384³ for the remainder of the box). The refinement of the initial conditions leads to a mass resolution for the high-resolution region of M_gas ≈ 1 × 10⁷ M_⊙ and M_DM ≈ 5.7 × 10⁷ M_⊙ for gas and DM particles, respectively. Here, we reached a maximum resolution of gas particles of h_sml, min ≈ 3.7 kpc. We list all simulations discussed in this work in Table 1.

All simulations were run using a Planck cosmology (Planck Collaboration XVI 2014) with matter densities of Ω_m = 0.307 and Ω_baryon = 0.048, a cosmological constant of Ω_Λ = 0.692, and the Hubble parameter H₀ = 67.77 km s⁻¹ Mpc⁻¹.

2.2. Simulation code

We employed OPENGADGET3 (Groth et al. 2023), an advanced version of the cosmological TREE-SPH code GADGET2 (Springel 2005). OPENGADGET3 uses a Barnes and Hut tree (Barnes & Hut 1986) to solve the gravity at short range and a particle-mesh (PM) grid for long-range forces. For the zoom-in simulations used in this work, the PM grid was split to add a higher-resolution grid to the central region to save computing time.

All simulations presented in this work are non-radiative CR-MHD simulations. OPENGADGET3 uses an updated smoothed particle hydrodynamics (SPH) implementation (Beck et al. 2016a) with a spatially and time-dependent high-resolution shock capturing scheme (Dolag et al. 2005b; Cullen & Dehnen 2010). The simulations presented here were run with a Wendland C₄ kernel with 200 neighbors and bias correction (as introduced in Dehnen & Aly 2012) to ensure stabilization against the tensile (pairing) instability that arises for kernels that do not have a positive definite Fourier transformation. Our SPH implementation was further stabilized by a high-resolution maximum entropy scheme, which in this work was realized by physical rather than artificial thermal conduction.

The MHD solver was presented in Dolag et al. (2009) and extended in Bonafede et al. (2011) to include nonideal MHD in the form of magnetic diffusion and a subgrid model for magnetic dissipation in the form of magnetic reconnection, which can heat the thermal gas. To enforce the ∇ ⋅ B = 0 constraint, we adapted the hyperbolic constrained divergence cleaning of Tricco et al. (2016) (for the implementation into OpenGadget3 see Steinwandel & Price, in prep.).

We employed an on-the-fly shock finder (Beck et al. 2016b) to compute the necessary shock quantities for CR acceleration. Cosmic rays were represented as populations of protons and electrons and evolved in time using an on-the-fly Fokker-Planck solver, which we discuss in more detail in the next section.

2.3. Cosmic ray model

To model CRs, we employed the on-the-fly Fokker-Planck solver CRESCENDO introduced in Böss et al. (2023a) and will only briefly outline the solver and the parameters used here in comparison to the previous work. CRESCENDO attaches populations of CR protons and electrons to all resolution elements of the simulation. These populations are assumed to be isotropic in momentum space and are represented as piece-wise power laws,

$$\begin{array}{c}\hfill f(p)=f_i{\left(\frac{p}{p_i}\right)}^{-q_i},\end{array}$$(1)

for a number of logarithmically spaced bins. Protons were considered in the dimensionless momentum range $\widehat{p}\equiv \frac{p}{m_{p}c}\in [0.1,{10}^5]$ and electrons in the range $\widehat{p}\equiv \frac{p}{m_{e}c}\in [1,{10}^5]$, where m is the mass of the individual particle species. We discretized the populations with 6 bins (1/dex) for protons and 20 bins (4/dex) for electrons. This choice of resolution was mainly driven by memory constraints, with the focus on still being able to model accurate high-momentum losses for electrons. These populations were then evolved in time by solving the diffusion-advection equation in the two-moment approach, following Miniati (2001) (see also Girichidis et al. 2020; Ogrodnik et al. 2021; Hopkins et al. 2022, for other recent on-the-fly spectral CR models).

We accounted for injection at shocks via DSA, adiabatic changes due to density changes in the surrounding thermal gas, and energy losses of electrons due to synchrotron emission and inverse Compton (IC) scattering off CMB photons. This reduced the Fokker-Planck equation, which describes the time evolution of our CR populations to

$$\begin{array}{cc}\hfill \frac{Df(p,\mathbf{x},t)}{\mathit{Dt}}& =(\frac{1}{3}\mathrm{\nabla }·\mathbf{u})p\frac{\partial f(p,\mathbf{x},t)}{\partial p}\hfill \end{array}$$(2)

$$\begin{array}{cc}& +\frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2{\displaystyle \sum _l}{b}_{l}f(p,\mathbf{x},t)\right)\hfill \end{array}$$(3)

$$\begin{array}{cc}& +j(\mathbf{x},p,t),\hfill \end{array}$$(4)

where Eq. (2) denotes changes due to adiabatic expansion or compression of surrounding gas, Eq. (3) denotes radiative changes due to synchrotron emission and IC scattering off CMB photons with ${\displaystyle \sum _l}{b}_l\equiv {\left|\frac{\mathrm{d}p}{\mathrm{d}t}\right|}_{\mathrm{syn}}+{\left|\frac{\mathrm{d}p}{\mathrm{d}t}\right|}_{\mathrm{IC}}$, and Eq. (4) denotes the source term of CRs. In this work, the source term j(x, p, t) describes CR acceleration at shocks, modeled following the DSA parametrization by Ryu et al. (2019). In addition to that, we used the obliquity-dependent model introduced by Pais et al. (2018) for protons that are primarily accelerated at quasi-parallel shocks, and shifted it by 90° for electrons (see Appendix A for more details). We injected the energy into the full momentum range as a single power law following a linear DSA slope in momentum space,

$$\begin{array}{c}\hfill q=\frac{3r}{r-1};r=\frac{(\gamma +1){\mathcal{M}}_s^2}{(\gamma -1){\mathcal{M}}_s^2+2},\end{array}$$(5)

where r is the shock compression ratio, γ is the adiabatic index of the gas, and ℳ_s is the sonic Mach number. The DSA parametrization by Ryu et al. (2019) imposes a critical Mach number, ℳ_s = 2.25, below which no efficient CR acceleration can take place. We employed a fixed electron-to-proton energy injection ratio of K_ep = 0.01.

For this run, we used closed boundary conditions at the lower end of the distribution function, as they provide more numerical stability and mimic low-momentum cooling on adiabatic compression, which is not explicitly included in this simulation. We solved Eq. (2) in the two-moment approach by accounting for CR number and energy changes per bin for both protons and electrons for reasons of numerical stability (see Girichidis et al. 2020, for a discussion of the benefit for protons).

Feedback from the CR component to the thermal gas in the form of comoving pressure was calculated by solving the integral

$$\begin{array}{c}\hfill P_{\mathrm{CR},c}=\frac{4\pi }{3}{a}^4{\displaystyle \underset{p_{\min }}{\overset{p_{\mathrm{cut}}}{\int }}}\mathrm{d}p{p}^{2}T(p)f(p)\end{array}$$(6)

under the approximation T(p)≈pc. We started the injection of CRs at z = 4. All these effects were computed on the fly for every gas particle active in the timestep. Constraints on computing memory made it impossible to also account for transport processes, besides advection, which we accounted for naturally, since the CR populations are attached to our Lagrangian SPH particles. We shall study the effect of transport processes on CRs in galaxy clusters in future zoom-in simulations instead.

2.4. Fermi-II reacceleration

We implemented an on-the-fly treatment of turbulent reacceleration of CRs and briefly outline the implementation here (for details, see Appendix C). To follow the turbulent reacceleration of CRs, we adopted the model by Cassano & Brunetti (2005). We used the unified-cooling approach, which extends the b_l(p) term in Eq. (3) to include the systematic component of the reacceleration,

$$\begin{array}{c}\hfill {\displaystyle \sum _l}{b}_l(p)\equiv {\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{cool}}^{\mathrm{synch}}+{\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{cool}}^{\mathrm{IC}}+{\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{acc}}^{\mathrm{sys}}.\end{array}$$(7)

Following Cassano & Brunetti (2005), the change in momentum due to turbulent reacceleration per timestep can be expressed as

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{acc}}^{\mathrm{sys}}=-\chi p\approx -2\frac{D_{\mathrm{pp}}}{p}=-2D_0(t)p,\end{array}$$(8)

where we used D_pp = D₀(t) p².

The computation of D₀(t) given in Eq. (C.10) is quite expensive and the way it is coupled to Eq. (2) puts a significant timestep constraint on the solver. We therefore sub-cycled the solver if D₀(t)×Δt surpassed a critical value. Additionally, we introduced a cap for D₀(t) to avoid overly strong impacts of individual particles on the total runtime of the simulation.

Due to the additional computational cost, we switched this mode off in SLOW-CR3072³ and ran comparison zoom-in simulations with the mode switched on. For COMA, we switched the mode off to obtain a comparison and see the impact of resolution. In COMA-D_pp, we ran the turbulent reacceleration on the fly, but capped the value of D₀, which limits the need to sub-cycle the solver. We capped the value at D_0, max = 10⁻¹⁸ s⁻¹, which lies at the lower end of the values reported in Appendix B in Donnert & Brunetti (2014) but is close to the maximum values found in the Coma filaments of our simulation (we discuss this in Sect. 4). However, this model will need to be revised to include recent results that indicate that preferentially the solenoidal component of turbulence can efficiently reaccelerate CRs (Brunetti & Lazarian 2016; Brunetti & Vazza 2020), or allow for other reacceleration mechanisms (e.g., Tran et al. 2023).

2.5. Synchrotron emission

As in the introductory paper for CRESCENDO, we computed the synchrotron emissivity, j_ν, in units of [erg s⁻¹ Hz⁻¹ cm⁻³] as

$$\begin{array}{c}\hfill j_{\nu }(t)=\frac{\sqrt{3}{e}^3}{m_e^2{c}^3}B(t){\displaystyle \sum _{i=0}^{N_{\mathrm{bins}}}}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}}\mathrm{d}\theta {sin}^2\theta {\displaystyle \underset{{\widehat{p}}_{\mathrm{i}}}{\overset{{\widehat{p}}_{\mathrm{i}+1}}{\int }}}\mathrm{d}\widehat{p}4\pi {\widehat{p}}^{2}f(\widehat{p},t)K(x),\end{array}$$(9)

where e is the elementary charge of an electron, m_e its mass, c the speed of light, $\widehat{p}$ the dimensionless momentum, and K(x) the first synchrotron function,

$$\begin{array}{c}\hfill K(x)=x{\displaystyle {\int }_x^{\infty }}\mathrm{d}z{K}_{5/3}(z),\end{array}$$(10)

using the Bessel function, K_5/3, at a ratio between the observation frequency, ν, and critical frequency, ν_c,

$$\begin{array}{c}\hfill x\equiv \frac{\nu }{{\nu }_c}=\frac{\nu }{C_{\mathrm{crit}}B(t)sin\theta {\widehat{p}}^2};C_{\mathrm{crit}}=\frac{3e}{4\pi m_{e}c}·\end{array}$$(11)

All scalings of the synchrotron emission with frequency and magnetic field strength arise directly from the evolved CR electron distribution function, $f(\widehat{p},t)$, in Eq. (9), without further assumptions or imposed limits (see Appendix A in Böss et al. 2023b, for a more detailed description).

3. Full-sky projections

Having a full box simulation of a constrained local Universe puts us in the unique position to study full-sky projections observed from our Milky Way replica. To obtain the projections, we mapped the SPH particles of our simulation onto a HEALP IX sphere (Górski et al. 2005) following the algorithm described in Dolag et al. (2005c). All figures in this section were obtained from the output at redshift z = 0.

3.1. Cosmic ray electrons

We show the full-sky projec|ion of our CR electron component in Fig. 1. The labeled clusters are our best matches for the observed counterparts (for details on the matching process, see Hernández-Martínez et al. 2024) and the circles indicate their projected virial radii.

Fig. 1. Full-sky projections of the CR electron component. Upper panel: Energy density of electrons in the range E ∈ [5 − 500] MeV as the mean along the line of sight in a radius range r = 5 − 300 Mpc. This traces the injected energy density of the CR electrons with the longest lifetimes. Circles indicate the projected rvir of the cross-identified clusters (see Hernández-Martínez et al. 2024). Lower panel: CR electrons with energies ECR, e > 1 GeV integrated along the line of sight in a radius range r = 5 − 300 Mpc. The short cooling times of these electrons lead to this map tracing recent injection events.

We show the two opposite ends of the CR electron energy range. The upper plot shows the mean CR electron energy density along the line of sight in a radius range, r = 5 − 300 Mpc, for CRs if the energy range E ∈ [5 − 500] MeV. This traces the electrons with the longest lifetimes, with a bias toward electrons that have been accelerated and advected into high-density regions. These electrons provide potential seed populations for reacceleration by subsequent shocks and/or turbulence.

We find a persistent population of CR electrons in galaxy clusters and the cosmic web. However, this can only be considered to be an upper limit, especially in the highest-density regions, as we do not explicitly treat low-momentum energy losses for CRs. Hence, this picture should be viewed more as a tracer for shock injection and adiabatic compression and with that a tracer of potential CR electron seed populations. A distinct feature and advantage of an on-the-fly approach for CR evolution is that we find large-scale halos of electron energy density around clusters and filaments. This is consistent with recent observations of radio mega-halos (Cuciti et al. 2022), indicating a significant population of CR electrons in the cluster volume. It has been shown by Beduzzi et al. (2023) that given a turbulent velocity of v_turb > ∼150 km s⁻¹ (which is only a factor of 1.5 larger than, for example, the measurements of the atmosphere in Perseus by Hitomi Collaboration 2018) these electrons have almost unlimited lifetimes and with that can remain synchrotron-bright at low frequencies, providing the basis for these radio mega-halos. We shall discuss the implications of this in future work focusing on radio halos and radio relics in our simulation.

In the lower panel of Fig. 1, we show the energy contained in electrons with energies above 1 GeV. These electrons are potentially synchrotron-bright at frequencies of 144 MHz. Since the cooling times of these CRs are of the order ∼10² − 10³ Myr, this predominantly traces recent injection. We find significant injection at the accretion shocks of most clusters, as well as the accretion shocks around cosmic web filaments. This shows that accretion shocks around filaments, together with internal shocks, provide a significant, volume-filling source for high energy CRs, available to emit radio emission, given sufficiently strong magnetic fields.

3.2. Magnetic field

Since the synchrotron emissivity scales strongly with the magnetic field strength (L_ν ∝ B^{α₀ + 1}, where α₀ is the spectral index of the synchrotron spectrum), the observability of synchrotron emission of the cosmic web is tightly linked to the magnetic field strength in filaments.

The magnetic field in our simulation is set up as B = (10⁻¹⁴, 0, 0) G at the starting redshift z = 120. This field then evolves and is amplified by adiabatic compression and dynamo processes (see e.g., Stasyszyn et al. 2010; Pakmor et al. 2014; Marinacci et al. 2015; Vazza et al. 2014, 2018, 2021; Steinwandel et al. 2022, 2024). As this is a non-radiative simulation without star formation or cooling, we cannot account for astrophysical sources of a magnetic field that can contribute to volume-filling magnetic fields in low-density regions (Beck et al. 2013; Garaldi et al. 2021). While we obtain reasonable values in galaxy clusters (of the order ∼μG, see e.g., Bonafede et al. 2010; van Weeren et al. 2019, the latter for a review), the magnetic field strength in filaments remains negligible. We attribute this to a lack of resolution in the filaments to drive an efficient turbulent dynamo to amplify the magnetic field (see Steinwandel et al. 2022, for an extensive study performed with our simulation code).

To test the implications of different magnetic field strength values in filaments on their potential synchrotron emission, we adopted five magnetic field models, two based on pressure scaling and three based on density scaling. We show the mean magnetic field along the line of sight in a radius range, r = 10 − 300 Mpc, for the simulation and the scaling models in Fig. 2.

Fig. 2. Full-sky projections of the mean magnetic field along the line of sight for the simulated magnetic field and post-processing models. Upper left panel: Simulated magnetic field obtained from the simulation. Upper right panel: Modeled magnetic field, assuming a plasma beta of $\beta \equiv \frac{P_{\mathrm{th}}}{P_{\mathrm{B}}}=50$. Center left panel: Modeled magnetic field based on a fraction ℱ of the turbulent pressure. Center right panel: Modeled magnetic field based on flux freezing B ∝ ρ2/3. Lower left panel: Modeled magnetic field, assuming a saturated turbulent dynamo scaling B ∝ ρ1/2 down to ne = 10−4 cm−3 and steeply declining beyond this density to match expectations for void magnetic fields. Lower right panel: Modeled magnetic field, assuming a saturated turbulent dynamo scaling for all density regimes.

First, we tested a simple model where the magnetic field pressure is a constant factor of the thermal pressure P_th, such that

$$\begin{array}{c}\hfill B_{\beta }=\sqrt{\frac{8\pi P_{\mathrm{th}}}{\beta }},\end{array}$$(12)

where we have assumed a constant plasma-β of β = 50. This is an optimistic model, given that such low values of β are typically only found in the central regions of clusters, while realistic values for filaments and accretion shocks are of the order β ∼ 10² − 10³ (see e.g., Ha et al. 2023).

Second, we calculated the magnetic field as a fraction, ℱ, of the turbulent pressure,

$$\begin{array}{c}\hfill B_{\mathcal{F}}=\mathcal{F}\sqrt{4\pi \rho v_{\mathrm{turb}}^2}.\end{array}$$(13)

ℱ can be estimated from the magnetic Reynolds number following Table 1 in Schober et al. (2015) and takes typical values of ℱ ≈ 0.01 − 0.1 in the systems we are interested in. However, to better match the central values of the magnetic field strength in our Coma analog, we adapted ℱ = 1. This assumes pressure equilibrium between turbulent and magnetic pressure, as it is expected in a regime where the turbulent dynamo is saturated (see e.g., Steinwandel et al. 2024, and references therein). Recent results by Zhou et al. (2024) show that shocks are potential regions of efficient magnetic field generation and amplification and can reach this saturation value relatively quickly. We therefore used the model B_ℱ as a maximum assumption of their work.

For the magnetic field models scaling with the density, we first used a simple flux-freezing model where initial magnetic fields are amplified by adiabatic compression in collapsing halos. Following Hoeft et al. (2008), we adapted

$$\begin{array}{c}\hfill B_{\mathrm{ff}}=\frac{B}{0.1\mathrm{\mu }\mathrm{G}}{\left(\frac{n_e}{{10}^{-4}{\mathrm{cm}}^{-3}}\right)}^{2/3},\end{array}$$(14)

scaling this model up by a factor of four to better match the central magnetic field strength in our Coma analog. We shall discuss this in more detail in Sect. 4.

In ultrahigh-resolution simulations of massive galaxy clusters, Steinwandel et al. (2024) find that the turbulent dynamo is the dominant amplification process, even in low-density regions down to n_e ∼ 10⁻⁴. From this, we extrapolated two turbulent dynamo models. The first assumes that the turbulent dynamo amplifies the magnetic field down to densities of n_e ∼ 10⁻⁴ and below that we applied a fit function to intersect the observations by Carretti et al. (2023) with B ∼ 30 nG at filaments with densities of n_e = 10⁻⁵ cm⁻³ and extrapolated to B ∼ 10⁻¹⁴ G at n_e = 10⁻⁶ cm⁻³ as a proxy for void magnetic fields (see e.g., Neronov & Vovk 2010; Caprini & Gabici 2015, for the large differences in the estimation of this value).

This gives a model dependent on n_e as

$$\begin{array}{c}\hfill B_{\mathrm{dyn},\downarrow }(\rho )=\{\begin{array}{cc}2.5\times {10}^{-6}{\left(\frac{n_e}{{10}^{-3}{\mathrm{cm}}^{-3}}\right)}^{1/2}\hfill & \text{[G], if}{n}_e>{10}^{-4}{\text{cm}}^{-3}\hfill \\ {\displaystyle \sum _{i=1}^5}{p}_i{log}_{10}{(n_e)}^{i-1}\hfill & [{log}_{10}\text{G], if}{n}_e<{10}^{-4}{\text{cm}}^{-3}\hfill \end{array},\end{array}$$(15)

where p ∈ [ − 16.38, −16.0, −8.07, −1.71, −0.13].

We also considered an upper limit, as a pure saturated dynamo scaling with density as

$$\begin{array}{c}\hfill B_{\mathrm{dyn},\uparrow }(\rho )=2.5\times {10}^{-6}{\left(\frac{n_e}{{10}^{-3}{\mathrm{cm}}^{-3}}\right)}^{1/2}\text{[G]},\end{array}$$(16)

where we fit the results from Steinwandel et al. (2024) and extended the scaling to lower densities with B ∼ 250 nG at n_e, 0 = 10⁻⁵ cm⁻³. Identical to the previous scaling, this matches magnetic field strengths in cluster centers very well, but overestimates the magnetic field strength in filaments by roughly one order of magnitude, compared to the Carretti et al. (2023). This is, however, consistent with previous works by Vernstrom et al. (2017), Brown et al. (2017), O’Sullivan et al. (2019), Locatelli et al. (2021), who allow for up to 250 nG. To illustrate the scaling of our models as a function of electron density, we refer to Fig. D.1.

3.3. Synchrotron emission

By folding the lower panel of Fig. 1 with the different panels of Fig. 2 according to Eq. (9), we can obtain the synchrotron emissivity. From there, we arrive at the intrinsic synchrotron surface brightness by integrating the emissivity of the SPH particles along the line of sight. We show the result of this for the previously introduced magnetic field models in Fig. 3.

Fig. 3. (Intrinsic) synchrotron intensity at 144 MHz integrated along the line of sight between r = 10 − 300 Mpc. The different panels show the emission calculated based on the magnetic field configurations shown in Fig. 2.

In the case of the simulated magnetic field, we only obtain synchrotron emission from galaxy clusters. The emission in cluster centers is of the order ∼ten lower compared to observations, which is mainly driven by the flatter radio spectra. Our simulated spectra in the box stem purely from direct injection and advection, so they lack additional steepening from on-the-fly turbulent reacceleration. We discuss this in more detail for the filaments of Coma in Sect. 4 and for the radio emission in Coma, Virgo, and Perseus in Dolag et al. (in prep.).

The lack of emission in the filaments stems from the steeper dropoff of the magnetic field as a function of density. As was discussed before, the simulated magnetic field strength in the filaments is well below the nano-Gauss level and with that would require substantial CR injection to be synchrotron-bright.

Assuming a magnetic field following a constant plasma-β (B_β) yields more interesting results. We still primarily get synchrotron emission from the clusters; however, we also see more emission by smaller clusters that do not have a strong enough simulated magnetic field to be synchrotron-bright. In addition, we also see more relics far from the cluster center, such as in our Perseus replica to the very left of the map. This stems from the magnetic field scaling with the thermal pressure and with that automatically tracing regions of over-pressure; namely, shocks. We can observe some filamentary structure in the synchrotron emission; however, this is mainly driven by the halos inside the filaments, rather than diffuse emission from filaments.

The magnetic field scaling with turbulent velocity (B_ℱ) produces morphologically similar features as the constant plasma-β case. This naturally originates from shocks inducing a substantial amount of turbulent energy (e.g., Ryu et al. 2003; Vazza et al. 2011), so this model again amplifies magnetic field strength and, with that, synchrotron emission at shocks. This provides even more synchrotron emission from accretion shocks around clusters than in the B_β-model and the first indication of synchrotron emission from accretion shocks around cosmic web filaments.

The density-dependent magnetic field models B_ff and B_dyn,↓ both show similar behavior. These models lead to preferential emission in the center of halos with a steep decline toward lower densities. This leads to filamentary emission driven by halos within the filaments, but no diffuse emission from the gas in filaments or accretion shocks.

B_dyn,↑ is the only model that produces considerable diffuse synchrotron emission from filaments. We find consistent diffuse emission within the filaments and in some regions in the central and top right parts of the map, the magnetic field is even strong enough to illuminate part of the accretion shocks onto the clusters. The emission by filaments and accretion shocks have intensities more than five orders of magnitude lower than the intensities in cluster centers, making prospects of direct detection very difficult, as was expected.

Further, from a numerical standpoint, it is difficult to capture the spurious shocks in the low-density environment around filaments which can lead to an underestimation of the synchrotron emission in these maps. This mainly stems from the numerical resolution in these low-density regions. We discuss the impact of resolution on our findings in Sect. 4.

4. The Coma filaments

To study the potential synchrotron emission from cosmic web filaments in more detail, we shall now focus on the Coma cluster replica found in our simulation and its surroundings. The cosmic web filaments around our Coma replica have recently been studied in Malavasi et al. (2023), where they find good agreement with indications from observations, which provides a reasonable basis for our analysis.

4.1. Impact of resolution

To test the impact of resolution on our results, we added to our analysis by presenting results from zoom-in simulations of the Coma cluster (COMA and COMA-D_pp) at eight times the resolution of the cosmological box (SLOW-CR3072³). We show the comparison between SLOW-CR3072³ and COMA in the top and bottom panels of Fig. 4. From left to right, we show surface density, the absolute strength of the simulated magnetic field, the CR electron energy of the electrons with energies above 1 GeV, and the turbulent reacceleration coefficient given in Eq. (C.10).

Fig. 4. Cutout around our Coma replica with a side-length of d = 20 h−1 c Mpc. The circles represent r200 and r500. From left to right, we show surface density, absolute value of the simulated magnetic field, CR electron energy of electrons with energies above 1 GeV, and turbulent reacceleration coefficient. The upper panels show the cluster in the full cosmological box; the lower panels show it in the zoomed-in re-simulation.

We find excellent morphological agreement between the different setups, as can be seen in the density plot in the left panels, providing us with a good basis for comparison. In Coma itself, we can also see reasonable agreement in the magnetic field strength (second panels), injected CR energy (third panels), and reacceleration coefficient (last panels).

However, as we move away from Coma, we can observe the impact of resolution. In the magnetic field strength, this is evident as the magnetic field retains higher values further from the cluster center in the case of the zoom-in simulation. As was discussed previously, we attribute this to the lack of resolution to drive an efficient amplification via the turbulent dynamo in SLOW-CR3072³, as is shown in Steinwandel et al. (2022).

Another impact of resolution can be seen in the CR electron panels, where the increased resolution in the zoom-in simulation leads to a more accurate shock capturing and with that improved CR injection. This can be seen both in the strong shock of the ongoing merger, as well as in the surrounding medium. In the zoom-in simulation, even the equatorial shock from the current merger is captured, with the cluster in the full box only showing a hint of this shock in direct comparison.

Both simulations show an envelope of a previous shock moving away from the center and colliding with the accretion shock (see e.g., Zhang et al. 2020a,b, for a high-resolution study of this process). We find a prominent acceleration of CRs at a filament to the top right of the cluster in both simulations. This is, however, better refined in the zoom-in simulation, where a clear filament can be seen in the CR component, while the box simulation shows more CR injection on one side of the filament than the other.

4.2. Fermi-II reacceleration

The turbulent reacceleration coefficient maps in the rightmost panels of Fig. 4 again agree reasonably in morphology and absolute value. We find strong potential for reacceleration in the post-shock regions of the merger and the center of the cluster. Here values reach up to D₀ ∼ 10⁻¹⁶ s⁻¹, in agreement with the maximum reported in Donnert & Brunetti (2014). This shows how conservative a cap at D₀ = 10⁻¹⁸ s⁻¹ is in the case of a massive galaxy cluster. However the filaments only reach values of D₀ ∼ 3 × 10⁻¹⁸ s⁻¹. Since we focus on filaments in this work we accept this limitation for the benefit of substantially lower run-time of the simulation.

4.3. Magnetic field

In Fig. 5, we show RM mock-observations and radial profiles of the magnetic field strength for our Coma replica. The top panels show the RM mock observations where we rescaled the simulated magnetic field to match our magnetic field models in Sect. 3.2. the bottom panels show the radial profiles of the total magnetic field strength. Different colors correspond to the three simulations, and dashed black lines show the magnetic profile derived from the observations of Coma presented in Bonafede et al. (2010).

Fig. 5. Magnetic field models applied to the COMA simulation. Top: Rotation measure (RM) mock observations of our Coma replica assuming different magnetic field models. Bottom: Radial profiles of the magnetic field. We show the profiles of all discussed simulations for the various magnetic field models in the solid colored lines. The dashed black lines show the observed profile by Bonafede et al. (2010).

For the simulated magnetic fields, all simulations are in reasonable agreement with the observed profile up to a radius of ∼1 Mpc. Beyond this radius, the magnetic field profile drops off too steeply, which we attribute to the lack of resolution in the cluster outskirts to resolve the turbulent dynamo amplification, as was mentioned above.

The magnetic field model based on a constant plasma-β (B_β) similarly agrees well with observations within the inner Mpc, after which the magnetic field is overestimated due to the unrealistically large energy density ratio of the magnetic field to the thermal pressure in the cluster periphery.

Our model based on the fixed ratio between turbulent and magnetic pressure (B_ℱ) overestimates the central magnetic field strength by roughly 50 percent, but agrees very well with the radial decline in the cluster outskirts. As mentioned above, this upscaled model is the maximum assumption with pressure equipartition between turbulent kinetic and magnetic pressure and should be considered more as an upper limit.

We applied a similar up-scaling to the B_ff model, which yields good matches between modeled and observed magnetic field strength.

Both turbulent dynamo scaling models (B_dyn,↑ and B_dyn,↓) agree well with the central value of the magnetic field strength, but differ in the outskirts of the cluster. B_dyn,↓ transitions to a steeper density scaling at roughly the virial radius and with that a steeper radial decline.

B_dyn,↑ retains its $\sqrt{\rho }$-scaling and with that diverges from the observed radial decline, leading to an overestimation of the magnetic field in regions with n_e < 10⁻⁴ cm⁻³.

4.4. Synchrotron flux

In Fig. 6, we show the expected synchrotron flux of our simulated Coma cluster and its environment, assuming a LOFAR beam size of θ = 60″ × 60″ (as in Fig. 2 of Bonafede et al. 2022). We convert the intrinsic radio power P_ν to an observed flux by placing our Coma replica at a redshift of z = 0.0231 and converting it via $S_{\nu }=\frac{P_{\nu }}{4\pi D_L}$ where we use our intrinsic Planck Collaboration XVI (2014) cosmology for the computation of the luminosity distance D_L. We find that our simulated relic is 1–2 orders of magnitude less luminous than the observed one. We attribute this first to the possible connection of the relic to a narrow-angle tail galax in the real Coma cluster (e.g., Bonafede et al. 2022; Churazov et al. 2022), which cannot be reproduced here since we do not include galaxy formation and AGN jet physics in this simulation. A second reason is the known discrepancy between current DSA models and the up to 10 percent efficiency required to reproduce radio relic brightness (see discussion in Botteon et al. 2020b). We discuss the impact of CR injection parameter choices in more detail in Sect. 6. For this reason, we rescale the simulated synchrotron flux by a factor of 100, to roughly match the observed brightness of the radio relic in the zoom-in simulations using the simulated magnetic field. As the acceleration efficiency of CR electrons is still under debate in the literature, this instead allows us to discuss the synchrotron flux in the filaments relative to the flux in the relics.

Fig. 6. Same cutout as in Fig. 4, but showing the synchrotron surface brightness at 144 MHz, assuming a LOFAR beam size of θ = 60″ × 60″ as in Bonafede et al. (2022). The color bar is split at 1 μJy beam−1, which is a factor two above the sensitivity limit in the stacking approach by Hoang et al. (2023). From top to bottom, we show the different simulation runs, and from left to right we vary the magnetic field models. The lowest row shows the result of painting on synchrotron emission with the Hoeft & Brüggen (2007) model in a post-processing approach. For this, we use the COMA simulation due to the better performance of the shock finder at this resolution, compared to SLOW-CR30723.

The top panels show the Coma replica in the cosmological box (SLOW-CR3072³), the second row shows the cluster in the zoom-in simulation without turbulent reacceleration (COMA) and in the third row we include turbulent reacceleration where we cap the value of D₀ at D_0, max = 10⁻¹⁸ s⁻¹ (COMA-D_pp). From left to right we show the flux from the simulated magnetic field and the different magnetic field models introduced in Sect. 3.2.

We shall focus on the results of our on-the-fly treatment of CRs here and discuss the comparison to the Hoeft & Brüggen (2007) model shown in the lowest panels of Fig. 6 in Sect. 6.1.

The relic replica in our simulation is only really visible in the zoom-in simulations with the simulated magnetic field. As is shown in Fig. 5, we find that the simulated magnetic field drops off too steeply as a function of radius and with that is not strong enough to illuminate the relic in the lower resolution simulation, even though the CR injection is virtually identical, as is seen in Fig. 4. This drop in magnetic field is even more pronounced in the filaments and with that does not provide diffuse synchrotron emission of any considerable value.

Applying the magnetic field model based on a constant plasma-β (B_β) yields more synchrotron emission in the surrounding medium of the cluster for all simulations. In the case of SLOW-CR3072³ this gives a glimpse at the synchrotron flux of the accretion shock onto the filament in the top right corner of the image. However, it is very dim due to our shock finder being resolution-limited. With the zoom-in simulations, this accretion shock is more pronounced and brighter by roughly one order of magnitude, due to the more accurate shock capturing and with that more consistent acceleration. The inclusion of turbulent reacceleration in COMA-D_pp does not have any considerable impact on the synchrotron flux, due to the small value of D_pp in the filament region. Remnants of a shock from a previous merger, perpendicular to the ongoing merger, can be seen as diffuse emission around the cluster in the zoom-in simulations visible as very light blue regions to the north and south of the cluster.

As was the case with Fig. 3 the model B_ℱ produces morphologically very similar results to B_β. Similar to B_β the model B_ℱ gives increased magnetic field strength at shocks, which boosts the synchrotron emission at ongoing merger- and accretion shocks, compared to the simulated magnetic field. With this model even the accretion shock is visible, enveloping the cluster at a radius of roughly 5 Mpc. Studying the time evolution of the cluster we find that the accretion shock has been pushed out by the previous merger activity, consistent with the findings by Zhang et al. (2020a,b). Further, we find a secondary shock at a radius of roughly 10 Mpc, visible to the right of the cluster in the COMA and COMA-D_pp simulations akin to a runaway shock described by the same authors.

With its steeper density-scaling B_ff does not significantly boost magnetic field strength in the density regimes of accretion shocks and filaments, hence the simulated emission is generally very low. We can again see the impact of resolution between SLOW-CR3072³ and COMA where the improved shockfinder performance increases the simulated synchrotron flux.

For B_dyn,↓ this is even more pronounced, leading to no considerable synchrotron emission in any of the simulations. Given that this model is the only one tuned to match the most recent observations by Carretti et al. (2023) for magnetic field strength in filaments, it makes the detection of diffuse emission very unlikely with current instruments.

B_dyn,↑ generally provides the strongest magnetic field in the filaments, which directly translates to the most synchrotron flux in the case of this model. In SLOW-CR3072³ we again see some of the accretion shock onto the top right filament, while in the zoom-in simulations, this is the magnetic field model that provides some tangible diffuse emission from the filaments connecting to Coma. For COMA, this flux reaches F_ν ∼ 0.01 − 1 μJy beam⁻¹ and with that lies below current observational capabilities by 0.5–1.5 orders of magnitude compared to the noise level achieved in Hoang et al. (2023). Including the modeling of turbulent reacceleration of CRs, with a cap at D₀ = 10⁻¹⁸ s⁻¹ does not change the synchrotron emission significantly in any of the magnetic field models.

4.5. Resolved synchrotron spectrum

The simulation of the CR electron spectra allows us to self-consistently model synchrotron spectra from the injected and aged electron populations. In Fig. 7, we show the resolved spectral index maps between 144 MHz and 1.4 GHz. We constructed these spectra by fitting a single power law to the surface brightness per pixel of the images in Fig. 6 and images with the same resolution obtained at a frequency of 1.4 GHz as

$$\begin{array}{c}\hfill {\alpha }_{1}44MH{z}^1.4GHz(\mathrm{pix})=\frac{{log}_{10}(I_{1}44MHz(\mathrm{pix}))-{log}_{10}(I_1.4GHz(\mathrm{pix}))}{{log}_{10}(144\mathrm{MHz})-{log}_{10}(1.4\mathrm{GHz})}·\end{array}$$(17)

Fig. 7. Resolved synchrotron spectrum maps accompanying Fig. 6. We constructed the synchrotron slope between the images of Fig. 6 and images with the same resolution at 1.4 GHz by fitting a single power law between the respective pixel values. We impose a surface brightness cutoff for the slope construction, similar to the one in Fig. 6 and also cut pixels with slopes α < −3. The lowest row shows histograms of the pixels in the spectral slope maps. Solid lines are simple histograms, dotted lines are synchrotron flux-weighted histograms. The grayed-out area indicates spectral slopes that are too flat for the standard DSA picture.

We impose a surface brightness cutoff similar to the one in Fig. 6, beyond which we do not consider pixels for the slope calculation. This helps avoid polluting the image with extremely low surface brightness emission that nonetheless can show a flat spectrum due to recent injection at high Mach number accretion shocks. In addition, we exclude spectral indices smaller than α < −3, also purely for visualization reasons.

The spectral index maps show two main behaviors. First, we find the typical spectral steepening in the wake of the ongoing merger shocks. This stems from the synchrotron and IC cooling the CR electrons experience after the shock passage. Second, in the most optimistic magnetic field models which show enough diffuse synchrotron emission from the filaments to model a synchrotron slope, we find a spectral index for diffuse emission in filaments of α ∼ −2.0 to –2.5. This is steeper than measurements from recent observations by Vernstrom et al. (2023). The main reason for this is a choice in numerical design: the injection of CRs into our model happens in the middle of the numerical timestep, in order to subtract the injected CR energy from the entropy change in the hydro solver. This is done to conserve the total energy and for numerical stability of the computation (see the discussion in Böss et al. 2023a). All other CR operations, adiabatic changes, turbulent reacceleration, and most relevant here, radiative losses are computed at the end of the timestep for the whole timestep. Snapshots are then written at the beginning of the subsequent timestep. This means that a spectrum that has cooled over at least one time step is written in a snapshot. Due to the adaptive nature of the time-stepping in OPENGADGET3, particles in low-density regions such as filaments are always evolved with larger timesteps than those in high-density regions, such as clusters and shocks in clusters. Due to the fast cooling by IC scattering off the CMB alone, spectra of CRs injected at accretions shocks in filaments will be artificially more cooled than those of CRs accelerated at merger shocks.

The second reason is the wide frequency range over which we consider the resolved synchrotron spectrum in this figure. Going to lower frequencies, below 200 MHz, where the synchrotron spectrum is less affected by cooling processes and does therefore probe the injection-spectrum, we find an excellent agreement with Vernstrom et al. (2023), as we discuss in Sect. 4.6. With our imposed limits for slope construction, however, we only find sufficient diffuse emission in the filament to perform this slope construction for the magnetic field models B_β, B_ℱ and B_dyn,↑.

The lowest panels show the histograms of the synchrotron slope per pixel in the above maps. Solid lines show the simple histograms, while dotted lines show synchrotron flux-weighted histograms. Line colors correspond to the different simulations, as is indicated below the figure. The grayed-out area indicates the range of synchrotron slopes that are too flat to be consistent with DSA. We note that even though some lines seem to be crossing into this area because of their width we find no pixels in the actual maps that show a synchrotron spectrum flatter than α = −0.5, the theoretical limit of DSA.

These histograms generally show a broad distribution of spectral slopes, with the peaks at α > −1 being driven by the ongoing injection in the relics, as can be seen from the synchrotron flux-weighted histograms. For the magnetic field models with strong fields in the filaments and simulations including turbulent reacceleration, we can see a large number of pixels with steeper synchrotron spectra, originating from the diffuse emission around the cluster and in filaments, as was discussed above.

4.6. Filament spectrum

To give a more quantitative impression of the synchrotron spectra arising from the CR electron spectra present in the filaments we plot these spectra in Fig. 8. The position of the selected filament region is sketched in Fig. E.1.

Fig. 8. Momentum and synchrotron spectra of CR electron populations in the filament segment indicated in Fig. E.1. Top panels: Sum of all momentum spectra for particles contained in the selected region. Bottom panels: Synchrotron spectra of the selected regions under the assumption of varying magnetic field models, as is indicated by the colors. In all cases, the simulated magnetic field strength is too low to produce synchrotron emission visible in the plotted range.

The upper panels show the sum of the momentum spectra within the selected region. These spectra are naturally injection-dominated by high Mach-number accretion shocks around the filament and therefore lie close to the maximum possible DSA slope of q = −4. We can see the impact of resolution and with that, shockfinder performance discussed above in the comparison between the left-most panel and the panel to the right of it. The larger normalization of the spectrum is driven by a larger number of SPH particles containing a CR population, leading to a larger sum of the CR number and energy per bin. In the upper right panel an indication of the impact of turbulent reacceleration can be found in the larger normalization at the low-momentum end of the distribution function. This difference is however very small and can also be explained by a slight difference in the evolution of the cluster.

The synchrotron spectra (shown in the lower panels) on the other hand are dominated by the magnetic field strength. In all cases, the simulated magnetic field strength is too low to result in synchrotron power large enough to appear in our plotted range. For our assumed models the magnetic field will always be stronger in the central regions of the filament, where the electrons will settle after some time. During this time the electrons will cool primarily due to IC scattering off CMB photons. Hence the high-frequency end of the synchrotron spectra is steeper than one might expect from the sum of the electron spectra. In the low-frequency range ν < 200 Mhz, the synchrotron spectra are in excellent agreement with the slope α ≈ −1 to –1.5, as is indicated by the dashed black reference lines. In this range, the CR spectra are not very affected by cooling, so it gives a closer representation of the injection spectra, as was discussed above.

For the SLOW-CR3072³ simulation we find the steepest spectra with the lowest synchrotron power, which we attribute to a lack of resolution to capture spurious shocks around the filament. This leads to this spectrum being driven by older CR populations which have advected downstream of the acceleration zone and have undergone considerable cooling. In the case of COMA we see similar behavior; however, the total emission is larger, due to the increased numerical performance of our shock finder and with that more injection of CRs.

In the COMA-D_pp runs we find very comparable spectra, with only a slightly higher total emission in the case of COMA-D_pp. This difference is however so small that it can also be attributed to differences in the evolution of the cluster between the two simulations. In general, we find that most spectra are better fit by two power laws, rather than one. Due to the complex nonlinearities arising from the inclusion of an on-the-fly treatment of CRs in a fully cosmological simulation, no two simulations will be perfectly identical. This leads to a discrepancy between runs of CR injection both in time and space, which with the short-lived nature of high-energy electrons can lead to variations in the synchrotron spectra.

4.7. Individual CR spectrum

To visualize the time evolution of the CR electron spectrum contained in an individual SPH particle we show one such evolution in Fig. 9. The SPH particle is taken from the COMA simulation, since it was run with the highest number of output snapshots. In the left panel, we see the momentum spectra with bin boundaries indicated as dashed lines.

Fig. 9. Momentum- and synchrotron spectrum of an arbitrarily chosen SPH particle with an active CR electron population. We show the spectra at every second available output time between z = 1 and z = 0 in the COMA simulation. Left: Momentum spectrum of the electron population in the momentum range relevant for synchrotron emission. Dashed lines indicate the momentum bin boundaries imposed in our CR scheme. Right: The resulting synchrotron emissivity spectrum of the same electron population.

The spectral evolution indicates that the particle has experienced a shock at least twice: Once around 10 Gyr, and a second time around 14 Gyr. This is in agreement with recent results by Smolinski et al. (2023) who find similar multi-shock scenarios can provide the seeds for bright radio relics. However, none of the output snapshots shown in this figure have a detected shock for the SPH particle, indicating that all acceleration happened between the available output times. This would lead to an underestimation of synchrotron emission in a pure post-processing approach. We address this again in Sect. 6.1.

Between the acceleration episodes, the spectrum is subject to adiabatic expansion, as is evidenced by the left shift of the spectrum at a constant spectral slope. Energy losses can be seen as the spectral cutoff moving to lower momenta and a steepening of the spectrum below this cutoff. We evolve this spectral cutoff as an explicit variable, allowing us to represent partially filled momentum bins, as is evident multiple times in the shown spectra.

In the right panel of Fig. 9, we show the synchrotron emissivity as a function of frequency for the corresponding CR electron spectra in the left panel. With readability in mind, we only show the synchrotron spectra based on the simulated magnetic field and omit the other magnetic field models here.

For the synchrotron spectra corresponding to recently injected momentum spectra, we find unbroken power laws up to the 10 GHz shown here, most evident by the spectrum injected at the second acceleration event. They show slopes close to the theoretical maximum of DSA, consistent with the momentum spectra approaching the same limit. As the momentum spectra cool from synchrotron and IC cooling we can see a corresponding steepening in the synchrotron spectra. This primarily affects the high-frequency end of the synchrotron spectra, with the low-frequency end being dominated by adiabatic changes in the gas and variations in magnetic field strength.

5. Full box

5.1. Cosmic ray electron injection

In Fig. 10, we show the density-temperature phase-space of all particles contained in SLOW-CR3072³. The background shows the total mass contained per phase-space element, indicated by the color bar on top. Contours show the total CR energy per phase-space element, contained in electrons with energies above 1 GeV. These contours show two main acceleration regions. One in the moderate density, high-temperature region associated with merger shocks in the cluster periphery. These CRs could be seen in the form of radio relics, given sufficient magnetic field strength. The other is in the low density and moderate temperature region of the warm-hot intergalactic medium (WHIM) with temperatures T_WHIM ∼ 10⁵ − 10⁷ K (Cen & Ostriker 1999; Davé et al. 2001) where we can expect accretion shocks onto clusters and filaments. This provides a basis for potential synchrotron emission in cosmic web filaments, provided the magnetic field is sufficiently strong.

Fig. 10. Density – temperature phase diagram of all gas particles in the simulation. The colored background shows the mass histogram in phase space. The contours show the CR electron energy of CRs with energies above 1 GeV contained per phase-space element.

5.2. Synchrotron emissivity

Figure 11 shows the density-synchrotron emissivity phase space for the simulated and modeled magnetic field. As in Fig. 10, the background shows the total mass contained per phase-space element, here for the emissivity calculated from the simulation SLOW-CR3072³. We only calculated the mean of the particles with emissivities above j_ν > 10⁻⁵² erg s⁻¹ Hz⁻¹ cm⁻³, since the majority of the SPH particles in the simulations do not contain a CRe population, so the mean including all particles is close to zero. Lines show the mean value per density interval for SLOW-CR3072³ and the zoom-in simulations COMA and COMA-D_pp. The arrow indicates the upper limit estimated in Hoang et al. (2023) for n_e = 10⁻⁵ cm⁻³. We list the values for j_ν we find in this density bin in Table 2.

Fig. 11. Individual panels showing the phase-space diagrams in the density-synchrotron emissivity phase-space, containing all gas particles in the simulation. The background shows mass histograms per phase-space element from the original simulation output. The colored lines indicate the mean synchrotron emissivity per density interval for the cosmological box simulation and the zoom-in simulations with and without turbulent reacceleration. All emissivities are calculated at a frequency of ν = 144 MHz. The upper limit of the synchrotron emissivity is taken from Hoang et al. (2023).

In the case of SLOW-CR3072³, we fall below the upper limit of Hoang et al. (2023) by at least one order of magnitude for all considered magnetic field configurations. As was discussed in the previous section, this can be attributed to the lack of resolution to capture all of the spurious accretion shocks. Lower limits in SLOW-CR3072³ can be found with the models B_ff and B_dyn,↓, which both show synchrotron emissivities three orders of magnitude below the observational limits. The more optimistic magnetic field models B_β, B_ℱ and B_dyn,↑ still lie significantly below the observational limit by 1.5–2 orders of magnitude.

The zoom-in simulations again show the impact of resolution. With the shock finder performing better at this resolution the spurious accretion shocks are found more consistently, which leads to a consistent increase in synchrotron emissivity in low-density regions. This can be seen in the bump of synchrotron emissivity below n_e ∼ 5 × 10⁻⁴ cm⁻³ which is where we expect to see the turnover from internal merger shocks to accretion shocks (e.g., Vazza et al. 2011). The inclusion of turbulent reacceleration does not have a noticeable impact on the emissivity in filaments, at least within the modeling present in this work, as was mentioned above. Future work will need to account for different reacceleration models to study the potential impact on filaments.

In the case of the simulated magnetic field model B_sim we find that the mean emissivity in COMA is in agreement with the emissivity limit. We note, however, that this comparison has its limits since we are comparing mean emissivities between a volume-filling model in Hoang et al. (2023) and spurious shocks in this work.

Compared with the other magnetic field models we find that the increase in resolution in COMA consistently increases the mean synchrotron emissivity by two orders of magnitude compared to SLOW-CR3072³. As this increase is magnetic field model independent we can conclude that it is injection driven, meaning that we capture more of the spurious accretion shocks in the zoom-in simulation COMA, compared to the full cosmological box SLOW-CR3072³.

The lower magnetic field models B_ff and B_dyn,↓ both result in emissivities orders of magnitude below the observational limit. Especially in the case of B_dyn,↓, which matches the 30 nG field inferred by Carretti et al. (2023) at n_e ≈ 10⁻⁵ cm⁻³, this underlines the point of a significant challenge for observations to detect diffuse emission, even by stacking.

For B_β and B_ℱ we find emissivities above the observational limit, we note however that this model is biased toward accretion shocks due to its scaling with the thermal and turbulent pressure. B_dyn,↑ again is roughly in agreement with the observational limit.

As we discuss in Sect. 6, one to two orders of magnitude of synchrotron emissivity can be covered by a variation in the choice of CR injection parameters, giving the possibility of diffuse emission just below the detection limit, even with more conservative magnetic field estimates.

6. Discussion

6.1. Comparison to post-processing

First, we want to compare the result from our on-the-fly solver to the post-processing model by Hoeft & Brüggen (2007). We use their acceleration efficiency function Ψ and ξ_e = 10⁻⁵, similar to Nuza et al. (2017). We applied their model to the COMA simulation to make use of the improved shock finder performance compared to SLOW-CR3072³ and show the results in the lowest panels of Fig. 6.

For the relics in our Coma replica, we find that the model by Hoeft & Brüggen (2007) provides reasonable results and most notably, sharper relics. While in our simulations the region between cluster center and relic is filled with electrons that have been accelerated in the shock passage and subsequently advected downstream of the shock, a post-processing model only accounts for the acceleration at the currently detected shock. This naturally leads to thinner shocks, as the downstream CRs are not accounted for.

However, when it comes to the diffuse model synchrotron emission of filaments a post-processing model faces a number of limitations. For one it is naturally coupled to the performance of the shock finder at the time of output. As can be seen in the panels for B_β, B_ℱ and B_dyn,↑, which produce volume-filling synchrotron emission with our on-the-fly treatment of CRs, this emission is comparable in absolute flux where it is present; however, it is confined only to the regions of active shock detection. This leads to very spurious emission only at the detected accretion shocks. Shock detection in simulations is numerically challenging and requires a number of filtering steps to distinguish between real shocks, cold fronts, and shearing flows (see e.g., Vazza et al. 2009; Schaal & Springel 2015; Beck et al. 2016b). This leads to highly transient shock surfaces and introduces difficulties in finding consistent shock surfaces in and around galaxy clusters and filaments at a single timestep. To add to the complication this naturally scales with numerical resolution, making it easier to detect shocks in the higher resolved intracluster medium (ICM) than in the lower resolved warm-hot intergalactic medium (WHIM) (see e.g., Vazza et al. 2011, for a comparison between different numerical models).

While this is true for every individual timestep of the simulation, we find that our on-the-fly treatment of CRs, and the fact that they can potentially be injected at every timestep in every particle, helps to alleviate this problem. As a SPH particle moves through the region of an accretion shock, it is highly likely that the shock will be detected at some point through the shock passage and we are not necessarily reliant on the detection happening at the output timestep. This leads to the smoother injection of CRs and with that a more diffuse emission in the filament volume, as can be seen in the direct comparison between the lower two rows of the rightmost column in Fig. 6.

We also find that the advection of CRs with the thermal gas plays a role in the volume-filling synchrotron emission of filaments. This can be seen in a direct comparison between, for example, the rightmost row for COMA, COMA-D_pp and the post-processing approach in Fig. 6 where we find more volume-filling emission in the case of the on-the-fly treatment of CRs. Naturally, a pure post-processing model cannot account for this process and is therefore entirely driven by in situ emission. This has been addressed by multiple groups and has led to several recent developments of Fokker-Planck solvers with or without tracer particles to account for CR electron injection and transport in large-scale simulations (e.g., Wittor et al. 2017; Yang & Ruszkowski 2017; Winner et al. 2019; Ogrodnik et al. 2021; Vazza et al. 2021; Hopkins et al. 2022).

6.2. Choice of cosmic ray injection parameters

In this work, we have used a fixed set of CR parameters and evolved the injected CR population over time. These parameters are, namely, the injection momentum, ${\widehat{p}}_{\mathrm{inj}}$, as the momentum where the nonthermal power law of the CR electron population starts, the efficiency, η(ℳ_s, X_cr, θ_B), as the fraction of shock energy that is available to accelerate CRs to relativistic energies, and finally the ratio between proton and electron injection, K_ep. Most of these parameters are poorly constrained in the high Mach number, high-β plasma shocks expected in and around cosmic web filaments, with only a small number of studies available (e.g., Kang & Ryu 2013; Guo et al. 2014; Ryu et al. 2019; Kang et al. 2019; Kobzar et al. 2021; Ha et al. 2023). We shall briefly discuss the impact of these parameters on the synchrotron emission in cosmic web filaments. An overview of our parameters and the ones used in previous work in the literature can be found in Table 3.

6.2.1. Injection momentum, ${\widehat{p}}_{\mathrm{inj}}$

The strongest impact is caused by the choice of the connection point between thermal and nonthermal electron populations, ${\widehat{p}}_{\mathrm{inj}}$, beyond which the power-law distribution starts. Typical fixed values used in previous work on this point lie around ${\widehat{p}}_{\mathrm{inj}}=0.01-0.1$ (e.g., Hong et al. 2014; Donnert et al. 2016). For the electrons to be efficiently accelerated by DSA, their gyro radius must be on the order of that of the protons, which means a dependence of the injection momentum on the magnetic field and the temperature downstream of the shock.

Kang & Ryu (2013), for example, model this as

$$\begin{array}{c}\hfill p_{\mathrm{inj}}\approx 1.17m_p{v}_2(1+\frac{1.07}{{ϵ}_B}),\end{array}$$(18)

where m_p is the proton mass, v₂ is the downstream velocity of the shock, and ϵ_B = B₀/B_⊥ is the ratio between the mean magnetic field along the shock normal and the amplitude of the post-shock MHD turbulence.

The Hoeft & Brüggen (2007) model omits the explicit magnetic field dependence and relates the injection momentum only to the temperature downstream of the shock (T₂). This leads to an injection momentum of

$$\begin{array}{c}\hfill {\widehat{p}}_{\mathrm{inj},\mathrm{e}}\sim \frac{10k_B{T}_2}{m_e{c}^2},\end{array}$$(19)

which for typical shocks in the ICM with T₂ ∼ 10⁸ K is of the order ${\widehat{p}}_{\mathrm{inj},\mathrm{e}}\sim 0.1-0.2$.

A similar parametrization is used, for example, in Kang et al. (2019), Vazza et al. (2021), Ha et al. (2023) with

$$\begin{array}{c}\hfill {\widehat{p}}_{\mathrm{inj},\mathrm{e}}\approx \frac{3\sqrt{\frac{m_e}{m_p}}{p}_{\mathrm{th},\mathrm{p}}}{m_{e}c}=\frac{3\sqrt{2k_B{m}_e{T}_2}}{m_{e}c},\end{array}$$(20)

which for the shocks in our simulated cosmic web filaments with typical downstream temperatures of T₂ ∼ 10⁷ K would be of the order ${\widehat{p}}_{\mathrm{inj}}\approx 0.17$. This makes our choice of a fixed ${\widehat{p}}_{\mathrm{inj}}=1.0$ an optimistic assumption, especially for cold accretion shocks. We show the impact of the choice of the injection momentum on the synchrotron emissivity in Fig. 12. Given the high sonic Mach number, we detect in the surrounding of filaments with ℳ_s ≫ 10 the standard RH jump conditions converge toward q = 4. For this value of q, even an injection momentum two orders of magnitude higher only introduces a synchrotron emission half an order of magnitude higher.

Fig. 12. Synchrotron emissivity of an electron spectrum with a power-law slope, q, extending between an injection momentum, ${\widehat{p}}_{\mathrm{inj}}$, and a maximum momentum of ${\widehat{p}}_{\max }={10}^6$. Each spectrum is normalized to contain a total energy density of ϵ = 1 erg cm−3 and emits in a magnetic field with B = 1 μG. The color shows emissivity and contour lines are spaced apart by half an order of magnitude.

6.2.2. Acceleration efficiency, η(ℳ_s, X_cr, θ_B)

Mach number-dependent efficiency models that have been used in previous simulations typically show a strong dependence on Mach number for weak shocks ℳ_s < 5 and saturate at a fixed value for strong shocks ℳ_s > 7 − 10 (e.g., Hoeft & Brüggen 2007; Kang et al. 2007; Kang & Ryu 2013; Caprioli & Spitkovsky 2014; Ryu et al. 2019). This saturation value varies by about one order of magnitude (see the left panel of Fig. A.1 for our parametrization of the different DSA models). For accretion shocks, whose Mach numbers lie within the range of the efficiency models above the saturation value this means that the injected energy and with that the synchrotron emission scales linearly with this saturation value. With our choice of the Ryu et al. (2019) efficiency model, we lie at the lower end of the range for the saturation values. A different model choice would therefore increase the potential synchrotron brightness by up to a factor of 25. We note, however, that DSA models for high sonic Mach number, high-β, and low-temperature shocks are scarce, and further investigation in this direction is needed.

In the case of magnetic field angle dependent acceleration efficiency (η(θ_B)), Vazza et al. (2017), Wittor et al. (2017) find that allowing CR electron acceleration only at quasi-perpendicular shocks with θ > 60° does not significantly alter the expected emission from accretion shocks onto clusters and cosmic web filaments. This stems from the alignment of the shock surface and magnetic field, due to the adiabatic compression and with that amplification of the magnetic field only working on the perpendicular component. In Fig. 13, we show the density–shock obliquity phase space of SLOW-CR3072³ and COMA. We find a similar trend of shock obliquity increasing for low-density regions.

Fig. 13. Density – shock obliquity phase map. The colors indicate the 2D mass histogram for SLOW-CR30723, while the lines indicate the maximum of the histograms for SLOW-CR30723 and COMA as a function of density.

However, we also see a shock obliquity of θ_B ≈ 60° for higher densities. This result is in agreement with the findings of Banfi et al. (2020, see their Fig. 17). Wittor et al. (2017) show that for any upstream shock obliquity, θ_pre, the shock passage will always align the obliquity toward θ = 90°. This process only depends on the upstream obliquity and the shock compression ratio, which makes it more efficient for higher Mach number shocks.

Hence, if accounting for a shock obliquity-dependent acceleration efficiency, we preferentially accelerate CR electrons over CR protons. This can help to ease the tension between CR acceleration efficiency models and the lack of diffuse γ-ray detection in galaxy clusters (see e.g., Wittor et al. 2019, for a discussion). We shall address this question as well in upcoming work, analyzing the proton component of our SLOW-CR3072³ simulation.

6.2.3. Electron-to-proton ratio, K_ep

The final free parameter is the electron-to-proton injection ratio, K_ep. This is poorly constrained and is typically taken as a canonical value from SNe as K_ep = 0.001 − 0.01 (e.g., Park et al. 2015). Kang (2020) follow an alternative approach by making this parameter injection slope dependent as

$$\begin{array}{c}\hfill K_{\mathrm{ep}}={\left(\frac{m_e}{m_p}\right)}^{(q_{\mathrm{inj}}-3)/2}.\end{array}$$(21)

At strong accretion shocks with an injection slope of q_inj ≈ 4, this gives a value of K_ep ≈ 0.023. This means we underestimate the injection of CRe by a factor of 2.3 by using our constant value. However, it is not clear what role nonlinear DSA effects can play in these shock environments and if this can lead to a steepening of the injection slope (see Caprioli et al. 2020; Diesing & Caprioli 2021). Given the sensitivity of K_ep on the injection slope and the uncertainty of the real slope, we note that further study is required.

Vazza et al. (2015a) test this value for radio relics and find that to explain observed radio emission they require K_ep ≥ 10 − 10⁻² to avoid γ-ray emission in conflict with Fermi-LAT observations. Recent results indicate that this over-injection of CR protons can also be remedied by the shock-obliquity dependent acceleration of CRs (see Wittor 2021). Given that the observed brightness of radio relics is still an open problem, this variance in values for K_ep introduces a large variance in potential synchrotron emission. The emission in our model scales linearly with K_ep, so larger values for K_ep can significantly boost the synchrotron emission in filaments. Recent PIC simulations by Gupta et al. (2024) indicate that electron acceleration may become more efficient compared to proton acceleration at high Mach number shocks, leading to K_ep ∼ 0.1. This would increase the synchrotron emission by an order of magnitude compared to our current modeling. We shall investigate the impact of different values for K_ep in future zoom-in simulations.

6.3. Impact of magnetic field models on synchrotron cooling

On the point of cooling, specifically cooling due to synchrotron losses, we note that the spectral evolution should in principle be performed again for each of our different magnetic field models. Since we change the magnetic field model only in post-processing, while keeping the CR spectra the same we can underestimate the cooling due to synchrotron emission that the CR spectra experience in cases where the modeled magnetic field is higher than the simulated one.

However, we find that in the vast majority of the particles containing a CR population cooling is IC scattering dominated since B < B_CMB ≈ 3.2 μG (1 + z)⁴. As can be seen in Fig. 5, in the example of Coma this concerns all particles beyond the inner ∼500 kpc. Since, for this work, we are primarily interested in the synchrotron emission of cosmic web filaments, whose magnetic field is expected to lie significantly below B_CMB, we aim to discuss the impact of varying magnetic field models on cooling in this section.

In Fig. 14, we show the histograms of magnetic field strength for all SPH particles contained in the selection used to obtain Fig. 9. In addition, we limit the histogram computation to particles that also contain a CR electron population. The panels correspond to the four simulations used in this work and colors correspond to magnetic field models. The dashed black line indicates the CMB equivalent magnetic field at z = 0. All histograms are normalized to the total number of SPH particles with a CR electron population.

Fig. 14. Histograms of magnetic field strength in all SPH particles containing a CR electron population in the cutout selected for Fig. 9. Colored lines correspond to the magnetic field models, the dashed black line shows the equivalence magnetic field of the CMB at z = 0.

In all cases, the relative number of CR spectra that would experience first-order effects from stronger synchrotron cooling in the case of stronger magnetic field models is of the order N/N_tot ∝ 𝒪(10⁻³). We further note that for these particles to have such high magnetic fields in our models they have to lie close to the center of the filament and should hence have undergone significant cooling due to IC scattering since their acceleration at the accretion shock and subsequent advection into the filament. This reduces the significance of the additional cooling due to synchrotron emission. We therefore do not expect a significant impact on our results from not rerunning the spectral evolution with the various magnetic field models.

6.4. Impact of AGN feedback and star formation

With our simulation being non-radiative, we do not account for AGN feedback or star formation. We therefore also cannot account for CR injection by these processes. In the following, we briefly discuss the potential impact of this on our results. An AGN can potentially provide a significant CR electron seed population that can advect and diffuse in the cluster volume (see e.g., ZuHone et al. 2021; Vazza et al. 2023; Domínguez-Fernández et al. 2024). However, this advection is expected to be confined to clusters, where the most powerful AGNs are expected to be located, not filaments, which are the primary focus of this work. There the injected electron can potentially be reenergized by cluster merger shocks and can be observed as the so-called radio phoenices (van Weeren et al. 2019). The jets of AGNs can also inject significant turbulence, as for example is observed in Perseus (Sanders et al. 2016). This turbulence can then reaccelerate fossil electrons and provide the basis for radio halos.

We aim to include AGN feedback in the near future, as the modeling of radio galaxies in our sample should help to estimate the background pollution by discreet sources as well as additional feedback on the cluster gas by AGN-driven shocks and their impact on our results. This is, however, beyond the scope of this work.

Cosmic ray proton injection at SNe has been found to drive and/or modify the structure of galactic winds and the circumgalactic medium (CGM) (e.g., Pfrommer et al. 2017; Butsky & Quinn 2018; Fielding et al. 2020; Ji et al. 2020; Girichidis et al. 2024). Cosmic ray protons generally have lifetimes exceeding the Hubble time in the low-density CGM and ICM (see e.g., Fig. 7 in Brunetti & Lazarian 2007, for an illustration of CR proton cooling times at different densities), which under the assumption of efficient transport can provide a volume-filling CR proton population. There they can nonetheless interact with the thermal background gas and be scattered into π-ons and, depending on the charge of the π-on, decay into γ-photons, or electrons or positrons. These electrons could then contribute to diffuse synchrotron emission. In the context of CR proton acceleration at merger shocks, this process has been proposed as a source of radio halos (e.g., Dolag & Enßlin 2000). However, the non-detection of diffuse γ-ray emission from galaxy clusters puts a strong upper limit on this process (e.g., Pinzke & Pfrommer 2010; Pinzke et al. 2011; Ackermann et al. 2014). For this same reason, we omitted including an explicit treatment of secondary electrons.

6.5. Observational prospects

Given the potential impact of so-far poorly constrained CR injection parameters on our results, we can make a best-case assumption of the observational requirements for direct observation of diffuse synchrotron emission from cosmic web filaments. We shall discuss this again for the filaments we find around our Coma cluster replica, as it provides the widest range of parameters in our set due to the zoom-in re-simulation. As is shown in Fig. 6, the most optimistic magnetic field models, normalized to the observed synchrotron flux from the Coma relic in the higher resolution simulations leads to a diffuse synchrotron flux from the filaments of S_{ν = 144 MHz} ≈ 0.1 μJy beam⁻¹ for a beam size of θ = 60″ × 60″. This is three orders of magnitude below the noise of σ_rms = 0.4 mJy beam⁻¹ given in Bonafede et al. (2022). However, the choice of other parameters for p_inj, η(ℳ_s) and K_ep could boost this emission by up to two orders of magnitude. This would indicate diffuse emission just a factor of ten below the current noise level. The stacking approach by Hoang et al. (2023) reduces the noise to σ_noise = 0.6 μJy beam⁻¹, which lies just a factor of ∼ten above the diffuse emission we find in the most optimistic scenario.

On top of surface brightness limitations, we further need to account for a suppression of observability, due to the large structures we consider and missing baselines of LOFAR. These missing baselines suppress the observed flux from synchrotron emission on large scales and increase it on small scales. To estimate the impact of missing baselines on the observability, we can follow the approach of, for example, Churazov et al. (2022), using Eq. (3.28) in Remijan et al. (2019). This gives an expression for the scale θ_MRS at which the flux at a given frequency, ν_obs, is suppressed by the minimum baseline, D_min, as θ_MRS ≈ 0.6c/(νD_min). Using a minimum baseline of D_min = 0.1 km from Bonafede et al. (2021), we obtain θ_MRS ≈ 42.94′≈1.24 Mpc for our Coma replica.

To calculate the contribution to the emission, we folded our non-smoothed images, I, of Fig. 6 with a Gaussian beam of size θ_MRS to obtain an image, I′. This gives a relative contribution (I − I′)/I and leads to more filamentary emission in the cluster periphery while suppressing diffuse emission in the cluster and filament centers.

The relative contribution can be seen in the top panels of Fig. 15 in the case of COMA. We chose this simulation as an example since it showed the strongest synchrotron emission in the filaments. To obtain a more realistic instrument response for LOFAR we added the missing baseline contribution to the original flux images and show the result in the bottom panels of Fig. 15.

Fig. 15. To illustrate the impact of missing baselines on our simulated flux images, mock observations of COMA as in Fig. 6 under these considerations. Top: Relative change to the flux due to missing baselines in LOFAR. Bottom: Synchrotron flux under the consideration of instrument response from missing baselines.

As was mentioned above, this addition has two notable effects. First, it strongly suppressed the diffuse flux from the central parts of filaments, and to a lesser extent the downstream regions of our relic analogs. Second, it increases flux originating from the collision between accretion and merger shocks we find close to the virial radius. These shock structures are morphologically similar to the accretion relic candidate reported in Bonafede et al. (2022); however, in our simulation we find that the emission to the upper right and lower left that most closely resembles the accretion relic originates from a previous, small merger along this axis. A detailed study of the origin of a possible accretion relic is however beyond the scope of this work. The inclusion of missing baselines in the consideration of observational prospects further reduces the chance of observing this emission, at least in nearby filaments.

7. Conclusions

We have presented results from the first constrained cosmological simulation with an on-the-fly Fokker-Planck solver to study the time evolution of CR electron populations on a cosmological scale. We used this simulation to study potential diffuse synchrotron emission from cosmic web filaments. Our results can be summarized as follows:

• We find that the acceleration of CR electrons occurs within clusters as well as in accretion shocks onto clusters and cosmic web filaments. The electrons with the longest lifetimes settle into large halos around clusters and in cosmic web filaments. Given enough magnetic field strength and stabilization against cooling losses, by for example Fermi-II reacceleration, these electrons can provide the basis for synchrotron emission in these regions.

• We tested six different magnetic field configurations to study the threshold for obtaining considerable diffuse emission from cosmic web filaments and found that with our choice of CR injection parameters we require fields of the order B ∼ 100 nG to obtain considerable emission.

• To test the impact of resolution in our simulation and the inclusion of on-the-fly turbulent reacceleration, we performed zoom-in simulations of our Coma cluster replica. We find that higher resolution, and with that a more accurate capture of accretion shocks, increases the synchrotron emission by an order of magnitude in low-density regions. This would in turn decrease the required magnetic field strength to obtain the same level of synchrotron emission as in the low-resolution full box by three orders of magnitude, bringing it into agreement with current observations with B ≈ 30 nG.

• The inclusion of turbulent reacceleration in our modeling does not increase the synchrotron emission in filaments, at least in the gyro-resonant flavor used for this work.

• With our chosen parameters for CR injection, the surface brightness of diffuse emission from cosmic web filaments lies three orders of magnitude below the current noise limit for direct detection with LOFAR.

• In the most optimistic model for magnetic field strength and CR injection parameters, we can estimate synchrotron emission only a factor ∼5 below the noise level currently achievable by stacking of radio images.

• In the models that provide significant synchrotron emission from the filaments, we find that the spectral index of diffuse emission in the LOFAR band is in the range of α ≈ −1.0 to –1.5, in agreement with recent observations by Vernstrom et al. (2023).

In future work, we shall present the results of SLOW-CR3072³ for the radio halos of our Coma, Perseus, and Virgo replicas, radio halo relations, and γ-ray emission from the CR proton population. We shall further increase our sample size of zoom-in simulations and vary the parameters for CR injection and turbulent reacceleration to study their impact on the time evolution of radio halos.

Acknowledgments

We want to thank the anonymous referee for helpful comments that improved the quality of this work. We thank Sebastian Nuza for providing us with the fit function for the Hoeft & Brüggen (2007) model. LMB would like to thank Kamlesh Rajpurohit, Jennifer Schober, Mathias Hoeft, Sebastian Nuza and Damiano Caprioli for helpful discussions. LMB, KD and IK acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311 as well as support for the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. KD, EHM, BS and JS acknowledge support by the grant agreements ANR-21-CE31-0019 / 490702358 from the French Agence Nationale de la Recherche / DFG for the LOCALIZATION project. UPS is supported by the Simons Foundation through a Flatiron Research Fellowship at the Center for Computational Astrophysics of the Flatiron Institute. The Flatiron Institute is supported by the Simons Foundation. The calculations for the hydro-dynamical simulations were carried out at the Leibniz Supercomputer Center (LRZ) under the project pn68na.

References

Appendix A: Acceleration models

We show the DSA parametrizations present in CRESCENDO in Fig. A.1. For this work, we used the sonic Mach number-dependent acceleration model by Ryu et al. (2019), which is shown as the orange line in the left panel. Acceleration efficiency from the thermal pool is shown in solid lines and reacceleration efficiency in dashed lines. We indicate typical ranges of sonic Mach numbers for merger- and accretion shocks with gray bands.

Fig. A.1. DSA parametrizations present in our model. Left: Sonic Mach number ℳs dependent acceleration models. Right: Shock obliquity (θB) dependent acceleration model.

The right panel shows the shock obliquity-dependent acceleration model. For protons, we use the parametrization by Pais et al. (2018), who fit the data by Caprioli & Spitkovsky (2014). The study in Caprioli & Spitkovsky (2014) was performed for protons, which are found to be preferentially accelerated at quasi-parallel shocks. Electrons are found to be preferentially accelerated by quasi-perpendicular shocks, so we take the simple approach of shifting the Pais et al. (2018) model by 90°. Further studies of the exact dependence of electron acceleration on θ_B are needed to improve this simple model.

Appendix B: Time integration of the CR spectra

We shall briefly recap the time integration of CRESCENDO, to provide the relevant equations for the next section (for details, see Böss et al. (2023a)). To update the distribution of the CRs in the two-moment approach we need to compute number- and energy changes per bin over a timestep. We can obtain the number of CRs per bin by performing a volume-integral in momentum space in shells of the with of the individual bins.

$$\begin{array}{c}\hfill N_i=\frac{1}{\rho }{\displaystyle \underset{p_i}{\overset{p_{i+1}}{\int }}}dp4\pi p^{2}f(p)=\frac{4\pi f_i{p}_i^3}{\rho }\frac{({\left(\frac{p_{i+1}}{p_i}\right)}^{3-q_i}-1)}{3-q_i}\end{array}$$(B.1)

We calculated the energy contained in each bin in the ultra-relativistic limit T(p)≈pc and attached this energy to every CR in our previous volume-integral to obtain

$$\begin{array}{c}\hfill E_i=\frac{1}{\rho }{\displaystyle \underset{p_i}{\overset{p_{i+1}}{\int }}}dp4\pi c{p}^{3}f(p)=\frac{4\pi c{f}_i{p}_i^4}{\rho }\frac{({\left(\frac{p_{i+1}}{p_i}\right)}^{4-q_i}-1)}{4-q_i}.\end{array}$$(B.2)

For the time integration, Eq. 2 was rewritten in terms of number- and energy changes per bin i as

$$\begin{array}{c}\hfill \frac{D{N}_i}{\mathit{Dt}}=\frac{1}{\rho }{\left[(\frac{1}{3}\frac{\partial u}{\partial x}p+b_l(p))4\pi p^{2}f(p)\right]}_{p_i}^{p_{i+1}}\end{array}$$(B.3)

$$\begin{array}{cc}\hfill \frac{D{E}_i}{\mathit{Dt}}=& \frac{1}{\rho }{\left[(\frac{1}{3}\frac{\partial u}{\partial x}p+b_l(p))4\pi c{p}^{3}f(p)\right]}_{p_i}^{p_{i+1}}\hfill \end{array}$$(B.4)

$$\begin{array}{cc}& -(\frac{4}{3}\frac{\partial u}{\partial x}{E}_i+\frac{1}{\rho }{\displaystyle \underset{p_i}{\overset{p_{i+1}}{\int }}}dp{b}_l(p)4\pi c{p}^{2}f(p))\hfill \end{array}$$(B.5)

where we introduced the quantity R_i(q_i, p_i) for the energy loss integral per bin

$$\begin{array}{c}\hfill R_i(q_i,p_i)=\frac{4-q_i}{p_{i+1}^{4-q_i}-p_i^{4-q_i}}{\displaystyle \underset{p_i}{\overset{p_{i+1}}{\int }}}dp{p}^{2-q_i}(\frac{1}{3}\frac{\partial u}{\partial x}+{\displaystyle \sum _l^{N_{\mathrm{losses}}}}{b}_l(p))\end{array}$$(B.6)

for simplification. Integrating this in times gives

$$\begin{array}{cc}\hfill N_i^{t+\mathrm{\Delta }t}& =N_i^t+\frac{1}{\overline{\rho }}(F_{N_{i+1}}^m-F_{N_i}^m)\hfill \end{array}$$(B.7)

$$\begin{array}{cc}\hfill E_i^{t+\mathrm{\Delta }t}& \approx E_i^t(1-R_i(q_i,p_i))+\frac{1}{\overline{\rho }}(F_{E_{i+1}}^m-F_{E_i}^m)\hfill \end{array}$$(B.8)

where we can associate

$$\begin{array}{cc}\hfill F_{N_i}^m& ={\displaystyle \underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}}d{t}^{\prime }{b}_l(p)4\pi p^{2}f(t^{\prime },p){|}_{p_i}\hfill \end{array}$$(B.9)

$$\begin{array}{cc}\hfill F_{E_i}^m& ={\displaystyle \underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}}d{t}^{\prime }{b}_l(p)4\pi c{p}^{3}f(t^{\prime },p){|}_{p_i}\hfill \end{array}$$(B.10)

with the time-averaged fluxes over the momentum bin boundaries. We then transfer the problem into momentum space by substituting dt′ with the momentum change over time. This gives equations for the fluxes that depend only on the momentum of a particle p_u at which it must have started out before a timestep Δt to arrive at p_i.

$$\begin{array}{cc}\hfill F_{N_i}^m& ={\displaystyle \underset{p_i}{\overset{p_u}{\int }}}dp4\pi p^2{f}^m(p)\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill F_{E_i}^m& ={\displaystyle \underset{p_i}{\overset{p_u}{\int }}}dp4\pi c{p}^3{f}^m(p),\hfill \end{array}$$(B.12)

with

$$\begin{array}{cc}\hfill f^m(p)& =\{\begin{array}{cc}{f}_i{\left(\frac{p}{p_i}\right)}^{-q_i}\hfill & \mathrm{if}{p}_u>p_i\hfill \\ f_{i-1}{\left(\frac{p}{p_{i-1}}\right)}^{-q_{i-1}}\hfill & \mathrm{if}{p}_u\le p_i\hfill \end{array}\hfill \end{array}$$(B.13)

After the time integration is finished we end up with new values for CR number and energy per bin. To reconstruct the distribution function we solve

$$\begin{array}{c}\hfill \frac{E_i}{N_i{p}_{i}c}=\frac{3-q_i}{4-q_i}\frac{{\left(\frac{p_{i+1}}{p_i}\right)}^{4-q_i}-1}{{\left(\frac{p_{i+1}}{p_i}\right)}^{3-q_i}-1}\end{array}$$(B.14)

for q_i numerically using Brent’s method and obtain the updated norm by solving Eq. B.1 for f_i as

$$\begin{array}{c}\hfill f_i=\frac{\rho N_i}{4\pi p_i^3}\frac{3-q_i}{{\left(\frac{p_{i+1}}{p_i}\right)}^{3-q_i}-1}\end{array}$$(B.15)

Appendix C: Fermi-II reacceleration

As was discussed above, to model the turbulent reacceleration of CRs we adopt the model by Cassano & Brunetti (2005)². We only account for the systematic component of reacceleration, which can be expressed in the unified-cooling approach, extending the b_l(p) term in Eq. 3 to include the momentum change due to systematic reacceleration

$$\begin{array}{c}\hfill {\displaystyle \sum _l}{b}_l(p)\equiv {\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{cool}}^{\mathrm{synch}}+{\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{cool}}^{\mathrm{IC}}+{\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{acc}}^{\mathrm{sys}}.\end{array}$$(C.1)

This momentum change over time can be expressed as

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{acc}}^{\mathrm{sys}}=-\chi p\approx -2\frac{D_{\mathrm{pp}}}{p}.\end{array}$$(C.2)

For an isotropic distribution of wave and particle momenta and v_A < v_M, with v_A being the Alfvén velocity and

$$\begin{array}{c}\hfill v_{\mathrm{M}}^2\approx \frac{4}{3}{v}_{\mathrm{ion}}^2+v_{\mathrm{A}}^2\end{array}$$(C.3)

the momentum diffusion coefficient, D_pp, is given by Eilek (1979) as

$$\begin{array}{cc}\hfill D_{\mathrm{pp}}(p,t)& \approx 4.45{\pi }^2\frac{v_{\mathrm{M}}^2}{c}\frac{p^2}{B^2}{\displaystyle \underset{k_{\min }}{\overset{k_{\max }}{\int }}}\mathrm{d}kk{\mathcal{W}}_k^B(t)\hfill \end{array}$$(C.4)

$$\begin{array}{cc}& =D_0(t)p^2.\hfill \end{array}$$(C.5)

This assumes that the infall of substructure causes plasma instabilities and turbulence which in turn drives magneto-sonic (MS) waves. These instabilities lead to the systematic second-order Fermi reacceleration process. The integral on the r.h.s. of Eq. C.4 is therefore the integral over the spectrum of the MS waves. The boundaries of this integral are defined by the minimum and maximum wave numbers of the underlying turbulence spectrum. Combining Eq. C.2 and Eq. C.5 gives the short form of

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}p}{\mathrm{d}t}\right)}_{\mathrm{acc}}^{\mathrm{sys}}=-2D_0(t)p.\end{array}$$(C.6)

Explicitly solving D₀(t) proves difficult due to the integral on the r.h.s. of Eq. C.4. The time evolution of W_k(t), the modified spectrum of the MS waves, depends on wave-wave interaction and wave-particle interaction and is therefore very expensive to solve. Cassano & Brunetti (2005) use an approximation that simplifies the form of W_k to only depend on the injection spectrum of the waves I(k), which for this purpose is assumed to be a single power law, I(k) = I₀k^−a, and the most prominent dampening mechanism, the dampening by thermal electrons, Γ_th, e(k),

$$\begin{array}{c}\hfill W_k\simeq \frac{I(k)}{{\mathrm{\Gamma }}_{\mathrm{th},e}(k)},\end{array}$$(C.7)

with

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\mathrm{th},e}(k)& =\sqrt{32{\pi }^3}\rho \sqrt{m_{\mathrm{e}}{k}_{\mathrm{B}}T}{\left(\frac{v_{\mathrm{M}}}{B}\right)}^2\frac{W_k^B}{W_k}\mathcal{I}(x)k,\hfill \end{array}$$(C.8)

and

$$\begin{array}{c}\hfill \mathcal{I}(x)=2{\displaystyle \underset{1}{\overset{\infty }{\int }}}dx(\frac{1}{x}-\frac{1}{x^3})exp[-x^2{\left(\frac{v_M}{v_{\mathrm{th}}}\right)}^2].\end{array}$$(C.9)

Here v_M² is given by Eq. C.3 and $v_{\mathrm{th}}^2=\frac{2k_{B}T}{m_e}$ is the thermal velocity of the electrons (see Sect. 4.2 of Cassano & Brunetti 2005, for more details). To explicitly solve D₀ at every timestep, we express it as

$$\begin{array}{c}\hfill D_0(t)=\frac{4.45{\pi }^2}{c\sqrt{32{\pi }^3{m}_{\mathrm{e}}{k}_{\mathrm{B}}}}\frac{{\eta }_t{a}_k{E}_{\mathrm{turb}}}{n_{\mathrm{e}}{V}_{\mathrm{p}}\sqrt{T}\mathrm{\Delta }t\mathcal{I}(x)},\end{array}$$(C.10)

where η_t is the free parameter introduced by Cassano & Brunetti (2005) with η_t = 0.2 in this model, E_turb is the turbulent energy, V_p is the kernel volume of a particle and a_k is the turbulent scale factor

$$\begin{array}{c}\hfill a_k=\frac{k_{\max }^{-2/3}-k_{\min }^{-2/3}}{k_{\mathrm{mps}}^{-2/3}-k_{\mathrm{h}}^{-2/3}}.\end{array}$$(C.11)

The different components of this are k_min and k_max as the integral limits for Eq. C.4, the wavelength of the mean particle separation $k_{\mathrm{mps}}=N_{\mathrm{ngb},i}^{1/3}/(2h_i)$ and the maximum wavelength within a kernel of smoothing length h_i, k_h = 1/(2h_i). k_min and k_max are free parameters in this model and can be expressed as the inverse of the injection and damping scale of the MS waves which for the purpose of this work take the values

$$\begin{array}{cc}\hfill k_{\min }& ={({\mathrm{\Lambda }}_{\mathrm{inj}})}^{-1}={(200\mathrm{kpc})}^{-1}\hfill \\ \hfill k_{\max }& ={({\mathrm{\Lambda }}_{\mathrm{damp}})}^{-1}={(0.1\mathrm{kpc})}^{-1}.\hfill \end{array}$$

To estimate the turbulent velocity and from that the turbulent energy we use the same approach as Donnert & Brunetti (2014). We calculated E_turb from the root-mean-square (RMS) velocity which we can estimate by calculating the bulk velocity in the kernel; for example, for particle i,

$$\begin{array}{c}\hfill v_{\mathrm{bulk},i}={\displaystyle \sum _{j=0}^{N_{\mathrm{ngb}}}}\frac{m_j}{{\rho }_j}{\mathbf{v}}_{j}W(|{\mathbf{r}}_i-{\mathrm{r}}_{\mathrm{j}}|,h_i)\end{array}$$(C.12)

and subtracting that from the velocity of particle i

$$\begin{array}{c}\hfill {\mathbf{v}}_{\mathrm{rms},i}={\mathbf{v}}_i-v_{\mathrm{bulk},i}\end{array}$$(C.13)

such that

$$\begin{array}{c}\hfill E_{\mathrm{turb}}=\frac{1}{2}m{|{\mathbf{v}}_{\mathrm{rms}}|}^2.\end{array}$$(C.14)

We find that this provides a reasonable proxy for the turbulent energy captured in SPH, even when compared to an estimate based on a full Helmholtz decomposition (Groth et al., in prep). This approach does however use the total available turbulent energy to reaccelerate electrons. As was mentioned in Sect. 2.4, recent work by Brunetti & Lazarian (2016), Brunetti & Vazza (2020) suggests that only the solenoidal component of the turbulent velocity can efficiently reaccelerate CRs, therefore our current model will need to be revised to test the impact of these results on our simulations.

To evolve our electron spectra we then transform the problem into momentum space as for the energy losses and adiabatic changes.

With Eq. C.6 this leads to an upper integration boundary of

$$\begin{array}{c}\hfill p_u=p_i·e^{-2D_0\mathrm{\Delta }t}\end{array}$$(C.15)

and a spectral cut update of

$$\begin{array}{c}\hfill p_{\mathrm{cut},t+\mathrm{\Delta }t}=p_{\mathrm{cut},t}·e^{2D_0\mathrm{\Delta }t}\end{array}$$(C.16)

which we can then use to compute the fluxes between momentum bins according to Eq. B.11-B.12. Solving Eq. B.6 with Eq. C.6 gives

$$\begin{array}{c}\hfill R_{D_{\mathrm{pp}}}=-2D_0,\end{array}$$(C.17)

as all prefactors cancel out. We note that this approach is a simplification and neglects the stochastic component of the turbulent reacceleration. It does however counter the radiative loss processes of high-momentum electrons and with that provides the basis for studying the stabilization of radio haloes in galaxy clusters against the cooling times of their electrons. The exponential dependency on D₀ and the timestep Δt puts a significant timestep constraint on the simulation and needs to be handled with care. We therefore subcycle the solver if D₀Δt > 10⁻⁵. This increases the runtime by roughly a factor of 2-3, which would make it unfeasible in the case of the full cosmological box. We therefore switched this mode off in the box and ran comparison zoom-in simulations with the mode switched on.

Appendix D: Magnetic field models

In Fig. D.1 we show the scaling of our magnetic field models with electron density. We compute the scaling as the mean magnetic field strength per bin, obtained from the SLOW-CR3072³ simulation.

Fig. D.1. Scaling of our magnetic field models with electron density. For Bsim, Bβ and Bℱ we plot the mean magnetic field strength per density bin for the SLOW-CR30723 simulation.

Appendix E: Filament selection

In Fig. E.1, we sketch the selected region for the synchrotron spectra in Fig. 8. We note that this is a simplified illustration where we omit projection effects, as the filament does not lie perfectly in the xy plane of this projection.

Fig. E.1. In red we sketch the region selected for the filament spectra in Fig. 8. For the background image, we chose the COMA simulation with the Bdyn,↑ magnetic field model.

All Tables

All Figures

Fig. 1. Full-sky projections of the CR electron component. Upper panel: Energy density of electrons in the range E ∈ [5 − 500] MeV as the mean along the line of sight in a radius range r = 5 − 300 Mpc. This traces the injected energy density of the CR electrons with the longest lifetimes. Circles indicate the projected rvir of the cross-identified clusters (see Hernández-Martínez et al. 2024). Lower panel: CR electrons with energies ECR, e > 1 GeV integrated along the line of sight in a radius range r = 5 − 300 Mpc. The short cooling times of these electrons lead to this map tracing recent injection events.

In the text

Fig. 2. Full-sky projections of the mean magnetic field along the line of sight for the simulated magnetic field and post-processing models. Upper left panel: Simulated magnetic field obtained from the simulation. Upper right panel: Modeled magnetic field, assuming a plasma beta of $\beta \equiv \frac{P_{\mathrm{th}}}{P_{\mathrm{B}}}=50$. Center left panel: Modeled magnetic field based on a fraction ℱ of the turbulent pressure. Center right panel: Modeled magnetic field based on flux freezing B ∝ ρ2/3. Lower left panel: Modeled magnetic field, assuming a saturated turbulent dynamo scaling B ∝ ρ1/2 down to ne = 10−4 cm−3 and steeply declining beyond this density to match expectations for void magnetic fields. Lower right panel: Modeled magnetic field, assuming a saturated turbulent dynamo scaling for all density regimes.

In the text

	Fig. 3. (Intrinsic) synchrotron intensity at 144 MHz integrated along the line of sight between r = 10 − 300 Mpc. The different panels show the emission calculated based on the magnetic field configurations shown in Fig. 2.
In the text

Fig. 4. Cutout around our Coma replica with a side-length of d = 20 h−1 c Mpc. The circles represent r200 and r500. From left to right, we show surface density, absolute value of the simulated magnetic field, CR electron energy of electrons with energies above 1 GeV, and turbulent reacceleration coefficient. The upper panels show the cluster in the full cosmological box; the lower panels show it in the zoomed-in re-simulation.

In the text

Fig. 5. Magnetic field models applied to the COMA simulation. Top: Rotation measure (RM) mock observations of our Coma replica assuming different magnetic field models. Bottom: Radial profiles of the magnetic field. We show the profiles of all discussed simulations for the various magnetic field models in the solid colored lines. The dashed black lines show the observed profile by Bonafede et al. (2010).

In the text

Fig. 6. Same cutout as in Fig. 4, but showing the synchrotron surface brightness at 144 MHz, assuming a LOFAR beam size of θ = 60″ × 60″ as in Bonafede et al. (2022). The color bar is split at 1 μJy beam−1, which is a factor two above the sensitivity limit in the stacking approach by Hoang et al. (2023). From top to bottom, we show the different simulation runs, and from left to right we vary the magnetic field models. The lowest row shows the result of painting on synchrotron emission with the Hoeft & Brüggen (2007) model in a post-processing approach. For this, we use the COMA simulation due to the better performance of the shock finder at this resolution, compared to SLOW-CR30723.

In the text

Fig. 7. Resolved synchrotron spectrum maps accompanying Fig. 6. We constructed the synchrotron slope between the images of Fig. 6 and images with the same resolution at 1.4 GHz by fitting a single power law between the respective pixel values. We impose a surface brightness cutoff for the slope construction, similar to the one in Fig. 6 and also cut pixels with slopes α < −3. The lowest row shows histograms of the pixels in the spectral slope maps. Solid lines are simple histograms, dotted lines are synchrotron flux-weighted histograms. The grayed-out area indicates spectral slopes that are too flat for the standard DSA picture.

In the text

Fig. 8. Momentum and synchrotron spectra of CR electron populations in the filament segment indicated in Fig. E.1. Top panels: Sum of all momentum spectra for particles contained in the selected region. Bottom panels: Synchrotron spectra of the selected regions under the assumption of varying magnetic field models, as is indicated by the colors. In all cases, the simulated magnetic field strength is too low to produce synchrotron emission visible in the plotted range.

In the text

Fig. 9. Momentum- and synchrotron spectrum of an arbitrarily chosen SPH particle with an active CR electron population. We show the spectra at every second available output time between z = 1 and z = 0 in the COMA simulation. Left: Momentum spectrum of the electron population in the momentum range relevant for synchrotron emission. Dashed lines indicate the momentum bin boundaries imposed in our CR scheme. Right: The resulting synchrotron emissivity spectrum of the same electron population.

In the text

	Fig. 10. Density – temperature phase diagram of all gas particles in the simulation. The colored background shows the mass histogram in phase space. The contours show the CR electron energy of CRs with energies above 1 GeV contained per phase-space element.
In the text

Fig. 11. Individual panels showing the phase-space diagrams in the density-synchrotron emissivity phase-space, containing all gas particles in the simulation. The background shows mass histograms per phase-space element from the original simulation output. The colored lines indicate the mean synchrotron emissivity per density interval for the cosmological box simulation and the zoom-in simulations with and without turbulent reacceleration. All emissivities are calculated at a frequency of ν = 144 MHz. The upper limit of the synchrotron emissivity is taken from Hoang et al. (2023).

In the text

Fig. 12. Synchrotron emissivity of an electron spectrum with a power-law slope, q, extending between an injection momentum, ${\widehat{p}}_{\mathrm{inj}}$, and a maximum momentum of ${\widehat{p}}_{\max }={10}^6$. Each spectrum is normalized to contain a total energy density of ϵ = 1 erg cm−3 and emits in a magnetic field with B = 1 μG. The color shows emissivity and contour lines are spaced apart by half an order of magnitude.

In the text

	Fig. 13. Density – shock obliquity phase map. The colors indicate the 2D mass histogram for SLOW-CR30723, while the lines indicate the maximum of the histograms for SLOW-CR30723 and COMA as a function of density.
In the text

	Fig. 14. Histograms of magnetic field strength in all SPH particles containing a CR electron population in the cutout selected for Fig. 9. Colored lines correspond to the magnetic field models, the dashed black line shows the equivalence magnetic field of the CMB at z = 0.
In the text

	Fig. 15. To illustrate the impact of missing baselines on our simulated flux images, mock observations of COMA as in Fig. 6 under these considerations. Top: Relative change to the flux due to missing baselines in LOFAR. Bottom: Synchrotron flux under the consideration of instrument response from missing baselines.
In the text

	Fig. A.1. DSA parametrizations present in our model. Left: Sonic Mach number ℳs dependent acceleration models. Right: Shock obliquity (θB) dependent acceleration model.
In the text

	Fig. D.1. Scaling of our magnetic field models with electron density. For Bsim, Bβ and Bℱ we plot the mean magnetic field strength per density bin for the SLOW-CR30723 simulation.
In the text

	Fig. E.1. In red we sketch the region selected for the filament spectra in Fig. 8. For the background image, we chose the COMA simulation with the Bdyn,↑ magnetic field model.
In the text