Magnetosheath plasma flow model around Mercury

The magnetosheath is defined as the plasma region between the bow shock, where the super-magnetosonic solar wind plasma is decelerated and heated, and the outer boundary of the intrinsic planetary magnetic field, the so-called magnetopause. Based on the Soucek–Escoubet magnetosheath flow model at the Earth, we present an analytical magnetosheath plasma flow model around Mercury. The model can be used to estimate the plasma flow magnitude and direction at any given point in the magnetosheath exclusively on the basis of the plasma parameters of the upstream solar wind. The model serves as a useful tool to trace the magnetosheath plasma along the streamline both in a forward sense (away from the shock) and a backward sense (toward the shock), offering the opportunity of studying the growth or damping rate of a particular wave mode or evolution of turbulence energy spectra along the streamline in view of upcoming arrival of BepiColombo at Mercury.

1 Introduction

The magnetosphere of a planet constitutes an obstacle to the super-magnetosonic solar wind. Upstream of the planet a bow shock emerges, because the interplanetary magnetic field (IMF), embedded in the solar wind, cannot simply penetrate the magnetosphere. At the bow shock the super-magnetosonic solar wind plasma is decelerated and heated. The region with the subsonic, heated plasma downstream of the bow shock is called magnetosheath. The magnetosheath plays an important role in the interaction between bow shock and magnetosphere as it conveys energy between the solar wind and the planetary magnetosphere.

One of the earliest magnetosheath plasma flow models is the hydrodynamic model introduced by Spreiter et al. (1966). Basically the model solves the gas-dynamic differential equations of an unmagnetized fluid around an obstacle, represented by the magnetosphere. It has successfully been tested against in situ spacecraft data (Song et al., 1999; Stahara et al., 1993; Spreiter and Alksne, 1968) and applied to model the magnetospheres of various planets in our solar system (see Stahara, 2002, for a review). A decisive drawback of this model, however, is the high complexity and computational demands to calculate numerically a set of differential equations.

To reduce the computational complexity, several analytical plasma flow models have alternatively been proposed (Russell et al., 1983; Kallio and Koskinen, 2000; Romashets et al., 2008). An analytical magnetosheath flow model, which has successfully been tested against spacecraft observations at Earth, has been implemented by Soucek and Escoubet (2012). This model is based on the magnetic field model developed by Kobel and Flückiger (1994) and later modified and extended by Génot et al. (2011) to obtain a magnetosheath plasma flow model. The essential advantage of this model is its compatibility with a wide range of bow shock and magnetopause models while retaining the simplicity and computational efficiency of the original magnetic field model. Furthermore, the model allows us to calculate the plasma flow velocity at any point in the magnetosheath using only the spacecraft position and solar wind parameter upstream of the bow shock.

In this work we follow the procedure proposed by Soucek and Escoubet (2012) and rescale their terrestrial magnetosheath flow model to the space environment at Mercury. First, we introduce the Hermean bow shock and magnetopause model, used to obtain the magnetosheath plasma flow model. Second, we revisit the magnetic field model of Kobel and Flückiger (1994) which Soucek and Escoubet (2012) used to determine the plasma velocity direction in the magnetosheath. Third, we extend the model by the Rankine–Hugoniot relations in a similar way as Génot et al. (2011) to determine the velocity magnitude downstream of the shock.

The aim of this paper is to provide a tool to estimate the plasma flow at a given point of spacecraft observation inside the Hermean magnetosheath on the basis of the solar wind conditions.

2 Bow shock and magnetopause model at Mercury

In the following we use an aberrated Mercury Solar Magnetospheric (MSM) coordinate system. This coordinate system is based on the Mercury Solar Orbital (MSO) coordinate system, but its origin is shifted northward by 479 km from the MSO origin to account for Mercury's dipole offset and rotated into the solar wind velocity direction. In the MSO coordinate system the X_MSO axis points sunward, the Y_MSO axis points antiparallel to Mercury's orbital velocity and $Z_{\mathrm{MSO}}=X_{\mathrm{MSO}}\times Y_{\mathrm{MSO}}$ completes the right-handed system. To compensate for the aberration of the solar wind direction due to the orbital motion of Mercury around the sun, the X_MSO axis is rotated antiparallel to the solar wind flow velocity direction. In the MSM coordinate system the bow shock and magnetopause models are considered to be cylindrically symmetric around the X_MSM axis, reducing the three dimensions $\mathit{\{}{X}_{\mathrm{MSM}},Y_{\mathrm{MSM}},Z_{\mathrm{MSM}}\mathit{\}}$ to two dimensions {X_MSM,ρ_MSM} with $\mathit{\rho }=\sqrt{Y_{\mathrm{MSM}}^{\mathrm{2}}+Z_{\mathrm{MSM}}^{\mathrm{2}}}$.

