A Fast Bow Shock Location Predictor-Estimator From 2D and 3D Analytical Models: Application to Mars and the MAVEN Mission.

We present fast algorithms to automatically estimate the statistical position of the bow shock from spacecraft data, using existing analytical two-dimensional (2D) and three-dimensional (3D) models of the shock surface. We derive expressions of the standoff distances in 2D and 3D and of the normal to the bow shock at any given point on it. Two simple bow shock detection algorithms are constructed, one solely based on a geometrical predictor from existing models, the other using this predicted position to further refine it with the help of magnetometer data, an instrument flown on many planetary missions. Both empirical techniques are applicable to any planetary environment with a defined shock structure. Applied to the Martian environment and the NASA/MAVEN mission, the predicted shock position is on average within 0.15 planetary radius R p of the bow shock crossing. Using the predictor-corrector algorithm, this estimate is further refined to within a few minutes of the true crossing (≈0.05R p). Between 2014 and 2021, we detect 14,929 clear bow shock crossings, predominantly quasi-perpendicular. Thanks to 2D conic and 3D quadratic fits, we investigate the variability of the shock surface with respect to Mars Years (MY), solar longitude (Ls), and solar EUV flux levels. Although asymmetry in Y and Z Mars Solar Orbital coordinates is on average small, we show that for MY32 and MY35, Ls = [135°-225°] and high solar flux, it can become particularly noticeable, and is superimposed to the usual North-South asymmetry due in part to the presence of crustal magnetic fields. ©2022. The Authors.

Introduction

Historically, planetary bow shocks, their position, size, and shape, have been characterized statistically with the use of (empirical) analytical fitting models in two‐dimensional (2D) or three‐dimensional (3D) spatial coordinates. A classical starting point for characterizing the Earth’s bow shock in 3D includes the seminal work of Formisano (1979), who investigated the asymmetry of the shock with respect to the apparent solar wind flow direction, with the use of quadratic surface fits with nine free parameters. In parallel, other studies, such as that of Slavin and Holzer (1981) relied on a simple polar equation assuming axisymmetry along the Sun‐planet line, corrected by the apparent motion of the solar wind in the rest frame of the planet, the so‐called aberrated X axis. The 2D approach has the merit of needing only three free parameters but ignores the potential asymmetries of the shock as, for example, seen at Earth’s bow shock (e.g., Formisano, 1979; Merka et al., 2005; Peredo et al., 1995, 1993).

More advanced physics‐based models have also been proposed as a complement to those empirical attempts. A good introduction into analytical models of the bow shock, based on gasdynamic theory and magnetohydrodynamics solutions, is given in Verigin, Slavin, Szabo, Gombosi, et al. (2003) and recently in Kotova et al. (2021). These studies present analytical functions describing the curvature, bluntness, and skewing angle of the shock structure, which are arguably better suited to the fitting of the shock flanks; they are applicable to many planetary bow shock conditions.

At Mars, due to the sparsity of early data and the non‐collisional nature of the shock, the tendency has been to use in priority the simplest fitting model available with least free parameters, that is, an empirical 2D polar equation (Edberg et al., 2008; Hall et al., 2019; Russell, 1977; Slavin et al., 1991; Trotignon et al., 2006). Only recently with the NASA/Mars Atmospheric and Volatile EvolutioN (MAVEN) mission were quadratic fits used to characterize the general structure of the Martian bow shock, with Gruesbeck et al. (2018) providing fits to a careful subset of identified crossings in the first year of operations of the MAVEN mission.

In recent years, many studies have attempted to characterize the Martian shock position and shape and its evolution under various solar wind and EUV conditions (Hall et al., 2019 and references therein). Two missions have been used for this goal, the ESA/Mars Express mission and the NASA/MAVEN mission. Mars Express (hereafter MEX for brevity) was launched in 2003 and has been orbiting Mars since 2004, whereas MAVEN was launched 10 yr later in November 2013 and has been orbiting the planet since 22 September 2014. MAVEN’s scientific payload includes among others a fluxgate magnetometer (MAG, Connerney et al., 2015), two ion spectrometers including the Solar Wind Ion Analyzer (SWIA, Halekas et al., 2015) and the Suprathermal and Thermal Ion Composition instrument (STATIC, McFadden et al., 2015), and an electron spectrometer (Solar Wind Electron Analyzer, SWEA, Mitchell et al., 2016). MEX unfortunately does not carry any magnetometer but includes a plasma suite (ions and electrons) as part of the ASPERA‐3 package (Barabash et al., 2006), which was used to investigate the plasma boundaries at Mars (Dubinin et al., 2006). Both missions aim at studying the upper atmosphere and the magnetospheric environment of Mars.

Detections of the bow shock in spacecraft data have relied on manual determinations using as many instruments (including plasma instruments and magnetometer) as available to avoid ambiguous detections (Gruesbeck et al., 2018). Recently, Němec et al. (2020) proposed a region identification scheme based on selected plasma parameters and applied their technique to the MAVEN dataset in order to identify upstream solar wind, magnetosheath, and magnetosphere regions crossed by the spacecraft. This method has the advantage of mapping these regions statistically, removing certain biases usually associated with manually picked individual boundary crossings which may be orbit‐dependent. However, it requires the reliable knowledge of flow speed, ion density and magnetic field, which may not all always be available. A parallel trend has also been to apply machine‐learning techniques to plasma data for the labeling and identification of the regions crossed by a given spacecraft (see Breuillard et al., 2020; Hall et al., 2016; for Mars and Earth, respectively), as part of online databases (Génot et al., 2021). However, these studies require extensive amounts of time and patience, as well as large computer resources to become efficient. Sometimes, a precise determination of the shock position in the data is not of paramount importance and simpler, faster, more straightforward approaches, such as the one presented in this study, can be advantageous. This may be the case in statistical studies where one of those regions needs to be systematically excluded, as in space weather databases monitoring the solar wind. This can also be of interest when areas around the predicted bow shock must be excluded for qualified reasons, as in wave studies focusing on regions outside of foreshock and shock wake structures, or when a first guess of the location and geometry of the shock is needed.

We present in this study new, simple analytical algorithms using two types of historical fitting techniques (2D and 3D), which make it possible to quickly estimate from spacecraft spatial coordinates the statistical geometrical position of the shock in planetary atmospheres. Special emphasis on the Martian environment and the MAVEN dataset is given throughout, however the method is applicable to any planetary environment and spacecraft dataset. This first crude estimator can be refined further by applying additional criteria, for example, on the magnetic field amplitude measured by the MAVEN spacecraft. This provides a fast means to approximately and quite reliably identify the position of the shock so that solar wind and magnetosheath regions can be studied on a statistical level in the data. Moreover, other characteristics of the shock crossing, such as the quasi‐parallel (q ‖) or quasi‐perpendicular (q ⊥) geometry of the shock can easily be obtained by deriving the perpendicular direction to the shock at any point on the surface.

After a review of 2D and 3D bow shock fitting models at Mars in Section 2, we present in a consistent manner the leading equations behind these models in 2D and 3D, give the analytical expressions for standoff distances, and propose a geometric calculation of the normal to the bow shock at any point on the surface. Starting in Section 3, we introduce a predictor algorithm for a fast estimation of the shock position in spacecraft orbital coordinates and its timing (Section 3.1). In Section 3.2, we propose a simple correction on this position and timing with the sole help of magnetometer data (predictor‐corrector algorithm). Application to the MAVEN MAG dataset is then given as validation on a few examples and then extended to the whole available dataset. Finally, as a result of the automatic detection proposed here, statistical analytical fits are given for the MAVEN mission between November 2014 and February 2021, with a discussion of the shock’s asymmetry based on terminator and standoff distances. Applications for space weather‐related databases are also mentioned.

Bow Shock Models at Mars

In this section, following a survey of existing fitting models at Mars, we present comprehensive formulae for analytical fits in 2D polar coordinates and 3D Cartesian coordinates, with a calculation of bow shock subsolar and terminator standoff distances. We also show how to calculate the normal to the surface at a given point in space, in order to estimate the q ⊥ − q ‖ shock conditions.

Coordinate Systems and Solar Wind Flow Aberration

All spacecraft coordinates in this study are in Mars Solar Orbital coordinates (MSO) for simplicity, in accordance with most previous studies. In the MSO system, identical to the Sun‐state coordinate system, the +X MSO axis points toward the Sun from the planet’s center, +Z MSO is toward Mars’ North pole and perpendicular to the orbital plane defined as the X MSO–Y MSO plane passing through the center of Mars, and Y MSO completes the orthogonal system.

Because of the orbital motion of Mars with respect to the average direction of the solar wind, the apparent direction of the solar wind in the rest frame of the planet deviates from the anti‐sunward direction. As a result, an anti‐clockwise rotation by an angle α around the Z axis must be applied so that the bow shock’s major axis is aligned with respect to the X axis along the solar wind flow. This aberration, first seen in cometary tails and at the origin of the hypothesis by Biermann of a stellar wind (Biermann, 1951), is taken into account in the so‐called aberrated MSO coordinates, denoted , , and (although Z is left unchanged by the transformation). To unclutter notations, the “MSO” subscript is now dropped. Following Formisano et al. (1979), Slavin and Holzer (1981) define the angle as where V p is Mars’ orbital velocity and is the solar wind velocity, for example expressed in km s−1. Mars’ average orbital velocity is  km s−1. For the maximum value (26.5 km s−1), the angle is α = 3.8° for a typical solar wind speed of 400 km s−1. The angle assumed by all studies except those of Slavin and Holzer (1981) and Slavin et al. (1991) is 4°. In Slavin and Holzer (1981), the aberration angle was chosen to be varying with solar wind speed conditions. In Slavin et al. (1991), α = 3.2°. Figure 1 shows the aberration angle with respect to orbital velocity (abscissa) and to solar wind velocity (color code).

Figure 1

Aberration angle (in °) with respect to the orbital velocity of Mars (in km s−1) and the solar wind mean speed (color‐coded isocontours, in km s−1).

Parametric Models

Table 1 chronologically lists past bow shock studies at Mars and their characteristics in terms of solar activity, solar cycle number, and number of observations, including the recent MEX and MAVEN surveys.

Table 1

Statistical Studies on the Martian Bow Shock Position Replaced Chronologically (With Respect to In‐Situ Observations) in the Context of Solar Activity and Mars Year (MY)

References	Spacecraft	Years	N	Solar activity	Cycle #	Start	Max.	MY
Slavin and Holzer (1981)	Mariner 4	1965	2	Low	20	1964	1968	6
Russell (1977) a	Mars 2, 3, 5	1965–1974	11	Low‐Medium	20	1964	1968	9–11
Slavin and Holzer (1981)	Mars 2, 3	1971–1972	10	Medium	20	1964	1968	9–10
Slavin and Holzer (1981)	Mars 5	1974	4	Low	20	1964	1968	11
Slavin et al. (1991)	Mariner 4, Mars 2, 3, 5	1965–1974	24	Low‐Medium	20	1964	1968	6–11
Schwingenschuh et al. (1990)	Phobos 2	1989	∼100	High	22	1986	1989	19
Slavin et al. (1991)	Phobos 2	1989	94	High	22	1986	1989	19
Zhang et al. (1991)	Phobos 2	1989	97	High	22	1986	1989	19
Trotignon et al. (1991)	Phobos 2	1989	26	High	22	1986	1989	19
Trotignon et al. (1993) b	Phobos 2	1989	126	High	22	1986	1989	19
Verigin et al. (1993), (1999) c	Phobos 2	1989	126	High	22	1986	1989	19
Vignes et al. (2000)	MGS	1997–1998	450	Low	23	1996	2001	23–24
Trotignon et al. (2006)	MGS	09/1997–02/1999	573	Low‐Medium	23	1996	2001	23–24
Edberg et al. (2008)	MGS	09/1997–02/1999	619	Low‐Medium	23	1996	2001	23–24
Hall et al. (2016),  (2019) d	MEX	2004–2008	4,422	Medium‐Low	23	1996	2001	27–29
Ramstad et al. (2017) e	MEX	11/2005–12/2016	2,166	High‐Medium	23, 24	1996, 2008	2001, 2014	27–33
Hall et al. (2016),  (2019) f	MEX	2008–2015	7,669	Low‐High	24	2008	2014	30–32
Hall et al. (2019) g	MEX	2015–12/2017	1,494	High‐Medium	24	2008	2014	33
Halekas et al. (2017) h	MAVEN	10/2014–05/2016	−	High‐Medium	24	2008	2014	32–33
Gruesbeck et al. (2018)	MAVEN	11/2014–04/2017	1,799	High‐Medium	24	2008	2014	32–34
Němec et al. (2020) i	MAVEN	11/2014–02/2020	−	High‐low	24	2008	2014	32–35

