© The Authors 2024

1. 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, which are 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 megaparsec-sized regions of radio emission, typically with an arc-like morphology. This emission shows strong polarization and steep integrated radio spectra (L_ν ∝ ν^α with α ∼ −1.1) that steepen even further from the cluster edge toward its center (van Weeren et al. 2019). The high polarization in radio relics is attributed to the compression of the perpendicular component of preshock magnetic fields across shocks. Conversely, radio halos, which are formed by electrons accelerated through turbulence, exhibit notably lower polarization. Furthermore, spectral steepening within these relics primarily arises from the postshock aging of cosmic ray electrons due to synchrotron and inverse Compton losses. This allows us to gain insights into the interaction between magnetic fields and high-energy particles.

Merger shocks are found to propagate in the intracluster medium (ICM), which is hot (with temperatures reaching 10⁸ K), dilute (plasma density 10⁻³ − 10⁻⁵ cm⁻³), and weakly magnetized, with the magnetic field strength on the order of 1 μG. Correspondingly, ICM plasma is characterized by a high plasma beta, β ≫ 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 first-order 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 the hot ICM, such mechanisms are 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 the 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, 2014b; 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 it 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 quasi-perpendicular subluminal shocks (Kobzar et al. 2021; Ha et al. 2021) have 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én ion cyclotron (AIC) instability, which is 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 nonrelativistic shocks. A combined hybrid kinetic and test-particle model was 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, which are much larger than investigated previously with fully kinetic PIC simulations by Guo et al. (2018), Kobzar et al. (2021), and 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 it 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 the shock obliquity angle. Section 5 compares the results from 2D and 3D simulations. Finally, Sect. 6 summarizes our findings and presents concluding remarks.

2. 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-to-mass 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₀ 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₀, 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₀ = −[(u₀/c)×B₀], is also initialized and points in the z-direction, ${\mathit{E}}_0=E_{0z}\widehat{z}$. 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 Sect. 5.

In the following, the velocities are normalized to the ion Alfvén speed, $u_{\mathrm{A}}=B_0/\sqrt{4\pi N_{0\mathrm{i}}{m}_{\mathrm{i}}}$, calculated using the upstream magnetic field and the ion density, N_0i. The ratio between the Alfvén speed and the speed of light is set to u_A/c = 10⁻⁴, 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₀. The simulation-frame Mach number is defined as:

$$\begin{array}{c}\hfill M_{\mathrm{s},0}\equiv \frac{u_0}{c_{\mathrm{s}}}=\frac{u_0}{\sqrt{2\mathrm{\Gamma }{k}_{\mathrm{B}}{T}_0/m_{\mathrm{i}}}},\end{array}$$(1)

where c_s is the sound velocity and k_B the Boltzmann constant. The shock Mach number in the shock rest frame is:

$$\begin{array}{c}\hfill M_{\mathrm{s}}=\frac{u_{\mathrm{sh}}}{c_{\mathrm{s}}}=M_{\mathrm{s},0}(1+\frac{1}{r(M_{\mathrm{s}})-1}),\end{array}$$(2)

and the density jumps across the shock, r(M_s), in the limit of weakly magnetized flows, is equal to:

$$\begin{array}{c}\hfill r(M_{\mathrm{s}})=\frac{\mathrm{\Gamma }+1}{\mathrm{\Gamma }-1+2/M_{\mathrm{s}}^2}·\end{array}$$(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:

$$\begin{array}{c}\hfill \beta ={\beta }_{\mathrm{e}}+{\beta }_{\mathrm{i}}=\frac{8\pi (N_{\mathrm{e}}+N_{\mathrm{i}})k_{\mathrm{B}}{T}_0}{B_0^2},\end{array}$$(4)

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:

$$\begin{array}{c}\hfill \beta =\frac{2M_{\mathrm{A}}^2}{\mathrm{\Gamma }{M}_{\mathrm{s}}^2},\end{array}$$(5)

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_{\mathrm{th},\mathrm{i}}/u_{\mathrm{A}}=\sqrt{{\beta }_{\mathrm{i}}/2}$, which sets up also the normalized ion sound velocity, $c_{\mathrm{s}}/u_{\mathrm{A}}=\sqrt{2\mathrm{\Gamma }{\beta }_{\mathrm{i}}/2}$, and the change of the normalized inflow velocity, that is, the Alfvénic Mach number of the upstream flow, M_A, 0 = u₀/u_A.

Our unit of space is the ion skin depth:

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{i}}=\frac{c}{\sqrt{4\pi e^2{N}_{\mathrm{i}}/m_{\mathrm{i}}}},\end{array}$$(6)

or the ion gyroradius