Slavin et al. (2009) modeled the bow shock at Mercury by a conic section of the form

$$\begin{array}{}\text{(1)}& \mathit{\xi }=\sqrt{(x_{\text{BS}}-x_{\mathrm{0}}{)}^{\mathrm{2}}+{\mathit{\rho }}_{\text{BS}}^{\mathrm{2}}}={\displaystyle \frac{p\mathit{ϵ}}{\mathrm{1}+\mathit{ϵ}\cos \mathit{\varphi }}},\end{array}$$

with x₀ being the distance of the focus of the conic section from the dipole center along X_MSM, ${\mathit{\rho }}_{\text{BS}}=\sqrt{y_{\text{BS}}^{\mathrm{2}}+z_{\text{BS}}^{\mathrm{2}}}$ being the distance from the X_MSM axis, p being the focal parameter and ϵ being the eccentricity. With the advent of the MESSENGER (MErcury Surface, Space ENvironment, GEophysics and Ranging;  Solomon et al., 2007) spacecraft in an orbit around Mercury, it was possible to characterize the spatial location of the bow shock and magnetopause statistically. Winslow et al. (2013) determined that the best-fit parameters to the bow shock are given by x₀=0.5 R_M, p=2.75 R_M and ϵ=1.04. With these parameters the extrapolated subsolar bow shock stand-off distance is R_BS=1.9 R_M (Mercury radii, 1 R_M∼2440 km). For this work, it is advantageous to transform Eq. (1) into the origin of the MSM coordinate system with

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{r}_{\mathrm{BS}}& =\sqrt{(\mathit{\xi }\cos \mathit{\varphi }+x_{\mathrm{0}}{)}^{\mathrm{2}}+(\mathit{\xi }\sin \mathit{\varphi }{)}^{\mathrm{2}}},\\ \mathit{\theta }& =\arccos \left({\displaystyle \frac{\mathit{\xi }\cos \mathit{\varphi }+x_{\mathrm{0}}}{r_{\mathrm{BS}}}}\right),\end{array}\end{array}$$

where r_BS is the distance from the dipole center to the bow shock, and θ is the angle between r_BS and the X_MSM axis. Figure 1 shows a schematic illustration of the parameters ξ, ϕ, r_BS and θ, which are used in the formulation of the bow shock.

Korth et al. (2015) used the magnetopause model proposed from Shue et al. (1997) and found that the MESSENGER observations of magnetopause crossing are best fit by

$$\begin{array}{}\text{(3)}& r_{\mathrm{MP}}=\sqrt{x_{\mathrm{MP}}^{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{MP}}^{\mathrm{2}}}=R_{\mathrm{MP}}{\left({\displaystyle \frac{\mathrm{2}}{\mathrm{1}+\cos \mathit{\theta }}}\right)}^{\mathit{\alpha }},\end{array}$$

with α=0.5 being the best-fit flaring parameter, and R_MP=1.42 R_M being the subsolar stand-off magnetopause distance.

Figure 1 shows the Slavin et al. (2009) bow shock model (S09-BS) and Korth et al. (2015) magnetopause model (K15-MP) evaluated from Eqs. (2) and (3), respectively.

Figure 1 Schematic representation of the parameters used in the formulation of the bow shock and magnetopause in the MSM equatorial plane. The solid red line is the bow shock evaluated from Eq. (2) (S09-BS;  Slavin et al., 2009). The solid black line represents the magnetopause determined by Eq. (3) (K15-MP;  Korth et al., 2015). The dashed blue lines are the bow shock and magnetopause determined by Eq. (4) from the KF94 model (Kobel and Flückiger, 1994).

3 The KF94 magnetic field model

To obtain the magnetosheath plasma flow direction, we follow the procedure proposed by Soucek and Escoubet (2012) and use the magnetic field model developed by Kobel and Flückiger (1994). In the following we denote this model as KF94 model and mark all quantities pertaining to the KF94 model by a tilde, e.g., $\stackrel{\mathrm{̃}}{r}$. In the KF94 model the bow shock (BS) and magnetopause (MP) at Mercury are modeled by parabolic surfaces at a common focus with

$$\begin{array}{}\text{(4)}& {\stackrel{\mathrm{̃}}{r}}_{\{\mathrm{BS},\mathrm{MP}\}}={\displaystyle \frac{-\cos \mathit{\theta }+\sqrt{{\cos }^{\mathrm{2}}\mathit{\theta }+\mathrm{4}{R}_{\{\mathrm{BS},\mathrm{MP}\}}{b}_{\{\mathrm{BS},\mathrm{MP}\}}{\sin }^{\mathrm{2}}\mathit{\theta }}}{\mathrm{2}{b}_{\{\mathrm{BS},\mathrm{MP}\}}{\sin }^{\mathrm{2}}\mathit{\theta }}},\end{array}$$

