Secondary Whistler and Ion-cyclotron Instabilities Driven by Mirror Modes in Galaxy Clusters

Electron cyclotron waves (whistlers) are commonly observed in plasmas near Earth and the solar wind. In the presence of nonlinear mirror modes, bursts of whistlers, usually called lion roars, have been observed within low magnetic field regions associated with these modes. In the intracluster medium (ICM) of galaxy clusters, the excitation of the mirror instability is expected, but it is not yet clear whether electron and ion cyclotron (IC) waves can also be present under conditions where gas pressure dominates over magnetic pressure (high β). In this work, we perform fully kinetic particle-in-cell simulations of a plasma subject to a continuous amplification of the mean magnetic field B (t) to study the nonlinear stages of the mirror instability and the ensuing excitation of whistler and IC waves under ICM conditions. Once mirror modes reach nonlinear amplitudes, both whistler and IC waves start to emerge simultaneously, with subdominant amplitudes, propagating in low- B regions, quasi-parallel to B (t). We show that the underlying source of excitation is the pressure anisotropy of electrons and ions trapped in mirror modes with loss-cone-type distributions. We also observe that IC waves play an essential role in regulating the ion pressure anisotropy at nonlinear stages. We argue that whistler and IC waves are a concomitant feature at late stages of the mirror instability even at high β, and therefore, expected to be present in astrophysical environments like the ICM. We discuss the implications of our results for collisionless heating and dissipation of turbulence in the ICM.


