© The Authors 2026

1. Introduction

Radio relics have presented a problem for theorists ever since Mills et al. (1960) reported an unusually strong spectral flux density of 57 μJy in the vicinity of Abell S753¹. Follow-up observations (see, in particular, Hoskins & Murdoch 1970; McAdam & Schilizzi 1977; Wall & Schilizzi 1979; Subrahmanyan et al. 2003) were able to progressively resolve this emission into two distinct regions: one roughly circular and centred on the cluster and one thin and arc-like, sited at the periphery of the cluster. Today, these are recognised as a radio halo and a radio relic², respectively.

In the 1970s, it was recognised that radio relics could not be fit by traditional models (Wall & Schilizzi 1979). Whilst their emission follows a power law, as is characteristic of synchrotron radiation, the spectral index α (where S ∝ ν^α, for flux density S at frequency ν) varies from the outer to the inner edge, typically starting at α ≲ −0.5 and steepening strongly thereafter (Feretti et al. 2012). Inactive radio galaxies were initially favoured as sources (see e.g. Goss et al. 1987; Harris et al. 1993), however, by the end of the 1990s, it was clear that not every relic could be spatially associated with one (Feretti & Giovannini 1996). Moreover, relics show no apparent spectral cutoff (Komissarov & Gubanov 1994; Rajpurohit et al. 2020b), which, combined with the relatively short cooling times of megahertz- and gigahertz-emitting electrons, implies an in-situ source (Enßlin & Brüggen 2002; Kang et al. 2012).

It was eventually recognised that relics actually trace cosmological shocks (Enßlin et al. 1998). In particular, they predominantly trace low Mach number shocks (ℳ ≲ 3 − 5) driven by cluster mergers (see e.g. Hoeft & Brüggen 2007; Markevitch & Vikhlinin 2007). The observation of such shocks in the X-ray band is now well established, as is their spatial association with radio relics (see e.g. Brunetti & Jones 2014; van Weeren et al. 2019), and references therein). This link has been strengthened by observed correlations between the orientation of the radio relic and the merger axis (van Weeren et al. 2011b), and scalings of the radio relic power with both the X-ray luminosity (Bonafede et al. 2012) and mass of the host cluster (de Gasperin et al. 2014).

Two main mechanisms have been proposed to link such shocks with the resultant highly polarised, megaparsec-sized emission. However, whilst adiabatic compression scenarios have not been wholly ruled out³ (Enßlin & Gopal-Krishna 2001; Button & Marchegiani 2020), it is questionable whether they can replicate the observed scaling relations (Colafrancesco et al. 2017) and spectral slopes (see e.g. evidence in van Weeren et al. 2017). Instead, the prevailing theory is that the merger shocks (re-)accelerate existing populations of electrons at the cluster outskirts via diffusive shock acceleration (DSA; Fermi 1949; Drury 1983; Blandford & Eichler 1987). This theory is not without problems either, however, as such shocks are expected to have very low acceleration efficiencies (≲1%) (Kang & Jones 2005; Kang & Ryu 2013; Mou et al. 2023). Acceleration from the thermal pool is consequently incompatible with the observed spectral intensities (Pinzke et al. 2013; Vazza & Brüggen 2014; Botteon et al. 2020a). To solve this issue, a population of pre-existing, semi-relativistic ‘fossil’ electrons are usually invoked (Markevitch et al. 2005; Kang & Ryu 2011; Pinzke et al. 2013; Vazza et al. 2015; Botteon et al. 2020b), although some competing theories have emerged, such as turbulent re-acceleration (Fujita et al. 2015) and the multi-shock scenario (Inchingolo et al. 2022).

Whilst first-order Fermi re-acceleration has emerged as the standard paradigm for explaining radio relics, our understanding of their origin remains far from complete. In particular, we highlight the following seven major problems:

1. What is the origin of the seed electrons needed for re-acceleration? Four competing schools of thought currently exist: (i) seed electrons are accelerated via diffusive shock acceleration during structure formation events and accretion shocks and therefore have an external origin (Pinzke et al. 2013; Vazza et al. 2023), (ii) electrons are injected by radio galaxies at the extremities of a cluster (Enßlin et al. 1998; de Gasperin et al. 2015; Johnston-Hollitt 2017; Botteon et al. 2020b; Vazza et al. 2023), (iii) seed electrons originate from active galactic nuclei (AGN) lobes emitted from the cluster core, (Enßlin & Gopal-Krishna 2001; Shulevski et al. 2015; Vazza et al. 2021; ZuHone et al. 2021), and (iv) the electron population is common with that of radio haloes, and hence turbulent re-acceleration is the explanation (Beduzzi et al. 2024). Whilst these options are not necessarily mutually exclusive, they each suffer from issues. For example, respectively, (i) has yet to be shown to work in a fully cosmological MHD simulation, where the Fokker-Planck equation is solved explicitly⁴, (ii) radio galaxies are not always observationally evident, (iii) it is unclear how AGN lobes survive thermal instabilities during buoyancy, and (iv) it is unclear whether radio haloes can extend as far as some radio relics.

2. What is the origin of relic morphology? Radio relics exhibit a wide range of morphologies including irregular forms, such as in the case of Abell 2256 (van Weeren et al. 2012b), more regular, arc-like forms, such as the ‘Sausage’ relic (Kocevski et al. 2007; van Weeren et al. 2010), and even relatively linear forms, such as in the case of the ‘Toothbrush’ relic (van Weeren et al. 2012a). Moreover, with increasing resolution, it is clear that most relics have a filamentary nature (Rajpurohit et al. 2020a, 2022). These are often attributed to magnetic fields (Rudnick et al. 2022) but shock geometry no doubt also plays a role. Numerical studies have already made some inroads towards answering this question, with simulations of relic morphology (Wittor et al. 2019, 2021b; Lee et al. 2024; Nuza et al. 2024), investigations into substructure and relic ‘patchiness’ (Domínguez-Fernández et al. 2021a, 2024), and a demonstration of a potential formation process for so-called ‘wrong way’ relics (Böss et al. 2023) all being published recently. However, we are still far from being able to explain the morphology of individual relics.

3. Mach numbers inferred from radio measurements are typically larger than those inferred from X-ray data. How do we explain this discrepancy? Temperature and density jumps can be measured using X-ray observations, allowing Mach numbers to be inferred through the standard jump conditions. Meanwhile, the spectral slope of the radio emission can also be used to infer Mach numbers, assuming DSA. However, while in principle they measure the same shock, the two measurements rarely agree (see Wittor et al. 2021a; Lee et al. 2025, and references therein). Potential explanations include the fact that X-ray measurements are subject to projection effects (Wittor et al. 2021a) and that radio emission may be biased to higher values as a result of a Mach number distribution (Hoeft et al. 2011; Skillman et al. 2013; Hong et al. 2015; Roh et al. 2019; Wittor et al. 2019; Domínguez-Fernández et al. 2021a; Wittor et al. 2021a). However, a simulation that shows a full ab-initio development of this effect is missing from the literature.

4. Observations appear to imply μG magnetic fields in radio relics. How is this possible given that the surrounding intracluster medium (ICM) has strengths at least an order of magnitude lower? Magnetic field estimates have been made in relics based on inverse Compton emission studies (Chen et al. 2008), constraints provided by the width of the radio relic (Markevitch et al. 2005; Hoeft & Brüggen 2007; van Weeren et al. 2010), and constraints using Faraday rotation (Bonafede et al. 2013). All of these models suggest magnetic field strengths on the order of μG, albeit with some scatter in their estimates. Observations (Brunetti et al. 2001; Govoni et al. 2017) and simulations (see e.g. Skillman et al. 2013; Domínguez-Fernández et al. 2019; Nelson et al. 2024) of the surrounding ICM, meanwhile, typically imply magnetic fields that are only a fraction of this. This cannot be covered by standard shock compression alone (Donnert et al. 2016), although there has been some evidence that radio observations bias values slightly high (Wittor et al. 2019).

5. Why are standard cooling models unable to match spectral variations? Four cooling models are typically applied when trying to analyse radio relics (see e.g. van Weeren et al. 2012a; Rajpurohit et al. 2020a). The most simple of these assumes a one-time injection⁵ with subsequent cooling that is dependent on (‘KP’ Kardashev 1962; Pacholczyk 1970) or independent of (‘JP’ Jaffe & Perola 1973) the cosmic ray (CR) electron pitch angle, respectively. Extensions include extended injection until the present time (‘CI’ Pacholczyk 1970) or until a fixed time in the past (‘KGJP’ Komissarov & Gubanov 1994). Each of these mechanisms generates a characteristic spectral shape, the evolution of which can be probed using ‘colour-colour’ diagrams (Katz-Stone et al. 1993). These are created by correlating spectral index maps, thereby producing a phase space diagram with α_ν1^ν₃(x, y) as a function of α_ν1^ν₂(x, y), where ν₁ < ν₂ < ν₃, and x and y are spatial coordinates. Spectral features then translate to distinctive trajectories in this plane. However, despite fine-tuning, models are typically unable to replicate the tracks produced by observed radio relics.

6. Recent studies imply that CR electron acceleration is only efficient above a critical Mach number of ℳ_crit ≈ 2.3. How do we reconcile this with reports of radio relics at shocks apparently weaker than this? There are now several claims in the literature that suggest that efficient CR electron acceleration can only take place above ℳ_crit  ≈  2.3. This is typically because shocks below this value do not excite sufficiently strong plasma waves via CR electron-driven instabilities⁶. These waves are necessary in order to confine the CR electrons during their acceleration to radio-emitting energies. This result has been observed in particle-in-cell (PIC) simulations at the pre-acceleration phase of CR electrons in both quasi-parallel (Shalaby et al. 2022; Gupta et al. 2024) and quasi-perpendicular shocks (Kang et al. 2019; Ha et al. 2021; Boula et al. 2024). However, efficient acceleration may also fail at quasi-perpendicular shocks due to a lack of CR proton generated plasma waves (Ha et al. 2018; Ryu et al. 2019) or due to conditions on the firehose instability (Guo et al. 2014). Whilst shocks derived from radio data are generally consistent with the above studies, shocks derived from X-rays are highly inconsistent. Indeed, there are now a substantial number of relics reported with X-ray-derived Mach numbers below ℳ = 2.3 (see Wittor et al. 2021a, and references therein). How can we explain this apparent inconsistency?

7. If electrons are most efficiently accelerated at quasi-parallel shocks, why does polarisation data appear to imply quasi-perpendicular acceleration? Both observations (Masters et al. 2013; Liu et al. 2019) and simulations (Crumley et al. 2019; Winner et al. 2020; Shalaby et al. 2021, 2022) imply that DSA is most effective for electrons at quasi-parallel shocks⁷ (see theory introduced in Bykov & Uvarov 1999). This, it seems, is in contradiction with observations of radio relics; for example, observations of the ‘Sausage’ relic appear to show the magnetic field vectors aligning almost perfectly with the shock surface (van Weeren et al. 2010; Di Gennaro et al. 2021). This behaviour was replicated in part by Wittor et al. (2019) and Domínguez-Fernández et al. (2021b), however neither included magnetic-obliquity-dependent shock acceleration in their simulations. It remains to be seen whether these results can still be replicated in such a scenario.

This paper marks the first in a series where we attempt to tackle the problems listed above. In this study, we tackle problems 3–5. Our strategy is to first run cosmological simulations in order to identify typical shock conditions, and then use these to inform significantly higher-resolution shock tube simulations. In addition, we implement upstream density turbulence, as inferred from both observations (Simionescu et al. 2011; Eckert et al. 2015; Ghirardini et al. 2018) and previous simulations (Nagai & Lau 2011; Zhuravleva et al. 2013; Battaglia et al. 2015; Angelinelli et al. 2021). By employing this method, we combine the strengths of previous numerical approaches, whilst also resolving relevant physics that is, as yet, out of the reach of current-generation cosmological simulations.

The paper is organised as follows: in Sect. 2, we introduce both the codes and the simulation setup used in this study. In Sect. 3, we present an analysis of a cosmological galaxy cluster merger simulation. This is then used to inform a series of shock tube simulations, the results of which we present in Sect. 4. In this section, we show that upstream density turbulence results in: shock corrugation and the generation of downstream velocity turbulence (Sect. 4.1), the formation of a Mach number distribution at the shock front (Sect. 4.2), flatter CR electron spectra (Sect. 4.3), breaking of the assumption that distance from the shock front is a reliable indicator of time cooled for individual electrons (Sect. 4.4), an amplification of the magnetic field to μG levels (Sect. 4.5), and synchrotron emission that is better able to replicate observed colour-colour diagrams (Sect. 4.6). We finish the paper with a discussion of the caveats of our model (Sect. 5) and a summary of our main conclusions (Sect. 6).

2. Methodology

2.1. AREPO

All simulations presented in this paper were run with the moving-mesh code AREPO (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020). This code employs a set of unfixed mesh-generating points to construct a volume-filling Voronoi tessellation. Upon this, the equations for an ideal magnetohydrodynamic (MHD) fluid are solved using a second-order finite-volume Godunov scheme (Pakmor et al. 2011; Pakmor & Springel 2013) based on an HLLD Riemann solver (Miyoshi & Kusano 2005).

This method inherits the advantages of both Eulerian and Lagrangian codes. For example, cells are kept within a factor of two of a predetermined target mass (see following sections for values), being refined and merged as necessary. This means that spatial resolution increases naturally with structural complexity, allowing computational power to be concentrated where the system is most dynamic. Meanwhile, mesh-generating points follow the motion of the gas, leading to quasi-Lagrangian behaviour, which has the effect of substantially reducing numerical diffusion (see e.g. Pfrommer et al. 2022). The result is increased accuracy at the same resolution compared to competing codes; in particular, AREPO outperforms standard smoothed-particle hydrodynamics (SPH) (see comparison studies in Vogelsberger et al. 2012; Sijacki et al. 2012; Kereš et al. 2012; Bauer & Springel 2012).

Especially germane to this study is the impact this has on the resolution of the turbulent cascade (Kolmogorov 1941). Several studies have shown that AREPO can reproduce this (see e.g. Bauer & Springel 2012; Whittingham et al. 2021), with strong evidence that it can also reproduce the closely related growth of the magnetic coherence scale in a fluctuating small-scale dynamo, given sufficiently high resolution⁸ (Pakmor et al. 2017; Whittingham et al. 2021, 2023; Pfrommer et al. 2022; Pakmor et al. 2024). The same studies have shown that the employed Powell 8-wave divergence cleaner (Powell et al. 1999) is highly robust even in strongly dynamic flows.

2.2. Cosmological simulations

To make sure that our shock tubes probe the most relevant physics, we first analysed the formation of shocks in a series of cluster mergers. We used the PICO-Clusters (Plasmas In COsmological Clusters) suite (Berlok et al., in prep.) for this purpose, which comprises of 25 individual zoom-in simulations. The simulated clusters vary in mass⁹ between 10^14.9 − 10^15.5 M_⊙ at z = 0. Cosmological volumes are periodic with a side-length of 1 co-moving Gpc h⁻¹ and use a Planck 2018 cosmology (Planck Collaboration VI 2020); i.e. Hubble’s constant is H₀ = 100 h km s⁻¹ Mpc⁻¹ = 67.3 km s⁻¹ Mpc⁻¹ with the density parameters for matter, baryons, and a cosmological constant being Ω_m = 0.316, Ω_b = 0.049, and Ω_Λ = 0.684, respectively.