with $b_{\mathrm{BS}}=\mathrm{1}/(\mathrm{4}{R}_{\mathrm{BS}}-\mathrm{2}{R}_{\mathrm{MP}})$ and $b_{\mathrm{MP}}=\mathrm{1}/\left(\mathrm{2}{R}_{\mathrm{MP}}\right)$ defined by the stand-off distances R_{BS,MP}.

Under the assumption that the IMF is parallel to the solar wind, the magnetic field lines of the KF94 model represent the flow lines of the solar wind and magnetosheath plasma. Using the magnetic field vector direction in the KF94 model, Soucek and Escoubet (2012) determined the flow velocity vector at a given position $\mathit{r}=(x,\mathit{\rho })$ by

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{v}}_{\mathrm{x}}=v_{\mathrm{m}}(C/\mathrm{2}d-C/R_{\mathrm{MP}}),\\ & {\stackrel{\mathrm{̃}}{v}}_{\mathit{\rho }}=v_{\mathrm{m}}(C\mathit{\rho }/[\mathrm{2}d(d+x-R_{\mathrm{MP}}/\mathrm{2})\left]\right),\end{array}\end{array}$$

where v_m corresponds to the flow velocity magnitude, $d=|\mathit{r}-{\mathit{r}}_{\mathrm{0}}|$ is the difference between the given position in the magnetosheath and the parabolic surface focus ${\mathit{r}}_{\mathrm{0}}=(R_{\mathrm{MP}}/\mathrm{2},\mathrm{0})$ and $C=R_{\mathrm{MP}}(\mathrm{2}{R}_{\mathrm{BS}}-R_{\mathrm{MP}})/(\mathrm{2}{R}_{\mathrm{BS}}-\mathrm{2}{R}_{\mathrm{MP}})$ is a constant defined by the bow shock and magnetopause stand-off distances.

4 The magnetosheath plasma flow model around Mercury

To obtain the magnetosheath plasma flow at a specific point ($\mathit{r}=(x,\mathit{\rho })$) in the magnetosheath at Mercury, we evaluate the plasma flow direction first and then determine the magnitude of the velocity vector from the Rankine–Hugoniot relation across the bow shock.

A magnetic field model is used to describe the plasma flow here. The reason for this is explained as follows. Assumptions are made such that there are no proton sinks or sources in the magnetosheath. Strictly speaking, this assumption is only weakly justified because of the neutral particles such as the hydrogen corona and the sodium exosphere around Mercury. A stationary case is taken in the fluid picture. The continuity equation reads then as

$$\begin{array}{}\text{(6)}& \mathrm{\nabla }\cdot \left(n\mathit{U}\right)=\mathrm{0},\end{array}$$

where n and U are the density and velocity of protons, respectively. Now we make an analogy such that

$$\begin{array}{}\text{(7)}& \left(n\mathit{U}\right)⟶\mathit{B},\end{array}$$

holds, and Eq. (6) is equivalent to the divergence-free condition of magnetic field:

$$\begin{array}{}\text{(8)}& \mathrm{\nabla }\cdot \mathit{B}=\mathrm{0}.\end{array}$$

For a more detailed discussion, see e.g., Génot et al. (2011).

4.1 Plasma flow direction

Following the procedure proposed by Soucek and Escoubet (2012), we rescale the plasma flow direction from the KF94 model to Mercury's space environment as follows:

1. As a first step, we calculate the angle $\mathit{\theta }=\arccos (x/\sqrt{x^{\mathrm{2}}+{\mathit{\rho }}^{\mathrm{2}}})$ between r and the X_MSM axis.

2. Then we estimate the fractional distance, ℱ, of r between the bow shock and magnetopause from Eqs. (2) and (3) with

   $$\begin{array}{}\text{(9)}& \mathcal{F}={\displaystyle \frac{r\left(\mathit{\theta }\right)-r_{\mathrm{BS}}\left(\mathit{\theta }\right)}{r_{\mathrm{BS}}\left(\mathit{\theta }\right)-r_{\mathrm{MP}}\left(\mathit{\theta }\right)}}.\end{array}$$

3. Now we change into the KF94 model and calculate ${\stackrel{\mathrm{̃}}{r}}_{\mathrm{BS}}\left(\mathit{\theta }\right)$ and ${\stackrel{\mathrm{̃}}{r}}_{\mathrm{MP}}\left(\mathit{\theta }\right)$ from Eq. (4) with the angle θ. Note that the stand-off distances (R_BS and R_MP) and the focus (${\mathit{r}}_{\mathrm{0}}=(R_{\mathrm{MP}}/\mathrm{2},\mathrm{0})$) in Eq. (4) are the best-fit values from Eqs. (2) and (3).

