Quasi-perpendicular shocks of galaxy clusters in hybrid kinetic simulations: The structure of the shocks

Context: Cosmic ray acceleration in galaxy clusters is still an ongoing puzzle, with relativistic electrons forming radio relics at merger shocks and emitting synchrotron radiation. In the present work we perform hybrid-kinetic simulations in a range of various quasi-perpendicular foreshock conditions, including plasma beta, magnetic obliquity, and the shock Mach number. Aims: We study the ion kinetic physics, responsible for the shock structure and wave turbulence, that in turn affects the particle acceleration processes. Methods: We apply a recently developed generalized fluid-particle hybrid numerical code that can combine fluid modeling for both kinetic ions and fluid electrons. The model utilizes the exact form of the generalized Ohm's law, allowing for arbitrary choice of mass and energy densities, as well as the charge-to-mass ratio of the kinetic species. Results: We show that the properties of ion-driven multi-scale magnetic turbulence in merger shocks are in agreement with the ion structures observed in PIC simulations. In typical shocks with the sonic Mach number $M_s=3$, the magnetic structures and shock front density ripples grow and saturate at wavelengths reaching approximately four ion Larmor radii. Finally, we note that the steady-state structure of $M_s=3$ shocks in high-beta plasmas shows evidence that there is little difference between 2D and 3D simulations. The turbulence near the shock front seems to be a 2D-like structure in 3D simulations.


Introduction
Galaxy clusters are vast collections of galaxies, hot gas, and dark matter bound together with gravitational forces.These objects often host numerous shock waves, including large-scale cluster merger shocks, generated during the collision and merging of galaxy clusters, when the gas clouds within them collide, creating shock waves that propagate through the gas and heat it to high temperatures.Observations of radio relics, which also emit X-rays, have provided compelling evidence for the acceleration of relativistic electrons at merger shocks (Brunetti & Jones 2014).Such radio relics in galaxy clusters are up-to Mpc-size regions of radio emission, typically with an arc-like morphology.This emission shows strong polarization and steep integrated ra-β ≫ 1, representing the ratio of thermal to magnetic pressures.Due to high temperatures, even the most energetic merger shocks typically have low sonic Mach numbers, usually M s ≲ 5.
Diffusive Shock Acceleration (DSA), also known as the firstorder Fermi process, is widely considered the most plausible mechanism responsible for particle acceleration to high energies (see, e.g., Bell 1978a,b;Blandford & Ostriker 1978;Drury 1983).In this process, particles gain energy by repeatedly crossing the shock front while being scattered between upstream and downstream regions.The essential element of DSA is the particle injection mechanism.DSA is effective only for particles with Larmor radii larger than the shock thickness, typically a few gyroradii of thermal ions.Therefore, thermal particles should be pre-accelerated to suprathermal momenta, an order of a few ion thermal momentum.The injection is more difficult for electrons than ions due to the lower mass and smaller Larmor radii of electrons.Thus, electron injection to DSA requires some interactions with the waves in the shock that provide their pre-acceleration.For the conditions of low-Mach-number shocks in hot ICM, such mechanisms remain not fully understood.
It has been found in numerous simulations that protons are typically more efficiently accelerated at quasi-parallel shocks, with θ Bn < 45 • , (e.g., Kang & Ryu 2013;Caprioli & Spitkovsky 2014;Ryu et al. 2019), while acceleration of electrons is more efficient at quasi-perpendicular shocks, θ Bn > 45 • , (Guo et al. 2014a,b;Kang et al. 2019;Ha et al. 2021;Amano & Hoshino 2022), where θ Bn is the angle between the shock normal and the large-scale upstream magnetic field.The mechanisms responsible for electron injection at quasi-perpendicular ICM shocks have been extensively debated in recent years.Shock Drift Acceleration (SDA) has gained attention as a plausible mechanism for electron (pre-)acceleration in these conditions, as suggested by numerous PIC simulations (e.g., Krauss-Varban & Wu 1989;Matsukiyo et al. 2011;Guo et al. 2014a,b).These simulations revealed that at subluminal shocks, the SDA-accelerated electrons can be reflected off the shock to propagate along the upstream mean magnetic field and excite magnetic turbulence in the foreshock region (Matsukiyo et al. 2011;Guo et al. 2014a,b;Kang et al. 2019).This turbulence is driven by electron firehose instability (EFI) and can scatter electrons back toward the shock, enabling repeating cycles of interactions with the shock front.It has been shown, that this multi-cycle SDA process leads to the generation of supra-thermal electron spectra upstream of the shock and is more efficient in shocks with higher plasma beta.Nevertheless, the maximum electron energy that can be achieved in this process is still below the estimated injection threshold for DSA.Recent large-scale PIC simulations of quasiperpendicular subluminal shocks (Kobzar et al. 2021;Ha et al. 2021) shown that stochastic Shock Drift Acceleration (SSDA, Katou & Amano 2019) is the primary mechanism capable of providing electron pre-acceleration up to the injection energy.During SSDA, electrons undergo stochastic pitch-angle scattering off multi-scale magnetic waves.This allows them to remain confined within the shock transition for an extended time, gaining energy through gradient drift.A crucial condition for SSDA is the presence of wideband turbulence, with wavelengths extending to the ion-scale long-wave shock rippling.Such turbulence can be driven by ion kinetic instabilities resulting from temperature anisotropies at the shock.One such instability is the Alfv'en ion cyclotron (AIC) instability, frequently associated with the generation of shock rippling (Matsukiyo & Matsumoto 2015).
Fully kinetic PIC simulations are often constrained by limitations in both space and time, restricting their applicability to large-scale systems due to significant computational demands.
On the contrary, hybrid kinetic simulations, which combine kinetic ions and fluid electrons, can handle considerably larger macroscopic systems, which makes them particularly valuable for advancing our understanding of the ion DSA process and the ion-scale turbulence in non-relativistic shocks.A combined hybrid kinetic and test-particle model has been recently developed by Trotta & Burgess (2019) to study 2D and 3D turbulent structures at low-Mach-number shocks in low-beta solar-wind plasmas and the electron acceleration therein.This approach allowed the authors to study the influence of the shock surface fluctuations on the acceleration of suprathermal electrons.
The aim of this work is to investigate the ion-scale structure of quasi-perpendicular merger shocks with the hybrid kinetic very large-scale 2D numerical experiments, much larger than investigated previously with fully kinetic PIC simulations by Guo et al. (2018); Kobzar et al. (2021); Ha et al. (2021).The simulations explore a range of shock conditions, including high plasma beta, sonic Mach numbers, and subluminal shock obliquity angles.They are designed to enable us to observe the system's development of multi-scale ion-driven turbulence.The simulations are augmented with smaller-scale 3D runs.This analysis lays the foundation for further exploration of the role of ion-scale turbulence in electron acceleration at high-beta shocks and provides the basis for large-scale 3D simulations to study shock physics.
The paper is organized as follows.Section 2 introduces the numerical setup and methods used in the simulations.Section 3 describes the 2D structure and evolution of a generic high-beta shock with the sonic Mach number M s = 3 and β = 20.We analyze the linear and nonlinear properties of the wave turbulence.Section 4 systematically explores a dependence of the shock physics on the shock and ambient plasma parameters, such as plasma beta, the sonic Mach number, and shock obliquity angle.Section 5 compares results from 2D and 3D simulations.Finally, Section 6 summarizes our findings and presents concluding remarks.

