Microinstabilities in the Transition Region of Weak Quasi-Perpendicular Intracluster Shocks

Microinstabilities play important roles in both entropy generation and particle acceleration in collisionless shocks. Recent studies have suggested that in the transition zone of quasi-perpendicular ($Q_{\perp}$) shocks in the high-beta ($\beta=P_{\rm gas}/P_{\rm B}$) intracluster medium (ICM), the ion temperature anisotropy due to the reflected-gyrating ions could trigger the Alfv\'en ion cyclotron (AIC) instability and the ion-mirror instability, while the electron temperature anisotropy induced by magnetic field compression could excite the whistler instability and the electron-mirror instability. Adopting the numerical estimates for ion and electron temperature anisotropies found in particle-in-cell (PIC) simulations of $Q_{\perp}$-shocks with sonic Mach numbers, $M_{\rm s}=2-3$, we carry out a linear stability analysis for these microinstabilities. The kinetic properties of the microinstabilities and the ensuing plasma waves on both ion and electron scales are described for wide ranges of parameters, including the dependence on $\beta$ and the ion-to-electron mass ratio. In addition, the nonlinear evolution of induced plasma waves are examined by performing 2D PIC simulations with periodic boundary conditions. We find that for $\beta\approx 20-100$, the AIC instability could induce ion-scale waves and generate shock surface ripples in supercritical shocks above the AIC critical Mach number, $M_{\rm AIC}^{*} \approx 2.3$. Also electron-scale waves are generated primarily by the whistler instability in these high-$\beta$ shocks. The resulting multi-scale waves from electron to ion scales are thought to be essential in electron injection to the diffusive shock acceleration mechanism in $Q_{\perp}$-shocks in the ICM.