Note. N is the number of bow shock crossings considered in each study. MGS, Mars Global Surveyor. MEX, Mars Express. Column "Start" marks the start year of the corresponding solar cycle, whereas column "Max." gives the year when the maximum of the solar cycle occurred.

Observations by Gringauz et al. (1976) and analyzed further by Russell (1977).

In Trotignon et al. (2006), 127 Phobos 2 crossings of the bow shock were reported, that is, one more than in Trotignon et al. (1993).

Number of observations from Phobos 2 used in this study was presumably the same as in Trotignon et al. (1993).

Data from Hall et al. (2019), MY27–29, Table 3  with respect to Mars Years (MY).

About 7,000 orbits were first manually examined, “out of which 1,083 orbit inbound and outbound segments with identified BS, IMB […] crossings were included.” Orbital coverage of MEX is shown in their Figure 9. No discrimination with solar cycle or MY is given, although EUV flux and solar‐wind parameter dependence are studied.

Data from Hall et al. (2019), MY30–32, Table 3 with respect to MY. MY32 runs from 31 July 2013 to 17 June 2015.

Data from Hall et al. (2019), MY33, Table 3 with respect to MY. MY33 runs from 18 June 2015 to 4 May 2017.

Bow shock variations are obtained by fitting 2D‐gridded datasets of average plasma density jumps through the shock location as measured with MAVEN/SWIA and are discriminated against magnetosonic Mach number M ms, EUV flux, and dynamic pressure.

Number of individual crossings not disclosed due to the nature of the region detection scheme used. A total of 2,040 full orbits reported.

Bow shock models at Mars have been proposed since the end of the 1970s, including (but not limited to): Russell (1977), Slavin and Holzer (1981), Slavin et al. (1991), Schwingenschuh et al. (1990), Trotignon et al. (1991), Zhang, Schwingenschuh, Russell, and Luhmann (1991), Zhang, Schwingenschuh, Lichtenegger, et al. (1991), Trotignon et al. (1993), Verigin et al. (1999), Vignes et al. (2000), Trotignon et al. (2006), and Edberg et al., 2008, 2010. These studies were performed with several spacecraft including Viking, Mars Global Surveyor (MGS) and Phobos‐2, and for varying solar conditions. Gringauz et al. (1976) quoted in Russell (1977) reported 11 crossings for the Russian Mars‐2, 3, 5 satellites, but Slavin and Holzer (1981) later reanalyzed the datasets and found 14 crossings in total. Slavin et al. (1991) reported 94 crossings for Phobos 2, upped to 126 by Trotignon et al. (1993), as later reported in Trotignon et al. (2006). Using the same dataset, Verigin et al. (1993) reported a weak dependence of terminator bow shock position on solar wind dynamic pressure P sw; additionally in an analogy with magnetized planets regarding the dependence of the magnetopause thickness to P sw, these authors anticipated the discovery of crustal magnetic field sources at Mars, later vindicated with MGS measurements (Acuña et al., 1998).

In contrast with what was found at Venus, Slavin and Holzer (1981) and later Vignes et al., 2000,  2002; using MGS and data from previous mission suggested that the mean bow shock standoff distance was likely independent of the solar activity. Slavin et al. (1991) showed that the terminator distance, which is a marker of the swelling of the cavity flanks, varied by as much as 11% between the Mars‐2, 3, 5 observations (low activity) and the Phobos‐2 observations (high activity), although the number of crossings for each mission largely differed. Mars EXpress (MEX), with its very long activity spanning the end of solar cycle 23 and cycle 24 up to now (toward the beginning of new cycle 25), had the best chance to conclusively solve this aspect: Hall et al. (2019) found that for the years 2004–2017, the terminator distance underwent variations up to ∼7%, in agreement with the preliminary results of Trotignon et al. (1993). Of note, Mazelle et al. (2004) made a review of all available measurements before MEX started observing, and discussed the solar cycle variations and the differences observed with Venus (for which up to 35% increases of the bow shock location in the terminator plane with increasing activity have been reported, see Edberg et al., 2010; Russell et al., 1988; Zhang et al., 2008).

Of great import, Edberg et al. (2009, 2010) used MGS and MEX data in combination with ACE data extrapolated to Mars to study the dependence of the bow shock location to solar EUV flux and magnetosonic Mach number (noted M ms). They pointed out that the shape of the magnetosonic shock wave depends on the ratio of the solar wind speed to the magnetosonic speed. Later, for the entire period 2004–2015, Hall et al. (2016) identified 11,861 crossings in the MEX database, using electron spectrometer data with machine‐learning algorithms. This study was extended by Hall et al. (2019) up to 2017 with a total of 13,585 crossings, totaling 13 yr of operations of MEX. Both works used a standard 2D conic fit depending on the Mars Year (MY), with observed variations up to a few percent in terms of standoff bow shock distance. In addition to solar wind cycle variations and hemispherical changes from a 2D perspective, these studies confirmed the dependence of the shock position with P sw and, most drastically, with solar EUV flux. Correspondingly, Ramstad et al. (2017) studied with coincident electron and ion data a subset of only 1,083 inbound and outbound MEX orbits for the period 2005–2016. They evaluated the dependence of the Induced Magnetospheric Boundary (IMB) and bow shock (BS) to EUV flux and the solar wind’s lowest moments (density and bulk velocity), showing that the BS mostly expands and contracts with the IMB. However, they also found that the BS swelling in the flank due to increased EUV fluxes cannot be solely explained by a corresponding swelling of the IMB. Simultaneously with MAVEN with both magnetometer and ion measurements, Halekas et al. (2017) investigated how the Martian magnetosphere and bow shock responded to EUV flux, M ms, and P sw between October 2014 and May 2016 (0.85 MY). In agreement with previous studies, they showed that the shock inflates with increasing EUV flux and contracts with increasing dynamic pressure and M ms; this in turn leads to EUV flux and dynamic pressure competing against one another because of their common dependence on heliocentric distance d .

Recently, Němec et al. (2020) used MAVEN plasma and magnetometer data to construct maps of the solar wind/magnetosheath regions. Their modeled bow shock locations, explicitly dependent on P sw, EUV flux, and crustal field intensities, were in good agreement with the average fits of Trotignon et al. (2006), with appreciable differences in terminator extensions as compared to the results of Ramstad et al. (2017).

Broadly speaking, two approaches fitting the shape of bow shocks have historically been employed, one using a simple 2D polar form (e.g., Slavin & Holzer, 1981), the other the 3D general Cartesian conic form (e.g., Formisano, 1979; Formisano et al., 1979). Semi‐empirical models based on gasdynamic and MHD predictions such as those of Verigin et al., 1993, 1999 are not discussed in the following.

The 2D Polar Form

Assuming that a full‐formed bow shock in aberrated coordinates is symmetric with respect to the Xʹ axis, the 3D shape of the bow shock can be reduced to a 2D problem in the plane. All Martian studies except that of Gruesbeck et al. (2018) have made this assumption. A simple 2D conic of revolution (usually a parabola or a hyperbola), symmetric around the aberrated MSO axis X′ and decentered from its focus x is shown in Figure 2. Such a 2D conic takes the parametric form (e.g., Hall et al., 2019): where ϵ is the conic’s eccentricity, L the semilatus rectum (called terminator crossing by Volwerk et al., 2016; Zhang et al., 2008, at Venus, because the focus is taken at the center of the planet) and θ the angle measured from the focus of the conic, typically within the [−π/2, π/2] range. This range of angles depends on the nature of the conic section: for a parabola (ϵ = 1) [(borders excluded), for an ellipse (ϵ < 1) ], and for a hyperbola (ϵ > 1) . The equivalent rectangular (Cartesian) form of this equation is (Trotignon et al., 2006):

Figure 2

Typical 2D conic bow shock shape in the aberrated coordinate system. For a point M on the shock surface, ρ is the Euclidean distance to the shock from the center of the planet of radius R p, and r is the distance to the shock surface from the focus x of the conic with semilatus rectum L and making an angle θ with the X′ direction, so that Equation (1) holds. ϑ is the usual polar angle, with respect to the center of Mars. R ss and R td are the standoff subsolar and terminator distances.

In this representation, Trotignon et al. (2006) derived two additional useful quantities, the standoff shock distance along the X axis, R ss (also called subsolar aerocentric distance in Trotignon et al., 2006) and the standoff terminator distance R td along the cylindrical coordinate (which is none other than the diameter of the tail at Xʹ = 0 divided by 2, or as in Russell, 1977, the “dawn radius”):

We can derive another parameter of interest, that is, the aperture of the Mach cone related to the shock structure—the limiting Mach cone angle (Verigin, Slavin, Szabo, Kotova, & Gombosi, 2003). In the gasdynamic approach, it is defined as , where is the sonic Mach number. The sonic speed is , with ρ the solar wind ion mass density, γ = 5/3 the ratio of specific heats, and P = n sw k (T  + T ) the solar wind thermal pressure, with T and T the electron and ion temperatures, respectively. For a hyperbola (ϵ   1), the limiting Mach cone angle is exactly the angle made by the asymptotes of the hyperbola (Slavin et al., 1984):

as uncertainty. In a canonical form for the hyperbola, with a the distance from the nose to the intersection of the asymptotes on the X′ axis, and b that from the shock nose to the asymptote on the Y′ axis, . Since, by definition , the substitution readily yields expression (7). It is noteworthy to remark that for ϵ close to 1, the uncertainty increases to infinity; any determination of ϱ is thus unreliable for quasi‐parabolic curves.

For the fit, Slavin and Holzer (1981) and Slavin et al. (1991) rewrote Equation (1) as y = ax + b (posing y = 1/r, x = cos θ, a = ϵ/L, and b = 1/L) and performed simple linear regressions for a range of focus locations. As pointed out by Vignes et al. (2000), this may result in fitting biases when observations are widely disparate in their location: in this case direct fitting methods to Equation (1) should be preferred. With a direct polar fit to MGS data, Edberg et al. (2008) gave for example, the following fitted values: ϵ = 1.05 ± 0.04, L = 2.10 ± 0.09, and x  = 0.55 ± 0.12. However, to match the results plotted in Figure 1 of Edberg et al. (2008), the value of ϵ must be modified down to ϵ = 1.03, a marginal difference likely due to rounding errors. For comparison, the corresponding values derived by Hall et al. (2019), for MEX data but with a larger sample, are ϵ = 0.998 ± 0.003, L = 1.802 ± 0.002, and x  = 0.76. For MGS and MEX data (Edberg et al., 2008; Hall et al., 2016, 2019), the subsolar standoff distance is R ss = 1.63 ± 0.04 R p, whereas the terminator standoff distance is R td = 2.50 ± 0.09 R p. By comparison, Halekas et al. (2017) found large variations of the bow shock standoff distances in the early MAVEN data, with R ss ∼1.6–1.9 R p and R td ∼2.5–3.1 R p depending on EUV flux levels and combined with either M ms or the solar wind dynamic pressure. Differences toward upper values with previous studies likely stem from different EUV levels encountered by the respective missions, the solar EUV flux being one of the main drivers, through ionosphere and exosphere variations, of the bow shock position (Halekas et al., 2017; Hall et al., 2019).