Numerical Setup
Our studies use a new generalized fluid-particle hybrid numerical code of Amano (2018) that treats the collisionless plasma under the assumption of quasi-neutrality.The model generally consists of fluid ions, electrons, and an arbitrary number of kinetic species.The dynamics of fluid and kinetic populations are coupled in a self-consistent manner.The code employs the exact form of the generalized Ohm's law, which allows us to consider arbitrary mass and energy densities, as well as the charge-tomass ratio of the kinetic species.Here, we restrict this general model to only fluid electrons and kinetic ions.Our approach thus resembles a standard hybrid model but with proper treatment of finite electron inertia effects.We assume the ion-to-electron mass ratio m i /m e = 100, where m i and m e are the ion and electron mass, respectively.We demonstrate in Appendix A that our results are independent of this specific choice of the mass ratio.
This work shows the results of large-scale 2D simulations supported with two smaller-scale 3D runs.The 2D simulation setup utilizes a computational grid in the x − y plane.The simulation code is fully 3D, and the 2D simulations are performed with two grid points in the z-direction.The simulations follow all three components of particle momenta and electromagnetic fields.An electron-ion plasma beam is injected at the right side of the computational box and flows with a bulk simulation-frame velocity u 0 in the −x-direction.After reflecting off the conductive wall at the left boundary, the beam interacts with the incoming plasma, forming a shock that travels in the positive x-direction at the speed u sh , as measured in the upstream rest frame.The right x-boundary of the box is open, and periodic boundary conditions are applied in transverse y (and z) directions.
The injected plasma carries a large-scale magnetic field, B 0 , which lies in the x − y plane (i.e., in-plane in 2D simulations) and is inclined at the angle θ Bn with respect to the shock normal.We investigate quasi-perpendicular shocks, for which θ Bn > 45 • .In conjunction with the magnetic field, a motional electric field, E 0 = −[(u 0 /c) × B 0 ], is also initialized and points in the zdirection, E 0 = E 0z ẑ.The results of previous PIC simulations drive the choice of the in-plane magnetic field configuration for 2D runs Kobzar et al. (2021), which showed its fidelity in reproducing 3D physics, which is also supported by our 3D simulations, see Section 5.
In the following, the velocities are normalized to the ion Alfvén speed, u A = B 0 / √ 4πN 0i m i , calculated using the upstream magnetic field and the ion density, N 0i .The ratio between Alfvén speed and the speed of light is set to u A /c = 10 −4 , the adiabatic index is Γ = 5/3, and the resistivity is η = 0.The electrons and ions are initially in thermal equilibrium and their temperatures are T i = T e = T 0 .The simulation-frame Mach number is defined as: where c s is the sound velocity and k B the Boltzmann constant.
The shock Mach number in the shock rest frame is: and the density jumps across the shock, r(M s ), in the limit of weakly magnetized flows, is equal to: (3) The simulation parameters were selected to reproduce typical physical conditions observed in ICM shocks.We simulate shocks of sonic Mach number M s of total plasma beta: where N e = N i is the electron number density.The plasma beta is initially equally carried by electrons and ions, β e = β i .If T e = T i , then the plasma beta can be expressed as: where M A = u sh /u A is the Alfvénic Mach number of the shock.
Setting up the shock in the simulations with the given physical parameters of M s and β involves the change of the normalized ion thermal velocity, u th,i /u A = β i /2, which sets up also the normalized ion sound velocity, c s /u A = 2Γβ i /2, and the change of the normalized inflow velocity, i.e., the Alfvénic Mach number of the upstream flow, M A,0 = u 0 /u A .
Our unit of space is the ion skin depth: or the ion gyroradius where e is the electric charge.The time unit is the ion gyrofrequency, Ω ci = eB 0 /m i c.The time step δt = 10 −3 Ω −1 ci .We resolve λ i with 4 cells and use N ppc = N 0i = 128 particles (ions) per cell.These numerical parameters were selected based on convergence tests (see Appendix B).The transverse size of the simulation box is several r g to resolve large-scale shock ripple modes.In the figures presented below, the spatial coordinates are normalized with the ion gyroradius r g (M s , β) of each model, unless specified otherwise.
Investigations of merger shock structures in similar system parameter ranges have been performed earlier with 2D PIC simulations by Guo et al. (2018); Kobzar et al. (2021); Ha et al. (2021).However, due to computational constraints these simulations used moderately large transverse sizes of numerical boxes, typically L y ≈ 21.5λ i ≈ (1.7 − 5.5)r g in Guo et al. (2018), L y = 32λ i ≈ 7.8r g in Kobzar et al. (2021), and L y = (44 − 62)λ i ≈ (3.4 − 5.1)r g in Ha et al. (2021), and simulation times typically t max ∼ (20 − 32)Ω −1 ci , with only two long-term runs with t max Ω −1 ci = 50 for M s = 3, β = 50, and θ Bn = 73 • shock in Ha et al. (2021) and t max Ω −1 ci = 78 for the case with M s = 3, β = 5, and θ Bn = 75 • in Kobzar et al. (2021).Compared to these studies, our hybrid simulations are very-large scale, with L y = (48 − 192)λ i , covering physical ranges from 5.9r g to even 47r g (see seventh and eighth columns in Table 1).Some of the simulations are also evolved to strongly nonlinear evolution stages, t max ≫ 50Ω −1 ci (ninth column in Table 1).We also systematically explore a wide range of physical parameters, including plasma beta, the shock Mach number, and the magnetic field obliquity.In particular, investigations of the θ Bndependence over such a wide range of angles have not been undertaken before.

Evolution of the shock structure
This section presents the shock structure characteristics typical of a shock with sonic Mach number M s = 3 and a quasiperpendicular obliquity θ Bn = 75 • .Our discussion is based on simulation run A4 with plasma beta β = 20, which in the following is referred to as the reference run.A general shock structure is presented in Section 3.1.Section 3.2 discusses the linear analysis of the magnetic wave modes, Section 3.3 treats density perturbations and shock ripples, and Section 3.4 describes the nonlinear evolution of the turbulence in the shock.