4. In a next step we determine the fractional position within the magnetosheath in the KF94 model with

   $$\begin{array}{}\text{(10)}& \stackrel{\mathrm{̃}}{r}\left(\mathit{\theta }\right)=\mathcal{F}\left[{\stackrel{\mathrm{̃}}{r}}_{\mathrm{BS}}\right(\mathit{\theta })-{\stackrel{\mathrm{̃}}{r}}_{\mathrm{MP}}(\mathit{\theta }\left)\right]+{\stackrel{\mathrm{̃}}{r}}_{\mathrm{BS}}\left(\mathit{\theta }\right),\end{array}$$

   according to Eq. (9).

5. Using Eq. (5) we evaluate the KF94 flow velocity vector, $\stackrel{\mathrm{̃}}{\mathit{v}}$ for $\stackrel{\mathrm{̃}}{\mathit{v}}=({\stackrel{\mathrm{̃}}{v}}_x,{\stackrel{\mathrm{̃}}{v}}_{\mathit{\rho }})$ at the position $\stackrel{\mathrm{̃}}{r}\left(\mathit{\theta }\right)$. Note that the velocity magnitude v_m is determined in a later step.

6. With the obtained flow velocity vector, v, we are able to estimate the new position of an adjacent point along the same flow line ${\stackrel{\mathrm{̃}}{r}}^{\prime }=\stackrel{\mathrm{̃}}{r}+\stackrel{\mathrm{̃}}{v}\mathrm{\Delta }t$, by choosing an infinitesimally small time increment Δt.

7. Next we determine the angle between the new position and the X_MSM axis, θ^′, and the fractional distance inside the KF94 magnetosheath, ℱ^′, using Eq. (9).

8. Finally we transform the new position ${\stackrel{\mathrm{̃}}{r}}^{\prime }\left({\mathit{\theta }}^{\prime }\right)$ back from the KF94 model to the original MSM reference frame where the magnetosheath is confined by Eqs. (2) and (3). The new position, r^′, inside this magnetosheath is then given by

   $$\begin{array}{}\text{(11)}& {\mathit{r}}^{\prime }={\mathcal{F}}^{\prime }\left[{\mathit{r}}_{\mathrm{BS}}\right({\mathit{\theta }}^{\prime })-{\mathit{r}}_{\mathrm{MP}}({\mathit{\theta }}^{\prime }\left)\right]+{\mathit{r}}_{\mathrm{BS}}\left({\mathit{\theta }}^{\prime }\right),\end{array}$$

   and thus the plasma flow direction can be determined by $\mathit{v}=({\mathit{r}}^{\prime }-\mathit{r})/\mathrm{\Delta }t$.

Applying recursively this procedure (steps 1–8) yields the plasma flow line within Mercury's magnetosheath. In Fig. 2 five examples of flow lines are shown.

Figure 2 Schematic representation of the flow lines (red) in the MSM equatorial plane. The solid black lines are the bow shock and magnetopause evaluated from Eqs. (2) and (3), respectively.

4.2 Plasma flow magnitude

To evaluate the magnetosheath plasma velocity magnitude, v_m, we apply the Rankine–Hugoniot (RH) equations, which relate the upstream (u) with the downstream (d) plasma conditions. The downstream plasma flow velocity directly behind the bow shock, v^d, is derived by the following procedure:

1. From the given spacecraft position in the magnetosheath, $\mathit{r}=(x,\mathit{\rho })$, we trace the flow line back to the bow shock. Thereto we iteratively apply steps (1)–(8) from above, with reversed increments ${\stackrel{\mathrm{̃}}{r}}^{\prime }=\stackrel{\mathrm{̃}}{r}-\stackrel{\mathrm{̃}}{v}\mathrm{\Delta }t$ in step (6), until the bow shock is reached (${\mathcal{F}}^{\prime }=\mathrm{0}$). Then we calculate the angle θ_BS between the X_MSM axis and the bow shock intersection at (x_BS,ρ_BS) with ${\mathit{\theta }}_{\mathrm{BS}}=\arccos (x_{\mathrm{BS}}/\sqrt{x_{\mathrm{BS}}^{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{BS}}^{\mathrm{2}}})$.