We present in Table 2 the fitted conic parameters of the main quoted references in Table 1, in chronological order. It is interesting to remark that most shapes fitted are stricto sensu hyperbolic (ϵ ≥ 1), but in practice can be considered quasi‐parabolic as eccentricity ϵ ∼ 1, which makes it possible to calculate the limiting Mach cone angle. For example, when fitting bow shocks from different MYs, Hall et al. (2019) showed that eccentricities varied around ϵ = 1 by <5% between MY27 and MY33, with a marked tendency toward ellipsoidal shapes (only two consecutive years, MY28 and MY29 had eccentricities above 1).

Table 2

Summary of Martian Bow Shock 2D Conic Parameters

Reference	ϵ	L [R p]	x F [R p]	R ss [R p]	R td [R p]	α	Nature	ϱ [°]
Russell (1977) a	0.99 ± 0.11	2.985	0	1.50 ± 0.15	3.00 ± 0.13	0	Ellipse	−
Slavin and Holzer (1981) b	0.94 ± 0.04	1.94 ± 0.02	0.5	1.50 ± 0.04	2.36 m	tan^−1V_pV_sw	Ellipse	−
Slavin et al. (1991) c	1.02	1.68	0.7	1.55	2.29 m	tan^−1V_pV_sw	Hyperbola	11.4 ± 2.9
Schwingenschuh et al. (1990)	0.85	2.72	0	1.47 ± 0.03	2.72 m	3.8°	Ellipse	−
Trotignon et al. (1991)	0.95 ± 0.10	2.17 ± 0.03	0.5	1.62 ± 0.07	2.60 m	4°	Ellipse	−
Trotignon et al. (1993)	1.02 ± 0.01	2.17 ± 0.03	0.5	1.57 ± 0.03	2.6	4°	Hyperbola	11.4 ± 2.8
Vignes et al. (2000) d	1.03 ± 0.01	2.04 ± 0.02	0.64 ± 0.02	1.64 ± 0.08	2.62 ± 0.09	4°	Hyperbola	13.9 ± 2.3
Vignes et al. (2000) e	1.02 ± 0.02	1.93 ± 0.01	0.72	1.67 ± 0.03	2.56 ± 0.06	4°	Hyperbola	11.4 ± 5.6
Trotignon et al. (2006)	1.026 ± 0.002	2.081 ± 0.006	0.6	1.63 ± 0.01	2.63 ± 0.01	4°	Hyperbola	12.9 ± 0.5
Edberg et al. (2008) f	1.05 ± 0.04	2.10 ± 0.09	0.55 ± 0.12	1.58 ± 0.18	2.69 m	4°	Hyperbola	17.8 ± 6.8
Hall et al. (2016)	1.01 ± 0.11	1.82 ± 0.08	0.74−0.10+0.03	1.65−0.18+0.13	2.46−0.22+0.20	4°	Hyperbola	− k
Halekas et al. (2017) g	1.0	2.01–2.54	0.6	1.6–1.9 l	2.5–3.1 m	4°	Parabola	−
Ramstad et al. (2017) h	1.022	1.48	0.85	1.58	2.19	(4°)	Hyperbola	11.9 ± 2.7
Hall et al. (2019) i	0.998 ± 0.001	1.802 ± 0.002	0.76	1.662 l	2.445 ± 0.003	4°	Ellipse	−
All (one per mission) j	1.016 ± 0.012	2.01 ± 0.25	0.61 ± 0.10	1.61 ± 0.08	2.56 ± 0.20	4°	Hyperbola	13 ± 4

Note. Pre‐Mars Express results were already summarized in Trotignon et al. (2006), Table 1. The aberration angle α is given for each reference. ϱ is the limiting Mach cone angle, calculated by formula (7) in the case of a hyperbolic shape. The mean value for each mission is also given, with Mars 2, 3, and 5 and Mariner 4 (Slavin et al., 1991), Phobos 2 (Trotignon et al., 1993), Mars Global Surveyor (MGS; Edberg et al., 2008), and Mars Express (MEX; with Hall et al., 2016, 2019). The planetary radius of Mars is by definition R p = 3,389.5 km.

Because the Mars 2, 3, and 5 measurements reported by Gringauz et al., 1976; in total 11 crossings did not specify local times, aberration angle α was assumed to be zero.

These authors use the full definition of the aberration angle, resulting in , in contrast to the more recent studies. See Section 2.1.

Mariner 4, and Mars 2, 3, 5 data only here. Uncertainties on ϵ fitted values assumed to be 1%.

“Direct fit” method with all three parameters varying simultaneously.

“Slavin’s method,” using a linear regression in (1/r, cos θ) space.

Note that ϵ = 1.03 matches better with Figure 1 of Edberg et al. (2008), for which the Mach cone aperture would instead be ϱ = 13.9.

Fits were performed on 2D‐gridded density data, co‐depending on M ms and EUV flux levels on the one hand, and solar wind dynamic pressure and EUV flux levels on the other.The coordinate system adopted by Halekas et al. (2017) was the Mars Solar Electric (MSE) system, with the X axis lying anti‐parallel to the solar wind flow.

Ramstad et al. (2017) use the following rectangular function (required to be cylindrically symmetric with respect to the solar wind direction): , with ρ the radial distance to the bow shock on the Yʹ axis from the center of Mars, R ss the subsolar standoff bow shock distance on the X′ axis, and ς the so‐called scale length. This function is valid ∀x ≠ R ss since ρ(y = 0) = R ss. By definition, . ς is a constant equal to 33.54 R p derived in Ramstad et al. (2017) from the bow shock model values for R ss, R td, and ϵ of Vignes et al. (2000) and can be calculated as . The original values of R ss and ϵ in their study were fitted to a function ; we have assumed here nominal conditions for simplicity. We calculate the semilatus rectum as from formulae (5) and (6). Uncertainties on ϵ fitted values assumed to be 1%.

Here, we only recall the results for all MYs (MY27–33). Individual MYs have eccentricities below 1 (ellipse), except for MY28–29 (hyperbola).

That is, Mariner 4 and Mars 2, 3, and 5 (Slavin et al., 1991), Phobos 2 (Trotignon et al., 1993), MGS (Edberg et al., 2008), MEX (Hall et al., 2019), and MAVEN (Halekas et al., 2017). The listed uncertainties are the standard deviations of the series. Accordingly mean angles ϱ are calculated only for 3 values and are only given for for completeness here.

Although this is a hyperbola with cone angle ϱ = 8.1°, the large eccentricity uncertainty leads to a cone angle uncertainty of 44°, hence no ϱ value is provided here.

Calculated from formula (5).

Calculated from formula (6).

The 3D Cartesian Form

The more general way of characterizing the bow shock shape does not assume any symmetry with respect to any axis. A 3D shape model can be constructed in the form of a quadratic equation (e.g., Formisano et al., 1979; Gruesbeck et al., 2018; Simon Wedlund et al., 2017, for Earth, Mars, and comets): Here and for clarity in the equations, (x, y, z) coordinates are by definition the unaberrated (X MSO, Y MSO, Z MSO) coordinates. With the MAVEN spacecraft including both magnetometer and ion spectrometer, Gruesbeck et al. (2018) used a limited subset of bow shock crossings when ignoring rapid spatial motions of the boundary across the spacecraft due to the variable solar wind (see Halekas et al., 2017), leaving a database of only 1,799 crossings spanning about 3 years of data (November 2014–April 2017). For all bow shock detections considered in their study, the best least‐squares ellipsoid fit was obtained with A = 0.049, B = 0.157, C = 0.153, D = 0.026, E = 0.012, F = 0.051, G = 0.566, H = − 0.031, and I = 0.019 and is valid only for the dayside bow shock up to a few 0.1 R p downstream of the terminator (X ≳ −0.5 R p) because of the poor MAVEN orbital coverage on the nightside flanks of the shock. Gruesbeck et al. (2018) concluded that asymmetry of the shock surface was particularly pronounced in the North‐South direction, likely due to the influence of crustal magnetic fields. Located predominantly in the southern hemisphere of Mars between 30° and 85°S (Acuña et al., 1998), they tend to increase the altitude of the induced magnetospheric boundary and hence increase the subsolar standoff distance of the shock (see Matsunaga et al., 2017 and references therein).

Although quadratic surfaces are not necessarily centered on the planet nor is their main axis directed along the X MSO axis (see Appendix A), a simple estimator of the shock’s position in the subsolar and terminator directions can be of interest. We propose here such an expression, based on Equation (9). We calculate thus the subsolar standoff distance along the X MSO axis at coordinates (x, y = 0, and z = 0) by finding the positive root of the simplified quadratic equation (i.e., the intersection of the surface with the X MSO axis): whereas the terminator standoff distances in the Y − Z plane (non‐aberrated MSO coordinates) are similarly given by:

Because of the small angles involved at Mars, non‐aberrated coordinates are rather accurate for the subsolar standoff distance. That said, rotating the MSO coordinate system by a certain small angle α about the Z axis does impact the terminator distances by a few 0.01 R p. For the parameters given above, R ss = 1.56 R p, R td, = 2.50 R p, and R td, = 2.62 R p. Because the shock is a 3D object, the exact position of the tip of the ellipsoid may vary with respect to the values taken at the origin.

Comparison of Historical Models

A comparison of a representative selection of historical bow shock models (some of them as listed in Vignes et al., 2000; Trotignon et al., 2006) in the Xʹ −Yʹ/Zʹ plane is given in Figure 3. The 3D quadratic model of Gruesbeck et al. (2018) was rotated anticlockwise by 4° around the Z axis to ease the comparison. It is noteworthy to remark that although the fit is not valid for Xʹ ≲ −0.5 R p, the figure displays the fits for Xʹ > −2.5 R p to illustrate the differences in shock surface swelling.

Figure 3

Bow shock fitted models to observations in MSO aberrated coordinates. The 3D quadratic model of Gruesbeck et al. (2018) fitted from MAVEN data was rotated anticlockwise 4° around the Z axis. All other models are obtained in cylindrical conic form from other missions, including Mars Express (MEX), Mars Global Surveyor (MGS), Phobos 2, and the Mars 2, 3, and 5 missions. The fits of Hall et al. (2016, 2019) and Gruesbeck et al. (2018) consider all shock detection points of their respective studies. Because the cylindrical models are symmetric about the Xʹ axis, the figure’s cylindrical y‐axis is equivalent to the Yʹ or to the Zʹ axis, regardless. The coordinates are normalized to the radius of Mars, R p = 3,389.5 km.

All models are in excellent agreement around the subsolar point, with a mean subsolar standoff distance value of R ss = 1.59 ± 0.05 R p. The terminator standoff distance is also in very good agreement—however, for Xʹ ≲ 0 R p, the difference between fits becomes substantial, especially (a) between the recent MEX investigations of Hall et al., 2016, 2019 and the other fits on the one hand and (b) between the MAVEN fits and all other fits on the other. For MAVEN, this is due, as previously mentioned, to the lack of orbital coverage by the spacecraft for Xʹ < 0 R p. In this sense, MEX has a much better antisolar spatial coverage. Of note, the North‐South asymmetry in the 3D fits of Gruesbeck et al. (2018), coinciding with the presence of crustal fields in the Southern hemisphere of Mars, can easily be seen, a characteristic which no axisymmetric model directly quantifies.

