Morphology of radio relics I: What causes the substructure of synchrotron emission?
P. Domínguez-Fernández, M. Brüggen, F. Vazza, W. E. Banda-Barragán, K. Rajpurohit, A. Mignone, D. Mukherjee, B. Vaidya
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 29 September 2020 (MN L A TEX style file v2.2)
Morphology of radio relics I: What causes the substructure ofsynchrotron emission?
P. Dom´ınguez-Fern´andez (cid:63) , M. Br ¨uggen , F. Vazza , , , W. E. Banda-Barrag´an ,K. Rajpurohit , , A. Mignone , D. Mukherjee , B. Vaidya Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany Dipartimento di Fisica e Astronomia, Universit´a di Bologna, Via Gobetti 92/3, 40121, Bologna, Italy Istituto di Radio Astronomia, INAF, Via Gobetti 101, 40121 Bologna, Italy Dipartimento di Fisica, Universit`a di Torino, via Pietro Giuria 1, 10125 Torino, Italy Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Pune - 411007, India Discipline of Astronomy, Astrophysics and Space Engineering, Indian Institute of Technology Indore, Khandwa Road, Simrol, Indore 453552, India
Received / Accepted
ABSTRACT
High-resolution radio observations of cluster radio relics often show complex spatial and spec-tral features. However, it is not clear what these features reveal about the underlying magneticfield properties. We performed three-dimensional magneto-hydrodynamical simulations ofmerger shock waves propagating through a magnetised, turbulent intracluster medium. Ourmodel includes the diffusive shock acceleration of cosmic-ray electrons, their spatial advec-tion and energy losses at run-time. With this set-up we can investigate the relation betweenradio substructure and pre-shock plasma conditions in the host cluster. We find that upstreamturbulence plays a major role in shaping the properties of radio relics produced downstream.Within the assumption of diffusive shock acceleration, we can reproduce the observed dis-crepancy between the X-ray derived Mach number of shocks, and the Mach number inferredfrom radio spectra. Our simulated spectral index maps and profiles across the radio relic alsosuggest that the standard deviation of the upstream magnetic field must be relatively small( σ B (cid:54) µ G) in order to reproduce observations and therefore, radio relics can potentiallyconstrain the distribution of magnetic fields in galaxy clusters outskirts.
Key words: galaxy: clusters, general – methods: numerical – intergalactic medium – accel-eration of particles
Radio observations of galaxy clusters reveal non-thermal processesin the intracluster medium (ICM). The size of diffuse radio sourcesis of the order of Mpc. Such large-scale emission is typicallygrouped into two categories: radio haloes and radio relics . Thefirst category refers to objects located at the cluster center and ofa relatively regular and spherical shape, with little to absent polari-sation (e.g. Brunetti & Jones 2014). The second category, refers toobjects located at the cluster periphery, with elongated shapes andtypically large degrees of polarisation (see Br¨uggen et al. 2012 andvan Weeren et al. 2019 for reviews). Radio relics are the main focusof this work.The particle acceleration mechanisms causing diffuse radiosources are not fully understood. However, it seems clear thatshocks generated during the assembly of galaxy clusters play akey role in accelerating the synchrotron-emitting cosmic-ray elec-trons (see van Weeren et al. 2011 or Bykov et al. 2019 for a re- (cid:63)
E-mail: [email protected] view). The radio emission observed in relics is compatible withsynchrotron emission from cosmic-ray (CR) electrons with Lorentzfactors of γ e ∼ – , which are believed to be accelerated viamild shocks ( M radio1 ∼ . – . ) crossing the ICM (e.g. Clarke &Ensslin 2006; van Weeren et al. 2010, 2012).The first-order Fermi acceleration process, commonly knownas diffusive shock acceleration (DSA), explains the acceleration ofrelativistic particles by the passage of a collisionless shock wave(e.g. Blandford & Ostriker 1978; Drury 1983). A small fraction( (cid:28) − ) of thermal particles can cross the shock front multi-ple times, and receive a boost in their momentum proportional to ∆ v/c , where ∆ v is the velocity jump (difference) across the shock.This acceleration mechanism is observed to be much more efficientthan what is expected from theory (see van Weeren et al. 2019; Bot-teon et al. 2020). Hence, it has been proposed that the electrons arepre-accelerated (e.g. Kang et al. 2012a; Pinzke et al. 2013) beforethey enter the DSA mechanism. M radio is the Mach number inferred from radio observations.c (cid:13) a r X i v : . [ a s t r o - ph . H E ] S e p P. Dom´ınguez-Fern´andez et al.
Moreover, the DSA process does not offer a straightforwardexplanation for the non-detection of γ -rays from clusters (see Ack-ermann et al. 2010, 2014, 2016 and Vazza & Br¨uggen 2014).Among the most relevant open challenges in our understandingof radio relics are: (i) the discrepancy between the Mach numbersdetected in X-rays and those inferred from radio observations as-suming DSA (e.g. Botteon et al. 2020), and (ii) the high electronacceleration efficiency of the order of several percent for the weakshocks commonly associated with radio relics. Up until now, onlya few radio relics can readily be explained by the DSA model (e.g.Locatelli et al. 2020).Recent high-resolution radio observations have shown aplethora of complex structures in radio relics (e.g. Rajpurohit et al.2018, 2020; Owen et al. 2014; van Weeren et al. 2017; Di Gennaroet al. 2018). Attempts to explain the observed features struggle withthe vast range of scales, from cosmological scales ( (cid:38) Mpc), downto turbulent scales ( ∼ kpc) (e.g. Egan et al. 2016; Dom´ınguez-Fern´andez et al. 2019), or even down to plasma scales where par-ticle acceleration occurs ( ∼ − kpc for the largest gyroradiusof relativistic protons). A possible choice is to study particle ac-celeration from the microphysical point of view using Particle InCell (PIC) simulations (e.g. Guo et al. 2014; Caprioli & Spitkovsky2014a; Park et al. 2015; Caprioli & Spitkovsky 2014b; Ryu et al.2019; Kang et al. 2019), or conversely on larger scales using themagneto-hydrodynamical (MHD) approximation, as customarilydone with cosmological simulations (the reader may refer to Don-nert et al. 2018 for a review).Radio emission from radio relics has been modelled in previ-ous works on larger scales (e.g. Skillman et al. 2013; Hong et al.2015; Nuza et al. 2017). Due to the discrepancy in scales, it is notpossible for MHD simulations to model the emission produced bysingle electrons, but rather from an ensemble of tracer particles ,representing a whole distribution of electrons. Previously, this ap-proach has been applied to cosmological MHD simulations in post-processing (see Wittor et al. 2019).In this paper, we model the synchrotron emission at run-timein a small fraction of the ICM by means of a new hybrid parti-cle and fluid framework using the MHD code PLUTO (Mignoneet al. 2007; Vaidya et al. 2018). Our aim is to study the substructureobserved in radio relics (e.g. Rajpurohit et al. 2020). This methoduses Lagrangian particles embedded in a large-scale MHD flow,each with its individual energy spectrum. Here we consider a sim-plified scenario: we set up a shock tube in a turbulent medium thatis representative of a small region of the ICM. We then assume thatCR electrons are injected instantly at the shock discontinuity andacquire spectral properties according to DSA theory.The paper is structured as follows: in Section 2, we describeour numerical set-up and initial conditions. In Section 3, we includea description of the particles’ initial spectral distribution and evo-lution. In Section 3.3, we explain how we obtain the emission andspectral index maps. Section 4 shows our results and we summarizein Section 5.
The turbulent ICM initial conditions were produced using the MHDFLASH code version 4.6.1 (Fryxell et al. 2000; Calder et al. 2002),with the goal of generating realistic pre-shock conditions. We used the unsplit staggered mesh (USM) MHD solver which uses a con-strained transport (CT) method at cell interfaces for preserving thedivergence-free magnetic field property on a staggered grid (e.g.Lee et al. 2009). The simulation domain is chosen to be a cubicbox of size L = L x = L y = L z = 200 kpc, uniformly spacedover a cells grid, with periodic boundary conditions. We as-sumed an ideal gas equation of state with γ = 5 / .Turbulence was created following a spectral forcing method,based on the stochastic Ornstein-Uhlenbeck (OU) process to modelthe turbulent acceleration f with a finite autocorrelation time (e.g.Eswaran & Pope 1988; Schmidt et al. 2006; Federrath et al. 2010).The OU process describes the evolution of the forcing term inFourier space, ˆ f , with a stochastic differential equation: d ˆ f ( k , t ) = f ( k ) P ζ ( k ) d W ( t ) − ˆ f ( k , t ) dtT , (1)where f is the forcing amplitude, W ( t ) is a random Wiener pro-cess, P ζ is a projection tensor in Fourier space, ζ is the forcing pa-rameter ( ζ = 0(1) for purely compressive(solenoidal) forcing), andT is the autocorrelation time scale of the forcing (the reader may re-fer to Federrath et al. 2010 for a detailed explanation on turbulenceforcing). In this work, we solely focus on solenoidal subsonic tur-bulence forcing ( ∇ · f = 0 ), since several authors have shownthat the most dominant type of turbulence in the ICM should besubsonic with a large ( (cid:62) ) predominance of solenoidal modes(e.g. Miniati 2014; Vazza et al. 2017).We have run two simulations with two different stirring scales.The forcing amplitude, f , was set to be a paraboloid in Fourierspace in both simulations only containing power on the largestscales in a small interval of wavenumbers. We chose two dif-ferent intervals: (cid:54) kL/ π (cid:54) for the first simulation and (cid:54) kL/ π (cid:54) for the second simulation. As the power peaksin 2/3 and 1/4 of the box, we will refer to each of them as L/ and L/ , representing injection scales of
133 kpc and
50 kpc , re-spectively. Furthermore, the maximum k allowed in the box cor-responds to the 2L/3 case. Conversely, the L/ case satisfies theminimum condition where its eddy turnover time is smaller thanthe time needed for our main set-up (where a shock sweeps thisturbulent medium).The autocorrelation timescale was set equal to the dynami-cal timescale (also called eddy turn-over timescale) on the scale ofenergy injection, t L/ = 2 L/ σ v and t L/ = L/ σ v , respec-tively, where σ v is the rms velocity to be achieved at saturation.Both simulations were set to have an amplitude of the fluctuationsof σ v = 125 km/s. The turbulence is fully developed after roughlytwo dynamical timescales when the magnetic energy, the plasmabeta, and the rms Mach number become stable (although some tran-sient fluctuations can still be present depending on the balance be-tween mechanical energy from the forcing term and the dissipationrate). At this point, the magnetic field reaches a saturated state sincewe start with a relatively strong magnetic field seed of 0.2 µ G (forthe L/ case) and 0.4 µ G (for the L/ case). This is shown inFig. 1 where we plot the evolution of the total kinetic and magneticenergy and in Fig. 2, where we show the evolution of the rms Machnumber and the plasma beta, β , for both runs. Even after the mag-netic saturation, the thermal pressure fluctuations will continue toincrease due to turbulent dissipation and, as a consequence, the rmssonic Mach number decreases.We selected one snapshot from each of these two runs to rep-resent a small region of the ICM and each act as an initial con-dition (see Sec. 2 and Fig. 5). The snapshots are taken at the re-spective saturation times, which is t = 2 . Gyr for the run with c (cid:13) , 000–000 orphology of radio relics Figure 1.
Energy evolution for the runs with injection at L/ ( blue ) andat L/ ( green ). The kinetic energy is shown with dashed lines and the mag-netic energy with solid lines . Saturation is reached at t = 2 t L/ = t = 2 t L/ = E M /E K , is ∼ . and ∼ . for the 2L/3 and L/ cases, respectively. M a c h nu m b e r Time [Gyr] P l a s m a b e t a Figure 2.
Top panel: Evolution of rms Mach number for the runs with in-jection at L/ ( blue ) and at L/ ( green ). The arrows point to the selectedsnapshots to be our initial conditions. Bottom panel: Corresponding evo-lution of the plasma beta, β . The selected snapshots have an rms Machnumber of M s ∼ β ∼ Figure 3.
PDFs of magnetic field strength at the selected time for the runswith injection at L/ ( blue ) and at L/ ( green ). Extra dashed line: PDFof a post-merger (PM) galaxy cluster previously analysed in Dom´ınguez-Fern´andez et al. 2019. injection scale L/ and t = 0 . Gyr for the run with injectionscale L/ . In Fig. 3, we show the probability distribution functionof the magnetic field strength at those times. The two snapshotshave a sufficiently different distribution of magnetic fields (due tothe difference in the initial magnetic field seed and the interval ofwavenumbers of the stirring), yet in terms of the rms sonic Machnumber and plasma- β they lie in the range of characteristic valuesof the ICM (e.g. Ryu et al. 2008). Note that the tail of the L/ dis-tribution extends to ∼ times larger magnetic field values. (seeFig. 3). For reference we include the PDFs of one galaxy clus-ter previously produced in a cosmological MHD simulation andanalysed in Vazza et al. 2018; Dom´ınguez-Fern´andez et al. 2019.This cluster is classified as a post-merger (PM) cluster with a mass M = 1 . × M (cid:12) . We selected a region at a distance of ∼ Mpc from the cluster’s center with an extent of ∼ kpc. As canbe seen in Fig. 3, the PDF of the magnetic field broadly agrees withprevious results from cosmological MHD simulations. The high-magnetic field tail of the distribution is slightly more extended thanin cosmological simulations, owing to the larger Reynolds num-ber in these new simulations. The magnetic field strength at theoutskirts of the clusters is most likely underestimated in the cos-mological simulation due to the limited resolution (see Vazza et al.2018). A more extensive survey of the interaction between mergershocks with a larger range of different initial turbulent conditionsfor the ICM will be the subject of future work. In order to study the synchrotron emission in an MHD shock tube,we use the code PLUTO (Mignone et al. 2007), which solves the The reader can find this galaxy cluster in Dom´ınguez-Fern´andez et al.2019 with ID E16Bc (cid:13) , 000–000
P. Dom´ınguez-Fern´andez et al.
Figure 4.
Initial magnetic field configuration in the PLUTO code. Thestreamlines are coloured according to the magnitude of the magnetic field:green, light colours denote higher values, while dark blue colour indicateslower values. I denotes the post-shock region, II denotes the uniform pre-shock region and III the turbulent pre-shock region (see Sec. 2.1). The leftside is a uniform medium with a B x component matching the mean valueof the B x of the turbulent medium. We have one Lagrangian particle percell placed in the whole regions II and III. following conservation laws for ideal MHD: ∂ρ∂t + ∇ · ( ρ v ) =0 , (2) ∂ m ∂t + ∇ · (cid:20) mv − BB + I (cid:18) p + B (cid:19)(cid:21) T = 0 , (3) ∂ ( E t ) ∂t + ∇ · (cid:20)(cid:18) ρv ρe + p (cid:19) v − ( v × B ) × B (cid:21) =0 , (4) ∂ B ∂t − ∇ × ( v × B ) =0 , (5) ∇ · B =0 , (6) ρe = pγ − , (7) where ρ is the gas mass density, m = ρ v is the momentum density, p is the thermal pressure, B is the magnetic field (where a factor of / √ π has been absorbed in its definition), E t is the total energydensity, e the specific internal energy and where we assumed anideal equation of state (EOS), that is γ = 5 / .Our computational domain is a rectangular box (400 kpc ×
200 kpc ×
200 kpc with × × cells, respectively),where x ∈ [-200,200] kpc, y ∈ [-100,100] kpc, and z ∈ [-100,100]kpc (see Fig. 4). The right-hand half of the domain is filled with aturbulent medium (see Sec. 2.1), representing a realistic ICM, whilethe left-hand half contains a uniform medium in which the shockis launched. We define a shock discontinuity at x = − kpc(see Fig. 4 for the initial configuration of the magnetic field). Thisdefines three regions in our simulation box: a post-shock uniformregion (I), a pre-shock uniform region (II) and a pre-shock turbulentregion (III).The turbulent medium is produced externally, with the pro-cedure outlined in Sec. 2.1. The turbulent fields are then readinto PLUTO and interpolated onto the numerical mesh used toevolve shocks with a bi- or tri-linear interpolation at the desiredcoordinate location at the initial time. The initial boundary condi-tions of the computational domain are outflow in x (zero gradientacross the boundary) and periodic in y and z . We used a piece-wise parabolic method (PPM) for the spatial integration, whereasa nd oder TVD Runge-Kutta method for the time stepping with a Courant-Friedichs-Lewy (CFL) condition of . . The Riemannsolver for the flux computation that we used is the Harten-Lax-vanLeer-Discontinuities (HLLD) solver (see Miyoshi & Kusano 2005).We control the ∇ · B = 0 condition with the hyperbolic divergencecleaning technique where the induction equation is coupled to ageneralized Lagrange multiplier (GLM) (e.g. Dedner et al. 2002).In Appendix A, we compare how the GLM and constrained trans-port (CT) methods control the divergence-free condition and findthat both methods work reliably in PLUTO.The initial conditions for the density, pressure and velocity in re-gion II ( pre-shock uniform region at [-100,0] kpc) are set to themean value of the corresponding turbulent fields. In the case ofthe magnetic field in region II, we set it to be the mean valueof the B x component of the turbulent medium. The initial condi-tions for region I ( post-shock region) are selected according to theMHD Rankine-Hugoniot conditions (e.g. Landau & Lifshitz 1987).We performed simulations with two different turbulent media (seeSec. 2.1) produced with the code FLASH (e.g. Fryxell et al. 2000;Calder et al. 2002) by varying also the strength of the shock andthe angle, θ bn , of the upstream magnetic field with respect to thenormal of the shock. Shocks can be classified as quasi-parallel andquasi-perpendicular if θ bn (cid:54) ◦ or θ bn > ◦ , respectively. Weconsider two limits, i.e. θ bn = 0 ◦ and θ bn = 90 ◦ . This sums upto a total of six runs for which all the parameters are summarizedin Table 1. We show the projection maps of the two different initialturbulent media considered for this work in Fig. 5 and clarify howthese were selected in Sec. 2.1.Finally, we fill the computational domain from the shock disconti-nuity up to the right side of the box with one Lagrangian particleper cell. This gives us a total number of 3,145,728 Lagrangian par-ticles for each run.
Each Lagrangian particle represents an ensemble of CR electronsand is characterized by an energy distribution function, χ ( E, τ ) = N ( E, τ ) n , (8)which gives the number of electrons per fluid number density.These particles evolve according to the cosmic-ray transport equa-tion, dχdτ + ∂∂E (cid:20)(cid:18) − E ∇ µ u µ − c r E (cid:19) χ (cid:21) = 0 , (9)where the first term in square brackets accounts for energy lossesdue to adiabatic expansion, while the second one accounts for thesynchrotron and inverse-Compton losses for electrons with isotrop-ically distributed velocity vectors, c r = 43 σ T cβ m e c (cid:20) B π + a rad T (1 + z ) (cid:21) , (10)where β = v e /c is the velocity of the electrons, m e their massand a rad the radiation constant. For the present work, we assumea constant redshift of z = 0 . The reader may refer to Vaidya et al.(2018) for a complete description of the numerical implementationand the semi-analytical scheme used for solving Eq. (9). Note that here we define the direction of the upstream magnetic field asthe direction of the mean magnetic field of the turbulent mediumc (cid:13) , 000–000 orphology of radio relics Run ID Turbulent medium M θ bn [ ◦ ] ρ II [10 − g / cm ] B x, II [ µG ] k1p5 M2 parallel L/ L/ L/ L/ L/ × − k4 M3 perpendicular L/ × − Table 1.
Initial conditions: We denote our regions I, II, III where I is the post-shock region ([-2,-1] in box coordinates), II is the uniform pre-shock region ([-1,0] in box coordinates) and III is the turbulent pre-shock region ([0,2]). The initial conditions for the left side of the shock (region I) depend on the pre-shockconditions (region II) and the initial Mach number of the shock M through the Rankine-Hugoniot jump conditions. L denotes the length of the turbulentregion, i.e. 200 kpc (see more details in Sec. 2.1). Note that the magnetic field in region II has only an x-component. y [ k p c ] t=0 Gyr 2L/3 t=0 Gyr L/4 y [ k p c ] x [100 kpc] y [ k p c ] x [100 kpc] D e n s i t y [ − m u ] V e l o c i t y [ k m / s ] T e m p e r a t u r e [ K / µ ] D e n s i t y [ − m u ] V e l o c i t y [ k m / s ] T e m p e r a t u r e [ K / µ ] Figure 5.
Projected maps along the z -axis of the different initial conditions(25 kpc slice), where m u is the atomic mass unit and µ is the mean molecu-lar weight. Here we show the right-hand half-size of the box containing theturbulent media. From top to bottom the reader can see the average density,velocity and temperature field for the L/ ( left column ) and L/ ( rightcolumn ) cases, where L/ and L/ refer to the integral scale of each typeof turbulence ( L = 200 kpc). We study a simplified scenario where the non-thermal spectraof the particles are activated instantly at the shock discontinuity.While the particles follow the fluid flow since t = 0 , the particle’senergy distribution starts to evolve only when the particles havepassed a shock. We implemented a shock finder based on converg-ing flows and pressure jumps that we describe in Appendix B.Once the Lagrangian particles are activated at the shock dis-continuity, they get assigned an initial energy distribution. The cor-responding particle momentum distribution function is a power-lawdistribution, i.e f ( p ) ∝ p − q . The power-law index is given by the diffusive shock acceleration (DSA) theory (Drury 1983). q = 3 rr − , (11)where r is the shock compression ratio defined as the ratio of theupstream and downstream densities. If we consider test particle ac-celeration at a shock of Mach number M , then it is possible tore-write Eq. (11) making use of the Rankine-Hugoniot jump equa-tions (see derivation in Blandford & Eichler 1987): q = 3( γ + 1) M M −
1) = 4 M M − , (12)where for the second equality we have considered an adiabatic in-dex of γ = 5 / . The corresponding power-law index for themacro-particle distribution function can be also obtained by assum-ing isotropy, i.e N ( p, τ ) = (cid:82) Ω τ p f ≈ πp f , where d Ω τ is thesolid angle around the direction τ . Since the particles are relativis-tic, we have that N ( E, τ ) dE = N ( p, τ ) dp . Moreover, from theisotropic condition we also have that N ( E ) = 4 πN ( E, τ ) .Therefore, the tracer particle energy distribution function atthe activation time is given by χ ( E ) = N ( E ) n = N n E − p , (13)where p = q − is what it is usually called the injection spectralindex , N is the normalisation factor and n is the fluid numberdensity. N is assigned according to the energy contained in theshock. That is, we considered that the total energy per fluid numberdensity is (cid:90) χ ( E ) E dE = E tot n , (14)where E tot = η E shock = η ρ post v and η is the accelera-tion efficiency. For which finally one can obtain the normalisationfactor: N = η E shock (4 − q ) [ E max4 − q − E min4 − q ] if q (cid:54) = 4 η E shock log (cid:16) E max E min (cid:17) if q = 4 (15)The PLUTO code allows us to compute the maximum and mini-mum energy at each time step. The maximum energy E max is im-posed considering the maximum allowed Larmor radius for eachparticle. The minimum energy E min is estimated balancing thesynchrotron and acceleration time scales (see full explanation inVaidya et al. 2018; B¨ottcher & Dermer 2010; Mimica & Aloy 2012)so it depends on the acceleration efficiency. Moreover, the accel-eration efficiency in collisionless relativistic shocks can also de-pend on the energy of the particles (see Sironi et al. 2013). In this c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al. work, for simplicity we assume a fixed acceleration efficiency of η = 10 − and fixed energy limits of γ min = 1 and γ max = 10 .This acceleration efficiency agrees with the expectations of DSAfor strong shocks (e.g. Kang et al. 2012b) and lies in the rangeof values required to explain observations of radio relics (see Bot-teon et al. 2020). Nevertheless, this means that the final synchrotronemission obtained can be re-scaled to the desired η since in our casethe energy limits remain constant. The synchrotron emission of a tracer particle in a local magneticfield B (cid:48) in the direction ˆn (cid:48) los , per unit solid angle, volume and fre-quency is given by J (cid:48) syn ( ν (cid:48) obs , ˆn (cid:48) los , B (cid:48) ) = (cid:90) N ( E (cid:48) ) P ( ν (cid:48) obs , E (cid:48) , φ (cid:48) ) dE (cid:48) d Ω (cid:48) , (16)where P ( ν (cid:48) obs , E (cid:48) , φ (cid:48) ) is the synchrotron power per unit frequencyand unit solid angle emitted by a single particle that has energy E (cid:48) and φ (cid:48) is the angle that the velocity vector of the particle makeswith the direction ˆn (cid:48) los . Following Ginzburg & Syrovatskii 1965we get J (cid:48) syn ( ν (cid:48) obs , ˆn (cid:48) los , B (cid:48) ) = √ e πm e c | B (cid:48) × ˆn (cid:48) los | (cid:90) N ( E (cid:48) ) F ( ξ ) dE (cid:48) , (17)where ˆn (cid:48) los is the unit vector in the direction of the line of sight inthe comoving frame and F ( ξ ) is a Bessel function integral givenby F ( ξ ) = ξ (cid:90) ∞ ξ K / ( z (cid:48) ) dz (cid:48) , (18)where ξ = ν (cid:48) obs ν (cid:48) c = 4 πm e c ν (cid:48) obs eE (cid:48) | B (cid:48) × ˆn (cid:48) los | , (19)where ν (cid:48) c is defined as the critical frequency at which the emissionpeaks. Note that only those particles with a pitch angle coincid-ing with the angle between B (cid:48) and ˆn (cid:48) los contribute to the emissionalong the line of sight in Eq. (17).The emissivity in Eq. (17) is measured in the local comovingframe with the emitting volume. If we want the emissivity in a fixedobserver’s frame then we have to apply a transformation : J syn ( ν obs , ˆn los , B ) = D J (cid:48) syn ( ν (cid:48) obs , ˆn (cid:48) los , B (cid:48) ) , (20)where D is a Doppler factor given by D = 1 γ (1 − ˆn los · v ) , (21)where γ is the Lorentz factor of the tracer particle. The unit vectorson the direction of the line of sight in the comoving and observer’sframe are related through ˆn (cid:48) los = D (cid:20) ˆn los + (cid:18) γ γ + 1 v · ˆn los − γ (cid:19) v (cid:21) . (22) The reader should note that the primed quantities in Eq. (20) refer to thecomoving frame, whereas standard notation refers to the observer’s frame.
In the preceding equations, the vector ˆn los can be selected accord-ing to an observing angle θ obs with respect to the line-of-sight. Inthis work, we only show results considering θ obs = 0 ◦ , that is weconsider an observer’s reference frame in which z lies along theline of sight ˆn los and x and y are in the plane of the sky. The spe-cific intensity (or surface brightness) maps can then be obtained byintegrating along a line of sight as I ν = (cid:90) J syn ( ν obs , x, y, z ) dz, (23)in units of [erg cm − s − Hz − str − ]. This is doable due to thefact that the emissivity information J syn of the Lagrangian tracerparticles is interpolated back onto the Eulerian grid. It is also pos-sible to obtain spectral index maps by means of − α ( x, y ) = log [ I ν ( x, y ) /I ν ( x, y )]log( ν /ν ) , (24)and the integrated spectra (or net flux) can be obtained by integrat-ing the specific intensity I ν over the area covered by the source inthe plane of the sky, that is F ν ( ν ) = (cid:90) I ν ( ν, x, y ) dxdy, (25)in units of [erg s − Hz − str − ]. In this subsection we describe some of the features of the fluid start-ing with the evolution of velocity, temperature and magnetic field,for the L/ , M = 3 and θ bn = 90 ◦ case (see Fig. 6). In all ofour runs, the initial shock sweeps region II maintaining a constantvelocity and planar shape due to the fact that this region is initiallyuniform. Next, the shock enters region III, where the shock will nolonger be uniform and be affected by the anisotropies of the fluid.The two first rows in Fig. 6 highlight how weak internal shocksand turbulence are generated in the downstream as the main shocktravels through the simulation box. The third row shows the stream-lines of the magnetic field. The streamlines are coloured accordingto the magnitude of the y -component of the magnetic field, B y , toillustrate how the shock compression amplifies the field. The fourthrow shows 1D profiles of the magnetic field obtained by integrat-ing along the LOS. As expected from the MHD Rankine-Hugoniotconditions, the component parallel to the shock normal, B x , is con-served, while the other components are amplified and stretched as aresult of shock compression. We show also the magnetic field pro-file weighted with the synchrotron emission at 150 MHz (blue solidlines).In Fig. 7 we show how the standard deviation (volumetricvalue) of the magnetic and velocity field evolve for all runs listedin Tab. 1. The evolution is characterized by two phases: 1) a firstphase in which the shock is crossing region II (lasting 50–150 Myrdepending on the type of run), and we have a purely decaying tur-bulence system on the right-hand side of the box ([0,200] kpc, seeFig. 4), and 2) a second phase in which the shock has alreadyentered region III and is compressing the turbulent medium. Thedashed gray vertical lines in Fig. 7 define the beginning of this sec-ond phase (the time differs according to the initial Mach number ofeach run). c (cid:13) , 000–000 orphology of radio relics Figure 6.
Evolution of the L/ , M = 3 and θ bn = 90 ◦ case. First row : 20 kpc slice along the z-axis of velocity field.
Second row : 20 kpc slice alongthe z-axis of temperature
Third row : streamlines of total magnetic field (note that the colour-code only denotes the strength of y -component). Fourth row :1D magnetic profiles corresponding. The component’s profiles have no weighting, whereas the magnetic field strength, B , is weighted with the 150 MHzsynchrotron emissivity.Run ID σ V [km / s] σ B [ µ G] t shock [Myr]k1p5 M2 parallel 217.7 1.016 282k1p5 M3 parallel 388.5 1.016 188k4 M2 parallel 133.1 0.659 439k4 M3 parallel 246.1 0.513 292k1p5 M3 perpendicular 388.5 1.099 188k4 M3 perpendicular 246.1 0.513 293 Table 2.
Values used in Fig. 7 and 8. The second and third column are theinitial standard deviation of the velocity and magnetic field in the wholesimulation box. The fourth column shows the total time needed for theshock to cross the entire simulation box.
The evolution of the standard deviation of the velocity andmagnetic field are in general not correlated (see third column ofFig. 7). We also analysed the standard deviation of the parallel andperpendicular components of both fields with respect to the shocknormal. In general, we find that the evolution of the parallel com-ponent dominates the evolution of σ V , whereas the perpendicularcomponent dominates the evolution of σ B . The standard deviationof the parallel component (see first column of Fig. 7) decreaseswith time for both fields. The standard deviations of the perpen-dicular components (see second column of Fig. 7) show a differentevolution: the perpendicular component, σ V ⊥ , follows the same be- haviour of σ V , while the perpendicular σ B always increases rightafter the shock crossing in all runs. This trend is persistent and getsstronger with higher resolution (see Fig. C1 in Appendix C).The reason for this is that most velocity fluctuations are drivenin the direction parallel to the shock normal, while magnetic fieldfluctuations are initially driven perpendicular to the shock normaldue to compression. With time, the dynamics gets more compli-cated due to the shock compression and possibly also stretching ofthe magnetic field. The increase in σ B induces a delay in the ve-locity field dissipation. The plateau observed in σ V in the secondphase (evolution to the right of the dashed gray lines in Fig. 7) in-dicates that the shock-induced turbulence could be maintained onlyfor some time before σ V decreases due to turbulent dissipation.For all our runs, we verified that the density and also the tem-perature evolve in the same fashion as the velocity field in Fig. 7and it is only the standard deviation of the magnetic field that has acharacteristic evolution in both phases.In the following, we summarize the observed effects for thissecond phase:i) The role of different turbulent injection scales : We find that σ B decreases faster in a system with smaller magnetic fluctuations( L/ ) than in one with larger fluctuations ( L/ ). This is expectedas the turbulence injection scale is smaller in the L/ case andtherefore the eddy-turnover time is shorter. In fact, in the L/ case, c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al.
Figure 7.
Evolution of the standard deviation of the magnetic field ( left axis ) and velocity field ( right axis ). Top panels : L/ case. Bottom panels : L/ case. First column : Parallel component to the shock, i.e B (cid:107) = B x and v (cid:107) = v x . Second column : Perpendicular component, i.e B ⊥ = (cid:113) B y + B z and v ⊥ = (cid:113) v y + v z . Third column : Standard deviation of the magnetic and velocity field strengths. σ B and σ V are normalised to their value at t = 0 forpurposes of comparison between different runs (see Table 2). the shock-induced turbulence seems to only have an increasing ef-fect on σ B whenever the shock is perpendicular, i.e. θ bn = 90 ◦ .ii) The role of the main shock’s Mach number : The largest im-pact on increasing σ B is due to a larger Mach number. For example,the M = 3 and θ bn = 0 ◦ case shows an increase of ∼ % for the L/ turbulence and ∼ % for the L/ turbulence. Conversely, the M = 2 and θ bn = 0 ◦ case shows an increase of ∼ % for the L/ turbulence and a decrease of ∼ % for the L/ turbulence.This suggests that weak shocks ( M =
2) are less likely to mod-ify the initial distribution of magnetic fields for smaller turbulentscales.iii)
The role of obliquity : A perpendicular shock has thestrongest effect on broadening the downstream magnetic field dis-tribution. This is expected as the perpendicular components of themagnetic field are the only ones affected by the shock compres-sion. The largest increase in both runs is ∼ % (with respect tothe corresponding dashed gray vertical line in Fig. 7).Finally, we show in Fig. 8 the profiles of the total magnetic fieldstrength for all runs for two snapshots. These 1D profiles are ob-tained integrating along the LOS and then taking an average. Eachprofile is again normalised to the pre-shock magnetic field of eachrun, B pre . The two times shown in Fig. 8 correspond to when theshock is starting to compress the turbulent medium (left panel ofFig. 8), and when the shock has crossed the whole simulation box(right panel of Fig. 8). The overall magnetic amplification in thedownstream region differs in each run. The discrepancies observedin Fig. 8 can be explained by the different distribution of incidence angles between the upstream magnetic field and the shock normalalong the shock front in the different runs. The downstream regionin the M = 3 case develops very similar profiles for both typesof turbulence. In particular, the L/ case leads to higher magneticamplification (compare blue and red lines), whereas the L/ caseleads to comparable downstream magnetic profiles (compare pur-ple and orange lines). We observe that a lower Mach number, suchas M = 2 , leads to less magnetic amplification, owing to the lowershock compression factor. In this case, the final extent of the regionwhere a magnetic field amplification occurs in the downstream is ∼ kpc larger than in the case with a M = 3 shock. This meansthat the strength of this shock is insufficient to further compress theturbulent medium, thus producing a more extended turbulent mag-netic region which will be also reflected in the synchrotron emis-sion. We present a few three-dimensional renderings of the synchrotronemission produced by our modelling in Fig. 9, as seen along differ-ent lines of sight. Although the radio emission seems fairly uniformwhen observed edge-on (see following Fig. 10), the emission is notspatially uniform, but concentrated into threads and filaments inthe shock plane. The combination of shock compression and turbu-lence introduces anisotropies and fluctuations in the flow. This, inturn, directly affects the advection properties and energy evolutionof the CR particles.We show the surface brightness maps at 150 MHz for all runs c (cid:13)000
The role of obliquity : A perpendicular shock has thestrongest effect on broadening the downstream magnetic field dis-tribution. This is expected as the perpendicular components of themagnetic field are the only ones affected by the shock compres-sion. The largest increase in both runs is ∼ % (with respect tothe corresponding dashed gray vertical line in Fig. 7).Finally, we show in Fig. 8 the profiles of the total magnetic fieldstrength for all runs for two snapshots. These 1D profiles are ob-tained integrating along the LOS and then taking an average. Eachprofile is again normalised to the pre-shock magnetic field of eachrun, B pre . The two times shown in Fig. 8 correspond to when theshock is starting to compress the turbulent medium (left panel ofFig. 8), and when the shock has crossed the whole simulation box(right panel of Fig. 8). The overall magnetic amplification in thedownstream region differs in each run. The discrepancies observedin Fig. 8 can be explained by the different distribution of incidence angles between the upstream magnetic field and the shock normalalong the shock front in the different runs. The downstream regionin the M = 3 case develops very similar profiles for both typesof turbulence. In particular, the L/ case leads to higher magneticamplification (compare blue and red lines), whereas the L/ caseleads to comparable downstream magnetic profiles (compare pur-ple and orange lines). We observe that a lower Mach number, suchas M = 2 , leads to less magnetic amplification, owing to the lowershock compression factor. In this case, the final extent of the regionwhere a magnetic field amplification occurs in the downstream is ∼ kpc larger than in the case with a M = 3 shock. This meansthat the strength of this shock is insufficient to further compress theturbulent medium, thus producing a more extended turbulent mag-netic region which will be also reflected in the synchrotron emis-sion. We present a few three-dimensional renderings of the synchrotronemission produced by our modelling in Fig. 9, as seen along differ-ent lines of sight. Although the radio emission seems fairly uniformwhen observed edge-on (see following Fig. 10), the emission is notspatially uniform, but concentrated into threads and filaments inthe shock plane. The combination of shock compression and turbu-lence introduces anisotropies and fluctuations in the flow. This, inturn, directly affects the advection properties and energy evolutionof the CR particles.We show the surface brightness maps at 150 MHz for all runs c (cid:13)000 , 000–000 orphology of radio relics Figure 8.
1D magnetic field profile of all runs.
Left panel : Shock front justentering the turbulent medium at t shock / . Right panel : Shock front is al-most at the right end of the x -axis at t shock / . Refer to Tab. 2 for thevalue of t shock of each run. at a time when the shock front has reached almost the right end ofthe simulation box in Fig. 10 (the different times are specified in theupper left corner of each panel) . Note that we only applied Gaus-sian smoothing for the surface brightness maps in Fig. 10 (meant tomimic the finite spatial resolution of a typical LOFAR-HBA obser-vation), while all of our following analysis was done without anysmoothing. In the 3D view presented in Fig. 9, we can distinguishthe complex substructure of the emission in the form of filaments,bristles, ribbons or other structures that cannot be classified in a sin-gle group. Nevertheless, since the emissivity is not volume-fillingand the strength of the emission varies from region to region, someof this structure vanishes in projection.The morphology observed in Figs. 9 and 10 depends on threefactors: the strength of the shock, the obliquity of the shock, andthe type of turbulence. A higher Mach number leads to strongeremission and elongated patterns due to a stronger compression ofthe magnetic field. On the other hand, the role of the upstream tur-bulence is more complex.One can observe the effect of more elongated patterns for thecases with a shock of M = 3 (comparing panels (a), (c), (d) and(f) of Fig. 10) due to the increased stretching of the eddies.On the other hand, a shock of strength M = 2 can producedisrupted patterns in our simulated emission. The shock front is lesslikely to look totally disrupted in the L/ case (see panel (e)), whilethis effect is particularly noticeable in the L/ case (see panel (b)).We notice that this is a direct consequence of the sound speed in themedia. For example, the L/ turbulence easily leads to regionswith Mach numbers lower than our threshold ( M min = L/ turbulence is slightly higher than that of the L/ case (seeSec. 2.1). Subsequently, turbulent dissipation leads to an increasein sound speed and therefore, lowering the Mach number. Such lowMach numbers are not expected to accelerate electrons via the DSAprocess and therefore, they are excluded from our modelling. Infact, Kang et al. 2019 found that quasi-perpendicular shocks with We have included the surface brightness maps along the x-axis in Ap-pendix D for completeness.
Figure 9.
Visualization of synchrotron emissivity J syn isocurves for the L/ , M = 3 and θ bn = 0 ◦ run at t = 178 Myr. The emissivity is shownin units of [erg cm − s − Hz − str − ]. The axis are shown in units of [100kpc]. M (cid:46) . may not efficiently accelerate electrons through DSA.Assessing the impact of different turbulent injection scales in M =2 shocks, would then require a tailored set-up which we leave forfuture work.Finally, obliquity produces more elongated emission becauseshock compression only amplifies the component parallel to theshock front (i.e. the B y and B z components in this study).In Fig. 11, we show 1D profiles of the emission presented inFig. 10, at two more observing frequencies: 1.5 GHz and 650 MHz. c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al. x [100 kpc] y [ k p c ] t=178 Myr (a) 2L/3, M =3 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] x [100 kpc] y [ k p c ] t=268 Myr (b) 2L/3, M =2 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] x [100 kpc] y [ k p c ] t=178 Myr (c) 2L/3, M =3 , θ bn =90 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] x [100 kpc] y [ k p c ] t=292 Myr (d) L/4, M =3 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] x [100 kpc] y [ k p c ] t=439 Myr (e) L/4, M =2 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] x [100 kpc] y [ k p c ] t=293 Myr (f) L/4, M =3 , θ bn =90 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] Figure 10.
Surface brightness at 150 MHz for all runs in Tab. 1 (see Equation 23). We considered a beam of θ = 15” × to get the surface brightness( θ I ν ) in units of µ Jy/beam. We smoothed the maps with a Gaussian kernel with
FWHM = 7 . kpc (assuming z = 0 . ). The most extended emission is found in the runs with M = 2 (seebottom panel of Fig. 11) in agreement with the 1D magnetic pro-files previously shown in Fig. 8. These low Mach number runs alsoshow the lowest values of surface brightness along the downstreamdue to the initial normalisation depending on the Mach number.This suggests that, in our survey of merger shocks, only those with M ∼ would require a higher acceleration efficiency to reach ob-servable values as it has been pointed out in previous works (e.g.Botteon et al. 2020). Finally, the steepness of the emission profiledepends on the magnetic morphology. Comparing all the M = 3 runs, the L/ turbulence case shows a shallower decline com-pared to the L/ case. This is a direct consequence of the initialdistribution of the magnetic field strength. In Fig. 3 of Sec. 2.1, weshowed that initially the L/ run reaches higher magnetic field val-ues in the tail of its PDF, which leads to larger synchrotron losses. In Fig. 12, we show the spectral index maps for the two differentturbulent media with the M = 3 shock using Eq. (24). We canobserve the expected spectral index gradient starting from the shockfront (red) to the end of the downstream region (blue). In the L/ case, the spectral index values range between α = − . and α = − . , while it goes from α = − . to α = − . in the L/ case. This agrees with the previous subsection, where we discussedthe emission profiles. In the same way, the turbulent medium withsmaller initial fluctuations L/ is more likely to produce a steepergradient because the initial magnetic field distribution has a largertail (see Fig. 3).In Fig. 13, we show the corresponding spectral index profilesalong the downstream region for completeness. In the top panelof Fig. 13 we show how the profile changes when taking into ac-count different frequencies for one specific run: L/ , M = 3 , θ bn = 0 ◦ . The profiles start to differ beyond ∼ kpc fromthe shock front where the emission at lower frequencies decreasesmore slowly. In the lower panel, we show the spectral index profilesfor all of our runs between 650 and 150 MHz. The selected snap- c (cid:13) , 000–000 orphology of radio relics Figure 11.
1D surface brightness profiles obtained from the emission mapsat 1.5 GHz ( top panel ), 650 MHz ( middle panel ) and 150 MHz ( bottompanel ) for all runs. shots correspond to those in Fig. 10. For the higher Mach number(i.e. M = 3 ), the L/ profiles agree more with observations thanthe L/ profiles. In the L/ case, the spectral index profiles aresteeper than observed. For example, the relic in the cluster MACSJ0717.5+3745 (see van Weeren et al. 2017; Bonafede et al. 2018)with a Mach number of M = 2 . (inferred from the injectionspectral index α ), shows a spectral index steepening up to valuesof ∼ − . over a region of (cid:46) kpc. Another example is the“Toothbrush” relic (see Rajpurohit et al. 2018, 2020) which steep-ens also up to values of ∼ − . over a region of ∼ kpc. Thissuggests that the initial turbulent magnetic field distribution beforea shock crossing is rather narrow. For the forced turbulence used inthis work, this means that the standard deviation must be smallerthan σ B ∼ µ G (see Fig. 3 in Sec. 2.1 for the initial magnetic fielddistribution). On the other hand, this also suggests that the injectionscale (and also the magnetic coherence scale) of the turbulence ingalaxy clusters outskirts could be L/ ( ∼
133 kpc) or even larger.
Figure 12.
Spectral index maps obtained between 150 MHz and 650 MHzat t shock using Equation (24). Top panel : L/ , M = 3 , θ bn = 0 ◦ run. Bottom panel : L/ , M = 3 , θ bn = 0 ◦ . It is fairly common in observations of radio relics, that only theintegrated spectral index can be computed. This is done by fittingthe total observed flux (see Eq. 25) against different available fre-quencies. We computed the integrated spectral index in this fashionusing 1.5 GHz, 650 MHz and 150 MHz frequencies. In Fig. 14 weshow how the integrated spectral index evolves as the shock sweepsacross the simulation box. In the first ∼
60 Myr the integrated spec-tral indices differ quite significantly, whereas after ∼
140 Myr thevalue of the integrated spectra seems to converge to the same valuefor all runs. The relation between the real spectral index α and theintegrated one is often assumed to be (e.g. Kardashev 1962; Heav-ens & Meisenheimer 1987): α int = α + 12 = M + 1 M − , (26)and therefore we also plot the expected α int for different Machnumbers as a reference in Fig. 14 with gray dashed horizontal lines.The integrated spectral index from our runs does not follow a char- c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al. acteristic pattern and it does not remain strictly constant through atime span of roughly 200–300 Myr. In addition, we include the cor-responding evolution of the integrated spectral index for two extraruns with a completely uniform medium (density, velocity, pressureand magnetic field) with shocks of strength M = 2 and M = 3 .In the uniform media, Eq. (26) holds after ∼ Myr when the en-ergy spectrum reaches a steady state at the shock (see Kang et al.2017 for a one-dimensional uniform media study), however this isnot the case for all the turbulent media. In particular, the effect ofturbulence on the 3D distribution of the synchrotron emissivity isthat it makes it patchy and not volume-filling (as can be observedin Fig. 9). This in turn has an effect on the integrated flux and there-fore, the integrated spectral index.Hence, it can be difficult to recover the real spectral index andMach number through this method, even if recent high-resolutionobservations of radio relics show consistent integrated and injectionspectra (e.g. Rajpurohit et al. 2020). The fact that Eq. (26) onlyholds for a planar shock where all the fields are uniform, suggeststhat one should be careful when making use of it. One good way toconfirm if it is applicable or not is by cross-checking with the resultfrom high-resolution spectral index maps.This has also been found in previous studies. For example,Kang 2015 showed that the relation in Eq. (26) only holds for pla-nar shocks, but not for spherical shocks. In the presence of a turbu-lent medium, the geometry of the shock is more complicated lead-ing to the evolution observed in Fig. 14.We find that the integrated spectral index in Fig. 14 is biasedtowards higher Mach numbers for all runs. The main reason forthis is that the brightest radio emitting regions correspond to thestrongest shock compression regions, which in average gives abias towards higher Mach numbers. In Sec. 4.4. we will show ananalysis of the Mach number distribution and compare it to the oneinferred from the thermal fluid.In Fig. 15 we show how the emission at the shock front corre-lates with the magnetic field strength in 3D and 2D. We show theevolution for the L/ case with M = 3 and θ bn = 0 ◦ run. Theemissivity computed from Eq. (20) scales with the magnetic fieldand frequency as J syn ∝ B ( p +1) / ν − ( p − / , (27)(see Engel 1979) where α = ( p − / is the spectral index and p is related to the Mach number through Eqs. (13) and (12), α = ( p − q − M + 32( M − . (28)In the top panel of Fig. 15, we show this relation for different Machnumbers with coloured lines and we add an additional black dashedline corresponding to a fit of all the data points. Overall, we find thatthere is always a systematic mismatch with the initial real Machnumber of the shock. During the first ∼ Myr, the 3D distribu-tion shows a sharper relation coinciding with what is expected fromEq. (28). Nevertheless, shortly after that there will be a spread in theMach number and magnetic field distribution along the shock front.This spread will keep changing as a consequence of all the turbu-lent motions in the medium. The black dashed line shows that therelation is biased towards larger Mach numbers. For the 2D caseshown in the bottom panel of Fig. 15, we considered a magneticfield weighted with the radio emissivity: B w = (cid:82) B J syn dZ (cid:82) J syn dZ . (29) Figure 13.
Top panel : Spectral index profile for the L/ , M = 3 , θ bn =0 ◦ case at t = 178 Myr.
Bottom panel : Spectral index profiles between 150MHz and 650 MHz for all the cases. The profiles are computed at the sametimes as in Fig. 10.
We show the same relation pointed out in the upper panel as a ref-erence. The dashed black line in this case corresponds to the fit ofall data points in the 2D map. While in this case we have less datapoints due to the integral along the LOS, it is interesting to see thatthe bias towards higher Mach numbers is still there. This suggeststhat the bias is not due to projection effects. We will further discussthe reason behind this bias in Sec. 4.4. c (cid:13)000
We show the same relation pointed out in the upper panel as a ref-erence. The dashed black line in this case corresponds to the fit ofall data points in the 2D map. While in this case we have less datapoints due to the integral along the LOS, it is interesting to see thatthe bias towards higher Mach numbers is still there. This suggeststhat the bias is not due to projection effects. We will further discussthe reason behind this bias in Sec. 4.4. c (cid:13)000 , 000–000 orphology of radio relics Figure 14.
Integrated spectral index evolution computed from fitting thetotal flux. The gray horizontal dashed lines show the expected integratedspectral index assuming that α int = α + 1 / for different Mach numbers(see Eq. 26). The shaded areas correspond to the uniform runs with M = 2 ( blue ) and M = 3 ( gray ) shocks. One important feature that characterizes the magnetic field is its power spectrum : P ij ( k ) = 1(2 π ) (cid:90) (cid:90) (cid:90) e − i k · x R ij ( k ) d k , (30)where R ij = (cid:104) b i ( x ) b j ( x + x ) (cid:105) is the two-point correlationfunction between the magnetic fields b i and b j (e.g. Batchelor1951). In the case of homogeneous and isotropic fields the relationbetween the spectral energy and the 1D power spectrum is found tobe E ( k ) = 2 πk P ii ( k ) . (31)We will refer to the 1D power spectrum P ii ( k ) simply as P ( k ) . Inorder to obtain the 1D power spectrum we averaged the 3D powerspectrum of the magnetic field over spherical shells: P ( k ) = 1 N k (cid:88) k − < | k | (cid:54) k + P ( k ) . (32)We computed the energy power spectrum in the right half part ofthe simulation (region III, [0,200] kpc) (see Fig. 4). In the top panelof Fig. 16 we show the final magnetic energy spectra for all our runsand in the bottom panel we show the final power spectra computedfor the emissivity at 150 MHz. In Fig. 17, we present the wholeevolution of the magnetic energy spectrum E B ( k ) condensed in theform maps. Each of these maps contains the following information:the y -axis shows the evolution and the x -axis shows the wavenum-ber k coloured with the amplitude of the magnetic energy spectrumat that k . In this way, the darker colours denote the regions wherethe power is peaking (see Fig. 16 as a reference). The dashed gray Figure 15.
Phase-plots of the magnetic field versus the radio emission at150 MHz at the cells where the shock front is located for the L/ , M = 3 , θ bn = 0 ◦ run. Top panel : values extracted out of the 3D distributions, i.e.the emissivity J syn . The coloured lines show the expected fit for differentMach numbers and the black dashed line shows a fit of the data at t = 188 Myr.
Bottom panel : values extracted out of the 2D maps, i.e here we havethe surface brightness and the values of the magnetic field weighted withthe emission (see Eq. 29). line is a reference for the reader to know when the shock enters thisturbulent region, i.e. region III. We are interested in understandingthe effect of shocks in the magnetic energy spectrum, in presenceof an ICM with decaying turbulence. We can see in Fig. 17 twoimportant features: i) the wavenumbers k (cid:38) (corresponding toscales (cid:46) kpc) are largely unaffected by shocks with strength M = θ bn = 90 ◦ cases (see the two bottom panels of c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al.
Figure 16.
Top panel : Final magnetic energy spectrum for all our runs. Thespectra are computed in a (200 kpc) volume through Eqs. (31) and (32).The final time step differs for each run and it is specified in the legend. Bottom panel : Final power spectrum of the synchrotron emissivity at 150MHz for all our runs (same times as in the top panel).
Fig. 17) at k (cid:46) (corresponding to scales (cid:38) kpc). In thesecases, power shifts towards smaller wavenumbers (larger scales)due to the new turbulence introduced after the shock passage. Thisagrees with previous results from cosmological MHD simulations.In Dom´ınguez-Fern´andez et al. 2019, we observed the same effectwhen analysing the evolution of a merging galaxy cluster over atime span of almost 10 Gyr. After every merger, shock waves arecreated and after every shock crossing, the magnetic power shiftsto smaller wavenumbers (or larger scales). In this subsection we make the same analysis as in subsection 4.3.1,but for the radio emission and we compare it to the results from themagnetic field. We compute the power spectrum of the 3D distri-butions and also of the 2D maps (computed by integrating alongthe LOS). The 1D power spectrum is obtained by averaging the 3Dspectrum over spherical shells (see Eq. (32)) as mentioned above;while for the 2D maps, we averaged the 2D spectrum over annuli.Afterwards we can compute the characteristic length of the powerspectrum for the emission as λ c = (cid:82) k − P ( k ) dk (cid:82) P ( k ) dk , (33)and for the magnetic field as, λ B = (cid:82) k − P B ( k ) dk (cid:82) P B ( k ) dk , (34)where P ( k ) and P B ( k ) correspond to the power spectrum of thesynchrotron emission (see Eq. (20)) and the power spectrum of themagnetic field, respectively.In the first two columns of Fig. 18 we show the results for the3D case for all of the runs including also the characteristic scale ofthe magnetic field λ B . The characteristic scale of the radio emis-sion is in general of the same order of the characteristic scale ofthe magnetic field, that is of the order of (cid:46)
100 kpc. There are somespecifics regarding 1) the Mach number: higher Mach numbers leadto larger emission scales, for example a M = 3 leads to a maxi-mum scale of ∼ kpc, whereas a M = 2 leads to a maximumof ∼ kpc; and 2) the injection scale of the turbulence in the θ bn = 0 ◦ case: the characteristic scale of the emission seems notto be affected by the injection scale, but it cannot be directly cor-related with the magnetic field’s scale. For example, the top panelsof the first two columns in Fig. 18 show how λ c is of the same or-der for both cases while the underlying turbulence is different. Incontrast, λ B is directly affected by the underlying turbulence andtherefore, it is of a different order. The type of injection scale of theturbulence in the θ bn = 90 ◦ case plays a big role for the evolutionof λ B , but this is not reflected in the evolution of λ c . This hap-pens only due to the fact that the acceleration efficiency, η , doesnot depend on θ bn in our modelling. The role of a varying η will besubject of future work.The characteristic scales of the integrated LOS variables are alsoshown in Fig. 18 (last two columns). In this case, λ B can be smallerthan in reality is by roughly 17–23% in some cases, while λ c is onlysmaller by ∼ While the shock front’s Mach number distribution peaks at its ini-tial Mach number for almost the whole evolution, it develops a tail c (cid:13)000
100 kpc. There are somespecifics regarding 1) the Mach number: higher Mach numbers leadto larger emission scales, for example a M = 3 leads to a maxi-mum scale of ∼ kpc, whereas a M = 2 leads to a maximumof ∼ kpc; and 2) the injection scale of the turbulence in the θ bn = 0 ◦ case: the characteristic scale of the emission seems notto be affected by the injection scale, but it cannot be directly cor-related with the magnetic field’s scale. For example, the top panelsof the first two columns in Fig. 18 show how λ c is of the same or-der for both cases while the underlying turbulence is different. Incontrast, λ B is directly affected by the underlying turbulence andtherefore, it is of a different order. The type of injection scale of theturbulence in the θ bn = 90 ◦ case plays a big role for the evolutionof λ B , but this is not reflected in the evolution of λ c . This hap-pens only due to the fact that the acceleration efficiency, η , doesnot depend on θ bn in our modelling. The role of a varying η will besubject of future work.The characteristic scales of the integrated LOS variables are alsoshown in Fig. 18 (last two columns). In this case, λ B can be smallerthan in reality is by roughly 17–23% in some cases, while λ c is onlysmaller by ∼ While the shock front’s Mach number distribution peaks at its ini-tial Mach number for almost the whole evolution, it develops a tail c (cid:13)000 , 000–000 orphology of radio relics Figure 17.
Evolution of the magnetic energy spectrum in the turbulent part of the box (region III [0,2] in Fig. 4 corresponding to a (200 kpc) cube). At eachtime-step ( y -axis), we show the magnetic power spectrum by colouring the x -axis with its amplitude (see colourbar in logarithmic scale).c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al.
80 100 120 140 160 180
Time [Myr] λ c [ k p c ] M =3 , θ bn =0 ◦ B
150 200 250 300
Time [Myr] λ c [ k p c ] L/4, M =3 , θ bn =0 ◦ B
150 200 250
Time [Myr] λ c [ k p c ] M =2 , θ bn =0 ◦ B
200 250 300 350 400 450
Time [Myr] λ c [ k p c ] L/4, M =2 , θ bn =0 ◦ B
80 100 120 140 160 180
Time [Myr] λ c [ k p c ] M =3 , θ bn =90 ◦ B
150 200 250 300
Time [Myr] λ c [ k p c ] L/4, M =3 , θ bn =90 ◦ B Figure 18.
First two columns : Characteristic scale of the 3D distributions according to Eq. (33).
Last two columns : Characteristic scale of the 2D distributions.Each panel shows the evolution of λ c for 1.5 GHz, 650 MHz and 150 MHz and also the evolution of the magnetic field’s characteristic scale λ B . towards higher Mach numbers owing to turbulent motions . In or-der to study its impact on the emission, we obtain the emissivity, J syn , corresponding to the cells tracking the shock front at eachtime step and correlated it with its Mach number. We can do thisbecause the particles are being activated whenever they detect theshock front. In Fig. 19, we show the binned statistics for the twoquantities obtained considering 16 time steps.In order to generate an X-ray alike estimate of the Mach num-ber, we computed the X-ray emissivity J X − ray = 1 . × − T / (cid:18) ρm p (cid:19) , (35)in units of [erg cm − s − Hz − str − ], where m p is the proton’smass and T and ρ are in cgs units. Thus, we selected the emission J syn and J X − ray only at the shocked cells at each time step inorder to compute the binned statistics. We show the distributionsat different times with dashed coloured lines and that of the wholeevolution with the solid blue lines in Fig. 19. The top panels showthe binned statistics considering the radio emissivity and the bottompanels show that of the X-ray emissivity.The discrepancy between the two statistical distributions is ev-ident. While the radio emissivity is always biased towards largerMach numbers, the X-ray emissivity is biased towards smaller see an example for the case L/ , M = 3 , θ bn = 0 ◦ in Fig. B2 of theAppendix B Mach numbers. In fact, the peak of the binned distribution of theradio emissivity (temporal envelope) is similar to the values in thetail of the real 3D Mach number distribution (i.e. see an examplein Fig. B2 in Appendix B). Out of the six cases analysed, only the L/ , M = 2 , θ bn = 0 ◦ case shows a partial match between bothdistributions. We also observe that the radio-Mach number statis-tics fluctuates with time such that the peak Mach number can varya fair bit.In addition, we analysed two runs with uniform media (allfields) but the same Mach numbers (i.e. M = 2 and M = 3 ).We verified that in this case the distributions of both emissivitiesand Mach numbers are the same. Hence, our results suggest thatthe difference in Mach numbers in radio and in X-ray is a result ofturbulence.In the presence of turbulence, the radio emissivity is notvolume-filling, while the X-ray emissivity comes from the wholeshock front. The fact that compression is different from region toregion leads to a patchier radio emissivity that can probe a limitedpart of the shock front (see Fig. 9). Apart from that, the radio emis-sion is biased towards higher Mach numbers because the initial CRenergy is ∝ M (which in our work corresponds to the normal-isation N defined in Eq.(15)). Finally, the magnetic field fluctu-ations created after the shock-crossing play an important role. Asdiscussed in Sec. 4.1, the amplitude of the magnetic field fluctu-ations (which directly affects the radio emission) decreases moreslowly than the velocity, density and temperature fluctuations. This c (cid:13)000
Last two columns : Characteristic scale of the 2D distributions.Each panel shows the evolution of λ c for 1.5 GHz, 650 MHz and 150 MHz and also the evolution of the magnetic field’s characteristic scale λ B . towards higher Mach numbers owing to turbulent motions . In or-der to study its impact on the emission, we obtain the emissivity, J syn , corresponding to the cells tracking the shock front at eachtime step and correlated it with its Mach number. We can do thisbecause the particles are being activated whenever they detect theshock front. In Fig. 19, we show the binned statistics for the twoquantities obtained considering 16 time steps.In order to generate an X-ray alike estimate of the Mach num-ber, we computed the X-ray emissivity J X − ray = 1 . × − T / (cid:18) ρm p (cid:19) , (35)in units of [erg cm − s − Hz − str − ], where m p is the proton’smass and T and ρ are in cgs units. Thus, we selected the emission J syn and J X − ray only at the shocked cells at each time step inorder to compute the binned statistics. We show the distributionsat different times with dashed coloured lines and that of the wholeevolution with the solid blue lines in Fig. 19. The top panels showthe binned statistics considering the radio emissivity and the bottompanels show that of the X-ray emissivity.The discrepancy between the two statistical distributions is ev-ident. While the radio emissivity is always biased towards largerMach numbers, the X-ray emissivity is biased towards smaller see an example for the case L/ , M = 3 , θ bn = 0 ◦ in Fig. B2 of theAppendix B Mach numbers. In fact, the peak of the binned distribution of theradio emissivity (temporal envelope) is similar to the values in thetail of the real 3D Mach number distribution (i.e. see an examplein Fig. B2 in Appendix B). Out of the six cases analysed, only the L/ , M = 2 , θ bn = 0 ◦ case shows a partial match between bothdistributions. We also observe that the radio-Mach number statis-tics fluctuates with time such that the peak Mach number can varya fair bit.In addition, we analysed two runs with uniform media (allfields) but the same Mach numbers (i.e. M = 2 and M = 3 ).We verified that in this case the distributions of both emissivitiesand Mach numbers are the same. Hence, our results suggest thatthe difference in Mach numbers in radio and in X-ray is a result ofturbulence.In the presence of turbulence, the radio emissivity is notvolume-filling, while the X-ray emissivity comes from the wholeshock front. The fact that compression is different from region toregion leads to a patchier radio emissivity that can probe a limitedpart of the shock front (see Fig. 9). Apart from that, the radio emis-sion is biased towards higher Mach numbers because the initial CRenergy is ∝ M (which in our work corresponds to the normal-isation N defined in Eq.(15)). Finally, the magnetic field fluctu-ations created after the shock-crossing play an important role. Asdiscussed in Sec. 4.1, the amplitude of the magnetic field fluctu-ations (which directly affects the radio emission) decreases moreslowly than the velocity, density and temperature fluctuations. This c (cid:13)000 , 000–000 orphology of radio relics Figure 19.
Binned statistics of radio emissivity at 1.5 GHz ( left axis ) and X-ray emissivity ( right axis ) with the Mach number of the shock front for dif-ferent time steps (dashed coloured lines). The dark blue line shows the statistics taken into account the whole evolution and the shadowed areas denote thecorresponding 75-th percentile of the distribution. The emissivity is shown in units of [ − erg cm − s − Hz − str − ].c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al. adds to the discrepancy because the X-ray emissivity depends onlyon the temperature and density fields.
We have presented a hybrid framework to compute the synchrotronemission from a shock wave propagating through a medium withdecaying turbulence representing a small fraction of the ICM. Inour framework, the MHD grid represents a thermal fluid, whereasLagrangian particles represent CR electrons. We injected CR elec-trons at the shock discontinuity assuming diffusive shock accelera-tion. Each CR electron evolves according to the cosmic-ray trans-port equation in the diffusion approximation.Our simulations explored shocks with Mach numbers charac-teristic of radio relics, i.e. M = 2 and M = 3 . Moreover, wevaried the downstream turbulence using turbulence-in-a-box simu-lations: a solenoidal subsonic turbulence with power peaking at 2/3of the box (case L/ ) and a solenoidal subsonic turbulence withpower peaking at 1/4 of the box (case L/ ). One snapshot of eachsimulation was selected as an initial condition for our shock-tubesimulation. Our results can be summarized as follows:i) Impact of a shock on decaying turbulence : We find thatmild shocks produce magnetic fluctuations in the downstreamregion that do not correlate with fluctuations in velocity, densityand temperature. In fact, we find that magnetic fluctuations canincrease even when velocity fluctuations decrease. This, in turn,can affect the final extent of the magnetic downstream. We findthe strongest effect in perpendicular shocks, as expected fromtheory. Shocks with M = 2 travelling in a medium with smallerfluctuations, such as our L/ case, cause the least effect and seemto hardly modify the initial magnetic field distribution.ii) Radio emission : The existence of substructure in the syn-chrotron emission is a direct consequence of a turbulent medium.We found that M = 2 shocks in our set-up are unlikely toreproduce observable radio relics. The physical reason behind thisis that M = 2 shocks are not strong enough to modify the initialpre-shock magnetic field. For example, the relic at Abell 2744(e.g. Govoni et al. 2001; Eckert et al. 2016; Pearce et al. 2017;Paul et al. 2020) reaches a surface brightness of the order of tensof mJy/beam at 1.4 GHz, which would require an accelerationefficiency of η ∼ in our L/ turbulence.iii) Discrepancies in the spectral index : Our spectral index pro-files suggest that in the case of M = 3 shocks, a turbulent injectionscale of L/ or even larger reproduces observations better thanthe L/ case. The L/ initial magnetic field distribution allowsfor higher values of the magnetic field strength reflected in the tailof the PDF which steepens the spectral index profiles more than inthe L/ case. We conclude that an initial turbulent, magnetic fielddistribution in the ICM must have a standard deviation smallerthan σ B = 1 µ G. We compare our results of the integrated spectralindex to the relation α int = α + 1 / and find that this relation doesnot seem to hold in the presence of a turbulent medium. The reasonfor this is the distribution of Mach numbers within a shock front ina turbulent medium. As a consequence, the injected electrons willhave different initial energy spectra.iv) Discrepancies in Mach numbers : We find that the synchrotronemission is biased towards larger Mach numbers when comparingto the X-ray emission. This agrees with previous numerical work (e.g. Hong et al. 2015) and a number of observations of radiorelics. For example, X-ray observations of the Toothbrush relic inthe cluster 1RXS J0603.3+4214 infer a Mach number of
M ∼ . (e.g. Ogrean et al. 2013; van Weeren et al. 2016), while radioobservations infer a higher Mach number of M ∼ . (e.g.Rajpurohit et al. 2018, 2020). The source of this discrepancy liesin 1) the stronger dependence of the synchrotron emission on thecompression in the shock and 2) the fact that the amplitude ofthe magnetic field fluctuations (which affect the radio emission)decreases more slowly than the density and temperature fluctua-tions. Hence, higher Mach numbers in the tail of the Mach numberdistribution bias the overall Mach number.v) Magnetic energy spectrum : We find that scales (cid:46) kpc arelargely unaffected by shocks with M = Characteristic length of the radio emission : The characteristiclengths derived from the power spectrum of the emission, λ c , andmagnetic field, λ B , are of the same order. We find that λ c is ingeneral of the order of (cid:46) kpc. Analysing the LOS variables,we do not observe strong projection effects and λ B and λ c are only17–23% and 10–15% smaller than in 3D.In summary, we could identify the most important featuresthat link the observable properties of radio relics with the dynam-ical properties of the upstream ICM. Our work confirms that theMach numbers inferred from the radio emission are likely to beoverestimates of the real Mach number of the thermal fluid in thepresence of turbulence. This has been previously pointed out as apossible solution that can alleviate the problem of acceleration ef-ficiencies and as a possible explanation for the non-detection of γ -rays from galaxy clusters (see Ackermann et al. 2010, 2014, 2016and Vazza & Br¨uggen 2014). While CRe and CRp are expected tobe accelerated at the same places, their acceleration mechanismsand efficiencies will differ (e.g. Caprioli & Spitkovsky 2014a).CRe are accelerated preferentially at quasi-perpendicular shocksand CRp at quasi-parallel shocks. Recent works by Wittor et al.2020 and Banfi et al. 2020 used cosmological MHD simulations toshow that indeed the predominance of quasi-perpendicular shocksin merger and accretion shocks might be enough to explain the ab-sence of CRp. Therefore, in future work we will include the obliq-uity dependence in our acceleration efficiency, η . Moreover, we in-tend to survey a larger range of possible parameters, both, for theICM conditions and for the shock properties. This will help us toassess whether the large variety of relic sources can be explained bythe model adopted here. We will also present a detailed study con-sidering polarisation (Q and U Stokes parameters, as well as Rota-tion Measure), different lines-of-sight, projection and beam effects. The analysis presented in this work made use of computationalresources on the JUWELS cluster at the Juelich SupercomputingCentre (JSC), under project ”stressicm” with F.V. as principal in-vestigator and P.D.F as co-principal investigator. The 3D visualiza-tion of the synchrotron emissivity was done with the VisIt software(see Childs et al. 2012).P.D.F, F.V. and M.B. acknowledge the financial support from theEuropean Union’s Horizon 2020 program under the ERC Starting c (cid:13)000
M ∼ . (e.g. Ogrean et al. 2013; van Weeren et al. 2016), while radioobservations infer a higher Mach number of M ∼ . (e.g.Rajpurohit et al. 2018, 2020). The source of this discrepancy liesin 1) the stronger dependence of the synchrotron emission on thecompression in the shock and 2) the fact that the amplitude ofthe magnetic field fluctuations (which affect the radio emission)decreases more slowly than the density and temperature fluctua-tions. Hence, higher Mach numbers in the tail of the Mach numberdistribution bias the overall Mach number.v) Magnetic energy spectrum : We find that scales (cid:46) kpc arelargely unaffected by shocks with M = Characteristic length of the radio emission : The characteristiclengths derived from the power spectrum of the emission, λ c , andmagnetic field, λ B , are of the same order. We find that λ c is ingeneral of the order of (cid:46) kpc. Analysing the LOS variables,we do not observe strong projection effects and λ B and λ c are only17–23% and 10–15% smaller than in 3D.In summary, we could identify the most important featuresthat link the observable properties of radio relics with the dynam-ical properties of the upstream ICM. Our work confirms that theMach numbers inferred from the radio emission are likely to beoverestimates of the real Mach number of the thermal fluid in thepresence of turbulence. This has been previously pointed out as apossible solution that can alleviate the problem of acceleration ef-ficiencies and as a possible explanation for the non-detection of γ -rays from galaxy clusters (see Ackermann et al. 2010, 2014, 2016and Vazza & Br¨uggen 2014). While CRe and CRp are expected tobe accelerated at the same places, their acceleration mechanismsand efficiencies will differ (e.g. Caprioli & Spitkovsky 2014a).CRe are accelerated preferentially at quasi-perpendicular shocksand CRp at quasi-parallel shocks. Recent works by Wittor et al.2020 and Banfi et al. 2020 used cosmological MHD simulations toshow that indeed the predominance of quasi-perpendicular shocksin merger and accretion shocks might be enough to explain the ab-sence of CRp. Therefore, in future work we will include the obliq-uity dependence in our acceleration efficiency, η . Moreover, we in-tend to survey a larger range of possible parameters, both, for theICM conditions and for the shock properties. This will help us toassess whether the large variety of relic sources can be explained bythe model adopted here. We will also present a detailed study con-sidering polarisation (Q and U Stokes parameters, as well as Rota-tion Measure), different lines-of-sight, projection and beam effects. The analysis presented in this work made use of computationalresources on the JUWELS cluster at the Juelich SupercomputingCentre (JSC), under project ”stressicm” with F.V. as principal in-vestigator and P.D.F as co-principal investigator. The 3D visualiza-tion of the synchrotron emissivity was done with the VisIt software(see Childs et al. 2012).P.D.F, F.V. and M.B. acknowledge the financial support from theEuropean Union’s Horizon 2020 program under the ERC Starting c (cid:13)000 , 000–000 orphology of radio relics Grant ”MAGCOW”, no. 714196. W.B.B. acknowledges the finan-cial support from the Deutsche Forschungsgemeinschaft (DFG) viagrant BR2026125.We acknowledge our anonymous reviewer for helpful commentson the first version of this manuscript. Finally, we thank D. Ryu, C.Federrath and R. Mohapatra for fruitful scientific discussions.
The data underlying this article will be shared on reasonable re-quest to the corresponding author.
REFERENCES
Ackermann M. et al., 2014, ApJ, 787, 18Ackermann M. et al., 2016, ApJ, 819, 149Ackermann M. et al., 2010, ApJL, 717, L71Banda-Barrag´an W. E., Federrath C., Crocker R. M., BicknellG. V., 2018, MNRAS, 473, 3454Banfi S., Vazza F., Wittor D., 2020, MNRAS, 496, 3648Batchelor G. K., 1951, Mathematical Proceedings of the Cam-bridge Philosophical Society, 47, 359374Birdsall C. K. N., Estacio E. T., Plasma Theory, Simulation Group(PTSG), 2004, Computer Physics Communications, 164, 189Blandford R., Eichler D., 1987, Physical Reports, 154, 1Blandford R. D., Ostriker J. P., 1978, ApJL, 221, L29Bonafede A. et al., 2018, MNRAS, 478, 2927B¨ottcher M., Dermer C. D., 2010, ApJ, 711, 445Botteon A., Brunetti G., Ryu D., Roh S., 2020, A & A, 634, A64Br¨uggen M., Bykov A., Ryu D., R¨ottgering H., 2012, Science &Space Review, 166, 187Brunetti G., Jones T. W., 2014, International Journal of ModernPhysics D, 23, 1430007Bykov A. M., Kaastra J. S., Br¨uggen M., Markevitch M., FalangaM., Paerels F. B. S., 2019, Science & Space Review, 215, 27Calder A. C. et al., 2002, ApJS, 143, 201Caprioli D., Spitkovsky A., 2014a, ApJ, 783, 91Caprioli D., Spitkovsky A., 2014b, ApJ, 794, 46Childs H. et al., 2012, in High Performance Visualization–Enabling Extreme-Scale Scientific Insight, pp. 357–372Clarke T. E., Ensslin T. A., 2006, ApJ, 131, 2900Dedner A., Kemm F., Kr¨oner D., Munz C. D., Schnitzer T., We-senberg M., 2002, Journal of Computational Physics, 175, 645Di Gennaro G. et al., 2018, ApJ, 865, 24Dom´ınguez-Fern´andez P., Vazza F., Br¨uggen M., Brunetti G.,2019, MNRAS, 486, 623Donnert J., Vazza F., Br¨uggen M., ZuHone J., 2018, ArXiv e-printsDrury L. O., 1983, Reports on Progress in Physics, 46, 973Eckert D., Jauzac M., Vazza F., Owers M., Kneib J.-P., TcherninC., Intema H., Knowles K., 2016, Mon. Not. Roy. Astron. Soc.,461, 1302Egan H., O’Shea B. W., Hallman E., Burns J., Xu H., Collins D.,Li H., Norman M. L., 2016, ArXiv e-printsEngel A. R., 1979, Physics Bulletin, 30, 158Eswaran V., Pope S. B., 1988, Physics of Fluids, 31, 506Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., MacLow M. M., 2010, A & A, 512, A81Fryxell B. et al., 2000, ApJS, 131, 273Ginzburg V. L., Syrovatskii S. I., 1965, ARAA, 3, 297 Govoni F., Feretti L., Giovannini G., B¨ohringer H., Reiprich T. H.,Murgia M., 2001, A & A, 376, 803Guo X., Sironi L., Narayan R., 2014, ApJ, 797, 47Heavens A. F., Meisenheimer K., 1987, MNRAS, 225, 335Hong S. E., Kang H., Ryu D., 2015, ApJ, 812, 49Kang H., 2015, Journal of Korean Astronomical Society, 48, 9Kang H., Ryu D., Ha J.-H., 2019, ApJ, 876, 79Kang H., Ryu D., Ha J.-H., 2019, The Astrophysical Journal, 876,79Kang H., Ryu D., Jones T. W., 2012a, ApJ, 756, 97Kang H., Ryu D., Jones T. W., 2012b, ApJ, 756, 97Kang H., Ryu D., Jones T. W., 2017, ApJ, 840, 42Kardashev N. S., 1962, Astronomicheskii Zhurnal, 39, 393Kritsuk A. G. et al., 2011a, ApJ, 737, 13Kritsuk A. G. et al., 2011b, ApJ, 737, 13Landau L. D., Lifshitz E. M., 1987, Fluid Mechanics, Second Edi-tion: Volume 6 (Course of Theoretical Physics), 2nd edn., Courseof theoretical physics / by L. D. Landau and E. M. Lifshitz, Vol.6. Butterworth-HeinemannLee D., Deane A. E., Federrath C., 2009, Astronomical Societyof the Pacific Conference Series, Vol. 406, A New Multidimen-sional Unsplit MHD Solver in FLASH3, p. 243Locatelli N. T. et al., 2020, MNRASMignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O.,Zanni C., Ferrari A., 2007, ApJS, 170, 228Mimica P., Aloy M. A., 2012, MNRAS, 421, 2635Miniati F., 2014, ApJ, 782, 21Miyoshi T., Kusano K., 2005, Journal of Computational Physics,208, 315Nuza S. E., Gelszinnis J., Hoeft M., Yepes G., 2017, MNRAS,470, 240Ogrean G. A., Brggen M., van Weeren R. J., Rttgering H., CrostonJ. H., Hoeft M., 2013, Monthly Notices of the Royal Astronom-ical Society, 433, 812824Owen F. N., Rudnick L., Eilek J., Rau U., Bhatnagar S., Kogan L.,2014, ApJ, 794, 24Park J., Caprioli D., Spitkovsky A., 2015, Physical Review Let-ters, 114, 085003Paul S., Salunkhe S., Sonkamble S., Gupta P., Mroczkowski T.,Raychaudhury S., 2020, Astronomy & Astrophysics, 633, A59Pearce C. J. J. et al., 2017, ApJ, 845, 81Pinzke A., Oh S. P., Pfrommer C., 2013, MNRAS, 435, 1061Rajpurohit K. et al., 2018, ApJ, 852, 65Rajpurohit K. et al., 2020, A & A, 636, A30Ryu D., Kang H., Cho J., Das S., 2008, Science, 320, 909Ryu D., Kang H., Ha J.-H., 2019, ApJ, 883, 60Ryu D., Kang H., Hallman E., Jones T. W., 2003, ApJ, 593, 599Schaal K., Springel V., 2015, MNRAS, 446, 3992Schmidt W., Niemeyer J. C., Hillebrandt W., 2006, A & A, 450,265Sironi L., Spitkovsky A., Arons J., 2013, ApJ, 771, 54Skillman S. W., Xu H., Hallman E. J., O’Shea B. W., Burns J. O.,Li H., Collins D. C., Norman M. L., 2013, ApJ, 765, 21Vaidya B., Mignone A., Bodo G., Rossi P., Massaglia S., 2018,ApJ, 865, 144van Weeren R. J., Br¨uggen M., R¨ottgering H. J. A., Hoeft M.,2011, Journal of Astrophysics and Astronomy, 32, 505van Weeren R. J. et al., 2016, The Astrophysical Journal, 818, 204van Weeren R. J., de Gasperin F., Akamatsu H., Br¨uggen M., Fer-etti L., Kang H., Stroe A., Zandanel F., 2019, Science & SpaceReview, 215, 16van Weeren R. J. et al., 2017, ApJ, 835, 197 c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al. van Weeren R. J. et al., 2017, The Astrophysical Journal, 835, 197van Weeren R. J., R¨ottgering H. J. A., Br¨uggen M., Hoeft M.,2010, Science, 330, 347van Weeren R. J., R¨ottgering H. J. A., Intema H. T., Rudnick L.,Br¨uggen M., Hoeft M., Oonk J. B. R., 2012, A & A, 546, A124Vazza F., Br¨uggen M., 2014, MNRAS, 437, 2291Vazza F., Brunetti G., Br¨uggen M., Bonafede A., 2018, MNRAS,474, 1672Vazza F., Dolag K., Ryu D., Brunetti G., Gheller C., Kang H.,Pfrommer C., 2011, MNRAS, 418, 960Vazza F., Jones T. W., Br¨uggen M., Brunetti G., Gheller C., PorterD., Ryu D., 2017, MNRAS, 464, 210Wittor D., Hoeft M., Vazza F., Br¨uggen M., Dom´ınguez-Fern´andez P., 2019, MNRAS, 490, 3987Wittor D., Vazza F., Ryu D., Kang H., 2020, MNRAS, 495, L112
APPENDIX A: NUMERICAL DIVERGENCE OF THEMAGNETIC FIELD
The ∇· B = 0 condition in the particle module of the PLUTO codeis maintained with the hyperbolic divergence cleaning techniquewhere the induction equation is coupled to a generalized Lagrangemultiplier (GLM) (e.g. Dedner et al. 2002). In Kritsuk et al. 2011a,it has been argued that the best results for the divergence-free evolu-tion of the magnetic field are achieved using the constrained trans-port (CT) method. Keeping this in mind, we tested the differencesof the ∇ · B = 0 condition with the GLM and CT techniques. InFig. A1, we show the evolution of the magnetic field divergenceusing the same set-up as our L/ , M = 3 and θ bn = 90 ◦ run(see Table 1) for the whole simulation box. Both methods are com-parable and keep the numerical magnetic field divergence under ∼ . % of the local magnetic field value. This demonstrates thatthe use of the GLM divergence cleaning technique is robust againstCT for our particular set-up. Thus, we do not expect our numericalscheme to have an impact in our final synchrotron emission maps.Additionally, we show how the divergence condition behaves forthe interface between regions II and III in our setup in Fig. A1. Theeffects of the initial interpolation of the external input get dimin-ished after a few steps and the numerical magnetic field divergencedrops again below ∼ . % before the shock enters the region ofinterest (i.e. region III); see an additional discussion about this inAppendix D of Banda-Barrag´an et al. 2018. APPENDIX B: SHOCK FINDER
The algorithm to find shocks is already implemented in the PLUTOcode (see Vaidya et al. 2018). In the first step, shock cells are taggedthrough the ∇ · v < condition. Then we implemented an extracondition for the activation time of the tracer particles in order tocompute the initial Mach number and compression ratio. This inturn is needed for assigning the initial energy spectra with an spec-tral index q (see Eq. (12)) to each tracer particle. The Mach numberat the shock center is computed from the Rankine-Hugoniot pres-sure jump condition: ∆ log p (cid:62) log p p (cid:12)(cid:12)(cid:12)(cid:12) M = M min , (B1)where the subscripts 1 and 2 denote the pre-shock and post-shock regions, respectively. The minimum Mach number is set to M min = 1 . and it acts as a threshold to filter out weaker shocks. Time [Myr] -11 -10 -9 -8 -7 -6 -5 h | ∇ · B / B | GLMCTInterface
Figure A1.
Test on the numerical conservation of the ∇ · B = 0 condition(where h = ∆ x ). We compare two runs only differing in their divergencecleaning method: GLM ( blue ) and CT ( purple ). The gray line shows thecorresponding evolution for the interface between regions II and III usingthe GLM method. We performed this test using the set-up L/ , M = 3 and θ bn = 90 ◦ . Figure B1.
Numerical test on the shock’s directionality. We show a 2Dparallel shock propagating with an angle of ◦ . The pressure jump is computed with the neighbouring cells alongthe three directions for which a Mach number distribution is finallyobtained: M = M x + M y + M z (see Ryu et al. 2003, Vazzaet al. 2011, Schaal & Springel 2015).The PLUTO code is able to compute the compression ratio for theupdate of the spectra once the particle has crossed a shock, never-theless we implemented this extra condition as it was necessary forcomputing the compression ratio at the time of activation and nota time step afterwards as in Vaidya et al. 2018. In this fashion, thetracer particles have an initial spectrum consistent with DSA the- c (cid:13)000
Numerical test on the shock’s directionality. We show a 2Dparallel shock propagating with an angle of ◦ . The pressure jump is computed with the neighbouring cells alongthe three directions for which a Mach number distribution is finallyobtained: M = M x + M y + M z (see Ryu et al. 2003, Vazzaet al. 2011, Schaal & Springel 2015).The PLUTO code is able to compute the compression ratio for theupdate of the spectra once the particle has crossed a shock, never-theless we implemented this extra condition as it was necessary forcomputing the compression ratio at the time of activation and nota time step afterwards as in Vaidya et al. 2018. In this fashion, thetracer particles have an initial spectrum consistent with DSA the- c (cid:13)000 , 000–000 orphology of radio relics Figure B2.
Evolution of Mach number distribution at the shock front forthe case L/ , M = 3 , θ bn = 0 ◦ . ory at the time of activation and after a time step their spectra willevolve subject to radiative losses.We tested the directionality of the shock finder by setting up a 2Dparallel shock propagating with an angle of ◦ with respect to the x -axis and placing one particle per cell on a quarter of the domain.The shock propagates with a Mach number of M = 10 and theinterpolation of the grid quantity at the particle position in PLUTOis done with standard techniques used for Particle-In-Cell (PIC)codes (see Birdsall et al. 2004). We used the Nearest Grid Point(NGP) method for the implementation of this test. In Fig. B1, weshow the Mach grid distribution as well as the particle’s interpo-lated Mach number for a snapshot of this set-up. We also show theevolution of the shock’s Mach number in our final set-up for one ofour studied cases in Fig. B2. APPENDIX C: RESOLUTION
We tested our set-up doubling the resolution to × × cells and number of Lagrangian particles to 25,165,824. Here wewill limit ourselves to show the run L/ , M = 3 and θ bn = 90 ◦ .In Fig. C1 we show the evolution of the standard deviation of themagnetic and velocity field for the resolution used in this work, i.e. × × cells, and for the higher resolution. This resultwas shown and explained in detail in Sec. 4.1. This figure confirmsthat the velocity field (and also the density and temperature field)dynamics is largely unaffected by a higher resolution. In fact, it isonly the magnetic field fluctuations that are ∼ % enhanced. Thismeans that our result on the non-correlation between magnetic andvelocity fluctuations is even more accentuated at higher resolution.This is expected as we are increasing the effective Reynolds num-ber of the simulation. An upper limit of the Reynolds number isgiven by: Re max ≈ (cid:18) l ∆ x (cid:19) / , (C1) Figure C1.
Evolution of the standard deviation of the magnetic field ( leftgreen axis ) and velocity field ( right blue axis ). The solid lines correspondto the lower resolution run ( × × cells) and the dashed linescorrespond to the higher resolution run ( × × cells). where l is the maximum correlation scale in the flow and ∆ x is theresolution. In this case, l = 2 L/ ≈ . kpc, ∆ x = 1 . kpc and ∆ x = 0 . kpc would lead to an upper limit of theeffective Reynolds number of ∼ and ∼ for the low andhigh resolution runs, respectively. A lower limit is given in contrastby: Re min ≈ (cid:18) lε ∆ x (cid:19) / , (C2)where ε is a factor depending on the diffusivity of the numericalmethod. For second order finite difference/volume codes such asour case with the PLUTO code, one can assume ε ≈ (e.g. Kritsuket al. 2011b). This leads to a lower limit of the effective Reynoldsnumber of ∼ and ∼ for the low and high resolution runs,respectively.Finally, we show in Fig. C2 a comparison between surfacebrightness maps at 150 MHz with both resolutions. As it can be ob-served, the broader features as well as the extent of the downstreamare consistent in both resolutions, while the higher resolution runhighlights smaller features ( ∼ kpc). APPENDIX D: SURFACE BRIGHTNESS ALONG THEX-AXIS
We present the surface brightness maps at 150 MHz for all our runsas viewed from the x -axis in Fig. D1. We want to point out thatthese maps were obtained by projecting the emissivity already usedfor this work along the z -axis. In order to compute self-consistentlythe surface brightness maps changing the LOS, one has to changethe observing angle, θ obs , (see Sec. 3.3) and the vector ˆn los . Thisis turn shall be used for computing the integral and Doppler factorin the emissivity equation (see Sec. 3.2). The maximum value ofsurface brightness in each panel of Fig. D1, is lower than its corre-spondent panel in Fig. 10. This is because in this case, we end up c (cid:13) , 000–000 P. Dom´ınguez-Fern´andez et al.
Figure C2.
Surface brightness at 150 MHz for the run L/ , M = 3 and θ bn = 90 ◦ . We considered a beam of θ = 15” × to get the surface brightness( θ I ν ) in units of µ Jy/beam. No smoothing is considered here.
Left panel : × × cells simulation. Right panel : × × cells simulation. y [100 kpc] z [ k p c ] t=178 Myr (a) 2L/3, M =3 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] y [100 kpc] z [ k p c ] t=268 Myr (b) 2L/3, M =2 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] y [100 kpc] z [ k p c ] t=178 Myr (c) 2L/3, M =3 , θ bn =90 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] y [100 kpc] z [ k p c ] t=292 Myr (d) L/4, M =3 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] y [100 kpc] z [ k p c ] t=439 Myr (e) L/4, M =2 , θ bn =0 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] x [100 kpc] y [ k p c ] t=293 Myr (f) L/4, M =3 , θ bn =90 ◦ S u r f a c e b r i g h t n e ss a t M H z [ µ J y / b e a m ] Figure D1.
Surface brightness at 150 MHz integrated along the x -axis (in correspondence to Fig. 10). We considered a beam of θ = 15” × to get thesurface brightness ( θ I ν ) in units of µ Jy/beam. We smoothed the maps with a Gaussian kernel with
FWHM = 7 . kpc (assuming z = 0 . ). summing up the contribution of the emissivity for a smaller region( ∼ kpc). c (cid:13)000