2. In a next step we determine the bow shock tangent $\widehat{\mathit{t}}$ and normal $\widehat{\mathit{n}}$ unit vectors where the back-traced flow line intersects the bow shock. For any point along the bow shock, the normal, n, and tangent, t, vector can easily be computed by

   $$\begin{array}{}\text{(12)}& \begin{array}{rl}& \mathit{t}=\left[{\displaystyle \frac{\mathrm{d}{r}_{\mathrm{BS}}}{\mathrm{d}\mathit{\theta }}}\cos \mathit{\theta }-r_{\mathrm{BS}}\sin \mathit{\theta }\right]{\widehat{e}}_x+\left[{\displaystyle \frac{\mathrm{d}{r}_{\mathrm{BS}}}{\mathrm{d}\mathit{\theta }}}\sin \mathit{\theta }+r_{\mathrm{BS}}\cos \mathit{\theta }\right]{\widehat{e}}_{\mathit{\rho }},\\ & \mathit{n}=\left[{\displaystyle \frac{\mathrm{d}{r}_{\mathrm{BS}}}{\mathrm{d}\mathit{\theta }}}\sin \mathit{\theta }+r_{\mathrm{BS}}\cos \mathit{\theta }\right]{\widehat{e}}_x-\left[{\displaystyle \frac{\mathrm{d}{r}_{\mathrm{BS}}}{\mathrm{d}\mathit{\theta }}}\cos \mathit{\theta }+r_{\mathrm{BS}}\sin \mathit{\theta }\right]{\widehat{e}}_{\mathit{\rho }},\end{array}\end{array}$$

   where $\frac{\mathrm{d}{r}_{\mathrm{BS}}}{\mathrm{d}\mathit{\theta }}$ is numerically calculated from two consecutive points along the bow shock given by Eq. (2). The shock frame of reference is then defined by the normalized normal and tangent vector at θ_BS.

3. In the shock reference frame the RH equations can be combined to determine the downstream velocity vector component parallel ${\mathit{v}}_{\text{n}}^{\text{d}}$ and perpendicular ${\mathit{v}}_{\text{t}}^{\text{d}}$ to the shock normal (see e.g., Génot, 2008):

   $$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\mathit{v}}_{\text{n}}^{\text{d}}={\mathit{v}}_{\text{n}}^{\text{u}}{\displaystyle \frac{\mathrm{1}}{\mathcal{R}}},\\ & {\mathit{v}}_{\text{t}}^{\text{d}}={\mathit{v}}_{\text{n}}^{\text{u}}\left(\tan {\mathit{\theta }}_{\mathrm{Vn}}+{\displaystyle \frac{\mathrm{1}-\mathrm{1}/\mathcal{R}}{\left[{M_{\mathrm{A}}^{\text{u}}}^{\mathrm{2}}/\left(\mathcal{R}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{Bn}}\right)\right]-\mathrm{1}}}\tan {\mathit{\theta }}_{\mathrm{Bn}}\right),\end{array}\end{array}$$

   where ${\mathit{v}}_{\text{n}}^{\text{u}}$ is the upstream velocity vector component parallel to the shock normal, $\mathcal{R}={\mathit{\rho }}^{\text{d}}/{\mathit{\rho }}^{\text{u}}$ is the compression ratio between the upstream and downstream mass density, ${\mathit{\theta }}_{\mathrm{Vn}}=\arctan (v_{\text{t}}^{\text{u}}/v_{\text{n}}^{\text{u}})$ is the angle between the upstream velocity vector and the shock normal, ${\mathit{\theta }}_{\mathrm{Bn}}=\arctan (B_{\text{t}}^{\text{u}}/B_{\text{n}}^{\text{u}})$ is the angle between upstream magnetic field vector and shock normal, and $M_{\mathrm{A}}^{\text{u}}=v_{\text{n}}^{\text{u}}\frac{\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{\rho }}^{\text{u}}}}{B_{\text{n}}^{\text{u}}}$ is the upstream Alfvén Mach number.

4. In Eq. (13) all parameters pertain to the upstream side, except the compression ratio ℛ. However, ℛ can also be expressed by exclusively upstream parameters with (see e.g., Anderson, 1963)

   $$\begin{array}{}\text{(14)}& \begin{array}{rl}& ({M_A^{\text{u}}}^{\mathrm{2}}-\mathcal{R}{)}^{\mathrm{2}}(\mathit{\gamma }{\mathit{\beta }}^{\text{u}}\mathcal{R}+{M_{\mathrm{A}}^{\text{u}}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{Bn}}\left[\right(\mathit{\gamma }-\mathrm{1})\\ & \mathcal{R}-(\mathit{\gamma }+\mathrm{1})\left]\right)+\mathcal{R}{M_{\mathrm{A}}^{\text{u}}}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{Bn}}\left(\right[\mathit{\gamma }+(\mathrm{2}-\mathit{\gamma })\mathcal{R}]\\ & {M_{\mathrm{A}}^{\text{u}}}^{\mathrm{2}}+\mathcal{R}\left[\right(\mathit{\gamma }-\mathrm{1})\mathcal{R}-(\mathit{\gamma }+\mathrm{1}\left)\right])=\mathrm{0}\end{array}\end{array}$$

   where ${\mathit{\beta }}^{\text{u}}=\left(\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}{p}^{\text{u}}\right)/{B^{\text{u}}}^{\mathrm{2}}$ is the ratio of the upstream thermal to magnetic pressures, and γ is the polytropic index which is typically assumed to be $\mathit{\gamma }=\mathrm{5}/\mathrm{3}$. Equation (14) is equivalent to Eq. (2.43) in Anderson (1963) with different notations. The solution is given in the form of compression ratio as a function of the shock angle θ_Bn. Solutions exist for a compression ratio in the range $\mathrm{1}\le \mathcal{R}\le \mathrm{4}$ when using Eq. (14). Another class of solutions also exists for the expansion (ℛ<1) with a decrease of entropy from the upstream onto the downstream side. The latter case is not physically relevant and is not considered here. The upper limit of compression ratio (ℛ=4) corresponds to the limit of high Alfvén Mach number (M_A→∞) under a polytropic index of $\mathit{\gamma }=\mathrm{5}/\mathrm{3}$.