Each zoom-in simulation was evolved from z = 127 and has a dark matter resolution of 5.9 × 10⁷ M_⊙ (target gas mass of 1.1 × 10⁷ M_⊙) in the highest resolution region, which extends without contamination of lower resolution particles to a minimum of 3 × R₂₀₀ throughout the simulation. This resolution is roughly comparable to that of the TNG300 (Nelson et al. 2019), TNG-Cluster (Nelson et al. 2024), and MilleniumTNG (Pakmor et al. 2023) simulations. Like these simulations we also used the IllustrisTNG galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018b). This includes a Springel & Hernquist (2003) ISM, stellar and AGN feedback (Springel et al. 2005; Weinberger et al. 2017), radiative cooling and metal enrichment (Wiersma et al. 2009; Vogelsberger et al. 2013), magnetic fields (Pakmor & Springel 2013; Marinacci et al. 2018), and the Schaal & Springel (2015) shock finder.

Cosmic ray protons and electrons were not included in this model. This is not a problem, however, as we focus at this stage purely on the study of shocks in low-redshift massive mergers; i.e. those most expected to create radio relics. Such shocks are only minimally affected by CR protons in clusters (Pfrommer et al. 2017) and are negligibly affected by CR electrons, owing to their small energy content (Winner et al. 2019).

The IllustrisTNG model has been shown to produce clusters that compare well with observations (Springel et al. 2018; Pillepich et al. 2018a). We note, in particular, that density and pressure profiles compare favourably (Pakmor et al. 2023). This is important for this study, as it is these properties specifically that we probed in order to create higher resolution shock tube simulations.

2.3. Shock tube simulations

2.3.1. Resolution

We used the up- and downstream gas properties extracted from the aforementioned cosmological simulations to create initial conditions for a series of shock tubes. These shock tubes are three-dimensional with a periodic volume of 1800 × 300 × 300 kpc³. The volume can be further separated into four distinct regions, each 450 kpc long:

− Region I: Low-resolution piston region

− Region II: Piston region with increasing resolution

− Region III: High-resolution upstream region

− Region IV: Upstream region with decreasing resolution

We illustrate these regions along with the initial conditions for a Mach 3 shock in dashed lines in Fig. 1. Regions I and IV are necessary as our shock finder necessitates periodic boundary conditions. We must therefore allow for the reverse shock, as seen in Fig. 1 in region IV. To mitigate for this, we varied the resolution across the volume, so that cells in region I have side-length 75 kpc, whilst cells in region III have a side-length of 2.5 kpc. By ramping the resolution up and down in this fashion, we aimed to focus our computational resources on region III, whilst also preventing the accumulation of truncation errors.

Fig. 1. Top: Mean pressure along the x axis at t = 0 Myr (dashed) and t = 100 Myr (solid) in our Mach 3 ‘Flat’ simulation (see Table 1). Dotted lines indicate the edges of various regions in the shock tube (see Sect. 2.3.1). We greyed out panels that we do not analyse. Bottom: As above, but lines show the mean density. The initial pressure discontinuity causes a shock to propagate into region III. As a result of shock compression, a dense sheet moves rightwards, bounded by the density jumps at the contact discontinuity and the shock front. This feature is critical to the mechanism we analyse in this paper.

We set the target gas mass resolution in the high-resolution region of our shock tubes to m_gas ≈ 1.5 × 10⁴ M_⊙, with the exact value based on our initial density (see following section). We allowed the cells here to refine during the course of the simulation, which provides sub-kpc resolution after shock-compression. We fixed the resolution outside of region III, which leads to some minor variations from the expected density and pressure profiles to be evident in Fig. 1. These variations do not impact the accuracy of the dynamics in region III.

2.3.2. Density and pressure values

As we show in Sect. 3, at the outskirts of clusters, merger shocks are able to generate dense, shock-compressed sheets. Indeed, this forms a critical part of the model for radio relics we introduce in this paper. This scenario can be replicated by setting the mean density of all four regions to be the same and driving the shock purely by setting a jump in pressure. Indeed, this turns out to be a good approximation to the scenario observed in our cosmological simulations. Following these simulations, we set the density in our shock tube to n_e = 3.5 × 10⁻⁵ cm⁻³ (or, equivalently, ρ ≈ 6.7 × 10⁻²⁹ g cm⁻³) and the initial upstream pressure to P₁ = 1 × 10⁻¹³ dyne cm⁻². Regions I and II then form the high-pressure piston, whilst regions III and IV form the low-pressure upstream region into which the shock runs. These initial conditions are shown as dashed lines in Fig. 1.

The Mach number of the shock is determined entirely in our case by the choice of pressure in the piston. Note, this differs from the classic Sod (1978) shock tube setup, in which density also plays a role. We chose to investigate the impact of ℳ = 2 and ℳ = 3 shocks, as these are typical values reported in radio relics (see e.g. van Weeren et al. 2017). This fixed the initial downstream pressure to P₂ = 1 × 10⁻¹² dyne cm⁻² and P₂ = 2.32 × 10⁻¹² dyne cm⁻², respectively, with exact values being calculated using a Riemann solver. This setup generates a narrow, shock-compressed region as desired (see solid lines in Fig. 1).

By initialising the shock in this way, we may test a variety of upstream conditions. Indeed, we examine a full range of variables in a companion paper (Whittingham et al., in prep.). In this current paper, however, we restrict ourselves simply to the impact of including density and magnetic turbulence at all.

2.3.3. Simulation variations

To unpin the impact of magnetic and density turbulence we ran four separate simulations. These are listed in Table 1 and cover the various possible permutations. Density turbulence was added only to region III, with the other regions given a constant density (see Sect. 2.3.2 for values). A short buffer region without density turbulence was also provided at the start of region III to allow for the initial formation of a planar shock.

In simulations with magnetic turbulence, on the other hand, this was added to all regions in order to keep magnetic divergence to a minimum. In simulations without magnetic turbulence, the magnetic field was oriented in the x-direction and given a fixed strength equal to the root-mean-square (RMS) value in the turbulent runs. We did not correlate the magnetic field and density fluctuations in this initial study, leaving a study of this improvement to future work.

2.3.4. Generating turbulence

We generated turbulence using the method given in Ruszkowski et al. (2007) and Ehlert et al. (2018). In this method, turbulence is first created in k-space before being converted to configuration space using a fast Fourier transform. For density turbulence, we applied Kolmogorov (1941) turbulence (P(k)∝k^−5/3) below an injection scale of 150 kpc, with white noise imposed above this. This is similar in magnitude to the scale used in previous shock tube studies of relics (see e.g. Domínguez-Fernández et al. 2021a). We additionally implemented log-normal variance, as inferred from both simulations and observations (see e.g. Kawahara et al. 2008), using a Box-Muller random variate method. Following Zhuravleva et al. (2013), we set the relative variance of the distribution to σ_ρ/μ_ρ = 0.4, where σ_ρ is the standard deviation and μ_ρ is the mean of the density distribution, respectively.

We expect the peak of the magnetic power spectra to lag behind that of the kinetic power spectra by a factor of a few (see, e.g. Tevlin et al. 2025), and so chose an injection scale here of 40 kpc. As before, we applied Kolmogorov turbulence below this scale and white noise above it. We set the RMS strength, such that in the upstream ⟨P_th/P_B⟩ = 100, where P_th and P_B are the thermal and magnetic pressures, respectively. We find this is typical of values in the ICM value in our cosmological simulations (Tevlin et al. 2025). The result is an RMS field strength of approximately 0.16 μG. Following, Ehlert et al. (2018), we also applied Gaussian variance, with a standard variation, σ_B, determined by Parseval’s formula:

$$\begin{array}{c}\hfill {\sigma }_B^2=\frac{1}{N^2}{\displaystyle \sum _{\mathit{k}}}{|{\mathit{B}}_{i,\mathit{k}}|}^2,\end{array}$$(1)

where i covers the three spatial components, N is the number of cells in the initial conditions, and the sum runs over all k-values. Magnetic divergence was cleaned by projecting it out in k-space, i.e. by subtracting $\widehat{\mathit{k}}(\widehat{\mathit{k}}\mathbf{·}{\mathit{B}}_{\mathit{k}})$.

Note that our method of creating turbulence differs fundamentally from the method used in the shock tubes of Domínguez-Fernández et al. (2021a,b, 2024), who first drive turbulence, before letting it then decay over time. In the shock tube simulations presented in this study, all cells are initially at rest; i.e. there is no initial velocity turbulence. We discuss the advantages and disadvantages of our method in Sect. 5.

2.3.5. CR physics and tracer particles

We used the standard AREPO code (see Sect. 2.1) in our shock tubes, augmented with the Pfrommer et al. (2017) CR proton module, in which CR protons are treated as a relativistic fluid with effective adiabatic index, γ_a = 4/3. We injected CR protons at shocks according to DSA theory (see Sect. 2.4) and advected them thereafter with the gas. We did not account for streaming and diffusion, as strong turbulence in the vicinity of the shock is believed to decrease the diffusion coefficient to values approaching Bohm diffusion (Caprioli & Spitkovsky 2014). This leaves advection the dominant transport process. Cosmic ray protons are dynamically unimportant in our simulations, but are intrinsically linked to the modelling of CR electrons (see Sect. 2.5).

As AREPO is only a quasi-Lagrangian code, we also employed the use of tracer particles (Genel et al. 2013) for CR electron modelling. One tracer particle was placed in each upstream cell, with additional tracers placed at the same resolution at the end of region II, just behind the initial pressure discontinuity. This allows us to take advantage of the contact discontinuity that forms (see Fig. 1), helping us to define the volume of cells in the injected region. In total, this results in approximately 2.9 × 10⁶ tracer particles. The tracer particle functionality has been adapted to save values on-the-fly in order to be used by CREST (Winner et al. 2019). Of particular importance is the recording of up- and downstream shock properties. For this, we used a numerical shock finder.

2.4. Shock finder

In both cosmological and shock tube simulations, we used the Schaal & Springel (2015) shock finder, with extensions for CRs, as presented in Pfrommer et al. (2017). This is based on the Rankine-Hugoniot jump conditions with modifications for the case of additional CR pressure and injection. These produce the formula

$$\begin{array}{c}\hfill {\mathcal{M}}^2=(\frac{P_2}{P_1}-1)\frac{x_{\mathrm{s}}}{{\gamma }_{\mathrm{a},\mathrm{eff}}(x_{\mathrm{s}}-1)},\end{array}$$(2)

where γ_a, eff is the effective adiabatic index, P₁ and P₂ denote the total pre- and post-shock pressures, and x_s is the density jump at the shock.

The dissipated energy flux at the shock can be expressed as

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{diss}}={\epsilon }_{\mathrm{diss}}{A}_{\mathrm{shock}}\frac{\mathcal{M}{c}_{\mathrm{s},1}}{x_{\mathrm{s}}},\end{array}$$(3)

where ε_diss is the post-shock energy density minus the adiabatically compressed pre-shock energy density, A_shock is the cell’s shock surface, and c_s, 1 is the upstream sound speed (see Pfrommer et al. 2017, for details).

We guarded against spurious shocks by further ensuring the following conditions:

• (i)

  $\mathbf{\nabla }\mathbf{·}\upsilon <0$

• (ii)

  ∇T⋅∇ρ > 0

• (iii)

  ℳ > ℳ_min

• (iv)

  x_s > x_s,  min,

where υ, T, and ρ, are the gas velocity, temperature, and density respectively. The guards act to: (i) ensure converging velocity flows, (ii) filter against tangential and contact discontinuities, (iii) ensure a minimum Mach number (here, ℳ_min = 1.3), and (iv) ensure a minimum density jump, respectively. This last criterion is especially important for ensuring the stability of Eq. (2) at weak shocks and is calculated using the hydrodynamic jump condition with ℳ_min.

Shocks are numerically broadened in our simulations; typically with a thickness of two to three cells. The cell with the minimum velocity divergence (i.e. maximum compression) is labelled the ‘shock surface’ cell, whilst the cell directly outside the shock zone is labelled the ‘post-shock’ cell. We ran the shock finder in the before-snapshot mode in the cosmological simulation and in the on-the-fly mode in our shock tubes.

2.5. CREST

Cosmic ray electron modelling in this study was performed using the post-processing code CREST (Winner et al. 2019). This code uses the previously discussed tracer functionality in AREPO to store relevant data on the MHD timestep, allowing us to evolve the Fokker-Planck equation in the Lagrangian frame. In particular, we accounted for adiabatic changes; cooling via Coulomb, bremsstrahlung, inverse Compton, and synchrotron losses; and DSA with magnetic obliquity (see formulae in Pais et al. 2018; Winner et al. 2020). In principle, CREST is also capable of Fermi I re-acceleration, Fermi II momentum diffusion, and Bell amplification, but we did not use these features in the following work.

On each tracer particle, CREST samples the underlying CR electron spectral density field. In practice, this means that all points closest to a given tracer are represented by that tracer’s spectrum. To compute a thermodynamically extensive quantity such as the CR electron energy, we therefore constructed the volume associated with a given tracer by using a Voronoi tessellation of space with the tracer position as the mesh-generating point. This gives the tracer an evolving volume, V_cell, where E_e = ε_eV_cell is the CR electron energy in that volume, and ε_e is the CR electron energy density obtained by taking the second moment of the CR electron spectrum.

Tracers were initially assigned a purely thermal spectrum. This is acceptable, as whilst a pre-existing non-relativistic population is likely required to reach radio relic luminosities (see Sect. 1), the addition of this feature should only affect the normalisation of the spectra, not its shape at radio-emitting frequencies (see e.g. Pinzke et al. 2013; Winner et al. 2019). Modelling the impact of such a population is outside the scope of this work, and hence left to a future study.

For our cooling modules, we assumed that the gas is primordial and hence that the mass fraction of hydrogen is X_H = 0.76. We also assumed that the gas is fully ionised, resulting in a mean molecular weight of μ = 0.588 and an ionisation fraction of x_e = 1.157. To aid comparison with observations we assumed a redshift of z = 0.2, which is the approximate redshift of the Toothbrush and Sausage radio relics (van Weeren et al. 2010, 2012a). The energy density of the cosmic microwave background (CMB) was therefore set to ϵ_CMB = 8.65 × 10⁻¹³ erg cm⁻³, or equivalently a magnetic field strength of B_CMB ≈ 4.7 μG.

Throughout this work, we use dimensionless momentum, $p=\stackrel{\sim }{p}/(m_{\mathrm{e}}c)$, where $\stackrel{\sim }{p}$ is the dimensional momentum, m_e is the electron rest mass, and c is the speed of light. We also convert the transport equation (see Winner et al. 2019) into one-dimensional momentum space, such that f^1D = 4πp²f^3D, where f^1D and f^3D are the one- and three-dimensional distribution functions, respectively. We let p range in all simulations from 10⁻¹ to 10⁸ and used 10 logarithmically spaced bins per decade. We find this is sufficient to produce converged spectra. In the remainder of this section, we summarise the DSA implementation, having the greatest impact on our study. We direct the reader, however, to Winner et al. (2019) for a more thorough discussion of all the physics represented in the code.

Cosmic ray electrons are accelerated in CREST in shock surface and post-shock cells by attaching to the one-dimensional distribution function a power-law slope¹⁰ with

$$\begin{array}{c}\hfill \alpha =-\frac{x_{\mathrm{s}}+2}{x_{\mathrm{s}}-1}.\end{array}$$(4)

This slope was limited to a maximum of –2.2, following Caprioli et al. (2020). This has a very minor impact on our results, however, given the weak shocks we simulate and the correspondingly steep spectral slopes. To avoid spurious heating before the shock due to numerical broadening, we stopped updating the recorded density on encountering a shocked cell. When the tracer reaches a shock surface or post-shock cell, the density was then updated to its post-shock value. This produces a discrete jump at the shock.