General structure of a reference shock in high-beta plasma
Figure 1 shows the profiles of various physical quantities across a well-evolved shock at time Ω ci t = 30, averaged over the y- coordinate.The x-axis is shown in units of the ion gyroradius (Eq.7), which for the specific parameters of run A4 is r g ≈ 4.1λ i .
The approximate location of the shock front at x ≈ 15.5r g is shown in Figure 1 with a vertical blue dotted line.The corresponding ion velocity phase-space is presented in Figure 2. The normalized ion density profile, shown with black line in Figure 1a, exhibits a characteristic pattern of overshootundershoot oscillations past the shock front.The density compression at the double-peaked overshoot (at x ≈ 15.2r g and x ≈ 14.5r g ) reaches N i /N 0 ≈ 3.5, decreases to N i /N 0 ≈ 2 in the undershoot at x ≈ 13.8r g , and increases again to N i /N 0 ≈ 3.1 in the second overshoot at x ≈ 13.0r g .Further downstream, the density compression gradually relaxes, eventually reaching the Rankine-Hugoniot value, approximately r ∼ 3.This observed structure is a consequence of the reflection of a portion of the incoming ions from the shock back upstream, a phenomenon evidenced by the presence of ions with large positive values of v xi and v zi in the phase-space ahead of the shock at x/r g ≈ 15.5 (Fig. 2a and 2c).The reflected ions gyrate in the upstream magnetic field and acquire energy as they drift along the motional electric field, E 0 .This energy increase enables the ions to surpass the potential drop across the shock and be transported downstream after a single reflection.The reflected ions turn around back to the shock within approximately 1r g and further gyrate in the downstream shock-compressed magnetic field.As the latter is nearly perpendicular to the shock normal (θ Bn ≈ 85 • , compare Fig. 1e), the ion gyration occurs primarily in the x − z plane, giving rise to characteristic patterns in downstream x − v xi and x − v zi phase-spaces.Such shock structures are a prevalent characteristic of quasi-perpendicular shocks (Treumann 2009).
The reflection and gyration of ions give rise to a marked ion temperature anisotropy at the shock ramp and overshoot.Figure 1b shows y-averaged profiles of the diagonal components of the ion temperature tensor, T x , T y , and T z , normalized to T 0 .The temperature tensor is derived from the pressure tensor, the second-order moment of the ion distribution function, calculated in the local plasma rest frame.In Figure 1c, we present the temperature components parallel, T ∥ , and perpendicular, T ⊥ , to the local magnetic field, and we show the temperature ratio, T ⊥ /T ∥ , in Figure 1d. 1 The anisotropy is most pronounced at the shock ramp and overshoot, T ⊥ /T ∥ ≫ 1, and its magnitude decreases with distance from the shock, while the amplitude of T ⊥ oscillations declines and T ∥ increases.This results from the pitch-angle scattering of ions off magnetic turbulence in the shock transition, which isotropizes the ion distribution, as one can discern in phase spaces at x/r g ≲ 13 in Figure 2.
The ion temperature anisotropy sources the magnetic turbulence in the shock.Figure 1a shows the profiles of B x , B y , and B z magnetic field components.The corresponding maps of the normalized ion density, N i /N 0i , and magnetic field fluctuations, δB x , δB y , δB z , at Ω ci t = 30 are presented in Figure 3b.The field fluctuations are defined with respect to the flux-frozen magnetic field, as proposed in Guo et al. (2017): where N i (x) is the y-averaged ion density.They are normalized to the upstream magnetic field, B 0 .The amplitudes of all components of the field fluctuations are similar and substantial enough in the shock front to induce departures from the frozen-in prediction.The latter can be discerned through deviations of the B y profile (red line) from the density profile (a proxy for the flux-freezing prediction) in Figure 1a.Deviations are visible in the first overshoot and are most pronounced around the second overshoot, x/r g ≈ 11−14, where the amplitudes of all field components reach their maximum values.Figures 3a-c display snapshots of the shock structure as it evolves at times Ω ci t = 10, 30, and 55, respectively.The ion density maps reveal corrugations in the shock transition present at the shock front, as well as in the region of the second overshoot.These shock ripples appear quite early in our simulation, approximately at time Ω ci t ≈ 6.Their amplitude and wavelength progressively increase over time.By the time Ω ci t = 30, the shock has transitioned into the strongly nonlinear stage.The amplitude of the shock ripples reaches δN i /N 0 ≈ 5 and remains at this level by the end of the simulation.Strong density and the magnetic field fluctuations in the shock are confined to the region of approximately 6r g in width from the shock front (compare Figs. 1-2).Farther downstream, the system relaxes to a weakly turbulent state.These characteristics of the ion-scale shock structure are consistent with results for subluminal shocks with M s = 3 and β ∼ 20 obtained with PIC simulations (Guo et al. 2017;Ha et al. 2021).