Shock Geometry

Quasi‐Perpendicular or Quasi‐Parallel Shock?

A collisionless shock may have different behaviors depending on the upstream solar wind magnetic field (the IMF), which conditions how the solar wind is losing its energy to the magnetosheath. Two main cases are conveniently studied for their varying properties: q ‖ and q ⊥ shocks. Additional important physical quantities driving the shock structure and dynamics are the magnetosonic Mach number (which defines the shock’s criticality) and the plasma‐β (Balogh & Treumann, 2013).

It is useful to recall that a q ‖ shock condition is defined so that the background IMF lines are intersecting normally the shock surface, whereas a q ⊥ shock describes an IMF that is in effect in the tangent plane to the surface shock. Thus, the angle of importance is the angle between the average IMF vector upstream of the shock and the shock normal. This angle is in the literature almost always named θ Bn, which is kept here for convenience. The geometry of the shock is defined as follows:

Starting in the magnetic field compression region in the solar wind, q ⊥ shocks have structures, from the point of view of B‐fields, almost always characterized by (a) a foot, (b) a fast ramp, and (c) a wider overshoot followed by a more gradual undershoot (see Kennel et al., 1985, Figure 11). This classic picture is a first approximation as fine electron‐scale structures in the foreshock, foot and ramp can be seen with high‐cadence magnetic field measurements. Q ⊥ shocks reflect particles back upstream to satisfy the shock conditions and are on average diffusive. Magnetic structures trapping particles such as mirror modes are observed to predominantly take place in the magnetosheath behind a q ⊥ shock (Gary, 1992). On the other hand, q ‖ shocks are on average resistive and are usually characterized by heavy turbulence. Their foreshock contains MHD turbulence that can give rise to first‐order Fermi acceleration. Also common in the foreshock region, highly compressive structures such as Short Large‐Amplitude Magnetic Structures (SLAMS) are associated with large density variations: they originate from the steepening of ULF waves and are of great importance in the shock reformation (Burgess et al., 2005; Burgess & Scholer, 2014).

In order to determine the q ‖ or q ⊥ geometry of the shock crossed by a spacecraft, the normal direction to the shock surface needs to be first estimated. For a single spacecraft, this can be done either with methods that take advantage of upstream and downstream magnetic field measurements (coplanarity method as in Horbury et al., 2002 and Schwartz, 1998, although prone to rather large uncertainties) or through geometrical considerations only, as we propose in Section 2.3.2. The accuracy of the geometrical method presented here is linked to the assumption that the shock surface is smooth and does not possess any kinks or local structures where the current curls on itself. In practice, we do not expect such a smooth surface as the shock may assume a more rippled shape which depends on the upstream solar wind condition and the turbulence at the boundary (Moullard et al., 2006). However, our geometric determination may still be a useful first approximation of the geometry of the shock.

Determination of the Shock Normal

The normal to the shock surface, at point (r 0, θ 0) in polar coordinates, (x 0, y 0, z 0) in Cartesian coordinates, or (r 0, ϑ 0, φ 0) in spherical coordinates, is simply defined as the gradient vector of the (assumed) smooth surface f at that point. Mathematically we can express this condition as: where v is a vector tangential to the surface at that point. This leads to the following expressions in the 2D and 3D cases.

The 2D Case

For the 2D polar coordinate case, let f be equal to f(r, θ) = r − L/(1 + ϵ cos θ) following Equation (1) where θ is the angle from the focus x on the Xʹ axis. The gradient of f depends on the two variables (r, θ):

At point (r 0, θ 0) vector ∇f = (R 0, Θ0) is perpendicular to the surface. Note that because of the peculiarity of a conic, values in x must always be corrected by the focus distance x , because the typical polar angle ϑ is not strictly the same as the conic angle θ used in Equation 1; see Figure 2.

In Cartesian coordinates, using Equation (4), the gradient will be against directions along Xʹ and Yʹ and equal to:

This expression circumvents the ambiguity on the angle direction of the polar formula, and as such should be preferred when calculating the normal direction. Figure 4a displays our estimates of the normal direction to several points on the shock surface applied to the 2D bow shock polar fit of Edberg et al. (2008), and converted to Cartesian coordinates.

Figure 4

Bow shock normal. (a) Polar 2D case (model of Edberg et al., 2008). (b) Cartesian 3D case (model of Gruesbeck et al., 2018). All coordinates are normalized to the radius of Mars, R p = 3,389.5 km.

The 3D case

For the 3D Cartesian quadratic equation, f is simply equal to the left member of Equation (9).

The gradient of f is then:

At point (x 0, y 0, z 0) vector ∇f = (X 0, Y 0, Z 0) is perpendicular to the surface. Additionally, the equation of the tangent plane to the smooth surface at that point is of the general form: where subscript “0” in the gradient components denotes the gradient taken at points (x 0, y 0, z 0) for brevity. Thus, knowing the 3D position of the spacecraft at the expected bow shock position, one can calculate the transverse and tangent directions to the bow shock surface. Figure 4b shows our results with this technique applied to the 3D bow shock fit of Gruesbeck et al. (2018), assuming a spacecraft situated on random points of the shock’s surface.

With the knowledge of the normal direction to the shock at the spacecraft location, it becomes possible, from the 2D and 3D model cases, to calculate the angle θ Bn from the average direction of the magnetic field at the spacecraft position near the shock, noted B bg:

For added robustness, the function arctan2 is recommended for the calculation of the inverse tangent, as it returns a value corresponding to the correct quadrant of the Euclidean plane. Because of the inherent 3D nature of a spacecraft orbit and of the local magnetic field, we prefer the 3D calculation with Equation (19) over the corresponding 2D case [Equation 18]. It is however important to recall here that the shock’s local shape may be assuming that of a “corrugated iron” section, as evidenced, for example at Earth, with the Cluster quartet of spacecraft (Moullard et al., 2006). No method is foolproof in estimating θ Bn: in case studies, the local normal to the shock can be more carefully checked, for which several complementary methods exist, such as the magnetic coplanarity method mentioned earlier.

When applied to the MAVEN dataset, the geometric calculation of θ Bn, assuming a smooth surface, is expected to reach an uncertainty of about ±5° depending on the upstream field determination. We obtained this estimate for a small sample of crossings, by extending over a few minutes the time spans used to calculate the upstream magnetic field direction.

Detecting the Bow Shock in Spacecraft Orbits

Estimating the bow shock position from spacecraft spatial coordinates can be achieved either empirically or theoretically, depending on the precision needed. Semi‐empirical but computationally intensive techniques using machine‐learning imaging algorithms are currently attempted to detect automatically and precisely the exact position of the shock. Such techniques make use of the full plasma instrumental payload on board planetary missions, when available. Other techniques like that of Němec et al. (2020) use ion and magnetic field data in combination to statistically identify the regions crossed by the spacecraft. However, a faster approach, based on a simple geometrical estimator using a static analytical bow shock model (see Section 2), may still prove valuable for statistical studies or for new datasets. We present such an approach and its possible refinements in 2D and 3D coordinate systems based on magnetic field‐only measurements. Because in the following we do not take into account plasma measurements, and because the true signature of the shock may be difficult to detect with magnetic field data only, a range of simultaneous criteria to identify solar wind and magnetosheath regions is required to mitigate this ambiguity. These criteria are presented in Section 3.2.

Predictor Algorithm for the Shock Position From Existing Analytical Models

The predictor algorithm is based on the calculation of polar (2D, θ) or spherical (3D, [ϑ, φ]) angles in the corresponding frame of reference centered on the planet in MSO coordinates. These angles unequivocally define the predicted distance to the bow shock at the position of the spacecraft (X sc, Y sc, Z sc). Comparing this bow shock distance with the Euclidean spacecraft distance to the center of the planet gives access to the region in which the spacecraft is located either in the solar wind or in the magnetosheath.

The 2D Case

Our algorithm is (see Figure 2 for definitions of angle and distances):

Choose conic model r(θ), with eccentricity ϵ, semilatus rectum L and focus’ position (x , 0, 0).

Calculate the spacecraft’s Euclidean distance r sc from the chosen conic model focus x F, in aberrated MSO coordinates, so that: .

Calculate the angle at the position of the spacecraft: .

Calculate the predicted bow shock distance R bs at the corresponding spacecraft θ sc from the focus x : R bs = r(θ sc) = L/(1 + ϵ cos θ sc) following Equation (1).

Calculate ΔR = r sc − R bs. If ΔR goes from negative to positive (respectively, from positive to negative) values, the spacecraft is expected to move from the magnetosheath to the solar wind (respectively, from solar wind to magnetosheath). At the temporal precision of the spacecraft ephemerides, the closest value to ΔR = 0 defines the shock position and crossing time t bs.

This purely geometrical approach was tested for polar coordinate models (such as Edberg et al., 2008; Hall et al., 2016, 2019), provided that all spacecraft coordinates are first rotated 4° into the aberrated MSO system.

The 3D Case

For 3D quadric models (Gruesbeck et al., 2018), the aberration is already taken into account and there is no need to correct the spacecraft coordinates for the position of the focus of the conic. Thus we only need generalize the approach above to spherical coordinates (ρ, ϑ, φ), where is the planetocentric distance, whereas and represent azimuth and elevation (measured from the X–Y plane) by convention. To compensate for the inherent ambiguity on azimuth depending on the quadrant and gain robustness, the function arctan2 is preferred throughout for simplicity, in a programing sense. Equation (9) becomes a second‐degree equation of the form: with:

In the case of an ellipsoid of revolution (as it is the case for the parametrization of Gruesbeck et al., 2018), the bow shock distance R bs at the angles (ϑ, φ) corresponds to the positive root of this equation:

for a ≠ 0. For other parametrizations such as a hyperboloid of two sheets, there may be two positive roots, in which case the smallest root should be chosen. The denominator 2a in expression (23) never reaches zero, no matter the combination of angles chosen, which makes it a robust formula throughout any orbit. For azimuth and elevation angles (|ϑ| ≳ 115°, |φ| ≲ 55°), the model of Gruesbeck et al. (2018) is not applicable any more (standoff distances above 4 R p) as these particular angular combinations correspond to a tail‐flank position, lying outside of MAVEN’s orbital range.

The next step is to determine whether the spacecraft is inside the bow shock surface or outside of it in the orbital sequence. Our algorithm follows a similar sequence as for the 2D case, but all variables are calculated with respect to the center of Mars, in unaberrated MSO coordinates:

Choose 3D model of the bow shock with A, B, C, D, E, F, G, H, and I parameters.

Calculate the spacecraft’s Euclidean distance in non‐aberrated coordinates, .

Calculate (azimuth, elevation) angles (ϑ, φ) at position of the spacecraft: and .

Calculate the bow shock distance R bs at the corresponding spacecraft spherical angles with Equation (23).

Calculate ΔR = R sc − R bs. If ΔR goes from negative to positive (respectively, from positive to negative) values, the spacecraft is expected to move from the magnetosheath to the solar wind (respectively, from solar wind to magnetosheath). At the temporal precision of the spacecraft ephemerides, the closest value to ΔR = 0 defines the shock position and crossing time t bs.

Application to MAVEN Orbits

Using 1‐min‐averaged MAVEN orbits to test our prediction algorithm, we perform the automatic detection of the bow shock as shown in Figure 5. Because of the relatively poor temporal resolution of this dataset, as well as the fast approach in the early stages of the orbit insertion, some points in the orbit yield false positive detections which disappear when increasing the orbital data resolution to 1 s.

Figure 5