Following Webb et al. (1984) and Pinzke & Pfrommer (2010), we set the maximum momentum up to which we accelerate to

$$\begin{array}{c}\hfill p_{\max }=\frac{{\upsilon }_{\mathrm{post}}}{c}\sqrt{\frac{6\pi e}{{\sigma }_{\mathrm{T}}}\frac{B}{(B^2+B_{\mathrm{CMB}}^2)}\frac{x_{\mathrm{s}}(x_{\mathrm{s}}-1)}{(x_{\mathrm{s}}+1)}},\end{array}$$(5)

where υ_post is the post-shock velocity in the shock rest-frame, e is the elementary charge, σ_T is the Thomson cross section, and B is the strength of the magnetic field.

We then required that p_min, the minimum acceleration momentum, fulfils

$$\begin{array}{c}\hfill {\displaystyle \underset{p_{\min }}{\overset{p_{\max }}{\int }}}{f}_{\mathrm{th}}(p_{\min })p_{\min }^{-\alpha }{p}^{\alpha }{E}_{\mathrm{e},\mathrm{kin}}(p)\mathrm{d}p=\mathrm{\Delta }{\epsilon }_{\mathrm{e}},\end{array}$$(6)

where f_th(p) is the thermal Maxwellian¹¹, $E_{\mathrm{e},\mathrm{kin}}(p)=(\sqrt{1+p^2}-1)m_{\mathrm{e}}{c}^2$ is the electron kinetic energy, and Δε_e is the increase in energy density of the CR electron population due to acceleration. We used 0.1% of the liberated thermal energy at the shock, converting this to an energy density by dividing by the cell volume (see Winner et al. 2019).

The normalisation of the injected spectrum is subsequently

$$\begin{array}{c}\hfill C=f_{\mathrm{th}}(p_{\min })p_{\min }^{-\alpha }.\end{array}$$(7)

At the high momentum end of the spectrum, we imposed a super-exponential cutoff (see Enßlin et al. 1998; Zirakashvili & Aharonian 2007; Pinzke & Pfrommer 2010), which modifies the spectrum thusly:

$$\begin{array}{c}\hfill \stackrel{\sim }{f}(p)=f(p){[1+0.66{\left(\frac{p}{p_{\max }}\right)}^{2.5}]}^{1.8}exp[-{\left(\frac{p}{p_{\max }}\right)}^2].\end{array}$$(8)

As explained in Sect. 1, there are now several studies that suggest that CR electron acceleration can only take place above a critical Mach number of ℳ_crit = 2.3. Subsequently, we did not accelerate electrons below this Mach number in the simulations that include upstream density turbulence¹². We show in Whittingham et al. (in prep.) that this has remarkably little impact on the subsequent emission, even when the Mach number distribution peaks at ℳ = 2. This is due to the fact that emission is dominated by the tail of the Mach number distribution.

2.6. CRAYON+

In the final step of our methodology, we post-processed the resultant CREST output with the CRAYON+ code (Werhahn et al. 2021b), which is able to convert spectra into instantaneous emission for a range of CR electron processes. Naturally, we are most interested in the radio synchrotron emission. We summarise the salient formulae here, but encourage the interested reader to see Section 2.4 and associated appendices of Werhahn et al. (2021a) for a more comprehensive overview.

Following Rybicki & Lightman (1986), the synchrotron emissivity j_ν of a tracer at frequency ν is

$$\begin{array}{c}\hfill j_{\nu }=\frac{\sqrt{3}{e}^3{B}_{\perp }}{m_{\mathrm{e}}{c}^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}}f(p)F(\nu /{\nu }_c)\mathrm{d}p\propto B_{\perp }^{1-{\alpha }_{\nu }}{\nu }^{{\alpha }_{\nu }},\end{array}$$(9)

where B_⊥ is the projection of the magnetic field onto the plane perpendicular to the line of sight, α_ν is the radio spectral index, and F(ν/ν_c) is the dimensionless synchrotron kernel, with ν_c the critical frequency. These last two variables are, in turn, defined as

$$\begin{array}{c}\hfill {\nu }_{\mathrm{c}}=\frac{3e{B}_{\perp }{\gamma }^2}{4\pi m_{\mathrm{e}}c},\end{array}$$(10)

where $\gamma =\sqrt{1+p^2}$ is the Lorentz factor and

$$\begin{array}{c}\hfill F(x)=x{\displaystyle \underset{x}{\overset{\infty }{\int }}}{K}_{5/3}(\xi )\mathrm{d}\xi ,\end{array}$$(11)

where K_5/3 is the modified Bessel function of order 5/3, with x = ν/ν_c. To get the right-hand side of Eq. (9), we have assumed a power-law spectrum over an appropriate range of frequencies.

The specific radio synchrotron intensity, I_ν, at frequency ν is then obtained by integrating j_ν (with units of erg s⁻¹ Hz⁻¹ cm⁻³) along the line of sight, s:

$$\begin{array}{c}\hfill I_{\nu }=\frac{1}{4\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}}{j}_{\nu }ds.\end{array}$$(12)

This gives I_ν in units¹³ of erg cm⁻² s⁻¹ Hz⁻¹ sterad⁻¹. In the remainder of the paper, where projections are shown we have set the upper limit of the integral in Eq. (12) to be the depth of the box, equal to 300 kpc.

3. Cosmological analysis

As explained in Sects. 1 and 2, the first step in our strategy is to analyse shock conditions in cosmological simulations. To this end, we present here a case study. We use a variation¹⁴ of Halo 0003 from the PICO-Cluster suite (Berlok et al., in prep) as it features a particularly energetic merger, making it especially likely that the generated shock wave will form a radio relic (see similar shocks in e.g. Böss et al. 2023; Lee et al. 2024). However, our analysis generalises to all the cluster mergers we investigated.

The major merger in Halo 0003 takes place between z ≈ 0.3 and z ≈ 0.1. At the start, the progenitors have M₂₀₀ masses of 10^14.8 M_⊙ and 10^15.2 M_⊙, giving a mass ratio of roughly 1:2.5. We show gas conditions during the merger in Fig. 2, focusing on a period when the energy dissipation in the shock is at its most intense. At this time, the shock has become detached from the merging cluster that drove it, known as the ‘runaway’ phase (see definitions in Zhang et al. 2019). It has reached a distance of a little under 1.8 × R₅₀₀, which is within the range expected for radio relics (see, e.g. Bagchi et al. 2011; Erler et al. 2015).

Fig. 2. Cosmological simulation of a galaxy cluster undergoing a major merger at z = 0.14. Four left-most panels, clockwise from top left: Projected shock-dissipated energy rate, gas pressure, gas density, and dissipation-weighted Mach number, respectively. Projections have a depth of ±7.5 Mpc from the cluster centre. Right-most panels: enlarged cut-outs showing slices through the midplane of the projection. The white, dashed box indicates the size of the upstream region in our idealised shock tube simulations. The merger drives a shock wave out to the cluster outskirts. The collision of this shock wave with an accretion shock leads to the production of a thin shell of compressed gas (see white arrow). A movie of this process can be found here.

In the top left panel of Fig. 2, we show the projected shock-dissipated energy, where this has a depth of 15 Mpc. The merger has led to the formation of two megaparsec-sized, arc-like shocks. These align with the merger axis and are relatively symmetric, as expected when the initial cluster masses are similar (van Weeren et al. 2011a; Lee et al. 2025). The shocks from the major merger are superimposed on a background of less energetic shocks. Some of these, such as those towards the lower-right are from ongoing minor mergers, whilst the majority result from accretion.

The bottom-left panel has a two-dimensional colourbar, where the colour indicates the dissipated energy-weighted Mach number, with saturation being set proportional to this energy. It can be seen that the merger shock is made up predominantly of low Mach numbers (ℳ ≲ 5). Shocks at the edge of the cluster, meanwhile, are much higher (ℳ ≳ 100)¹⁵. These shocks define the filamentary structure, within which the cluster sits. Indeed, the merging cluster has arrived from the filament evident at the top right of the panel (see movie).

In the middle column, we show the projected pressure (top) and density (bottom) respectively. The brightest central galaxy (BCG) is located at the peak value in both panels. Infalling galaxies can also be seen as bright spots at the edge of the cluster. The merger has led to a region of increased pressure centred on the cluster. This region is well-delineated, in contrast to the more diffuse density distribution. It is also clear that the region of high pressure is bounded by the merger shocks.

To remove the effects of projection, in the right-hand column we show slices of the same quantities. We zoom in on a region of size 2 × 2 Mpc² in order to better explore the shock details. At this time and position, the outwards-moving merger shock encounters an inwards-moving accretion shock¹⁶¹⁷. These are evident as ridges in the density and pressure slices in the corresponding movie. The time at which the two shocks meet can be considered a Riemann problem with density, pressure and velocity specified at both sides of the discontinuity. The result of this is two-fold: firstly, an inwards travelling rarefaction wave is generated, and secondly, a thin, shock-compressed density shell forms behind the merger shock. This feature is marked by the white arrow in the bottom right panel, and is consistent with previous studies (Zhang et al. 2019, 2020).

It can be seen that the pressure downstream of this feature is approximately constant, whilst the density jumps. This indicates that the edge closest to the arrow is actually a contact discontinuity. The jump itself is exacerbated by the rarefaction wave, which lowers the density downstream of this discontinuity. Such waves have also been observed in previous studies of cluster mergers (Shi et al. 2020). This feature is key to several results shown in the remainder of the paper.

The subsequent evolution of the shock wave can be well-approximated as a pressure-driven shock wave on a flat density background (see Sect. 2.3). This lends itself to shock tube simulations, where we can approach the problem with significantly higher resolution. This is necessary, as whilst the overall pressure and density profiles of our simulated clusters are converged (Berlok et al., in prep; Tevlin et al. 2025), density turbulence is not sufficiently resolved at these radii. In both of the right-most panels, we have added a white, dashed box with dimensions 450 kpc × 300 kpc, which represents the size of the high-resolution region in our shock tube initial conditions. This region is traversed by only 𝒪(10) cells, which is clearly insufficient. The shock tube method also allows us to precisely define the various turbulent parameters, allowing for a clearer investigation into their impact on downstream properties.

4. Shock tube analysis

Our shock tube setup, including initial pressure and density values, is informed by the above cosmological analysis, with values given in Sect. 2.3. We note that the pressure jump dominates the white, dashed box shown in Fig. 2, and so chose to make the shock exclusively pressure-driven, i.e. the downstream density is set equal to the mean upstream density. The mass resolution in the high-resolution upstream region of these simulations (i.e. region III) is approximately 730 times greater than in the cosmological simulation, allowing for the resolution of density turbulence.

As previously mentioned, multiple studies have shown that the ICM is clumpy. We implemented this in our shock tube simulations using the method described in Sect. 2.3.4. The injection scale of the density turbulence was set to half the smallest box dimension, i.e. 150 kpc, and the amplitude of the fluctuations was taken from Zhuravleva et al. (2013). We remind the reader that there is no initial velocity turbulence; that is, the density fluctuations are initially at rest.

We show the state of our fiducial Mach 3 simulation at t = 180 Myr in Fig. 3, where the shock front is travelling from left to right. Each panel shows a 2D slice (except for panel iii, which is a thin projection) and x = 0 marks the beginning of ‘region III’ (see Sect. 2.3). The time was chosen such that several characteristics of the shock are on display simultaneously. In particular, we draw attention to the following three consequences of including upstream density turbulence: (i) the formation of downstream turbulence, (ii) shock corrugation, (iii) the formation of a Mach number distribution.

Fig. 3. Mach 3 shock is driven into a turbulent upstream density field. The initial pressure and density values for the shock tube correspond to those shown in Fig. 2. Panels: (i) gas pressure, (ii) gas density, (iii) dissipated-energy-weighted Mach number, (iv) the fraction of the gas speed in the x-component, (v) the x-component of gas velocity, (vi) gas speed divided by sound speed. Data are shown in slices except for (iii), which is a thin projection of 35 kpc. All values are shown in the shock rest frame at t = 180 Myr. Upstream density turbulence directly leads to shock corrugation, a distribution of Mach numbers at the shock front, and downstream velocity turbulence. An animated version of this figure can be found here.

4.1. Downstream turbulence and shock corrugation

Perhaps the most striking of these three is the generation of downstream turbulence¹⁸. The impact of this can be seen in almost all panels, but we start with panel (iv). Here we show the fraction of the shock-frame gas velocity in the x-component. To calculate the shock-frame, we took the median x-coordinate for all shock surface cells at each snapshot. We then converted this to an average shock velocity by dividing by the snapshot cadence. Finally, we subtracted the subsequent value from the lab-frame gas velocities. In a purely laminar flow aligned with the shock normal, there would be no deviations from |υ_x|/|υ| = 1, as is the case in the simulations without density turbulence.

This panel is supported by panel (v), which shows the x-component of the velocity. Blue shades indicate gas moving from right to left, whilst red shades indicate gas moving the opposite way. These values are, once again, both given in the frame of the shock. The reduction in speed after crossing the shock is a standard theoretical result, however, panel (v) shows that including density turbulence causes gas to move away from the shock front at varying speeds. Indeed, in some cases, gas even moves towards the shock front from the left. Whilst we directly model CR electrons later in the paper, this already shows us that distance from the shock front is not a reliable indicator of time since injection.

The turbulence results predominantly from a Rayleigh-Taylor instability, which leads to the formation of density ‘fingers’. These are visible, for example, in panel (ii) at y ≈ 120 kpc and y ≈ 290 kpc. This spacing is no coincidence, as such fingers form at positions where the shock reaches higher density clumps, which in turn is set by the injection scale. In the non-linear regime of this instability, the resulting velocity shear causes the Kelvin-Helmholtz instability to form as a parasitic instability. The effect of this can be seen particularly clearly in the upper half of panel (ii) and the associated movie.

We explain the emergence of the Rayleigh-Taylor instability as follows: density fluctuations result in the shock speed varying along its surface, with the shock accelerating into lower density regions. An initially planar shock will subsequently warp. This leads to the pressure dropping behind the most advanced part of the shock, as can be seen in panel (i). The consequence of this is that the pressure gradient starts to point towards the far downstream. The density jump from the downstream behind the contact discontinuity towards the shock-compressed region, however, provides a gradient pointing the opposite way. What is more, as higher density fluctuations are less easily accelerated by the shock, the back of the compressed region also begins to warp. Indeed, as seen in Fig. 3, this happens to a greater extent than at the shock front.

We show a schematic that illustrates this scenario in Fig. 4, with orange colours representing pressure and blue colours representing density, respectively. Pressure and density gradients have become misaligned, which, due to the inviscid vorticity equation, results in a baroclinic torque, as

$$\begin{array}{c}\hfill \frac{\mathrm{D}\mathit{\omega }}{\mathrm{D}t}=\frac{1}{{\rho }^2}\mathbf{\nabla }\rho \mathbf{\times }\mathbf{\nabla }P,\end{array}$$(13)

where the left-hand term is the Lagrangian derivative of vorticity ω. This mechanism rapidly generates downstream velocity turbulence, and leads to a reinforcement of the finger-like structures, with eddies directed as shown by the white, dashed arrows.

Fig. 4. Schematic showing how upstream density turbulence causes velocity turbulence to be generated downstream. The corrugation of the shock front naturally leads to a misalignment of pressure and density gradients (shown here by red and green arrows, respectively). This results in a baroclinic term, which in turn induces vorticity, causing the contact discontinuity to become Rayleigh-Taylor unstable.