Linear analysis of the wave turbulence
The ion temperature anisotropy T ⊥ /T ∥ > 1 at the shock and downstream should excite the Alfvén Ion Cyclotron (AIC) instability and the mirror instability (e.g., Winske & Quest 1988;McKean et al. 1995;Lowe & Burgess 2003).Figures 3a2-3a4 show magnetic fluctuations at time Ω ci t = 10, when the shock is at the end of the linear stage.Dominant waves in B x and B z have their wavectors primarily in the y-direction and so nearly aligned with the background compressed magnetic field.They are circularly polarized and move downward along the mean field direction with velocity equal to ∼ 3.7 of the Alfven velocity in the overshoot (see below that this is compatible with the linear analysis).These waves are then consistent with AIC modes.Fluctuations in B y have oblique wave vectors that are slightly weaker, which are the characteristics of the mirror mode.The mirror modes should have zero real frequency, but in our simulation, they move in the same direction as the AIC modes.The asymmetric wave propagation in directions parallel and anti-parallel to the magnetic field has been observed except for strictly perpendicular shocks, indicating that a finite parallel drift of the reflected ions associated with obliquity plays a role.Although the linear analysis presented below assuming a symmetric ion distribution gives a reasonable agreement with the simulations in terms of wavelength, we may need to consider the finite asymmetry to understand the downward wave propagation of AIC and mirror modes.This requires a more detailed analysis of the obliquity angle dependence and is left for future study (see Section 4.3).
Figure 4 shows Fourier power spectra of the magnetic waves in the region from the shock ramp to behind the second overshoot, 1 ≲ x/r g ≲ 5.5 at Ω ci t = 10.The maximum wave power in B x and B z waves is at k ∼ k y ≈ 0.9/λ i , which corresponds to λ ≈ 7λ i .These waves are present in the undershoot and also in the first overshoot, where they co-exist with smaller-scale fluctuations with k y λ i ≈ 1.1 and 1.3 (λ ≈ 6.3λ i and λ ≈ 4.4λ i , respectively), and the (k x , k y )λ i ≈ (0.8, 0.9) mode in B z with wavelength λ ≈ 5.2λ i .The B y waves in the same region also have k y λ i ≈ 0.9, 1.1 and 1.3 wavevectors, though Fig. 3: Structure of the reference shock with M s = 3, β = 20, and θ Bn = 75 • (run A4) at different times Ω ci t = 10 (panels a * ), Ω ci t = 30 (panels b * ), and Ω ci t = 55 (panels c * ).From top to bottom, shown are the maps of the normalized ion density, N i /N 0i (panels *1), and the normalized magnetic field fluctuations (see Equations 8-10), δB x (panels *2), δB y (panels *3), and δB z (panels *4).The scaling of all distributions is logarithmic and for the magnetic fields it is sign-preserving, and, e.g., for δB z it is: sgn(δB z ) • {2 + log[max(|δB z |/B 0 , 10 −2 )]}.The level of "0" on the color scale hence corresponds to |δB|/B 0 ≤ 10 −2 (see, e.g., Kobzar et al. 2021).with small k x λ i < 0.7 components.These waves are well visible in the δB y map in Figure 3a3, but not in the spectrum in Fig- ure 4b, because they are much weaker than the waves present at x/r g ≈ 3, that dominate the Fourier power spectrum.The peak signal in B y waves in this second overshoot region is at (k x , k y )λ i ≈ (0.75, 0.6), so that λ ≈ 6.8λ i .Significant wave power is also in oblique modes with (k x , k y )λ i ≈ (0.75, 0.75), (0.8, 0.9), and (1.0, 1.0), with the respective wavelengths of λ/λ i ≈ 5.6, 5.2 and 4.4.Finally, in the second overshoot, slightly longer wavelength modes exist in B x and B z with k y λ i ≈ 0.75 (λ ≈ 8.4λ i ).
We have performed a linear dispersion analysis to identify the magnetic waves in the shock transition with the AIC and mirror modes.The model is similar to that used in Kim et al. (2021) and assumes a bi-Maxwellian ion velocity distribution function.Because in the shock, the ion density and velocity distributions are highly nonuniform, we adopt the numerical values of the model parameters -the parallel ion beta β i∥ ≡ β ∥ and the temperature ratio, T ⊥ /T ∥ -measured at two locations across the shock: the shock ramp at x/r g ≈ 5.5 and the region between the undershoot and the second overshoot, 3 ≲ x/r g ≲ 3.5, at time Ω ci t = 10 (compare Fig. 3a).In the shock ramp, the temperature ratio is high, T ⊥ /T ∥ ≈ 7, β ∥ ≈ 9.3 and plasma is weakly compressed, r ≈ 1.5.On the other hand, in the undershoot-second overshoot region, the temperature ratio is moderate, T ⊥ /T ∥ ≈ 2.7, β ∥ ≈ 23 and plasma compression r ≈ 2.7.By calculating the linear properties in these two regions, we estimate the range of wavevectors and growth rates that should be observed in our simulation.In Appendix C we compare the linear theory predictions for the parameters in the shock ramp with a simulation in the periodic-box.We demonstrate that linear theory results are satisfactorily reproduced.
Figure 5 shows the linear analysis results for both regions.Parallel, k ∥ , and perpendicular, k ⊥ , wave vectors are defined with respect to the average magnetic field so that k ∥ corresponds to Fig. 4: Fourier power spectra of the magnetic waves in the region from the shock ramp to behind the second overshoot, 1 ≲ x/r g ≲ 5.5 at Ω ci t = 10 for run A4 (compare Fig. 3a).Note different dynamic scales in each panel.also has a longer wavelength of λ ≈ 7.6λ i , corresponding to (k ∥ , k ⊥ )λ ′ i ≈ (0.35, 0.35).We conclude that the predicted range of the AIC and mirror modes wavelengths is consistent with waves with λ/λ i ≈ 4 − 8 observed in the shock transition.The estimated growth time scales, τ ≡ 1/γ, of the AIC modes, τ AIC ≈ (0.5 − 0.9)/Ω ci , and the mirror modes, τ mirror ≈ (1.0 − 1.9)/Ω ci , are also in good agreement with the fact that both modes develop very early in our simulation.
We note finally that the ion-scale magnetic field turbulence structure observed in our hybrid simulations agrees with the fully kinetic PIC simulations performed for the same physical parameters.A detailed comparison is provided in Appendix D. The alignment found in the results derived from the hybrid-kinetic method and the PIC simulations confirms the reliability of hybrid simulations for investigating ion-scale shock physics.

Plasma density perturbations and shock ripples
The lack of electron-scale waves in hybrid simulations in the linear stage enables us to witness the initial evolution of plasma density fluctuations in the shock front (compare Appendix D).One can see in Figure 3a and also in Figure D.1d, that perturbations in ion density on top of the shock plasma compression appear very early in the reference shock simulation.This is also true for other plasma betas investigated here for shocks with M s = 3 (compare, e.g., Fig. D.1b for β = 5).Previous work suggested that shock rippling is associated with AIC instability.However, since this instability is purely electromagnetic and transverse, generating density compressions must involve magnetic perturbations with components parallel to the mean magnetic field, such as the oblique mirror modes.The characteristics of the ripple wave structure and evolution observed in our simulations are in agreement with other studies suggesting that density fluctuations at the shock front are the effects of the nonlinear evolution of AIC and mirror modes, along with nonlinear couplings between these modes (e.g., Winske & Quest 1988;Lowe & Burgess 1999;Shoji et al. 2009;Kim et al. 2021).
Earlier PIC simulations of shocks in high-beta plasmas by Kobzar et al. (2021) and Ha et al. (2021) reported much slower growth rates and much longer wavelengths of the initial shock front ripple modes than observed in our hybrid simulation.For example, Kobzar et al. (2021) for M s = 3 and β = 5 shock observed the ripples to emerge at Ω ci t ≈ 25 with wavelength λ ripple ≈ 16λ i , in agreement with their linear dispersion analysis that predicted the peak growth rate of the AIC modes of γ/Ω ci ≈ 0.076.As we demonstrate in Section 3.2, in our hybrid reference run, the initial growth rates are much faster, and AIC mode wavelengths are much shorter than these values, which is also in line with our linear analysis that approximates numerical values at the overshoot.Moreover, we show in Appendix D that results of A4 and A2 simulations reproduce the early ionscale structures of the magnetic field fluctuations and plasma density in PIC simulations.However, note that there is no disagreement between these findings.The linear analysis by Kobzar et al. (2021) is based on a different model that takes into account the estimates of the number of ions transmitted and reflected at the shock, and the temperature ratio T ⊥ /T ∥ = 4.7 at an already evolved shock.On the other hand, the linear analysis model used by Kim et al. (2021) is similar to ours but takes the parameters from shock simulations averaged over an extended shock transition zone, e.g., T ⊥ /T ∥ = 2 for M s = 3 shocks in β = 5−100 plasmas.Consequently, these linear estimates provide growth rates and wavelengths matching the shock conditions at later nonlinear evolution stages, at which our hybrid simulations show comparable turbulence characteristics.