$$\begin{array}{c}\hfill r_{\mathrm{g}}\equiv \frac{u_0}{{\mathrm{\Omega }}_{\mathrm{ci}}}=\frac{u_0}{u_{\mathrm{sh}}}{M}_{\mathrm{s}}{\lambda }_{\mathrm{i}}\sqrt{\beta \mathrm{\Gamma }/2},\end{array}$$(7)

where e is the electric charge. The time unit is the ion gyrofrequency, Ω_ci = eB₀/m_ic. The time step $\delta t={10}^{-3}{\mathrm{\Omega }}_{\mathrm{ci}}^{-1}$. 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.

The list of simulation runs presented in this paper is summarized in Table 1. Runs A1–A6 explore M_s = 3 shocks in the plasma beta range β = 3 − 100 and for θ_Bn = 75°. The Mach number M_s = 3 represents the typical merger shock conditions and was used in earlier PIC simulations studies (Kobzar et al. 2021). The choice of θ_Bn = 75° enables us to directly compare the results for ion physics of hybrid-kinetic and fully kinetic modelings. Runs B1–B6 study the quasi-perpendicular shock obliquity angle range, θ_Bn = 55 ° −89°, for M_s = 3 and β = 10. Runs C1–C6 scan the range of sonic Mach numbers, M_s = 2.0 − 5.5, for β = 20 and θ_Bn = 75°. Runs D1–D2 are performed to compare the θ_Bn-dependence for M_s = 3 and β = 20. Finally, runs E1–3D/E2–3D are the 3D simulations for M_s = 3, θ_Bn = 75° and β = 5 and 10.

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), and 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.8 r_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 }\sim (20-32){\mathrm{\Omega }}_{\mathrm{ci}}^{-1}$, with only two long-term runs with $t_{\max }{\mathrm{\Omega }}_{\mathrm{ci}}^{-1}=50$ for M_s = 3, β = 50, and θ_Bn = 73° shock in Ha et al. (2021) and $t_{\max }{\mathrm{\Omega }}_{\mathrm{ci}}^{-1}=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.9 r_g to even 47 r_g (see seventh and eighth columns in Table 1). Some of the simulations are also evolved to strongly nonlinear evolution stages, $t_{\max }\gg 50{\mathrm{\Omega }}_{\mathrm{ci}}^{-1}$ (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 θ_Bn-dependence over such a wide range of angles have not been undertaken before.

3. 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 quasi-perpendicular 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 Sect. 3.1. Section 3.2 discusses the linear analysis of the magnetic wave modes, Sect. 3.3 treats density perturbations and shock ripples, and Sect. 3.4 describes the nonlinear evolution of the turbulence in the shock.

3.1. 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 Ω_cit = 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.5 r_g is shown in Fig. 1 with a vertical blue dotted line. The corresponding ion velocity phase-space is presented in Fig. 2.

Fig. 1. Structure of the reference shock with Ms = 3, β = 20, and θBn = 75° (run A4) at time Ωcit = 30. Displayed are the y-averaged profiles of (a) the ion density normalized to the far upstream density N0i (black line), Bx, By and Bz magnetic field components normalized to their upstream values (blue, red, and green lines, respectively), (b) profiles of Tx, Ty, and Tz 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 T0, (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 rg.

Fig. 2. Ion velocity phase-space distributions, x − vi (i = x,  y,  z), averaged in the y-direction for Ms = 3, β = 20, and θBn = 75° shock (run A4) at time Ωcit = 30. The vertical dotted line marks the approximate shock location.

The normalized ion density profile, shown with black line in Fig. 1a, exhibits a characteristic pattern of overshoot-undershoot oscillations past the shock front. The density compression at the double-peaked overshoot (at x ≈ 15.2 r_g and x ≈ 14.5 r_g) reaches N_i/N₀ ≈ 3.5, decreases to N_i/N₀ ≈ 2 in the undershoot at x ≈ 13.8 r_g, and increases again to N_i/N₀ ≈ 3.1 in the second overshoot at x ≈ 13.0 r_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 (Figs. 2a and c). The reflected ions gyrate in the upstream magnetic field and acquire energy as they drift along the motional electric field, E₀. 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 1 r_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₀. 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 Fig. 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 Fig. 1d¹. 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 Fig. 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 Ω_cit = 30 are presented in Fig. 3b. The field fluctuations are defined with respect to the flux-frozen magnetic field, as proposed in Guo et al. (2017):

$$\begin{array}{cc}& \delta B_x=B_x-B_{0}cos{\theta }_{B\mathrm{n}},\hfill \end{array}$$(8)

$$\begin{array}{cc}& \delta B_y=B_y-B_{0}sin{\theta }_{B\mathrm{n}}\frac{N_{\mathrm{i}}(x)}{N_{0\mathrm{i}}},\hfill \end{array}$$(9)

$$\begin{array}{cc}& \delta B_z=B_z,\hfill \end{array}$$(10)

Fig. 3. Structure of the reference shock with Ms = 3, β = 20, and θBn = 75° (run A4) at different times Ωcit = 10 (panels a*), Ωcit = 30 (panels b*), and Ωcit = 55 (panels c*). From top to bottom, shown are the maps of the normalized ion density, Ni/N0i (panels *1), and the normalized magnetic field fluctuations (see Eqs. (8)–(10)), δBx (panels *2), δBy (panels *3), and δBz (panels *4). The scaling of all distributions is logarithmic and for the magnetic fields it is sign-preserving, and, e.g., for δBz it is: sgn(δBz)⋅{2 + log[max(|δBz|/B0, 10−2)]}. The level of “0” on the color scale hence corresponds to |δB|/B0 ≤ 10−2 (see, e.g., Kobzar et al. 2021).

where N_i(x) is the y-averaged ion density. They are normalized to the upstream magnetic field, B₀. 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 Fig. 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 Ω_cit = 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 Ω_cit ≈ 6. Their amplitude and wavelength progressively increase over time. By the time Ω_cit = 30, the shock has transitioned into the strongly nonlinear stage. The amplitude of the shock ripples reaches δN_i/N₀ ≈ 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 6 r_g in width from the shock front (compare Figs. 1 and 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).

3.2. 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–a4 show magnetic fluctuations at time Ω_cit = 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 Alfvén 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 antiparallel 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 Sect. 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 Ω_cit = 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 coexist 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 with small k_xλ_i < 0.7 components. These waves are well visible in the δB_y map in Fig. 3a3, but not in the spectrum in Fig. 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).

Fig. 4. Fourier power spectra of the magnetic waves in the region from the shock ramp to behind the second overshoot, 1 ≲ x/rg ≲ 5.5 at Ωcit = 10 for run A4 (compare Fig. 3a). We note the different dynamic scales in each panel.

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 Ω_cit = 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 k_∥ ≈ k_y in our simulation, and k_⊥ ≈ k_x. Note that primed units of the wave vectors and the growth rates are expressed in terms of the local parameters so that ${\lambda }_{\mathrm{i}}^{\prime }={\lambda }_{\mathrm{i}}/\sqrt{r}$ and ${\mathrm{\Omega }}_{\mathrm{ci}}^{\prime }={\mathrm{\Omega }}_{\mathrm{ci}}r$.

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.

The area of the strongest growth at $k_{\Vert }{\lambda }_{\mathrm{i}}^{\prime }\approx 1$ and k_⊥ ≈ 0 in panel a for the shock ramp pertains to the AIC mode whose wavelength in the simulation units is λ ≈ 5.1λ_i. The growth rate and the oscillation frequency of this mode are very high, $\gamma /{\mathrm{\Omega }}_{\mathrm{ci}}^{\prime }\approx 1.2$ and $\omega /{\mathrm{\Omega }}_{\mathrm{ci}}^{\prime }\approx 1.5$ so that the phase velocity of the AIC wave is v_ph = ω/k ≈ 1.8 u_A (in units of the upstream Alfvén velocity). The oblique mirror mode at $(k_{\Vert },k_{\perp }){\lambda }_{\mathrm{i}}^{\prime }\approx (1.35,0.45)$ has wavelength of λ ≈ 3.6λ_i and $\gamma /{\mathrm{\Omega }}_{\mathrm{ci}}^{\prime }\approx 0.67$. The AIC mode in the undershoot-second overshoot region in panel b peaks at much smaller linear growth of $\gamma /{\mathrm{\Omega }}_{\mathrm{ci}}^{\prime }\approx 0.4$ and much smaller $k_{\Vert }{\lambda }_{\mathrm{i}}^{\prime }\approx 0.4$, which corresponds to λ ≈ 9.6λ_i. With the real frequency of $\omega /{\mathrm{\Omega }}_{\mathrm{ci}}^{\prime }\approx 0.8$, the phase velocity of this mode is v_ph ≈ 3.3 u_A. A weakly growing mirror mode with $\gamma /{\mathrm{\Omega }}_{\mathrm{ci}}^{\prime }\approx 0.2$ also has a longer wavelength of λ ≈ 7.6λ_i, corresponding to $(k_{\Vert },k_{\perp }){\lambda }_{\mathrm{i}}^{\prime }\approx (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.

3.3. 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 Fig. 3a and also in Fig. 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 Ω_cit ≈ 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 Sect. 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 ion-scale 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.

3.4. 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 Ω_cit ≈ 6. One can see in Fig. 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 Ω_cit = 10).

Fig. 6. Temporal evolution of the dominant ky wavevector component of the magnetic field and ion density fluctuations in the shock front for the reference simulation with Ms = 3, β = 20, and θBn = 75° (run A4). The error bars reflect resolution in the Fourier space (see text).

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 Ω_cit ≈ 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 Ω_cit ≈ 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 larger-scale ones (see Figs. 3b and c). However, the structures do not evolve toward even longer wavelengths by the end of our reference run, A4 at Ω_cit = 56. On the other hand, weak downstream turbulence beyond the second overshoot remains at a quasi-steady state at times Ω_cit ≳ 30. Note, that we observe similar behavior for M_s = 3 shock propagating in β = 5 plasma (run A2), which we follow up to Ω_cit = 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.

4. Parameter dependence

In this section we discuss the dependence of the results obtained for the reference run A4 with M_s = 3, β = 20, and θ_Bn = 75° on the shock parameters: plasma beta in the range β = 5 − 100 in Sect. 4.1, the shock Mach number in the range M_s = 2 − 5.5 in Sect. 4.2, and the range of the upstream magnetic field inclination angle, θ_Bn = 55 ° −90°, in Sect. 4.3.

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):