Automatic detection of bow shock in MSO coordinates normalized to the planet’s radius, using the general quadric formula of Gruesbeck et al. (2018). The bow shock surface is in brown, the orbit of MAVEN between 1 November 2014 and 7 February 2021 is in blue. Detections of the crossings from inside the shock surface to outside of it are shown as orange circles, whereas outside‐to‐inside crossings are depicted by yellow circles. Coordinates are normalized to the radius of Mars, R p = 3,389.5 km.

Thanks to the simple algorithms presented above, we can statistically predict bow shock crossings in a given spacecraft orbit. To help identify the solar wind region, we can distinguish between trajectories moving from the magnetosheath to the solar wind region, and vice‐versa. For each orbit intersecting the bow shock model, two points per orbit will be identified. At 1 s resolution, we predicted a total of 16,515 bow shock crossings using the 3D analytical model of Gruesbeck et al. (2018) for the MAVEN dataset between November 2014 and February 2021, including 8,256 crossings from the sheath to the solar wind and 8,259 crossings from the solar wind to the sheath.

Actual crossings may in practice be different than the predictions and the algorithms may fail to pinpoint the location of the shock sometimes by several tens of minutes, mostly because of solar wind varying conditions, or due to the geometry of the shock at the point of passage for an individual orbit (Halekas et al., 2017). We estimate thus the precision of these automatic estimates to be of the order of ±0.08R p (±270 km) around the “true” bow shock location, from a representative subset of data. Because of variable shock position from orbit to orbit and the geometric average nature of the detection, some orbits that may have experienced shock crossings but lie inside the average shock position will not be tested for potential detection. It is estimated that only a few hundred potential crossings were ignored in the process. Consequently, the true shock structure location should be checked directly in the magnetic field and plasma data. Moreover, since the bow shock is a dynamical object, it may experience fast forward and backward motions, crossing the spacecraft trajectory several times more per orbit. This can be seen, for example, in Figure 6c, where the total magnetic field undergoes sharp intermittent jumps in the foreshock area (23 June 2018, around 01:00 UT). As in Gruesbeck et al. (2018), these multiple crossings, usually the first one in the temporal sequence for crossings into the sheath, and the last one for crossings into the solar wind, are identified as one with our algorithm.

Figure 6

Examples of total magnetic field amplitudes |B| at 1 s resolution measured by the MAVEN/MAG instrument throughout the mission (left y‐axes), and calculated running median absolute deviations ς mad,B (right y‐axes in green). (a) 25 December 2014 (beginning of mission). The first crossing is quite oblique (θ Bn ≈ 45°) followed by two highly q ⊥ shock crossings (θ Bn   85°). (b) 23–24 September 2016. The two detected crossings are q ⊥, the first one with θ Bn ≈  58°, the second with θ Bn  ≈  78°. (c) 23 June 2018, with two detected q ‖ crossings (θ Bn ≈ 8, 25°). (d) 28 March 2020, with five crossings all oblique toward q ⊥ conditions, with θ Bn ≈ 45°, 49°, 80°, 82°, and 88°, successively. The predictor geometric detections (Section 3.1) are in dashed lines and labeled “automatic,” whereas the predictor‐corrector detections proposed in Section 3.2 are in solid lines. Highlighted in different colors are crossings from solar wind to magnetosheath (labeled SW  →  MS, red) and from magnetosheath to solar wind (labeled MS  →  SW, blue). Calculations of θ Bn angles were performed using median averages of B over the color‐highlighted regions (blue for MS  →  SW crossings, red for SW  →  MS crossings). The threshold ς th = 0.5 is shown as a horizontal dashed line (right y‐axis, green).

We successfully applied this 3D algorithm to the retrieval of undisturbed solar wind density and velocity moments in the MAVEN/SWIA data with 1 min resolution, as part of the Helio4Cast solar wind in‐situ data catalog, enabling the statistical study of interplanetary coronal mass ejections and high speed streams (Möstl et al., 2020).

Refining the Position of the Shock: A Predictor‐Corrector Algorithm

Because of variations in the shock position, the automatic detection may give inaccurate predictions. We present here a fast method to correct to a certain extent for these discrepancies. It makes use of the magnitude of B to identify the position of the shock structure, either when crossing from the magnetosheath into the solar wind or vice‐versa. As before, the main assumption is that the shock is crossed twice per orbit at maximum, although in practice the shock structure may be crossed several times due to the fast motion of the boundary across the spacecraft trajectory (Halekas et al., 2017). As we are interested in the statistical position of the shock, this assumption nonetheless provides a valuable estimate of the average position of the shock during those times. All magnetic field data are assumed here to be of the order of 1 s resolution.

From the point of view of a single spacecraft’s magnetic field measurements across the shock boundary, the total magnetic field intensity increases sharply, from typically 5 to 10 nT in the solar wind at Mars to about twice that level on average when moving into the magnetosheath. Additionally, fluctuations increase, going from small standard deviations around the mean in the solar wind, to comparatively larger fluctuations in the magnetosheath. Typical crossings of the Martian shock illustrating these behaviors are shown, for example, in Figure 6 which presents examples of bow shock crossings around Mars as seen with the MAG instrument on board MAVEN throughout the mission. As discussed in Section 2.3, q ⊥ shocks usually display distinct features (foot, sharp ramp, and overshoot) as in Figure 6a, 6b, and 6d, whereas q ‖ geometries have less clear signatures in magnetic field data, making it difficult to detect reliably (Figure 6c). Consequently, the corrector method presented here is biased toward the detection of q ⊥ crossings.

Our predictor‐corrector detection algorithm attempts to consistently identify solar wind undisturbed regions in the chosen dataset, characterized by relatively low magnetic field intensities in combination with small fluctuations. For each orbit of a spacecraft, it proceeds as follows, with user‐defined values (Δt bs, B th, T, and ς th) discussed afterward:

Calculate predictor estimate of the shock crossing time t bs with the automatic algorithm in 2D aberrated polar coordinates or in 3D (as in Section 3.1).

If t bs exists in the considered orbit, choose a small time interval symmetric around the estimated shock so that [t bs − Δt bs, t bs + Δt bs] (user‐defined).

Calculate a robust estimate of the average magnetic field, for example, the median of the magnetic field in half of this interval, noted |B sw|1/2. By definition |B sw|1/2 corresponds to the assumed solar wind region magnetic field. For crossings from the solar wind into the sheath, the first half‐interval is selected. For crossings from the sheath to the solar wind, the second half‐interval is selected. There are two possibilities at this junction:

If |B sw|1/2 > B th, where B th is the solar wind‐to‐magnetosheath threshold (user‐defined), the position of the solar wind region is difficult to assess. In that case:

Make the time interval float around the estimated location of the shock, by increments of Δt bs/3 in one direction or the other, so that t bs = t bs ± Δt bs/3,

Repeat interval shift until |B sw|1/2 ≤ B th or until a maximum shift of 2Δt bs is effected from the estimated shock timing in either direction. If |B sw|1/2 > B th still, either the spacecraft is always in the sheath or in an usually high solar wind B‐field region, in which case the crossing is altogether ignored and removed from the database. If not, go to step (b):

If |B sw|1/2 ≤ B th, this half‐interval is a good candidate for undisturbed solar wind conditions. In that case:

Calculate the running Median Absolute Deviation (MAD) ς mad,B of the total B‐field signal over a temporal window of duration T (user‐defined) in the chosen half‐interval, and smooth further the result with a running median over a time span enclosing the shock structure in its entirety (e.g., 2T or 3T). This also helps remove potentially abrupt but temporally isolated changes in the signal. Note that this particular choice of ς mad,B is somewhat arbitrary. After several tests including running standard deviations, normalized or not to the “solar wind” signal, the choice of a smoothed running MAD was empirically found to work consistently well with the MAVEN dataset at 1s resolution. For all times t at which the total magnetic field B tot is measured over a running interval [t , t  + T]:

Compare ς mad,B to threshold value ς th (user‐defined). A jump above a certain threshold ς th indicates a transition between a less turbulent region (ς mad,B < ς th) to a more turbulent one (ς mad,B > ς th), an indicator of the presence of a shock‐like structure. If this threshold is reached, take the first (respectively, last) time this happens in the chosen interval for solar wind‐to‐sheath (respectively, sheath‐to‐solar wind) crossings, and correct the original timing t bs of step 1 to new t bs, corr. If not, discard crossing.

Repeat for each orbit.

At Mars, we tested step 1 (predictor) in the previous section with the MAVEN dataset: on average, the detected shock was within ±0.08R p (±270 km) of the true shock crossing, corresponding to about ±30 min of data along the orbit. We thus set the interval of study for the corrector algorithm to Δt bs = 30 min. This value depends on the orbit inclination with respect to the bow shock surface and is thus mission‐dependent.

We determined the typical user‐defined threshold values for the MAVEN mission manually, for simplicity, on a reduced dataset. We found a good compromise by trial and error with B th = 11 nT, because the undisturbed solar wind magnetic field in Mars’ vicinity is of the order of 2–6 nT on average (Slavin & Holzer, 1981), but can reach up to about 10 nT or more when solar transient effects such as coronal mass ejections or co‐rotating interaction regions are involved (Liu et al., 2021). In future studies, a more dynamic criterion in step (3a) may be preferred, that is, where the amplitude of the magnetic field is normalized to the assumed upstream solar wind value. The criterion for being in the magnetosheath could, for example, become , where γ th is an adequately chosen threshold (γ  ∼  1.5 for a clear increase of magnetic field when moving into the magnetosheath, with nominal solar wind levels γ  ∼  1).

Because the shock appears in measurements as a turbulent structure whereas the solar wind is on average less so, step (3b–i) calculates a measure of the variability of the magnetic field in the vicinity of the shock. The value of T is also mission‐ and instrument‐dependent; for MAVEN/MAG data at a resolution of 1 s, a value T = 120 s was chosen, which adequately captures the mean magnetic field variations. When in the solar wind, we found that the intrinsic variability of the signal ς mad,B, calculated over a running window of duration T = 120 s, was on average <0.5, which is adopted as the threshold ς th. This makes it possible to detect the very first perturbations in the solar wind leading to the creation of the shock structure, a point which is identified here as the position of the shock proper at time t bs. Again, in future studies, the threshold can be normalized to the magnetic field level, as in Halekas et al. (2017) where, together with constraints on plasma parameters, the normalized root‐sum‐squared value of the magnetic field was chosen so that RSS(B)/|B| < 0.15 to identify undisturbed solar wind intervals.

Application to MAVEN Dataset

Applied to the MAVEN dataset at 1 s orbital and magnetic field resolution, the corrector algorithm reaches an accuracy of ±0.02R p (±70 km) around the “true” shock (manually picked on a reduced dataset for comparison), a factor 4 increase in accuracy with respect to step 1. In the temporal datasets, this corresponds to only a few minutes of continuous data along MAVEN’s orbit. This is epitomized in Figure 6, which shows examples of shock crossing predictions from step 1 (dashed lines) in magnetic field data (left axis) as compared to the predictions from the predictor‐corrector algorithm (unbroken lines, red for crossings into the magnetosheath, and blue for crossings into the solar wind). Right axes in green show ς mad,B and the associated threshold of 0.5, where ς mad,B < 0.5 mainly occurs for upstream solar wind intervals. With the addition of criterion B th on the total magnetic field, only clear solar wind regions are captured by our algorithm, whereas ambiguous regions from the point of view of the magnetic field data are rejected.