Nonlinear evolution of the shock structure
As noted above, in our reference simulation, we observe a transition of the shock system into a nonlinear phase already at Ω ci t ≈ 6.One can see in Figure 3 that with time, the magnetic field and plasma density fluctuations in the shock front grow in wavelength, and large-scale shock ripples develop.Figure 6 illustrates the temporal evolution of the magnetic field and density turbulence in the region of the forward peak of the first overshoot.Shown are the values of the k y wavevector component that is approximately aligned with the compressed magnetic field direction for individual components of the magnetic field as well as the ion density.Although the wave power spectra of turbulence reveal strong signals at multiple distinct wavevectors, we display only the values of the dominant k y (compare Fig. 4 for magnetic field data at Ω ci t = 10).One can note that with time the spectrum of magnetic AIC (primarily B x and B z components) and mirror (mainly B y component) waves shifts to lower wave numbers.Such energy transfer from higher to lower wave numbers is due to nonlinear effects, such as the parametric decay instability (Shoji et al. 2009).As time progresses beyond Ω ci t ≈ 30, the field components exhibit saturation at k y ≈ 0.2, corresponding to a wavelength of approximately 30λ i .
Similar evolution from short to long wavelengths followed by saturation is observed for density fluctuations.Past Ω ci t ≈ 20 the ripple wavelength is λ ripple ≈ 30λ i .The ripple structure is neither steady nor regular, though, as at locations along the shock surface, smaller ripples appear at times in addition to the largerscale ones (see Fig. 3b-c.However, the structures do not evolve toward even longer wavelengths by the end of our reference run, A4 at Ω ci t = 56.On the other hand, weak downstream turbulence beyond the second overshoot remains at a quasi-steady state at times Ω ci t ≳ 30.Note, that we observe similar behavior for M s = 3 shock propagating in β = 5 plasma (run A2), which we follow up to Ω ci t = 87.To the best of our knowledge, these characteristics of the nonlinear evolution and saturation are reported here for the first time in the context of weak shocks in high-beta plasmas.This is facilitated in the hybrid approach, which enables extended, large-scale simulations of shocks, a capability beyond the current constraints of PIC simulations.

Parameter dependence
In this section we discuss the dependence of the results obtained for the reference run A4 with For all shocks analyzed here, the motions of shock-reflected ions give rise to the ion temperature anisotropies that catalyze the development of the ion cyclotron and mirror modes.These waves scatter the ions in pitch angle and reduce their anisotropy, bringing it closer to the upper limit provided by the marginal stability condition discussed in Gary et al. (1997): where β ∥ = 8πN i k B T ∥ /B 2 is the proton plasma beta at a given location, calculated using the proton temperature in the parallel direction.with lower marginal stability limits, resulting in faster growth rates of ion-scale instabilities.Consequently, shocks in higherβ plasmas initially generate more intense turbulence that provides more efficient ion scattering, leading to a quicker decline in anisotropy behind the shock.Therefore, in high-β shocks, turbulence attains its peak and becomes localized closer to the shock front (Fig. 8b-c for β = 20 and 50).Conversely, lower-β shocks sustain a finite level of ion isotropy further downstream.The turbulence exists in a more extended downstream region, where the magnetic energy also reaches its maximum (Fig. 8a for β = 5).
The relationship between the plasma beta and the location of the peak magnetic energy is most evident in Figures 9 and 10 that present profiles of the magnetic field energy density at times Ω ci t = 26 and Ω ci t = 55, respectively.The maximum magnetic energy density is observed approximately 1r g behind the shock with β = 100, about 2r g for β = 50 and 20, to around 3r g for β = 5.Lower-β shocks develop weaker turbulence, Figure 8.By the early nonlinear simulation phase at Ω ci t = 26, a substantial difference in the total magnetic energy, spanning nearly one order of magnitude, is observed when the plasma beta varies by a factor of ∼ 30, Figure 9.The amplitude levels do not change in the strongly nonlinear stage; compare Figure 10a.From the early stage of the evolution, the magnetic turbulence is dominated by the B z field component, Figure 10b-d.The amplitude difference between B z and B x fluctuations is an artifact of the 2D geometry, while the smaller strength of B y waves is attributed to the lower growth rate of the mirror mode instability.These findings are consistent with PIC simulation results (Guo et al. 2017;Guo et al. 2018).
Runs A2 and A5, have been followed to late times Ω ci t ≥ 62.As in the reference run, we observe a saturation of the ripple growth.The nonlinear wavelengths of the ripples are λ ripple ≈ 16λ i and λ ripple ≈ 48λ i for β = 5 and β = 50, respectively.In units of the ion skindepth, the wavelengths thus grow with plasma beta, as expected.However, note that in units of the ion gyroradius, the nonlinear ripple wavelength is approximately independent of beta, and for the three cases with β = 5, 20 and 50 it is λ ripple ≈ 3.7 − 3.9r g .

Dependence on the sonic Mach number M s
Figure 11 shows maps of the ion density and B z magnetic field fluctuations for the example cases of shocks with the shock Mach number M s = 2, 2.25, and 5.5 in β = 20 plasmas with θ B n = 75 • .Figure 12 presents the profiles of the total energy density of the magnetic field turbulence for all simulation runs C for shocks in the Mach number range M s = 2 − 5.5.One can see in Figure 7b that the amplitude of the temperature anisotropy peak at the shock front is a factor of approximately 4 larger for a strong shock with M s = 5.5 compared to a weak shock with M s = 2.0.At the same time, the marginal stability limits at the shock are similar in each case.In the case of the M s = 5.5 shock, the anisotropy drops right at the shock front, whereas for lower-M s shocks, it decreases in gradually wider downstream regions.In effect, since larger temperature anisotropies lead to faster growth of instabilities and stronger turbulence, in a high-M s = 5.5 shock, the turbulence attains its peak within less than 1r g downstream from the shock front, Fig. 11c and Fig. 12 (magenta line).With decreasing Mach number, the waves grow slower and hence the magnetic energy density reaches maximum progressively farther from the shock front, from about 2.5r g for M s = 2.9 shock, through approximately 3r g , 4r g , and 6r g , respectively for shocks with M s = 2.7, 2.45, and 2.25.At the same time, the turbulence extends in wider downstream regions and   becomes weaker with lower M s .One can see in Fig. 12, that the amplitude of the magnetic field is two orders of magnitude smaller for the M s = 2 shock compared to the M s ≈ 3 shock and more than three orders of magnitude smaller compared to the strong M s = 5 shock (compare also Fig. 11).
One can see in the upper panels of Figure 11 and in Figure 3, that higher-M s shocks, M s = 3, 5.5, develop distinct shock front ripples.Density corrugations in the shock front, in particular in the backward peak of the first overshoot, are also visible in the M s = 2.25 shock (Fig. 11b1), whereas the M s = 2 shock surface  of the second overshoot (∼ 5r g behind the shock) and further downstream.
Based on similar observations from PIC simulations, Ha et al. (2021) suggested that shock surface ripples appear only in shocks with the sonic Mach number above the critical Mach number M * AIC ≈ 2.3.At and above this M s value, the temperature anisotropy averaged over a fixed width in the shock transition, encompassing the first and the second overshoots, significantly exceeded the marginal stability condition, as indicated on the right-hand side of Equation 11.This coincided with an enhanced ion reflection efficiency from the shock at shocks with M s ≳ 2.3.At shocks with M s < M * AIC , the temperature anisotropy was T ⊥ /T ∥ ≈ 1 on average, and hence below the stability limit.
In effect, the AIC instability was not excited in these shocks.The sonic Mach number at which the shock ripples appear in our hybrid simulations is consistent with these findings.However, we note that in all our simulation runs C, the temperature anisotropy is markedly above the marginal stability condition (e.g., T ⊥ /T ∥ ≳ 2 within ∼ 5r g of the M s = 2 shock front) and relaxes via pitch angle scattering only far downstream of the second overshoot for very weak shocks (compare Fig. 7b).Because the ion temperature anisotropy in shocks with M s ≲ 2.3 is small, the resulting growth time scales of the AIC and mirror instabilities are slower than the plasma advection time scales through these shocks.Consequently, the turbulence is built far from the shock front and the shock surface remains laminar.Small amplitude of magnetic field fluctuations in very weak shocks results from a low amount of free energy stored in temperature anisotropies at these shocks.reflection from the oblique shock, the ions acquire a finite drift in a direction parallel to the mean magnetic field that contributes to T ∥ .The apparent temperature anisotropy decreases as the parallel drift speed increases for smaller θ Bn .An additional factor contributing to lower anisotropy at less oblique shocks is the ion drift in the motional electric field, the strength of which decreases with sin(θ Bn ).

Dependence on the obliquity angle θ Bn
The patterns in the temperature anisotropy profiles observed with an increasing θ Bn resemble those seen in temperature anisotropy profiles as Mach numbers progress from small to large values.Consequently, more oblique shocks should generate faster-growing and stronger waves, whose amplitude peaks closer to the shock front, compared to shocks with smaller θ Bn .One can see such a trend in the B z maps in Figure 13 and in the magnetic energy density profiles in Figure 14.However, despite a temperature anisotropy difference of approximately 3.5 between shocks with θ B n = 55 • and θ B n = 89 • , there is only a slight difference, by a factor of a few, in the strength of the magnetic field fluctuations in these cases.This might be because the free energy is provided not only by the temperature anisotropy but also by parallel ion drift.Consequently, a smaller anisotropy may not necessarily lead to a much lower level of turbulence.
All shocks we analyze here develop density corrugations in the shock surface.Due to the slower growth of instabilities, the shock ripples emerge at later times at lower-θ Bn shocks.The ripple wavelength becomes longer as the obliquity decreases.At the final simulation time Ω ci t = 56 for β = 20 shocks, Figure 13, the ripple wave in the strictly-perpendicular shock, θ Bn = 90 • , has a wavelength λ ripple ≈ 25λ i , in θ Bn = 75 • shock λ ripple ≈ 30λ i (see also Section 3.4), and in the shock inclined at angle θ Bn = 55 • the ripple wavelength is λ ripple ≈ 50λ i .This dependence of λ ripple on the field obliquity possibly results from the nonlinear evolution of the Alfvenic turbulence in the overshoot that couples the ion temperature anisotropy-driven modes with waves generated due to the parallel ion drift.Such dependence has not been reported before.We leave a detailed analysis of these effects to a future study.
We have demonstrated in Section 3.4 that the wavelength of the ripple mode saturates past the simulation time Ω ci t ≈ 20.We see a similar trend in the simulation of the strictly-perpendicular shock.As in the case with θ Bn = 75 • shock, the ripple wave structure is quasi-steady with smaller-scale density fluctuations occasionally appearing on top of the dominant mode.In this case, we also observe symmetric wave propagation parallel and anti-parallel to the mean magnetic field (see Section 3.2).In the case of the oblique shock with θ Bn = 75 • , the ripple wavelength evolves slowly.Confirming the saturation in this case necessitates an extended simulation time and possibly a larger computational box.

Comparison with 3D simulations
Our studies in this work are mostly based on 2D simulations, since running 3D simulations encompassing large ion-scales of the shock-front ripples is computationally demanding, even in the hybrid approach.Therefore, to study the dimensionality effects, we performed two runs for shocks with Mach number M s = 3, obliquity angle θ Bn = 75 • and plasma beta β = 5 and 10, runs E1-3D and E2-3D, respectively.As in 2D simulations, the upstream magnetic field lies in the x − y plane.The size of the simulation box along the y and z directions is equal, L y = L z .All other simulation parameters are the same as in 2D runs.The E1-3D run serves as a reference for studies of large-scale effects, while smaller-size E2-3D simulation provides supportive insights.
Figure 15 compares maps of the ion density and magnetic field fluctuations for the 2D run A2 (panels a*) and run E1-3D (panels b*) with β = 5 at time Ω ci t = 36.Maps for run E1-3D show a 2D cut across the computational box in the plasma flow direction in the middle of the box at z = 48λ i .One can note a very good correspondence in the density and wave mode structures between 2D and 3D.In particular, the wavelength λ ripple ≈ 12 − 14λ i of the shock ripples in E1-3D run is in agreement with the ripple wavelength measured in 2D run A2 at the same simulation stage.
The cut through the box in the direction transverse to the plasma flow at x = x sh − 0.12r g presented in Figure 16, shows the structure of the shock surface ripples and associated nonlinear AIC modes.It is apparent that the turbulence structure is 2D-like, and the waves have only weak long-wave components along the z-direction.We note that in the early-stage of the system evolution in the shock front region bounded by the first and the second overshoots, there is a significant wave component to δB x , δB y and δB z fluctuations with k z ≈ k y , and the waves are oblique in general (not shown).However, a steady-state structure in this region is 2D-like, and only far downstream plasma hosts quasi-isotropic turbulence.Similar characteristics are observed in run E2-3D with β = 10.
Finally, Figures 17a and 17b compare the profiles of the normalized energy density of the magnetic fluctuations between 2D (solid lines) and 3D (dash-dotted lines) simulations at Ω ci t = 24 and Ω ci t = 36, respectively.At an earlier stage, Figure 17a, the magnetic field fluctuation levels in 2D and 3D are comparable in the case with β = 10.At the same time, shocks with β = 5 differ between 2D and 3D by nearly an order of magnitude in the average magnetic energy density in the region close to the shock front.Similar magnetic turbulence levels are nevertheless reached by Ω ci t = 36 in this case, Figure 17b, because the wave growth rates are lower for smaller β = 5, compared to β = 10 (see Section 4.1).
The results presented in this section indicate that the steadystate structure of weak shocks in high-beta plasmas does not considerably differ between fully 3D simulations and simulations with dimensionality restricted to 2D with an in-plane magnetic field geometry.While this conjecture may not hold for strong shocks, like the case with M s = 5.5 considered here in 2D (see Section 4.2), it is expected to be valid for shocks with M s ≲ 3, that are considered more representative of the merger shocks in ICM.Further studies are necessary to demonstrate conclusively that the effects of higher plasma beta or different magnetic field obliquity do not alter this conclusion.

Summary and concluding remarks
In this study we employ hybrid-kinetic numerical simulations to investigate the ion-scale kinetic instabilities at merger shocks of galaxy clusters.Our results can be summarized as follows: • • For the first time we also investigate high-β quasiperpendicular shocks in a wide range of the magnetic field obliquity.Despite significant variations in temperature anisotropy at the shock, we only observe a moderate differences in the strength of magnetic field fluctuations at various θ Bn .This indicates that, in addition to temperature anisotropy, parallel ion drift plays a significant role in driving magnetic turbulence.We also observe that shock ripples have longer wavelengths as the field obliquity decreases.Our results indicate that the ripple growth saturation occurs independent of θ Bn .However, to confirm this conjecture, longer simulation times and potentially larger computational boxes are required for shocks with low θ Bn .• The comparison between 2D and 3D simulations suggests that the steady-state structure of merger shocks with M s ≲ 3  3).The shock position, x sh , is marked with the white dotted line.
in high-beta plasmas does not significantly differ between fully 3D and 2D simulations with an in-plane mean magnetic field.The shock surface ripples and associated nonlinear Alfven ion cyclotron (AIC) modes exhibit a 2D-like structure at the shock front and immediate downstream, with quasi-isotropic turbulence emerging only in the far downstream region.Steady-state shocks in both 2D and 3D develop comparable levels of turbulence.
The studies in this work provide the fundamental framework for conducting comprehensive large-scale 3D simulations aimed at understanding shock physics in a wide range of physical conditions.They may also guide parameter choice for large-scale PIC simulations.Equally important is that they lay the foundation for further investigations of the impact of ion-scale turbulence on electron acceleration in high-beta shocks by means of test-particle simulations (e.g., Trotta & Burgess 2019) or similar methods.As emphasized in Section 1, multi-scale turbulence is vital for efficient electron acceleration via the SSDA process.Our recent PIC simulations (Kobzar et al. 2021) indicate that electron-scale turbulence has less impact on electron energization to supra-thermal energies.Hence, the investigation of electron acceleration can be reliably conducted through test particle simulations.Our generalized fluid particle hybrid code is particularly well-suited for this purpose, as it allows for the incorpo- appearance of temperature anisotropy.Subsequently, the system further relaxes toward complete ion distribution isotropy.In this regard, periodic box simulation results significantly differ from the shock simulations.In the latter, a continuous inflow of fresh ions from the upstream sustains high-temperature anisotropy at the shock, leading to the excitation of instabilities.Consequently, meaningful comparisons between the two simulations are viable only during the early evolution of the periodic-box system, for Ω ci t ≲ 10.
In Figure C.1, it is evident that the growth rate of transverse magnetic field fluctuations, δB y and δB z , associated primarily with AIC modes, significantly exceeds that of δB x waves, corresponding to mirror modes.During the time interval Ω ′ ci t = 2.5 − 5, the measured growth rates are γ AIC ≈ 0.8Ω ′ ci for δB y , δB z components and γ m ≈ 0.4Ω ′ ci for δB x .The slower growth of mirror modes aligns with linear theory.However, both growth rates are notably smaller than the linear predictions of γ AIC ≈ 1.2Ω ′ ci and γ m ≈ 0.6Ω ′ ci .This suggests that, even at this early stage, nonlinear effects are at play, which is unsurprising given the challenge of capturing a truly linear evolution when the wave growth rate is very high.Nonetheless, we believe comparing early-stage periodic-box turbulence with its characteristics observed at the shock ramp is viable.ci t ≤ 5.The Fourier power in B y and B z waves is maximum at k ∥ λ ′ i ≈ 0.7 − 1.4, in reasonable agreement with linear theory (compare Figure 4a).The corresponding range of wavelengths calculated using the shock simulation values is λ ≈ 3.7 − 7.3λ i , which is consistent with the AIC waves observed at the first overshoot (compare Figure 4a and 4c).Slightly oblique waves in B x with wave vectors k ⊥ λ ′ i ≲ 0.25 and similar k ∥ λ ′ i ≈ 0.6 − 1.4 have also their counterparts in the waves observed at overshoot in δB y fluctuations (see description in Section 3.2).Although these waves are most likely due to nonlinear effect of AIC and mirror waves coupling, the oblique modes with k ∥ λ ′ i ≈ 0.6 − 1.4 and k ⊥ λ ′ i ≈ 0.4 − 0.7 seen in δB x and also δB y are consistent with the predictions of the linear analysis.We conclude that periodic-box simulations satisfactorily reproduce both the linear theory predictions and the properties of waves observed in the first overshoot in the early-stage evolution of the shock.
During subsequent evolution, we observed waves in all magnetic field components growing in wavelengths (not shown).This is in line with the known nonlinear evolution of AIC and mirror waves (e.g., Shoji et al. 2009).