5. By solving Eq. (14) for ℛ the downstream velocity magnitude $v_{\mathrm{d}}=\sqrt{{v_{\text{n}}^{\text{d}}}^{\mathrm{2}}+{v_{\text{t}}^{\text{d}}}^{\mathrm{2}}}$ is therefore entirely determined by only the upstream plasma parameters.

Since the detailed density profile along the flow line is unknown, we assume in a first approximation a constant plasma density and thus a constant velocity magnitude along the flow line. Therefore, the velocity magnitude at a given point r directly corresponds to the velocity magnitude downstream of the shock, and v_m=v_d. Although this assumption has the tendency to underestimate the velocity close to the magnetopause, it yields satisfactory results in a first approach (Génot et al., 2011).

The entire procedure from above is implemented in an IDL computer program which can be retrieved from OSF (Schmid, 2020). The program is designed to evaluate the plasma flow velocity vector at a given observation point of a spacecraft inside the Hermean magnetosheath exclusively on the basis of the upstream solar wind conditions. As the solar wind input parameters, we use the solar wind propagation model of Tao et al. (2005), which is modified by the orbital motion of Mercury. The model is a one-dimensional magnetohydrodynamic model and takes the OMNI dataset as input to compute propagation at all solar system bodies including Mercury. The correction for the orbital motion is achieved by adding the solar wind velocity vector, V_SW, (which is radially away from the Sun) and the orbital motion velocity vector of Mercury, V_mercury, with $\mathit{V}={\mathit{V}}_{\mathrm{SW}}+{\mathit{V}}_{\text{mercury}}$. To obtain V_mercury, we use the dataset provided from the Navigation and Ancillary Information Facility (NAIF; Acton, 1996), which provide the distance between Mercury and Sun, D, and absolute velocity of Mercury, V_mercury. To determine aberration velocity vector, we first calculate the aberration angle ϕ on the basis of Mercury's elliptical orbit with

$$\begin{array}{}\text{(15)}& \mathit{\varphi }=\arctan \left(b/\sqrt{a^{\mathrm{2}}-b^{\mathrm{2}}}\sin \left[\text{arccos}\left(\right(\mathrm{1}-p/D)/e)\right]\right),\end{array}$$

where a is the semimajor axis, b the semiminor axis, p the semi-latus rectum and e the eccentricity of the ellipse. With the aberration angle ϕ and the consideration whether Mercury moves towards or away from the sun, we subsequently obtain the aberration velocity vector ${\mathit{V}}_{\text{mercury,x}}=\pm V_{\text{mercury}}\sin \left(\mathit{\varphi }\right)$ and ${\mathit{V}}_{\text{mercury,y}}=\pm V_{\text{mercury}}\cos \left(\mathit{\varphi }\right)$. The transformation due to this abberation effect is made by applying a two-dimensional rotation matrix to the spatial coordinates (spanning the x and ρ coordinates):

$$\begin{array}{}\text{(16)}& \left(\begin{array}{l}{x}^{\prime }\\ y^{\prime }\end{array}\right)=\left(\begin{array}{ll}\cos {\mathit{\theta }}_{\mathrm{a}}& -\sin {\mathit{\theta }}_{\mathrm{a}}\\ \sin {\mathit{\theta }}_{\mathrm{a}}& \cos {\mathit{\theta }}_{\mathrm{a}}\end{array}\right)\left(\begin{array}{l}x\\ \mathit{\rho }\end{array}\right),\end{array}$$