At the beginning of the mission (Figure 6a, December 2014), the automatic predictor algorithm gives a reliable estimate of this q ⊥ shock’s position: this is expected since the prediction is based on the 3D quadric model of Gruesbeck et al. (2018) who specifically performed shock fits on the first years of the mission. In this case, the predictor‐corrector algorithm only corrects the shock’s estimated location by a few minutes. Figure 6b (September 2016) displays a case where the shock position is hard to ascertain from magnetic‐field data only: the predictor estimate is off by up to 20 min for all crossings. Because of the constraints on the magnetic field amplitude and the lack of significant variations in |B| between 22:00 and 01:00 UT, the predictor‐corrector algorithm ignores the two first expected crossings but corrects well for the two next crossings (around 03:30 and 04:45 UT). Figure 6c shows a more complex mix between sheath and solar wind conditions, and even though the boundaries are more subtle and the overall B‐field magnitude below 10 nT, the predictor‐corrector algorithm manages to estimate the position of the q ‖ shock well, ignoring possibly unclear crossings which do not fulfill the combined threshold conditions on ς mad,B and B th (around 01:00 UT). Figure 6d shows yet another example of the superiority of the predictor‐corrector algorithm for some very clear bow shock crossings in 2020, after the orbit of MAVEN had been altered into a different orbit than at the beginning of the mission.

The final corrected timings for the detections yield with this algorithm a lower estimate of the total actual number of crossings encountered by a spacecraft throughout its mission. Events occurring when |B sw|1/2 > B th, even after shifting the temporal window significantly, or when ς mad,B < ς th, were discarded in the final selection as can be seen in Figure 6. They may indicate that the magnetic field was either too turbulent or too complex in its structure (e.g., multiple crossings as is regularly the case with q ‖ crossings) for the corrector algorithm to capture. Moreover, as discussed previously, the analytical approximation model used for the determination in step 1 (either 2D or 3D) is likely to underestimate the true number of crossings due to the planetary bow shock variability (Halekas et al., 2017). Because our study is primarily interested in the statistical position of the bow shock throughout the mission, this loss of potential detections may be compensated by the large number of orbits of the considered spacecraft.

For MAVEN, from the original 16,515 candidate detections from step 1 (predictor) using the 3D model of Gruesbeck et al. (2018) as a first approximation, our predictor‐corrector algorithm selected 14,929 events (7,494 detections from the solar wind to the magnetosheath and 7,435 from the sheath to the solar wind) for the period 01 November 2014–07 February 2021. This is a 10% decrease in number of crossings, leaving out the less ambiguous events from the predictor algorithm. On average, the correction to the original timing is within about Δt ± 20 min along the spacecraft orbit, which corresponds to a percentage difference in radial distance of about ±25% (not shown). The calculated average difference |ΔR bs| over the entire dataset between predictor and predictor‐corrector algorithm is ∼0.11R p, that is, ∼380 km or about 5% difference. The complete list of MAVEN crossings (times and spatial coordinates) derived with the predictor‐corrector algorithm is compiled in Simon Wedlund et al. (2021).

Out of these 14,929 crossings, we found that an overwhelming number (11,967) was predominantly q ⊥ crossings with θ Bn > 45° as calculated with formula (19), the remaining 2,962 events being classified as more reminiscent of q ‖ shock conditions. It is important to recall that the method is biased toward clear signatures of q ⊥‐like bow shock crossings, in practice filtering out many q ‖ crossings, which may explain part of the discrepancy. About 22% of the q ⊥ crossings are highly perpendicular shocks (θ Bn ≥ 80°, 2,674 crossings), whereas only ∼1% of the q ‖ shocks are highly parallel ones (θ Bn ≤ 10°, 30 crossings). This is shown in Figure 7, where the highest number of crossings occurs for θ Bn ∼ 80°. This is also in qualitative agreement with the results of Vignes et al. (2002) for the MGS mission, when they investigated a proportion of 93 q ⊥ shocks for only 23 q ‖ shocks.

Figure 7

Statistical distribution of crossings with respect to θ Bn angles, the angle between the normal to the shock and the average magnetic field direction. The limit between q ‖ and q ⊥ conditions is for θ Bn = 45°.

Application to the Study of the Martian Bow Shock Variability

In this section, as an application and statistical test of our new predictor‐corrector algorithm, we present 2D and 3D fits of the average bow shock position with the MAVEN spacecraft for the 1 November 2014–7 February 2021 period. First, we sort out the detected shock positions (compiled in Simon Wedlund et al., 2021) by MY (MY32–MY35), aerocentric solar longitude Ls range (four seasons centered on equinoxes and solstices), EUV flux (two regimes, one for higher solar flux and one for lower solar flux), and shock geometry (q ⊥ or q ‖). These cases correspond to:

Mars Years: MY32 (incomplete), 33, 34, and 35, inspired by the work of Hall et al. (2019) on MEX datasets.

Solar longitude Ls ranges from [315° to 45°] (centered on Northern Hemisphere [NH] spring equinox), [45°–135°] (NH summer solstice), [135°–225°] (NH autumn equinox), [225°–315°] (NH winter solstice). Ls defines the geographic Martian season, with Ls = 251° (Ls = 71°, respectively), marking perihelion (aphelion) conditions.

Two EUV flux levels, inspired by the works of Halekas et al. (2017) and Gruesbeck et al. (2018) on the early MAVEN datasets. The EUV flux at Mars is obtained using the FISM‐IUVS daily irradiances at 121.5 nm calculated from the Mars EUVM model (Thiemann et al., 2017). The median of the EUV flux in the 2014–2021 period is 0.002 8 W m−2 and defines two EUV flux levels, one “high” for fluxes strictly above that limit, one “low” for fluxes below.

In Section 4.1, we perform 2D and 3D fits to the found Mars bow shock positions using the spacecraft ephemerides. A discussion of these fits and what they imply is given in Section 4.2.

Statistical Position of the Martian Bow Shock

In 2D MSO aberrated coordinates, the polar equation, Equation (1), can be rewritten in the linear form y = ax + b (see Trotignon et al., 2006): with a linear regression in the (r, Xʹ − x F) space performed for a chosen focus location x F. First we chose a focus location randomly between 0 and 1R p, and for each linear fit performed, the residuals are calculated. The adopted focus point is the one that minimizes the residuals. Because MAVEN’s orbits are not suited to bow shock detections for Xʹ < −0.5 R p, additional constraints on the tail distributions are necessary to obtain a more realistic conic fit. This can be achieved, for example, by using the predictions from a chosen pre‐existing model for deeply negative Xʹ values, such as those of Edberg et al. (2008, noted “E08” in the following) or Hall et al. (2019, noted “H19”) where bow shock detections were reported downstream to and to , respectively. First, additional “ghost” points (representing 10% of the total number of detections for the considered case) are calculated for for the chosen model and randomized spatially around this result to give a more realistic tail spread. The linear regressions are then performed on the new constrained dataset. Tests were performed on the robustness of this method using different analytical models: E08 and H19 fits are essentially the same around the nose of the shock downstream to about −1 R p where patent differences start to appear. Because of the added cloud of ghost points, this is expected and thus 2D fits presented below are only valid in practice in the range [− 0.5 R p, R ss]. Incidentally, differences on the terminator and subsolar standoff distances are less than in each case. It is noteworthy to add that the determination of the nature of the conic section found from the fits can be significantly altered when using tail models either from E08 (hyperbola) or from H19 (fits’ nature given by these authors, ranging from ellipse to hyperbola, depend on the MY considered): in that case, the fit's nature will naturally be biased toward matching that of their respective parent tail model.

Table 3 and Figure 8 display our 2D fits when the E08 tail model supplements the MAVEN dataset for additional constraints on the tail. Candidate bow shock detections are also drawn as semi‐transparent circles. Despite the constraints on the predictor‐corrector algorithm, several detected points appear to fall well into the magnetosheath of Mars, and are false detections. Because of their relative scarcity and thanks to the large statistical database, these points do not significantly impact the final fits, which remain robust.

Table 3

Martian Bow Shock 2D Conic Parameters in Aberrated MSO Coordinates From Linear Regression Fits Applied to Equation (25) and the MAVEN Orbits and Magnetic Field Data (Predictor‐Corrector Algorithm)

Case	ɛ	L [R p]	x F [R p]	R ss [R p]	R td [R p]	R 2	Nature	ϱ [°]	# detections
All points, this work	1.00	1.75	0.86	1.74	2.46	0.98	Parabola	−	14 ,929
MY32, this work	0.83	2.24	0.65	1.87	2.70	0.97	Ellipse	−	1,196
MY33, this work	0.99	1.88	0.75	1.69	2.51	0.98	Ellipse	−	4,586
MY34, this work	1.02	1.72	0.84	1.69	2.44	0.96	Hyperbola	11	5,073
MY35, this work	1.02	1.63	0.91	1.72	2.39	0.98	Hyperbola	11	4,074
Ls = [315°–45°], this work	1.01	1.73	0.86	1.72	2.45	0.98	Hyperbola	8	3,793
Ls = [45°–135°], this work	1.00	1.81	0.71	1.61	2.42	0.99	Parabola	−	3,746
Ls = [135°–225°], this work	0.99	1.82	0.71	1.62	2.42	0.98	Ellipse	−	3,134
Ls = [225°–315°], this work	0.98	1.91	0.86	1.82	2.62	0.98	Ellipse	−	4,256
EUV flux ≥0.0028 W m−2	1.00	1.79	0.91	1.80	2.54	0.98	Parabola	−	6,502
EUV flux <0.0028 W m−2	1.00	1.75	0.79	1.67	2.41	0.98	Parabola	−	8,427
Quasi‐⊥	1.00	1.79	0.82	1.72	2.48	0.98	Parabola	−	11 ,967
Quasi‐‖	1.06	1.47	1.07	1.78	2.37	0.94	Hyperbola	19	2,962

Note. Subsolar and terminator standoff distances R ss and R td are calculated with Equations (5) and (6). For hyperbolae, the Mach cone aperture ϱ is also given as calculated by Equation (5), (6). For each fit, the coefficient of determination R 2 gives a measure of the goodness of the linear regression. Due to the large data spread, uncertainties on R ss and R td are of the order of 5% and of the order of 2% for the other quantities. R p = 3,389.5 km is the radius of Mars.

Figure 8

The 2D fits performed on the predictor‐corrector algorithm for the detection of bow shock crossings in the MAVEN dataset, 2014–2021 in aberrated MSO coordinates (), and parametrized in Table 3. (a) and (b) Versus MY 32–35. (c) and (d) Versus Ls (season) ranges. (e) Versus EUV flux levels. (f) Versus shock geometry (q ⊥ and q ‖). (g) All detected points in the current database color‐coded by year, and comparison to the analytical quadric fit of Hall et al. (2019). All coordinates are expressed in units of the planet’s radius, that is, R p = 3,389.5 km. Superimposed on all panels are the corresponding analytical models of Hall et al. (2019) for MY 27–33, except for MYs 32 and 33, where their corresponding yearly fits are plotted. Candidate detections points for each case are also drawn as filled circles of varying colors, with the opacity giving a measure of the density of points in that area, giving more or less weight to the fitting method.

A 3D Case

In 3D, we perform quadric fits using the method first put forward by Taubin (1991), adapted to the quadratic surface of Equation (5), (6). This fitting method constructs scatter matrices from local gradients S of tested model T and finds the diagonal matrix of the generalized eigenvalue problem so that Tv = λ Sv, where v is the generalized eigenvector of T and S, and λ are the eigenvalues. Because of the scatter of points in the database, uncertainties on the found parameters A to I are of the order of 1%, in a least‐squares sense.

Table 4 collects all 3D fit parameters for each case; all fitted surfaces are ellipsoids of revolution. For completeness, we present and give the physical interpretation of these parameters in Appendix A, in terms of principal axes, their direction and lengths and the centering of the ellipsoids. Figure 9 shows the corresponding fits and their sections onto the X MSO − Y MSO (dawn‐dusk hemispheres), X MSO − Z MSO (South‐North hemispheres) and Y MSO − Z MSO (at the terminator, i.e., X MSO = 0) non‐aberrated MSO coordinates.

Table 4