Appendix D: Comparison with 2D PIC simulations
In Figure D.1, we present a comparison of 2D results obtained for the same shock parameters using two different simulation approaches: particle-in-cell (PIC) simulations in the left panels (a) and (c) and hybrid simulations in the right panels (b) and (d)).We compare the early stage structures of shocks with Mach number M s = 3, mean magnetic field inclination of θ Bn = 75 • , and two different plasma betas: β = 5 in top panels and β = 20 in bottom panels.PIC simulations results for the case of β = 5 are taken from Kobzar et al. (2021).The results for β = 20 have not yet been published (Kobzar et al. (2024), in preparation).We present hybrid simulations results from runs A2 and A4, though for a clearer one-to-one comparison, we display only portions of the computational boxes for hybrid runs in  results are considered to be contemporaneous.PIC simulations have been performed with 20 particles per cell per species, ionto-electron mass ratio m i /m e = 100, and the electron skin depth λ e = c/ω pe = 15 cells, where ω pe = 4πe 2 N e /m e is the electron plasma frequency.
In the linear phase of PIC simulations, the prevailing turbulence in the first overshoot is the electron whistler waves, evident as short-scale magnetic field fluctuations at 18 ≲ x/λ i ≲ 24 in As the PIC and hybrid simulations progress, the ion-scale waves evolve into multi-scale turbulence.This includes the development of long-wave shock-front rippling modes, whose nonlinear wavelengths closely match.For instance, the wavelength of λ ripple ≈ 16λ i of the ripples observed at time Ω ci t = 36 in run A2, closely fits the shock corrugation wavelength reported in Kobzar et al. (2021).The agreement between the results obtained with the hybrid-kinetic approach and our PIC simulations, as well as other PIC simulation studies (e.g., Ha et al. 2021;Guo et al. 2017;Guo et al. 2018) attests to the validity of the hybrid simulations in exploring ion-scale shock physics and the role of ion-scale turbulence in particle acceleration.