Upstream density fluctuations are critical for the formation of this scenario. Without them, there is only a very limited pressure gradient and no corrugation at the back of the shock-compressed region. The later condition in particular would mean that ∇ρ and ∇P were approximately anti-parallel, leading to a vanishing cross-product in Eq. (13) and hence little to no growth of the instability. We show this to be true in our simulations in Appendix B.

We note that Hu et al. (2022) identified the related Richtmyer-Meshkov instability as causing turbulence in their shock tube simulations. We may discount this here as the Richtmyer-Meshkov instability takes place at the shock front, whilst the instability we observe takes place at the contact discontinuity. Vorticity may also be generated by curved shock surfaces according to the theorem by Crocco (1937). We discount this as the primary driver in this case, however, as this mechanism would produce vorticity oriented in the opposite direction. This would reduce the size of the fingers rather than increase them.

4.2. Mach number distribution

4.2.1. Form

We now turn our attention to the third consequence of including density fluctuations: the formation of a Mach number distribution. It has long been known that shocks in clusters vary in strength (Ryu et al. 2003; Kang & Jones 2005; Pfrommer et al. 2006; Hoeft et al. 2008; Skillman et al. 2008; Vazza et al. 2009). It is also inferred from synchrotron fluctuations that this is true within radio relics as well (Rajpurohit et al. 2020a, 2021). This was supported by the shock tube simulations of Domínguez-Fernández et al. (2021a). Like them, we find that an initially well-defined Mach number decomposes into a distribution upon encountering turbulence. This can be seen in the top right panel of Fig. 3, which shows a thin projection¹⁹ of depth 35 kpc, where colours indicate the Mach number at the shock front.

The shock in Fig. 3 was initially set up with ℳ = 3. However, whilst many cells continue to show similar values, it is clear that that there is now a distribution, with values lying approximately between ℳ = 1.9 and ℳ = 4.8. We quantify the development of this distribution in Fig. 5. For this, we calculated the distribution at each snapshot, where each cell has been weighted by its contribution to the overall shock surface. We show the median of these distributions in Fig. 5 as solid lines and provide the 25−75% percentile range as shaded values.

Fig. 5. Mach number distributions for all models, where each cell has been weighted by its normalised contribution to the shock surface. Lines indicate the median taken over all snapshots, whilst the shaded values indicate the interquartile range. The addition of magnetic turbulence (Flat) to a homogeneous density distribution (Flat-ConstB) broadens the distribution only very mildly. By contrast, the addition of upstream density fluctuations in the simulation (Turb and Turb-ConstB) turns an extremely narrow distribution into a much broader one. The width of the distribution is proportional to the peak Mach number.

The dotted grey line indicates the values for simulation Flat-ConstB, which has neither density nor magnetic turbulence in the initial conditions. It can be seen that the distribution forms an approximate delta function. This indicates the base-line accuracy of our shock finder. When we add magnetic fluctuations, as in Flat, the distribution broadens slightly. This is due to the additional pressure fluctuations caused by the magnetic field (see Appendix E). However, this effect is totally dominated by the impact of including density fluctuations. Indeed, there is essentially no difference between the Mach number distributions in Turb-ConstB and Turb. The initial conditions for these simulations differ only by their inclusion of a uniform or turbulent magnetic field, respectively.

The peak of the distribution meanwhile remains at the initial Mach number (ℳ = 2 in the top panel and ℳ = 3 in the bottom). Finally, we note that the distribution remains remarkably stable over time; the difference between the 25th and 75th percentile is very low. This is in contrast to Domínguez-Fernández et al. (2021a) who, whilst also finding a tail towards higher Mach numbers, do not see such stability (see their Fig. B2). We attribute this to our different methods of generating turbulence (see Sect. 5).

4.2.2. Origin

Returning to Fig. 3, it can be seen from panel (iii) that the higher Mach numbers are typically found towards the back of the shock front, whilst lower Mach numbers are found at its most advanced position. An intuitive understanding of this can be gained by inspecting panel (vi). Here we show a slice, where colours indicate the shock-frame gas speed as a fraction of the sound speed. We expect this quantity to provide a fair approximation to the more rigorous result given by Eq. (2). We remind the reader that we have defined the shock-frame based on the median of the local shock velocities. This is reasonable as the local velocity fluctuations along the shock front are small compared to its median shock speed.

Gas is stationary in the lab-frame and so, under our approximation, approaches the shock at the same speed, as can be seen in the right-hand side of panel (v). Additionally, the gas is in approximate pressure equilibrium, as can be seen by inspecting panel (i). Consequently, |υ|/c_s in the upstream is determined purely by the density fluctuations, as $c_{\mathrm{s}}=\sqrt{{\gamma }_{\mathrm{a}}P/\rho }$, where the adiabatic index is set to γ_a = 5/3. This means that denser (more rarefied) gas has a lower (higher) sound speed, respectively. The consequence of this is that we expect denser gas to produce higher Mach numbers, and more rarefied gas to produce lower Mach numbers. This is indeed the case, as can be seen by comparing panels (ii) and (iii). As discussed in Sect. 4.1, lower density gas also causes the shock front to accelerate locally. As a result, lower Mach numbers are found further forward. We show that this is generically true across the shock front in Appendix C.

4.3. Radio versus X-ray Mach number discrepancy

It has been argued previously that the highest Mach numbers in a distribution are responsible for the radio emission (see e.g. Hoeft et al. 2011; Wittor et al. 2021a; Domínguez-Fernández et al. 2021a). This is currently the leading theory for relieving the tension in the so-called ℳ_radio − ℳ_{X − ray} discrepancy, where Mach numbers derived from radio observations are found to be higher than those derived from X-ray observations. The effect has yet to be modelled end-to-end, however, using a full Fokker-Planck solver. We resolve this issue here, using the tracers discussed in Sect. 2.3.5 and the CR electron post-processing code CREST, introduced in Sect. 2.5.

We show the resulting volume-weighted non-thermal electron spectra at t = 250 Myr in the upper panels of Fig. 6. At this time, the shock in the Mach 3 simulations has reached the end of the high-resolution upstream region (region III). We remind the reader that we present dimensionless momenta $p=\stackrel{\sim }{p}/(m_{\mathrm{e}}c)$ and that we use the 1D distribution function, f^1D = 4πp²f^3D. We have further multiplied the spectra by p^2.2, which is the maximum slope of injection we allow (see Sect. 2.5). Using this convention, such a slope would appear as a horizontal line.

The grey, dashed lines show spectra from our Flat simulations, which do not include density fluctuations. We have calculated their spectral slope, a_e, flat, in each panel by fitting a straight line between 10 < p < 10³ using the method of least squares. This is approximately the region unaffected by cooling.

If we assume a single Mach number is responsible for the CR electron spectra slope, α_e, we may derive the following formula from the jump conditions (Enßlin et al. 1998; Clarke et al. 2011):

$$\begin{array}{c}\hfill \mathcal{M}=\sqrt{\frac{2(2-{\alpha }_{\mathrm{e}})}{1-2{\alpha }_{\mathrm{e}}-3{\gamma }_{\mathrm{a}}}},\end{array}$$(14)

or equivalently

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{e}}=-\frac{2({\mathcal{M}}^2+1)}{{\mathcal{M}}^2-1},\end{array}$$(15)

where we set the adiabatic index γ_a = 5/3. These formulae predict that ℳ = 2 and ℳ = 3 shocks should produce slopes of α_e = −3.33 and α_e = −2.5, respectively. The actual values of α_e, flat = −3.36 and α_e, flat = −2.49, on the other hand, imply ℳ = 1.98 and ℳ = 3.03. Whilst some of this variation is due to the distribution of Mach numbers shown in Fig. 5, at higher Mach numbers the error arises predominantly due to the fact that Eq. (14) has a pole at α_e = −2 and consequently shallow slopes must be measured with exponentially more accuracy. With this in mind, we consider the Flat runs to be good approximations of single Mach number shocks. That is, their slopes are what would be expected if the radio-derived Mach number matched that derived from X-ray observations, which tend to average over upstream and downstream properties (Wittor et al. 2021a).

The solid black lines in Fig. 6 show spectra from our Turb runs, which include density fluctuations. It can be seen that the resultant spectral slopes, labelled α_turb, are systematically shallower than in the Flat case. Indeed, Mach numbers derived from these slopes would imply shocks with strength ℳ = 2.54 and ℳ = 3.25 for our two simulations with initial Mach numbers of ℳ = 2 and 3, respectively²⁰. This is despite the fact that the shock strength continues to peak at ℳ = 2 and ℳ = 3, respectively, as seen in Fig. 5.

Fig. 6. Top row: Volume-weighted non-thermal electron spectra generated from our Turb (solid black) and Flat (dashed grey) simulations at t = 250 Myr. Blue lines show contributions to the Turb spectrum, where tracers have been binned by time since injection. Bottom row: Histograms indicating the total number of injected tracers relative to the shock front, where bins have a width of 2 kpc. Blue and grey colours represent our Turb and Flat simulations respectively, whilst colour saturations indicate the same time bins as above. The broadened Mach number distribution seen in Fig. 5 causes a shallower slope, αe, in the Turb runs compared to the theoretical expectation. This is especially noticeable for low Mach number shocks. In the case of upstream density fluctuations, substantial mixing takes place downstream, such that distance from the shock front is no longer a good indication of cooling time.

It is notable that the discrepancy is greater in weaker shocks. This might seem to be in contrast to our finding in Fig. 5, which shows that the tail of the Mach number distribution actually extends further in stronger shocks. This effect is outweighed, however, by the strong inverse dependence of α_e on ℳ, with α_e having a strong dependence on ℳ at ℳ = 2, but a substantially weaker dependence at ℳ = 3. For example, from Eq. (15), it can be seen that a shock with ℳ = 2.5 would produce a spectral slope of α_e = −2.76, which is 0.57 shallower than that of a ℳ = 2 shock. On the other hand, a ℳ = 3.5 shock will produce a spectral slope of α_e = −2.36, which is only 0.14 shallower than an ℳ = 3 shock. The tail of the Mach number distribution is therefore more influential in weaker shocks, as the shallower slopes from these fluctuations can more easily affect the integrated spectra at higher frequencies, thus increasing the discrepancy between the expected and inferred Mach number.

Of course, in radio observations, only the momenta range 10⁴ ≲ p ≲ 10⁵ is available for observation, as it is this range that produces megahertz- and gigahertz-emitting electrons. Standard theory dictates that the slope of the CR electron spectrum in this range should be a_e − 1 due to cooling²¹ (Condon & Ransom 2016). Indeed, it is this fact that is used to infer the spectral slope in radio relic observations (see, e.g. van Weeren et al. 2012a; Rajpurohit et al. 2020a). This result relies, however, upon the spectrum being produced by a single Mach number shock. We have added the theoretical slope to Fig. 6 for the Turb runs to guide the eye. It can be seen that, when the shock is made up of a distribution of Mach numbers, the slope in this region is actually slightly shallower than a_e − 1. Once again, this is especially evident in the ℳ = 2 run. Indeed, if we measure the slope in this range for the Mach 2 Turb run, we would infer a spectral slope of α_e, turb = −2.5. This in turn produces an estimated Mach number of ℳ = 3.01, which is ∼50% higher compared to the true peak value. For comparison, we show the slope measurements and inferred Mach numbers for all of our runs in Table 2.

As previously, this effect originates due to the shallower slopes of higher Mach number shocks. At the higher momenta end, where f(p) is lower, individual tracers can have a larger impact on the volume-weighted mean. The effect can be more clearly seen if we plot the spectra binned by the time elapsed since their injection at the shock front. These bins are given in the legend in the top right panel of Fig. 6. We have chosen the time bins such that the number of tracers in them approximately doubles each time, which produces an approximate doubling of the amplitude of the binned spectra in regions where the cooling time exceeds the simulation run-time. The maximum time since injection is 244 Myr, as we ignore tracers injected during the initial buffer region (see Sect. 2.3.3). As the electrons cool rapidly at p ≳ 10⁴, binning in this fashion allows us to pick out regions of momenta which are dominated by a specific age. Indeed, it is this layering that produces the overall volume-weighted slope. It can be seen that the binned spectra bend slightly, becoming slightly shallower at higher momenta. This results from the layering of different injected slopes as tracers sample the Mach number distribution. It is this bending effect that ultimately produces a slope of less than a_e − 1.

4.4. Breaking of the laminar flow assumption

It is frequently assumed that distance from the shock front is a reliable measure for time cooled (see e.g. van Weeren et al. 2012a; Rajpurohit et al. 2020a). We now investigate whether this assumption is valid, particularly in light of the turbulence we identified in Sect. 4.1. For this purpose, we reuse the time bins just discussed, making histograms of the distance from the shock median for each binned population. These histograms are presented in the bottom row of Fig. 6. The grey shaded regions indicate data from the Flat simulations. It can be seen that whilst there is some overlap between time bins for the Flat simulation, owing to the aforementioned magnetic pressure fluctuations, the amount is relatively minimal. This indicates that the flow is predominantly laminar. In contrast, tracers from the Turb simulations, represented by the blue shaded regions, overlap strongly. Indeed, there are tracers in the last time bin that are still close to the shock front.

Two further aspects further distinguish the Turb simulations from their Flat counterparts: firstly, there are well-populated regions of negative distance due to the corrugation of the shock front (see Sect. 4.1), and secondly, the maximum distance from the shock is substantially larger in the Turb runs. Indeed, the relative maximum distance increases proportionally with the strength of the shock: tracers reach 32% and 67% further in Mach 2 and Mach 3 shocks, respectively. This comes about primarily due to the inertia of the high-density fluctuations; higher density clumps experience less acceleration and can therefore better penetrate downstream. Higher Mach number shocks, meanwhile, have faster speeds, resulting in a greater distance opening up between the shock front and such clumps. This effect is further reinforced by the Rayleigh-Taylor instability as analysed in Sect. 4.1, which itself is more effective at faster shock speeds, when the downstream is more dynamic.

The consequence of all of these effects is that, once turbulent density fluctuations are taken into account, distance from the shock front is no longer a good indicator for time since injection for individual electrons. Indeed, as we show in Appendix D, it is not possible to make an analogous figure to Fig. 6 where the tracers have been binned by distance.

These effects can be seen particularly well in Fig. 7, where we show slices through the injected region in our Mach 3 Turb simulation at t = 180 Myr. The slices presented here may be directly compared with those presented in Fig. 3; data have simply been taken from the tracers here rather than the gas cells. In panel (i) of Fig. 7, we show the time since injection, with lighter colours indicating more recent injection. Regions without colour indicate a lack of injection. Within the shock-compressed region, this results from the shock front falling below the critical Mach number (see discussion in Whittingham et al., in prep.). The impact of the range of post-shock gas speeds, as already shown in panel (v) of Fig. 3, can clearly be seen, with relatively freshly injected CR electrons evident at varying distances from the shock front. With this said, the overall distribution of CR electron ages mostly traces the shape of the Rayleigh-Taylor fingers, with extra perturbations formed by the additionally generated turbulence. As this instability acts at the contact discontinuity, the resultant eddies are generally populated by older material. This effect is also responsible, however, for bringing the same aged CR electrons closer to the shock front.

Fig. 7. Slices through our fiducial Mach 3 simulation at t = 180 Myr showing: (i) time since injection, (ii) CR electron number density at p = 2 × 104 (see text), (iii) magnetic field strength, and (iv) synchrotron emissivity at 150 MHz. Only tracers that have undergone DSA are shown. A Rayleigh-Taylor instability induces a counter-streaming plume (in the shock rest frame) that brings aged, and therefore cooled, electrons close to the shock front. However, these electrons are still relatively bright at 150 MHz due to magnetic field amplification up to μG levels. An animated version of this figure can be found here.