INTRODUCTION
Major mergers of galaxy clusters are known to drive weak shocks with sonic Mach numbers, M s 3, in the hot intracluster medium (ICM) of high β (e.g., Ryu et al. 2003;Skillman et al. 2008;Vazza et al. 2009;Hong et al. 2014;Ha et al. 2018). Here, the plasma beta, β = P gas /P B , is the ratio of the gas pressure to the magnetic pressure. The radiative signatures of such shocks have been detected in X-ray and radio observations (e.g., Brüggen et al. 2012;Brunetti & Jones 2014). In the case of so-called radio relics, the radio emission has been interpreted as the synchrotron radiation from relativis-tic electrons accelerated via diffusive shock acceleration (DSA) in merger-driven shocks (see van Weeren et al. 2019, for a review).
To explain the origin of radio relics, this DSA model requires an electron preacceleration mechanism, because postshock thermal electrons do not have momenta large enough to participate in the standard DSA process, in which cosmic ray (CR) electrons diffuse across the shock (Kang et al. 2012). Since the width of the shock transition layer is comparable to the gyro-radius of postshock thermal ions, thermal electrons need to be energized to the so-called injection momentum, p inj ∼ a few × p i,th . Here, p i,th = √ 2m i k B T i2 is the ion thermal momentum in the postshock gas of temperature T i2 , m i is the ion mass, and k B is the Boltzmann constant. For shocks in the solar wind, the electron injection is observed preferentially at the quasi-perpendicular (Q ⊥ ) configuration Table 1. Linear Properties of the Instabilities driven by Perpendicular Temperature Anisotropies instability AIC whistler ion-mirror electron-mirror free energy source T i⊥ > T i T e⊥ > T e T i⊥ > T i T e⊥ > T e propagation angle with γm a parallel parallel oblique oblique wavenumber ck/ωpi ≤ 1 ck/ωpe ≤ 1 ck/ωpi ≤ 1 ck/ωpe ≤ 1 wave frequency 0 < ωr < Ωci Ωci < ωr < Ωce ωr = 0 ωr = 0 wave polarization LHCP b RHCP b Non-propagating Non-propagating a γm is the maximum growth rate.
with θ Bn 45 • , where θ Bn is the shock obliquity angle between the shock normal and the upstream magnetic field direction (e.g., Gosling et al. 1989;Oka et al. 2006;Burgess 2006). The electron preacceleration has been a key outstanding problem in understanding the production of CR electrons in weak ICM shocks. Previous studies have shown that, in low-M s , high-β, Q ⊥ -shocks, thermal electrons could be preaccelerated primarily through the Fermilike acceleration in the shock foot (Matsukiyo et al. 2011;Guo et al. 2014a,b;Kang et al. 2019) and the stochastic shock drift acceleration (SSDA) in the shock transition region (Katou & Amano 2019;Niemiec et al. 2019;. Although it has been shown that the electron preaccleration would be be enhanced by preexisting strong magnetic fluctuations in the low-β (∼ 1) regime (e.g. Guo & Giacalone 2015;Trotta et al. 2020), the effect of such turbulence on highβ ICM shocks has yet to be investigated and will not be considered here.
Both the Fermi-like acceleration and SSDA mechanisms rely on the various microinstabilities triggered by the ion and electron temperature anisotropies in the shock structure (Gary 1993). If T e > T e⊥ , for example, the electron firehose instability (EFI) can grow with the following two branches: the nonresonant, parallelpropagating mode with left-hand circular polarization, and the resonant, non-propagating, oblique mode (Gary & Nishimura 2003). Hereafter, the subscripts and ⊥ denote the parallel and perpendicular directions to the background magnetic field, B 0 , respectively. Under the condition of T e⊥ > T e , by contrast, the whistler instability and the electron-mirror (e-mirror) instability can be triggered (Scharer & Trivelpiece 1967;Gary 1992;Hellinger &Štverák 2018). The most unstable whistler mode propagates in the direction parallel to B 0 with right-hand circular polarization, while the e-mirror mode is non-propagating and has the maximum growth rate at the wavevector direction oblique with respect to B 0 . In the case of T i⊥ > T i , the Alfvén ion cyclotron instability (AIC, or the proton cyclotron instability) and the ion-mirror (i-mirror) instability may become unstable (Winske & Quest 1988;Gary 1993;Gary et al. 1997;Burgess 2006). The fastest-growing mode of the AIC instability propagates in the direction parallel to B 0 with left-hand circular polarization, while the i-mirror mode is non-propagating and has the maximum growth rate at the wavevector direction oblique with respect to B 0 . Table 1 summarizes these linear properties of the instabilities driven by perpendicular temperature anisotropies, which are relevant for the present study.
In the foot of Q ⊥ -shocks, the shock-reflected electrons backstream mainly along the upstream magnetic field and induce an electron parallel anisotropy (T e > T e⊥ ), which could trigger the EFI and facilitate the Fermi-like preacceleration (Guo et al. 2014b;Kang et al. 2019;Kim et al. 2020). In the transition region behind the shock ramp, on the other hand, the AIC and i-mirror instabilities can be triggered by the ion perpendicular anisotropy (T i⊥ > T i ) mainly due to the shock-reflected ions advected downstream, while the whistler and e-mirror instabilities can be excited by the electron perpendicular anisotropy (T e⊥ > T e ) mainly due to magnetic field compression at the shock ramp (Guo et al. 2017;Katou & Amano 2019). Such multi-scale waves from electron to ion scales are essential in the electron preacceleration via the SSDA (Matsukiyo & Matsumoto 2015;Niemiec et al. 2019;. Using two-dimensional (2D) particle-in-cell (PIC) simulations for β ≈ 20 − 100, Q ⊥ -shocks, Kang et al. (2019) showed that the Fermi-like preacceleration involving multiple cycles of shock drift acceleration (SDA) in the shock foot could be effective only in supercritical shocks with M s greater than the EFI critical Mach number, M * ef ≈ 2.3. However, they argued that the electron preacceleration may not proceed all the way to p inj , be-cause the EFI-driven waves are limited to electron scales. Niemiec et al. (2019), on the other hand, performed a PIC simulation for M s = 3 shock with β = 5 in a 2D domain large enough to include ion-scale fluctuations, and suggested that electrons could be energized beyond p inj via the SSDA due to stochastic pitch-angle scattering off the multi-scale waves excited in the shock transition zone.
Furthermore,  found that in β ≈ 1 plasmas, the AIC instability is triggered and the ensuing electron preacceleration operates only in supercritical shocks with the Alfvénic Mach number greater than the critical Mach number, M * AIC ≈ 3.5. In a separate paper (Ha et al. 2021, we report a similar study of β ≈ 20 − 100 shocks, which is design to explore through 2D PIC simulations how the multi-scale waves excited mainly by the AIC and whistler instabilities in the shock transition can assist the electron injection to DSA in ICM shocks. In this paper, adopting the numerical estimates for the temperature anisotropies in the transition region of the simulated shocks of HKRK2021 1 , we first perform a linear stability analysis for microinstabilities for wide ranges of parameters such as M s = 2 − 3, β = 1 − 100, and the ion-to-electron mass ratio, m i /m e = 50 − 1836. This approach allows us to identify the most dominant modes of possible microinstabilities and to evaluate their linear properties for the set of realistic parameters pertaining to ICM shocks. Hence, this kind of analyses on kinetic scales can provide crucial insights for theoretical modelings and/or larger scale simulations of particle acceleration at weak high-β shocks. However, one of the limitations of PIC simulations is that it can follow kinetic plasma processes mainly at the low end of temporal and spatial scales owing to severe requirements of computational resources (e.g., Pohl et al. 2020).
In addition, adopting the same setup as in the linear analysis but only for the models with β = 50 and m i /m e = 50, we carry out 2D PIC simulations with periodic boundary conditions (periodic-box simulations, hereafter) to study the nonlinear evolution of the plasma waves excited by such microinstabilities. Note that throughout the paper we refer two different sets 1 In HKRK2021 and hereafter, the transition zone is defined as the downstream region of r L,i behind the shock ramp, where r L,i ≈ u 0 /Ω up ci is the gyroradius of incoming ions; u 0 is the preshock flow speed defined in the downstream rest frame and Ω up ci is the gyro-frequency in the upstream. Both the first and second overshoots and the accompanying undershoot are included in this zone, beyond which the downstream states satisfy the canonical Rankine Hugoniot relation (see Figure 1 of HKRK2021). of PIC simulations: (1) The 'periodic-box simulations' are designed to study the nonlinear evolution of the excited plasma waves in the same set-up as in the linear analysis, and will be presented in Section 3. (2) The 'shock simulations' reported in HKRK2021 provide the numerical estimates for the ion and electron temperature anisotropies in the shock transition zone.
The paper is organized as follows. Section 2 describes the linear analysis of the AIC, whistler, and mirror instabilities. In Section 3, we present the evolution of the waves driven by these instabilities in 2D periodic-box simulations. In Section 4, the implication of our work on the shock criticality and shock surface ripples is discussed. A brief summary is given in Section 5.

Basic Equations
We consider a homogeneous, collisionless, magnetized plasma, which is specified by the density and temperature of ions and electrons, n i , n e , T i , T e , and the background magnetic field B 0 . The linear dispersion relation of general electromagnetic (EM) modes is given as where the dielectric tensor, ij , is determined by the plasma parameters and the velocity distribution functions (VDFs) of particles. Here, k i and k j are the components of the wavevector k. Then, the complex frequency, ω = ω r + iγ, 2 can be calculated as a function of the wavenumber, k, and the propagating angle, θ, between k and B 0 . Without loss of generality, we set B 0 = B 0ẑ along the +z direction and k = k xx + k zẑ in the x-z plane, as schematically illustrated in Figure  1(a). In order to compute ij , we adopt the VDFs with bi-Maxwellian distributions for ions and electrons: where v ⊥ = v 2 x + v 2 y and v = v z . The subscript a denotes e or i defined as the electron or ion species, respectively. Here, n 0 is the number density of electrons or ions, which satisfies the charge neutrality condition, i.e., n 0 = n e = n i . The parallel and perpendicular (to B 0 ) thermal velocities are v T a = 2k B T a /m a and v T a⊥ = 2k B T a⊥ /m a , respectively. Then, the perpendicular temperature anisotropy of each particle species is given as A a ≡ T a⊥ /T a = v 2 T a⊥ /v 2 T a . The schematic configuration of the thermal velocity ellipsoid for a bi-Maxwellian VDF with the temperature anisotropy A a is shown in Figure 1(b). As A a increases, the thermal velocity surface in the velocity space deviates further away from the spherical shape. Under these considerations, ij is given as Equation (3) in Kim et al. (2020) without the bulk drift velocities. We note that in the shock simulations of HKRK2021, the VDFs of ions and electrons in the transition zone are likely non-gyrotropic and non-Maxwellian due to the SDA-reflected ions and electrons accelerated via the gradient-B drift (see Figure  4 of Guo et al. (2014a)). However, we expect the effects of non-Maxwellian VDFs on the linear predictions would be only marginal, since the fractions of particles in the suprathermal tail are order of 10 −2 for ions and 10 −2 for electrons in the downstream region of the fiducial M s = 3 shock model (KRH2019).
For n 0 , B 0 , T a , and T a⊥ of the homogeneous background plasma, we adopt the numerical values, averaged over the transition zone of the simulated shocks of HKRK2021, where the preshock conditions are specified with the typical parameters of high-β ICM plasmas, n ICM = 10 −4 cm −3 , k B T ICM = (k B T i + k B T e )/2 = 8.6 KeV, and β ICM = 20 − 100. Again, in the shock simulations, both the ion and electron distributions are spatially nonuniform in the transition zone, where the flow structure oscillates with overshoots and undershoots in the longitudinal direction and ripples in the transverse direction. Hence, we focus on qualitative analyses of the instabilities rather than making precise quantitative predictions.
Throughout the paper, the plasma beta, β a = 8πn a k B T a /B 2 0 , the plasma frequency, ω 2 pa = 4πn a e 2 /m a , and the gyro-frequency, Ω ca = eB 0 /m a c, for electrons and ions are used. Note that in HKRK2021 the results are expressed in terms of the upstream parameters, n up 0 ≈ n 0 /r and B up 0 ≈ B 0 /r, where r is the density compression ratio across the shock ramp. So for example ω up pa ≈ ω pa / √ r and Ω up ca ≈ Ω ca /r. Plasma waves are characterized with the growth rate, γ, and the real frequency, ω r , which are calculated by solving the dispersion relation in Equation (1) for wavevector k. If the propagation angle of the wave with the maximum growth rate, γ m , is θ m ≈ 0 • , the wave mode is called 'parallel-propagating'. If θ m 0 • , it is 'oblique-propagating'. If the wave frequency, ω r ≈ 0, the mode is 'non-propagating'. The wave polarization, P , can be estimated also using the solution of the dispersion relation as follows: where δE ± ≡ δE x k,ω ∓ iδE y k,ω (Verscharen & Chandran 2013). The left-hand circular polarization (LHCP) corresponds to P = −1, whereas the right-hand circular polarization (RHCP) corresponds to P = +1. Waves are in general elliptically polarized with P = ±1. In the case of non-propagating mode (ω r = 0), P = 0 (see Table 1).

Linear Properties of AIC, Whistler and Mirror Instabilities
In this section, we report the results of the linear stability analysis for the microinstabilities triggered by the ion and electron temperature anisotropies in the transition region of high-β, Q ⊥ -shocks. The first column of Table 2 lists the model name, which is assigned with the two parameters, the shock Mach number, M s , and β (≈ β e + β i ) in the shock transition region. For example, the LM3.0β50 model has M s = 3.0 and β ≈ 50. The values of the parameters, β e , β i , A e , and A i , are listed in the 3 − 6 columns of the table. For the models of β ≈ 20 − 100 and the mass ration m i /m e = 50, they are obtained with n 0 , B 0 , T a , and T a⊥ estimated by averaging the numerical values over the transition region in the simulated shocks with M s = 2 − 3 and β up = 20 − 100 of HKRK2021. 3 Considering the uncertainties in averaging over nonlinear structures with a The quantities, βe, βi, Ae, and Ai, are obtained by the averaging numerical values over the transition zone in the simulated Q ⊥ -shocks presented in HKRK2021. b Linear predictions for the fastest growing mode of the ion-driven instabilities, (γm/Ωci,ckm/ωpi,θm), normalized with the ion gyro and plasma frequencies. θm is given in units of degree. c Linear predictions for the fastest growing mode of the electron-driven instabilities, (γm/Ωce,ckm/ωpe,θm), normalized with the electron gyro and plasma frequencies. θm is given in units of degree.
overshoot/undershoot oscillations, they are given only up to two significant figures.
For the models with higher mass ratios, LM3.0β50-m100 with m i /m e = 100 and LM3.0β50-m1836 with m i /m e = 1836, the parameters for the LM3.0β50 model (β e = 19, β i = 31, A e = 1.2, and A i = 2.0) are used only for the linear analysis. Also we carried out two additional shock simulations for M3.0β1 with β = 1 and M3.0β5 with β = 5, which were not considered in HKRK2021, in order to obtain the parameters to be used for LM3.0β1 and LM3.0β5. Our fiducial models have m i /m e = 50, which is adopted in order to ease the requirements of computational resources for the periodicbox PIC simulations that will be described in Section 3.
The linear predictions for the AIC, whistler, i-mirror, and e-mirror instabilities are given in the 8−11 columns of Table 2. The three numbers inside each parenthesis present the linear properties of the fastest growing mode: (γ m /Ω ci ,ck m /ω pi ,θ m ) for the AIC and i-mirror instabilities, and (γ m /Ω ce ,ck m /ω pe ,θ m ) for the whistler and e-mirror instabilities. Here, k m is the wavenumber that has the maximum growth rate γ m at θ m , and θ m is given in units of degree. For a clear distinction between the ion and electron mirror modes, in the 10 − 11 columns, γ m of each mirror mode, obtained with either isotropic electrons (i.e., A e = 1, A i > 1) or isotropic ions (i.e., A i = 1, A e > 1), is shown. Note that 'stable' means that waves cannot grow because γ m < 0, and 'quasi-stable' corresponds to γ m /Ω ci < 10 −2 . Figure 2 shows the linear analysis results for the LM3.0β50-m1836 model. For the adopted parameters, ω pe /Ω ce = 26. Panels (a)-(c) display the growth (or damping) rate at θ m of the AIC, whistler, and mirror instabilities, respectively, as a function of the wavenumber. To make a uniform comparison, γ and k are normalized with Ω ci and ω pi /c, respectively, for both the ion-driven and electron-driven instabilities. Note that in panel (c) both γ and k are given in the logarithmic scales, in order to show both the i-mirror and e-mirror modes in the same panel. To examine the effects of A i and A e separately and also their combination, we present the black dashed lines for the case with both the ion and electron anisotropies, the red solid lines with the ion anisotropy only, and the blue solid lines with the electron anisotropy only.
The AIC instability induces quasi-parallel modes with θ m = 0 • . Although A i > 1 is the main free energy source which drives the AIC instability, we find that A e > 1 reduces the growth rate (see the red and black lines in panel (a) and also Ahmadi et al. (2016)). By contrast, the whistler instability is unstable for A e > 1, Note that γ and ωr are normalized with Ωci and k is normalized with ωpi/c, uniformly for both the ion and electron modes. and the growth rate is independent of A i . The whistler mode is also quasi-parallel propagating with θ m = 0 • . The mirror modes, on the other hand, are highly oblique with θ m = 64 • for LM3.0β50-m1836. The e-mirror mode (blue) at high-k (ck/ω pi > 0.3) grows much faster than the i-mirror mode (red) at low-k (ck/ω pi < 0.3). With both A i > 1 and A e > 1, a mixture of the two mirror modes appears in the intermediate-k regime (ck/ω pi ∼ 0.3).
In the LM3.0β50-m1836 model, the maximum growth rates are given in the following order: where γ WI , γ EM , γ AIC and γ IM are the maximum growth rates of the whistler, e-mirror, AIC and i-mirror instabilities, respectively. Note that in general γ WI > γ EM (Gary & Karimabadi 2006), and γ AIC > γ IM under space plasma conditions with low-β and large temperature anisotropies (Gary 1992(Gary , 1993. The real frequency, ω r /Ω ci , at θ m for the mixed case (A e = 1.2 and A i = 2.0) are shown in panels (d)-(f) of Figure 2. The AIC-driven mode has ω r /Ω ci ∼ 0.25 − 0.5 for ck/ω pi ∼ 0.1 − 0.4, while the whistler mode has ω r /Ω ci ∼ 80 − 350 for ck/ω pi ∼ 5 − 20. The mirror modes are non-propagating or purely growing with ω r = 0. Moreover, the polarization, calculated using the solutions of the dispersion relation, is P = −1, +1, and 0 for the AIC, whistler, and mirror instabilities, respectively, as expected.

Parameter Dependence of Linear Properties
As listed in Table 2, we consider a number of models to explore the dependence on m i /m e and β. The upper panels of Figure 3 show the linear predictions for the models with M s = 3, β = 50, and m i /m e = 50 − 1836, while the lower panels are for the models with M s = 3, m i /m e = 50, and β = 1 − 100. For a higher mass ratio, electrons go through more gyro-motions per the ion gyro-time, Ω −1 ci . Nevertheless, γ AIC /Ω ci and γ EM /Ω ce are almost independent of m i /m e . In the case of the whistler and i-mirror instabilities, on the other hand, overall, the normalized growth rates are slightly lower for smaller m i /m e . Also the damping rate for the whistler instability is slightly higher for smaller m i /m e in the small wavenumber regime (ck/ω pe ∼ 0.1). As a result, the growth of the whistler and i-mirror instabilities may be somewhat suppressed in the shock simulations with reduced mass ratios. However, even in the case of m i /m e = 50, this effect is expected to be only minor, because the inequality in Equation (4) Table 2. Note that for the mirror modes, θm depends on β, but not on mi/me. and the changes of k m and θ m are negligible (see Table  2).
The plasma beta is another important parameter that affects the stability of the system. Note that the anisotropy parameters, A e and A i , are almost independent of β for β ≈ 20 − 100, the range relevant for ICM shocks (see Table 2), although they tend to increase slightly with increasing β in the second digit to the right of the decimal point. In the low-β case (LM3.0β1), A i = 1.2 is significantly smaller than those of other highβ models due to the strong magnetization of ions. This is because A i in the shock transition is closely related to the fraction of reflected ions. On the other hand, A e in the shock transition is not substantially affected by β, because it is mainly determined by the magnetic field compression rather than the fraction of reflected electrons. Given the same temperature anisotropies, the growth of the instabilities tends to be suppressed by strong magnetic fields at low-β plasmas. As can be seen in the lower panels of Figure 3, the peak values of either γ m /Ω ci for the ion-driven modes or γ m /Ω ce for the electron-driven modes increase with increasing β. For the AIC, whistler, and i-mirror instabilities, γ m /Ω ci or γ m /Ω ce occurs at smaller ck/ω pi or ck/ω pe , for higher β. But such a trend is not obvious in the case of the e-mirror mode.
In the high-β cases (β ≈ 20 − 100) with M s = 3, all the AIC, whistler, i-mirror and e-mirror waves can be triggered, as shown in the lower panels of Figure 3, leading to the generation of multi-scale waves from elec-tron to ion scales. On the other hand, in the LM3.0β5 model (red solid lines), the e-mirror mode is stable, but other modes are unstable. In the LM3.0β1 model (gray solid lines), all the instabilities are stable with negative growth rates.
The sonic Mach number, M s , is the key parameter that determines the temperature anisotropies in the transition of high-β ICM shocks (β ≈ 20 − 100), since the ion reflection fraction and the magnetic field compression are closely related to M s . Figure 4 shows the growth rates of the instabilities for M s = 2.0 (black), 2.3 (red), and 3.0 (blue), in the cases of β = 20 (top), 50 (middle) and 100 (bottom). As M s increases, both A e and A i increase, so all the modes grow faster and k m shifts towards larger k, regardless of β.
Note that the AIC and whistler modes have γ m at θ m = 0 independent of M s , whereas θ m decreases with increasing M s for the i-mirror and e-mirror modes (see also Table 2). In LM2.0β50 and LM2.0β100, the AIC instability is stable or quasi-stable, while the whistler and mirror modes can grow. In the case of LM2.0β20, all the instabilities are stable (see black lines in top panels). In the models with M s = 2.3 − 3 (red and blue lines), on the other hand, the four instabilities are unstable, and hence multi-scale plasma waves can be generated.
The parameter dependence can be summarized as follows. (1) For the AIC mode, the maximum normalized growth rate, γ m /Ω ca , and the corresponding normalized wavenumber, ck m /ω pa , are almost independent of m i /m e . For the whistler, i-mirror, and e-mirror modes,   on the other hand, γ m /Ω ca is only slightly enhanced for larger m i /m e , whereas ck m /ω pa exhibits almost no dependence.
(2) For all the modes, the overall trend shows that γ m /Ω ca is higher and ck m /ω pa is smaller for higher β.
(3) For all the modes, γ m /Ω ca is higher and ck m /ω pa is larger for higher M s cases with larger A e and larger A i .

Numerical Setup
To investigate the development and nonlinear evolution of the instabilities, we performed 2D PIC simulations with periodic boundary conditions for the three fiducial models, LM2.0β50, LM2.3β50 and LM3.0β50, with the same setup described in Section 2.1. Electrons and ions are prescribed with bi-Maxwellian VDFs with β e , β i , A e , A i given in Table 2. As noted before, here m i /m e = 50 is employed due to the computational limitations, but at least the early, linear-stage develop-ment of the plasma instabilities under consideration is expected to depend rather weakly on the mass ratio.
We point that the setup for these 2D periodic-box simulations should intrinsically differ from the condition in the transition zone of shocks in the following aspects.
(1) The initial distributions of ions and electrons are allowed to relax in the periodic-box simulations. By contrast, the shock-reflected ions and electrons are continuously convected into the transition zone behind the shock ramp, leading to the continuous excitation of the instabilities. (2) Homogeneous spatial distributions and bi-Maxwellian VDFs are assumed for the periodic-box simulations. On the other hand, as noted in Section 2.1, both the ion and electron distributions are spatially nonuniform and the VDFs of ions and electrons contain suprathermal tails in the shock transition region. Nevertheless, the periodic-box simulations like ours are often used to investigate the nonlinear evolution and properties of microinstabilities in the either upstream or downstream region near the shock (e.g. Scholer et al. 2000;Guo et al. 2014b;Trotta et al. 2020).  The simulations were carried out using a parallelized EM PIC code, TRISTAN-MP (Buneman 1993;Spitkovsky 2005). The simulation domain is a square of box size L x = L z = 84.8c/ω pi = 600c/ω pe in the x-z plane, which consists of the grid cells of ∆x = ∆z = 0.1c/ω pe . In each cell, 32 particles (16 for ions and 16 for electrons) are placed. The time step of the simulations is ∆t = 0.045ω −1 pe , and the simulations ran up to t end = 130Ω −1 ci . Interpreting the results of our PIC simulations could be limited by numerical noises and aliases due to a finite number of macroparticles on discretized grids (e.g. Pohl et al. 2020). However, we expect that the overall results of the PIC simulations are reasonably converged, judging from the previous studies (Guo et al. 2014b;Kang et al. 2019;Kim et al. 2020).

Results of Periodic-Box Simulations
With the inequality in Equation (4), we expect that the whistler mode grows much faster than other modes, resulting in the relaxation of A e during the early stage. As the whistler and e-mirror modes grow and then decay on the time scale of τ WI ≡ 1/γ WI , the AIC and i-mirror modes become dominant later on the time scale of τ AIC ≡ 1/γ AIC . Figure 5 shows the magnetic field fluctuations, δB y (upper panels) and δB z (lower panels), in the x-z plane (simulation plane) at three different times in the LM3.0β50 model. Here, the growth time scales, τ WI and τ AIC , are estimated with γ m of each mode in Table 2. At t ∼ τ WI , the transverse component, δB y , appears on electron scales and the waves containing it propagate parallel to B 0 in panel (a), but the longitudinal component, δB z , does not grow significantly in panel (d). In this early stage, the dominant mode is the whistler mode, while the e-mirror mode is much weak to be clearly manifested. As A e decreases in time due to the electron scattering off the excited waves, the whistler waves decay as shown in panel (b). On the time scale of τ AIC , both the AIC and i-mirror instabilities grow and become dominant. It is clear that the AIC-driven waves, shown in panel (c), are parallel-propagating, while the i-mirror-driven waves, shown in panel (f), are obliquepropagating; the blue arrow in the bottom-left corner of panel (f) denotes the wavevector of the i-mirror-driven mode with the maximum growth rate. Figure 6 shows the time evolution of the power spectrum for the magnetic field fluctuations, δB 2 y (k), for LM2.0β50, LM2.3β50, and LM3.0β50 at t ∼ τ WI , t ∼ . Power spectra of the magnetic field fluctuations, δB 2 y (k), in the period-box simulations for LM3.0β50 (top), LM2.3β50 (middle), and LM2.0β50 (bottom), plotted in the k -k ⊥ (that is, kz −kx) plane. The results are shown at t ∼ τWI (left), t ∼ 3τWI (middle), and t ∼ τAIC (right). See the text for the remarks on τAIC for LM2.0β50. The gray star symbol marks the location of the maximum linear growth rate, γm, estimated from the linear analysis. In the models with Ms ≥ 2.3, AIC, whistler and i-mirror waves appear, while those waves do not grow substantially in the model with Ms = 2.
3τ WI , and t ∼ τ AIC . Again, the growth time scale of each mode is estimated with γ m listed in Table 2, except for the LM2.0β50 model, in which the AIC instability is stable, and so the output time of panel (i) is chosen at the evolutionary stage similar to that of LM2.3β50. In the cases of M s = 2.3 and 3, whistler waves are excited dominantly at quasi-parallel propagating angles at t ∼ τ WI . After the initial linear stage, the energy of the whistler waves is transferred to smaller wavenumbers and the waves gradually decay, as shown in panels (b) and (e). On the time scale of ∼ τ AIC , AIC waves and i-mirror waves appear dominantly at quasi-parallel and highly oblique angles, respectively, as shown in panels (c) and (f). This is consistent with the evolutionary behavior which we have described with Figure 5. For the AIC and whistler instabilities, the linear predictions for k m with the maximum growth rate (gray star symbols) agree reasonably well with the peak locations of the magnetic power spectrum realized in the PIC simulations. But the linear estimates for the i-mirror mode are slightly off, because γ m is obtained without the electron anisotropy, as stated in Section 2.2. In summary, the results of the periodic-box simulations are quite consistent with the linear predictions described earlier. Also we note that the results of our PIC simulations are in good agreement with those of Ahmadi et al. (2016), in which PIC simulations were carried out to explore the evolution of the instabilities due to the temperature anisotropies in space plasmas with β ∼ 1. The bottom panels of Figure 6 confirm that waves do not grow noticeably in the LM2.0β50 model. In these periodic-box simulations, the electron-scale waves develop first and then decay as A e is relaxed in the early stage, followed by the growth of the ion-scale waves due to A i . In the shock transition region, by contrast, temperature anisotropies are to be supplied continuously by newly reflected-gyrating ions and magnetic field compression, hence multi-scale plasma waves from electron to ion scales are expected to be simultaneously present.

Shock Criticality
As mentioned in the introduction, the Fermi-like acceleration, which relies on the upstream waves excited by the EFI, is effective only for supercritical shocks with M s ≥ M * EFI ≈ 2.3 in β ≈ 20 − 100 plasma (Guo et al. 2014b;Kang et al. 2019). The SSDA, which depends on the multi-scale waves excited mainly by the AIC and whistler instabilities, is thought to occur in supercritical shocks with M s ≥ M * AIC ≈ 3.5 in β ≈ 1 plasmas  and M s ≥ M * AIC ≈ 2.3 in β ≈ 20 − 100 plasmas (HKRK2021). We suggest that both M * EFI and M * AIC are related to the sonic critical Mach number, M * s , for ion reflection, since the structure of collisionless shocks is governed primarily by the dynamics of shockreflected ions. Table 3 summarizes the shock criticality of the simulated shock models and the stability of the linear analysis models. The first column lists the name of the simulated shock models considered in HKRK2021, and the two additional models for low-β shocks performed for this study. The shock criticality of each model is given in the second column. The name of the corresponding linear analysis models is given in the third column, while the last four columns show the stability for the four instabilities (see also Table 2). We note that the name of the shock models includes β up in the preshock, upstream plasmas, while that of the linear analysis models includes β in the shock transition zone given in Table 2.
As discussed in Sections 2 and 3, in β ≈ 20 − 100 plasmas, the AIC instability operates for M s 2.3, while whistler waves are induced regardless of M s . In the M3.0β5 model, the e-mirror mode is stable, while the other three modes are unstable. This is in good agreement with the 2D simulation of a M s = 5 and β = 5 shock reported earlier by Niemiec et al. (2019). In the M3.0β1 model, by contrast, A i is smaller than that of high-β models, and all the four instability modes are suppressed by strong magnetization. This is consistent with the results of M * AIC ≈ 3.5 for shocks with β ≈ 1 presented by .

Shock Surface Rippling
Another important feature of supercritical shocks above M * AIC is the shock surface rippling. According to previous shock simulations (e.g., Lowe & Burgess 2003;Matsukiyo & Matsumoto 2015;Niemiec et al. 2019;, the rippling has the characters of AIC waves with the fastest growing mode at θ m ∼ 0, the propagation speed close to the local Alfvén speed, and the wavelengths of ∼ λ AIC (≈ 30c/ω pi ). In fact, shock ripples have been observed at the Earth's bow shock and the interplanetary shocks inside the heliosphere, and investigated extensively in space physics (e.g. Winske & Quest 1988;Moullard et al. 2006;Johlander et al. 2016). At the interplanetary shocks, ripples on scales even larger than λ AIC have been detected as well, and are thought to be triggered by the upstream magnetic structures produced by backstreaming ions (Kajdič et al. 2019).
The parallel-propagating AIC and whistler waves in homogeneous plasmas are purely electromagnetic and incompressible with both the electric and magnetic wave vectors pointing normal to B 0 . The fluctuating magnetic fields of oblique mirror modes, on the other hand, have a substantial longitudinal component, that is, δB has a significant component parallel to B 0 (Gary 1993). Since the density fluctuations are proportional to the parallel electric and magnetic field fluctuations (Hojo et al. 1993), we expect to see only weak ion density fluctuations due to the i-mirror mode in our 2D periodic-box simulations.
Panel ( However, the ion density fluctuations of the rippling waves propagating along the shock surface behind the shock ramp are rather significant in the shock simulation for the M3.0β50 model in HKRK2021. Panel (c) shows that both the variations of [ n i −n 0 x,avg /n 0 ] ≈ ±0.1 and [ B y −B 0 x,avg /B 0 ] ≈ ±0.1 have similar amplitudes; the fluctuations of n i are much larger than those of the linear prediction expected for the parallel-propagating AIC mode. Note that here the quantities are averaged along the x-direction over the shock transition zone including the first and second overshoot oscillations behind the ramp. Hence, the basic assumptions of the linear theory, such as the homogeneous background, charge neutrality, and zero net-current, are likely to be violated in this region.
We point that such large-amplitude fluctuations of n i , comparable to the fluctuations of B y , were previously recognized in the 2D hybrid simulations of supercritical, perpendicular shocks presented in Winske & Quest (1988). The authors suggested that the large compressive waves might result from nonlinear effects in addition to oblique i-mirror modes. The effects due to nonlinear couplings between various wave modes could be significant as well (e.g. Shukla & Stenflo 1985;. Therefore, the pure AIC-driven waves in the shock transition could have been modified by such possible nonlinearities, leading to the enhancement of ion density fluctuations. Variations in the transverse component of magnetic field, By − B0 x,avg/B0, and the iondensity, ni − n0 x,avg/n0, averaged over the x-domain in the 2D periodic-box simulation for the LM3.0β50 model, plotted along B0 (z-direction) at three different times. (c): Variations in the longitudinal component of magnetic field, By−B0 x,avg/B0 (red), and the ion-density, ni−n0 x,avg/n0 (black), averaged along the x-direction over the shock transition zone in the 2D shock simulation for the M3.0β50 model in HKRK2021, plotted along the y-direction. Note that the preshock magnetic field, B up 0 , lies in the x-y plane, and the obliquity angle between B up 0 and the y-axis is θBn = 63 • in the shock simulation.

SUMMARY
In supercritical Q ⊥ -shocks, a substantial fraction of incoming ions and electrons are reflected, and the transverse components of magnetic fields are amplified at the shock ramp. The reflected-gyrating ions and the amplified magnetic fields induce the ion and electron perpendicular temperature anisotropies, A i and A e , respectively, in the shock transition region (Guo et al. 2017). They in turn trigger various microinstabilities including the AIC, whistler, i-mirror, and e-mirror instabili-ties (e.g., Winske & Quest 1988;Lowe & Burgess 2003;Guo et al. 2017). The kinetic properties of these four instabilities are summarized in Table 1. The multi-scale plasma waves generated by these microinstabilities are thought to be crucial for the electron preacceleration via the SSDA (e.g., Katou & Amano 2019;Niemiec et al. 2019;. In this work, adopting the numerical estimates for the ion and electron temperature anisotropies found in the 2D PIC simulations of Q ⊥ -shocks with M s = 2 − 3 (see Table 2), we have carried out the kinetic linear analysis of the microinstabilities for wide ranges of parameters, β = 1 − 100 and m i /m e = 50 − 1836. The linear predictions for the fastest growing mode, γ m , k m , θ m , of each instability are given in Table 2. In addition, in order to investigate the development and nonlinear evolution of the waves induced by the microinstabilities, we have performed 2D PIC simulations with periodic boundary conditions for the three fiducial models, LM2.0β50, LM2.3β50, and LM3.0β50. Finally, the results have been also compared with the 2D PIC simulations for ICM shocks reported in HKRK2021.
The main results can be summarized as follows: 1. In the LM3.0β50-m1836 model with the real mass ratio, which represents a typical supercritical ICM shock, the maximum growth rates of the four instabilities have the following order: γ WI γ EM γ AIC > γ IM (Fig. 2). Hence, the parallel-propagating AIC and whistler waves are expected to be more dominant than the obliquepropagating mirror waves. 2. In the LM2.0β50 model, which represents a subcritical ICM shock, by contrast, the AIC mode is stable (Table 2), so mainly the electron-scale whistler waves are generated. 3. The maximum normalized growth rates for the AIC and e-mirror instabilities, γ AIC /Ω ci and γ EM /Ω ce , are almost independent of m i /m e , while γ WI /Ω ce for the whistler instability and γ IM /Ω ci for the i-mirror instability are slightly lower for smaller m i /m e (Fig. 3). 4. For all the four instabilities, the maximum normalized growth rates increase with increasing β (Fig. 3). 5. As the sonic Mach number M s increases, both A e and A i increase, all the modes grow faster, and k m of each mode shifts towards larger k, regardless of β in the range of β ≈ 20 − 100 (Fig. 4). 6. The critical sonic Mach number to trigger the AIC instability in the shock transition is M * AIC ∼ 2.3 for β ≈ 20−100. It is similar to the EFI critical number suggested in our previous study, that is, M * AIC ≈ M * ef ≈ 2.3 (Kang et al. 2019). For β = 1, on the other hand, M * AIC 3 is slightly higher, because the AIC mode is suppressed by the strong magnetization of ions (Hellinger & Mangeney 1997;. 7. The 2D periodic-box simulations confirm the linear predictions. In the early stage of ∼ τ WI , electron-scale waves develop and then decay as A e is relaxed, followed by the growth of ion-scale waves on the time scale of ∼ τ AIC (Figs. 5 and 6). 8. The rippling waves propagating along the shock surface have the characteristics of AIC waves. Although the AIC waves are parallel-propagating, electromagnetic, incompressible in the linear regime, the amplitudes of the longitudinal magnetic field and ion-density fluctuations associated with the overshoots in the shock transition are similar and of the order of 10% according to the shock simulation for the M3.0β50 model (Fig. 7). It is expected that the inhomogeneity in the shock transition and the nonlinear effects could lead to the generation of such large-amplitude fluctuations of the ion-density along the shock surface.
In conclusion, our results well support the suggestion for the generation of multi-scale plasma waves via various microinstabilities in the transition region of highβ, supercritical, Q ⊥ -shocks (Guo et al. 2017;Katou & Amano 2019;Niemiec et al. 2019;. A detailed description of the shock structure and the electron preacceleration in such ICM shocks, realized in 2D PIC simulations, is reported in HKRK2021.