Fig. 1 :
Fig. 1: Structure of the reference shock with M s = 3, β = 20, and θ Bn = 75 • (run A4) at time Ω ci t = 30.Displayed are the y-averaged profiles of (a) the ion density normalized to the far upstream density N 0i (black line), B x , B y and B z magnetic field components normalized to their upstream values (blue, red, and green lines, respectively), (b) profiles of T x , T y , and T z temperature components (red, blue, and black lines, respectively), and (c) the components parallel, T ∥ (red line), and perpendicular, T ⊥ (cyan line), to the local magnetic field, normalized to T 0 , (d) the T ⊥ /T ∥ temperature ratio, and (e) the mean magnetic field inclination angle with respect to the shock normal, θ Bn .The vertical dotted line in each panel marks the shock location.The x coordinate is normalized to the ion gyroradius r g .
Fig. 5: Linear growth rate for the AIC and mirror modes for conditions representing (a) the shock ramp and (b) the region between the shock undershoot and the second overshoot in run A4.Prime denotes quantities measured in the local plasma frame.

Fig. 6 :
Fig. 6: Temporal evolution of the dominant k y wavevector component of the magnetic field and ion density fluctuations in the shock front for the reference simulation with M s = 3, β = 20, and θ Bn = 75 • (run A4).The error bars reflect resolution in the Fourier space (see text).
M s = 3, β = 20, and θ Bn = 75 • on the shock parameters: plasma beta in the range β = 5 − 100 in Section 4.1, the shock Mach number in the range M s = 2 − 5.5 in Section 4.2, and the range of the upstream magnetic field inclination angle, θ B n = 55 • − 90 • , in Section 4.3.
Figure 8 shows maps of the ion density and B z magnetic field fluctuations for the selected cases of β = 5, 20, and 50, for shock Mach number M s = 3, and θ B n = 75 • .One can see in Figure 7a, that the amplitude of the temperature anisotropy peak at the shock front remains relatively unchanged, regardless of the plasma beta.At the same time, higher β values are associated

