# Seasonal and Mach Number Variation of the Martian Bow Shock Structure

Jacob Fruchtman, Jasper Halekas, Jacob Gruesbeck, David Mitchell,

Christian Mazelle

► To cite this version:

Jacob Fruchtman, Jasper Halekas, Jacob Gruesbeck, David Mitchell, Christian Mazelle. Seasonal and Mach Number Variation of the Martian Bow Shock Structure. Journal of Geophysical Research Space Physics, 2023, 128 (8), pp.30. 10.1029/2023JA031759. [hal-04939514](https://hal.science/hal-04939514v1)

HAL Id: hal-04939514

https://hal.science/hal-04939514v1

Submitted on 11 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
JGR Space Physics

RESEARCH ARTICLE

10.1029/2023JA031759

**Key Points:**

* We found that the bow shock magnetic field structure has at most a very weak seasonal variation in the magnetic jump and the overshoot

* The Martian bow shock overshoot varies much more significantly with the solar wind flow Mach number than with any local shock parameters

* The Martian shock responds to solar wind variations in the same way as the Earth's shock

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

J. Fruchtman,
jfruchtman@uiowa.edu

**Citation:**

Fruchtman, J., Halekas, J., Gruesbeck, J., Mitchell, D., & Mazelle, C. (2023). Seasonal and Mach number variation of the Martian bow shock structure.
*Journal of Geophysical Research: Space Physics, 128*, e2023JA031759. https://doi.org/10.1029/2023JA031759

Received 8 JUN 2023
Accepted 8 AUG 2023
Corrected 2 SEP 2023

This article was corrected on 2 SEP 2023.
See the end of the full text for details.

©2023. The Authors.

This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

Seasonal and Mach Number Variation of the Martian Bow
Shock Structure

Jacob Fruchtman¹ ID, Jasper Halekas¹ ID, Jacob Gruesbeck² ID, David Mitchell³, and Christian Mazelle⁴ ID

¹Department of Physics and Astronomy, University of Iowa, Iowa City, IA, USA, ²Solar System Exploration Division, NASA Goddard Space Flight Center, Greenbelt, MD, USA, ³Space Sciences Laboratory, University of California, Berkeley, Berkeley, CA, USA, ⁴Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, UPS, CNES, Toulouse, France

**Abstract** Like its terrestrial counterpart, the Martian bow shock can exhibit foot, ramp, and overshoot substructures. The shock extent is known to vary with solar wind and season, but it remains an open question whether other features of the bow shock also vary seasonally. By comparing parameters of the shock during different seasons and solar wind conditions, we investigated whether the presence of seasonal planetary ions has any effect on the shock structure. We utilized a series of algorithms to automatically collect and analyze 7056 shock crossings from MAVEN spacecraft data between November 2014 and 2019. We compared the magnetic shock jump and overshoot amplitude as functions of upstream parameters and season. We found that the magnetic shock jump agrees to first order with Rankine-Hugoniot predictions, and that the overshoot amplitude dependence on critical ratio, beta and Alfvén Mach number agree with the results of previous studies. These trends are identical for all seasons. We also found that the shock jump and fast Mach number correlate primarily with solar zenith angle, whereas overshoot amplitude correlates primarily with shock extent. We found that the most likely cause of the latter was an unexpected strong correlation between overshoot amplitude and the solar wind flow fast Mach number. The seasonal variability of the Sun-Mars system does not appear to be strongly reflected in the resulting shock structure. The bow shock structure is controlled mainly by the properties of the solar wind, rather than by Mars' location in its orbit around the Sun.

**Plain Language Summary** The size of the Martian bow shock is known to vary with solar wind and season, but it remains an open question whether other features of the bow shock also vary seasonally. By comparing the shock during different seasons and solar wind conditions, we investigated whether the presence of seasonal planetary ions has any effect on the shock structure. We wrote an algorithm to identify whenever the MAVEN spacecraft crossed the bow shock, and measure relevant quantities in the vicinity to characterize the shock's behavior. We found that the shocks' measured magnetic structures' dependencies on solar wind parameters roughly agree with theoretical predictions and results of previous studies. We found that the primary and secondary substructures have negligible seasonal variation, while the latter is most strongly correlated with an entirely different solar wind parameter than we predicted—the solar wind flow directed magnetosonic Mach number—which is independent of the direction of the shock surface. Thus, the bow shock structure is controlled mainly by the properties of the solar wind, rather than by Mars' location in its orbit around the Sun.

**1. Introduction**

Mars does not have its own intrinsic global magnetic field (Acuña et al., 1992). Much like Venus (T. Zhang et al., 2008), its (non-crustal) magnetic field is induced by atmospheric interaction with extrinsic solar phenomena: First, solar EUV ionizes atoms in the dayside upper exosphere. Henceforth, mass loading of the solar wind by the Martian corona slows down the heliospheric plasma close to the planet in order to conserve momentum. The frozen-in condition of MHD plasmas therefore requires the magnetic field lines to drape around the planet (Bertucci et al., 2011; Szegő et al., 2000). Temporal variations in the solar wind magnetic field generate the electric field which can drive a current in the ionosphere that in turn produces an opposing magnetic field. Both processes lead to magnetic field pileup in front of the ionosphere (Bertucci et al., 2005). The resulting magnetic pressure slows and deflects the solar wind plasma and forces it to flow around the pile-up region. Together, these processes generate an obstacle to the flow from which fast magnetosonic waves are counter-propagated to form a “bow shock.” This bow shock separates the solar wind plasma into two regions—the shocked downstream
"magnetosheath" and the unshocked upstream solar wind (henceforth we will refer to the shocked and unshocked solar wind regions respectively by "downstream"/"magnetosheath" and "upstream"/"solar wind" and denote variables $Q$ within each region with subscripts $Q_2/Q_D$ and $Q_1/Q_U$).

## 1.1. Bow Shock Structure

Collisionless bow shocks are discontinuities in the structure of a flowing magnetized plasma (such as the solar wind) approaching an obstacle (such as Mars), where dissipative processes convert sufficient amounts of upstream kinetic energy into downstream heat in order to slow down the downstream plasma below its wave speed (in this case, the solar wind's fast magnetosonic wave speed) (Kennel et al., 1985). If the bow shock is static, locally planar (i.e., has a single local surface normal $\hat{n}$) and structurally varies only along this $\hat{n}$ direction, then its structure can be described by the MHD Rankine-Hugoniot equations (C. Russell et al., 1982). It has been noted (Kennel et al., 1985; Marshall, 1955) that when the upstream fast magnetosonic Mach number $M_{fms}$ is above some [first] critical Mach number $M_{crit}$, resistivity cannot account for all the dissipation demanded by Rankine-Hugoniot equations. Further, it has long been observed that a controlling parameter of the structure and dynamics of a [region of a] shock is the shock normal angle $\theta_{BN} = \arccos(\frac{\hat{n}\cdot B_U}{|B_U|})$: Quasiperpendicular ($\theta_{BN} > 45^{\circ}$) shocks can be highly dissimilar from quasi-parallel shocks (Gosling et al., 1989). As with much of the extant literature, each of these two will be introduced separately.

For a quasi-perpendicular subcritical shock, the magnetic profile can be well modeled to lowest order as a step or sigmoid function where the magnitude of the field jumps when crossing from upstream of the shock to downstream (Kennel et al., 1985). As such, we can quantify the magnetic profile in terms of the magnetic "jump condition" $B_D/B_U$. We will define this ratio as

$$J = B_D/B_U \qquad (1)$$

for compactness. This does not mean magnetic shock profiles are perfectly laminar, however. It has been long understood that the upstream parameters which primarily control the Rankine-Hugoniot jump conditions ($M_{fms}$, Alfven Mach number $M_A$, $\theta_{BN}$, and plasma $\beta$) also control other MHD scale shock substructures. The resistive dissipation mechanism for collisionless shocks can involve the emission of magnetosonic whistler waves (Kennel et al., 1985; Wilson III et al., 2017), which in turn appear as turbulence in the magnetic structure. Wilson III et al. (2017) showed that these whistler waves can be common and large in amplitude around quasi-perpendicular shocks even at low Mach numbers and plasma betas.

In supercritical quasi-perpendicular shocks, additional structures form in order to allow conservation laws to be satisfied. These additional substructures are within a convected ion gyroradius of the discontinuity, where the now significantly more prominent kinetic scale physics comes into play (Bale et al., 2005). Ions which do not travel fast enough in the shock normal direction when passing the discontinuity can be reflected back upstream by the Lorentz force over the course of a single gyration (Gedalin, 1996). This additional population of reflecting ions causes the region immediately downstream of the shock to become more densely populated with ions, in turn inducing a stronger magnetic field. This immediately downstream magnetic substructure is called the "overshoot" as it takes the form of a peak where the magnetic "ramp" across the discontinuity continues to increase beyond the MHD scale asymptotic downstream field before eventually dropping back down to it (C. Russell et al., 1982). This overshoot can be quantified in terms of the "overshoot amplitude"

$$A = (B_{\max} - B_D)/B_D \qquad (2)$$

where $B_{\max}$ is the maximum measured strength of the magnetic field in the overshoot. The top panel of Figure 1 (which will be discussed in more detail throughout Section 2) highlights the overshoot regions in red. After the reflected particles have crossed back upstream of the "ramp," they continue to generate stronger electromagnetic fields by the same processes that generate the overshoot (Gedalin, 1996). These new fields plus the original plasma population's extant magnetic and motional electric fields accelerate the backstreaming particles while continuing to reflect them -this time back downstream toward the shock again (Gedalin, 1996; Woods, 1971). This denser immediately upstream region with slightly increased magnetic field is called the "foot" (Wong, 1968). The foot is much more difficult to accurately characterize automatically due to its small intermediate height nature, and so this study will do nothing with it other than attempt to identify regions upstream from the foot.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example magnetic field shock profile at 2015-01-29T17:10:12, successive partitioning of interval into sub-intervals, and panels for interval identification. (a) Magnetic field profile during this interval (black), hyperbolic tangent fit (yellow), overshoot region (red), “upstream measuring interval” (green), and “downstream measuring interval” (cyan). (b) Color-coded REGion IDentification (Halekas et al., 2017) of the local plasma: pristine solar wind (blue), disturbed solar wind/magnetosheath (cyan). (c) Initial division of time series into relevant sub-intervals, partitioned by “A” and “B” flags. (d) Automatically calculated refined division of relevant time range, with “F” denoting a conservative zeroth order estimate of the bottom of the shock foot (orange). (e) The boxcar averages $\langle B \rangle_{2\tau_{ci}}$ (green), $\langle B \rangle_{3\tau_{ci}}$ (red) and $\langle B \rangle_{4\tau_{ci}}$ (cyan), as well as horizontal lines of heights $B_D^{Fit}$ and $B_U^{Fit}$. (f) “B_meanMode”: the result of averaging together the components of each of the boxcar averages in (e) then taking a $\tau_{ci}$ width boxcar average of the result; Horizontal lines as above, and the orange highlighted region denotes the zeroth order foot estimate. (g) $(B_meanMode - B_U^{Fit})/B_D^{Fit}$ within the fitted interval upstream of the foot. (h) The standard deviation between the boxcar averages of (e), normalized to the mean of the components of (e), which is then boxcar averaged over $\tau_{ci}$ to smooth, and then boxcar averaged again over measuring interval width $T = \tau_{ci}[60s/\tau_{ci}]$ to find the center (yellow vertical line) of the interval (green) with the minimum average value.</figcaption></figure>

Quasiparallel shocks, by contrast, are much less well studied on account of being more difficult to characterize with a limited number of parameters. As $\theta_{BN}$ decreases into the quasi-parallel regime, the shock transition can become not so much a uniform, coherent, step-like structure with a single shock surface as a patchwork of smaller turbulent structures (Schwartz & Burgess, 1991). Due to the closer alignment between the “nominal” shock normal and the upstream magnetic field lines, suprathermal charged particles may more easily backstream along these field lines—having either been reflected at the shock transition (Gosling et al., 1989) or even diffused from
the downstream region following heating (Kennel et al., 1985). The backstreaming of suprathermal ions generates instabilities in the “upstream” plasma, in turn inducing an ultra-low frequency wavefield (τ ~ 10 s) “upstream” of the shock (Schwartz, 1991; Schwartz & Burgess, 1991). There are a variety of non-linear instabilities that can form structures in the foreshock. Some waves in this wavefield can strongly increase in amplitude whereupon they are referred to as Short Large Amplitude Magnetic Structures (Burgess et al., 2005; Schwartz, 1991; Schwartz & Burgess, 1991).

At the most macroscopic level, the Martian bow shock can be roughly characterized as an axially symmetric (i.e., rotated about the abberated $X_{MSO}$ axis) conic section—the calculation of average conic parameters which we follow are given in Trotignon et al. (2006). While this may be accurate to lowest order, more precise measurements have shown the presence of asymmetries. Mars has crustal magnetic fields—hypothesized to be vestigial remnants of a former intrinsic magnetic field (Acuña et al., 1998). Previous studies have shown that these fields are stronger in Mars' southern hemisphere (Acuña et al., 1998; Connerney et al., 2015), and that the bow shock is “pushed” to higher altitudes in the vicinity of these fields (Edberg et al., 2008; Garnier et al., 2022b; Gruesbeck et al., 2018; Hall et al., 2016)—thereby generating a southern asymmetry. The standoff distance has also been found to expand with solar EUV (Edberg et al., 2008; Garnier et al., 2022a, 2022b; Gruesbeck et al., 2018; Halekas et al., 2017; Hall et al., 2016), but decrease with solar wind dynamic pressure and solar wind magnetosonic Mach number $M_{ms} = |V_{SW}|/(C_S^2 + V_A^2)^{1/2}$ (Garnier et al., 2022a; Gruesbeck et al., 2018; Halekas et al., 2017; Hall et al., 2016). The bow shock standoff distance is smaller in the quasi-parallel region (Ma et al., 2004; Vignes et al., 2000; T.-L. Zhang et al., 1991), corresponding to “dawn” side (i.e., at negative $Y_{MSO}$) of the bow shock (Gruesbeck et al., 2018)- a result of the predominantly Parker spiral nature of the solar wind magnetic field. A much lower timescale bow shock variation -connected to the EUV dependence-arises from seasonal variation. However, it has not been as well studied how these various parameters affect the structure of the Martian bow shock.

## 1.2. Martian Atmospheric and Seasonal Variation

Like Earth, the rotational axis of Mars is tilted, and thus the planet has seasons. Mars seasons can be defined in terms of Mars Solar Longitude $L_S$ (with $L_S = 90°$ corresponding to southern winter solstice), and such that each season covers a 90° interval with southern (northern) autumn (spring) starting at $L_S = 0°$ (Zurek & Martin, 1993). On the other hand, Mars has a much more elliptical orbit than Earth, with $\epsilon_M = 0.0935$ compared to Earth's $\epsilon_E = 0.0167$ (Williams, 2021). Considering that solar radiation drops off as $r^{-2}(L_S) \propto (1 + \epsilon_M \cos(L_S - 251°))^2 \approx (1 + 0.1 \cos(L_S - 251°))^2$—such that the solar radiation can vary by ±20%—one might expect that radiation driven atmospheric phenomena should vary not insignificantly with season. In fact, previous studies have found exactly this.

Zurek and Martin (1993) observed that Martian global dust storms peak following perihelion ($L_S = 251°$), during southern spring and summer. The rate of hydrogen escape through the hydrogen exosphere has also been observed to peak during approximately this same period, increasing from the minimum rate by over an order of magnitude (Bhattacharyya et al., 2015; M. S. Chaffin et al., 2014; Halekas, 2017; Rahmati et al., 2018). It has been theorized (M. Chaffin et al., 2017; Heavens et al., 2018) that these phenomena are connected: dust storms transport water molecules to higher altitudes, whereupon they are photodissociated to further populate the hydrogen exosphere.

These are not the only observed seasonal variations in the Martian system: Halekas (2017) found that the trend of the escape rate and upstream column density has a maximum at approximately southern summer solstice ($L_S = 270°$), an order of magnitude greater than the minimum column density near aphelion. The spatial extent of the bow shock itself varies with season: As Mars approaches the sun, enhanced EUV heats the atmosphere, thereby inducing thermal expansion and increasing the exosphere thickness—and thus bow shock height. Increased mass loading during higher EUV may also contribute to the increase in bow shock height (Hall et al., 2016).

The exospheric neutral hydrogen atoms discussed above can be ionized, whether by solar EUV or charge exchange with solar wind ions, and then mass load the solar wind as pick-up ions. Once picked-up, these planetary ions can be accelerated by the solar wind's motional electric field $E = -V_{SW} \times B_{SW}$ and begin to undergo cycloidal motion. Those pickup ions produced close enough to the shock to reach the discontinuity early in their gyration will have insufficient normal velocity to avoid being reflected (Dubinin et al., 2006). Such reflected protons may repeatedly bounce within the foot-ramp-overshoot structure (Gedalin, 1996) until they are sufficiently electrically
accelerated to pass fully downstream of the shock. Given that the extensive seasonal changes in the Martian environment are known to affect the size of the bow shock, we seek to determine whether the introduction of additional planetary ions into the upstream flow changes the magnetic structure of the bow shock.

## 1.3. MAVEN

The primary source of data for studying plasma in the Martian environment at this time is the Mars Atmosphere and Volatile Evolution (MAVEN) spacecraft (Jakosky et al., 2015), which arrived in September 2014, and began collecting data on 15 November 2014. MAVEN has a comprehensive suite of instruments for studying the magnetic field, as well as the ion and electron plasma, and other properties of the local environment. Magnetic field data is collected with the magnetometer (MAG) instrument (Connerney et al., 2015). For our purposes, we use the l2 data, at 1 s cadence.

The primary source of ion plasma data is the Solar Wind Ion Analyzer (SWIA) instrument (Halekas et al., 2015). SWIA data has two primary modes, Fine (SWIF) and Coarse (SWIC) for use in the solar wind and Martian magnetosphere, respectively. From the energy distribution data, 8 s cadence ion moments (n, v, T, P) can be calculated. As the energy distributions possess two peaks for H⁺ and He⁺⁺ species, separate temperature moments Tproton and Talpha are calculated for each (Halekas et al., 2017). As we are primarily interested in the effect of protons (as the dominant species in the solar wind), we henceforth use only Tproton for temperature measurements. The effect of alpha particles on the bow shock structure is neglected.

Along with ion data, MAVEN also has a module for electron measurements. The Solar Wind Electron Analyzer (Mitchell et al., 2016) instrument measures electron energy and angular distributions. In particular, we utilize electron number density n_e and temperature T_e moments calculated from SWEA's energy spectrogram (SPEC) data, with spacecraft potential corrected for. In April 2022, it was determined that the SWIA density moments were uniformly offset by a factor of 0.7, and so the entire SWIA moment data was re-calibrated. The SWEA calibration was initially set based on SWIA measurements, and was not yet re-calibrated at the time our study was conducted. In order to recalibrate our SWEA data to that of SWIA, we manually divided the electron density moments by 0.7 to match the ion counterparts. There also existed occasional issues with the electron moment data. For instance, occasions exist where the spacecraft potential is poorly measured or not measured at all, leading to non-existent or NAN moments. Additionally, on 16 November 2019, due to blanking out of low energy measurements due to an instrumental issue, the overall moments began to be computed to have values reduced by a factor of approximately 10. We will handle these issues by rejecting any data point which measures NAN for upstream electron parameters, and by choosing 15 November 2019 as the final day of our run of MAVEN data.

In this study, we attempted to determine what, if any, seasonal dependence the Martian bow shock structure may have by analyzing MAVEN data. Such variations are necessarily statistical in nature, requiring a very large number of data points per season. As such, we built an algorithm to run through the data, 1 day at a time, and try to locate, fit, measure, and characterize the bow shock crossings in each. This process is repeated for the first three and a half Martian years of MAVEN data, from 15 November 2014 to 15 November 2019.

## 2. Methodology

In any given day, we attempt to partition the time-series data into solar wind-magnetosheath (on the inbound side—magnetosheath-solar wind on the outbound) regions, fit each shock interval to hyperbolic tangent functions, locate and measure the overshoot, find small representative upstream and downstream intervals in-which average measurements can be made, and then calculate a number of representative quantities. This takes anywhere from 650 to 1100 s, and thus to save computation time, we start by making a number of checks to determine whether any measurements will be possible.

1. Halekas et al. (2017) provided an algorithm for identifying the pristine solar wind, “mvn_swia_regid” (hereby “REGID”)—publicly available in the SPEDAS code base (Angelopoulos et al., 2019), which assigns Region ID's given time averaged ion velocity magnitude $V_{MSO}^{ion}$ [km/s], ion temperature magnitude $T_{ion}^{e}$ [eV], number density $n^{ion} [cm^{-3}]$, altitude and MSO coordinates [km], and magnetic field magnitude $B[nT]$, as well as the temporal magnetic standard deviation $\sigma_B[nT]$. It assigns IDs for pristine solar wind (“REGID = 1” where altitude >500, $V_{ion}^{e} > 200$, $\sqrt{T_{ion}^{e}/V_{ion}^{e}} < 0.012$, and $\sigma_B/B < 0.15$), disturbed solar wind/magnetosphere
("REGID = 2" where $V^{\text{ion}} > 200$, altitude >300 and either $\sqrt{T^{\text{ion}}}/V^{\text{ion}} > .015$ and/or $\sigma_B/B > 0.25$), the "ionosphere" (REGID = 3 and 4, see code), and "magnetotail lobe" (REGID = 5, not used here) based on MAVEN moment data. If this algorithm detects less than 40 s (10 array elements) where the spacecraft is in the pristine solar wind, then the day is skipped.

2. If SWEA moment data cannot be generated, or there is less than 400 s (100 array elements) where the electron density is finite or the minimum computed electron number density is greater than 1,000 cm⁻³ (In which case it failed to calculate moments in the upstream solar wind), then we skip to the next day.

## 2.1. Identifying Precise Solar Wind-Shock-Sheath Intervals

In order to identify a day's bow shock crossings and respective local solar wind and magnetosheath intervals, we must start by partitioning it in such a way as to locate intervals of adjacent solar wind-magnetosheath regions. To lowest order, the REGID algorithm does exactly this, designating temporal regions as being within (pristine) solar wind, magnetosheath, or other sub-sheath regions. Being designed to identify pristine solar wind, the REGID algorithm does not however switch between regions perfectly cleanly: As can be seen in Figure 1b, numerous instances exist where it determines that MAVEN has left the pristine solar wind and entered the magnetosheath/disturbed solar wind for a few time-series elements before returning, or vice versa. It may also oscillate between the two at the boundary. Similar phenomena occur at the low-altitude sheath boundary. Further, there exist array elements where no region ID has been chosen. In order to use the REGID data, we must “clean” it up by filling in null values and smoothing out perturbations from the overall behavior (for compactness, the details are given in Text S3 in Supporting Information S1). Once this data has been cleaned, we can assign pointers for initial zeroth order estimates of when the spacecraft passes the shock and a rough timestamp for a downstream edge of consideration. We denote these “A” and “B” respectively, and their positions for the 2015-01-29T17 shock crossing are shown in Figure 1c. Not shown is the initial upstream boundary for this solar wind-magnetosheath interval, corresponding to the adjacent “A” flag. Note that outbound shocks, where magnetosheath-solar wind intervals are considered, are done separately by means of time reversing our data before conducting identical operations. For simplicity, the rest of this methodology will be written in terms of inbound shocks only.

As has been previously noted, REGID makes no distinction between the magnetosheath and non-pristine solar wind—thus the boundary between the calculated “solar wind” and “sheath” periods will be upstream of the shock itself. Further, no IDs are provided for the pileup region—our cleaning procedure places “B0” halfway between the closest “magnetosheath” and “ionosphere” identifiers as an initial anchor—thus we may be including the pileup region within our interval. As such, it is necessary to shrink our interval to something representative of the local solar wind and sheath at either side of the shock crossing. The flag refining procedure is an iterative process. The zeroth order magnetic field $B_0$ and suprathermal ion density $n_0^{\text{ion}}$ profiles obey $\text{sign}\left(\frac{dB_0}{dt}\right) \times \text{sign}\left(\frac{dn_0^{\text{ion}}}{dt}\right) = 1$ in the sheath and $\text{sign}\left(\frac{dB_0}{dt}\right) \times \text{sign}\left(\frac{dn_0^{\text{ion}}}{dt}\right) = -1$ in the pileup region, so we first shift our “B” flag upstream while the latter is true. We next refine our guessed shock location “A.” While the shock ramp should nominally correspond to the global maximum in $\frac{dB_0}{dt}$ within our interval, there exist certain scenarios in which our interval still extends far enough downstream to include induced magnetic fields where a larger maximum in $\frac{dB_0}{dt}$ may be present. As a work around, we update “A” by iterating downstream over large local maxima (greater than 25% of the global maximum of our interval) in $\frac{dB_0}{dt}$ until we reach a local maximum greater than the subsequent maximum. Using this information we can identify a zeroth order guess for the bottom of the foot, to which we assign a flag we denote “F.” Note that this flag is a conservative estimate of where the upstream is roughly quiescent, and will therefore tend to include turbulent foreshock regions as part of the “foot.” We go into more detail within our Supporting Information S1, in Text S4. Figure 1d shows the final placement of these flags. It should be noted that whether or not “Foot” is further upstream than the initial “A” flag varies with the shock in question.

## 2.2. Fitting the Shock

To lowest order, a well-behaved shock can be expected to have a shock profile of a step function plus a number of higher order perturbations near the discontinuity. Depending on the upstream parameters, this discontinuity's ramp can be broader or steeper. As such, to lowest order, the profile can be modeled as a sigmoid function.

In particular, we attempt to model our shock profiles to hyper-tangent functions $f(t) = m_0 \tanh(m_1(t - m_2)) + m_3$. Here, $B_D^{\text{FIT}} = m_3 + m_0$ and $B_U^{\text{FIT}} = m_3 - m_0$ are the best estimates of the downstream and upstream equilibrium
magnetic field strengths respectively, $1/m_1$ controls the steepness of the sigmoid ramp, and $m_2$ is the best estimate for the time of the center of the ramp. IDL's curvefit function (Curvefit, 2013) utilizes a $\chi^2$ gradient descent algorithm based on initial input fit parameters. As such, it is possible for the curvefit function to generate a false result with ramp centered at the wrong location due to reducing the $\chi^2$ to a local (false) minimum. Thus, it is imperative to ensure that our initial fit parameters are as close to correct as possible. This is done by pulling downstream our upstream boundary until it is at most 10 min upstream of our calculated foot—then attempting to more precisely set our “A” flag to the center of our ramp, and then pulling upstream our “B” flag until it is no farther than 40 min downstream of our “A” flag.

Once this is done we can calculate initial fit parameters $m_0, m_1, m_2$, and $m_3$. From our output fit—should it be verified as a valid result, we can calculate the time of the center of the ramp exactly. Once we've done this, we can attempt to refit our shock within a more local region of 2 min upstream of the foot to 7 min downstream of the center of the ramp. Unfortunately, our above “A” refinement algorithm is not always sufficient, and there is not a single optimal algorithm for correctly refining the “A” flag in all scenarios. As such we also try multiple combinations of backup subroutines and choose at each shock the outputted fit with the lowest reduced $\chi^2$. We leave the descriptions of the various fitting routines to our Supporting Information S1. With this new set of fitted shock profiles we have the center of our ramps (hereby the “shock location”), asymptotic “upstream” and “downstream” magnetic field magnitude values (hereafter notated $B_U^{Fit}$ and $B_D^{Fit}$), and temporal ramp widths. Such values can be used for more detailed measurements. Figure 1a shows the automatically generated fit for 2015-01-29T17:10 in yellow.

## 2.3. Determining the Overshoot and Representative Upstream and Downstream Intervals

The behavior of the plasma near shocks is not entirely discontinuous: neither $B, V_s, T_s, n_s$ or any other parameter (with a few exceptions) actually follow step or sigmoid functions. Further characterization requires identification of the overshoot (for purposes of measuring $B_{max}$), and where the overshoot and subsequent undershoot end for the purpose of identifying a metastable sheath interval. The structure of collisionless shocks depends on the upstream densities and pressures, as well as the full vectorial magnetic field. To consistently measure these values, it is necessary to identify downstream and upstream intervals over which average values can be calculated. In order to ensure a long enough sample size to average out high frequency perturbations—in particular to wash out those associated with ion cyclotron motion, but not too long to avoid including isolated high amplitude perturbations uncharacteristic of the overall environment, we choose measuring intervals of duration $T = \tau_{ci}[60s/\tau_{ci}]$, the smallest integer multiple of the local upstream ion cyclotron period (calculated from a minute average immediately upstream of “F”) that is greater than a minute. We will use the shock profile in Figure 1 as an example case. The downstream portion of this—finding the end of the overshoot and identifying an interval—are conducted together, and performed twice in case the first time fails.

When locating the overshoot, undershoot, and metastable sheath region where an interval can be chosen with which to average over downstream data, our first attempt at doing so involves automating the methodology for manually locating these intervals as laid out in Burne et al. (2021). We identify the start of the overshoot as where “B_meanMode”—the $\tau_{ci}$ (as calculated from a 1 min interval just upstream of the foot) boxcar average of the element wise average (panel f) of the $2\tau_{ci}^{upstream}$, $3\tau_{ci}^{upstream}$, and $4\tau_{ci}^{upstream}$ width boxcar averages (panel e) of the field—is greater than $B_D^{Fit}$ for the first time following the shock, the end of the overshoot as when the array passes below again and then the undershoot ends when it returns back above. We then set our downstream measuring interval immediately afterward. If any of these fail, we attempt a number of other similar algorithms based upon which flag failed to trigger. If all of the above fail, we resort to finding when $B_meanMode$ reaches a peak and subsequent valley—setting the overshoot in between the shock and valley with $B_{max}$ equal to the maximum 1 s resolution magnitude in the overshoot in either case, and setting our downstream interval to after our valley has ended. Figure 1a shows the automatically determined overshoot region, colored in red, and automatically determined downstream representative interval in cyan. We go into more detail within Text S6 in Supporting Information S1.

The identification of a representative upstream interval (of the same duration as the downstream counterpart) is a somewhat complex procedure. Figure 1 gives an illustrative example of this. For the purposes of ensuring our upstream measurements characterize local properties, we require that our measuring interval must be within the interval which has been fitted. To ensure our measurement will be of roughly pristine solar wind, we also demand
that the chosen interval must be upstream of the “Foot” flag. To zeroth order, and assuming no other structure in the plasma, B should asymptotically decrease with distance from the shock and approach a value of $B_D^{Fit}$. Thus our first check is whether the time derivative of a zeroth order estimate of B (not shown) reaches a minimum at the upstream boundary. As this has failed, we next utilize the $\langle B \rangle_{2\tau_{ci}}$, $\langle B \rangle_{3\tau_{ci}}$ and $\langle B \rangle_{4\tau_{ci}}$ boxcar averages (panel e) to search for when the field is most quiescent. Upstream of the shock, pick-up protons generate a high amplitude wavefield with period $\tau_{ci}^{upstream}$ on top of the already extant upstream magnetic field. This wavefield is a common feature of the near-Mars upstream environment (Mazelle et al., 2004). Where $\langle B \rangle_{2\tau_{ci}}$, $\langle B \rangle_{3\tau_{ci}}$ and $\langle B \rangle_{4\tau_{ci}}$ are in closest agreement for the longest interval should correspond to when the plasma is most quiescent. Conversely, they will begin to noticeably diverge in the presence of instabilities as well as immediately before the foot. We therefore calculate the standard deviation between $\langle B \rangle_{2\tau_{ci}}$, $\langle B \rangle_{3\tau_{ci}}$ and $\langle B \rangle_{4\tau_{ci}}$ and normalize by their mean in order to determine the amount of disagreement, then boxcar average over $\tau_{ci}^{upstream}$ to smooth over isolated one-off disagreements. To find where this reaches a minimum for the entire duration $\Delta t = T$—the length of our measuring interval—we boxcar average again over width $T$ (panel h). Since $(B_{meanMode} - B_U^{Fit})/B_U^{Fit}$ (panel g) is on average less than 25% within this interval (and therefore agrees sufficiently with our fit) and the downstream edge of the interval has an averaged standard deviation (panel h) below our chosen threshold (as an additional check that the interval does not include the foot) in this case, we set this as our upstream measuring interval. In Figure 1a, the interval is highlighted in green. Further detail can be found in Text S7 in Supporting Information S1.

## 2.4. Calculating Representative Quantities

Once representative upstream and downstream intervals have been identified, we can calculate the shock normal vectors and average upstream parameters. For compactness we will henceforth write the time average of a scalar quantity $q$ measured in the upstream (downstream) interval of interest as $\langle q \rangle_U$ ($\langle q \rangle_D$), and the standard deviation in the interval as $\sigma_{qU}$ ($\sigma_{qD}$). For a vectorial time-series quantity $A$, we denote the component wise time-average over the upstream (downstream) interval of interest as $\langle A \rangle_U$ ($\langle A \rangle_D$).

Schwartz (1998) provides methods of calculating a shock normal vector from single spacecraft measurements of the magnetic field **B** and bulk flow velocity **V** in an arbitrary coordinate system. We compute the three mixed-mode methods $\hat{n}_{MX1}$, $\hat{n}_{MX2}$, and $\hat{n}_{MX3}$—using the fine and coarse ion velocity $\mathbf{v}^{ion}$ for $\langle V \rangle_U$ and $\langle V \rangle_D$—and then take the normalized component-wise average of the three in order to obtain a more consistently accurate normal vector $\hat{n}_{AVG}$. Once the average shock normal vector has been found we can calculate the shock normal angle $\theta_{BN}$ in terms of $\hat{n}_{AVG}$ and $\langle \mathbf{B} \rangle_U$. In the example case in Figure 1 the calculated shock normal angle is $\theta_{NB}^{calc} = 82.855904^\circ$. This methodology does a very good job of locating the global scale shock normal—that is, a normal to the large-scale average shock surface calculated by Trotignon et al. (2006), since it relies on average quantities well upstream and downstream of the shock. Figure 2a shows the shock normal vectors projected onto the locations of their respective shock crossing locations in cylindrical MSO coordinates, along with the shock conic section and spacecraft trajectory, for 2015-01-29. Figure 2b aggregates all of our data set's shock crossing locations together in a cylindrical projection, coloring the points based on the dot product between the measured (average) normal vector and the normal vector calculated from the Trotignon et al. (2006) shock conic at the point on the surface closest to the measured shock crossing location.

Once this is done, we can calculate other derived upstream parameters such as the plasma beta (ion, electron, and total), total sound speed, Alfven speed, and fast magnetosonic Mach number $M_{fms}$. For consistency, all derived parameters are defined as functions of time-series averaged raw quantities $\langle n_{ion} \rangle_U$, $\langle n_e \rangle_U$, $\langle \mathbf{v}^{ion} \rangle_U$, $\langle \mathbf{B} \rangle_U \equiv \langle |\mathbf{B}| \rangle_U$, $\langle \mathbf{B} \rangle_D \equiv \langle |\mathbf{B}| \rangle_D$, $\langle \mathbf{B} \rangle_U$, $\langle \mathbf{B} \rangle_U$, $\langle T_{proton} \rangle_U$, and $\langle T_e \rangle_U$. Thus for instance the fast magnetosonic Mach number $M_{fms}$ is calculated as:

$$M_{fms} = \frac{\hat{n}_{AVG} \cdot \langle \mathbf{v}^{ion} \rangle_U}{\left[ 0.5 \left( C_s^2 + V_A^2 + \sqrt{(C_s^2 + V_A^2)^2 - 4C_s^2 V_A^2 \cos^2 \theta_{NB}} \right) \right]^{1/2}} \quad (3)$$

where $V_A^2 = \langle \mathbf{B} \rangle_U^2 / (\mu_0 m_p \langle n^{ion} \rangle_U)$ and $C_s^2 = \langle (T_e)_U + 5/3 \langle T_p \rangle_U \rangle / m_p$. From the corrected equations in chapter 6 of (C. T. Russell et al., 2016), we can calculate the first critical Mach number $M_{crit}$ from our measured $\theta_{NB}$ and $\beta = \beta_i + \beta_e$. We make the necessary simplification of assuming that the MSO frame is the bow shock rest frame:
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) The calculated normal vectors (arrows) for detected shock crossings on 2015-01-29 projected onto their locations within the MAVEN spacecraft trajectory (black) for that day. The brown arc is the Trotignon et al. (2006) shock conic section. (b) Every approved data point projected into their crossing location in cylindrical MSO coordinates. Points are colored by dot product of $\hat{n}_{AVG}$ with the conic surface normal $\hat{n}_{conic}$ at the closest point on the conic to the shock crossing location.</figcaption></figure>
calculations of the frame velocity of the shock assume a fixed downstream density, but the ion density in the Martian magnetosheath is rarely sufficiently stable that $\langle n^{ion} \rangle_D$ can be justifiably considered characteristic of the entire magnetosheath.

## 2.5. Data Point Verification

The fits and interval identification routines do not always work correctly. More specifically, downstream of the overshoot the sheath does not always converge to a relatively stable magnitude for sufficient lengths of time within our window of consideration such that a characteristic region of measurement can be accurately identified. If the mean downstream field $\langle B \rangle_D$ in our chosen interval varies significantly from the asymptotic fit value $B_D^{Fit}$, then we throw the data point out. In other cases, primarily at low $\theta_{NB_1}$, there does not exist a well-defined shock or magnetosheath. Here too, we throw out data points. The shock is not entirely stationary: it can ripple back and forth across the spacecraft thus causing the latter to briefly cross between the solar wind and sheath multiple times before settling on the other side. If our overall fitted interval covers multiple shock crossings such that the upstream fit, upstream measuring interval, downstream fit, and downstream measuring interval do not all correspond to the same crossing, then the normal vector is not being precisely measured and thus we have to throw the point out. Other causes for throwing out points are entirely the result of failure of the algorithm: The overshoot (and $B_{max}$) can be improperly identified. Related to this are scenarios where the downstream measuring interval has failed to pass downstream of the overshoot or has passed too far and entered the beginning of a crustal field or a non-local sheath region. On the upstream side, we manually reject the data point if our algorithm fails to position our representative.

We also throw out points automatically based on sanity check based triggers: If SWE moment data cannot be computed within our upstream measuring interval, then the point is automatically thrown out. Similarly, we require $0 < \beta \le 20$, $1 \le M_{fms} \le 10$, $7 < \langle n_e \rangle_U / \langle n_i \rangle_U < 2$, $1.5 \le B_D^{Fit} / B_U^{Fit}$, $\langle B \rangle_D / \langle B \rangle_U \le 7$, $\hat{n}_{AVG} \cdot \hat{n}_{CONIC} \ge .5$, and $(\langle V^{ion} \rangle_D \cdot \hat{n}_{AVG}) / ((\langle V^{ion} \rangle_U \cdot \hat{n}_{AVG}) > 0$. In order to ensure our fits and measuring intervals do not output values that deviate too far, we further demand $|B_U^{Fit} - \langle B \rangle_U| \le 1.3$, $2|B_U^{Fit} - \langle B \rangle_U| / (B_U^{Fit} + \langle B \rangle_U) \le .25$, $|B_U^{Fit} - \langle B \rangle_U| \le \max(3\sigma_{B_U}, 1.3)$, and $|B_D^{Fit} - \langle B \rangle_D| \le 3\sigma_{B_D}$.

The SWIA instrument is calibrated once a month within a magnetosheath crossing. SWIA measurements during these events are entirely non-physical. In order to throw away any such events, we reject all crossings where the measured density jump $\langle n_{ion} \rangle_D / \langle n_{ion} \rangle_U > 50$.

It should be noted that the attached Supporting Information S1 describes each of the above steps in our algorithm in more extensive detail, including supplying the relevant filenames to be found in Fruchtman (2023a).

# 3. Results

After every unreliable data point has been thrown out, 3846 remain out of the 7056 total fitted crossings which reached the point of needing to be thrown out as described above. Using our aggregated data, we can look at variations in the bow shock structure at multiple system scales. In particular, we will look at the global scale spatial variation in the bow shock, the MHD scale magnetic jump, and the kinetic scale overshoot. For ease of reference, all parameter definitions and correlation coefficients discussed below are tabulated in the supplement.

## 3.1. Global Scale Spatial Variation in the Bow Shock

Starting at the global scale, we may observe spatial variations in the bow shock by plotting the locations of our identified shock crossings as functions of altitude and solar zenith angle SZA, and coloring each point based on the associated value at the shock for various parameters of interest. We do so in Figure 3. We briefly define the shock extent (the shock’s “latus rectum”) as

$$ \Delta = \sqrt{(Altitude^2 - X_{MSO}^2) + (X_{MSO} - X_0)^2 + \epsilon(X_{MSO} - X_0)} \quad (4) $$

where $X_0 = .6R_{mars}$ and $\epsilon = 1.026$ are conic parameters calculated in Trotignon et al. (2006). A $\Delta$ contour line is drawn for the $\Delta_{Trotignon} = 2.081R_M$ surface. In panel (a), we show that the shock normal upstream ion velocity $V_N \equiv \langle V^{ion} \rangle_U \cdot \hat{n}_{AVG}$ strongly varies with SZA ($r(V_N, SZA) = 0.525$) as should be expected: the solar wind flows approximately anti-sunward, therefore crossing the shock at more oblique angles toward the flank.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** We plot our data as functions of altitude and SZA. Points are colored by (a) $\langle V^{ion} \rangle_U \cdot \hat{n}_{AVG}$, (b) Degrees from the southern summer solstice ($L_s = 270^\circ$) when shock is maximally expanded. (c) The solar wind Magnetosonic Mach number, (d) $\theta_{BN}$ (e) $M_{fms}$, (f) $M_A$, (g) critical ratio, (h) Magnetic Jump ratio, (i) $\log_{10}(A)$. The black contour line shows the Trotignon et al. (2006) conic.</figcaption></figure>

In panels (b)–(d), we show successful replication of results from previous studies: Points in panel b are colored by degrees $\Delta L_s$ from the $L_s = 270^\circ$ southern summer solstice (hereby $\Delta L_s(270)$). We see that the shock extent does in fact increase as Mars approaches the southern summer solstice ($r(\delta L_s(270), \Delta) = -0.357$). This is entirely expected: the southern summer solstice is shortly after perihelion ($L_s = 250^\circ$) and the total daytime EUV entering the atmosphere reaches a maximum. Increasing the temperature of the atmosphere causes the volume of the atmosphere and the rate of mass loading (and thus the shock) to expand. Panel (c)'s plot of the solar wind magnetosonic Mach number

$$M_{ms} = |\langle V^{ion} \rangle_U| / \sqrt{C_S^2 + V_A^2} \quad (5)$$
shows an extremely strong anti-correlation with shock extent ($r(M_{ms}, \Delta) = -0.471$), as previously reported (Farris & Russell, 1994; Peredo et al., 1995). Consistent with the results of T.-L. Zhang et al. (1991), we see in panel (d) that $\theta_{NB}$ possesses a just barely visible variation with $\Delta$ ($r(\theta_{NB}, \Delta) = 0.254$): at lowest $\Delta$, there is a thin layer that is much more quasi-parallel on average than the rest of the data set.

In panels (e)-(i) we show our new findings. We see in panel (e) that $M_{fms}$ decreases with SZA ($r(M_{fms}, \text{SZA}) = -0.477$). This should not be surprising, as its numerator is simply $V_N$. Unexpectedly, this does not hold true for $M_A = V_N/V_A$, which has a much stronger correlation with shock extent $r(M_A, \Delta) = -0.372$ than it does SZA ($r(M_A, \text{SZA}) = -0.187$) or altitude ($r(M_A, \text{alt}) = -0.333$). Similarly, we see in panel (g) a stronger anticorrelation between $M_{fms}/M_{crit}$ and shock extent ($r(M_{fms}/M_{crit}, \Delta) = -0.409$) or altitude ($r(M_{fms}/M_{crit}, \text{alt}) = -0.476$) than with the also visible SZA ($r(M_{fms}/M_{crit}, \text{SZA}) = -0.349$) dependence.

In terms of the magnetic structure itself, we see first in panel (h) that the jump $J^{\text{Fit}} = B_D^{\text{Fit}}/B_U^{\text{Fit}}$ is clearly anti-correlated with SZA ($r(J^{\text{Fit}}, \text{SZA}) = -0.552$). This should not be surprising: the jump is a consequence of conservation law requirements for deflecting the magnetized solar wind. Less deflection is necessary toward the flanks of the magnetosheath (high SZA) than near the subsolar point. Unexpectedly, we see in panel (i) that overshoot amplitude is clearly anticorrelated with shock extent ($r(\log_{10}(A), \Delta) = -0.437$) but negligibly with SZA ($r(A, \text{SZA}) = 0.088$) or altitude ($r(\log_{10}(A), \text{alt}) = -0.193$).

## 3.2. Deviation of the Magnetic Jump From Rankine-Hugoniot Predictions

Before discussing the seasonal independence of the magnetic jump, we start by considering the level of agreement of $J^{\text{Fit}}$ with the Rankine-Hugoniot predictions for the jump $J^{\text{RH}} = (B_D/B_U)^{\text{RH}} = F(M_{fms}, \theta_{BN}, \beta)$ as calculated using the (corrected) equations in C. T. Russell et al. (2016). Far from the kinetic scale behavior across the discontinuity, the solar wind plasma should reduce to MHD behavior and thus to the lowest order of ideal MHD should obey the RH equations. We do see from Figure 4 that $J^{\text{Fit}}$'s median dependence on upstream conditions obeys the predictions to approximately first order, though it diverges as $\theta_{NB}$ decreases through the quasi-parallel regime and $\beta$ approaches 10. In particular, $J^{\text{Fit}}$'s dependence on $M_{fms}$ is on average in agreement with $J^{\text{RH}}$ to at least first order throughout its range. The median lines for the $M_A$ dependence too are in very good agreement, for $M_A < 10$. The disagreements with Rankine Hugoniot predictions may be of some concern, so we will discuss the various apparent inconsistencies in turn.

It is evident that our data's $\theta_{BN}$ dependence disagrees with the Rankine-Hugoniot predictions in the quasi-parallel ($Q_\|, \theta_{BN} < 45^\circ$) regime. We hypothesized that this could be explained by survivorship bias: $Q_\|$ shocks tend to be much less well structured than quasi-perpendicular ($Q_\perp, \theta_{BN} > 45^\circ$) shocks, and thus are less likely to be successfully fitted. Of those that are, poorly structured shocks are typically rejected during the (automatic and manual) quality assurance phase due to issues with properly identifying the shock location, overshoot, $B_{\text{max}}$ and/or local approximate sheath. As such, it would stand to reason that the only surviving shocks calculated to be quasi-parallel are those more similar in structure to quasi-perpendicular shocks, and thus have greater jumps. Counting quality assurance rejected shocks however, we find 1–2835/5026 = 2191/5026 = 43.6% of quasi-perpendicular shocks were rejected compared to 1–1010/2039 = 1029/2039 = 50.5% of quasi-parallel shocks, suggesting that this may not account for the discrepancy.

Another possible set of explanations for the quasi-parallel divergence can be considered by noting that $\hat{n}_{AVG} \cdot \hat{n}_{conic}$ itself breaks down at low $\theta_{BN}$. For most cases our calculated shock normal is very close in direction to the closest positioned global shock conic surface normals: $\hat{n}_{AVG} \cdot \hat{n}_{conic} > .94$ for most shocks. However not only is this not uniformly true, but- as can be seen in Figure 2b—there is a clear dependence of the shock normal dot product on standoff distance ($r(\hat{n}_{AVG} \cdot \hat{n}_{conic}, \Delta) = 0.247$), with a clear drop in $\hat{n}_{AVG} \cdot \hat{n}_{conic}$ at the inner edge. This is the same spatial trend as $\theta_{BN}$, and in fact there is a small but significant (0.000000) Spearman rank correlation between these terms: $\rho(\theta_{BN}, \hat{n}_{AVG} \cdot \hat{n}_{conic}) = 0.431$. The quasi-parallel breakdown may potentially be explained by whatever has caused this correlation.

The first possible explanation for the $\hat{n}_{AVG} \cdot \hat{n}_{conic}$ correlation is that many of the shocks determined to have low $\theta_{BN}$ have been incorrectly calculated and are in fact quasi-perpendicular. This discrepancy has more wide reaching implications as it invalidates our calculations for quasi-parallel $M_{crit}$ and $M_{fms}$ (to second order) or in worst case scenario anything derived from $\hat{n}_{AVG}$, including $M_{fms}$ (to lowest order) and $M_A$. An argument for this can be made by considering an alternative set of derived quantities $\theta_{BN_{conic}}$, $M_{fms}^{conic}$, $M_A^{conic}$, and $M_{crit}^{conic}$ defined such that $\hat{n}_{AVG}$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** *J* = *B*<sub>D</sub>/*B*<sub>U</sub> versus *M*<sub>fms</sub> (top), *β* (row 2), *θ*<sub>BN</sub> (row 3) and *M*<sub>A</sub> (bottom). On the left we superimpose the predicted Magnetic jump calculated via Rankine-Hugoniot Jump Conditions *J*<sup>RH</sup> = *f*(*M*<sub>fms</sub>, *θ*<sub>BN</sub>, *β*) (orange, red median line) along with the ratio of asymptotic downstream and upstream magnitudes calculated from the fit *J*<sup>Fit</sup> = *B*<sub>D</sub><sup>Fit</sup>/*B*<sub>U</sub><sup>Fit</sup> (blue, black median line). On the right, we split each by southern season: spring (green), summer (red), autumn (gold), winter (blue), and total median (gray). Median lines use bin sizes *b*<sub>M<sub>fms</sub></sub> = 1, *b*<sub>β</sub> = 0.25, *b*<sub>log10β</sub> = 0.2, *b*<sub>θ<sub>BN</sub></sub> = 5°, *b*<sub>M<sub>A</sub></sub> = 1, with error bars showing the 0.25 and 0.75 quartile values in each.</figcaption></figure>

is replaced with $\hat{n}_{conic}$. The correlation coefficients of $M_{fms}^{conic}$ and $M_A^{conic}$ with $J^{Fit}$ increase significantly (see supplement for table of values) compared to the empirically calculated values, as do the correlation coefficients of $M_A^{conic}$ and $M_{fms}^{conic}/M_{crit}^{conic}$ with A. This is even true when considering only the subset where $\theta_{BN_{conic}} > 45^\circ$. This could be interpreted as suggesting that the conic is in some ways more accurate than our calculated normal vector: this alternative method allows us to more accurately predict the jump from these parameters. On the other hand, the correlation between shock normal angle and $J^{Fit}$ drops when calculated in this way. Furthermore, the Rankine-Hugoniot predicted jump is significantly less correlated with $J^{Fit}$ when using the conic normal. Thus, our disagreement with the conic normal at low $\theta_{BN}$ cannot entirely solve the problem of RH diverging at low shock normal angle.
A more plausible explanation is one that arises from physical processes: Low $\theta_{BN}$ shocks are very often non-planar and non-stationary (C. T. Russell et al., 2016). That is, they ripple and reform. As the Rankine-Hugoniot equations assume both planarity and stationarity, their resulting jump conditions no longer hold. Furthermore, as the non-planar shock ripples, the local surface normals $\hat{n}_{local}(x)$ fluctuate variably in quasi-perpendicular or even more quasi-parallel directions (Hao et al., 2016). As a result, the shock begins to take on a mix of quasi-perpendicular and quasi-parallel features. This includes the upstream and downstream fields in the vicinity. Our calculation of $\hat{n}_{AVG}$ involves time-averaging over ranges relatively far upstream and downstream of the shock crossing, to calculate $\hat{n}$ in three different methods before averaging together the components. Thus we have ultimately averaged over some upstream and downstream segments of the instability. Resulting $\theta_{BN}$ which are measured to be quasi-parallel may therefore correspond to shocks with features which conspire to result in overall quasi-perpendicular-like structures, and in general are non-planar and thus non RH-compliant. Similarly, shocks which are genuinely quasi-parallel at the MHD scale may have segments of their non-planar instabilities averaged over in such a way that they conspire to produce a calculated $\hat{n}_{AVG}$ which is much more quasi-perpendicular than is indicated by the overall magnetic structure.

These latter two explanations both unfortunately suggest significant breakdowns within our “Quasi-parallel” data for those parameters dependent on $\hat{n}_{AVG}$ ($M_{fms}$, $M_{crit}$, $M_A$, $\theta_{BN}$, and $J^{RH}$) and potentially also within our quasi-perpendicular data set. While we will continue to report our results as measured, it will be important to keep these caveats in mind.

Looking at more than just our median lines, another inconsistency becomes apparent at high $J^{Fit}$. In our calculation of $J^{RH}$ and other Rankine-Hugoniot predicted parameters, we make the simplifying assumption of a single fluid with specific heat in $d=3$ spatial dimensions $\gamma = (d+2)/d = 5/3$. The jump conditions predict that in this case that when $\theta_{BN} \to 90^\circ$, $\lim_{M_{fms} \to \infty} J^{RH} = \lim_{M_A \to \infty} J^{RH} = 4$. Our individual median lines for $M_{fms}$, $\theta_{BN}$, and $\beta$ do asymptote to 4. The median line for $M_A$ fails to asymptotically increase to 4, but just barely. Points in a median line merely represent the value of the 50% quantile in a given bin, however. Our individual shocks, let alone the $M_A$ median line, quite clearly do not possess this asymptotic maximum of $J^{Fit} \le 4$.

The failure of the asymptotic limit to the single fluid $\gamma = 5/3$ prediction for our overall and $M_A$ median data, despite the agreement for the median lines of our input parameters ($M_{fms}$, $\beta$, and $\theta_{BN}$), hints at a more significant systematic disagreement with RH. We calculated $\beta = \beta_i + \beta_e = \frac{n_i T_p}{B^2/2\mu_0}$ and $C_S^2 = (\gamma_i T_p + \gamma_e T_e)/m_p = (5/3 T_p + T_e)/m_p$. The Rankine-Hugoniot Equations and resulting derivations however assume a single fluid that does not distinguish between ions and electrons. Thus, there are only a single temperature $T$, number density $n$, and specific heat $\gamma$ such that $\beta = \frac{n_T}{B^2/2\mu_0}$ and $C_S^2 = \gamma T/m_p = \gamma \beta V_A^2/2$. Such a system can be modeled as a limiting case of two-fluid MHD where $\gamma_e \to \gamma_i$. Empirically, we utilize separate species temperatures $T_{proton}$ and $T_e$ and species specific polytropic indices $\gamma_i = 5/3$ and $\gamma_e = 1$. Thus, quantities such as $J^{RH}$ derived from the RH equations may not be correctly calculated.

In particular, $\lim_{M_{fms} \to \infty} (\lim_{\theta_{BN} \to \infty} J^{RH}) = 4$ is only true for a single fluid system where $\gamma = 5/3$ or a multiple fluid system where ions and electrons have equal densities $n = n_i = n_e$, temperatures $T = T_p + T_e$ and identical specific polytropic indices $\gamma = 5/3$. In either of these cases, $C_S^2 = \gamma T/m_p = \gamma \beta V_A^2/2$ holds. For a 2-fluid system where $\gamma_i \neq \gamma_e$, this simplified relationship between $C_S^2$, $V_A^2$ and $\beta$ used in the derivation of the RH predictions is no longer valid, and thus the predictions, including the asymptotic limit, may no longer be entirely accurate.

### 3.3. Seasonal Variation of the Magnetic Shock Jump

Having covered the discrepancies from the RH predictions, we can now return to discuss how the asymptotic fit magnetic jump $J^{Fit}$ varies with season. Figure 4's right hand side splits the $L_S$ dependence of the jump into seasonal bins of our upstream parameters. Through each are drawn separate seasonal median lines through bins of width $\Delta M_{fms} = 1$, $\Delta(\log_{10}\beta) = 0.2$, $\Delta\theta_{BN} = 5^\circ$, and $\Delta M_A = 1$. Close inspection reveals a very slight seasonal variation throughout: In any given bin other than the large $M_A$ outlying points, the largest median jump is in either winter or autumn. Just as consistently, the lowest median jump is in either summer or spring. The overall symmetry in median jumps within each of the Autumn/Winter or Spring/Summer pairs suggests that correlation exists with distance from the summer solstice.

Quantitatively, we find that the correlation of $J^{Fit}$ with angular distance from $L_S = 270^\circ$ is $r = 0.13$. This suggests that there is minimal direct correlation between $J^{Fit}$ and $\Delta L_S$ from the southern summer solstice. For sake of
image clarity, we do not draw quartile bars for the individual seasons except where the seasonal deviation diverge enough for error bars to be necessary. Where not drawn, the widths of the envelopes allow the reader to imagine them easily. Given that all error bars for each season contain the median lines of each other season, it is therefore likely that this seasonal variation is little more than a statistical artifact. We further note that the same statistically negligible results are obtained when Lyman alpha EUV flux quartiles are used instead of seasons. We conducted further statistical analysis using Mann-Whitney and Kolmogorov-Smirnov tests and found that each season's distribution of $J^{Fit}$ was statistically distinct. However, tests on variation of any given season's distribution from year to year shows significant difference from 1 year to another and in orderings which are not consistent between seasons. This variance can therefore not be attributed to season.

## 3.4. Earth Versus Mars Overshoot

Mellott and Livesey (1987) previously found that for quasi-perpendicular shock overshoots at Earth, A is significantly correlated with critical ratio, electron $\beta (r_{Mellott}^{Q\perp}(A, \beta_e) = .79)$, and Alfven mach number $(r_{Mellott}^{Q\perp}(A, M_A) = .91)$. For purposes of comparing Earth's overshoot amplitude A with Mars', we overlaid their (4 Hz resolution $B_{max}$) plots over our (1 Hz resolution $B_{max}$) counterpart results in the left hand side of Figure 5. Their plots match our own almost exactly. Furthermore, they explain the pairs of high $M_{fms}/M_{crit}$ and high A outliers in Mellott's A versus $M_{fms}/M_{crit}$ plot. When compared with Mars data, they appear entirely in-line with a now readily apparent dependence of A on $f(\log(\beta_e))$. This demonstrates that the overshoot dependence on solar wind parameters is very comparable between Earth and Mars.

As their data consisted entirely of shocks measured to be quasi-perpendicular, we have colored those of our shocks measured to be quasi-parallel separately to differentiate between them. We see no distinction between our quasi-parallel and quasi-perpendicular overshoots, however. Like Mellott and Livesey (1987), we found that the jump is uncorrelated with $\theta_{BN}$ for quasi-perpendicular shocks, and this lack of correlation reoccurs for our data in the quasi-parallel range. This indicates either that overshoots in more highly structured quasi-parallel shocks behave as if they were quasi-perpendicular, and in general the $\theta_{BN}$ independence in the quasi-perpendicular regime can be extended to the entire range, or that our $\theta_{BN}$ are poorly measured in the quasi-parallel regime.

## 3.5. Seasonal and Mach Number Variation in the Overshoot

It would be logical to presume that A should spatially vary in the same sense as (or at least the superpositions of) their most strongly correlated parameters. However, we saw in Section 3.1 that while A and the Mellott parameters are all clearly anticorrelated with shock extent (note $r(\log(\beta_e), \Delta) = -0.376$), the latter are more strongly anticorrelated with altitude and/or SZA than A is (note $r(\log(\beta_e), \text{SZA}) = 0.119$, $r(\log(\beta_e), \text{alt}) = -0.0745$), and less anticorrelated with $\Delta$ than A. A naive hypothesis for the source of this unexpected anticorrelation between A and $\Delta$ is that there may be a transitive anticorrelation with the seasonal shock expansion.

Figure 5 (right) and 6a show that the hypothesis that A will increase with distance from the sun is not supported. While it is true that a minimum exists at approximately the $L_s = 250^\circ$ perihelion, broader trends do not appear to agree. In particular, when considering our individual upstream parameter dependence A plots, we see the opposite trend entirely. In any given bin, summer or spring has the largest median A, while autumn or winter has the lowest. We easily see that the differences between season in each bin are small. Even at large critical ratio (where it is worthwhile to draw quantile bars), we see that our median lines are close enough to differ only by statistical artifacts.

We see from Figure 5 that A is significantly correlated with $M_{ms}$ ($r = 0.634$). A is in fact much more highly correlated with $M_{ms}$ than $M_{fms}/M_{crit}$ ($r = 0.437$), $M_A$ ($r = 0.521$), or $\beta_e$ ($r = 0.465$). That the overshoot amplitude correlates more with $M_{ms}$ than $M_{fms}$ or critical ratio (or any other expected parameter) is peculiar. $M_{ms}$ is a global parameter, the fast Mach number of bulk flow in the direction of the bulk flow. $M_{fms}$ is a local parameter differing from $M_{ms}$ in that it decreases as the inbound plasma arrives at the shock at an oblique angle $\theta_{NV}$, and increases as the angle between the upstream magnetic field and the shock normal decreases.

Returning to the plot of A versus $M_{ms}$, we note that due to the narrowness of the envelope, we can try to normalize the overshoot amplitude to a function of $M_{ms}$ in order to verify if it continues to exhibit any significant variation
<figure><img src="image_7.png" /><figcaption>**Figure 5.** A versus $M_{fms}/M_{crit}$ (top), $\beta_e$ (row 2), $M_A$ (row 3) and $M_{ms}$ (bottom). On the left, we superimpose the 4 Hz resolution data from Mellott and Livesey (1987) (black) over our own $Q_\perp$ (blue) and $Q_\parallel$ (red) 1 Hz overshoot amplitude measurements. On the right we separately plot by southern season: spring (green), summer (red), autumn (gold), and winter (blue). All lines shown are median lines (for binsizes $b_{M_{fms}/M_{crit}} = 1, b_{\beta_e} = .25, b_{\log_{10}\beta_e} = .25, b_{M_A} = 1, b_{M_{ms}} = 1$), with error bars showing the 0.25 and 0.75 quartile values in each. The electronic scanned figures in Mellott are distorted, and thus we attempted to un-distort and properly size to line up with our own axes as accurately as possible.</figcaption></figure>

due specifically to $L_s$. We accomplish this by boxcar averaging our A data over an interval of $\Delta M_{ms} = 0.5$, which we fit to the function $A_{fit} = m_0e^{m_1(M_{ms}-m_2)} + m_3$ using a $\chi^2$ gradient descent algorithm, and then dividing A by the result $A_{fit} = 0.13(M_{ms} - 0.50)^{1.30} - 0.28$. The $M_{ms}$ normalized overshoot amplitude is shown in Figure 6b, for a subset that excludes large $A/A_{fit} > 10$ outliers arising from very small $A_{fit}$. We see that any dependence on $L_s$ has been flattened, such that one can draw a straight line through our error bars. It therefore follows that any variation with $L_s$ is entirely insignificant. The apparent $L_s$ dependence is therefore likely due to some covariance between $M_{ms}$ and $L_s$ and may actually represent a hidden Mach number dependence. Interestingly, we find -as shown in panel (c)—that A decreases with Lyman EUV flux, but this relationship entirely vanishes once A has been normalized (panel (d)).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) *A* as a function of *L*<sub>s</sub> (b) *A* normalized to an empirical function of *M*<sub>ms</sub> within reduced range [0,4]. (c) *A* versus Lyman α Irradiance (d) normalized *A* versus Lyman α Irradiance. Error bars here are 25th and 75th quartile marks. Horizontal purple bands show the range of values which are within all error bars.</figcaption></figure>

**3.6. Solar Cycle Variation**

The orbit of Mars around the Sun—and the seasons thus produced—is not the only periodic phenomenon in the Mars-Sun system. In the above discussion, we neglected to consider the possibility of temporal variation across the 11-year-period solar cycle (Hathaway, 2015). Solar Cycle 24 reached solar maximum in April 2014 and ended at solar minimum in December 2019. The temporal range of our study is just within this interval. Plotting various parameters as functions of time should allow us to investigate possible solar cycle dependencies: a parameter that depends on the solar cycle should have an envelope that monotonically increases or decreases with time, and approach a slope of zero toward either edge. In fact, we found the unexpected result that *A* appears to hold such a solar cycle dependence. We show a plot of such in Figure 7(top).

We can see that *A* versus *time* appears to take the form of the second half of a cosine-like function, with period of approximately 11 years. Solar cycle temporal variation is primarily measured in terms of sunspot number by empirically derived periodic functions which are not (Hathaway, 2015), but can be loosely approximated as, cosine functions. This is almost exactly what one would expect from a solar cycle dependent parameter. Drawing median lines with 180 days bins makes this clearer. A number of empirical functions have been derived to model the asymmetrical variation across a solar period. For sake of simplicity—and opposing maximum and minimum, a (half period) cosine function can be sufficient to lowest order. To confirm our hypothesis of solar cycle
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Overshoot time-series data. Top: A versus time with median lines and quartile bars over 180 days bins, and colored by season. Note the increase in slope toward the middle before leveling off toward 2019. Highly resembles half of a 10–11 year periodic function. middle: A versus time with best fit to a cosine function. Previous maximum and minimum within a few months of Solar Cycle 24’s minimum and maximum respectively. bottom: Time-series data for A normalized by a function of $M_{ms}$ previously shown in Figure 6. Here, the overall temporal trend vanishes entirely.</figcaption></figure>
dependence, we fit our A data to the function $y = -a_0 \cos(\frac{2\pi}{a_1}(t - a_2)) + a_3$ using a $\chi^2$ gradient descent algorithm.
The resulting function is overlaid on top of our data in the middle panel of Figure 7. As can be seen, our hypothesis of sinusoidal (or at least periodic) solar-cycle dependence seems validated.

The large envelope of A at any given time in comparison to the solar cycle variation suggests that this variation is controlled by another parameter entirely. Given the strong correlation of A with the global solar wind parameter $M_{ms}$, this would be a reasonable first choice to verify. In the bottom panel of Figure 7, we plot the time series data of A normalized by the previously discussed function of $M_{ms}$, with median line and quartile error bars. As can be clearly seen, the 11 years sinusoidal variation mostly vanishes. Thus, the empirical solar cycle variation of the overshoot too is due to a hidden $M_{ms}$ dependence.

## 4. Conclusions

It has been well understood that the Martian bow shock expands and contracts with season, reaching its largest distance during the southern summer. One might expect that other features of the bow shock may also vary with season. In order to investigate this we collected data for thousands of MAVEN shock crossings in order to conduct statistical analysis. The structure of a supercritical shock includes two magnetic substructures of interest which arise from different physical scales and processes: The primary substructure is the magnetic jump $J = B_D / B_U$ itself which arises from macroscopic fluid scale physics. The second substructure is the ion gyroradius thick overshoot on the downstream side of the shock discontinuity, which arises from kinetic scale reflective processes. Both of these structures were measured.

We found that $J$ generally agrees with the Rankine-Hugoniot jump conditions in quasi-perpendicular shocks to lowest order for a $\gamma = 5/3$ plasma, but breaks down in the quasi-parallel regime. We also found a large number of high Mach number shocks where the magnetic jump exceeds the theoretical asymptotic limit $J_{\gamma=5/3}^{RH} \le 4$. This is a potential sign that the single fluid approximation from which the Rankine-Hugoniot jump conditions are derived is insufficient to model the shock at Mars and that a 2-fluid equivalent may be necessary.

In terms of the overshoot substructure, we found that the dependence of the overshoot amplitude at Mars on upstream parameters $M_{fms}/M_{crit}$, $M_A$ and $\beta_e$, when $B_{max}$ is measured at 1 s resolution, is in very good agreement with corresponding trends at Earth found in prior literature in which $B_{max}$ is measured at 1/4 s resolution. These results are in sufficiently high agreement that the outlying points in Mellott and Livesey (1987), for which many less shocks were measured, are firmly within the envelope of our own results. The most highly correlating upstream parameters, as determined by Mellott and Livesey (1987), were not the most highly correlating parameters in our own data set. We found that the overshoot amplitude is most significantly correlated with the solar wind flow Mach number $M_{ms}$, which is the primary controller of the shock extent. We hypothesized that the increased number of slower upstream pickup ions during summer dust storm season may have some effect on the size of the overshoot, which would be reflected in a variation of the overshoot amplitude with season. We found that no such significant trend with season exists. We did however find a very small and on average almost monotonic temporal variation arising from $M_{ms}$ which may be due to the solar cycle. Together, these results suggest that the Martian overshoot is much more significantly controlled by properties of the solar wind than any variable feature of the planet itself.

The coverage of our study was of insufficient duration to include solar maximum or minimum. A larger study may identify additional solar cycle trends. We also neglected to properly measure the size of the foot except insofar as attempting to locate times sufficiently upstream of the shock to also be upstream of the foot. Shock feet are generated by the same process that generates the overshoot. A counterpart study focused on the foot may be capable of better identifying variations induced by increased numbers of pickup ions undergoing shock drift acceleration. Given the similarity between our Mars overshoot results and prior Earth literature, it may be informative to determine whether a single empirical trend exists at all bow shocks regardless of how the bow shock originated, or whether this paper's result is simply a coincidence. This may be achieved in a follow up study focused on other planets, such as Venus.

## Data Availability Statement

The aggregated data used for our results in the study, as well as the codebase developed for this study, are available at Zenodo via https://doi.org/10.5281/zenodo.8161685 with the MIT License (Fruchtman, 2023a). Due to the file size of the code base, https://doi.org/10.5281/zenodo.8161688 contains just the aggregated data, with the Creative Commons Attribution 4.0 International license (Fruchtman, 2023b).
Acknowledgments

The MAVEN project is supported by NASA through the Mars Exploration Program, and MAVEN data are publicly available through the Planetary Data System. This work was funded by Grant 80NSSC20K0571 from the NASA Solar System Workings program.

References

Acuña, M., Connerney, J., Wasilewski, P., Lin, R., Anderson, K., Carlson, C., et al. (1992). Mars observer magnetic fields investigation. *Journal of Geophysical Research*, **97**(E5), 7799–7814. https://doi.org/10.1029/92je00344

Acuña, M., Connerney, J., Wasilewski, P. a., Lin, R., Anderson, K., Carlson, C., et al. (1998). Magnetic field and plasma observations at mars: Initial results of the mars global surveyor mission. *Science*, **279**(5357), 1676–1680. https://doi.org/10.1126/science.279.5357.1676

Angelopoulos, V., Cruce, P., Drozdov, A., Grimes, E. W., Hatzigeorgiu, N., King, D. A., et al. (2019). The space physics environment data analysis system (SPEDAS) [Software]. Space Science Reviews, **215**, 9. https://doi.org/10.1007/s11214-018-0576-4

Bale, S., Balikhin, M., Horbury, T., Krasnoselskikh, V., Kucharek, H., Möbius, E., et al. (2005). Quasi-perpendicular shock structure and processes. *Space Science Reviews*, **118**(1), 161–203. https://doi.org/10.1007/s11214-005-3827-0

Bertucci, C., Duru, F., Edberg, N., Fraenz, M., Martínez, C., Szego, K., & Vaisberg, O. (2011). The induced magnetospheres of mars, venus, and titan. *Space Science Reviews*, **162**(1), 113–171. https://doi.org/10.1007/s11214-011-9845-1

Bertucci, C., Mazelle, C., Acuna, M., Russell, C., & Slavin, J. (2005). Structure of the magnetic pileup boundary at mars and Venus. *Journal of Geophysical Research*, **110**(A1), A01209. https://doi.org/10.1029/2004ja010592

Bhattacharyya, D., Clarke, J. T., Bertaux, J.-L., Chaufray, J.-Y., & Mayyasi, M. (2015). A strong seasonal dependence in the Martian hydrogen exosphere. *Geophysical Research Letters*, **42**(20), 8678–8685. https://doi.org/10.1002/2015gl065804

Burgess, D., Lucek, E., Scholer, M., Bale, S., Balikhin, M., Balogh, A., et al. (2005). Quasi-parallel shock structure and processes. *Space Science Reviews*, **118**(1), 205–222. https://doi.org/10.1007/s11214-005-3832-3

Burne, S., Bertucci, C., Mazelle, C., Morales, L. F., Meziane, K., Halekas, J., et al. (2021). The structure of the Martian Quasi-Perpendicular Super-critical shock as seen by MAVEN. *Journal of Geophysical Research (Space Physics)*, **126**(9), e28938. https://doi.org/10.1029/2020JA028938

Chaffin, M., Deighan, J., Schneider, N., & Stewart, A. (2017). Elevated atmospheric escape of atomic hydrogen from mars induced by high-altitude water. *Nature Geoscience*, **10**(3), 174–178. https://doi.org/10.1038/ngeo2887

Chaffin, M. S., Chaufray, J.-Y., Stewart, I., Montmessin, F., Schneider, N. M., & Bertaux, J.-L. (2014). Unexpected variability of Martian hydrogen escape. *Geophysical Research Letters*, **41**(2), 314–320. https://doi.org/10.1002/2013GL058578

Connerney, J., Espley, J., Lawton, P., Murphy, S., Odom, J., Oliversen, R., & Sheppard, D. (2015). The maven magnetic field investigation. *Space Science Reviews*, **195**(1), 257–291. https://doi.org/10.1007/s11214-015-0169-4

Curvefit. (2013). Curvefit [Software]. L3Harris Geospatial Solutions, Inc. Retrieved from https://www.l3harrisgeospatial.com/docs/curvefit.html

Dubinin, E., Fraenz, M., Woch, J., Barabash, S., Lundin, R., & Yamauchi, M. (2006). Hydrogen exosphere at mars: Pickup protons and their acceleration at the bow shock. *Geophysical Research Letters*, **33**(22), L22103. https://doi.org/10.1029/2006GL027799

Edberg, N., Lester, M., Cowley, S., & Eriksson, A. (2008). Statistical analysis of the location of the Martian magnetic pileup boundary and bow shock and the influence of crustal magnetic fields. *Journal of Geophysical Research*, **113**(A8), A08206. https://doi.org/10.1029/2008ja013096

Farris, M. H., & Russell, C. T. (1994). Determining the standoff distance of the bow shock: Mach number dependence and use of models. *Journal of Geophysical Research*, **99**(A9), 17681–17689. https://doi.org/10.1029/94JA01020

Fruchtman, J. (2023a). jacobfruchtman/MAVEN-shock-analyzer: MAVEN Shock analyzer (with expanded dataset) [Software]. Zenodo. https://doi.org/10.5281/zenodo.8161685

Fruchtman, J. (2023b). Maven shock crossing data [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.8161688

Garnier, P., Jacquey, C., Gendre, X., Génot, V., Mazelle, C., Fang, X., et al. (2022a). The drivers of the Martian bow shock location: A statistical analysis of mars atmosphere and volatile evolution and mars express observations. *Journal of Geophysical Research: Space Physics*, **127**(5), e2021JA030147. https://doi.org/10.1029/2021JA030147

Garnier, P., Jacquey, C., Gendre, X., Génot, V., Mazelle, C., Fang, X., et al. (2022b). The influence of crustal magnetic fields on the Martian bow shock location: A statistical analysis of maven and mars express observations. *Journal of Geophysical Research: Space Physics*, **127**(5), e2021JA030146. https://doi.org/10.1029/2021JA030146

Gedalin, M. (1996). Ion reflection at the shock front revisited. *Journal of Geophysical Research*, **101**(A3), 4871–4878. https://doi.org/10.1029/95ja03669

Gosling, J., Thomsen, M., Bame, S., & Russell, C. (1989). Ion reflection and downstream thermalization at the quasi-parallel bow shock. *Journal of Geophysical Research*, **94**(A8), 10027–10037. https://doi.org/10.1029/ja094ia08p10027

Gruesbeck, J. R., Espley, J. R., Connerney, J. E. P., DiBraccio, G. A., Soobiah, Y. I., Brain, D., et al. (2018). The three-dimensional bow shock of mars as observed by maven. *Journal of Geophysical Research: Space Physics*, **123**(6), 4542–4555. https://doi.org/10.1029/2018JA025366

Halekas, J. (2017). Seasonal variability of the hydrogen exosphere of mars. *Journal of Geophysical Research: Planets*, **122**(5), 901–911. https://doi.org/10.1002/2017JE005306

Halekas, J., Ruhunusiri, S., Harada, Y., Collinson, G., Mitchell, D., Mazelle, C., et al. (2017). Structure, dynamics, and seasonal variability of the mars-solar wind interaction: Maven solar wind ion analyzer in-flight performance and science results. *Journal of Geophysical Research: Space Physics*, **122**(1), 547–578. https://doi.org/10.1002/2016ja023167

Halekas, J., Taylor, E., Dalton, G., Johnson, G., Curtis, D., McFadden, J., et al. (2015). The solar wind ion analyzer for maven. *Space Science Reviews*, **195**(1), 125–151. https://doi.org/10.1007/s11214-013-0029-z

Hall, B. E. S., Lester, M., Sánchez-Cano, B., Nichols, J. D., Andrews, D. J., Edberg, N. J. T., et al. (2016). Annual variations in the Martian bow shock location as observed by the mars express mission. *Journal of Geophysical Research: Space Physics*, **121**(11), 11474–11494. https://doi.org/10.1002/2016JA023316

Hao, Y., Lu, Q., Gao, X., & Wang, S. (2016). Ion dynamics at a rippled quasi-parallel shock: 2d hybrid simulations. *The Astrophysical Journal*, **823**(1), 7. https://doi.org/10.3847/0004-637x/823/1/7

Hathaway, D. H. (2015). The solar cycle. *Living Reviews in Solar Physics*, **12**(1), 4. https://doi.org/10.1007/lrp2-2015-4

Heavens, N. G., Kleinböhl, A., Chaffin, M. S., Halekas, J. S., Kass, D. M., Hayne, P. O., et al. (2018). Hydrogen escape from mars enhanced by deep convection in dust storms. *Nature Astronomy*, **2**(2), 126–132. https://doi.org/10.1038/s41550-017-0353-4

Jakosky, B. M., Lin, R. P., Grebowsky, J. M., Luhmann, J. G., Mitchell, D., Beutelschies, G., et al. (2015). The mars atmosphere and volatile evolution (maven) mission. *Space Science Reviews*, **195**(1), 3–48. https://doi.org/10.1007/s11214-015-0139-x

Kennel, C. F., Edmiston, J. P., & Hada, T. (1985). A quarter century of collisionless shock research (Vol. 34, pp. 1–36). Washington DC American Geophysical Union Geophysical Monograph Series. https://doi.org/10.1029/GM034p0001

Ma, Y., Nagy, A. F., Sokolov, I. V., & Hansen, K. C. (2004). Three-dimensional, multispecies, high spatial resolution MHD studies of the solar wind interaction with mars. *Journal of Geophysical Research*, **109**(A7), A07211. https://doi.org/10.1029/2003JA010367

Marshall, W. (1955). The structure of Magneto-Hydrodynamic shock waves. *Proceedings of the Royal Society of London Series A*, **233**(1194), 367–376. https://doi.org/10.1098/rspa.1955.0272
Mazelle, C., Winterhalter, D., Sauer, K., Trotignon, J., Acuña, M., Baumgärtel, K., et al. (2004). Bow shock and upstream phenomena at Mars. *Space Science Reviews*, 111(1), 115-181. https://doi.org/10.1023/b:spac.0000032717.98679.d0

Mellott, M., & Livesey, W. (1987). Shock overshoots revisited. *Journal of Geophysical Research*, 92(A12), 13661-13665. https://doi.org/10.1029/ja092ia12p13661

Mitchell, D., Mazelle, C., Sauvaud, J.-A., Thocaven, J.-J., Rouzaud, J., Fedorov, A., et al. (2016). The maven solar wind electron analyzer. *Space Science Reviews*, 200(1), 495-528. https://doi.org/10.1007/s11214-015-0232-1

Peredo, M., Slavin, J. A., Mazur, E., & Curtis, S. A. (1995). Three-dimensional position and shape of the bow shock and their variation with alfvénic, sonic and magnetosonic mach numbers and interplanetary magnetic field orientation. *Journal of Geophysical Research*, 100(A5), 7907-7916. https://doi.org/10.1029/94JA02545

Rahmati, A., Larson, D., Cravens, T., Lillis, R., Halekas, J., McFadden, J., et al. (2018). Seasonal variability of neutral escape from mars as derived from maven pickup ion observations. *Journal of Geophysical Research: Planets*, 123(5), 1192-1202. https://doi.org/10.1029/2018je005560

Russell, C., Hoppe, M., & Livesey, W. (1982). Overshoots in planetary bow shocks. *Nature*, 296(5852), 45-48. https://doi.org/10.1038/296045a0

Russell, C. T., Luhmann, J. G., & Strangeway, R. J. (2016). *Space physics: An introduction*. Cambridge University Press.

Schwartz, S. J. (1991). Magnetic field structures and related phenomena at quasi-parallel shocks. *Advances in Space Research*, 11(9), 231-240. https://doi.org/10.1016/0273-1177(91)90039-m

Schwartz, S. J. (1998). Shock and discontinuity normals, mach numbers, and related parameters. *ISSI Scientific Reports Series*, 1, 249-270.

Schwartz, S. J., & Burgess, D. (1991). Quasi-parallel shocks: A patchwork of three-dimensional structures. *Geophysical Research Letters*, 18(3), 373-376. https://doi.org/10.1029/91GL00138

Szegö, K., Glassmeier, K.-H., Bingham, R., Bogdanov, A., Fischer, C., Haerendel, G., & others (2000). Physics of mass loaded plasmas. *Space Science Reviews*, 94(3), 429-671. https://doi.org/10.1023/a:1026568530975

Trotignon, J., Mazelle, C., Bertucci, C., & Acuña, M. (2006). Martian shock and magnetic pile-up boundary positions and shapes determined from the phobos 2 and mars global surveyor data sets. *Planetary and Space Science*, 54(4), 357-369. https://doi.org/10.1016/j.pss.2006.01.003

Vignes, D., Mazelle, C., Rme, H., Acuña, M. H., Connerney, J. E. P., Lin, R. P., et al. (2000). The solar wind interaction with mars: Locations and shapes of the bow shock and the magnetic pile-up boundary from the observations of the MAG/ER experiment onboard mars global surveyor. *Geophysical Research Letters*, 27(1), 49-52. https://doi.org/10.1029/1999GL010703

Williams, D. R. (2021). *Mars fact sheet*—NASA. NASA. Retrieved from https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html

Wilson, L., III., Koval, A., Szabo, A., Stevens, M., Kasper, J., Cattell, C., & Krasnoselskikh, V. (2017). Revisiting the structure of low-Mach number, low-beta, quasi-perpendicular shocks. *Journal of Geophysical Research: Space Physics*, 122(9), 9115-9133. https://doi.org/10.1002/2017ja024352

Wong, H. (1968). Foot of normal hydromagnetic shock. *Bulletin of the American Physical Society*, 13, 1518.

Woods, L. (1971). On double-structured, perpendicular, magneto-plasma shock waves. *Plasma Physics*, 13(4), 289-302. https://doi.org/10.1088/0032-1028/13/4/302

Zhang, T., Delva, M., Baumjohann, W., Volwerk, M., Russell, C., Wei, H., et al. (2008). Induced magnetosphere and its outer boundary at Venus. *Journal of Geophysical Research*, 113(E9), E00B20. https://doi.org/10.1029/2008je003215

Zhang, T.-L., Schwingenschuh, K., Russell, C. T., & Luhmann, J. G. (1991). Asymmetries in the location of the Venus and mars bow shock. *Geophysical Research Letters*, 18(2), 127-129. https://doi.org/10.1029/90GL02723

Zurek, R. W., & Martin, L. J. (1993). Interannual variability of planet-encircling dust storms on Mars. *Journal of Geophysical Research*, 98(E2), 3247-3259. https://doi.org/10.1029/92je02936

Erratum

The originally published version of this article contained an incorrect version of Figure 3. The figure has since been replaced, and this may be considered the authoritative version of record.