In panel (ii) of Fig. 7, we show the CR electron distribution function f(p), where we have set p = 2 × 10⁴. We choose this momentum in particular as it is the one which contributes most to the 150 MHz emission²². Older electrons have had more time to cool, and so f(p) at a fixed momentum generally reflects the distribution shown in panel (i). Some noticeable differences are evident at the shock front, however. Here it can be seen that the most advanced part of the shock has lower f(p) values compared to neighbouring regions. This is due to the weaker shocks that take place here, as may be seen by comparing with panel (iii) in Fig. 3. Such shocks produce spectra with steeper slopes and lower normalisations²³, which act to reduce f(p) values. Both of these effects are also evident in Fig. 6.

Observationally, f(p) is not available to us; instead, we must use synchrotron emission to infer it. However, this is modulated by the magnetic field strength of the emitting region, as can be seen by inspection of Eq. (9). We therefore focus now on the magnetic field downstream.

4.5. Magnetic field amplification

4.5.1. Strength

As discussed in Sect. 2.3.4, the upstream turbulent magnetic field in our simulations has an RMS field strength of B₁ ≈ 0.16 μG. In the ℳ = 3 case without density fluctuations we expect a shock compression ratio of x_s = 3. Assuming that the magnetic components perpendicular to the direction of shock propagation are amplified by this same ratio, for an isotropic magnetic field we should expect the RMS strength immediately behind the shock to be

$$\begin{array}{c}\hfill B_2=B_1\sqrt{\frac{1+2x_{\mathrm{s}}^2}{3}},\end{array}$$(16)

and hence B₂ ≈ 2.5B₁ = 0.4 μG. This is problematic, as multiple studies now indicate that magnetic fields in radio relics are at least a factor 3 higher than this (see citations given in Sect. 1).

In panel (iii) of Fig. 7, we show the magnetic field strength in the injected region. It is evident that values indeed start around B = 0.4 μG. However, many regions of the downstream actually exceed μG strengths, in accordance with observed radio relics. Even at its highest strengths, however, the magnetic field remains dynamically subdominant. As can be seen in Fig. 8, minimum values reach P_th/P_B = 10, with the majority having P_th/P_B > 100²⁴. Moreover, peak magnetic field strengths also remain below the equivalent value for the CMB magnetic field strength at redshift z = 0.2, B_CMB ≈ 4.7 μG. This means that inverse Compton cooling still dominates over synchrotron cooling in the downstream.

Fig. 8. Slices through our fiducial Mach 3 simulation at t = 180 Myr showing plasma beta values. The grey lines mark the corrugated shock surface and the contact discontinuity (i.e. the region within which tracers have been injected). Despite significant amplification, beta values are typically well above 10, limiting the ability of the magnetic field to affect dynamics.

We quantify the range of magnetic field strengths reached in Fig. 9, where we show the field strength as a function of density for the initial upstream and final downstream regions. To do this we take all gas cells in the high-resolution upstream region (region III) at t = 0 and all gas cells at t = 250 Myr that are in the shock-compressed region. Contours of these distributions are shown in blue and red, respectively.

Fig. 9. Left: Phase space diagram of magnetic field strength versus gas number density for our Mach 3 Flat-TurbB simulation. Contours cover 25%, 50%, 75%, and 95% of the probability density (represented by darker to lighter shades, respectively). Blue colours indicate the initial upstream distribution, whilst red colours show the final state in the injected region at t = 250 Myr. A cross marks the median of the distribution. Right: As previous, except data are taken from the Mach 3 fiducial simulation. Dashed lines indicate the empirical scaling of the 75% boundaries. Amplification is weaker than the naive expectation for Flat-TurbB due to magnetic decay. Meanwhile, for the fiducial simulation, amplification is significantly stronger than that possible from compression alone.

It can be seen that the fiducial simulation, Turb, indeed regularly reaches μG strengths and generally significantly exceeds those achieved in the Flat-TurbB simulation. Moreover, the Flat-TurbB simulation fails to match the naively expected amplification: with x_s = 3 and an amplification of the magnetic field by a factor of 2.5 according to Eq. (16), we should expect a scaling of B ∝ n^α_B, with α_B = log(2.5)/log(3)≈0.83. Empirically, however, we recover B ∝ n^0.5. This is a direct result of turbulent magnetic decay caused by magnetic tension forces. Indeed, we show in Appendix F that the B ∝ n^0.5 scaling can even be predicted, once this process is taken into account.

The decay of the magnetic field significantly exacerbates the problem of μG magnetic fields in radio relics; without subsequent amplification in the downstream, 0.4 μG becomes a maximum value for our simulation²⁵. This makes the empirical scalings of the fiducial simulation all the more remarkable; compression alone cannot produce scalings greater than B ∝ n, and magnetic decay acts to reduce the field strengths over time as well.

4.5.2. Origin

To understand the scalings we observe, we move to analysis of the individual magnetic components, which we show in Fig. 10. The x-component, which is parallel to the direction of shock propagation, should initially experience no changes across the shock, due to the ∇⋅B = 0 constraint. On the other hand, following flux conservation, the y- and z- components should be amplified by the shock compression ratio, x_s. This would lead to scalings of B = const. and B ∝ n, respectively. Indeed, by looking at individual tracer trajectories, we find that the magnetic field is indeed initially amplified in this way²⁶. It can be seen on the left-hand side of Fig. 10, however, that for the Flat-TurbB simulation we ultimately recover scalings of B_x ∝ n^−0.1 and B_y, z ∝ n^0.6, respectively. This is a result of the aforementioned magnetic decay: the red contours, taken at t = 250 Myr, average over regions of freshly amplified and decayed magnetic fields.

Fig. 10. Left: As Fig. 9 but now showing individual components of the magnetic field. All components are more highly amplified in the fiducial simulation owing to time-varying compression. It is particularly notable, however, that the x component is the most highly amplified. This is indicative of amplification by shearing.

Inspecting the right-hand side of Fig. 10, we find that there are several differences compared to the Flat-TurbB simulations. Firstly, the top and bottom of the contours no longer scale the same way. In fact, the y- and z-components have a more extended tail towards low B-values, compared to the Flat-TurbB. This is likely both a result of the extended Mach number distribution, which produces a distribution of x_s values, but also due to regions of rarefaction taking place behind the shock front, as are visible in Fig. 3 and the associated animated figure.

The upper contours also show much stronger scaling. Indeed, the y- and z-components are very close to B ∝ n. This is possible if compression happens at late times, when the magnetic field has less time to decay. Regions of compression can be seen in Figs. 7 and 8; they are the filament-like shapes, which are predominantly sheets seen in cross-section. Indeed, this compression now happens across all three components, as the Rayleigh-Taylor instability also leads to the gas being compressed in the y- and z-directions.

However, whilst a significant amount of amplification arises from compression, some decay will still naturally occur, and compression cannot account for the peak scaling observed for the x-component. Indeed, it is notable that this component experiences the strongest amplification. By applying similar analysis to our lower resolution simulations (see Appendix G), we find that any difference in magnetic field strength is explained by a corresponding change in density. As a magnetic dynamo should also produce higher strengths for the same density, this implies that a dynamo is either not at all, or only weakly, active in our simulations. This is supported by the animated version of Fig. 3, where it can be seen that eddies experience few if any turnovers, and a full turbulent cascade does not have time to form.

The varying speeds of the downstream flow (see panel v of Fig. 3) does, however, lead naturally to the stretching and shearing of the magnetic field. This process is easily able to produce B ∝ n^α_B scalings with α_B > 1. Moreover, magnetic filaments pointing in the x-direction are clearly seen in panel (iii) of Fig. 7 and in Fig. 8, as would be formed by shearing motions. We hence conclude that shearing is likely critical for the production of μG magnetic fields in radio relics²⁷.

4.5.3. Impact on observational inferences

Returning to Fig. 7, we may now inspect how the magnetic field strength affects synchrotron emissivity, which is shown in panel (iv). This is calculated using CRAYON+, as explained in Sect. 2.6, and is shown in this figure at 150 MHz. The colourbar for this panel is given the same overall range as was given to f(p = 2 × 10⁴). It can be seen that the dynamic range, however, is substantially smaller in panel (iv). This is due to the significant amplification that takes place towards the back of the shock-compressed region. There are two ways in which this impacts the synchrotron emissivity. Firstly, as can be seen by inspecting Eq. (9), there is a linear scaling of the emission with the component perpendicular to the line of sight. Secondly, through Eq. (10), a change in the magnetic field strength shifts the momenta at which most emission takes place. Specifically, as ν_c ∝ Bγ², an increased magnetic field strength means that lower momenta, which have higher values of f(p), contribute more at the same frequency of emission. Together these effects boost emission from the rear of the shock-compressed zone, thereby evening out some of the effects of cooling. Indeed, in some regions where f(p) is low, such as towards the back of the shock at y ≈ 210 kpc, the synchrotron emissivity is now higher than at the front of the shock. In general, however, the emissivity still decreases from front to back. We also note that the brightest regions are still at the front of the shock, at the parts where we previously identified the highest Mach number shocks. The front of the shock undergoes milder magnetic field amplification and hence milder modulation, and so our analysis in Sect. 4.4 regarding features close to the shock front holds for the synchrotron emission too.

We now consider how this modulation affects how we infer variables from observations of radio relics. As stated earlier, synchrotron emission is our primary window into CR electron properties in radio relics, and so it is critical to understand how observations are biased. In Fig. 11, we investigate the impact of synchrotron-weighting on the projected time since injection and projected magnetic field strength. On the left-hand side, we show the volume-weighted values, whilst on the right-hand side we show the synchrotron-weighted projections. We use the 150 MHz channel for this, as previously. We only include values where the line of sight crosses at least one injected particle. The colourbars in each row have the same scale, so panels may be directly compared.

Fig. 11. Projections with a depth of 300 kpc through our fiducial Mach 3 simulation at t = 180 Myr. These show: i) volume-weighted time since injection, ii) synchrotron-weighted time since injection, iii) volume-weighted magnetic field strength, and iv) synchrotron-weighted magnetic field strength. Regions with no injected tracers have been masked (see text for further details). Synchrotron emission is calculated at ν = 150 MHz and is biased towards fresher injection and stronger magnetic fields. Observations made through this channel consequently appear to show a relatively laminar flow and μG-strength magnetic fields, even though this is not typical of the downstream.

In panel (i), we project the time since injection across all tracers. We set tracers that were not injected to have values equal to the simulation run-time (t = 180 Myr). This leads to a gradual fade-out towards the left-hand side, as tracers age and as we project through non-injected tracers as well. This effect also highlights the Rayleigh-Taylor fingers, where the majority of injected tracers are. These are not as clear as in Figs. 3 and 7 as we are now projecting through significant amounts of turbulence. On the right-hand side of panel (i), the non-injected tracers reduce the projected average in the most advanced parts of the shock, leading to these regions being darker than expected. This is a useful feature, however, as it shows us that the back of the shock is predominantly planar; it is the front of the shock that corrugates, as it advances into regions of lower density. Indeed, we could see this effect in cross-section in panel (iii) of Fig. 3. We should consequently expect the brightest part of the simulated radio relic to also be relatively flat, as it is here that the shock-dissipated energy is highest (see Sect. 4.2.2 and Appendix C). We revisit this analysis in Sect. 4.6.1.

When we weight the projection by the 150 MHz channel, as in panel (ii), we start to recover the relationship between distance from the shock and time since injection; on large enough scales, distance once again becomes a reasonable proxy for age (see Appendix D). On kiloparsec-scales, however, there is still substantial variance. This typically takes the form of horizontal striations. An exception to this is towards the very front of the shock, where two darker, curved regions are visible behind the shock corrugations. If we compare this panel with Fig. 3, which shows the same simulation at the same time, we can see that this is due to the turbulent transport of older material towards the shock front. As we showed in Sect. 4.1, turbulence is most effective behind the regions where the shock corrugates.

In panel (iii), we show the volume-weighted projected magnetic field strength. This reaches typical values of around 0.4 μG, with peak values reaching closer to 0.65 μG, as expected from Fig. 9. When we weight the projection by synchrotron emission at 150 MHz, as shown in panel (iv), we find a significant shift towards higher field strengths. Now, a substantial amount of the projection is above 1 μG, more closely in line with observational inferences from radio relics. Once again, the emission-weighted projection predominantly shows horizontal striations, especially towards the rear of the shock-compressed zone. Indeed, this explains the patterns seen in panel (ii): the synchrotron emission is proportional to B_⊥^{1 − α_ν} (see Eq. (9)), hence CR electrons moving along these features dominate the projected quantities.

Magnetic field features stretched parallel to the x-axis could also be seen in Figs. 7 and 8. By comparing these plots with the velocity structures evident in Fig. 3, we reinforce our previous conclusion that these features originate from the shearing of gas. Curved filaments of high magnetic strength are also evident, however. Through similar comparisons, it can be seen that these form in the troughs between Rayleigh-Taylor fingers, which may point to adiabatic compression (see Fig. 7 again). All of these mechanisms require turbulence, and this is mainly generated at the contact discontinuity. Consequently, the strongest magnetic field values in Fig. 11 are seen towards the back of the shock-compressed region, with minimal extra amplification at the very front of it.

In Fig. 12, we quantify just how much the synchrotron-weighted projections are biased relative to their volume-weighted counterparts. We show the initial state in the high-resolution region (region III) in dotted lines, and the final state in the shock-compressed region at t = 250 Myr in solid lines, with blue indicating volume-weighting and orange indicating synchrotron-weighting. The distribution is initially narrow and peaks at approximately 0.16 μG, matching the chosen RMS value (see Sect. 2.3.4). By construction, the initial distribution is the same in both Mach 2 and Mach 3 simulations. At the end of each simulation run, the volume-weighted distribution has broadened and now shows a double peak. This is a line-of-sight effect; with the left-hand peak representing lines of sight that predominantly pass through regions outside the shock-compressed region, and the right-hand peak representing the opposite scenario²⁸. This effect biases these projections slightly low relative to the true volume-averages in the injected region (see values in Fig. 9). The projected volume-weighted RMS magnetic field strengths in the Mach 2 and Mach 3 simulations are 0.23 μG and 0.33 μG, respectively.

Fig. 12. Left: Probability density distributions of the projected magnetic field strength in our fiducial Mach 2 simulation. Solid (dashed) lines represent projections through the shock-compressed region at t = 250 Myr (high-resolution region at t = 0 Myr). Blues lines represent volume-weighting, whilst orange lines represent synchrotron-weighting, where the synchrotron emission is calculated at ν = 150 MHz. Arrows are placed at the RMS magnetic field strength calculated for each distribution. Right: As previous, except data come from the fiducial Mach 3 simulation. Synchrotron-weighting significantly overestimates the average magnetic field strength. Magnetic field values are able to reach μG values even in weaker shocks, although amplification is more effective at higher Mach number shocks.

The synchrotron-weighted distribution, however, shows significantly higher RMS values of 0.48 μG and 0.78 μG, respectively, which is more than double the volume-weighted values. This is due to the strong tail towards higher magnetic field strengths, which reach approximately 1.07 μG and 1.73 μG in the Mach 2 and Mach 3 simulations. Figure 12 shows clearly that strong caution should be exercised when using CR electron observations to infer the average magnetic field strength in radio relics.

4.6. Impact on emission

4.6.1. Spectral index and intensity maps