Martian Bow Shock 3D Conic Parameters From Quadric Surface Fits Applied to the MAVEN Orbits and Magnetic Field Data (Predictor‐Corrector Algorithm)

Case	A	B	C	D	E	F	G	H	I	R ss	R td,y	R td,z
Gruesbeck et al. (2018)  a	0.049 0	0.157 0	0.153 0	0.026 0	0.012 0	0.051 0	0.566 0	−0.031 0	0.019 0	1.557	2.624	2.495
All points, this work	0.176 9	0.160 9	0.155 9	0.005 7	0.004 4	0.028 1	0.377 3	−0.032 3	0.014 3	1.539	2.595	2.487
MY32, this work	0.136 9	0.141 9	0.140 0	0.038 1	0.017 8	0.054 7	0.378 3	0.004 4	0.027 9	1.654	2.639	2.575
MY33, this work	0.166 0	0.151 6	0.150 1	0.023 5	−0.006 1	0.020 3	0.380 7	−0.028 4	−0.001 8	1.562	2.664	2.587
MY34, this work	0.171 9	0.174 2	0.155 1	−0.011 2	−0.006 8	0.016 5	0.395 5	−0.025 6	0.012 7	1.522	2.470	2.499
MY35, this work	0.557 7	0.224 5	0.200 0	−0.050 9	−0.010 3	0.096 3	−0.142 1	0.001 1	−0.011 9	1.472	2.108	2.266
Ls = [315°–45°], this work	0.155 4	0.162 5	0.158 7	−0.026 5	−0.008 1	0.002 3	0.426 0	0.001 3	0.004 4	1.513	2.477	2.496
Ls = [45°–135°], this work	0.171 9	0.176 1	0.155 5	−0.027 5	0.004 8	0.056 4	0.410 7	−0.028 7	0.005 7	1.497	2.466	2.518
Ls = [135°–225°], this work	0.249 0	0.147 3	0.262 4	0.040 9	0.041 5	0.109 8	0.236 8	−0.094 1	−0.151 0	1.584	2.945	2.261
Ls = [225°–315°], this work	0.155 9	0.148 4	0.140 0	0.004 7	0.022 7	0.040 0	0.358 3	−0.007 2	−0.019 7	1.632	2.620	2.744
EUV flux ≥ 0.002 8 W m−2	0.109 6	0.148 0	0.150 0	0.027 4	0.003 1	0.035 5	0.432 9	−0.029 3	0.004 5	1.634	2.700	2.567
EUV flux <0.0028 W m−2	0.213 8	0.180 7	0.157 7	−0.023 1	−0.005 1	0.031 4	0.347 3	−0.020 7	0.018 5	1.498	2.410	2.460
Quasi‐⊥	0.179 8	0.160 7	0.153 9	0.001 6	0.004 0	0.033 0	0.377 7	−0.034 8	0.012 4	1.531	2.605	2.509
Quasi‐‖	0.142 7	0.167 5	0.166 6	−0.000 4	0.005 1	0.005 0	0.399 2	−0.009 8	0.023 0	1.595	2.473	2.382

Note. See Equation (9) for the definition of parameters A to I and Equations (11), (12), and (13) for those of the subsolar standoff distance R ss along the X MSO axis and the terminator standoff distances R td along the Y MSO and Z MSO axes. Uncertainties on the parameters are of the order of 1% in a least squares sense. All quadrics here are ellipsoids. The domain of validity for each fit is shown in Figure 9: fits are valid for X MSO ≥ −0.5 R p on average. The number of fitting points used for each case is the same as for the 2D fits, see Table 3 (last column). Also, see Appendix A for a physical interpretation of the tabulated parameters. R ss and R td are expressed in units of Mars radius R p = 3,389.5 km.

For all points considered in their data subset.

Figure 9

The 3D fits performed on the refined predictor‐corrector algorithm for the detection of bow shock crossings in the MAVEN dataset, 2014–2021 in the X MSO − Y MSO, X MSO − Y MSO, and Y MSO − Z MSO planes (traces of ellipsoids of revolution parametrized in Table 4). (a) Versus Mars Years (MY). (b) Versus Ls (season) ranges. (c) Versus EUV flux levels, with their corresponding subset of detected points (blue and orange dots). (d) Versus bow shock geometry, q ⊥ (blue dots) and q ‖ (orange dots). (e) All detected points in the current database with a comparison of present fit (black line) to the analytical quadric fit of Gruesbeck et al. (2018, orange dashed line). On each figure, superimposed crosses show where the nose of the shock is located, in the plane of projection (see Appendix B). All coordinates are expressed in units of the planet’s radius, that is, R p = 3,389.5 km.

Because of the spacecraft changing orbits during the mission, some of the ellipsoid fits appear anomalous in their orientation. This is especially obvious for MY35 when MAVEN, as of March 2020, decreased its apogee to ∼4,500 km and hence its revolution period to 3.5 hr to accommodate Mars 2020 rover operations on the ground. Consequently, MAVEN only seldom explored regions below X < 0 R p for half of MY35: this makes it difficult to constrain the fit, and we end up with an ellipsoid having its longest principal axis unphysically tilted almost 90° in the X‐Z plane (purple curves of Figure 9a). A similar issue is found for Ls = 135°–225° (yellow, Figure 9b), which has the lowest number of detections among Ls ranges and for which the orbit was never favorable for detections due to orbit precession. Thus, in these two cases, no physical interpretation should be drawn from the axes orientations of the ellipsoid and the fit should only be valid for near‐subsolar crossings. More robust physics‐based analytical models could be used to overcome these fitting issues (Kotova et al., 2021).

Discussion

Our 2D and 3D fits give some insight on how the Martian bow shock is moving globally for different conditions of MY, Ls, and EUV flux and complements previous studies with the MAVEN mission (Gruesbeck et al., 2018; Halekas et al., 2017; Němec et al., 2020). As the bow shock position is connected to the balance between thermal pressure from the plasma in the ionosphere and the dynamic pressure from the solar wind, any variation of these two quantities will have repercussions on the position of the shock.

When assuming axisymmetry around the aberrated axis X′ in the 2D polar rectangular coordinates case, Table 2 and the average subsolar and terminator distances can be a first guide for our interpretation. Our new results with MAVEN (all points, Table 3) agree rather well with past measurements (Table 2) considering the data spread and estimated uncertainties: R ss = 1.74 ± 0.09 compared to 1.61 ± 0.08 R p and R td = 2.46 ± 0.13 compared to 2.56 ± 0.20 R p. More specifically, for:

Mars Years: The subsolar standoff distance decreases by as much as 10% between MY32 and MY33–MY34, from 1.87 to 1.69 R p, although some of this variation may be stemming from the relatively lower statistics for the first year (1,196 points for MY32 compared to for all other years), due to the MAVEN mission starting toward the end of MY32. A similar tendency is seen for terminator standoff distances, with a 11% decrease seen between MY32 and MY35. Following Hall et al. (2019), these variations may be connected through solar EUV irradiance to the solar cycle itself, when descending from the maximum of solar cycle 24 (encompassed by MY32) toward a minimum of activity (MY34) and the start of solar cycle 25 (MY35). A variation in R ss of similar magnitude ( between minimum and maximum of activity) was shown by Hall et al. (2019) using MEX data for the previous solar cycles (23–24).

Seasonal variations: In contrast, the Ls ranges have a much more even statistics throughout, with more than 3,000 detections per season. Arguably, this makes comparing results between seasons statistically more significant than for the previous case. Overall, for northern spring equinox (Ls = [315°–45°]) and winter solstice (Ls = [225°–315°]) conditions, the bow shock appears to expand in the subsolar direction by about 7%–13% from its summer and autumn position (R ss ≥ 1.72 R p compared to R ss ≈ 1.61 R p). Simultaneously, the area encompassed by the bow shock conic is also increased during those two instances. One possible driving factor behind these changes may be in turn linked to changes in Mars’ dayside upper atmosphere and extended exosphere, and how they expand and contract with seasons, increasing or decreasing the size of the obstacle to the solar wind flow (Hall et al., 2016 and references therein). A denser lower atmosphere around perihelion (Ls° ∼ 251° where the EUV flux is highest on average) and during the dust storm season in the autumn (Trainer et al., 2019) may drive the ionosphere to expand significantly at constant EUV flux (Dubinin et al., 2019; Sánchez‐Cano et al., 2016), offering a more efficient obstacle to the solar wind. Similarly the expansion of Mars’ extended exosphere (notably modulated by the solar wind flux) increases the efficiency of the solar wind charge‐exchange process (with a net conversion of fast solar wind ions to slow‐moving heavy ions of planetary origin, effectively slowing down the solar wind, see Edberg et al., 2009; Halekas et al., 2017). Both aspects result in the standoff distance moving outwards. The opposite effect is expected when upper atmosphere densities are lower in the deep summer and in the beginning of the autumn and the bow shock surface shrinks. The fits and characteristics of the shock appear consistent with this picture.

EUV flux variations: The effect of a relatively larger flux on the shock position is twofold, globally increasing the ionization rates in the ionosphere and through photoionization of the extended exosphere as well as heating up and expanding the neutral atmosphere‐exosphere system (Edberg et al., 2009; Forbes et al., 2008; Hall et al., 2016). Photoionization of exospheric neutrals creates newly born ions that are picked up by the solar wind convective electric field, resulting in mass‐loading and slowing down of the solar wind flow (Yamauchi et al., 2015, with the presence of pickup ions in the foreshock region). Such combined effects have been shown to expand the bow shock in the solar wind direction (Mazelle et al., 2004). The two fits we present here, one for higher and one for lower EUV fluxes (more than 6,500 points each), display the expected behavior, with a larger standoff distance by 7% and a noticeably larger flaring of the fitted conic for the higher EUV fluxes (terminator distances increasing from 2.41 to 2.54 R p, i.e., 5%).

Shock conditions: q ‖ and q ⊥ bow shock crossings are related to the average interplanetary magnetic field’s (IMF) direction and the spacecraft’s orbit (more precisely, the spherical quadrant in which the spacecraft emerges into the solar wind). Because the predictor‐corrector algorithm favors q ⊥ detections (Section 3.2), the statistics between the two cases is heavily unbalanced (see Vignes et al., 2002, for a similar result). On average, we find no significant difference between the two conditions, with the shock surface slightly contracting and flaring up in q ⊥ conditions with respect to q ‖ conditions (|ΔR ss,td| ∼4%). Such a tendency is marginal considering that these percentages are at the precision limit obtained with the fits.

Let us now look at our 3D fit results. Figure 9 clearly shows several asymmetries depending on the MY, Ls, EUV flux’ and shock condition. The usual pronounced North‐South asymmetry (X MSO − Z MSO plane, second column, and also Y MSO − Z MSO plane, third column), mostly ascribed to the presence of crustal magnetic fields in the southern hemisphere (Gruesbeck et al., 2018), is clearly seen for all cases with the standoff subsolar distances being skewed toward that hemisphere. This is shown by crosses representing the tip of the projected ellipsoid (calculated by the formulae in Appendix B) located all in the fourth quadrant in the X MSO − Z MSO plane. A similar tendency is sometimes marginally observed in the X MSO − Y MSO plane (first column), when the shock surface is skewed toward the dawn hemisphere (−Y MSO), with standoff subsolar distances on average larger than on the dusk hemisphere. This is true for MY32 and MY33 (Figure 9a) and for larger EUV fluxes (Figure 9c). For lower EUV fluxes (and, incidentally, all other cases), the opposite seems to be taking place with the position of the maximum standoff distance being in the dusk +Y MSO hemisphere. It is difficult at this stage to tell if these latter (small) effects may stem mainly or not from the dawn‐dusk asymmetry of the atmosphere and hence of the ionosphere (Gupta et al., 2019). Likewise, as in the 2D case, it is important to note that many drivers of the shock position (EUV flux, atmospheric seasons, etc.) all act in combination at any given time: our fits do not discriminate precisely between these effects. A finer characterization of each driver separately is left for another study.

