Making Waves: Mirror Mode Structures Around Mars Observed by the MAVEN Spacecraft

Abstract We present an in‐depth analysis of a time interval when quasi‐linear mirror mode structures were detected by magnetic field and plasma measurements as observed by the NASA/Mars Atmosphere and Volatile EvolutioN spacecraft. We employ ion and electron spectrometers in tandem to support the magnetic field measurements and confirm that the signatures are indeed mirror modes. Wedged against the magnetic pile‐up boundary, the low‐frequency signatures last on average ∼10 s with corresponding sizes of the order of 15–30 upstream solar wind proton thermal gyroradii, or 10–20 proton gyroradii in the immediate wake of the quasi‐perpendicular bow shock. Their peak‐to‐peak amplitudes are of the order of 30–35 nT with respect to the background field, and appear as a mixture of dips and peaks, suggesting that they may have been at different stages in their evolution. Situated in a marginally stable plasma with β ‖ ∼ 1, we hypothesize that these so‐called magnetic bottles, containing a relatively higher energy and denser ion population with respect to the background plasma, are formed upstream of the spacecraft behind the quasi‐perpendicular shock. These signatures are very reminiscent of magnetic bottles found at other unmagnetized objects such as Venus and comets, also interpreted as mirror modes. Our case study constitutes the first unmistakable identification and characterization of mirror modes at Mars from the joint points of view of magnetic field, electron and ion measurements. Up until now, the lack of high‐temporal resolution plasma measurements has prevented such an in‐depth study.