INTRODUCTION
Several classes of astrophysical plasmas display fully developed turbulent states and a weak collisionality, in the sense that the particles' mean free path is several orders of magnitude larger than the typical radius at which they gyrate around the ambient magnetic field. These two characteristics alone can make the transport properties and global evolution of the astrophysical environment in question challenging and dependent on the local evolution at particles' scales. Therefore a detailed study of the behavior of these plasmas at the kinetic level becomes a necessity.
That is the case of the intracluster medium of galaxy clusters (ICM). The ICM is a hot, magnetized (Bonafede, A. et al. (2010)), weakly collisional and turbulent (Schuecker, P. et al. (2004); Zhuravleva et al. (2014);Hitomi Collaboration et al. (2016)) gas in the plasma state where the thermal pressure greatly exceeds the magnetic pressure (β ≡ 8πP/B 2 ∼ 10 − 100, P is the isotropic thermal pressure and B the magnetic field strength). In these conditions, departures from thermodynamic equilibrium, such as pressure anisotropies, are easy to achieve. For example, slow compression of the magnetic field increases particle kinetic energy perpendicular to the magnetic field such that the magnetic moment (or, the magnetic flux through the particle gyro-orbit) remains con-stant, leading to an excess of perpendicular pressure P ⊥ over parallel pressure P ∥ . However, pressure anisotropy cannot grow unchecked. Pressure anisotropies can easily destabilize microinstabilities such as mirror, firehose, ion-cyclotron and whistler (Schekochihin et al. (2005); Schekochihin & Cowley (2006)). The back reaction of these instabilities on the particles can maintain pressure anisotropy near its marginally unstable value, and are thought to play an important role in several aspects of ICM transport and heating (Kunz et al. (2011); Berlok et al. (2021); Drake et al. (2021); Perrone & Latter (2022a,b); Ley et al. (2023); Tran et al. (2023)).
In a similar vein, the solar wind and some regions of the Earth's magnetosheath and magnetosphere host plasmas that are also collisionless and turbulent. Even when the plasma β is lower than in the ICM (β i ∼ 1 − 10, β e ∼ 1), we can encounter some similarities. In particular, the plasma is also pressure anisotropic, and the same microinstabilities above mentioned are found to be present, usually in their fully developed, nonlinear stage (Bale et al. (2009)). Particularly important to this work is the presence of the mirror instability (Chandrasekhar et al. (1958); Rudakov & Sagdeev (1961); Hasegawa (1969); Southwood & Kivelson (1993); Kivelson & Southwood (1996); Pokhotelov et al. (2002Pokhotelov et al. ( , 2004) and its interplay with the whistler and (potentially) ion-cyclotron in-stabilities (Gary (1992), Gary & Wang (1996)). An example of this has been observed in these space plasmas, and termed whistler lion roars.
Whistler lion roars are short bursts of right-hand polarized waves, with frequencies below the electron cyclotron frequency (ω c,e ) commonly observed in the Earth's magnetosheath and magnetosphere (Smith et al. (1969); Tsurutani et al. (1982); Baumjohann et al. (1999); Breuillard et al. (2018); Giagkiozis et al. (2018); Kitamura et al. (2020); Zhang et al. (2021)), therefore identified as whistler waves. They have also been observed in Saturn's magnetosheath (Píša et al. (2018)) and the solar wind. They are observed in regions of locally low magnetic field strength (magnetic troughs, or magnetic holes) of magnetic fluctuations. These magnetic troughs are usually identified as structures produced by mirror instability modes, which are able to trap electrons with low parallel velocity within these regions due to the aforementioned invariance of magnetic moment (Southwood & Kivelson (1993)).
Several mechanisms have been proposed to explain the excitation of whistler lion roars. They usually invoke the pressure anisotropy P ⊥,e > P ∥,e that electrons generate while trapped inside the magnetic troughs (P ⊥,e and P ∥,e are, respectively, the electron pressure perpendicular and parallel with respect to the local magnetic field B). Other mechanisms have also been proposed involving counterpropagating electron beams inside these regions, and butterfly distributions in pitch-angle (Zhang et al. (2021); Jiang et al. (2022)). As the waves propagate out from the magnetic troughs, they are thought to interact with electrons, regulating the number of trapped electron inside magnetic troughs and also the global anisotropy of electrons in the magnetosheath. This way, there would be a causal connection between an ion-scale mirror instability with an electron scale whistler instability at nonlinear stages, providing valuable insight into the interaction of mirror modes with electrons.
The question arises as to whether a similar interplay can be expected in the ICM. Such behavior would imply a more complex scenario in which several microinstabilities would be causally connected and coexisting with each other, and several channels of turbulent energy dissipation would open, leading to a much richer dynamics.
Mirror instability and its consequences have been extensively studied using particle-in-cell (PIC) simulations of moderately and high-β plasmas, both hybrid (Kunz et al. (2014); Melville et al. (2016); Arzamasskiy et al. (2023)) and fully kinetic (Sironi & Narayan (2015); Riquelme et al. (2015Riquelme et al. ( , 2016; Ley et al. (2023)), up to nonlinear stages. Consistent with early theoretical works (Southwood & Kivelson (1993); Kivelson & Southwood (1996)), it has been demonstrated that mirror modes are efficient in trapping ions inside regions of low magnetic field strength during their secular growth (Kunz et al. (2014)). When mirror modes reach amplitudes of order δB/B ∼ 1, they reach a saturated stage and the ions eventually undergo scattering, allowing them to escape. This trapping process is similar for electrons, and it has been shown to have important consequences in the electron viscosity and thermal conduction of the plasma (Riquelme et al. (2016); Roberg-Clark et al. (2016). Interestingly, Riquelme et al. (2016) reported the observation of whistler waves in the nonlinear, saturated stages of mirror modes in their simulations, along with ion-cyclotron (IC) waves, although they did not pinpoint the cause of the excitation.
In this work, we use PIC simulations to investigate the nonlinear stages of the mirror instability at moderate and high-β, focusing on the abovementioned excitation of whistler and IC waves. We observe that, indeed, both right hand and left hand polarized, quasi parallel-propagating waves are excited at the end of mirror's secular growth and during its saturated stage, and provide evidence for their excitation mechanism associated to the pressure anisotropy electrons and ions within magnetic troughs of mirror modes. The right-and left-handed circular polarization of these waves lead to their identification as electron-cyclotron (i.e. whistlers) and ioncyclotron (IC) waves. We also provide some additional discussion about their nature. We describe the interaction of these waves with electrons and ions, and their effect on the regulation of the pressure anisotropy at late stages.
This paper is organized as follows. Section §2 describes our simulation setup and the runs we perform. Section §3 shows our simulation results starting from the excitation of the mirror instability, an early whistler burst and then the late excitation of the electron and ion cyclotron waves at nonlinear stages of the mirror instability. We also detail the mechanism by which these cyclotron waves are excited during the saturated stage of mirror modes, by tracking ions and electrons throughout the simulations. We also describe the subsequent interaction of these waves with the ions and electrons at late stages. In section §4 we discuss the dependence of our results on the mass ratio used in our simulations and show that they are fairly insensitive to it. In section §5 we present results of simulations at different initial ion plasma beta, and show these cyclotron waves are also present at lower and higher betas as well. Finally, we discuss the implication of our work in the context of galaxy clusters and present our conclusions in section §6.

SIMULATION SETUP
We perform fully kinetic, 2.5D particle-in-cell (PIC) simulations using TRISTAN-MP (Buneman (1993);Spitkovsky (2005)), in which we continuously shear a collisionless, magnetized plasma composed of ions and electrons (Riquelme et al. (2012)). The magnetic field is initially spatially uniform and starts pointing along the x-axis. A shear velocity field is imposed with v = −sxŷ (red arrows in fig. 1), where x is the distance along the x-axis and s is a constant shear rate. We solve the PIC system of equations using shearing coordinates, as implemented in Riquelme et al. (2012) (The suitability of this approach to studying ion Larmor scale phenomena is also discussed in Riquelme et al. (2015)). The conservation of magnetic flux implies that the y-component of the magnetic field B evolves as dB y /dt = −sB 0 , whereas dB x /dt = 0 and dB z /dt = 0. The action of the shear then Figure 1. The evolution of the simulation domain. Panel a: Initially, the box is straight, the magnetic field is initialized pointing in thex direction and a shear velocity field v = −sxŷ is imposed in the y-direction (red arrows). Panel b: The velocity field shears the box continuously throughout the simulation, amplifying the magnetic field and changing its direction in the process due to magnetic flux conservation. continuously amplifies the magnetic field strength such that its magnitude evolves as B(t) = B 0 √ 1 + s 2 t 2 . In our simulations, ions and electrons are initialized with Maxwell-Jüttner distributions (the relativistic generalization of the Maxwell-Boltzmann distribution, Jüttner (1911)) with equal initial temperatures T init i = T init e , and k B T init i /m i c 2 between 0.01 and 0.02. The physical parameters of our simulations are the initial temperature of ions and electrons (T init i = T init e ), the initial ion plasma beta, β init i , the mass ratio between ions and electrons m i /m e , and the ratio between the initial ion cyclotron frequency and the shear frequency, ω init c,i /s, that we call the "scale-separation ratio". The numerical parameters in our simulations are the number of macroparticles per cell, N ppc , the plasma skin depth in terms of grid point spacing, c/ ω 2 p,e + ω 2 p,i /∆x, and the domain size in terms of the initial ion Larmor radius, L/R init L,i , where R init L,i = v th,i /ω init c,i and v 2 th,i = k B T i /m i . These physical and numerical parameters are listed in Table 1. We fix c/ ω 2 p,e + ω 2 p,i /∆x = 3.5 in the simulations presented in Table 1.
In the bulk of the paper we discuss a representative, fiducial simulation with m i /m e = 8, β init i = 20 (thus β init = β init i + β init e = 40) and ω init c,i = 800 (simulation b20m8w200 in Table  1, highlighted in boldface). We vary the above parameters in a series of simulations, all listed in Table 1. Importantly, given the available computational capabilities, performing a simulation with realistic mass ratio m i /m e = 1836 becomes prohibitively expensive. Therefore, a range of values of ionto-electron mass ratio are presented in order to ensure that our results do not strongly depend on this parameter. The effects of varying these parameters are discussed in § §4 & 5.
In the absence of a scattering mechanism and/or collisions, the ion and electron magnetic moments µ j ≡ p 2 ⊥,j /(2m j B) and longitudinal action J j ≡ p j,∥ dℓ are adiabatic invariants (p ⊥,j and p ∥,j are the components of the momentum of a particle of species j perpendicular and parallel to the local magnetic field, respectively, and j = i, e), and therefore are conserved as the system evolves, provided that the variation of B is sufficiently slow compared to the particle cyclotron frequencies; in our case, s ≪ ω c,j , where ω c,j = eB/m j c is the cyclotron frequency of particles of species j, c is the speed of light, and e is the magnitude of the electric charge.
The continuous amplification of the magnetic field B implies that the particles' adiabatic invariance drives a pressure anisotropy in the plasma such that P ⊥,j > P ∥,j . In the very early stages of the simulation, we expect the evolution of P ⊥,j and P ∥,j to be dictated by the double-adiabatic scalings (Chew et al. (1956)). Soon after this stage, however, the pressure anisotropy acts as a free energy source in the plasma and is able to excite several kinetic microinstabilities after surpassing their excitation thresholds, which are proportional to β −α , (α ∼ 0.5 − 1) (Hasegawa (1969); Gary & Lee (1994); Gary & Wang (1996)). These microinstabilities break the adiabatic invariants and act upon the pressure anisotropy to regulate the anisotropy growth in the nonlinear stages.
In our simulations, and given our initial physical parameters (namely, β init i ≡ 8πP init i /B 2init = 20), we expect the dominant instability to be the mirror instability. Mirror modes are purely growing (i.e. zero real frequency), with the fastest growing modes propagating highly obliquely with respect to the mean magnetic field. Their most unstable wavenumbers satisfy k ⊥ R L,i ∼ 1, where R L,i is the ion Larmor radius. This instability presents Landau resonances with particles of very small parallel momentum, p ∥ ≈ 0, that become trapped in between mirror modes, and contribute to regulating the pressure anisotropy.
In addition to the mirror instability, we also observe wave activity that we associate with the ion-cyclotron (Gary (1992)) and whistler (Gary & Wang (1996)) instabilities at ion and electron scales, respectively, during the late stages of our simulations. Ion cyclotron (IC) modes are left circularly polarized and have real frequency below the ioncyclotron frequency ω c,i , with modes of maximum growth rate propagating parallel to the mean magnetic field B. Similarly, whistler modes are right circularly polarized and have real frequency below the electron cyclotron frequency ω c,e , with modes of maximum growth rate also propagating parallel to B. As we will see, this wave activity is associated with the ion and electron trapping processes that mirror modes generate.
3. RESULTS Figures 2 and 3 summarize the evolution of magnetic field fluctuations and particle pressure anisotropy over time. Figure 2 shows the fluctuations in the magnetic field δB ≡ B − ⟨B⟩ (where ⟨·⟩ denotes a volume average over the entire simulation domain) in its three different components at two different times: t · s = 0.4 (first row, panels a,b and c) and at t · s = 1.4 (second row, panels d, e and f ). The black arrows in panels a-f denote the direction of the mean magnetic field ⟨B⟩ at those particular times. The components of δB are defined as parallel with respect to the main field ⟨B⟩ (δB ∥ , panels b and e), perpendicular to ⟨B⟩ in the plane of the simulation (δB ⊥,xy , panels a and d) and perpendicular to ⟨B⟩ in the direction out of the simulation plane (δB z , panels c and f ). Additionally, figure 2g shows the evolution of the energy in each of the three components of δB, normalized by B(t) 2 ; δB 2 ∥ (blue line), δB 2 ⊥,xy (red line), and δB 2 z (green line). Figure 3a shows the evolution of the ion pressure anisotropy ∆P i ≡ P ⊥,i − P ∥,i for run b20m8w800, and the dashed gray line shows the approximate instability threshold for the mirror instability (Hasegawa (1969);Hellinger (2007)). We can see that the ion anisotropy surpasses the mirror threshold very early in the simulation, and reaches its maximum value at t · s ≈ 0.5 (we will call this stage the anisotropy overshoot hereafter). We will show that this is consistent with the beginning of the secular growth of mirror modes (Kunz et al. (2014), Riquelme et al. (2016)). Figure  3b shows the same for the electron pressure anisotropy, which we will show relaxes by efficient scattering.

Mirror Instability Evolution
Since mirror modes are highly oblique, their evolution is well represented by the time trace of δB 2 ∥ shown in fig. 2g. We identify both a linear, exponentially growing phase until t · s ≈ 0.45, and a subsequent nonlinear, slower growing secular phase, consistent with the different evolutionary phases of the ion and electron pressure anisotropies described above. Besides the break in the mirror mode's evolution at t · s ≈ 0.45, a second break in the secular growth occurs around t · s = 0.6 followed by a shallower slope of growth. We will show that this break coincides with the excitation of both whistler and IC waves in δB 2 ⊥,xy and δB 2 z , implying that whistler and IC waves, albeit smaller in amplitude, modulate the evolution of mirror modes during nonlinear stages.
3.1.1. Linear, exponentially growing mirror phase After an early CGL phase of the pressure anisotropy ∆P j (j = i, e, see fig. 3), fig. 2g shows the excitation of the mirror instability starting at t · s ≈ 0.35, mainly in the parallel component of the magnetic fluctuations, δB ∥ (blue line), consistent with theoretical expectations (Southwood & Kivelson (1993); Pokhotelov et al. (2004)). Figure 2g also shows that δB ∥ grows first and it has the largest amplitude throughout the entire simulation, meaning that the mirror instability is indeed the dominant instability. Figure 2b (i.e. δB 2 ∥ ) shows the linear, exponentially growing phase of mirror modes at t · s = 0.4, where small filamentary structures of high local magnetic field amplitude start to emerge and slowly grow, in between wider regions of low local magnetic field amplitude. The obliqueness of the modes is readily apparent, as well as the fact that the mirror generated magnetic fluctuations lie mainly in the (k,B) plane (they can be seen in δB 2 ⊥,xy too, but not in δB 2 z , as expected from linear theory (Pokhotelov et al. (2004))). The oblique nature of mirror modes can also be seen in fig. 4a, where we show the power spectrum in space of δB ∥ at t · s = 0.4. The solid and dashed lines represent the directions parallel and perpendicular to the mean magnetic field ⟨B⟩, respectively. Therefore, we can see that at t · s = 0.4, the power is mostly concentrated between wavevectors 0.44 ≲ kR init L,i ≲ 1.35 and angles of 52 • ≲ θ k ≲ 77 • , where θ k ≡ cos −1 (k · ⟨B⟩/kB) is the angle between mirror modes' wavevector and the mean magnetic field ⟨B⟩.
It should be emphasized that the ion-cyclotron wave activity only starts at t · s = 0.6, and not before. There is no sign of an early excitation of the ion-cyclotron instability competing with the mirror instability for the available free energy in ∆P i . Instead, at earlier stages, only the mirror instability is excited, consistent with our initial conditions of high-beta (β init i = 20), where the mirror instability is expected to dominate (e.g. Riquelme et al. (2015)).
The absence of ion-cyclotron waves early in the simulation (0 < t · s < 0.6) is clearly seen in fig. 5a, where we show the power spectrum in time and space of δB z (ω, k ∥ ) + iδB ⊥,xy (ω, k ∥ ) at early stages: 0.3 < t·s < 0.5. This particular combination of the two perpendicular components of δB allows us to disentangle the parallel-propagating waves (with respect to the main magnetic field ⟨B⟩, e.g. ion-cyclotron and whistlers), and also their left-handed and right-handed circular polarizations (Ley et al. (2019); Tran et al. (2023)). In this case, the left-hand circularly polarized wave activity is shown for ω > 0, whereas right-hand circularly polarized wave activity is shown for ω < 0. We readily see that, apart from the ω ≈ 0 power consistent with mirror modes appearing in δB ⊥,xy , there is no left-handed polarized wave activity throughout 0.3 < t · s < 0.5, only right-handed polarized waves, which corresponds to an early excitation of the whistler instability, as we will see in section 3.2.

Nonlinear, secular mirror phase
At t · s ≈ 0.45, we can clearly see the beginning of the secular growth of the mirror instability, where the modes reach nonlinear amplitudes, and keep growing but at a slower rate. This evolution is consistent with previous works (Kunz et al. (2014); Riquelme et al. (2016)).  The evolution of the ion pressure anisotropy ∆Pi/P ∥,i for run b20m8w800 is shown as a solid green line. The dashed green line shows the double-adiabatic evolution of ∆Pi/P ∥,i (Chew et al. (1956)). The dashed gray line shows the approximate threshold for the mirror instability: 1/β ∥,i (Hasegawa (1969)). The dotted-dashed orange line shows the threshold for the IC instability from Gary & Lee (1994) for γIC /ωc,i = 10 −2 (γIC is the IC growth rate). The red dashed line shows the bestfit to ∆Pi/P ∥,i = Aiβ α i ∥,i from t · s = 0.7 to t · s = 2.0, with Ai = 0.544 ± 0.003 and αi = 0.445 ± 0.003. Panel b: The evolution of the electron pressure anisotropy ∆Pe/P ∥,e is shown as solid orange line. The dashed orange line shows the double-adiabatic evolution of ∆Pe/P ∥,e . The dashed blue line shows the best-fit to ∆Pe/P ∥,e = Aeβ αe ∥,e from t · s = 0.7 to t · s = 2.0, with Ae = 0.036 ± 0.0002 and αe = 0.341 ± 0.003. The dashed gray line shows the linear threshold for the anisotropic whistler instability from (Gary & Wang (1996)) for growth rate γW /ωc,e = 0.01. (γW is the whistler growth rate).
Interestingly, the mirror secular growth is interrupted at t · s ≈ 0.6, and the slope of δB 2 ∥ breaks. This is also approximately where the ion pressure anisotropy experiences its fastest decline ( fig. 3). Mirror modes continue to grow, but at a much slower rate. This is consistent with the saturation of energy in the subdominant components δB 2 ⊥,xy and δB 2 z (solid red and green line in fig. 2g, respectively), which also presents a distinct pattern of oscillations. This activity is a clear evidence of a new burst of waves with components mainly in the direction perpendicular to δB, and we will see that they are consistent with both electron cyclotron waves (whistlers) and ion cyclotron waves excited by electron and ion populations, respectively, that become trapped within mirror modes (see sec. 3.3). Figure 2e shows a late, nonlinear stage of the mirror instability, at t·s = 1.4. At this time, the regions of high magnetic field of mirror modes (e.g. red filamentary structures seen in fig. 2b) have grown significantly and merged with neighboring structures to form wider and sharper regions of high local amplitudes (δB ∥ /B ∼ 0.9), whose sizes are comparable to regions of low magnetic field. At this stage, most of the power is concentrated in wavevectors 0.2 ≲ kR init L,i ≲ 1.1, and angles 57 • ≲ θ k ≲ 85 • (see fig. 4b).
After reaching its overshoot, the ion anisotropy starts to decrease towards marginal stability. However, this decrease stops around t · s ≈ 0.65 at ∆P i /P ∥,i ≈ 0.18, well above the approximate mirror threshold (dashed gray line, (Hasegawa (1969); Hellinger (2007))). The anisotropy then reaches a marginal stability level that is above the mirror threshold, similar to some previous works using both hybrid and fully kinetic simulations (Sironi & Narayan (2015); Melville et al. (2016); Ley et al. (2023)).
In order to better characterize the evolution of ∆P i , we fit a relation ∆P i = A i β αi ∥,i from 0.7 ≤ t · s ≤ 2 (In our simulations, the shear motion continuously amplifies B, therefore β ∥,i also evolves.). As shown in fig. 3a, our best-fit parameters are A i = 0.544 ± 0.003 and α i = −0.445 ± 0.003. The obtained exponent is consistent with marginal stability threshold given by the ion-cyclotron instability for lower β i (Gary & Lee (1994)). Indeed, the threshold for the IC instability, ∆P i = 0.53β −0.4 ∥,i , is plotted as dotted-dashed orange line in fig. 3a for γ IC /ω c,i = 10 −2 (Gary & Lee (1994)), and we can clearly see the similarity with our best-fit threshold, even at this higher value of initial β init ∥,i . This observation was also reported in Sironi & Narayan (2015), and we will see that, indeed, we do observe ion-cyclotron waves as part of the saturated phase of the mirror instability that starts at t · s = 0.6. The presence of ion and electron cyclotron waves coexisting with mirror modes at late, nonlinear stages of the mirror instability has been reported in previous works (Riquelme et al. (2016); Sironi & Narayan (2015); Ahmadi et al. (2018)). In §3.3, we argue that a natural explanation of the source of these cyclotron waves is pressure anisotropy of ions trapped within nonlinear mirror modes.   Figure 3b shows the evolution of the electron pressure anisotropy ∆P e ≡ P ⊥,e − P ∥,e for run b20m8w800. Initially, the electrons develop their own pressure anisotropy alongside ions and for the same reasons. The anisotropy follows double-adiabatic (CGL) scaling (dashed orange line) until t · s ≈ 0.4, when it has already reached a value significantly larger than the theoretical threshold for the growth of whistler modes, marked by grey-dashed lines (Gary & Wang (1996)). Around this time, the whistler instability starts to grow, as seen by the time trace of δB 2 z in fig. 2g, which is

First Whistler
. Panel a: The power spectrum of δBz(ω, k ∥ ) + iδB ∥,xy (ω, k ∥ ) in the entire simulation domain and between 0.3 < t · s < 0.5. The frequency is normalized by the initial electron cyclotron frequency ωc,e, and the wavevector is normalized by the plasma frequency ωp,e over the speed of light c. The solid black line shows the linear dispersion relation ωr(k) for the whistler instability according to our linear dispersion solver, whereas the dashed black line shows its growth rate γ. Panel b: The power spectrum in space of δBz(kx, ky) at t · s = 0.4. The wavenumbers kx, ky are normalized to the initial Larmor radius of the electrons, R init L,e . The solid and dashed black lines represent the direction parallel and perpendicular to the main magnetic field at that time. a rough proxy for whistler waves, and also because there are no left-handed IC waves as shown in fig. 5a. At t · s ≈ 0.45 the whistler modes saturate and enter a regime of quasisteady amplitude, which lasts until t · s ≈ 0.53. During this t · s ≈ 0.4 − 0.53 period, ∆P e is rapidly drawn down by frequent scattering, reaching a more slowly decreasing regime between t · s ≈ 0.53 and 0.6. The draw down of electron anisotropy happens at a time when the ion anisotropy is still growing. This lasts until mirror modes are sufficiently high amplitude to start trapping the electrons (t · s = 0.6).
The presence of whistler modes at t · s = 0.4 can be seen mainly in the perpendicular components of δB, namely, δB ⊥,xy and δB z , figures 2a and 2c, respectively. They propagate quasi-parallel to the main magnetic field B in a fairly homogeneous way inside the simulation domain. This quasi-parallel propagation can also be seen in fig. 5b, where we show the power spectrum in space of δB z (k x , k y ) at t · s = 0.4 for run b20m8w800, and the solid and dashed black lines indicate the directions parallel and perpendicular to the main magnetic field ⟨B⟩ at t · s = 0.4. The power of δB z (k x , k y ) is concentrated at parallel propagation and wavevectors 0.6 < kR init L,e < 1. We show the whistler wave frequencies in the power spec- fig. 5a. We can see that the power is localized in the region ω < 0, i.e. right-handed circularly polarized waves, consistent with the whistler polarization, and within frequencies 0.02 < ω/ω c,e < 0.05. As mentioned above, no IC activity is present during this time period.
We also calculated the theoretical dispersion relation of the anisotropic whistler instability using a linear dispersion solver assuming an initial bi-maxwellian distribution of electrons (Tran et al. (2023)), using the initial parameters and values of T ⊥,e , T ∥,e directly from the simulations. The dispersion relation ω(k) is shown as a solid black line in fig. 5a, whereas the instability growth rate is shown in dashed black lines. We can see that the power in right-hand circularly polarized waves is consistent with the whistler dispersion relation.
This way, the early evolution of the electrons is determined by an early burst of whistler modes associated to the initial electron pressure anisotropy growth. We will see that, once electrons start to become trapped in between mirror modes at t · s ≈ 0.6, another burst of whistler activity happens, this time associated with the trapping process within mirror modes during their secular and saturated phase.

Whistler and Ion-cyclotron Excitations
At the end of its secular growth, when mirror modes have reached sufficiently high-amplitudes, we simultaneously observe right-hand and left-hand circularly polarized wave activity, which we identify as whistler and ion-cyclotron waves, respectively. We will see below ( §3.3) that these whistler and ion-cyclotron waves propagate mainly in regions of locally low magnetic field (magnetic troughs). The source of this wave activity is identified to be the pressure anisotropic population of ions and electrons mainly due to trapped parti-cles inside the magnetic troughs. The whistlers and ion cyclotron waves then pitch-angle scatter both the trapped and untrapped particles, contributing to regulation of the global anisotropy. Figure 6 shows different spectral properties of the late burst of waves excited from t · s ≈ 0.6 onwards. Figure 6a shows the power spectrum in time of δB z (ω) + iδB ⊥,xy (ω) between 0.5 < t · s < 1.1, so we can see both left-hand (solid blue line) and right-hand (solid orange line) circular polarizations. The power spectrum peaks at low-frequencies, consistent with the nature of the dominant mirror modes (mainly appearing in δB ⊥,xy ). Additionally, we can clearly see a secondary peak at around ω ∼ 0.2ω c,i , with a spread that goes from ω ∼ 0.1ω c,i to ω ∼ 0.3ω c,i , in both left and right hand circular polarizations. This constitutes the characteristic feature informing the late burst of wave activity. This peak resembles observations of whistler lion roars in the Earth's Magnetosheath (see e.g. figs. 1 and 2 of Giagkiozis et al. (2018), fig. 3 of Zhang et al. (2021) for right-hand polarized waves.). Figure 6b shows the spectrogram of δB z (ω) + iδB ⊥,xy (ω) in frequency and time, ranging 0.4 < t · s < 1.3, with positive frequencies representing left-hand circularly polarized waves, and negative frequencies denoting right-hand circularly polarized waves. Here we can also see the early burst of whistler waves starting at t·s ≈ 0.4 and peaking at t·s ≈ 0.45 (see section §3.2), followed by the burst of both left-hand and right-hand circularly polarized waves at t·s ≈ 0.53 and peaking at t · s ≈ 0.65. This coincides with the rise in amplitude of δB 2 z and δB ⊥,xy (see fig. 2)g, and the waves are continuously maintained throughout the simulation at around the same frequencies.
Finally, figure 6c shows the power spectrum of δB z (ω, k ∥ ) + iδB ⊥,xy (ω, k ∥ ) in time and space, at 0.5 < t · s < 1.1. Frequencies and wavenumbers are normalized by ω c,i and ω p,i /c, respectively. Here we can also see the power at low frequencies consistent with the dominance of mirror modes appearing in δB ⊥,xy . The burst of left and right hand circularly polarized waves can be seen concentrated around frequencies ω ≈ 0.2ω c,i and ω ≈ −0.15ω c,i , respectively. Their range in wavenumbers is 0.2 ≲ ck ∥ /ω p,i ≲ 0.5. Overall, the power spectra of both left and right hand polarized waves are very similar to those of ion-cyclotron and electron cyclotron whistlers, and we will identify these waves as such from now on. In the next section, we will confirm that the population of particles that excites these waves have anisotropic distributions that are IC and whistler unstable.
The morphology of IC and whistler waves can also be seen in figures 2d and 2f . The short wavelength, wavepacket-like structures are identified with whistler modes, which propagate mainly through regions of low magnetic field strength of mirror modes, as we can see from δB ⊥,xy ( blue shaded regions in fig. 2d). The IC modes, on the other hand, are identified as the longer wavelength, extended modes that can be seen in δB z . The IC modes seem to propagate through the entire simulation box, given their ion-scale wavelength, whereas whistler modes clearly propagate within mirrors'  Figure 6. Panel a: The power spectrum of δBz(ω) + iδB ⊥,xy (ω) as a function of frequency. The frequencies are normalized by the initial ion-cyclotron frequency. The power spectrum of lefthanded circularly polarized waves (ω > 0) is shown as a solid blue line, whereas the power spectrum corresponding to righthanded circularly polarized waves (ω < 0) is shown as an orange line folded into positive frequencies. Panel b: Spectrogram of δBz(ω) + iδB ⊥,xy (ω) in frequency and time, at 0.4 < t · s < 1.3. The frequency is normalized by the initial ion-cyclotron frequency. Positive and negatives frequencies corresponds to left-hand and right-hand circularly polarized waves, respectively. Panel c: The power spectrum of δBz(ω, k ∥ ) + iδB ⊥ (ω, k ∥ ) at 0.5 < t · s < 1.1. Frequencies are normalized by the initial ion gyrofequency, and wavenumbers are normalized by the initial ion skin depth. Here also, positive and negative frequencies show left-hand and righthand polarized waves, respectively. Figure 7. The power spectrum in space of δB ⊥,xy (kx, ky) at t · s = 0.9. The wavenumbers kx, ky are normalized by the initial ion Larmor radius R init L,i . The solid and dashed white lines represent, respectively, the direction parallel and perpendicular to the main magnetic field at that time. magnetic troughs. This also resembles magnetosheath's observations of whistler waves within magnetic troughs (e.g. Kitamura et al. (2020)).
The peak frequencies observed in figure 6 for both ioncyclotron and whistler waves can be understood in terms of their dispersion relations. At high-β and kR L,e ∼ 1, and for quasi-parallel propagation, the dispersion relation for whistler waves can be written as (Stix (1992); Drake et al. (2021)) where d e = c/ω p,e and d i = c/ω p,i are the electron and ion skin depths, respectively. Knowing that d 2 i = R 2 L,i /β i , we can also write Similarly, at high-β and kR L,i ∼ 1, and for quasi-parallel propagation, the ion-cyclotron wave dispersion relation is approximately (Stix (1992)) and we can also write We can estimate k W , k IC by looking at the power spectrum of any of the perpendicular components of the magnetic field fluctuations. Figure 7 shows the power spectrum of δB ⊥,xy (k x , k y ) at t · s = 0.9, where the solid and dashed white lines denote the direction parallel and perpendicular to the mean magnetic field B at that time, respectively. Apart from the power in the perpendicular direction corresponding to the mirror modes, in the power parallel to B (i.e. along the solid black line in fig. 7) we can distinguish large wavenumbers centered at (k y R init L,i , k x R init L,i ) ≈ (0.75, −1.5) (and also at (−1.5, 0.75)), corresponding to whistlers, and also smaller wavenumbers centered at (k x R init L,i , k y R init L,i ) ≈ (0.5, 0.7), corresponding to ion-cyclotron waves.
The large wavenumber extent in k x , k y observed in fig.  7 gives us an approximate range of wavenumbers 1.5 ≲ k W R init L,i ≲ 3.2 for whistlers, implying frequencies 0 , consistent with the frequencies observed in the negative half of fig. 6c, corresponding to right-hand polarized waves. Similarly, the small wavenumber extent in k x , k y gives us a range of wavenumbers 0.4 ≲ k W R init L,i ≲ 1.1, implying frequencies 0.1 ≲ ω IC /ω init c,i ≲ 0.25, also consistent with the frequencies in the positive half of fig. 6c, corresponding to left-hand polarized waves.

2D Particle Distributions
The specific time at which ion and electron cyclotron wave activity saturates, which coincides with the end of mirror instability's secular growth (t · s ≈ 0.6), and the propagation of whistler waves within regions of low-magnetic field strength, give a hint towards uncovering the mechanism by which the whistler and IC waves are excited.
As a first step, we explore the evolution of the pressure anisotropy of ions and electrons at the time at which the IC and whistler waves are excited. At this time, mirror modes have achieved high amplitudes, and created sharp regions of high and low magnetic field strength, making the plasma spatially inhomogeneous. This implies that, in general, the plasma β of ions and electrons would not be the same at different locations in the simulation domain, making the anisotropy thresholds for the growth of the modes different in different regions. For this reason, a more appropriate method would be to measure the 2D distribution of pressure anisotropy, β ∥ and δB ∥ /B in the simulation domain. Figure 8 shows the distribution of ion and electron pressure anisotropy as a function of ion β ∥,i (panels a, b, c) and electron β ∥,e (panels g, h, i), respectively, and the distribution of δB ∥ /B versus β ∥,i (panels d, e, f ) and electron β ∥,e (panels j, k, l), respectively. These distributions are shown at three different times: beginning of the simulation (t · s ≈ 0, left column); end of mirror's secular growth and beginning of ion and electron cyclotron wave activity (t · s = 0.6, middle column), and a late stage well into the saturated regime of mirror instability (t · s = 1.4, right column). In the top row of fig. 8 (i.e. panels a, b, and c), the dashed gray line corresponds to the approximate mirror instability thresh-old 1/β ∥,i (Hasegawa (1969)), the dashed-dotted orange line corresponds to the theoretical IC threshold 0.53/β 0.4 ∥,i from Gary & Lee (1994) for γ IC /ω c,i = 10 −2 , and the solid black line is the best-fit to the global ion anisotropy derived in section 3.1 (see fig. 3a). In the third row of fig. 8 (panels g, h, i), the dotted-dashed black line shows the whistler instability threshold 0.36/β 0.55 ∥,e from Gary & Wang (1996), for γ W /ω c,e = 10 −2 .
Starting with the ions, we can see that, from a stable, isotropic distribution at the very beginning of the simulation ( fig. 8a), the ions become anisotropic enough to surpass both the mirror and the theoretical IC threshold from Gary & Lee (1994), as well as our best-fit instability threshold, as shown in fig. 8b. At this point (t · s = 0.6), we start to observe the excitation of ion-cyclotron waves that seem to interact with the ions and start driving them towards a marginally stable state. This can be seen in fig. 8c, where the distribution becomes bimodal, with one population of ions under both the IC-threshold and our best-fit threshold (centered at β ∥,i ∼ 5 and P ⊥,i /P ∥,i ∼ 1.2), meaning that they are driven towards marginal stability with respect to the IC threshold. Interestingly, there exists another ion population that is still unstable (centered at β ∥,i ∼ 18 and P ⊥,i /P ∥,i ∼ 1.4), therefore IC waves could then continue being excited even at this late stages. This could explain the sustained amplitude observed in δB 2 z and δB 2 ⊥,xy in figure 2g. Therefore, we can see that the unstable population has a higher β ∥,i , and the marginally stable population moves to lower β ∥,i .
For a similar value of P ∥,i , the difference in the values of β ∥,i between the unstable and marginally stable populations should imply a difference in the local magnetic field strength (recall β ∥,i = 8πP ∥,i /B 2 ). This gives us a hint on the location of the unstable and marginally stable populations in the domain, as mirror modes generate distinct regions of low and high magnetic field strength.
As we can see in figs. 8d, 8e, and 8f , the ions also separate into two populations now in δB ∥ /B. Starting from zero magnetic field fluctuations at the beginning (t · s ≈ 0, fig.  8d), we see how δB ∥ /B starts to grow at t · s = 0.6 ( fig. 8e), until we clearly see the bimodal distribution at t · s = 1.4, separating the two ion populations: the high-β ∥,i population located in regions of δB ∥ /B < 0 (i.e. low-B strength), and the low-β ∥,i population located in regions of δB ∥ /B > 0 (i.e. high-B strength).
We can therefore conclude that, after mirror modes develop and the IC waves are excited (t · s ≳ 0.6), the ions separate into two populations: one of low-β ∥,i , located mainly in high-B strength regions, and marginally stable to IC waves, and the second population with high-β ∥,i , low-B strength regions, and still unstable to IC waves. This suggests that the IC wave are excited by the unstable ion populations in regions of low magnetic field strength, and then interact with the ions in such a way that the ions move to regions of high-B strength and low β ∥,i . In sections 3.5 and 3.6 we will see that the population of ions that contribute most to the , t · s = 0.6 (middle column), and t · s = 1.4 (right column). The dashed gray line represents the approximate mirror instability threshold 1/β ∥,i (Hasegawa (1969)), the dotted-dashed orange line represents the IC instability threshold from Gary & Lee (1994) for γIC /ωc,i = 10 −2 (γIC is the IC instability growth rate), and the solid black line represents our best-fit threshold from section 3.1 (see fig. 3a). Second row: The distribution of δB ∥ /B versus ion β ∥,i for the same three times as in the top row. Third row: The distribution of electron P ⊥,e /P ∥,e versus β ∥,e in the simulation domain at the same three times as in the top row. The dotted-dashed black line represents the whistler instability threshold from Gary & Wang (1996). Fourth row: The distribution of δB ∥ /B versus electron β ∥,e for the same three times as in the top row. An animated version of this plot is available in the online version.
anisotropy that destabilize the IC waves are the ones that become trapped within mirror troughs. In the case of the electrons, we can see a similar evolution. From a stable, isotropic distribution at t · s ≈ 0 ( fig.  8d), we can see how part of it becomes now whistler unstable at t · s = 0.6 ( fig. 8e), after which the excited whistler waves interact with the electrons driving again part of the distribution gradually towards marginal stability, also generating a bimodal distribution similar to that of the ions. At t · s = 1.4 ( fig. 8f ), we can see that the electron population with low β ∥,e (centered at β ∥,e ∼ 5 and P ⊥,e /P ∥,e ∼ 1) is marginally stable with respect to the whistler threshold, whereas the electron population with high β ∥,e (centered at β ∥,e ∼ 18 and P ⊥,e /P ∥,e ∼ 1.2) is still unstable with respect to the whistler threshold. This also implies that whistler waves can still be excited at late stages in the simulation.
Analogously, the electrons also separate into two populations with respect to δB ∥ /B. Similarly to ions, we also see that the population with low-β ∥,e is located in regions of δB ∥ /B < 0 (low B strength), whereas the high-β ∥,e population is located in regions of δB ∥ /B > 0 (high B strength). In this sense, we also conclude that in the case of electrons, the unstable population is located mainly in regions of low-B strength and high-β ∥,e , where whistler waves are being excited, and the marginally stable population is located mainly in regions of high-B field and low-β ∥,e . This also suggests that whistler waves interact with electrons so they move to regions of high-B strength. We will also see in sections 3.5 and 3.6 that the electrons that contributes the most to the pressure anisotropy that destabilizes whistler waves are the ones that become trapped within mirror modes.

Physical Mechanism of Secondary IC/Whistler Excitation: Trapped and Passing Particles
In this section, we study the evolution of the ions and electrons that become trapped within mirror modes as part of the mirror instability's interaction with the particles. We characterize the pressure anisotropy and distribution functions of these populations at the moment of trapping, and provide evidence that they are able to destabilize parallel propagating modes that ultimately allow them to escape the mirrors and regulate the overall anisotropy.
As part of their evolution, and after reaching secular growth, mirror modes start to trap particles of low parallel momentum p ∥,j (j = i, e) in regions of low local magnetic field strength. The trapped particles bounce between these regions and conserve their magnetic moment in the process (Southwood & Kivelson (1993); Kunz et al. (2014)). In order to investigate the relation between this trapping process and the excitation of the these late IC and whistler waves, we select and track a population of ions and electrons throughout the evolution of the simulation, and study the trapped and passing (i.e. untrapped) subpopulations separately.
We select and track two populations of ions and two populations of electron having relatively small and large parallel momentum at a particular time in the simulation. This way, we make sure that we can capture particles that eventually be- p ,e /p ,e0 come trapped and others that remained passing. In our fiducial simulation b20m8w800, the two populations of ions that we track have parallel momentum −0.12 < p ∥,i /m i c < 0.12 and 0.3395 < p ∥,i /m i c < 0.3405 at t·s = 0.4. Similarly, the two populations of electrons have −0.2 < p ∥,e /m e c < 0.2 and 0.4599 < p ∥,i /m i c < 0.4601 at t · s = 0.4.
In order to study the behavior of the tracked particles when the IC and whistler activity starts, we ask how many particles become trapped and how many become passing during the interval of time at which this activity happens, which we denote by ∆τ LR . To answer this, we look at fig. 2g and define ∆τ LR as the interval of time 0.52 < t · s < 0.62, which covers the exponential growth that δB 2 z and δB 2 ⊥,xy undergo before saturating. This interval of time also covers the majority of the secular growth of mirror modes (see δB 2 ∥ ). Having this time interval well defined, we now must define the criterion by which we consider a particle to become trapped and passing during ∆τ LR , and for this we look at the evolution of their parallel momentum. Similarly to Ley et al. (2023), we define a particle as trapped during ∆τ LR if the median of its parallel momentum over ∆τ LR is smaller than one standard deviation over ∆τ LR . We then define a particle as passing if the median of its parallel momentum over ∆τ LR is greater than or equal than one standard deviation over ∆τ LR . This is a statement of small variation of p ∥,j over ∆τ LR , which in turn is a proxy for an oscillatory be-havior of p ∥,j , characteristic of a bouncing particle between mirror points. We confirm that this simple criterion gives excellent results separating trapped from passing particles. Figure 9 shows the evolution of the parallel momentum of a trapped and a passing ion (panels a) and a trapped and a passing electron (panels b), where the dashed vertical gray lines indicate ∆τ LR . We can see the the oscillation pattern in the evolution of the parallel momentum of the trapped ion during ∆τ LR and until t · s ≈ 0.7, when it escapes. The parallel momentum of the passing ion evolves without major changes as the ion streams through the simulation box. This behavior is consistent with previous works using hybrid and fully kinetic simulations Kunz et al. (2014); Riquelme et al. (2016).
In figure 9d we can also see the oscillating pattern of the parallel momentum of the trapped electron, indicating bouncing inside mirror modes, which ends at t · s ≈ 1.1, when it escapes. The parallel momentum of the passing electron does not vary significantly during ∆τ LR , confirming that it was streaming along field lines at least at that interval.
It is worth noting, however, what happens after ∆τ LR . Our criterion for identifying particles as trapped and passing was only within ∆τ LR , and after that period of time particles can continue evolving into the saturated stage of mirror modes, where they can escape, be trapped again or continue streaming unperturbed. Indeed, by looking at its parallel momentum, we can see that after escaping and streaming for a while, the trapped ion shown in figure 9a gets trapped again at t·s ≈ 1.1, bounces inside a mirror mode and escapes again at t · s ≈ 1.4. Similarly, we can also see that the trapped electron shown in figure 9b gets trapped again at t · s ≈ 1.2 and seems to stay trapped until the end of the simulation. Interestingly, the passing electron also gets trapped at around t · s ≈ 0.7, by looking at its parallel momentum, and then escapes again at t · s ≈ 1.2. Therefore, in a statistical sense, we can consider the particles as trapped and passing only over the particular period of time ∆τ LR that we chose, after which they can continue evolving and turn into passing or trapped again, as long as the mirror saturation persists in the simulation.

Physical Mechanism of Secondary IC/Whistler Excitation: Distribution Functions
In this section, we look at the evolution of the pressure anisotropy and distribution functions of trapped and passing ions and electrons defined according to the criterion described in section §3.5. We see that during ∆τ LR , both trapped ions and trapped electrons contribute most of the pressure anisotropy necessary to destabilize IC and whistler modes. We show that these IC and whistler waves interact in a quasilinear fashion with ions and electrons, respectively, and quickly regulate their pressure anisotropy such that their distributions evolve to a more isotropic state. Figure 10a shows the evolution of the pressure anisotropy of trapped and passing ions. We can see that the anisotropy of trapped ions initially follows a double-adiabatic (CGL, dotted blue line) evolution until t · s ≈ 0.5 (i.e. just start- ing ∆τ LR ), when the mirror modes start to trap them. We can readily see that during ∆τ LR , the trapped ions develop a significant anisotropy, peaking at around t · s ≈ 0.55. The anisotropy is quickly regulated and converges to the best-fit threshold that we derived in section 3.1 and show in figure 3a. Similarly, the pressure anisotropy of passing ions evolves in a relatively unperturbed fashion following CGL evolution (dotted red line) through the majority of ∆τ LR , until t · s ≈ 0.6, where it passes from negative values (consistent with passing ions having preferentially large parallel momentum) to a positive but, more isotropic state consistent with the best-fit threshold from fig. 3a.
The behavior of the pressure anisotropy of trapped and passing particles can be understood as follows. Mirror modes interact resonantly with ions and electrons according to the resonance condition ω M − k ∥,M v ∥ = 0, where ω M and k ∥,M are the frequency and parallel wavenumber of mirror modes, respectively, and v ∥ is the parallel velocity of the particle. The very low frequency of mirror modes, ω M ∼ 0, implies that the resonant particles are the ones having very low v ∥ (v ∥ < γ M /k ∥,M , where γ M is the mirror growth rate, Southwood & Kivelson (1993); Pokhotelov et al. (2002)). These are the particles that become trapped within mirror modes (Kivelson & Southwood (1996)). Consequently, all trapped particles have very low parallel velocity and, as a whole, they should naturally have a pressure anisotropy P ⊥,j > P ∥,j (j = i, e). Similarly, all passing particles have large v ∥ , and therefore they have a pressure anisotropy P ∥,j > P ⊥,j . In this sense, fig. 10 is consistent with the trapping argument described in Kivelson & Southwood (1996) (see their fig. 1).
The fact that both trapped and passing ions evolve into the average level of ion anisotropy shown in fig 3a shows that their trapped or passing condition corresponds to a transient state, that passes after a time comparable to ∆τ LR . Also, notice that the anisotropy of the two populations (and for the whole population for that matter) is significant enough to drive IC waves unstable (see section 3.3), and therefore this can provide evidence for the source of the IC waves that we see. If this is the case, their interaction with ions is the source of the quick regulation of the anisotropy that we see in fig. 10a. Interestingly, under this scenario, the regulation of the pressure anisotropy of passing ions, which happens at the same time as that of the trapped ions, should also be due to the interaction with these IC waves, meaning that the IC waves interact with both populations of trapped and passing ions simultaneously, and therefore regulate the global ion anisotropy. We confirm that this is the case by looking at the evolution of the distribution functions of trapped and passing ions.
In the case of electrons, we observe a similar evolution in figure 10b. Initially, both trapped and passing electrons detach from their respective CGL evolution (dotted blue and red lines, respectively), and develop a significant anisotropy ∆P e > 0, that peaks at t · s ≈ 0.4. We also see that trapped electrons detach from their CGL evolution much earlier than passing electrons. This evolution then leads to the early burst of whistler waves, which also quickly regulates and drives anisotropies of both trapped and passing electrons towards a more isotropic state (see section 3.2). As expected, the anisotropy of trapped electrons is higher than the one of the passing electrons. After this process, and during ∆τ LR , the anisotropy of trapped electrons increases again, while that of passing electrons continues to decrease. This way, we see that trapped electrons build up a pressure anisotropy ∆P e > 0 that is also quickly regulated after ∆τ LR , converging to an anisotropy level similar to the one of the general electron populations. The anisotropy ∆P e < 0 of the passing electrons also gets regulated towards a similar anisotropy level during the same time. This evolution of trapped electrons also suggests that they become anisotropic enough to destabilize whistler waves, and therefore could be the source of the whistler activity observed at t · s > 0.6. We provide evidence of this by showing the evolution of the distribution function of electrons. Figure 11 shows the distribution functions of trapped and passing ions and electrons at three different times t·s = 0.57, t · s = 0.61, and t · s = 0.75, spanning ∆τ LR and also part of mirror's saturated stage. In the following we describe the evolution of each population: The distribution of trapped ions (figs. 11a, 11b, and 11c) shows a clear loss-cone like form at t · s = 0.57 (all outside the loss-cone), meaning that all trapped ions are effectively trapped in mirror troughs. At this time, trapped ions have reached its maximum pressure anisotropy according to figure 10a.
Once IC waves are excited, they interact with both trapped and passing ions via pitch-angle scattering in a quasilinear fashion (Kennel & Engelmann (1966)). This diffusion process happens along paths of constant particle's energy in the frame moving with the waves (see e.g. Squire et al. (2022) We plot these contours in solid white lines in each plot of figure 11 as v 2 ⊥,j + (v ∥,j − ω/k ∥ ) 2 ≈ v 2 ⊥,j + v 2 ∥,j = const., as in a high-β scenario, the phase velocity of an IC wave offers a small correction of order v A /v th,i = 1/β. Additionally, the IC waves in our simulations are destabilized in both parallel and anti-parallel directions to B. We can see that the relaxation of the distribution function of trapped ions by the quasi-linear interaction with IC waves agrees very well with these paths, by looking at t · s = 0.61 and t · s = 0.75.
The distribution of passing ions (figs. 11d, 11e, and 11f ) shows, on the one hand, a concentration of ions at low perpendicular velocities and relatively large parallel velocities, and it looks fairly symmetric in v ∥ . This is consistent with having untrapped ions mainly streaming along the mean magnetic field in both directions. On the other hand, the population of large parallel velocity is also shown at v ∥ /c ≈ 0.3 (see section 3.5). Interestingly, the passing ions also interact quasilinearly with IC waves, and this is particularly evident in the evolution of passing ions. Indeed, we can clearly see how the large parallel velocity population of passing ions evolves along the contours of of constant particle energy with  Figure 11. The distribution function f (v ∥,j , v ⊥,j ) of trapped and passing ions and electrons at three different times: t · s = 0.57 (first column), t · s = 0.61 (second column), and t · s = 0.75 (third column). The distribution function ftrapped(v ∥,i , v ⊥,i ) of the trapped ions is shown in the first row, fpassing(v ∥,i , v ⊥,i ) for the passing ions are shown in the second row, ftrapped(v ∥,e , v ⊥,e ) for the trapped electrons are shown in the third row, and fpassing(v ∥,e , v ⊥,e ) for the passing electrons are shown in the fourth row. In all the plots, the solid white curves denote contours of constant particle energy in the frame moving with the waves: . An animation is available. excellent agreement at t · s = 0.61 and t · s = 0.75. We can understand the evolution of this population by looking at the gyroresonance condition If we look at the peak power at positive frequencies in the power spectrum shown in fig. 6c, we can estimate the frequency and wavenumber at which most of the power of IC waves resides: ω/ω init c,i ≈ 0.2, and ck ∥ /ω init p,i ≈ ±0.15. From eq. (6) we can estimate then the parallel velocity of the ions interacting gyroresonantly with these IC waves: which gives v ∥,i /c ≈ 0.36 and v ∥ /c ≈ −0.24, which falls in the range of the large parallel velocity population. The quasilinear evolution also happens with the population with smaller parallel velocity.
The population of trapped electrons (figs. 11g, 11h, and 11i) shows a very similar evolution to that of trapped ions; the loss-cone like distribution is also apparent. The evolution of this distribution is also consistent with a quasilinear interaction now between the electron and whistler waves, driving the distribution towards isotropy along paths of constant particle energy, as can be seen at later times in figure 11.
Finally, the population of passing electrons (figs 11j, 11k, and 11l) also shows a very similar evolution to that of the ions. The populated loss-cone shape of the distribution is also apparent, and we can see the quasilinear evolution of the distribution function along constant particle energy contours at later times.
This way, we have provided evidence for the source of both IC and whistler waves observed in our simulations. Once ions and electrons get trapped in regions of low magnetic field strength of mirror modes, they get significantly anisotropic with a loss-cone like distribution, which is able to destabilize parallel-propagating IC and whistler waves, respectively. These waves then interact with both population of trapped and passing particles in a quasilinear fashion, driving both populations of trapped and passing ions and electrons towards a more isotropic state. Consequently, this mechanism can contribute to regulate the global anisotropy of ions and electrons, and can thus be a pathway for particle escape and consequent saturation of mirror modes (Kunz et al. (2014)).

MASS-RATIO DEPENDENCE
In this section, we compare simulations with different mass ratios: m i /m e = 8, m i /m e = 32, but with the same initial conditions for ions, as shown for runs b20m8w800, b20m32w800,and b20m64w800 in Table 1, although with somewhat different temperatures. We see that IC and whistler waves' signatures do appear in all three simulations, and thus they do not seem to present a strong dependence on mass ratio. Figure 12 shows the evolution of δB 2 ∥ (panel a) and δB 2 z (panel b) for the three runs with mass ratios: m i /m e = 8, 32, and 64 (runs b20m8w800, b20m32w800, and b20m64w800 in table 1). We can see a very consistent evolution of δB 2 ∥ in all three runs, meaning that m i /m e does not play a significant role on the early evolution and saturation of the mirror instability. Similarly, δB 2 z shows the same features in all three runs, especially during mirrors' secular growth and saturated stages (t · s ≈ 0.5 onwards). The early peak in δB 2 ∥ at t · s ≈ 0.4 corresponding to the early whistler burst is also seen in the three runs, but more prominently in the simulation with m i /m e = 8. This is possibly due to an enhancement of this wave activity by the ions, which are able to weakly feel the presence of whistlers, as the mass separation is not very large. This effect disappears as the mass ratio increases, and the early whistlers only affect the electrons. More importantly, for t · s > 0.5, all three runs show a very similar evolution of δB 2 ∥ . Figure 13 shows the evolution of the pressure anisotropy of ions (panel a) and electrons (panel b) for the same three runs. In the case of the ions, we can see an overall evolution that is very consistent in all three runs, both in early and late stages. We can see a smaller anisotropy overshoot for the simulation with m i /m e = 8 at t · s ≈ 0.4, coincident with the enhancement seen in δB 2 z , during the early whistler burst, suggesting that ions can weakly interact with the whistlers at this mass ratio, and consequently their anisotropy does not reach the same overshoot as the rest of the runs. Notwithstanding the foregoing, we can see how all three runs display a very similar pressure anisotropy evolution afterwards, which is also well described by the best-fit threshold ∆P i ∝ β −0.45 i shown in fig. 3.
In the case of the electron pressure anisotropy ∆P e , we can also see a similar evolution overall in fig. 13b. The overshoot at t·s ≈ 0.4 is larger for decreasing mass ratios, possibly due to the fact that the whistler amplitude required for efficient scattering decreases as m i /m e increases, as explained above. This means that, after ∆P e /P e,∥ has surpassed the threshold for efficient growth of the whistler modes, the simulations with larger m i /m e take shorter times to reach the necessary whistler amplitude to efficiently scatter the electrons. This implies that the overshoot decreases for higher mass ratios. During late stages, we can see a very similar evolution of ∆P e in all three runs, that is even more evident for m i /m e = 32 and m i /m e = 64 (orange and green curves in fig. 13a), which essentially lie on top of each other.
Here we also see a very similar power distribution at both mass ratios, showing both left-hand and right-hand polarized waves (positive and negative frequencies, respectively). The peak power is also observed at the same frequencies and wavenumbers as in fig. 6 for both polarizations.  Figure 14. The power spectrum of δBz(ω, k ∥ ) + iδB ⊥ (ω, k ∥ ) at 0.5 < t · s < 0.7 for mi/me = 32 (run b20m32w800, left panel) and mi/me = 64 (run b20m64w800, right panel). Positive and negatives frequencies show the power in left-hand and right-hand polarized waves, respectively. This way, we can see that the linear and nonlinear evolution of the mirror instability and the late IC and whistler evolution are well captured in our simulations, and it does not strongly depend on mass ratio.

DEPENDENCE ON INITIAL PLASMA β
We tested whether the IC and whistler waves' activity is present in simulations with β init i = 2 (i.e, total β init = 4), and β init i = 40 (i.e. total β init = 80), and compare them with our fiducial simulation at β init i = 20. We confirm that the mirror instability can develop in all simulations, and both IC and whistler waves do appear at nonlinear stages.
The power spectrum of δB z (ω, k ∥ ) + iδB ⊥,xy (ω, k ∥ ) is shown in figure 15, and we can see that it is similar among the three β i cases. In all three cases we see the power concentrated at ω ∼ 0 corresponding to mirror modes. In addition, we also see a concentration of power in right and left polarized waves, so both IC and whistler waves are also present, although their peak frequency changes. For the β init i = 2 case we see that the peak frequency is at ω/ω init c,i ≈ 0.5, whereas in the β init i = 40 case it shifts to smaller values, ω/ω init c,i ≈ 0.1. This shift in peak frequency can also be explained by the IC and whistler dispersion relations analogous to our discussion in section 3.3. Figure 16 compares the evolution of δB 2 ∥ (i.e., mainly the development of the mirror instability) for the three runs with different initial β init (the other phyiscal parameters are the same, see table 1). In all three cases we can see an exponential phase followed by the secular and saturated stages characteristic of the mirror instability, which develops earlier for higher initial β init , consistent with the smaller anisotropy threshold for the growth of the mirror instability at larger beta. The amplitude of δB 2 ∥ at the saturated stage is comparable for both β init = 20 and β init = 40 runs, and is smaller for the β init = 2 run, as also seen by previous works (e.g. Riquelme et al. (2015)).
Indeed, when we look at the evolution of δB 2 z , we can see that for both β init = 20 and β init = 40 runs, the evolution is similar: both display an early whistler burst at t·s ≈ 0.4, and a IC/whistler excitation stage (t · s ≈ 0.5 onwards) at almost the same amplitude. In the case of the β init = 2 run, we can see that the first exponential growth in δB 2 z at t · s ≈ 0.6 is consistent with an IC burst (see e.g. Ley et al. (2019)), after which we see the typical oscillation pattern that the excitation of late IC and whistler waves produces, from t · s ≈ 0.8 onwards, saturating at a similar amplitude than the rest of the runs, and displaying a very high-frequency oscillation.
In figure 17, we compare the evolution of the ion and electron pressure anisotropy plotted as a function of their parallel plasma β i for the three simulations with different initial β i (As in all our simulations the mean magnetic field strength is continuously increasing, so the particles' β i decreases over time, and therefore the simulation evolves towards the left in fig. 17.).
In the case of the ions ( fig. 17a), we can see a similar overshoot and subsequent regulation, but the overshoot occurs at a lower anisotropy value for increasing β i . This is consistent with the inverse β i dependence of the mirror instability threshold: mirror modes are excited earlier at higher β i , and therefore have relatively more time to regulate the anisotropy before it reaches a higher overshoot. Interestingly, the saturated stage of the ion pressure anisotropy is consistent with the theoretical IC threshold from Gary & Lee (1994): for γ IC /ω c,i = 10 −2 (see fig. 3a) in all three runs, suggesting a universality in the threshold that ∆P i /P ∥,i follows, as a consequence of the excitation of IC waves during mirrors' saturated stage. (In the case of the β init i = 40 run, however, it is more unclear whether it can follow the above mentioned threshold at late stages, given the short duration of this run.) In the case of electrons ( fig. 17b), we can also see that the overshoot is reached at lower values of the pressure anisotropy ∆P e /P ∥,e for increasing initial beta, consistent with an inverse-β i dependence now of the whistler instability anisotropy threshold. It is interesting to note that after the anisotropy overshoot, and during these late stages, the electron pressure anisotropy tends to be significantly smaller than the expectation from the threshold for the whistler instability in the higher initial β i runs (β init i = 20 and β init i = 40), irrespective of the generation of pressure anisotropy that the continuous amplification of the magnetic field produces as a consequence of the shear motion in the simulation. Notice, however, that in low magnetic field regions the electron pressure anisotropy is larger than the whistler threshold for growth rate γ = 0.01ωc,e, from Gary & Wang (1996). and, therefore, enough to excite whistlers (fig 8). This shows the key role played by mirror-generated magnetic troughs in creating the conditions to excite whistlers despite the fact that, globally, the pressure anisotropy may be not be enough to make these waves unstable. On the other hand, in the β init i = 2 run, ∆P e /P ∥,e continues to weakly grow because of the continuous B amplification, and this is done following a marginal stability state well described by the threshold of the whistler instability ∆P e /P ∥,e ∝ β −0.55 (Gary & Wang (1996)), consistent with previous works at lower β ∥,e (Ahmadi et al. (2018)).
The persistence of the late IC and whistler activity at different initial plasma β i suggests that this phenomenon is a natural consequence of the excitation of the mirror instability. In other words, in a weakly collisional plasma with an initial plasma β i sufficiently high to effectively excite the mirror instability, the excitation of IC and whistler waves at its late, saturated stages seems to be ubiquitous. ror instability, and provide an interesting physical connection between ion-scale instabilities and electron-scale physics.
In this work, we did not vary the scale-separation ratio ω c,i /s. In an environment like the ICM, turbulent eddies could drive the plasma locally through shear motions at kinetic scales with a wide range of frequencies s, and we typically expect larger kinetic energy at low frequencies (i.e., higher ω c,i /s). For larger values of ω c,i /s, previous works have shown that mirror modes can develop comparatively earlier in the simulations, therefore having relatively more time to saturate, and reaching similar amplitudes (Kunz et al. (2014); Melville et al. (2016); Riquelme et al. (2016); Ley et al. (2023)). In this sense, we would expect a similar late excitation of IC and whistler waves once mirror modes have reached a saturated stage.
The excitation of IC and whistler waves at saturated stages of the mirror instability modulates its nonlinear evolution, and therefore could affect transport processes in the ICM in which mirror modes come into play.
Particularly important is the pressure anisotropy regulation in the context of collisionless heating and dissipation via magnetic pumping in the ICM (Kunz et al. (2011);Ley et al. (2023)). The marginal stability level that the ion pressure anisotropy reaches at the saturated stage, ∆P i ∝ β 0.45 ∥,i (see fig. 3a, also correctly pointed out by Sironi & Narayan (2015)) is larger than the usual mirror threshold 1/β ∥,i by a factor ∼ β 0.55 ∥,i . which directly translates into an excess heating of the same order. Indeed, given that β is estimated to be β ∼ 10 − 100, and that the heating rate is directly proportional to the pressure anisotropy, this could imply a heating rate several times larger than predicted from the mirror threshold, enhancing the efficiency of the mechanism by draining more energy from the turbulent motions that drive the pumping.
The structures of high and low magnetic field that mirror modes produce in the saturated stage seem to be persistent in time, and its energy δB 2 ∥ does not decrease as long as the amplification of the mean magnetic field B is maintained (see fig. 2g). Even when this amplification is halted or reversed, the decaying timescales of mirror modes are large compared to the typical ion gyroperiod (Melville et al. (2016); Ley et al. (2023)). This implies that the trapping process of ions and electrons also persists, along with the excitation of secondary IC and whistlers. This source of whistler waves can have interesting implications in the context of ICM thermal conduction models like whistler-regulated MHD (Drake et al. (2021)), as they can dominate the electron scattering in the presence of mirror modes.
This source of whistler waves associated to mirror modes can also influence the suppression of the effective heat conductivity in the plasma even in the absence of heat-fluxes (Komarov et al. (2016); Riquelme et al. (2016); Roberg-Clark et al. (2016), and this can have consequences in larger-scale instabilities such as the Magneto-thermal instability (MTI, Balbus (2000); Berlok et al. (2021); Perrone & Latter (2022a,b)).
Future work aimed towards 3D fully kinetic PIC simulations would be required to have a full understanding of the consequences of the mirror instability and secondary IC/whistler excitation in these high-β plasmas.
We thank Aaron Tran for providing the dispersion solver used in this work, and we thank Lorenzo Sironi, Jonathan Squire and Alexander Schekochihin for useful comments and discussion. F.L. acknowledges support from NSF Grant PHY-2010189. M.R. thanks support from ANID Fondecyt Regular grant No. 119167. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant No. ACI-1548562. This work used the XSEDE supercomputer Stampede2 at the Texas Advanced Computer Center (TACC) through allocation TG-AST190019 (Towns et al. (2014)). This research was performed using the compute resources and assistance of the UW-Madison Center For High Throughput Computing (CHTC) in the Department of Computer Sciences. This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).