We now turn our attention to the projected synchrotron emission itself, and how features in intensity and spectral index maps can be explained with the knowledge we have gained over the last sections. We first focus on the intensity maps, which we show for our Mach 3 simulation at t = 180 Myr in the upper panels of Fig. 13. The intensity of the emission is calculated for each pixel using Eq. (12). In the left-hand panel, we focus on emission at 150 MHz. The maximum values here reach approximately 6 × 10⁻³ μJy arcsec⁻². Typical radio relic surface brightnesses, meanwhile, are 0.1 − 1 μJy arcsec⁻² (van Weeren et al. 2019). This is not problematic, however, as we are not modelling first-order Fermi re-acceleration of a fossil relativistic electron population in our simulations (see Sect. 2.5). Modelling this is expected to increase the emission at Mach 3 by a factor of 10 to 100 (Pinzke et al. 2013), helping bring our simulations into line with observations.

Fig. 13. Left column: Synchrotron intensity at 150 MHz (top) and 1.5 GHz (bottom), respectively, for our fiducial Mach 3 simulation at t = 180 Myr. The projection depth is 300 kpc. Right column: Spectral index maps between 325 MHz and 150 MHz (top) and 1.5 GHz and 150 MHz (bottom), respectively. Intensity fluctuations in this region are a consequence of Mach number variations, whilst spectral index variations towards the shock front are predominantly a result of the corrugated shock front. Both intensity and spectral index maps show long striations predominantly orientated along the x axis, resulting from the magnetic filaments observed in Fig. 11. The white arrow indicates a particularly clear example caused by adiabatic compression of the post-shock magnetic field. An animated version of this figure can be found here.

It can be seen that the peak of the emission traces an approximately straight line, with the shock corrugation being visible as promontories ahead of it. This follows from our observation in the previous section that the back of the shock is relatively flat, especially when seen in projection (see Figs. 11 and C.1). The most advanced regions of the shock in Fig. 13 show emission that is a factor of 10 − 100 weaker than the main body of the relic, indicating that they could be trickier to locate in observations; if the shock is not observed edge-on, this feature could be overpowered by the emission behind it. We will investigate the full impact of orientation on such features in future work.

Regions of higher intensity can be seen at the shock front at y ≈ 130 kpc and y ≈ 270 kpc. These regions are especially evident in the 1.5 GHz observations, where cooling acts faster. Such regions are produced by the more intense emission created when the shock is strongest, as can be seen by comparing these features with the position of the highest projected Mach numbers in Fig. C.1. We note that the spatial separation of these ‘hotspots’ is approximately equivalent to the turbulent injection scale (l_inj = 150 kpc). We postpone a more detailed investigation of this correlation to the next paper in this series.

Finally, we observe that the striations towards the back of the injected region align well with the patterns seen in the synchrotron-weighted magnetic field in panel (iv) of Fig. 11. Indeed, individual features in both plots can be matched up with one another. This implies that such synchrotron ‘striations’ are the result of particularly strong magnetic field lines. Evidence for the origin of these striations can be seen in the spectral index maps, which we show in the right-hand panels of Fig. 13. We calculate the spectral index in each pixel for these maps as

$$\begin{array}{c}\hfill {\alpha }_{{\nu }_1}^{{\nu }_2}=\frac{{log}_{10}(I_{{\nu }_2}/I_{{\nu }_1})}{{log}_{10}({\nu }_2/{\nu }_1)},\end{array}$$(17)

where ν₂ > ν₁ so that α_ν1^ν₂ < 0 for a spectrum that decreases with increasing frequency.

The brightest filaments in the left-hand panels are evident as red-orange colours in the spectral index maps, indicating relatively fresh injection. Filaments with weaker emission are also evident, however, in greenish colours. It is notable that each filament shows a relatively uniform colour; the strong magnetic fields in them should lead to rapid cooling, creating a colour gradient if tracers travelled along them. Indeed, the lack of colour gradient indicates that this is not the case. We find instead that the magnetic field is amplified simultaneously along its length, temporarily brightening CR electrons of the same age. This happens in one of two main ways: (i) the filament is formed close to the shock front and is then advected downstream, typically being rotated vertically due to the gas flow, or (ii) compression and shearing at the contact discontinuity amplifies a region of the magnetic field simultaneously, as is the case for the horseshoe-like filament evident where the white arrow is located in the top left panel. Both formation mechanisms can be witnessed in the animated version of Fig. 13, linked here.

Towards the front of the injected region, it can be seen that the spectral index undergoes fluctuations. Specifically, behind the most advanced part of the shock, the spectral index temporarily reduces in value before increasing again. This feature is best seen in the ${\alpha }_{150\mathrm{MHz}}^{1.5\mathrm{GHz}}$ data, in which the effects of cooling are particularly strong. This indicates that more aged material temporarily dominates the emission, as could also be seen in the synchrotron-weighted projection of the time since injection in Fig. 11. As explained in Sect. 4.5.3, this feature is produced by the transport of older material towards the shock front due to the Rayleigh-Taylor instability (see Fig. 3). As previously, the observation of this feature will be easiest when the relic is seen edge-on.

The spectral index values in Fig. 13 approximately match the values and trends observed for the Toothbrush relic (see e.g. van Weeren et al. 2012a; Rajpurohit et al. 2020a). To quantify this statement, however, we turn to the use of colour-colour diagrams.

4.6.2. Colour-colour diagrams

Colour-colour diagrams are often used in radio observations to understand the underlying acceleration and cooling processes. They are created by plotting two spectral index maps against each other. Points on this plane then measure the spectral curvature, i.e. the functional behaviour of dlog₁₀j_ν/dlog₁₀ν in the given frequency range, which, in turn, is a measure of the spectral shape. As the CR electrons age, they trace a trajectory in the plane, indicating how the spectral shape evolves with time. This can then be compared against the trajectories produced by various models.

As discussed in Sect. 1, four models are typically used. These are abbreviated: KP, JP, CI, and KGJP²⁹. The first two of these model a one-time injection, assuming that the pitch angle distribution in each CR electron population is either fixed or continuously re-isotropised, respectively. The difference, here, is that the KP model leads to an increasingly inhomogeneously cooled population along the line of sight, leading to flatter spectral curvature at later times. This generates a curved trajectory in the colour-colour plane. In contrast, the JP model produces homogeneous populations along the line of sight, generating relatively straight lines in the colour-colour plane, indicative of approximately constant spectral curvature. The CI and KGJP models are extensions of the JP model, and assume continuous injection and injection for a limited time only, respectively. These models also produce (initially) curved trajectories in the colour-colour plane, for the same reasons as above.

In theory, such tracks are insensitive to the rate of cooling, as this only affects the rate at which the spectra develops, not its overall shape. Varying the density and magnetic field strength³⁰ should therefore only change the distance of the points from the shock, not their locus on the plane. The key assumption underpinning each model, however, is that lines of sight intersect CR electron populations that have aged for the same amount of time. This implicitly assumes an edge-on perspective and that distance from the shock front, d, is a reliable measure of time cooled, t; specifically, that d = υ_postt, where υ_post is the post-shock velocity given by the jump conditions.

In Fig. 14, we show a colour-colour diagram for our fiducial Mach 3 simulation at t = 180 Myr, created from the data presented in Fig. 13. In the left-hand panel, we show a scatter plot, where each point matches a pixel in the projected emission maps. We colour each point by its distance from the shock median. For reference, we include the cooling models fitted to the B1 region of the Toothbrush relic by Rajpurohit et al. (2020a). We do not specifically aim to replicate this radio relic in the current work, but by coincidence our simulation approximately matches the same maximum spectral index values, with both observed and simulated values reaching α ≈ −0.65. We should therefore expect such models to also fit the simulation. It is evident, however, that none of them do.

Fig. 14. Left: Colour-colour diagram created from our fiducial Mach 3 simulation at t = 180 Myr using spectral indices between 1.5 GHz and 150 MHz, and 650 MHz and 150 MHz, respectively. Points are coloured based on their distance to the shock median. The diagonal grey line indicates zero spectral curvature. The remaining curves indicate spectral cooling models fitted to the B1 region of the Toothbrush radio relic by Rajpurohit et al. (2020a, see text for details). All models begin at ${\alpha }_{150\mathrm{MHz}}^{1.5\mathrm{GHz}}={\alpha }_{150\mathrm{MHz}}^{650\mathrm{MHz}}=-0.65$. Right: As previous, but points have been binned over such that colour now represents the phase space density. Contours show the analogous data observed for the B1 region by Rajpurohit et al. (2020a) and contain 50%, 75%, and 95% of the scatter points from their figure 10 (represented by darker to lighter shades, respectively). Data from the simulation contradict traditional cooling models but are in rough agreement with the observed colour-colour diagram for a real radio relic.

In particular, the KP and CI models produce shapes fundamentally incompatible with the simulations. The KGJP model appears to fit within the scatter, but this, too, shows an initially curving trajectory, not evident in the simulated data, and is too steep once continuous injection is turned off. Moreover, as the models start at the data point with the highest spectral index, their trajectories should be compared with the bottom of the scatter, not its mid-point. With this in mind, it can be seen that the JP model is the closest to describing the data. This is perhaps to be expected, as our setup is a one-time injection with cooling independent of pitch angle (see Sect. 2.5), which matches the theoretical basis of the model. However, even the JP model has flaws; firstly the trajectory produced is too steep, and secondly, this model is unable to replicate the subtle flattening of the simulated trajectory, which begins at ${\alpha }_{150\mathrm{MHz}}^{1.5\mathrm{GHz}}\approx -1.3$.

This ‘flattening’ feature is perhaps best seen in the right-hand panel of Fig. 14, where we show the phase space density. This feature appears in the Sausage relic (see figure 19 of Stroe et al. 2013) and in the Toothbrush relic. To show this, we overlay contours indicating the phase space density observed for the B1 region of the Toothbrush relic, as given in figure 10 of Rajpurohit et al. (2020a). To produce the contours, we have binned the scatter points provided and smoothed the resulting 2D histogram with a Gaussian kernel. It can be seen that the trajectory of the observed data also flattens at approximately the same point. To help explain the origin of the feature, we present a schematic in Fig. 15.

Fig. 15. Schematic showing why the trajectory in the colour-colour plane flattens, as shown by the break in the dashed black lines in the left-hand panel, rather than following the standard JP model (dashed grey lines). Lines of sight initially intersect homogeneous CR electron populations (top right panel), with the same steep spectral curvature. Due to turbulence, later lines of sight intersect inhomogeneous populations (bottom right panel), made up of superpositions of cooled and fresher spectra (dashed purple and lilac lines, respectively). The resultant weaker spectral curvature correspondingly flattens the trajectory in the colour-colour plane. The effect is exaggerated here for explanatory purposes.

Here, the distribution of points in the colour-colour plane is represented by the two dashed, black lines in the left-hand panel. These lines represent the approximate boundaries of the distribution of spectral indices resulting from the Mach number distribution. The trajectory is initially steep, but constant, as lines of sight intersect homogeneous populations of CR electrons, which have all cooled for the same amount of time. This scenario is shown in the top right panel. As time increases, the ‘break’ frequency moves leftwards, as shown by the dashed grey line in the same panel. The spectrum thereby becomes increasingly steep when measured at the same frequencies. However, the spectral curvature, which is the relation of the logarithmic slope of the radio surface brightness measured between ν₃ and ν₁, and ν₂ and ν₁, stays the same. If this evolution continued, it would result in the continuation of the straight trajectory in the left-hand panel, shown by the light grey, dashed lines.

Instead, however, the trajectory flattens. This is due to the turbulence in our simulations, which results in lines of sight intersecting inhomogeneous populations that have cooled for different amounts of time. The result of this scenario is shown in the bottom right panel. Lines of sight are now contributed to by both cooled spectra (purple line) and by more freshly injected spectra (lilac line). This weakens the curvature of the integrated spectra, reducing the difference between α_ν1^ν₃ and α_ν1^ν₂ and hence flattening the trajectory in the colour-colour plane as well. This effect can be seen by comparing the binned spectra in Fig. 6 and Fig. D.1, which bin by age and by distance respectively. It is clear that binning by distance produces weaker spectral curvature, as we no longer bin over homogeneously cooled populations.

Under this interpretation, the trajectory will flatten when the intensity of the more freshly injected spectra contributes significantly to the overall intensity of the pixel. Whilst this effect begins in Fig. 14 at approximately the same point for both the simulated and the observed data, we note that the overall trajectory of the observed data at ${\alpha }_{150\mathrm{MHz}}^{1.5\mathrm{GHz}}<-1.3$ is even flatter than that seen in our simulation. There are multiple possibilities to explain this, including smoothing and projection effects (Rajpurohit et al. 2020a), or potentially still missing physics, such as Mach number dependent electron acceleration efficiencies, insufficient mixing due to post-shock turbulence not being strong enough, pitch-angle dependent cooling in combination with an incomplete pitch-angle isotropisation process, or turbulent re-acceleration (Fujita et al. 2015; Kang 2024). We defer a more detailed analysis of the differences to future work, where we may compare with a more like-for-like simulation.

Finally, in the right-hand panel of Fig. 14, it can be seen that both the observed and simulated data produce trajectories with similar widths. This remains fairly constant over the length of the trajectory, and will initially be set by the range of Mach numbers at the shock front, as these produce different spectral indices (see Eq. (14)). Points lying close to the dashed, grey line, which marks zero spectral curvature, have experienced minimal cooling, and should therefore be able to constrain the overall distribution. We will investigate the use of this as a diagnostic in future work as well.

5. Discussion

5.1. What affects the formation of the Rayleigh-Taylor instability?

The Rayleigh-Taylor instability has played a large role in the results shown in this paper. Specifically, in Sect. 4.5 we showed that it plays a fundamental role in amplifying the magnetic field through compression and shearing, whilst in Sect. 4.6 we showed that it breaks the commonly used assumption that d = υ_postt, thereby helping to explain observed trends in spectral features. Accordingly, it is important to understand under what circumstances the instability forms.

As explained in Sect. 4.1, in order to generate a Rayleigh-Taylor instability, we require the gradients of density and pressure to be misaligned with one another. Specifically, they must be misaligned such that the induced vorticity amplifies the perturbation. In our case, this means that the density gradient must point across the contact discontinuity towards the shock-compressed region. Furthermore, the amount of vorticity generated is proportional to the size of the gradients. Indeed, the growth rate of the instability follows

$$\begin{array}{c}\hfill \sigma =\sqrt{|\mathbf{a}|\left(\frac{{\rho }_2-{\rho }_3}{{\rho }_2+{\rho }_3}\right)k},\end{array}$$(18)

where a ∝ −∇P is the acceleration of the lighter material into the denser material, the density terms in brackets constitute the Atwood number, the subscripts ‘2’ and ‘3’ indicate post-shock gas and gas to the left of the contact discontinuity, respectively, and k is the wave number of the perturbation (Chandrasekhar 1961). In the scenarios presented in this paper, the most influential variable in this equation is the density.

The difference between ρ₂ and ρ₃ is amplified in our simulations due to the rarefaction wave produced (see Fig. 1). This wave exists in cosmological simulations as well, as can be seen in the movie of Fig. 2, but is somewhat masked by the overall radial density gradient found in the cluster. We would therefore expect a slightly slower growth of the instability in a real galaxy cluster; for example, if all other variables remained the same, but ρ₃ increased to 1 × 10⁻²⁸ g cm⁻³, we would expect the growth rate to be halved.