where the aberration angle θ_a is given by the radial solar wind velocity, V_SW, and the apparent solar wind velocity, V with ${\mathit{\theta }}_{\mathrm{a}}=\arccos \left({\mathit{V}}_{\mathrm{SW}}\cdot \mathit{V}\right)$. In Fig. 3 the results of the model are shown for the average solar wind plasma parameters during the entire MESSENGER operation service between 2011 and 2015. After modifying the solar wind velocity vector of the Tao et al. (2005) model by the orbital motion of Mercury (Acton, 1996), the average input solar wind plasma parameters for our model are density of $n^{\text{u}}\approx \mathrm{40}{\mathrm{cm}}^{-\mathrm{3}}$, temperature of T≈18 eV, flow speed of $\left|{\mathit{V}}^{\text{u}}\right|\approx -\mathrm{400}\mathrm{km}/\mathrm{s}$, magnetic field magnitude of $\left|{\mathit{B}}^{\text{u}}\right|\approx \mathrm{20}\mathrm{nT}$ with the radial component B_r≈18 nT (ignoring the sign) and the tangential component B_t≈16 nT (ignoring the sign). The mean values of density, temperature, flow speed and magnetic field are valid for nearly 1500 d of observations of MESSENGER (confirmed by one of the reviewers). It is worthwhile to note that one finds an angle of 23^∘ from the Tao dataset, which is consistent to the spiral angle at Mercury at an average position of 0.4 astronomical units (AU) from the Sun for a solar wind speed of $\mathrm{400}\mathrm{km}/\mathrm{s}$. The B_x component is computed from the B_y component using the Parker spiral field in the Tao model. The Alfvén Mach number in the upstream region is $M_{\mathrm{A}}=v^{\text{u}}/V_{\mathrm{A}}\approx \mathrm{5.8}$ (with an Alfvén speed of $V_{\mathrm{A}}\approx \mathrm{69}\mathrm{km}/\mathrm{s}$), and the plasma parameter beta (upstream) is $\mathit{\beta }=\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}n{k}_{\mathrm{B}}T/B^{\mathrm{2}}\approx \mathrm{0.72}$ (where μ₀ is the permeability of free space, and k_B is the Boltzmann constant) in our setup. Color coded is the obtained velocity magnitude v_m. Additionally, the back-traced flow line from a virtual spacecraft located at $x_{\mathrm{MSM}}=-\mathrm{3}{\mathrm{R}}_{\mathrm{M}}$ and ρ_MSM=3 R_M is plotted (green line). At the bow shock intersection the calculated shock normal $\widehat{\mathit{n}}$ and tangent $\widehat{\mathit{t}}$ are illustrated in thin black lines. The thick black line downstream of the bow shock shows the velocity vector direction evaluated from the RH relations, which is in good agreement with the streamline direction determined by the KF94 model.

Figure 3 Color-coded flow speed calculated from Eq. (13) with the averaged upstream parameters from the Tao et al. (2005) solar wind propagation model at Mercury between 2011 and 2015. Also plotted is a schematic representation of back-traced flow line (green) from a virtual spacecraft (S/C) to the bow shock (BS) with the respective shock normal $\widehat{\mathit{n}}$ and tangent $\widehat{\mathit{t}}$ obtained from Eq. (12). The thick black line downstream of the BS is the downstream velocity vector direction determined by Eq. (13). Alfvén Mach number is 5.8, and plasma beta is 0.79.

At the virtual spacecraft position the model predicts a magnetosheath plasma flow velocity of $v_{\mathrm{x}}\approx -\mathrm{200}\mathrm{km}/\mathrm{s}$ and $v_{\mathit{\rho }}\approx \mathrm{17}\mathrm{km}/\mathrm{s}$.

A naive picture of the compression by a factor of 4 (in the limit of high Mach number) is not realistic, because the interplanetary magnetic field around Mercury reaches a magnitude of 20 to 50 nT, and the Alfvén Mach number is correspondingly smaller than that around the Earth by 20 to 50, respectively. A picture of the constant density and a reduced flow speed to $\mathrm{1}/\mathrm{4}$ of the solar wind speed at the magnetopause (along the streamline tangential to the magnetopause) is not valid, either, since the adiabatic expansion breaks down, and the model is not applicable to the flow in the subsolar region and at the magnetopause.

5 Discussion and conclusions

Here we present the first analytical magnetosheath plasma flow model for the space environment around Mercury. The model is based on the magnetosheath model by Soucek and Escoubet (2012), which has successfully been tested against spacecraft observations at Earth. The proposed model is relatively simple to implement and provides the possibility to trace the flow lines inside the Hermean magnetosheath.

The model presented in this paper is generally robust and easy to implement for its analytic expression using upstream parameters. It can help to determine the local plasma conditions of a spacecraft in the magnetosheath exclusively on the basis of the upstream solar wind parameters. Two applications are in mind in view of the BepiColombo mission, where the Mercury Magnetospheric Orbiter (MMO also referred to as Mio) will probe Mercury's magnetosheath and solar wind with unprecedented fast measurements of the particle distribution functions: (1) the Tátrallyay method to observationally determine the growth rate or damping rate of specific mode such as the mirror mode along the streamline (Tátrallyay and Erdös, 2002; Tátrallyay et al., 2008) and (2) the Guicking method to observationally track the spectral evolution of turbulent fluctuations along the streamline (Guicking et al., 2010, 2012).