$$\begin{array}{c}\hfill \frac{T_{\perp }}{T_{\Vert }}-1\simeq \frac{1.6}{{\beta }_{\Vert }^{0.72}},\end{array}$$(11)

where β_∥ = 8πN_ik_BT_∥/B² is the proton plasma beta at a given location, calculated using the proton temperature in the parallel direction. Figure 7 shows average profiles of the ion temperature anisotropy (solid lines) for runs A2–A6 for shocks with different plasma β and M_s = 3, θ_Bn = 75° (panel a), runs C1–C6 for different shock sonic Mach numbers and β = 20, θ_Bn = 75° (panel b), and runs A3 and B1–B7 for a range of obliquity angles θ_Bn = 55 ° −89° for shocks with M_s = 3 and β = 10 (panel c). According to Eq. (11), marginal stability limits are shown in these plots with dotted lines. The profiles are presented at an early stage of the shock evolution at time Ω_cit = 8.

Fig. 7. Profiles (y-averaged) of the ion temperature anisotropy (solid lines); the marginal stability condition of Eq. (11) (dotted lines) at time Ωcit = 8 for shocks with different plasma β and Ms = 3 and θBn = 75° (runs A2–A6, panel a); different sonic Mach numbers and β = 20 and θBn = 75° (runs C1–C6, panel b); and a range of obliquity angles for Ms = 3 and β = 10 (runs A3 and B1–B7, panel c). The x coordinate was normalized to the ion gyroradius rg(Ms, β) of each model.