Slower shocks will also reduce the growth rate, as this affects a in Eq. (18). These can be generated by increasing the upstream average density or weakening the Mach number. Indeed, we have predominantly shown results from our Mach 3 simulations throughout this paper, as the growth rate here is 1.5 times faster than in the corresponding Mach 2 simulation. The turbulence is subsequently better developed in these simulations at the end of the 250 Myr runtime. Nonetheless, as we show in Fig. 12, even the Mach 2 simulations are able to produce sufficient turbulence to generate μG strength magnetic fields.

The instability will be totally suppressed when ρ₃ > ρ₂; i.e. when the density achieved due to shock compression is lower than the initial density downstream. As can be seen in the movie of Fig. 2, this is the case for the propagating merger shock wave during the initial stages of the cluster merger, before it reaches an accretion shock. Along with geometric and shock-dissipated energy arguments (Vazza et al. 2012), as well as the density argument given above, this may help explain why radio relics are so rarely observed close to the cluster centre: the mechanism presented in this paper requires conditions only available at the cluster outskirts.

5.2. How does our method for generating turbulence impact the results?

Many authors have used shock tubes to probe resolution-dependent physics. In particular, they have frequently been used to investigate how upstream turbulence affects downstream magnetic field amplification (see e.g. Ji et al. 2016; Hu et al. 2022; Hew & Federrath 2023). Few papers, however, have focused on the gas and shock properties relevant to radio relics, with the notable exception of Domínguez-Fernández et al. (2021a). In contrast to our own work, the authors of this paper used the Ornstein–Uhlenbeck method to generate upstream fluctuations. This process creates density, velocity, and magnetic turbulence simultaneously (see Federrath et al. 2010, and references therein). The turbulence then decays over time as the shock progresses through it.

The key advantage of the Ornstein–Uhlenbeck method is that turbulent properties are self-consistently modelled. Hence fluctuations are correlated, unlike in the simulations we present. Moreover, the magnetic field will display intermittency, which is not true of the Gaussian-random fields used in our initial conditions. Such intermittency is produced naturally by small-scale dynamos (Sur et al. 2021), which are expected to generate the magnetic field in the ICM (Tevlin et al. 2025). Although density turbulence dominates the effects shown in this paper, magnetic fluctuations could be influential for the striations analysed in Sect. 4.6. We will investigate this effect in a future study.

As already stated, in our own simulations, there is no initial velocity turbulence (see Sect. 2.3.4). Such turbulence should exist in the ICM, however, which will also result in Mach number variations. The velocity fluctuations would need to be a substantial fraction of the shock speed, however, before they could outweigh the effect of the density fluctuations shown in this paper. We leave a precise parameter study of this effect to future work as well.

6. Conclusions

A standard theory of radio relics has emerged in the past two decades. This states that radio relics are formed by the megaparsec-sized shocks generated during cluster mergers and that such shocks (re-)accelerate electrons, enabling them to emit in the radio band. This framework has been largely successful in explaining observational data, however, several significant challenges remain. This study has focused on three of those problems:

• (i)

  What is the origin of the ℳ_radio − ℳ_{X − ray} discrepancy if both measurement methods trace the same physical shock?

• (ii)

  How do magnetic fields reach μG strength in radio relics (as inferred from cooling length arguments) if the surrounding ICM has field strengths an order of magnitude smaller?

• (iii)

  Why do spectral variations fail to match standard cooling models in colour-colour diagrams?

To answer these problems, we took a hybrid approach. First we identified typical conditions in cosmological zoom-in simulations of cluster mergers before applying our findings to significantly higher-resolution idealised shock tubes. In this last step, we used the CR electron spectral code CREST and the emission code CRAYON+, thereby producing synchrotron data ab-initio.

From our cosmological simulations, we have identified that the merger shock typically meets an accretion shock at distances commonly recorded for radio relics. This leads to the production of a narrow, shock-compressed density sheet (Fig. 2). This scenario can be modelled as a shock tube problem (see Sect. 2.3). In addition, we included upstream density turbulence in our shock tubes, which is not well-resolved at the outskirts of our zoom-in simulations but is evident in many previous studies and observations. We chose a log-normal distribution and set the variance of the distribution according to results given in Zhuravleva et al. (2013).

We find that density fluctuations are directly responsible for the following effects:

• (i)

  A Mach number distribution forming at the shock front (Figs. 3 and 5),

• (ii)

  A Rayleigh-Taylor instability forming at the contact discontinuity (Figs. 3 and 4), and

• (iii)

  Shock corrugation (Fig. 3).

We show how these features form in Sect. 4.1. These effects provide answers to our three main questions:

• X-ray versus radio discrepancy: In Sect. 4.3, we show that the tail of the Mach number distribution dominates the higher frequencies. In particular, it flattens the integrated injected spectra and results in the slope at radio emitting frequencies being shallower than the theoretical value of α_e − 1 (Fig. 6). Both of these effects lead to an over-estimation of the Mach number with respect to the actual peak of the distribution (Table 2). The effect is particularly strong in weak shocks (ℳ ≲ 2).

• Magnetic field amplification: The Rayleigh-Taylor instability causes higher levels of amplification than would be expected due to the standard theory of shock-compression alone. This boosts the magnetic field strength from ICM-like conditions up to μG levels, although the volume-average remains significantly below this (Figs. 7 and 8). We find that the amplification is predominantly driven by the time-varying adiabatic compression and shearing that results from the instability (Figs. 9 and 10). Additionally, we find that the projected synchrotron emission is dominated by the strongest magnetic field values encountered along the line of sight. This means that estimations of the magnetic field made from emission are heavily biased by these regions; values inferred from cooling length arguments are only representative of the tail of the distribution, not the volume-average (Figs. 11 and 12).

• Cooling models: The turbulence generated by the Rayleigh-Taylor instability breaks the laminar flow assumption. This is especially true on an electron-by-electron basis, where distance from the shock front is a poor indicator of time since injection (Figs. 6 and 7). The relationship is somewhat better recovered when projecting the synchrotron-weighted time since injection, but even here d = υ_postt no longer holds true (Fig. 11). Turbulence leads to the mixing of spectra along the line of sight, which acts to flatten the spectral curvature and explains why standard cooling models do not work for radio relics (Figs. 14 and 15).

In addition, we find that features in surface brightness and spectral index maps can be attributed to specific mechanisms (Fig. 13). Specifically, the corrugation of the shock front leads to filamentary emission in projection. Mach number fluctuations, on the other hand, cause fluctuations in intensity along the shock front that spatially correlate with the injection scale. The Rayleigh-Taylor instability can produce spectral index fluctuations towards the front of the shock. Finally, magnetic fluctuations, either advected downstream, sheared, or compressed at the contact discontinuity, produce striated emission aligned with the shock normal. In the next papers in the series, we will show how these effects can be used to determine ICM conditions and will address the critical Mach number question.

Data availability

Movies are available on youtube via URLs indicated in the text. The data underlying this article will be shared on reasonable request to the corresponding author.

Acknowledgments

The authors thank Kamlesh Rajpurohit for providing data used in Fig. 13 of this paper. They thank Thomas Berlok for help with the analysis tools package PAICOS (Berlok et al. 2024) and support in accessing the PICO-CLUSTER simulations. They would also like to thank Rüdiger Pakmor, Anatoly Spitkovsky, and the organisers and participants of Galaxy Clusters & Radio Relics II for stimulating discussions and presentations, which helped improve this paper. JW acknowledges support by the German Science Foundation (DFG) under grant 444932369. CP and JW acknowledge support by the European Research Council under ERC-AdG grant PICOGAL-101019746. CP and LJ acknowledge support by the DFG Research Unit FOR-5195. PG acknowledges financial support from the European Research Council via the ERC Synergy Grant ‘ECOGAL’ (project ID 855130). The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time on the GCS Supercomputer SuperMUC-NG at the Leibniz Supercomputing Centre (LRZ) in Garching, Germany, under project pn68cu.

References

Appendix A: Numerical stability

In order to trust our results, it is critical that the simulations are numerically converged. For this purpose, we re-ran our fiducial simulation, Turb, using a target gas mass of m_gas ≈ 1.2 × 10⁵ M_⊙. This is eight times lower in mass resolution than the fiducial run or, equivalently, two times lower in spatial resolution. We find that all of our results continue to hold, although the magnetic field amplification is slightly weaker, as discussed in Sect. 4.5.2 and in App. G.

As evidence for convergence, we present a slice through the simulation at t = 180 Myr in Fig. A.1, with colours indicating density. This may be directly compared with panel (ii) in Fig. 3. It can be seen that all structures continue to be evident, except with lower resolution.

Fig. A.1. Version of our fiducial simulation, run with 8 times lower mass resolution. Colours indicate density, with the same limits and bounds as shown in Fig. 3. Higher resolution only leads to smaller-scale structure, indicating the numerical stability of our simulations.

Appendix B: The lack of a Rayleigh-Taylor instability in the Flat run

In Sect. 4.1, we claimed that density perturbations were necessary for the formation of the Rayleigh-Taylor instability in our simulations. This instability forms a critical part of our mechanism for solving the magnetic field and cooling model problems outlined in Sect. 1. To show that perturbations are necessary, we show slices through our Flat shock tube simulation in Fig. B.1. This simulation has no density turbulence in the upstream initial conditions, but includes magnetic turbulence (see Sect. 2.3.3). The panels shown are analogous to the panels in the top row of Fig. 3 and are shown at the same time (i.e. t = 180 Myr).

Fig. B.1. As Fig. 3, but only the top row is shown, and data are taken from the Flat simulation. Without upstream density fluctuations, the Rayleigh-Taylor instability is restricted to a very low growth rate.

It can be seen that without upstream density turbulence, the shock front shows no curvature. Magnetic pressure fluctuations, however, lead to mild variations in density and pressure gradients at the contact discontinuity. Consequently, the Rayleigh-Taylor instability still takes place here; indeed, this forms the bumps seen in the middle panel at the trailing edge of the shock-compressed zone. The growth-rate of this instability has been severely restricted, however. This is partially due to the significantly smaller pressure gradient, and partially due to the lack of corrugation at the contact discontinuity, which helped to seed the Rayleigh-Taylor instability in Fig. 3.

Appendix C: Mach numbers in projection

In Sect. 4.2, we argued that lower Mach numbers should be found typically in a more advanced position, and that stronger Mach numbers should be found towards the back of the shock. In panel (iii) of Fig. 3, specifically, we analysed a thin projection and found the result to be consistent. In Fig. C.1, we show the projected dissipated-energy weighted Mach number in our Mach 3 Turb simulation at t = 180 Myr. This is analogous to the previous projection, except the projection depth has increased from 35 kpc to 300 kpc (i.e. the box size). It can be seen that our argument continues to hold.

Fig. D.1. Projected dissipation-weighted Mach number, shown for the same simulation and time as in Fig. 3. The projection depth is 300 kpc. Weaker Mach numbers tend to form where the shock advances ahead of the median position, whilst stronger ones form where it lags behind.

We have already shown in Sect. 4.2 that the properties of the Mach distribution remain relatively stable over time. With this said, however, we must also emphasise that our argument only holds true in a statistical sense; Mach numbers do not decrease monotonically in value towards the most advanced part of the shock. This can also be seen in Fig. C.1.

Appendix D: Spectra binned by distance

In Sect. 4.4 we presented evidence that distance from the shock is not a reliable proxy for time since injection, at least on a tracer-by-tracer basis (see caveats given in Sect. 4.5.3). We showed, in particular, that the volume-averaged spectrum is created through the layering of tracers binned by age, and that tracers in these bins overlap significantly when their spatial distributions are overplotted (see Fig. 6). In Fig. D.1, we plot the inverse; that is, we bin tracers by their distance and make histograms of the temporal distributions. To remove some degree of arbitrariness, we use the same time bins as previously, but now multiply out by the theoretical post-shock speed υ_post. We have calculated this by taking the actual shock speeds of approximately 1, 000 km s⁻¹ and 1, 500 km s⁻¹ in the Mach 2 and Mach 3 simulations, respectively, and applying the standard jump conditions, giving υ_post ≈ 440 km s⁻¹ and υ_post ≈ 500 km s⁻¹. This equates to distance intervals of approximately 0, 0.04, 4, 18, 46, and 107 kpc in the Mach 2 simulations, and intervals 1.14 times higher in the Mach 3 simulations.

The layering in Fig. D.1 is clearly less well-defined than in Fig. 6. Nonetheless, it still remains true that CR electrons further from the shock are, on average, more cooled. Consequently, it is still possible to pick distance intervals that dominate a particular momentum range. This is especially true at later ages, when the cooled part of the spectra is steeper. This is despite the fact that the tracers are strongly mixed, as can be seen by inspecting the bottom row of Fig. D.1.

The layering of spectra can be partially improved by removing the assumption that the distance bins in Fig. D.1 must be proportional to the time bins given in Fig. 6. However, in order to see an improvement, the spacing of bins must be picked independently for Mach 2 and Mach 3 simulations. Moreover, we find that, regardless of binning, it is not possible to produce a figure entirely analogous to Fig. 6. This implies that the standard formula of d = υ_postt cannot be adapted in a linear fashion, and that laminar flow is only a good assumption for simulations without density turbulence.

Appendix E: Plasma beta distribution

In Fig. E.1 we quantify our earlier statement in Sect. 4.5.1 that despite significant amplification, the magnetic field does not play a major dynamical role. To do this, we show histograms of the plasma beta values in the initial upstream region (region III, in the nomenclature given in Sect. 2.3.1) at t = 0 Myr and in the shock-compressed region at t = 250 Myr. As discussed in Sect. 2.3.4, we implemented density and magnetic field turbulence in our simulations independently of one another. This leads to a wide initial distribution with a strong tail towards higher plasma beta values. The RMS magnetic field strength was picked, however, to produce a peak at P_th/P_B = 100, as can be seen in the figure.

After the gas has undergone shock compression, the width of the distribution in Fig. E.1 increases. A minority of cells reach values P_th/P_B < 10, but no cell reaches below P_th/P_B = 4. We may therefore safely treat the downstream region as being in the kinetic regime, meaning that hydrodynamic models of the downstream remain valid.

Fig. E.1. As Fig. 6, except tracers are now binned by distance from the median shock position, rather than time since injection. Distance bins are equal to the time bins in Fig. 6 multiplied by the relevant post-shock velocity (approximately 440 and 500 km/s for Mach 2 and Mach 3 shocks, respectively). Contributions binned by distance do not match those binned by time since injection due to the highly non-laminar flow. This can also be seen in the strongly overlapping distributions of time since injection within each distance bin (bottom row).

Appendix F: Magnetic field decay

In Sect. 4.5, we claimed that the reduced scaling of the magnetic field with density in the Flat-TurbB simulation was due to the decay of the magnetic field. We also claimed that, in the absence of additional amplification mechanisms, this was an expected physical process in the post-shock region. We justify this claim here.

In Fig. F.1, we show the RMS magnetic field strength over time, calculating the average using all tracers that passed the shock before t = 10 Myr. It can be seen that the solid lilac line starts at the upstream average value of 0.16 μG before being strongly amplified by the next snapshot³¹. Shortly thereafter, however, the magnetic field strength enters a prolonged decline. This decline is naturally explained by the untangling of the magnetic field: a tangled magnetic field feels the effects of magnetic tension, which is communicated at the Alfvén speed, υ_A. As a result, the magnetic field will untangle itself over a timescale of τ = l_inj,  B/υ_A, where l_inj,  B is the outer scale of the magnetic turbulence. Without a process to keep the field amplified, this leads to the overall magnetic energy decaying as E_B ∝ exp( − t/τ), and hence the magnetic field strength decays as B ∝ exp( − 0.5 t/τ).