At the moment a comparison with plasma data is technically not feasible for our model. Due to the limited particle measurements on board MESSENGER, it is not possible to obtain the plasma parameters properly in the solar wind and magnetosheath, that is, by covering the full velocity distributions and to compare with the model velocities. Above all, the plasma instrument is located behind the heat shield and has just a limited field of view. Due to this fact, the majority (thermal core part) of solar wind particles cannot be detected. Comparison with the numerical simulations would be another possibility to test for the model, but a quantitative comparison remains a challenge for the reason that there are large discrepancies in the density and flow velocity among various simulation models on the dayside from the subsolar point to the northern terminator (Aizawa et al., 2021).

An approximation that the magnetic field is more aligned with the solar wind flow direction is more justified at Mercury than at the Earth because of the Parker spiral nature. One of the possible consequences of our assumption is that the magnetic field magnitude would change or evolve in the same sense as that of the plasma (or particle number flux) in the magnetosheath. The correlation between the magnetic field and that of the particle flux in the magnetosheath would ideally be tested against the plasma and magnetic field data on the arrival of BepiColombo at Mercury. The change in the number density can be interpreted as the change in the cross-sectional area of a flux tube across which the plasma streams. The change in the flow velocity can then be compared with that from the adiabatic expansion and that from the measurement.

The following remarks are drawn as scientific message. First, Fig. 1 visually demonstrates that different models predict different shapes of the tail and magnetosheath, which is an overlooked issue in the Mercury magnetosphere community. Second, Fig. 2 shows that the flow lines near the subsolar region (Sun-to-planet line if neglecting the planetary orbital motion) expand abruptly so that the adiabatic expansion may break down. In particular, the adiabatic expansion plays an important role in predicting the flow in the magnetosheath. Logical continuation of our model construction would be to evaluate also the density and velocity profile along the flow line and to test under which conditions the adiabatic expansion breaks down.

Although the proposed model has a good performance overall for a wide range of upstream conditions, the accuracy strongly depends on the used bow shock and magnetopause model. Here we utilize the bow shock and magnetopause model from Slavin et al. (2009) and Korth et al. (2015), which were adopted from MESSENGER boundary crossing observations.

The presented model is cylindrically symmetric around the X_MSM axis. In reality, however, non-radial IMF conditions will lead to a spatially asymmetric magnetosheath (Nishino et al., 2008; Dimmock and Nykyri, 2013; Dimmock et al., 2016). On the quasi-perpendicular side, where the shock-normal angles θ_Bn are greater than 45^∘, the magnetosheath is known to be thicker with larger plasma flow velocities than on the quasi-parallel side, where ${\mathit{\theta }}_{\mathrm{Bn}}<\mathrm{45}{}^{\circ }$. Such asymmetries cannot be reproduced by the simple model presented here but should be addressed in future work.

Furthermore, the method used to determine the flow velocity magnitude can possibly be improved. Here we assumed a constant plasma density and velocity along the flow line which has the tendency to underestimate the plasma velocity in regions with lower densities, e.g., close to the magnetopause. Génot et al. (2011) proposed a simple ad hoc model of a plasma density profile which has been implemented by Soucek and Escoubet (2012). While this ad hoc density model showed good correspondence with in situ spacecraft plasma observation at Earth, the solar wind and magnetospheric conditions at other planets can be very different (like at Mercury) and thus might give a worse prediction (Soucek and Escoubet, 2012). Our model inherits the properties from the Soucek–Escoubet model by scaling the Kobel–Flückinger model of the near-Earth environment: (1) time stationary flow and (2) axially symmetric around the axis of (apparent) solar wind penetrating the planet. Assumption of time stationary flow may break down when the change in the solar wind state is not negligible. Assumption of the axisymmetric magnetosheath may also break down when the magnetopause location is not symmetric between the northern and the southern hemisphere (in particular, in the tail region).

At this stage we decided not to include an ad hoc density profile, also because it can hardly be tested due to the limited plasma observations around Mercury. The assumption of constant density implies a constant velocity along a given flow line. The velocity profile may vary considerably from that estimated in the earlier models such as the Spreiter model (Spreiter et al., 1966), the Genot model (Génot et al., 2011), and the Soucek–Escoubet model (Soucek and Escoubet, 2012). As mentioned above, the constant density will likely underestimate the propagation timing in the magnetosheath. In future work such a density profile should be evaluated and included.