Fig. 8 :
Fig. 8: Structure of shocks with M s = 3 and θ Bn = 75 • at Ω ci t = 55 for different values of plasma beta, β = 5, 20, 50, runs A2, A4, and A5 in panels a-c, respectively.Upper panels (*1) show distributions of normalized ion density, and lower panels (*2) the normalized B z magnetic field fluctuations.The scaling is linear (compare Fig. 3c for the case with β = 20, run A4).White dotted lines depict the approximate location of the shock, x sh , with respect to which distance is calculated in Figures 9 and 10.
Fig. 8: Structure of shocks with M s = 3 and θ Bn = 75 • at Ω ci t = 55 for different values of plasma beta, β = 5, 20, 50, runs A2, A4, and A5 in panels a-c, respectively.Upper panels (*1) show distributions of normalized ion density, and lower panels (*2) the normalized B z magnetic field fluctuations.The scaling is linear (compare Fig. 3c for the case with β = 20, run A4).White dotted lines depict the approximate location of the shock, x sh , with respect to which distance is calculated in Figures 9 and 10.

Fig. 9 :
Fig. 9: The y-averaged profiles at time Ω ci t = 26 of the normalized total energy density of the magnetic field fluctuations across the shocks with M s = 3 and θ Bn = 75 • and plasma beta in the range β = 3 − 100, runs A1-A6.The x coordinate is normalized to the ion gyroradius r g and measured relative to the shock location x sh .

Fig. 10 :
Fig. 10: The y-averaged profiles at time Ω ci t = 55 of the normalized energy density of the magnetic field fluctuations (total and in B x , B y , and B z components) across shocks with M s = 3 and θ Bn = 75 • and plasma β = 5, 20, and 50, runs A2, A4, A5.

Figure 13
Figure 13 presents maps of the ion density and B z magnetic field fluctuations for shocks with M s = 3 in β = 20 plasmas, and different θ Bn angles, 55 • , 75 • , and 90 • .Figure 14 shows profiles of the total energy density of the magnetic field turbulence for M s = 3 shocks in the range θ Bn = 55 • − 89 • and for a lower value of plasma beta, β = 10.We note, that all shocks are quasiperpendicular, θ Bn > 45 • , and all but the strictly-perpendicular one, θ Bn = 90 • , are subluminal shocks with θ Bn < θ Bn,crit , where the critical superluminality angle is θ Bn,crit = arccos(u sh /c) ≈ 89.93 • for shocks in β = 20 and θ Bn,crit ≈ 89.95 • for β = 10.The variations in the magnitude of the temperature anisotropy at the shock front with the field obliquity seen in Figure 7c are related to the geometry of the ion reflection from the shock and their gyration in the upstream magnetic field.Upon
The structure of a subluminal large-scale 2D shock with fiducial merger shock parameters, M s = 3, β = 20 and θ Bn = 75 • , is governed by the ion temperature anisotropies in the shock due to the ion reflection from the shock and their gyration in the magnetic field in the shock transition.These anisotropies drive multi-scale magnetic turbulence whose properties in the linear stage are consistent with AIC and mirror modes and agree very well with characteristics observed in fully kinetic PIC simulations.This confirms the credibility of hybrid simulations in exploring the ion-scale physics of shocks.•With growing plasma β shocks generate stronger turbulence that is confined closer to the shock front.Shocks in very high-beta plasmas, β ≈ 50 − 100, produce an order of magnitude more intense turbulence than lower-β shocks.• Nonlinear evolution of AIC and mirror modes, and the nonlinear couplings between these modes are responsible for the formation of plasma density fluctuations in the shock front.Initial short-wave ripples evolve toward long wavelengths.We show for the first time that in M s = 3 shocks in high-β plasmas, β = 5 − 50, the growth of the shock front corrugations and magnetic waves saturates.The shock ripples reach wavelengths of λ ripple ≈ 3.7 − 3.9r g .This result has important consequences for electron acceleration processes at cluster shocks.The saturated wavelengths of the waves impose constraints on the maximum energy of particles that are confined at the shock through resonant scattering off the waves and undergo the SSDA process.• At very weak shocks with M s ≲ 2.3 the AIC growth time scale is slower than the plasma advection time scale through the shock.Therefore, such shocks do not develop shock surface ripples and form weak turbulence in the shock downstream.This is consistent with the notion of the critical Mach number M * AIC ≈ 2.3 to trigger multi-scale waves in the shock front, proposed inHa et al. (2021).

Fig. 15 :
Fig. 15: Structure of the M s = 3, β = 5, θ Bn = 75 • shock at Ω ci t = 36 in 2D and 3D simulation runs A2 (panels a*) and E1-3D (panels b*).Two-dimensional maps of the normalized ion density and magnetic field fluctuations are shown, and for run E1-3D the cut at z = 48λ i is displayed (see Fig.3).The shock position, x sh , is marked with the white dotted line.

Fig. 16 :
Fig. 16: Maps of density and B x magnetic field fluctuations in the y − z plane at x = x sh − 0.12r g (see Fig. 15) at Ω ci t = 36 for E1-3D run.
Figure C.2 shows distributions of δB x , δB y , and δB z magnetic field fluctuations at time Ω ′ ci t = 4. Figure C.3 displays corresponding Fourier power spectra averaged in the time interval 3 ≤ Ω ′ Figure D.1.The shock structures are compared at slightly different times due to the distinct nature of the two simulation codes.Nevertheless, the

Fig
Fig. C.1: Time evolution of the magnetic field fluctuation energy densities (top), the temperature anisotropy (center), and the parallel and perpendicular temperature components (bottom) for the periodic-box simulation.
Figure D.1a and at 41 ≲ x/λ i ≲ 46 in Figure D.1c (middle and bottom panels).However, strong ion-scale fluctuations grow in parallel and emerge downstream of the first overshoot.The characteristics of these waves in δB x (not shown) and δB z components are consistent with AIC instability modes, and in the δB y component with the mirror modes.The structure of these fluctuations closely resembles the waves observed in hybrid simulations in the same shock regions.In particular, in the β = 5 shock, the dominant wavelength of the nearly transverse AIC fluctuations is approximately λ AIC ≈ 5λ i , and the field-aligned component of the oblique mirror modes, λ m∥ ≈ 6λ i , in both cases.Similarly, the wavelengths of both the AIC and mirror waves are λ AIC ≈ λ m∥ ≈ 7.5λ i in the shocks in β = 20 plasmas.

Fig
Fig. C.2: Maps of the normalized magnetic field fluctuations for the periodic box simulation at time Ω ′ ci t = 4. From left to right shown are (a) δB x = B x − B 0 , (b) δB y , and (c) δB z .

Fig
Fig. C.3: Fourier power spectra of magnetic field fluctuations in the periodic box simulation averaged in the time interval 3 ≤ Ω ′ ci t ≤ 5. Note that the dynamic range is different in each panel.

Table 1 :
Parameters of the hybrid-kinetic simulation runs discussed in this work.