Fig. F.1. Histograms of the plasma beta values for all gas cells upstream at t = 0 Myr (dashed) and in the shock-injected region at t = 250 Myr (solid) for our Mach 3 fiducial simulation, weighted by the cell volume. The distribution initially peaks at 100 (marked by a dotted, grey line). Amplification downstream is able to increase the magnetic field pressure up to ∼10% of the gas pressure in a fraction of the cells.

In our initial conditions, we set l_inj,  B = 40 kpc. However, post-shock, the field is compressed along the x-direction by the shock compression ratio, x_s. For a Mach 3 shock, the effective injection scale hence becomes ${\stackrel{\sim }{l}}_{\mathrm{inj},B}=40/3$ kpc. This leads to a model for the observed decay of the field strength:

$$\begin{array}{c}\hfill B=B_{2}exp(-0.5\frac{{\upsilon }_{\mathrm{A}}(t)}{{\stackrel{\sim }{l}}_{\mathrm{inj},B}}t),\end{array}$$(F.1)

where B₂ is the peak magnetic field strength reached post-shock.

We plot this as the dotted lilac line in Fig. F.1. To calculate it, we keep ${\stackrel{\sim }{l}}_{\mathrm{inj},B}$ fixed, but set υ_A(t) to the magnetic-energy weighted mean Alfvén speed calculated for the same group of tracers as used for the solid line. We also renormalise B₂ to the time-averaged post-shock value. It can be seen that this simple model does a remarkably good job of predicting the decay of the field strength.

As the magnetic field becomes weaker, the decay slows down. This is because ${\upsilon }_{\mathrm{A}}=|B|/\sqrt{4\pi \rho }$, where ρ is the gas density. Since the post-shock gas density is effectively a constant in the Flat simulations, a weaker magnetic field results in a slower Alfvén speed³² and hence a longer Alfvén timescale on the outer scale of magnetic turbulence. We plot this quantity in orange in Fig. F.1.

We note that the timescale ranges from ∼150−250 Myr, whilst the system evolves for a comparable amount of time. This explains why we see magnetic decay in the first place; were the timescale much greater, there would be no time for the decay to set in. In contrast, in the upstream, the Alfvén timescale is closer to 450 Myr, owing predominantly to the larger value of l_inj,  B.

Finally, we compute the average post-amplification RMS magnetic field strength in Fig. F.1 to be 0.27 μG. This means that the effective magnetic field increase is by only a factor of 1.7, and hence we recover a scaling of B ∝ n^α_B, with α_B = log(1.7)/log(3)≈0.5, as observed in Fig. 9.

Fig. F.2. RMS magnetic field strength (solid lilac) and the effective Alfvén timescale (solid orange) as a function of time for tracers in the Flat-TurbB simulation. Averages are calculated using all tracers that passed the shock before t = 10 Myr. The magnetic field initially experiences the expected level of amplification, but then subsequently decays. This decay is well-described by a simple model (Eq. F.1), which replicates the untangling of the turbulent magnetic field (dotted lilac).

Fig. F.3. As Fig. 9, except solid lines now show data from our lower resolution Mach 3 fiducial simulation. The faint, dashed lines behind indicate data from the original, higher resolution version. Dashed, black lines show how the 75% percentile contour scales. It can be seen that the scalings remain the same as in Fig. 9. This is evidence against the existence of a dynamo in our simulations.

Appendix G: Magnetic amplification at lower resolution

In Sect. 4.5.2, we claimed that a magnetic dynamo was unlikely to be active in our simulations. We show evidence for this in Fig. G.1. Here, we show how density and magnetic field strength scales for the lower resolution version of our Mach 3 fiducial simulation, as initially introduced in App. A. For comparison, we show data from the original, higher resolution simulation behind in fainter, dashed lines.

It is evident that, whilst the lower resolution simulation does not quite reach the same peak magnetic field strength values, the empirical scalings are identical. In particular, any increase in magnetic field strength is compensated for by an additional increase in density. For subsonic turbulence, however, a magnetic dynamo should also result in additional amplification at the same density, thereby increasing α_B, where B ∝ n^α_B. We have repeated this analysis on the individual magnetic field components and find the same behaviour here as well. We therefore conclude that a dynamo is not relevant for magnetic field amplification in our simulations.

All Tables

All Figures

Fig. 1. Top: Mean pressure along the x axis at t = 0 Myr (dashed) and t = 100 Myr (solid) in our Mach 3 ‘Flat’ simulation (see Table 1). Dotted lines indicate the edges of various regions in the shock tube (see Sect. 2.3.1). We greyed out panels that we do not analyse. Bottom: As above, but lines show the mean density. The initial pressure discontinuity causes a shock to propagate into region III. As a result of shock compression, a dense sheet moves rightwards, bounded by the density jumps at the contact discontinuity and the shock front. This feature is critical to the mechanism we analyse in this paper.

In the text

Fig. 2. Cosmological simulation of a galaxy cluster undergoing a major merger at z = 0.14. Four left-most panels, clockwise from top left: Projected shock-dissipated energy rate, gas pressure, gas density, and dissipation-weighted Mach number, respectively. Projections have a depth of ±7.5 Mpc from the cluster centre. Right-most panels: enlarged cut-outs showing slices through the midplane of the projection. The white, dashed box indicates the size of the upstream region in our idealised shock tube simulations. The merger drives a shock wave out to the cluster outskirts. The collision of this shock wave with an accretion shock leads to the production of a thin shell of compressed gas (see white arrow). A movie of this process can be found here.

In the text

Fig. 3. Mach 3 shock is driven into a turbulent upstream density field. The initial pressure and density values for the shock tube correspond to those shown in Fig. 2. Panels: (i) gas pressure, (ii) gas density, (iii) dissipated-energy-weighted Mach number, (iv) the fraction of the gas speed in the x-component, (v) the x-component of gas velocity, (vi) gas speed divided by sound speed. Data are shown in slices except for (iii), which is a thin projection of 35 kpc. All values are shown in the shock rest frame at t = 180 Myr. Upstream density turbulence directly leads to shock corrugation, a distribution of Mach numbers at the shock front, and downstream velocity turbulence. An animated version of this figure can be found here.

In the text

Fig. 4. Schematic showing how upstream density turbulence causes velocity turbulence to be generated downstream. The corrugation of the shock front naturally leads to a misalignment of pressure and density gradients (shown here by red and green arrows, respectively). This results in a baroclinic term, which in turn induces vorticity, causing the contact discontinuity to become Rayleigh-Taylor unstable.

In the text

Fig. 5. Mach number distributions for all models, where each cell has been weighted by its normalised contribution to the shock surface. Lines indicate the median taken over all snapshots, whilst the shaded values indicate the interquartile range. The addition of magnetic turbulence (Flat) to a homogeneous density distribution (Flat-ConstB) broadens the distribution only very mildly. By contrast, the addition of upstream density fluctuations in the simulation (Turb and Turb-ConstB) turns an extremely narrow distribution into a much broader one. The width of the distribution is proportional to the peak Mach number.

In the text

Fig. 6. Top row: Volume-weighted non-thermal electron spectra generated from our Turb (solid black) and Flat (dashed grey) simulations at t = 250 Myr. Blue lines show contributions to the Turb spectrum, where tracers have been binned by time since injection. Bottom row: Histograms indicating the total number of injected tracers relative to the shock front, where bins have a width of 2 kpc. Blue and grey colours represent our Turb and Flat simulations respectively, whilst colour saturations indicate the same time bins as above. The broadened Mach number distribution seen in Fig. 5 causes a shallower slope, αe, in the Turb runs compared to the theoretical expectation. This is especially noticeable for low Mach number shocks. In the case of upstream density fluctuations, substantial mixing takes place downstream, such that distance from the shock front is no longer a good indication of cooling time.

In the text

Fig. 7. Slices through our fiducial Mach 3 simulation at t = 180 Myr showing: (i) time since injection, (ii) CR electron number density at p = 2 × 104 (see text), (iii) magnetic field strength, and (iv) synchrotron emissivity at 150 MHz. Only tracers that have undergone DSA are shown. A Rayleigh-Taylor instability induces a counter-streaming plume (in the shock rest frame) that brings aged, and therefore cooled, electrons close to the shock front. However, these electrons are still relatively bright at 150 MHz due to magnetic field amplification up to μG levels. An animated version of this figure can be found here.

In the text

	Fig. 8. Slices through our fiducial Mach 3 simulation at t = 180 Myr showing plasma beta values. The grey lines mark the corrugated shock surface and the contact discontinuity (i.e. the region within which tracers have been injected). Despite significant amplification, beta values are typically well above 10, limiting the ability of the magnetic field to affect dynamics.
In the text

Fig. 9. Left: Phase space diagram of magnetic field strength versus gas number density for our Mach 3 Flat-TurbB simulation. Contours cover 25%, 50%, 75%, and 95% of the probability density (represented by darker to lighter shades, respectively). Blue colours indicate the initial upstream distribution, whilst red colours show the final state in the injected region at t = 250 Myr. A cross marks the median of the distribution. Right: As previous, except data are taken from the Mach 3 fiducial simulation. Dashed lines indicate the empirical scaling of the 75% boundaries. Amplification is weaker than the naive expectation for Flat-TurbB due to magnetic decay. Meanwhile, for the fiducial simulation, amplification is significantly stronger than that possible from compression alone.

In the text

	Fig. 10. Left: As Fig. 9 but now showing individual components of the magnetic field. All components are more highly amplified in the fiducial simulation owing to time-varying compression. It is particularly notable, however, that the x component is the most highly amplified. This is indicative of amplification by shearing.
In the text

Fig. 11. Projections with a depth of 300 kpc through our fiducial Mach 3 simulation at t = 180 Myr. These show: i) volume-weighted time since injection, ii) synchrotron-weighted time since injection, iii) volume-weighted magnetic field strength, and iv) synchrotron-weighted magnetic field strength. Regions with no injected tracers have been masked (see text for further details). Synchrotron emission is calculated at ν = 150 MHz and is biased towards fresher injection and stronger magnetic fields. Observations made through this channel consequently appear to show a relatively laminar flow and μG-strength magnetic fields, even though this is not typical of the downstream.

In the text

Fig. 12. Left: Probability density distributions of the projected magnetic field strength in our fiducial Mach 2 simulation. Solid (dashed) lines represent projections through the shock-compressed region at t = 250 Myr (high-resolution region at t = 0 Myr). Blues lines represent volume-weighting, whilst orange lines represent synchrotron-weighting, where the synchrotron emission is calculated at ν = 150 MHz. Arrows are placed at the RMS magnetic field strength calculated for each distribution. Right: As previous, except data come from the fiducial Mach 3 simulation. Synchrotron-weighting significantly overestimates the average magnetic field strength. Magnetic field values are able to reach μG values even in weaker shocks, although amplification is more effective at higher Mach number shocks.

In the text

Fig. 13. Left column: Synchrotron intensity at 150 MHz (top) and 1.5 GHz (bottom), respectively, for our fiducial Mach 3 simulation at t = 180 Myr. The projection depth is 300 kpc. Right column: Spectral index maps between 325 MHz and 150 MHz (top) and 1.5 GHz and 150 MHz (bottom), respectively. Intensity fluctuations in this region are a consequence of Mach number variations, whilst spectral index variations towards the shock front are predominantly a result of the corrugated shock front. Both intensity and spectral index maps show long striations predominantly orientated along the x axis, resulting from the magnetic filaments observed in Fig. 11. The white arrow indicates a particularly clear example caused by adiabatic compression of the post-shock magnetic field. An animated version of this figure can be found here.

In the text

Fig. 14. Left: Colour-colour diagram created from our fiducial Mach 3 simulation at t = 180 Myr using spectral indices between 1.5 GHz and 150 MHz, and 650 MHz and 150 MHz, respectively. Points are coloured based on their distance to the shock median. The diagonal grey line indicates zero spectral curvature. The remaining curves indicate spectral cooling models fitted to the B1 region of the Toothbrush radio relic by Rajpurohit et al. (2020a, see text for details). All models begin at ${\alpha }_{150\mathrm{MHz}}^{1.5\mathrm{GHz}}={\alpha }_{150\mathrm{MHz}}^{650\mathrm{MHz}}=-0.65$. Right: As previous, but points have been binned over such that colour now represents the phase space density. Contours show the analogous data observed for the B1 region by Rajpurohit et al. (2020a) and contain 50%, 75%, and 95% of the scatter points from their figure 10 (represented by darker to lighter shades, respectively). Data from the simulation contradict traditional cooling models but are in rough agreement with the observed colour-colour diagram for a real radio relic.

In the text

Fig. 15. Schematic showing why the trajectory in the colour-colour plane flattens, as shown by the break in the dashed black lines in the left-hand panel, rather than following the standard JP model (dashed grey lines). Lines of sight initially intersect homogeneous CR electron populations (top right panel), with the same steep spectral curvature. Due to turbulence, later lines of sight intersect inhomogeneous populations (bottom right panel), made up of superpositions of cooled and fresher spectra (dashed purple and lilac lines, respectively). The resultant weaker spectral curvature correspondingly flattens the trajectory in the colour-colour plane. The effect is exaggerated here for explanatory purposes.

In the text

	Fig. A.1. Version of our fiducial simulation, run with 8 times lower mass resolution. Colours indicate density, with the same limits and bounds as shown in Fig. 3. Higher resolution only leads to smaller-scale structure, indicating the numerical stability of our simulations.
In the text

	Fig. B.1. As Fig. 3, but only the top row is shown, and data are taken from the Flat simulation. Without upstream density fluctuations, the Rayleigh-Taylor instability is restricted to a very low growth rate.
In the text

	Fig. D.1. Projected dissipation-weighted Mach number, shown for the same simulation and time as in Fig. 3. The projection depth is 300 kpc. Weaker Mach numbers tend to form where the shock advances ahead of the median position, whilst stronger ones form where it lags behind.
In the text

Fig. E.1. As Fig. 6, except tracers are now binned by distance from the median shock position, rather than time since injection. Distance bins are equal to the time bins in Fig. 6 multiplied by the relevant post-shock velocity (approximately 440 and 500 km/s for Mach 2 and Mach 3 shocks, respectively). Contributions binned by distance do not match those binned by time since injection due to the highly non-laminar flow. This can also be seen in the strongly overlapping distributions of time since injection within each distance bin (bottom row).

In the text

Fig. F.1. Histograms of the plasma beta values for all gas cells upstream at t = 0 Myr (dashed) and in the shock-injected region at t = 250 Myr (solid) for our Mach 3 fiducial simulation, weighted by the cell volume. The distribution initially peaks at 100 (marked by a dotted, grey line). Amplification downstream is able to increase the magnetic field pressure up to ∼10% of the gas pressure in a fraction of the cells.

In the text

Fig. F.2. RMS magnetic field strength (solid lilac) and the effective Alfvén timescale (solid orange) as a function of time for tracers in the Flat-TurbB simulation. Averages are calculated using all tracers that passed the shock before t = 10 Myr. The magnetic field initially experiences the expected level of amplification, but then subsequently decays. This decay is well-described by a simple model (Eq. F.1), which replicates the untangling of the turbulent magnetic field (dotted lilac).

In the text

	Fig. F.3. As Fig. 9, except solid lines now show data from our lower resolution Mach 3 fiducial simulation. The faint, dashed lines behind indicate data from the original, higher resolution version. Dashed, black lines show how the 75% percentile contour scales. It can be seen that the scalings remain the same as in Fig. 9. This is evidence against the existence of a dynamo in our simulations.
In the text