A comparison between all standoff distances, subsolar and terminator alike, and calculated by our 2D and 3D algorithms is shown in Figure 10, and based on Tables 3 and 4. The standoff distances calculated from the 3D fits show the same tendencies as their 2D counterparts, although since the X and Y coordinates are not solar‐wind aberrated and hence no axisymmetry is considered, the comparison between the 2D and 3D cases can only be that of general trends. From MY32 to MY35, a general decrease of standoff distances can be seen. Excluding R td,, the other standoff distances first steadily decrease from Ls centered on 0° (labeled “Ls1”) to 180° (“Ls3”) but then increase significantly toward Ls values around 270° (“Ls4”), which may be linked to the EUV flux becoming maximum at perihelion Ls = 251°. This result is arguably in contrast to those presented in Vignes et al. (2002) although our statistics with MAVEN is much larger than in their study. In an identical way to the 2D fits, larger EUV fluxes result in a bow shock surface significantly expanding in the solar wind toward the subsolar direction. With respect to the geometry of the shock, the subsolar standoff distances R ss appear to marginally increase from q ⊥ to q ‖ conditions, although the inverse trend is seen for the terminator distances. Again, these differences are slight, which may reflect in part the bias against q ‖ conditions of our bow shock estimator (thus yielding a low amount of q ‖ shock detections).

Figure 10

Comparison of standoff distances, both at the subsolar point and at the terminator, calculated from the 2D and 3D fits, and for each case as in Tables 3 and 4. Terminator standoff values are in blue, orange, and yellow (wider bars), whereas subsolar standoff values are in violet and green (thinner bars). For brevity in the axis labeling, Ls1 = [315°–45°], Ls2 = [45°–135°], Ls3 = [135°–225°], and Ls4 = [225°–315°]. All distances are expressed in units of the planet’s radius, that is, R p = 3,389.5 km.

As a preliminary conclusion, we note that:

The X MSO − Y MSO and X MSO − Z MSO asymmetry seems particularly marked for Ls = [135°–225°] (labeled “Ls3” on the figure), MY32, MY35, and higher EUV fluxes: it can readily be seen by comparing the length of the blue and red bars. As explained earlier, the number of points used for fits for MY32 is the lowest of all the cases because MAVEN arrived at Mars late in MY32. This is in qualitative agreement with the conclusions of Gruesbeck et al. (2018). On average and outside of those special cases, the shock’s shape stays rather symmetric about the X MSO axis: the terminator distances R td, (3D fits) and R td (2D fits) indeed seem to match rather well most of the time. This axisymmetric tendency can be further amplified by aligning the X MSO − Y MSO plane with the solar wind aberration system, rotating the 3D quadric surface 4° anticlockwise around the Z MSO axis; new standoff distances for the 3D fits (R td, and R ss) differ by less than 5% with their corresponding 2D fits values (not shown).

Although the 3D and 2D conic fits retain strong similarities in their behavior, the 3D fits (seemingly paradoxically) appear more robust and less affected by external assumptions. It is recalled here that not only do the 2D fits assume axisymmetry around , but certain 2D fits had to also be constrained at larger Euclidean distances from the center of the planet due to the poor coverage of MAVEN for . This superiority of the 3D fitting algorithm is due to: (a) the number of fitting variables (A to I, allowing more flexibility despite risking over‐determination of the linear system of equations), (b) the natural asymmetry of the shock (however small), and (c) the fitting points being statistically better distributed over a larger space (both in X MSO − Y MSO and X MSO − Z MSO planes instead of a single polar plane) and thus optimizing the fits.

Because the Martian seasons (monitored by Ls ranges) to a degree and the EUV flux both depend on Mars’ heliocentric distance, correlations between these fits are to be expected. For example, similar fits for low Ls values () and low EUV flux can be seen in Figures 8c and 8e (orange curves) and Figure 9b and 9c (red curves).

Because the solar cycle is a continuous underlying driver of the shock’s position regardless of the binnings adopted here (Hall et al., 2019), we expect also correlations between EUV flux and Mars year results. This effect is most clearly exemplified with the shock standoff decrease when going from the declining phase of solar cycle 24 (MY32 and MY33, high fluxes) to the next solar minimum (MY34, low fluxes)

Conclusions

In this study, we presented a fast method to estimate automatically the position of the bow shock in real spacecraft orbits, as well as analytical expressions for the normal direction to the shock surface at any point in its close vicinity. After a survey of existing analytical smooth models of the bow shock surface at the planet Mars based on 2D and 3D fits, we used these models as a first prediction of the shock position in the data and refined this prediction further with a predictor‐corrector algorithm based on the median absolute deviation of the magnetic field around the predicted shock. This method, biased toward the detection of q ⊥ shocks but not entirely limited to them, does not substitute for a detailed analysis of the crossing or for machine‐learning techniques currently developed for space missions. It however finds a useful application when it is necessary to quickly determine the position of the spacecraft, or at least an estimate thereof, with respect to the bow shock.

As part of the solar wind and space weather database Helio4Cast, our technique was successfully used to retrieve solar wind undisturbed parameters from the MAVEN mission (Möstl et al., 2020). We also successfully applied the predictor‐corrector method to the MAVEN orbit and magnetic field data between November 2014 and February 2021 (see list compiled by Simon Wedlund et al., 2021), and performed a series of fits, in 2D and in 3D, to test our method and investigate statistically the shape of the shock depending on MY, solar longitude Ls, and two solar EUV flux levels. The 3D fitting has obvious advantages over the 2D polar axisymmetric geometry usually used to describe the shock structure, namely, a more accurate estimate of asymmetries in the global structure, and taking full advantage of the 3D distribution of bow shock detections in space. This is especially important for bodies such as Mars with large orbital eccentricities and axial tilts to the ecliptic, and for which the heliocentric distance is a strong driver of the EUV flux input and seasonal changes on the planet.

Expectedly, we found the Martian shock to be highly asymmetric with respect to the North‐South hemispheres, in agreement with previous studies (see, e.g., Gruesbeck et al., 2018; Halekas et al., 2017; Hall et al., 2016). Such an asymmetry is in part linked to the presence of crustal magnetic fields at Mars; a specific study taking into account the planet’s rotation and the location of these crustal magnetic sources on the nightside or on the dayside is left for the future. Bow shock fits for quasi‐perpendicular and parallel shock conditions were, to the precision of our approach, almost identical. In addition, the shock appeared noticeably asymmetric with respect to Y MSO and Z MSO directions in specific conditions, namely, for MY32 and MY35, Ls = [135°–225°] and larger EUV fluxes. Despite this observed asymmetry, solar‐aberrated axisymmetric models may still provide a worthy first approximation of the shock’s shape and position.

To investigate further the conditions of the shock’s asymmetry throughout different solar cycles, solar drivers, and internal drivers, such as crustal magnetic fields, and isolate their respective contribution, a full exploitation of MAVEN’s continuously growing datasets is warranted; likewise, a reanalysis of past encounters at Mars using 3D quadric fits would be a welcome addition. These are left for future studies. Applications of these methods, especially in 3D, to other bodies with large orbital eccentricities (such as Mercury) may also prove of interest.

Table A1

Characteristics of the 3D Martian Bow Shock as Derived From the MAVEN Dataset, See Table 4 for the Parameters of the 3D Surfaces Considered

Case	L 1	L 2	L 3	V	V_1 [R p]			V_2 [R p]			V_3 [R p]			P center [R p]			P nose [R p]
		[R p ]		[Rp3]	X MSO	Y MSO	Z MSO	X MSO	Y MSO	Z MSO	X MSO	Y MSO	Z MSO	X MSO	Y MSO	Z MSO	X MSO	Y MSO	Z MSO
All points, Gruesbeck et al. (2018)	8.20	4.34	4.11	613.30	−0.97	0.10	0.22	0.08	−0.72	0.69	0.22	0.69	0.69	−6.45	0.59	0.99	1.52	−0.21	−0.80
All points, this work	2.84	2.74	2.55	83.37	−0.44	−0.06	0.90	−0.16	0.99	−0.01	−0.88	−0.15	−0.45	−1.07	0.12	0.05	1.18	0.50	1.19
MY32, this work	3.42	3.09	2.68	118.38	0.76	−0.28	−0.59	0.13	−0.82	0.55	0.64	0.50	0.59	−1.44	0.17	0.17	1.16	−0.79	−1.83
MY33, this work	2.97	2.82	2.65	92.83	−0.50	0.61	0.62	−0.06	0.69	−0.72	−0.87	−0.40	−0.31	−1.17	0.19	0.09	1.13	1.23	0.90
MY34, this work	2.85	2.70	2.60	83.79	−0.37	0.05	0.93	−0.65	−0.73	−0.21	0.67	−0.68	0.30	−1.15	0.04	0.02	0.60	−1.73	0.80
MY35, this work	2.28	2.13	1.34	27.19	−0.13	0.06	0.99	0.08	1.00	−0.05	−0.99	0.08	−0.13	0.13	0.01	0.00	1.45	−0.09	0.17
Ls = [315°–45°], this work	2.99	2.86	2.73	97.79	0.78	0.62	0.12	−0.25	0.13	0.96	0.58	−0.78	0.25	−1.38	−0.12	−0.01	0.94	1.73	0.34
Ls = [45°–135°], this work	3.08	2.68	2.52	87.10	−0.61	−0.23	0.76	0.21	0.87	0.44	−0.76	0.43	−0.48	−1.23	−0.02	0.20	0.69	−1.10	1.42
Ls = [135°–225°], this work	2.80	2.36	1.88	51.95	0.13	−0.99	0.11	0.74	0.03	−0.67	0.65	0.17	0.74	−0.59	0.35	0.38	1.17	0.41	−1.19
Ls = [225°–315°], this work	3.13	2.86	2.65	99.65	−0.48	−0.33	0.81	−0.47	0.88	0.08	−0.74	−0.34	−0.58	−1.18	0.02	0.24	0.79	0.94	1.76
EUV flux ≥0.002 8 W m−2	3.81	3.13	3.00	150.08	−0.92	0.25	0.32	−0.02	−0.81	0.58	−0.40	−0.53	−0.75	−2.05	0.29	0.22	1.44	−0.67	−0.98
EUV flux <0.0028 W m−2	2.73	2.54	2.27	65.77	−0.26	−0.02	0.97	−0.27	−0.96	−0.09	0.93	−0.28	0.24	−0.81	0.01	0.02	1.30	−0.63	0.57
Quasi‐ ⊥	2.87	2.73	2.53	83.11	−0.43	−0.10	0.90	−0.09	0.99	0.06	−0.90	−0.06	−0.44	−1.06	0.11	0.07	1.21	0.26	1.18
Quasi‐‖	3.00	2.79	2.75	96.12	−0.99	−0.02	0.10	0.09	−0.66	0.75	0.06	0.75	0.66	−1.40	0.03	−0.05	1.58	0.08	−0.36

Note. All quadrics are ellipsoids. L (i = 1, 2, 3) are the lengths in units of R p of the principal axes of the ellipsoids and V their volume. are the eigenvectors of matrix M, that is, the normalized directions of the principal axes in MSO Cartesian coordinates (because the values are normalized to R p and rounded down, a value of 1.00 or 0.00 is not stricto sensu 1 or 0). P center and P nose are the positions of the center of the ellipsoid and its sunward nose, in MSO Cartesian coordinates. The domain of validity for each fit is shown in Figure 9: Fits are valid for X MSO ≥ −0.5 R p on average. Mars' radius is R p = 3,389.5 km.