Introduction
Mirror mode-like (MM) structures have been found everywhere in solar system plasmas (Tsurutani, Lakhina, et al., 2011), from the solar wind (Bale et al., 2009;Kaufmann et al., 1970;Winterhalter et al., 1995) to Earth (Lucek et al., 1999b;Tsurutani et al., 1984), Mars (Bertucci et al., 2004;Espley et al., 2004) and Venus (Volwerk et al., 2008a(Volwerk et al., , 2008b, Jupiter (Erdős & Balogh, 1996) and Saturn (Violante et al., 1995), as well as comets (Glassmeier et al., 1993;Mazelle et al., 1991;Volwerk, Richter, et al., 2016). MM waves are low-frequency long-wavelength transverse waves, usually linearly polarised and non-propagating in the plasma rest frame. Their non-linear evolution has been discussed both for electrons and protons (Balikhin et al., 2010;Kivelson & Southwood, 1996;Soucek & Escoubet, 2011). An MM structure classically takes the spatial shape of a magnetic bottle imprisoning pockets of high-density plasma drifting with the ambient plasma. Thus, in the data, MMs commonly appear as sudden dips or peaks in the magnetic field intensity, anti-correlated with plasma density variations and with only little magnetic field angular variation across the structure (Tsurutani, Lakhina, et al., 2011). In competition with Alfvén ion cyclotron waves, with which they are co-generated in the plasma, MMs typically grow in a high-β plasma (β ‖ > 1) from an ion temperature anisotropy, itself triggered by any asymmetry upstream in the solar wind flow (Gary, 1992). The anisotropy either arises from the ring-beam distribution of locally new picked-up ions, as seen at Mars, Venus or comets, or as the downstream result of a quasi-perpendicular bow shock crossing, as seen on all objects, magnetized or unmagnetized (Tsurutani, Echer, et al., 2011). The mirror mode instability criterion (MMIC) is usually expressed as (Hasegawa, 1969): SIMON WEDLUND ET AL.
10.1029/2021JA029811 3 of 28 and/or temporal resolution. Our event took place during the early part of the MAVEN mission on 25 December 2014 around 11:30 UT, in the receding phase of solar cycle 24. MMs were found wedged against the MPB behind a quasi-perpendicular bow shock crossing-a classical empirical picture reminiscent of Earth-bound observations. After describing our detection methods based on magnetometer data only at 1 -Hz sampling frequency and validated against calculated ion and electron plasma moments (Section 2), we present an in-depth analysis of these MM structures in Section 3.2, providing recommendations for a new set of B-field-only detection criteria (Section 3.3). Finally, we discuss their possible origin in Section 4. One of the purposes of this study is to prepare for a statistical analysis of MM occurrences using the whole MAVEN data set.

Instrumentation
The MAG package consists of two tri-axial fluxgate magnetometers mounted at the extremity of two boomlets on MAVEN's solar panels . It measures 3-component magnetic fields with a nominal frequency of 32 Hz and an accuracy better than 0.05%. Such high temporal resolution is not necessary to investigate the global behavior of MM structures and hence a 1-s resolution is adopted throughout.
As part of the Particles and Fields package on board MAVEN, the Solar Wind Ion Analyzer (SWIA) is an electrostatic ion analyzer measuring ion differential fluxes with a maximum temporal resolution of 4 s (Halekas, Taylor, et al., 2015). Sensitivity is automatically switched to accommodate for a large dynamic range of measured ion fluxes (magnetosphere and solar wind). SWIA has thus two three-dimensional (3D) scanning modes, one coarse (360° × 90°, Δϕ = Δθ = 22.5° in anode/deflection angles, energy 25 eV/q to 25 keV/q , suitable for magnetosphere and pickup ions), one fine (45° × 45°, Δϕ = 4.5°, Δθ = 3.75°, 10% energy windows, suited to solar wind ions), both with 48 energy steps. Additionally, two separate 3D data packet telemetry modes are used on board which result in the Survey (lower cadence, maximum availability) and Archive (higher cadence but lower availability) modes. The combination of the scanning and telemetry modes results in four different data products, labeled in the following SWIFA (SWIA Fine Archive), SWIFS (Fine Survey), SWICA (Coarse Archive) and SWICS (Coarse Survey). The last two modes are thus those with lowest angular resolution and are usually selected in the magnetosheath.
Each of these four data products yields different moments at different temporal resolution. Onboard moments are automatically calculated either in coarse or fine modes, making use of the highest cadence available: in practice, they are usually provided at the highest possible temporal resolution of 4 s. Because SWIA does not discriminate between ion masses, the plasma composition needs to be assumed first (Halekas, Taylor, et al., 2015). Discrimination in the solar wind between H + and He 2+ ions can be done on an energy-per-charge basis. Depending on composition assumptions, errors of the order of √ are introduced on the onboard-calculated moments, M being the mass of the plasma (for example, M = 1 amu for protons) . The proton-dominated composition assumption is usually fulfilled in the not-so-dense parts of the Martian plasma environment, that is, the solar wind and magnetosheath down to the outer edge of the Magnetic Pile-up Boundary (MPB, sometimes also referred to as Induced Magnetosphere Boundary or IMB), which the present study concentrates on. In regions where heavy ions dominate (as in the ionosphere and during deep-dip campaigns), it is however necessary to call in supporting instruments with mass-resolving capabilities, such as the Suprathermal and Thermal Ion Composition analyser (STATIC) . In the following, ion moments of order 0 (density N i in cm −3 ), 1 (velocity vector V in km s −1 ) and 2 (pressure tensor ̄ in eV cm −3 and associated temperatures T in eV) are used. Calculations of these moments depend also on the energy thresholds and extreme care must be achieved in their derivation .
The Solar Wind Electron Analyzer (SWEA), located on a 1.5-m long boom, is a hemispheric electrostatic electron analyser designed to measure the energy and angular distribution of photo-and solar-wind electrons in the 3-4,600 eV energy range in 64 bins, with a resolution of ΔE/E ∼ 17% and a maximum cadence of 2 s (Mitchell et al., 2016). It swipes almost all 4π sr of solid angle with a 22.5° angular resolution in the azimuth direction and 20° along the direction of the elevation angle. As MAVEN is a 3-axis stabilised spacecraft, the field of view (FOV) is broadened up to 360° × 120° by the use of electrostatic deflectors theoretically covering 87% of the sky, but reduced to 79% due to spacecraft obstruction. Omnidirectional energy spectra (product 'SPEC') are calculated onboard at the highest temporal resolution by integrating the 3D distribution over angular sectors. Densities SIMON WEDLUND ET AL.
10.1029/2021JA029811 4 of 28 presented in the following are calculated from these onboard energy fluxes by integrating over the 64 energy bins, assuming isotropy over SWEA's blind spots. With its temporal resolution of ∼2 s, SWEA complements SWIA for total plasma density measurements when structures lasting less than about 8 s (2 × Δt SWIA ) are identified in magnetometer datasets, which is often the case for MMs. Because of quasi-neutrality, the electron density N e must theoretically be equal to the ion density, that is, N e ∼ N i . The absolute determination of the electron density depends on different experimental parameters (such as the product used) and on the accuracy in the determination of the spacecraft potential. The electron density from the SPEC product has the better available time resolution but its absolute value can differ from the moments computed on the ground from the 3D velocity distribution. Hence, only absolute variations in electron densities may be used in those cases to corroborate ion density variations.
With the help of MAG, SWEA and SWIA, several temporal and spatial scales are probed, incrementally accessing plasma features at 1 s and below, to 2 s, to 4-8 s, respectively.
The coordinate system adopted in the following for all vectorial and tensorial quantities such as B, ̄ or V, is that of Mars Solar Orbital coordinates (also known as Sun-state coordinates), abbreviated MSO: centred on Mars, the X axis points towards the Sun from the planet's centre, Z is in the direction of the North pole from the orbital plane, and Y completes the right-hand triad, so that the X-Y plane is Mars' orbital plane around the Sun. All quantities are thus expressed in MSO coordinates.

Detection Criteria: B-Field Only
Because of the ubiquity of magnetometers on board space missions (at the notable exception of Mars Express) and their usually high sampling frequency, mirror mode structures have traditionally been located by B-field only measurements (Joy et al., 2006;Lucek et al., 1999b;Soucek et al., 2008;Volwerk, Schmid, et al., 2016), relying on compressional and linearly polarised signatures. However, structures such as foreshock waves or fast mode waves linked to pickup ions are compressional in nature and may also fulfill these criteria. In those cases, plasma measurements are necessary to lift the ambiguity; one way is to check for the anticorrelation between density and magnetic field variations that is expected for MM structures (Hasegawa, 1969). Another is to evaluate, when possible, the MM-unstable condition with the help of the MMIC, assuming that the spacecraft crosses the source region of the MMs. In the following, B-field-only criteria are first used to pinpoint promising structures in the magnetosheath. They are then validated against plasma measurements (Section 2.3). Table 1 lists the MM detection criteria chosen, with references that inspired them. Following previous studies at Venus, only quasi-linearly polarised waves are sought. As in Volwerk, Schmid, et al. (2016), magnetic field measurements are first low-pass filtered using a 2-min-wide Butterworth filter (f band = 1/120 Hz, passband ripple 2 dB and stopband attenuation 20 dB) to approximate the total background field |B bg | = B bg . A compressional signature consistent with MMs is defined by large fluctuations around the background field, that is, Δ| |∕ bg = | | |bg |∕ bg , with threshold values typically ranging from 0.1 to 0.2. For example, Volwerk, Schmid, et al. (2016) used a threshold of Δ|B|/B bg ≥ 0.20. When dip and peak mirror structures are both present in the studied interval, Joy et al. (2006) advocated the use of B-field lower and upper quartiles to determine the 'true' background level: this may be an issue especially when long portions of sinusoidal-like oscillations lead to an under-or over-estimation of the ambient field intensity and hence a correspondingly smaller Δ|B|/B bg . We have checked that this issue does not significantly affect the results presented here by lowering the absolute fluctuation Δ|B|/B bg threshold from 0.20 to 0.15 in order to potentially capture more events.
A minimum variance analysis (Sonnerup & Scheible, 1998) (MVA) on the magnetic field directions is then performed with a 15-s moving window. The aim of MVA for MM detection is to constrain the wave polarisation of the detected structures (Lucek et al., 1999b). Ideally, linear structures are sought, adopting a cigar-shaped variance ellipsoid with maximum and intermediate eigenvalues λ max and λ int clearly distinct from one another (Génot et al., 2001;Tátrallyay & Erdös, 2002). In practice, non-linear structures have also been observed and studied theoretically (Califano et al., 2008), with elliptically polarised MMs found more often than stricto sensu linearly polarised ones (Génot, 2008;Génot et al., 2001). Consequently, criteria based on eigenvalue ratios should not be too constraining if one wants to capture candidate events straying somewhat from the pure linear wave polarisation (Sergeev et al., 2006), hence the "quasi-linear" nomenclature adopted here. In parallel, the quasi-linear nature of the sought-after waves can be reinforced by calculating the angle Θ maxV (Φ minV ) between the background field 5 of 28 direction and the maximum (minimum) variance directions (i.e., eigenvectors). Threshold values on these angles were adopted from Volwerk et al. (2008b) and Volwerk, Schmid, et al. (2016).
Quasi-linear MM wave candidates are finally automatically detected from B-field measurements only if the criteria listed in Table 1 are fulfilled: (a) compressional structures (Criterion 1), (b) quasi-linear wave polarisation (Criteria 2 − 4) and (c) sufficiently large average B-field intensities in a ±2 min interval around the MM candidate (magnetosheath conditions, Criterion 5). Values for each criteria have been determined empirically starting from theoretical considerations (Price et al., 1986), past observations at Earth and Venus, and tweaked to the MAVEN Mars data set to capture at least once most structures seen in a test interval. Criterion 5 (〈B bg 〉 ≥ 11 nT), inspired by a similar one from Dimmock et al. (2015) for observations at Earth, aims at removing potential outliers outside of typical magnetosheath B-field levels in an interval around the MM candidate.
Non-Gaussian statistics may additionally help with the classification of MM candidates . The propensity of the signal distribution in a given time interval (containing, for instance, m successive measurements of B) toward peaks or dips can be measured by its skewness, expressed as Skew where σ is the standard deviation of the module |B j | of B at each time step j. A positive skewness denotes the presence of peaks whereas a negative skewness points to the presence of dips (Criterion 6 in Table 1). Similarly, the excess kurtosis, defined as eKurt gives a measure of the fluctuations of the signal (presence of outliers in the tail of the probability distribution function): a negative kurtosis, significantly different from 0 (pure Gaussian statistics), implies a B-field signal with large oscillations, a characteristic that trains of MMs share . This is given in Criterion 6.
Finally, magnetic azimuth and elevation angles, defined as = arctan ( ∕ ) and = arctan can be calculated. The magnetic field is expected to rotate by less than about 15° for linear MMs (Treumann et al., 2004), a property that can be checked after the detection is performed.
Note that the eigenvalue criteria (Criteria 3 and 4) aim at reinforcing the quasi-linear polarised nature of the wave modes, ensuring that the maximum variance direction (the tangential component of the eigenvector triad) is well defined and that quasi-degeneracy of the covariance matrix is kept to the two minimum eigenvalues. Although the direction of the wave vector k is in general ill-determined, Criteria 4-5 imply that the variance ellipsoid is in  Note. The absolute B-field fluctuation is estimated by Δ|B|/B bg , with B bg the background intensity (low-pass Butterworth filter). A moving MVA analysis with a 15-s window is performed. Θ maxV (Φ minV ) is the angle between the maximum (minimum) variance direction and that of the background field. λ max /λ int is the ratio of the maximum and intermediate eigenvalues, whereas λ int /λ min is the ratio of the intermediate and minimum eigenvalues, as determined from the MVA analysis. Criterion 5, with average B-field 〈B bg 〉 performed on an interval ±2 min around the potential candidate MM, further checks for potential outliers that are tentatively in the solar wind. Criterion 6 classifies MMs into predominantly dips and peaks in the time interval considered (for example ±2 min around the detected candidate). Criterion 7 checks for changes of magnetic azimuth and elevation angles across an individual MM structure Δ(az, el). The initial values are those used for the first detection of the structures (Section 3.2), whereas the revised values were manually obtained after validation of the detections (Section 3.3). 10.1029/2021JA029811 6 of 28 general cigar-shaped, and elliptically polarised structures are mostly filtered out to the benefit of linearly polarised waves.

Validation With Plasma Measurements
In parallel, SWIA plasma densities N i for the available modes are calculated from the moments of order 0 of the ion distribution function. Moreover, the background-normalised amplitude ΔN/N = (N i − N bg )/N bg is calculated where N bg , the background plasma density, is obtained by applying a low-pass Butterworth filter to the SWIA N i densities, with the same filter parameters as for the magnetic field. In this way, ΔN/N has the advantage of changing signs for peaks or dips in the signal. For comparison, total B-field intensities are downsampled to SWIA's mode resolution and normalised amplitudes are calculated in the same way, that is, Δ ∕ bg = | | − bg ∕ bg . When ΔB/B and ΔN/N are of opposite signs at a candidate MM, they are in effect anticorrelated. Those candidates that fulfill the B-field-only criteria and the anticorrelation signature are considered MM signatures proper.
For further confirmation of the nature of the MM candidates, two separate criteria are further checked against: 1. Direction (also, Criterion 7 in Table 1): the magnetic field direction in a region containing MMs should not change by more than 15° on average. This stems from theoretical considerations (due to a small parallel wave vector component, see Treumann et al., 2004;Tsurutani, Echer, et al., 2011) and the fact that MM, as linearly polarised waves, are predominantly growing at large angles to B. Such directional change can be measured by the azimuth and elevation angles of the magnetic field, as shown in Figure 4. 2. Instability criterion MMIC (Equation 1): because large uncertainties in deriving reliable temperatures from the diagonalization of the higher-order moment (the pressure tensor) may arise due to composition assumptions, caution must be exercised when calculating MMIC in the inner magnetosheath and, especially, close to the MPB, where heavier ions may be present. In the time span where MM-like structures were found in the present study (11:26-11:30 UT), STATIC energy-mass spectrograms agreed with SWIA and additionally showed that both He 2+ and heavier ions such as O + were present in the higher energy-per-charge levels. However, choosing a reduced energy threshold to restrain the contamination did not significantly impact the determination of the moments of the ion distribution function with SWIA. Consequently, in this study, the moments were calculated in the full range 25−25 × 10 3 eV/q.
In order to calculate the MMIC criterion from Equation 1, ion pressures P ‖ and P ⊥ in the parallel and perpendicular direction with respect to the ambient B -field must be first estimated. This is done using the full pressure tensor calculated as the moment of order 2 of the ion distribution function from the appropriate SWICA mode (maximum angular resolution available). Temporal resolution of the SWICA mode for this study is 8 s. The pressure tensor is symmetric by construction and consists of nine elements, that is, 3 diagonal elements P xx , P yy and P zz , and a total of 6 off-diagonal elements P xy = P yx , P xz = P zx and P yz = P zy .
Two methods are usually adopted to obtain P ‖ and P ⊥ . They are presented in more detail in Appendix B. In short, Method 1 uses a direct diagonalization of the pressure tensor, whereas Method 2 first rotates the pressure tensor into a Mean Field-Aligned coordinate system (MFA, see Appendix A) with one direction along the mean B-field and diagonalizes the remaining 2 × 2 tensor to obtain the two perpendicular directions . Method 1 may provide a quick first estimate of directions but implicitly assumes that one of the principal axes of the tensor ellipsoid is aligned with the magnetic field direction, an hypothesis that is usually false. Moreover, off-diagonal terms still present in the MFA coordinate system have physical significance: indicative of non-gyrotropy (Che et al., 2018;Swisdak, 2016), they correspond to shear stresses in the ambient plasma, which Method 1 ignores. Thus, Method 2 is preferred in the following. The MFA direction is evaluated by taking the low-pass filtered B bg (as calculated in Section 2.2) since the average field over the SWICA measuring time of 8 s (taken as ±4 s around the middle of the scanning interval) gives almost identical results in the time span considered in this study (Appendix A).
As discussed in Halekas, Brain, et al. (2017), the limited field of view (FOV) of the SWIA instrument (360° × 90° for the SWICA mode) impacts the quality of the retrieval. At times when the ambient magnetic field direction is in the FOV of SWIA, one of the perpendicular components may be difficult to estimate (the maximum perpendicular eigenvalue of the pressure tensor taken as an upper limit). At times when it is not in the instrument's FOV, the derivation of the parallel pressure, and hence parallel temperature, will suffer from comparatively larger uncertainties. In the interval of study where MM-like structures were detected (2014-12-25 11:26-11:30 UT,), 7 of 28 the background magnetic field direction was unfortunately not in the FOV of the instrument (with angles between the Z axis of SWIA and the magnetic field direction exceeding 135°, see Figure 2). This made the determination of the parallel component of the pressure tensor difficult. As a result P ‖ was likely underestimated, as the flux integration over all angles could not capture the full 3D ion velocity distribution function (VDF) in that direction. However, magnetosheath ion VDFs are not beam-like and may spread over several angular sectors of SWIA. A cursory examination of the ion VDFs during this interval showed that the peak of the 3D distribution (in the bulk plasma direction) was within SWIA's FOV, so that the retrieved parallel component of the pressure tensor was more representative of the core contribution of the ions. As such, they represent a good proxy of P ‖ and T ‖ .
The final derivation from the pressure tensor of the parallel and perpendicular temperatures (T ‖ = P ‖ /N i and T ⊥ = P ⊥ /N i with the temperatures in eV), as well as the plasma-β ⊥ with respect to the background field B bg assumes only one ion species (protons), enabling the estimation of the MMIC. For a faster and automatic detection bypassing the MMIC criterion, the normalised (electron and ion) density variations in anticorrelation with |B|field variations proved to be an adequate and sufficient tool to confirm the presence of MM structures during the time intervals considered in this study.

General Description
On 25 December 2014 around 09:46 UT and around 11:40 UT, two prominently quasi-perpendicular bow shock crossings took place. The first of these crossings was extensively studied by Burne et al. (2021) who noted the quiet upstream solar wind conditions with no coronal mass ejection (CME), or other solar wind transients. These authors concluded that the mass-loaded shock around those times exhibited well-defined supercritical features, such as foot, ramp and overshoot, with scales reminiscent of those found at Earth, making it a baseline example at Mars. The magnetosheath behind these two quasi-perpendicular shocks is thus expected to be in a so-called 'quiet' state, with no external disturbance in the solar wind flow.
The corresponding orbit of MAVEN during those two crossings is shown in MSO coordinates in Figure 1 between 09:00 and 12:00 UT, and shown together with the statistical position of the bow shock (BS), the direction of the background magnetic field (green) and that of the normal to the shock surface (blue). The normal to the BS was calculated using the analytical model of Gruesbeck et al. (2018) assuming the shock surface to be smooth. Moving first from the solar wind to the magnetosheath, MAVEN resided ∼1 hr in the magnetosphere, circling the planet and exiting again into the upstream solar wind around 11:40:16 UT in the (+Y MSO , +Z MSO ) dayside quadrant. The angle between the average magnetic field and the normal to the bow shock was found for this crossing to be θ Bn = 76°, assuming a smooth bow shock surface. Figure 2 presents an overview of magnetic field and plasma measurements between 11:20 and 11:50 UT. Panels (a) and (b) show the total magnetic field intensity |B| and magnetic cone and clock angles as a measure of the field's direction, defined as cone = arctan and clock = arctan ( ∕ ) . A cone angle of 0° (180°) implies a sunward (antisunward) B-field direction, whereas a clock angle of 0° (90°) indicates a direction in the +Y MSO (+Z MSO ) direction. Several regions encountered by MAVEN in its exit toward the upstream solar wind are highlighted: after spending some time in the magnetosphere (labeled 'MSp') with large rotations of the field first directed roughly antisunward and then abruptly rotating toward the sunward direction, magnetic field fluctuations become attenuated indicating the presence of the MPB (depicted as a gray zone) when transiting toward the magnetosheath (labeled 'MSh') around 11:24:45 UT. In the magnetosheath, several quasi-periodic dips and peaks are observed in sequence in the magnetic field data, about ten minutes before the quasi-perpendicular BS crossing into the solar wind ('SW'). Field directions remain on average fairly constant in the magnetosheath, but tend to fluctuate more and more when the spacecraft closes in on the shock structure. In the upstream solar wind, no indication of MHs was found.
Panel (c) of Figure 2 shows for reference the ion velocity measured by SWIA in SWICA mode, which is valid mostly in the magnetosheath. The plasma is moving in the magnetosheath at velocities of about 150 km/s on average, with little change in direction and amplitude until the spacecraft exits into the upstream solar wind around 11:40 UT. Panel (d) calculates the angle, noted α B−V , between the background magnetic field direction (downsampled to SWIA's 8-s resolution) and the ion velocity vector measured by SWIA. The background magnetic SIMON WEDLUND ET AL.
10.1029/2021JA029811 8 of 28 field B bg is obtained from a low-pass Butterworth filter of the 1-s B-field data, as explained in Section 2.2. The angle remained close to about 100° with little variation around that value (±6°) during the entire interval shown here, that is, the plasma flow and the magnetic field stayed quite perpendicular to each other throughout. This could result in ion pickup ring distributions favoring the emergence of instabilities in the solar wind as well as downstream of the shock (Price, 1989).
Panels (e) and (f) display the ion and electron differential flux spectra from SWIA and SWEA (omnidirectional), respectively, with Φ j the differential particle flux of charged species j (with j = [ions,e − ]). It is important to remark that SWIA's field of view during this interval does not encompass the direction of the magnetic field, resulting in a likely underestimate of the calculated moments along the parallel direction to the field (See Section 3.2). In the magnetosphere ('MSp'), two separate ion populations of different energies but with rather constant fluxes can be seen, one below 100 eV/q, the other above 400 eV/q, likely mixing fast protons and heavier species produced by EUV photoionization of Mars' ionosphere. This coincides with lower electron fluxes with energies around 20-25 eV, typical of photoelectron energies. In the magnetosheath ('MSh'), the ion spectra broaden in energy with a main population of protons at a mean energy of about 200 eV/q and large flux fluctuations that appear related to the large magnetic field fluctuations seen then. Similar effects are seen on the electron spectra, with an intensification of electron energies and fluxes together with a high flux variability. Ion differential fluxes reach about 4 × 10 5 cm −2 s −1 sr −1 eV −1 whereas mean electron fluxes are around 10 8 cm −2 s −1 sr −1 eV −1 .
Closer to the bow shock around 11:37 UT, in the turbulent wake and transiting into the solar wind, the magnetic field varies a lot and a steady drop in ion energies with an increase in ion fluxes can be seen. This ion energy drop is concomitant with a broadening of the electron energy distribution, with electron fluxes above 10.1029/2021JA029811 9 of 28 2 × 10 8 cm −2 s −1 sr −1 eV −1 around 25 eV energy. In the solar wind, two stable main populations of ions are present, interpreted as protons (yellow-orange broad line) and He 2+ particles (blue-green thinner line) at about 600 and 1,200 eV/q, equivalent to a bulk speed of ∼340 km/s.
The upstream solar wind conditions after 11:41 UT were nominal in the sense of Slavin and Holzer (1981), with |B| ≈ 6.6 nT, V sw ≈ 340 km/s, T i ≈ 7 eV, T e ≈ 12 eV and n i ∼ 20 cm −3 , all in agreement with the previous orbit of MAVEN (Burne et al., 2021), except for slightly higher plasma temperatures and ion density. The magnetosonic Mach number ms = sw∕ √ 2 + 2 (with Alfvén and sound velocities V a and V s ) had a value around 6 (supercritical shock), whereas the solar wind ion plasma-β reached to about 1.3. This indicates that the solar wind conditions remained quite steady over the 3 hr between the two quasi-perpendicular crossings, with no CME or other solar transients effects disturbing the Martian plasma environment. The wave content for the full interval is shown in Figure 3, where the power spectral densities (PSD) for the compressional, right-handed and left-handed polarisation terms are plotted using the 32-Hz magnetic field data. First, the matrix of passage M MFA from MSO coordinates to the MFA reference frame is calculated as in Equation A1 (Appendix A), using the background magnetic field B bg direction as the z direction of the MFA system. Then the MSO B-field vector is rotated into MFA so that the two first directions are perpendicular to B bg , noted B ⊥1 and B ⊥2 , and the third one is parallel to B bg , noted B ‖ . The compressional component, left-handed and right-handed polarisations are thus defined as: where i is the imaginary unit. Cross-spectrum analysis is performed with a FFT Welch method to evaluate the PSD of each component (Welch, 1967), using a moving window of 1,024 consecutive points (32 s) for a shift of 64 points (2 s).
Most of the wave power is observed at frequencies much lower than the local cyclotron frequency, assuming only protons throughout, and, additionally, at frequencies lower than the solar wind proton cyclotron frequency (at about 0.1 Hz). In the upstream region between 11:43 and 11:50 UT, the signal appears mostly compressive with no right-handed polarisation component, which is consistent with the presence of proton cyclotron waves 10.1029/2021JA029811 11 of 28 (Romanelli et al., 2016). Downstream of the shock, the wave power is more evenly redistributed which suggests no clear relation between the upstream turbulence and the turbulence in the magnetosheath.

In-Depth Analysis
In Figures 2 and 3, two vertical dashed lines highlight one interval between 11:26 and 11:30 UT where fluctuations lodged in against the MPB appear to contain MM-like structures. This interval is shown as a zoom-in in Figure 4, which displays ion/electron energy spectra (panels a and b), magnetic field |B| (panel c), and ion densities 10.1029/2021JA029811 12 of 28 N i and electron density variations ΔN e /N e (panels e and f). The resolution of the ion spectrum and calculated moments for the SWICA mode shown here is 8 s, whereas the onboard moments have a resolution twice as good, that is, 4 s. SWICA and onboard moments agree rather well during this interval, with fluctuations closely following each other. The 2-s resolution electron densities are unavailable after 11:26:45 UT; electron density estimates calculated from the electron distribution function have afterward a reduced cadence of 8 s, which is too large to be of use for the typical width of the dips found in the magnetic field data. Hence, only the maximum resolution available is displayed; moreover, because these electron densities are not fully calibrated, only relative variations, calculated here as a forward difference, are shown.
Candidate MM detections using the B-field only criteria of Table 1 are highlighted as blue-shaded rectangles. Skewness and excess kurtosis of the interval 11:26-11:30 UT are −0.20 ± 0.16 and −1.09 ± 0.32, respectively, indicating an interval containing mostly dips and a heavy tail distribution with large fluctuations. The first large dip in B-field around 11:26:45 UT coincides with an increase in both electron and ion fluxes, flux energies and densities, an observation that seems accurate for all other detected dips in the interval, except the last one around 11:29:40 UT for which the SWICA poor spectral resolution is inconclusive, possibly indicating the slow transition from dips to peaks in the MM B-field structure (see Joy et al. (2006) for a clearer example of this effect). This is confirmed by the fact that when shifting the interval window by increments of +30 s (11:26:30-11:31:30 UT, etc.), the skewness, initially negative, steadily increases to become positive for a shift of +1 min. Throughout these periods, both magnetic field azimuth and elevation appear very stable (panel d), which is consistent with the presence of MM structures. Looking more closely into the magnetic field measurement (panel c), it appears oscillations in the total B-field take place starting around 11:27:30 UT, with a quasi-period of ∼9.3 s, as determined by a Fourier analysis of the smoothed detrended signal. Several of these peaks or dips are not captured by the automatic detection algorithm. The gain in temporal resolution (4 s) from the onboard ion datasets allows the ion density variations to follow even more closely those of the magnetic field throughout the interval. A supplementary indicator is given for the first MM detection (dip) around 11:26:45 UT by SWEA's ΔN e /N e at a resolution of 2 s, which allows this 8-9 s-long structure to be probed with almost 4 electron density data points. The B − N i and B − N e anticorrelations for this event are striking, and represent the clearest and most unambiguous MM signature detected to date at Mars. The consistency found between onboard and manually calculated moments, as well as the electron density data for this event, gives confidence for the rest of the interval. Moreover, it suggests that possibly every period with increased electron and ion fluxes, energies and densities with respect to the background levels in Figure 4 are MMs proper. The only exceptions are maybe toward the end of the interval, where the two last structures display a very rough anticorrelation between ion density and magnetic fields, with eigenvalue criteria only just fulfilled, in contrast to the previous structures. When studying in more detail the electron spectrograms at those times (panel b), electron fluxes and magnetic field intensities do appear anticorrelated, giving more confidence that these are indeed MMs too.
Out of a measurement window of 240 s, 21 s fulfilled the MM criteria above, resulting in 4-5 separate structures if two structures are separated by at least 8 s (that is, half of the MVA window and, also, the temporal resolution of the SWICA moments in this region), as in Volwerk, Schmid, et al. (2016). The B-field-only algorithm captured many prominent MM structures that can be identified by eye, but discarded those that are deemed ambiguous. Moreover, because the magnetic field background level seems to cut through the middle of the oscillations in the second part of the interval, several promising structures were missed altogether.
Uncertainties on the MVA directions and magnetic field components were evaluated following the prescriptions of Sonnerup & Scheible (1998). During the period concerned, the mean statistical error on the B-field magnitude along the maximum variance direction was less than 5 nT, that is, about 10% and below of the background field. Hodograms of the B-field variation in the MVA orthogonal planes (not shown) displayed clear linearly polarised signatures. This is giving additional confidence on the validity of the detections performed. For a discussion of SIMON WEDLUND ET AL.
10.1029/2021JA029811 13 of 28 MVA and how it is applied to the detection of ULF waves, the reader is for example referred to the review of Volwerk (2006).
Let us now investigate further the points raised in Section 2.3, that is, magnetic field direction and instability criterion MMIC. First, during the interval, the magnetic field azimuth and elevation angles do not change significantly, with less than 15° rotation (standard deviations ≤ 7.5 • ). More precisely, angular changes across each MMlike structures are less than 5° on average, in excellent agreement with the theoretical and observational expectations of MMs (Tsurutani, Lakhina, et al., 2011). Second, the MMIC makes use of the parallel and perpendicular temperatures which are calculated from the pressure tensor as explained in Section 2.3 and Appendix B, with the MFA system based on B bg as defined in Appendix A. Results for P ‖ , the two perpendicular directions P ⊥1 and P ⊥2 , as well as corresponding temperatures and plasma-β ‖ and β ⊥ are shown in Figure 5. The diagonalization method yields off-diagonal components of the pressure tensor that are small in comparison to the diagonal elements, with larger fluctuations seen when entering the candidate MM structures, especially the one around 11:26:45 UT. Clear anisotropies in pressure (and hence temperatures, although their variations are in comparison somewhat dampened by the abrupt variations seen in the densities) are seen, with the two perpendicular components larger 10.1029/2021JA029811 14 of 28 than the parallel component by a factor 2 or so. Because the mean magnetic field vector was in the blind sector of SWIA's limited FOV, the parallel components in pressure, temperature and plasma-β are all expected to be lower estimates. This in turn implies that our calculated anisotropies T ⊥ /T ‖ averaging 2.2 over the studied interval are upper estimates (see Figure B1c). However, our retrieved anisotropies are certainly not anomalous as values of 2 or more are common in the Martian magnetosheath . Figure 6 shows the interplay between magnetic and ion pressures (parallel and perpendicular) (panel a). Magnetic pressure P B = |B 2 |/2μ 0 and ion pressures, especially the perpendicular pressure P ⊥ , appear to be rather well anticorrelated throughout the interval, at the resolution of SWIA. This is in agreement with results at comets (Mazelle, 1990;Mazelle et al., 1991).
In order to give a measure of the non-gyrotropy of the pressure tensor, the quantity ζ (Equation B8, Appendix B), which we refer to as "gyrotropy index" based on the off-diagonal elements of the pressure tensor following Swisdak (2016), is shown in panel (b), left axis. Any value of ζ ≠ 0 is a departure from the normal distribution (for which ζ = 0), which is observed here throughout the interval, with ζ = 0.05-0.12 (mild departure from perfect gyrotropy). Early in the interval the gyrotropy index is much larger: this can be ascribed to the presence of heavy pickup ions, confirmed by a cursory look at the mass spectra from STATIC (see also Figure 2), essentially unmagnetized early in their trajectory. It is interesting to note that when entering the most prominent MM dip structure around 11:26:45 UT, ζ decreases together with the magnetic field intensity; consequently, gyrotropy index and perpendicular pressures appear also anticorrelated, as can be seen later in the interval. This implies that in an MM structure, the off-diagonal elements of the pressure tensor tend to diminish with respect to the diagonal (parallel and perpendicular) elements.
How does that compare with spatial scales probed by MAVEN throughout the interval? One way of evaluating this aspect is to define a quantity commensurable with a "fluid-like" gyroradius for core-distribution protons, that is, an equivalent spatial scale calculated as = ⟂∕ | | , assuming most ions (the core ions) have a perpendicular velocity equal to V ⊥ , the perpendicular bulk plasma velocity. In this way, would be a rough estimate of the mean ion gyroradius for ions populating the maximum of the non-Maxwellian ion distribution function. The decrease in gyrotropy index ζ coincides thus with (panel b, right axis) increasing by almost a factor 2 in mirror mode 'dip' structures, because of a corresponding decrease of the magnetic field intensity. One tentative interpretation in a kinetic sense could be that ions with larger gyroradii (and potentially seeping into the mirror is given as a comparison to previous works, theoretical (Hellinger et al., 2006) or observational .
SIMON WEDLUND ET AL.

10.1029/2021JA029811
15 of 28 loss cone of the 'magnetic bottle') escape, eventually tending toward a more gyrotropic distribution. This effect could become especially important for smaller MMs or MHs. We note also that and the local thermal proton gyroradius, determined from T ⊥ = P ⊥ /N i , have very similar values during the interval considered.
The MM instability criterion MMIC is shown in Figure 6c. The region crossed by the spacecraft and containing the train of MMs is interpreted as being marginally unstable to the generation of MM structures, with MMIC moving from positive values before 11:26 UT (0 ≤ MMIC < 1) to negative values afterward. Significantly negative values (MMIC ≤ −1) start to occur starting 11:28:30 UT and onwards, coinciding with the start of the quasi-periodic B-field signal. However, because of the likely overestimate of T ⊥ /T ‖ , the exact time where the MMIC may cross the threshold line MMIC = 0 is likely inaccurate. In order to investigate how changes in T ‖ may modify this interpretation, we decided to increase by increments of 10% the value of T ‖ up to a 50% overall increase, which we expect to be much larger than the uncertainty on the temperature retrievals (see discussion in Section 2.3), hence creating an "envelope of confidence" for our interpretation. This is illustrated as a light blue region above the nominal level in black in Figure 6c. This increase shifts the curve upwards toward the stability line at MMIC ∼ 0 − , making the plasma reach marginal stability almost all the time during the interval. At 11:26:45 UT, in the middle of the most prominent MM structure found (vertical light gray area on the figure), the derived MMIC remains very stable and close to the threshold. Combined with predominantly MM-unstable conditions at play for times after around 11:29 UT, this would imply that the region crossed by MAVEN at that time is not the source region of these MMs and should be found upstream of the spacecraft. Owing to their duration (around 10 s, i.e., of the order of 1,000 km in size for a lower-estimate plasma velocity of 100 km/s) and possibly advanced development stage (mixed presence of dips and peaks), this source region may be up to several hundreds of kilometres upstream in size, and not necessarily in the path of the spacecraft.
In Figure 6d, the temperature anisotropy T ⊥ /T ‖ is shown as a function of proton plasma-β ‖ and time (color code), with a Levenberg-Marquardt fit of the form: as in Fuselier et al. (1994) and Hellinger et al. (2006) to investigate the marginal stability condition with respect to MM and ion cyclotron instabilities. The modified MMIC criterion translates to [a = 0.77, b = 0.76 c = 0.016] (Hellinger et al., 2006), with any point on the right of that line denoting MM-unstable conditions. In the interval of study, a = 1.26 and b = 0.12 with c = 0, which is in contrast with previous studies at Earth of the proton anisotropy: Anderson et al. (1994) obtained values a = 0.85 and b = 0.48 whereas Génot, Budnik, Hellinger, et al. (2009) found a = 0.47 and b = 0.56 (c = 0 in both cases). A flatter inverse correlation with β ⊥ (a = 1.41, b = 0.09, not shown) was found in the study interval. This in turn is closer to the values found by Fuselier et al. (1994) (a = 1.40, b = 0.26) for He 2+ ions, also in the Earth's magnetosheath. Such differences are not surprising since (a) our present analysis does not discriminate between protons and other ions (including alpha particles) in the ion distributions, and (b) parallel temperatures are likely underestimated because of the FOV of SWIA. Moreover, we limited ourselves here to one event spanning a few minutes only. A statistical study of the temperature anisotropy behavior with respect to the plasma-β, using the full MAVEN data set, is left for another study.

Associated Scales and Validation of B-Field Criteria
A reanalysis of the current interval is now possible on the basis of the plasma measurements and the high confidence level in the density variations. The temporal width of each structure is evaluated by remarking that the first crossing between the background magnetic field level on either side of the detection constitutes the start and end points of the structure. It is clear that each event has a different scale in the parallel and perpendicular directions, as remarked by Zhang et al. (2008): this is seen in Figure 7a for the clearest event around 11:26:45 UT, when comparing the parallel and cross-field magnetic fluctuation components normalised to the background field, δ|B ‖ | = |B ‖ − B bg |/B bg and δ|B ⊥ |.
For this particular event, one can specifically check for the species-by-species phase relationship between density and the magnetic perturbations, which we apply here to the ion measurements. It is approximately, in the hypothesis of a cold component of the plasma, in this case the cold photoelectrons (Hasegawa, 1969;Mazelle et al., 1991): Equation 6 formalizes that for a temperature anisotropy in the perpendicular direction (T ⊥ /T ‖ > 1), the density is expected to be out of phase with the magnetic field. In Figure 7b, we show the comparison between the left-hand side and the right-hand side of the equation, with δN i = ΔN i /N i . In order to evaluate the uncertainty in the parallel temperature (derived from SWIA's limited field of view), the shaded green region denotes a 50% variation of T ‖ (1-1.5 × T ‖ ). Although the δB term has larger variations in the time span considered, the agreement, both in shapes and amplitudes of variations, is striking: it quantitatively confirms the anticorrelation between density and magnetic field fluctuations and the ion density. Similarly, Mazelle et al. (1991) found that Equation 6 was also accurate for the low-frequency MM-like signatures found at comet 1P/Halley.
To try to further characterize the particle content of this MM, the mean electron and ion omnidirectional fluxes measured by SWEA and SWIA are shown in Figure 7c where mean fluxes were derived successively inside and outside the structure (see details in the caption). As already remarked from Figure 4, the inside of the MM (with δB = (|B| − B bg )/B bg ) were calculated using a low-pass Butterworth filter to estimate the background contributions N bg and B bg . The light green-shaded region illustrates an increase in T ‖ of up to 50% to evaluate how sensitive the results are to an underestimation of T ‖ (see text for explanation). (c) Particle differential electron and ion fluxes inside and outside the MM structure. Black lines correspond to the mean fluxes inside of the structure (gray zone on panels a and b) whereas the blue and red lines correspond to mean fluxes 8 s (= resolution of SWICA) outside of the structure (light blue and red zones on panels a and b). Notice that the zones do not start exactly at a full second since the original temporal width of the MM was determined by an interpolation of the 1 -Hz total B-field intersecting the background field level B bg . 17 of 28 structure corresponds to a distinct increase in both peak energy and flux with respect to the outside, both for electrons and for ions. Interestingly, the fluxes on either side of the MM are similar, suggesting that the surrounding plasma may well have very similar characteristics. This paints the picture of a magnetic bottle containing a relatively energetic and dense ion population, and drifting with the surrounding cooler and diluted plasma.
Reanalysis of all the MM structures present in the 11:26-11:30 UT interval proceeded as follows. First, we manually select the location and temporal width of 9 "dip" MMs for simplicity, assuming that their peak counterparts are part of the same structure, by comparing the |B|-field variations with the onboard moment and SWICA density variations. Their average width is 8.7 ± 2.1 s, with an average dip-to-dip period of ∼9.3 s. Four structures last up to 10 s whereas the shortest around 11:28:28 UT last only about 5 s. Assuming a plasma bulk speed of about 150 km/s (mostly in the perpendicular direction to the magnetic field as shown by SWIA), and that the spacecraft is at rest with respect to the plasma, the size of the MM structures varies between 750 and 1,500 km (a large fraction of the planetary radius, 0.22-0.44 R p ), which is much larger than the equivalent spatial scale of the peak-distribution ions in this region ( ∼ 35 km , i.e., MM sizes of the order of ∼20-40̃ , see Figure 6b). By comparison, we can estimate the thermal proton gyroradius, noted = √ 2 ⟂∕ ∕ | | , in the possible source regions, either in the solar wind or in the immediate wake of the quasi-perpendicular shock. In the solar wind upstream of the shock, around 11:45 UT, the perpendicular proton temperature is 6.5 × 10 4 K (i.e., T ⊥ = P ⊥ / N i ∼ 5.6 eV as measured on average from the SWIFA mode, a textbook value, see Slavin & Holzer, 1981) and the upstream B-field magnitude was ∼6.6 nT. This corresponds to a solar wind proton thermal gyroradius sw of the order of 50 km (∼0.015 R p ), implying that MM sizes vary from ∼15 to 30 sw throughout this interval. Such an estimate, albeit significantly smaller, is in line with similar considerations at Venus made by Zhang et al. (2008), who deduced MHs to have sizes of 40-100 r p in the shape of a prolate spheroid. Similarly, just behind the shock around 11:37 UT where the magnetic field intensity is lower than at the location of the MMs, T ⊥ ∼ 150 eV, and |B| ∼ 24 nT, resulting in a magnetosheath thermal gyroradius ms ∼ 75 km, that is, the observed MMs deeper in  Table  2 of Tsurutani, Lakhina, et al. (2011)).
We also look into the magnetic field power spectral density (PSD) estimates of the compressional, right-handed and left-handed components during the full interval, computed with Welch's algorithm and a 50% Hamming window overlap (Welch, 1967). This is shown in Figure 8, with the mean MM period of 9.3 s (0.017 Hz) appearing to almost coincide with a series of larger peaks in the compressional PSD, with the two circular polarizations being muted in comparison. During the interval considered, the mean proton gyrofrequency ̄ was equal to 0.64 Hz: above this limit, sub-ion scales start. Two main trends in the overall PSDs can be seen, with a Kolmogorov-type turbulence (spectral index − 1.61 ∼ −5/3 in the fluid limit) steepening above the proton gyrofrequency (spectral index − 2.66 = −8/3). The mean period of the MMs lies thus in the Kolmogorov-like region of the spectrum, as expected from studies of the low-frequency wave turbulence in the Earth's magnetosheath (Sahraoui et al., 2003(Sahraoui et al., , 2004 or from simulations (Hellinger et al., 2017). This is also in agreement with the study of Ruhunusiri et al. (2017) who calculated turbulence spectra in the environment of Mars with MAVEN magnetometer data. They found on average spectral indices of the order of −1.43 and −2.50 close to the MPB, values which are close to our estimates. This has one main consequence for the evolution of our detected MMs: the environment in which these MMs are embedded has had time to evolve into a fully developed energy cascade. (g) angles between background field and maximum (respectively, minimum) variance Θ maxV (Φ minV ) directions. Initially retained ranges from Table 1 are highlighted in transparent blue, whereas the revised criteria, determined from the PDFs within 1σ of the mean, yield the ranges in orange. Normal distribution fits are shown as blue lines. Out of 240 points in the interval, the total number of revised measurement points was 33 (initial criteria: 21 points), thus equal to a total residence time in a MM structure of 33 s.
SIMON WEDLUND ET AL.

10.1029/2021JA029811
19 of 28 Finally, manually picking the location and duration of MMs in our interval results in the statistics shown in Figure 9, with 9 structures totaling 78 s, from which a revised set of B-field only criteria was extracted: this new set of criteria has the advantage of being unbiased with respect to other choices of criteria adopted at Venus and Earth. Criteria were on average relaxed (Δ|B|/|B bg |, λ max and λ min , Θ maxV ) to better take into account the characteristics of the detected MMs. To detect as many MMs as possible and limit false positive detections, we decided to keep within 1σ of the mean of the PDFs shown in Figure 9, depending on the spread in values within the range. The revised detections are highlighted in orange, whereas the detections using the initial criteria are in light blue; common periods are in gray. The revised set of criteria is given in Table 1, yielding 7 'dip' structures and 3 'peak' structures, with 1 potential outlier, totaling 33 s out of a 240 -s interval. Some of these peak structures appear to be nothing more than counterparts of dip structures, depending on the choice of the background field. Their location with respect to detrended B-field values, where MM structures have typical peak-to-peak amplitudes of 20-35 nT, are also emphasized as areas in light orange on the figure.
Robust detections across an entire orbit may be difficult as small adjustments in angles and eigenvalue ratios can change the results, increasing or decreasing the number of potential candidates. One difficulty lies in the aforementioned mixed presence of peaks and dips in the study interval considered. Automatic detection applied to larger datasets may in this way benefit from considering smaller intervals of a few minutes on which to calculate the background field levels, depending on the skewness of the B-field distribution in that interval (Ala-Lahti et al., 2018;Génot, Budnik, Hellinger, et al., 2009). However, we stress that any attempt at detecting MMs based on B-field only measurements will always result in false positive detections and that plasma data is necessary for a positive identification.

Discussion
Although magnetic field signatures in the magnetosheath of Mars were previously found consistent with the probable presence of MM waves (Bertucci et al., 2004;Espley et al., 2004), no dedicated study has been possible until the advent of MAVEN and its full magnetic field and high-time resolution plasma suite of instruments. Consequently, we presented in this study the first in-depth characterization of MMs in the Martian environment.
Where do these structures originate? The question of the origin and development of MM structures around an unmagnetized object such as Mars, Venus or around comets is of particular interest to the space physics community. Through competition with the left-hand Alfvén ion cyclotron mode (which heats the plasma), MMs consume the locally available free energy of the temperature or pressure anisotropy and may help in the thermalization of the magnetosheath (Shoji et al., 2009;Soucek et al., 2008). The source of the anisotropy can be twofold: (a) it may come from the intrinsic properties of the plasma in the wake of the quasi-perpendicular shock which heats and deflects the solar wind (Bale et al., 2005;Peng et al., 2015) or (b) from the pickup ion process several proton gyroradii in the upstream solar wind (Bader et al., 2019;Price, 1989). The first mechanism occurs at magnetized (Earth, Jupiter) and unmagnetized/weakly magnetized planets (Mars, Venus) alike, and implies a relatively local MM generation in the magnetosheath. The second mechanism, the pickup ion process, is specific to unmagnetized objects with an extended exosphere such as Mars, Venus or comets and creates ring/ring-beam ion velocity distributions which in a high-β plasma are unstable to the generation of MMs (Price, 1989): it can drive their growth everywhere where the exosphere extends, from the upstream solar wind and in the magnetosheath down to the MPB. In the solar wind, MM-like structures, called MHs or magnetic depressions, have been detected upstream of Mars (Madanian et al., 2020) and at comets (Volwerk et al., 2014). This mechanism may become prominent when the bow shock is either quasi-parallel, that is, less conducing to the rise of anisotropies and more permeable to wave transmission, or weak (for example at comets, and in a lesser measure, at Mars); it may then constitute an important remote source of MMs. Indeed, if the source of the anisotropy is upstream of the shock, MM waves found downstream of the shock but originating upstream of it may consequently have been transported over large distances, almost unchanged, through the quasi-parallel shock. In this latter scenario, if we recall that MH structures in the solar wind are somewhat larger than the MMs observed in the sheath, a damping or attenuation process could be at play when crossing the shock into the magnetosheath, with the trapped particles within the magnetic bottles losing energy. This scenario is reminiscent of the one found by Plaschke et al. (2018) at comet 67P/Churuymov-Gerasimenko for low outgassing rates (and at about the time a nascent bow shock was expected to take take place), with relatively smaller MH structures deep in the coma found to be the end result of upstream larger MHs.

of 28
In our case however, and although contributions from local pickup ions cannot be ruled out entirely, the first mechanism (source in the wake of the quasi-perpendicular shock) seems the likeliest. 1. First, as we approach closer and closer to the quasi-perpendicular bow shock, the plasma becomes more and more unstable to the generation of MMs, with the detected structures close to the MPB mostly embedded in a plasma marginally stable/unstable to the generation of MMs. 2. Second, the average angle α B−V between the magnetic field direction and the plasma ion velocity is extremely stable around 103° ± 5° (see in Figure 2). This implies in velocity phase space that the ion velocity distribution functions (VDF) of newborn pickup ions are closer to a ring distribution, which preferentially contributes to a large increase in P ⊥ of these newly picked-up ions as compared to the parallel direction (Price, 1989), in agreement with our derived anisotropies. A cursory examination of the VDFs measured by SWIA during the interval (not shown) reveals indeed a non-Maxwellian (or rather non-bi-Maxwellian) behavior, with free energy available for light particles, such as protons, to drive the growth of instabilities. This may in turn favor the growth of the MM instability over that of the proton cyclotron instability owing to large plasma-β and the presence of a small quantity of heavier ions in the plasma (Price et al., 1986), in keeping with the observations of Russell et al. (1987) at comet 1P/Halley. Moreover, because the α B−V angle is large also in the solar wind (≈ 98 ± 5 • ) as displayed in Figure 2, an extra source of anisotropy upstream of the shock may be present, driving MM unstable conditions in those remote locations, with potential mode conversions through the boundary into the magnetosheath. 3. Third, the train of MMs found here shows a variety in shape (dips and peaks), duration and size that point toward different stages in their evolution. 4. Fourth, the turbulence spectrum during our events and shown in Figure 8 displays a typical fully developed energy cascade with Kolmogorov-like behaviors, implying that the plasma turbulence has had time to evolve from its generation source. 5. Finally, no magnetic holes are found immediately upstream of the shock, at least in the path of the spacecraft.
All five points above are consistent with the scenario of a nonlocal source region for these MMs likely situated in the bow shock's immediate wake where a large part of the anisotropy grows.

Conclusions
For the first time, we established the unmistakable presence of a train of linear MMs in the magnetosheath of Mars using in combination MAVEN's magnetic field and high-resolution plasma measurements-ions and electrons-in the early part of the mission. Their characteristics have a textbook flair to them: with wave frequencies below the local proton gyrofrequency and lasting approximately 5-10 s (as previously found at Venus by Volwerk, Richter, et al. 2016), these magnetic islands follow a clear B − N i and B − N e anticorrelation. Trapping distinctly denser and more energetic ions on the inside, they appear to float in a relatively less energetic and more diluted ambient plasma.
These MMs are situated on the dayside in the wake of a highly quasi-perpendicular supercritical shock (as studied by Burne et al., 2021), deep in the magnetosheath (R sc ∼ 1.35 R p , solar zenith angle SZA ∼ 60°) and close to the transition into the MPB: this is precisely this region where Ruhunusiri et al. (2015) statistically found waves matching MM characteristics. The upstream solar wind conditions are quiet (|B| ≈ 6.6 nT, V sw ≈ 340 km/s, T i ∼ 7 eV, T e ∼ 12 eV, n i ∼ 20 cm −3 ), with no CME or other transient solar effects detected at that time. Because of the nature of the shock, magnetic holes created upstream in the solar wind may not easily cross the quasi-perpendicular barrier, which could imply a more local generation mechanism behind the bow shock. Moreover, the magnetic field and plasma signatures of these MMs seem to be in constant evolution: they vary in size (from a few hundreds of kilometres to a few thousands of kilometres) and in peak-to-peak amplitude (from ≈ 20 to 35 nT with respect to the background field), possibly implying that these structures are in various stages of evolution. There is also indication that the MMs transit from peaks to dips in a region where the plasma is marginally unstable to MM generation. This could in turn imply that the generator region of these particular MMs, if assumed to be immediately downstream of the shock region, was never probed by the MAVEN spacecraft. Although the quasi-perpendicular shock and its turbulent wake may be the main source of the free energy in the plasma, our results do not exclude two possibilities for the generation of these MMs: (a) local source of anisotropy through pickup ion processes upstream and downstream of the shock, and (b) MHs forming in the upstream solar wind and be-21 of 28 ing transported through the quasi-perpendicular bow shock into the magnetosheath. However, such a crossing would in practice be difficult and raise questions as to the energy/momentum transfer and wave mode conversion through the shock. Finally, the question of their transport from the remote source region behind the shock down to the magnetic pileup boundary where they are detected appears reminiscent of the apparent accumulation of MMs at the Earth's magnetopause (Erkaev et al., 2001;Omidi et al., 1994;Shoji et al., 2009) and at the Venus induced magnetospheric boundary (Volwerk, Schmid, et al., 2016).
Further tests of this idea could be done with the use of dedicated plasma flow models in the magnetosheath, as done in Guicking et al. (2010) for Venus and Soucek et al. (2015) for Earth, to reconstruct the possible history of MMs along streamlines up to the shock region. Using the magnetic field-only detection criteria developed in the present work, statistical studies of the location of MMs in Mars' magnetosheath (currently under way) should shed some more light on the place(s) of origin of mirror mode structures, under which conditions they predominantly appear and if the accumulation toward the magnetic pileup boundary is systematic or not.

Appendix A: Magnetic Field-Aligned Coordinate System
The Magnetic Field-Aligned system (MFA) can be determined by first calculating the ambient magnetic field direction B ambient . In the case of an ion instrument such as MAVEN/SWIA, which takes a finite time in order to scan through angles and energies, one way of estimating B ambient is to calculate the average magnetic field <B> during a full measurement scan of SWIA. In the SWICA mode during the time interval chosen in this study (11:25-11:30 UT on 2014-12-25), one measurement takes for instance about Δt = 8 s. Another way of estimating B ambient is to calculate the so-called background field B bg from a low-pass filtering of B as in Section 2.2.
Assuming the latter for now, the matrix of passage M MFA between MSO coordinates (base ) is usually chosen so that the third axis in the MFA system is aligned with the background B-field direction, that is, ̂′ = bg∕| bg| . Because in three dimensions there is an infinite number of perpendicular vectors to a given vector, the perpendicular directions are arbitrary: a choice with respect to one perpendicular axis must be made to complete the right-hand rule. One possibility used in the literature is to assume that ̂′ is perpendicular to the position vector R of the spacecraft (Laakso et al., 2010) so that ̂′ = ×̂′ . A common choice outside of Earth studies chooses a plane containing the ̂′ direction and performs a counterclockwise rotation of 90° of ̂′ so that the rotation matrix between MSO and MFA frames becomes: This is the convention adopted in the present study. By construction, M MFA is orthogonal, so that the transpose of the matrix is its inverse, that is, ⊺ MFA = −1 MFA , and its determinant obeys det ( MFA) = 1 . The instantaneous magnetic field vector in MFA coordinates is simply: for column vectors. In the convention above, B ‖ is thus along the z direction and contains most of the signal, whereas the perpendicular directions oscillate closely around zero.
Tests (not shown) were performed in the interval 11:26-11:30 UT on 2014-12-25 to determine the rotation matrix M MFA using as ambient field either the ±4 s averaged B-field over the SWICA measurement time <B> (which can be argued to be more physical with respect to the ion measurements), or the low-pass filtered field B bg (used in the B-field only approach, a measure of the macroscopic ambient field over a 2-min span were also similar, with only occasional spikes when using ⟨ ⟩ linked to the temporal resolution of SWIA (see Figure B1). In the main study, the low-pass filtered background field, noted bg , was used for simplicity and to keep consistent with the MM detection algorithm.

Appendix B: Deriving Parallel and Perpendicular Pressures
As explained in Section 2.3, the anisotropy in the plasma can be estimated by studying the behaviour of the ion pressure tensor in the Mean Field-Aligned (MFA) coordinate system. The pressure tensor ̄ is a symmetric second-rank tensor with 9 elements of the form:̄= in any {̂̂̂} normalised system of coordinates. It is good to recall also for later check-ups that, under any change of coordinate system, three main invariants  exist for a tensor (Cayley-Hamilton's theorem). For a symmetric matrix, they reduce to: First the pressure tensor, originally expressed in the SWIA instrument coordinate system in eV cm −3 (̄swia ), is rotated into the MSO coordinate system using a rotation matrix from the appropriate SPICE kernel to yield ̄M SO . Ideally, in a coordinate system aligned with the mean magnetic field, off-diagonal terms will tend to be small with respect to diagonal terms. Notwithstanding the system of reference, and because the trace of a tensor is invariant with respect to coordinate transforms (rotations), the scalar pressure p = P xx + P yy + P zz remains constant. It follows from this simple statement that when diagonalizing the 3 × 3 pressure tensor, the eigenvalues found will stay the same no matter which coordinate system is used.
To obtain the pressure components parallel and perpendicular to the background or average magnetic field, two methods are usually adopted, the first one based on the direct diagonalization of the 3 × 3 tensor ̄ , the other relatively more involved with a rotation first into the MFA coordinate system and a diagonalization of the remaining 2 × 2 tensor matrix.

B1. Method 1: Direct 3 × 3 Matrix Diagonalization
Becausd that it can be performed without a supporting magnetometer, this method is usually the one used for onboard moment calculations. In any coordinate system in which the initial pressure tensor is expressed, the pressure tensor can always be diagonalized so that: SIMON WEDLUND ET AL.
10.1029/2021JA029811 23 of 28 with S a non-singular matrix containing the right eigenvectors of the system, giving the directions of the principal axes of the diagonalized tensor in the chosen coordinate system. Here the parallel direction to the ambient magnetic field is assumed to be along the z direction, whereas the two other components, P ⊥1 and P ⊥2 , are assumed to be identical in the gyrotropic assumption, by definition P ⊥1 ∼ P ⊥2 ∼ P ⊥ .
In order to forego any assumption on the direction of the ambient magnetic field, one classic way to sort the eigenvalues found in Equation B5 is to isolate one that is most different from the two others. This component is by inference the parallel scalar pressure P ‖ , whereas the two remaining components are assumed to be the perpendicular pressures. The total perpendicular pressure is then given by ⟂ = 1 2 ( ⟂1 + ⟂2) . The situation with negligible off-diagonal terms in the pressure tensor is mostly encountered in the textbook cases of a spherically symmetric distribution function or a gyrotropic one, when the velocity distribution function is cylindrically symmetric about the mean field direction, here assumed along z (Paschmann et al., 1998). However, this situation seldom occurs in practice and determining the parallel and perpendicular pressures from this technique may yield false estimates, sometimes large.

B2. Method 2: Alignment to MFA System and 2 × 2 Matrix Diagonalization
When the magnetic field direction is known, this method is arguably more physical than Method 1, since the pressure tensor is analyzed directly in the MFA coordinate system itself instead of an arbitrary coordinate system (Swisdak, 2016). It thus requires the simultaneous knowledge of the magnetic field direction and amplitude. Calculations can be readily made in instrument coordinates (in which case the B-field is transformed into SWIA instrument coordinates) or in MSO coordinates (in which case the pressure tensor is first rotated into the MSO coordinate system). The latter is chosen to keep with the conventions of the main text.
Following the calculation of the matrix of passage from MSO to MFA coordinates (see Appendix A), the MSO pressure tensor can be rotated into MFA coordinates so that: is also a symmetric tensor.
By definition, ′ = ‖ , whereas ′ , ′ and ′ are non-zero off-diagonal elements indicative of shear stresses, that is, the flow of momentum in the x (respectively, y) direction by motion of the plasma in the y (z) direction. Keeping the third column of the tensor aside, and since the determination of perpendicular directions is always arbitrary (see Appendix A), the two perpendicular directions are simply obtained by diagonalizing the remaining 2 × 2 matrix (which is always diagonalizable, by construction): where S′ is a singular matrix containing the right eigenvectors of the system. By convention, the perpendicular component eigenvalues are sorted by increasing value.

B3. Comparison of Methods
A comparison of the two methods is presented in Figure B1, with a (direct diagonalization) shown in panels (b-d) as a blue line, whereas Method 2's results are given in orange (and, for comparison, yellow, when using SWIA's mode scanning resolution of 8 s). Although local differences are clearly seen, the two methods agree rather well on average. At the time when a clear MM structure is detected (gray zone), P ⊥ ∼ 2.5P ‖ for both methods. However, at the beginning of the interval, Method 1 exhibits abrupt variations of a factor 3-4 over a time scale comparable with the SWICA scanning temporal resolution (∼8 s): this is highly suspicious and linked to the rather arbitrary sorting of the eigenvalues, with the two closest eigenvalues assumed to be the two perpendicular components, in keeping with the gyrotropic assumption. Method 2, on the contrary, provides a more gradual and smoother evolution of the pressures, both parallel and perpendicular, as well as the temperature ratios. This appears more consistent with the measurements themselves and the physics at play. Except from a few less gradual point-to-point variations, a similar conclusion can be made when using method 2 and estimating the ambient field direction from the mean field <B> over the measurement time of SWIA (yellow) instead of using the macroscopic low-pass filtered B bg .
Because of the arbitrary way the parallel component is chosen in the diagonalized pressure tensor in Method 1, only Method 2 should be safely used in this particular case especially at the beginning of the interval, which has Figure B1. Comparison of retrievals for the pressure tensor parallel and perpendicular components. (a) Magnetic field |B| from MAVEN/MAG at 1 s resolution, low-pass-Butterworth filtered field |B bg | as in Section 2.2, and average field <|B|> over SWICA's mode resolution (8 s, centered), during the interval 11:26-11:30 UT on 2014-12-25. (b) Parallel P ‖ (continuous lines) and total perpendicular ion pressures (dashed line) ⟂ = 1 2 ( ⟂1 + ⟂2) , as retrieved from MAVEN/SWIA (SWICA mode). (c) Temperature (or pressure) anisotropy. (d) Mirror mode instability criterion (MMIC), with the zero line separating MM-unstable and MM-stable conditions. Throughout panels 2-4, method 1 (direct diagonalization) is in blue whereas method 2 is in orange and yellow (tensor expressed in MFA, using either B bg or <B>, respectively, to estimate the ambient B-field direction). The gray-shaded zone represents the clearest MM structure in this interval as detected by a combination of plasma and magnetic field measurements.
SIMON WEDLUND ET AL.
10.1029/2021JA029811 25 of 28 repercussions on any other quantities derived from the ion pressure tensor: T ‖ and T ⊥ , plasma-β ‖ and β ⊥ . Ultimately, this may alter the interpretation of the MMIC of Equation 1. This is shown in Figure B1d, where Method 1 would imply sharp oscillations around the MM-stable line at the beginning of the interval, which is misleading. In contrast, Method 2 predicts marginally MM-stable conditions (MMIC≈ 0 ) in this time span. After 11:28 UT, both methods agree rather well.
It is important to note here that the validity of any of those methods depends on the quality of the ion velocity distribution measured in the first place, and the field of view (FOV) of the instrument . If the direction of the magnetic field for example lies in the blind sectors of the plasma instrument due its limited FOV, the parallel estimate of the pressure tensor will become difficult to assess and likely underestimated. Likewise, if the instrument does scan through the direction where the magnetic field points, the parallel direction and only one perpendicular direction (corresponding to the largest eigenvalue) will be well-defined. Assumptions on either P ‖ or P ⊥ must thus be made depending on the case. For SWIA, the coarse mode usually adopted in the magnetosheath has an angular FOV of 360° × 90° (SWICA mode), which, despite its broad coverage in comparison to other more FOV-limited modes, needs to be checked against the magnetic field direction for each event, separately.
In conclusion, Method 2, being more physical, should always be preferred over Method 1 when magnetic field measurements are available, unless the velocity distribution function is shown to be spherically or cylindrically symmetric. Vigilant care in the interpretation is strongly recommended in all cases.