4.1. Dependence on plasma β

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 θ_Bn = 75°. One can see in Fig. 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 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 (Figs. 8b and 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).

Fig. 8. Structure of shocks with Ms = 3 and θBn = 75° at Ωcit = 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 Bz 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, xsh, with respect to which distance is calculated in Figs. 9 and 10.

The relationship between the plasma beta and the location of the peak magnetic energy is most evident in Figs. 9 and 10 that present profiles of the magnetic field energy density at times Ω_cit = 26 and Ω_cit = 55, respectively. The maximum magnetic energy density is observed approximately 1 r_g behind the shock with β = 100, about 2 r_g for β = 50 and 20, to around 3 r_g for β = 5. Lower-β shocks develop weaker turbulence, Fig. 8. By the early nonlinear simulation phase at Ω_cit = 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, Fig. 9. The amplitude levels do not change in the strongly nonlinear stage; compare Fig. 10a. From the early stage of the evolution, the magnetic turbulence is dominated by the B_z field component, Figs. 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, 2018).

Fig. 9. Profiles (y-averaged) at time Ωcit = 26 of the normalized total energy density of the magnetic field fluctuations across the shocks with Ms = 3 and θBn = 75° and plasma beta in the range β = 3 − 100, runs A1–A6. The x coordinate was normalized to the ion gyroradius rg and measured relative to the shock location xsh.

Fig. 10. Profiles (y-averaged) at time Ωcit = 55 of the normalized energy density of the magnetic field fluctuations (total and in Bx, By, and Bz components) across shocks with Ms = 3 and θBn = 75° and plasma β = 5, 20, and 50, and runs A2, A4, and A5.

Runs A2 and A5, have been followed to late times Ω_cit ≥ 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.9 r_g.

4.2. 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 θ_Bn = 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 Fig. 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 1 r_g downstream from the shock front, Figs. 11c and 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.5 r_g for M_s = 2.9 shock, through approximately 3 r_g, 4 r_g, and 6 r_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).

Fig. 11. Structure of shock with β = 20 and θBn = 75° at time Ωcit = 20 for different values of the sonic Mach number, Ms = 2.00, 2.25, 5.50, runs C1, C2, and C6 in panels a–c, respectively (compare Fig. 10).

Fig. 12. Profiles (y-averaged) at time Ωcit = 20 of the normalized total energy density of the magnetic field fluctuations across shocks with β = 20 and θBn = 75° and Mach numbers in a range Ms = 2 − 5.5 and runs C1–C6 (compare Fig. 9).

One can see in the upper panels of Fig. 11 and in Fig. 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 is smooth and the ion-scale turbulence is built only in the region of the second overshoot (∼5 r_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_{\text{AIC}}^{*}\approx 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 Eq. (11). This coincided with an enhanced ion reflection efficiency from the shock at shocks with M_s ≳ 2.3. At shocks with $M_{\text{s}}<M_{\text{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 ∼5 r_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.

4.3. Dependence on the obliquity angle θ_Bn

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 quasi-perpendicular, θ_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.

Fig. 13. Structure of shocks with β = 20 and Ms = 3 at time Ωcit = 56 for different values of the obliquity angle, θBn = 55°, 75°, and 90°, and runs D1, A4, and D2 in panels a–c, respectively (compare Fig. 10).

The variations in the magnitude of the temperature anisotropy at the shock front with the field obliquity seen in Fig. 7c are related to the geometry of the ion reflection from the shock and their gyration in the upstream magnetic field. Upon 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).

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 Fig. 13 and in the magnetic energy density profiles in Fig. 14. However, despite a temperature anisotropy difference of approximately 3.5 between shocks with θ_Bn = 55° and θ_Bn = 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.

Fig. 14. Profiles (y-averaged) at time Ωcit = 30 of the normalized total energy density of the magnetic field fluctuations across the shocks with β = 10 and Ms = 3 and obliquity angles in the range θBn = 55 ° −89°, and runs B1–B4, A3, and B5–B7 (compare Fig. 9).

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 Ω_cit = 56 for β = 20 shocks, Fig. 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 Sect. 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 Alfvénic 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 Sect. 3.4 that the wavelength of the ripple mode saturates past the simulation time Ω_cit ≈ 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 antiparallel to the mean magnetic field (see Sect. 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.

5. 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 Ω_cit = 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.

Fig. 15. Structure of the Ms = 3, β = 5, θBn = 75° shock at Ωcit = 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, xsh, is marked with the white dotted line.

The cut through the box in the direction transverse to the plasma flow at x = x_sh − 0.12 r_g presented in Fig. 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.

Fig. 16. Maps of density and Bx magnetic field fluctuations in the y − z plane at x = xsh − 0.12 rg (see Fig. 15) at Ωcit = 36 for the E1–3D run.

Finally, Figs. 17a and b compare the profiles of the normalized energy density of the magnetic fluctuations between 2D (solid lines) and 3D (dash-dotted lines) simulations at Ω_cit = 24 and Ω_cit = 36, respectively. At an earlier stage, Fig. 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 Ω_cit = 36 in this case, Fig. 17b, because the wave growth rates are lower for smaller β = 5, compared to β = 10 (see Sect. 4.1).

Fig. 17. Profiles (y-averaged) of the normalized energy density of the magnetic field fluctuations across shocks in 2D (solid lines) and 3D (dash-dotted lines) simulations, (a) at time Ωcit = 24 for shocks with β = 5 and 10, Ms = 3, θBn = 75°, and runs A2, A3, E1–3D, and E2–3D, and (b) at Ωcit = 36 for runs A2 and E1–3D with β = 5.

The results presented in this section indicate that the steady-state 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 Sect. 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.

6. Summary and concluding remarks

For this study we employed hybrid-kinetic numerical simulations to investigate the ion-scale kinetic instabilities at merger shocks of galaxy clusters. Our results can be summarized as follows:

• 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.

• The nonlinear evolution of AIC and mirror modes in addition to 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.9 r_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 they form weak turbulence in the shock downstream. This is consistent with the notion of the critical Mach number $M_{\text{AIC}}^{*}\approx 2.3$ to trigger multi-scale waves in the shock front, which has been proposed in Ha et al. (2021).

• For the first time, we also investigate high-β quasi-perpendicular shocks in a wide range of the magnetic field obliquity. Despite significant variations in temperature anisotropy at the shock, we only observed 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 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 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 Sect. 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 of an 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 incorporation of an additional energetic electron component and the exploration of electron acceleration within turbulent shock structures. Such investigations are currently ongoing.

Acknowledgments

We thank the anonymous referee for useful comments. The authors acknowledge the support of Narodowe Centrum Nauki through research projects no. 2019/33/B/ST9/02569 (S.B. and J.N.) and no. 2016/22/E/ST9/00061 (O.K.). We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2022/015967. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #520 (Energy Partition across Collisionless Shocks).

References

Appendix A: Mass ratio dependence

The hybrid simulation code adopted in this study includes the finite electron inertia effect. The model thus correctly considers the dispersive effect on the linear dispersion relation associated with the finite electron mass, which appears at the electron skin depth scale. This indicates that the finite mass ratio should not affect the results as long as the electron skin depth is unresolved. In this case, we can rather consider the finite inertia effect as a numerical stabilization trick to suppress high-frequency whistler noise (Amano et al. 2014).

As an example, Fig. A.1 compares the snapshots of the temperature anisotropy at Ω_cit = 10 for β = 5, M_s = 3.0, θ_Bn = 75° with two different mass ratios m_i/m_e = 100 and 1836. As we can see, the results are almost identical. This is natural since, in both cases, the spatial resolution 0.25λ_i is larger than the electron skin depth (0.1λ_i and 0.02λ_i for m_i/m_e = 100 and 1836, respectively).

Fig. A.1. Temperature anisotropy profiles across the shock transition at Ωcit = 10, showing a dependence on the ion-to-electron mass ratio, mi/me, for the test case of β = 5, Ms = 3.0, and θBn = 75°.

Note that, in general, attempting to resolve the electron skin depth with the hybrid code is not advisable. In fact, the electron kinetic effect, such as the cyclotron damping ignored in the model, starts to play a role at this scale. Therefore, one has to be careful in interpreting the physics at the electron scale in the hybrid model.

Appendix B: Dependence on the number of particles per cell

In Fig. B.1, we explore the impact of varying the number of computational particles per cell, N_ppc, with values of N_ppc = 32 (red), 64 (blue), 128 (green), and 256 (purple). Noticeable distinctions are not particularly evident within the shock front, where the temperature anisotropy is most pronounced, especially for N_ppc > 64. However, these differences become noteworthy past the shock front, where strong magnetic turbulence develops. For N_ppc > 128 the results are convergent. We use N_ppc = 128 in our simulations.

Fig. B.1. Temperature anisotropy profiles across the shock transition at Ωcit = 10, showing a dependence on the number of particles per cell, Nppc, for the test case of β = 10, Ms = 3.0, and θBn = 75°.

Appendix C: Periodic box simulation of the wave turbulence evolution

In Figure 4 in Section 3.2, we present the results of the linear analysis using the local macroscopic properties of the plasma at two regions of the shock transition – the shock ramp and the region between the undershoot and the second overshoot. To confirm the compatibility of early-stage instabilities with linear calculations and explore the nonlinear evolution of these instabilities, we conducted 2D simulations with periodic boundary conditions for the values at the shock ramp. Hence, we assume the temperature anisotropy T_⊥/T_∥ = 7, parallel ion beta β_∥ = 9.33, and the electron plasma beta β_e = 10. Note that while these parameters are measured in the simulation run A4 at time Ω_cit = 10, they show no significant differences from parameters measured at the shock ramp at earlier times. Therefore, they should be considered representative of the early-stage shock evolution at this location.

The simulation domain is a square box of size $L_x=L_z=60{\lambda }_i^{\prime }$ in the x − y plane, and the mean magnetic field B₀ is along the x-direction (note that in shock simulations, the compressed downstream magnetic field is essentially aligned with the y-direction). Here, we denote the unit length scale ${\lambda }_i^{\prime }$ to emphasize that wavelengths and wave vectors are given in terms of the local values measured in the shock simulation A4. All other simulation parameters are the same as in the shock simulations.

For the given simulation setup, the early evolution of the instabilities should correspond to the linear theory results presented in Figure 4a. We expect that the AIC waves are formed in B_y and B_z magnetic field components, and the oblique compressive mirror mode is seen primarily in the B_x component. Figure C.1 shows the time evolution of the magnetic field fluctuations and temperature profiles. Field turbulence rapidly grows and saturates around Ω_cit ≈ 10, coinciding with the near disappearance 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 Ω_cit ≲ 10.

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.

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 ${\mathrm{\Omega }}_{ci}^{\prime }t=2.5-5$, the measured growth rates are ${\gamma }_{\text{AIC}}\approx 0.8{\mathrm{\Omega }}_{ci}^{\prime }$ for δB_y, δB_z components and ${\gamma }_m\approx 0.4{\mathrm{\Omega }}_{ci}^{\prime }$ 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 ${\gamma }_{\text{AIC}}\approx 1.2{\Omega }_{ci}^{\prime }$ and ${\gamma }_m\approx 0.6{\mathrm{\Omega }}_{ci}^{\prime }$. 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.

Figure C.2 shows distributions of δB_x, δB_y, and δB_z magnetic field fluctuations at time ${\mathrm{\Omega }}_{ci}^{\prime }t=4$. Figure C.3 displays corresponding Fourier power spectra averaged in the time interval $3\le {\mathrm{\Omega }}_{ci}^{\prime }t\le 5$. The Fourier power in B_y and B_z waves is maximum at $k_{\parallel }{\lambda }_i^{\prime }\approx 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_{\perp }{\lambda }_i^{\prime }\lesssim 0.25$ and similar $k_{\parallel }{\lambda }_i^{\prime }\approx 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_{\parallel }{\lambda }_i^{\prime }\approx 0.6-1.4$ and $k_{\perp }{\lambda }_i^{\prime }\approx 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.

Fig. C.2. Maps of the normalized magnetic field fluctuations for the periodic box simulation at time ${\mathrm{\Omega }}_{ci}^{\prime }t=4$. From left to right, (a) δBx = Bx − B0, (b) δBy, and (c) δBz are shown.

Fig. C.3. Fourier power spectra of magnetic field fluctuations in the periodic box simulation averaged in the time interval $3\le {\mathrm{\Omega }}_{ci}^{\prime }t\le 5$. We note that the dynamic range is different in each panel.

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., in prep.). 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 Figure D.1. The shock structures are compared at slightly different times due to the distinct nature of the two simulation codes. Nevertheless, the results are considered to be contemporaneous. PIC simulations have been performed with 20 particles per cell per species, ion-to-electron mass ratio m_i/m_e = 100, and the electron skin depth λ_e = c/ω_pe = 15 cells, where ${\omega }_{\mathit{pe}}=\sqrt{4\pi e^2{N}_e/m_e}$ is the electron plasma frequency.

Fig. D.1. Comparison of the shock structures obtained with PIC (left panels) and hybrid (right panels) simulations for shocks with Ms = 3 and θBn = 75°. The top row shows results for β = 5 at time Ωcit = 10 for the PIC (a) and Ωcit = 13 for the hybrid (b) simulation (run A2). The bottom row displays the results for β = 20 at time Ωcit = 10 for the PIC (c) and Ωcit = 11 for the hybrid (d) simulation (run A4). In each panel, from top to bottom, normalized distributions of ion density, and δBy and δBz magnetic field fluctuations are shown. We note that the scaling of all distributions is logarithmic (see Fig. 3). Only portions of the computational boxes are shown for hybrid simulation runs A2 and A4.

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 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.

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 Ω_cit = 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, 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.

All Tables

All Figures

Fig. 1. Structure of the reference shock with Ms = 3, β = 20, and θBn = 75° (run A4) at time Ωcit = 30. Displayed are the y-averaged profiles of (a) the ion density normalized to the far upstream density N0i (black line), Bx, By and Bz magnetic field components normalized to their upstream values (blue, red, and green lines, respectively), (b) profiles of Tx, Ty, and Tz 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 T0, (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 rg.

In the text

	Fig. 2. Ion velocity phase-space distributions, x − vi (i = x,  y,  z), averaged in the y-direction for Ms = 3, β = 20, and θBn = 75° shock (run A4) at time Ωcit = 30. The vertical dotted line marks the approximate shock location.
In the text

Fig. 3. Structure of the reference shock with Ms = 3, β = 20, and θBn = 75° (run A4) at different times Ωcit = 10 (panels a*), Ωcit = 30 (panels b*), and Ωcit = 55 (panels c*). From top to bottom, shown are the maps of the normalized ion density, Ni/N0i (panels *1), and the normalized magnetic field fluctuations (see Eqs. (8)–(10)), δBx (panels *2), δBy (panels *3), and δBz (panels *4). The scaling of all distributions is logarithmic and for the magnetic fields it is sign-preserving, and, e.g., for δBz it is: sgn(δBz)⋅{2 + log[max(|δBz|/B0, 10−2)]}. The level of “0” on the color scale hence corresponds to |δB|/B0 ≤ 10−2 (see, e.g., Kobzar et al. 2021).

In the text

	Fig. 4. Fourier power spectra of the magnetic waves in the region from the shock ramp to behind the second overshoot, 1 ≲ x/rg ≲ 5.5 at Ωcit = 10 for run A4 (compare Fig. 3a). We note the different dynamic scales in each panel.
In the text

	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.
In the text

	Fig. 6. Temporal evolution of the dominant ky wavevector component of the magnetic field and ion density fluctuations in the shock front for the reference simulation with Ms = 3, β = 20, and θBn = 75° (run A4). The error bars reflect resolution in the Fourier space (see text).
In the text

Fig. 7. Profiles (y-averaged) of the ion temperature anisotropy (solid lines); the marginal stability condition of Eq. (11) (dotted lines) at time Ωcit = 8 for shocks with different plasma β and Ms = 3 and θBn = 75° (runs A2–A6, panel a); different sonic Mach numbers and β = 20 and θBn = 75° (runs C1–C6, panel b); and a range of obliquity angles for Ms = 3 and β = 10 (runs A3 and B1–B7, panel c). The x coordinate was normalized to the ion gyroradius rg(Ms, β) of each model.

In the text

Fig. 8. Structure of shocks with Ms = 3 and θBn = 75° at Ωcit = 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 Bz 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, xsh, with respect to which distance is calculated in Figs. 9 and 10.

In the text

	Fig. 9. Profiles (y-averaged) at time Ωcit = 26 of the normalized total energy density of the magnetic field fluctuations across the shocks with Ms = 3 and θBn = 75° and plasma beta in the range β = 3 − 100, runs A1–A6. The x coordinate was normalized to the ion gyroradius rg and measured relative to the shock location xsh.
In the text

	Fig. 10. Profiles (y-averaged) at time Ωcit = 55 of the normalized energy density of the magnetic field fluctuations (total and in Bx, By, and Bz components) across shocks with Ms = 3 and θBn = 75° and plasma β = 5, 20, and 50, and runs A2, A4, and A5.
In the text

	Fig. 11. Structure of shock with β = 20 and θBn = 75° at time Ωcit = 20 for different values of the sonic Mach number, Ms = 2.00, 2.25, 5.50, runs C1, C2, and C6 in panels a–c, respectively (compare Fig. 10).
In the text

	Fig. 12. Profiles (y-averaged) at time Ωcit = 20 of the normalized total energy density of the magnetic field fluctuations across shocks with β = 20 and θBn = 75° and Mach numbers in a range Ms = 2 − 5.5 and runs C1–C6 (compare Fig. 9).
In the text

	Fig. 13. Structure of shocks with β = 20 and Ms = 3 at time Ωcit = 56 for different values of the obliquity angle, θBn = 55°, 75°, and 90°, and runs D1, A4, and D2 in panels a–c, respectively (compare Fig. 10).
In the text

	Fig. 14. Profiles (y-averaged) at time Ωcit = 30 of the normalized total energy density of the magnetic field fluctuations across the shocks with β = 10 and Ms = 3 and obliquity angles in the range θBn = 55 ° −89°, and runs B1–B4, A3, and B5–B7 (compare Fig. 9).
In the text

	Fig. 15. Structure of the Ms = 3, β = 5, θBn = 75° shock at Ωcit = 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, xsh, is marked with the white dotted line.
In the text

	Fig. 16. Maps of density and Bx magnetic field fluctuations in the y − z plane at x = xsh − 0.12 rg (see Fig. 15) at Ωcit = 36 for the E1–3D run.
In the text

	Fig. 17. Profiles (y-averaged) of the normalized energy density of the magnetic field fluctuations across shocks in 2D (solid lines) and 3D (dash-dotted lines) simulations, (a) at time Ωcit = 24 for shocks with β = 5 and 10, Ms = 3, θBn = 75°, and runs A2, A3, E1–3D, and E2–3D, and (b) at Ωcit = 36 for runs A2 and E1–3D with β = 5.
In the text

	Fig. A.1. Temperature anisotropy profiles across the shock transition at Ωcit = 10, showing a dependence on the ion-to-electron mass ratio, mi/me, for the test case of β = 5, Ms = 3.0, and θBn = 75°.
In the text

	Fig. B.1. Temperature anisotropy profiles across the shock transition at Ωcit = 10, showing a dependence on the number of particles per cell, Nppc, for the test case of β = 10, Ms = 3.0, and θBn = 75°.
In the text

	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.
In the text

	Fig. C.2. Maps of the normalized magnetic field fluctuations for the periodic box simulation at time ${\mathrm{\Omega }}_{ci}^{\prime }t=4$. From left to right, (a) δBx = Bx − B0, (b) δBy, and (c) δBz are shown.
In the text

	Fig. C.3. Fourier power spectra of magnetic field fluctuations in the periodic box simulation averaged in the time interval $3\le {\mathrm{\Omega }}_{ci}^{\prime }t\le 5$. We note that the dynamic range is different in each panel.
In the text

Fig. D.1. Comparison of the shock structures obtained with PIC (left panels) and hybrid (right panels) simulations for shocks with Ms = 3 and θBn = 75°. The top row shows results for β = 5 at time Ωcit = 10 for the PIC (a) and Ωcit = 13 for the hybrid (b) simulation (run A2). The bottom row displays the results for β = 20 at time Ωcit = 10 for the PIC (c) and Ωcit = 11 for the hybrid (d) simulation (run A4). In each panel, from top to bottom, normalized distributions of ion density, and δBy and δBz magnetic field fluctuations are shown. We note that the scaling of all distributions is logarithmic (see Fig. 3). Only portions of the computational boxes are shown for hybrid simulation runs A2 and A4.

In the text