Scalar-potential mapping of the steady-state magnetosheath model

The steady-state magnetosheath model has various applications for studying the plasma and magnetic field profile around the planetary magnetospheres. In particular, the magnetosheath model is analytically obtained by solving the Laplace equation for parabolic boundaries (bow shock and magnetopause). We address the question, how can we utilize the magnetosheath model by transforming into a more general, empirical, non-parabolic magnetosheath geometry? To achieve the goal, we develop the scalar-potential mapping method which provides a semi-analytic estimate of steady-state flow velocity and magnetic field in the empirical magnetosheath domain. The method makes use of a coordinate transformation from the empirical magnetosheath domain into the parabolic magnetosheath domain and evaluates a set of variables (shell variable and connector variable) to utilize the solutions of the Laplace equation obtained for the parabolic magnetosheath domain. Our model uses two invariants of transformation: the zenith angle within the magnetosheath with respect to the direction to the Sun and the ratio of the distance to the magnetopause and the thickness of magnetosheath along the magnetopause-normal direction. The use of magnetopause-normal direction makes a marked difference from the earlier model construction using the radial direction as reference. The plasma flow and magnetic field can be determined as a function of the upstream condition (flow velocity or magnetic field) in a wide range of zenith angles. The scalar-potential mapping method is computationally inexpensive, using analytic expressions as much as possible, and is applicable to various planetary magnetosheath domains.

1 Introduction

Steady-state plasma flow and magnetic field can be regarded as a realization of potential field in the planetary magnetosheath region when the vorticity and the electric current are ignored. In such a case, the potential is obtained by solving the Laplace equation, which was elegantly and analytically solved by Kobel and Flückiger (1994) for a parabolic shape of magnetosheath (hereafter KF). The KF potential was further extended to the stream function in the magnetosheath by Guicking et al. (2012). The KF solution made a series of breakthroughs in magnetosheath research. One of the most successful applications is the ability to track the plasma parcel along the streamline in the modeled magnetosheath. The tracking method was extensively used to observationally study the mirror mode growth (e.g., Tatrallyay et al., 2002; Génot et al., 2011) and the streamwise turbulence evolution in the magnetosheath (Guicking et al., 2012). Predictive models of plasma flow and magnetic field serve as a useful tool when combined with the numerical simulation or the observational data.

The KF potential is obtained using the assumption that the planetary bow shock and magnetopause have a parabolic shape sharing the same focal point. Empirical models of the bow shock and magnetopause (fitted to the spacecraft data), on the other hand, are not necessarily parabolically or co-focally shaped. For example, the empirical Earth bow shock model by Farris et al. (1991) and Cairns et al. (1995) has a parabolic shape, but the focal point differs from that of the KF solution. The empirical magnetopause model by Shue et al. (1997) applies a power-law scaling to the parabolic shape such that the magnetic field lines appear stretched in the tail region. The gap between the KF parabolic magnetosheath and the empirical magnetosheath needs to be filled when applying the KF potential in the empirical magnetosheath.

Naively speaking, one wishes to find a conformal mapping (angle-preserving mapping) from the KF parabolic magnetosheath onto a non-parabolic empirical magnetosheath shape such as the analytic extension of magnetopause shape (Narita et al., 2023). However, no general mathematical algorithm is known so far to obtain the conformal mapping when the spatial domain is not properly bounded. The problem lies in the fact that the magnetosheath is bounded only by two sides, i.e., the standing shock and the magnetopause in the radial direction to the planet, but it is not bounded along the flow in the tail region. The algorithms of numerical conformal mapping are so far proposed for spatially bounded domains (Papamichael and Whiteman, 1973; Chakravarthy and Anderson, 1979; Fornberg, 1980; Karageorghis et al., 1996) or domains with a closed shape of internal boundaries (Wei et al., 2014).

Here we address the question, how can we utilize the KF magnetosheath model by transforming into a more general, empirical, non-parabolic magnetosheath geometry? To achieve the goal, we develop a mapping method which provides a semi-analytic estimate of steady-state flow velocity and magnetic field in the empirical magnetosheath domain. Our scalar-potential mapping method is computationally inexpensive by using the analytic expression as much as possible, and it is applicable to various planetary magnetosheath domains.

This work is organized in the following fashion. After reviewing the magnetosheath model constructed by Kobel and Flückiger (1994) (Sect. 2) and discussing different mapping methods (Sect. 3), the detailed procedure of the magnetopause-normal mapping is presented (Sect. 4) and discussed (Sect. 5), with concluding remarks (Sect. 6).

2 Revisiting the magnetosheath scalar potential

2.1 Parabolic coordinates

In the KF parabolic coordinates, the shell variable v (iso-contour lines enveloping the magnetosphere) and the connector u (iso-contour lines connecting from the bow shock to the magnetopause) play an important role in computing the flow velocity and magnetic field in the magnetosheath. These variables are explicitly evaluated using Cartesian coordinates and the radial distance from the focal point as

$$\begin{array}{}\text{(1)}& {\displaystyle }v& {\displaystyle }=\sqrt{r_{\mathrm{0}}+(x_{\mathrm{k}}-x_{\mathrm{0}})},\text{(2)}& {\displaystyle }u& {\displaystyle }=\sqrt{r_{\mathrm{0}}-(x_{\mathrm{k}}-x_{\mathrm{0}})},\end{array}$$

where r₀ is the distance to the focus at x₀:

$$\begin{array}{}\text{(3)}& r_{\mathrm{0}}=\sqrt{(x_{\mathrm{k}}-x_{\mathrm{0}}{)}^{\mathrm{2}}+y_{\mathrm{k}}^{\mathrm{2}}+z_{\mathrm{k}}^{\mathrm{2}}}.\end{array}$$

The focus is along the x axis and is defined as

$$\begin{array}{}\text{(4)}& x_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{R}_{\mathrm{mp}}.\end{array}$$

x_k, y_k, and z_k are the Cartesian representation of the KF magnetosheath model (i.e., with the pre-fixed bow shock and magnetopause shapes) obtained by projecting the position vector onto the unit vectors e_x, e_y, and e_z:

$$\begin{array}{}\text{(5)}& {\displaystyle }{x}_{\mathrm{k}}& {\displaystyle }={\mathit{r}}^{\left(\mathrm{k}\right)}\cdot {\mathit{e}}_x,\text{(6)}& {\displaystyle }{y}_{\mathrm{k}}& {\displaystyle }={\mathit{r}}^{\left(\mathrm{k}\right)}\cdot {\mathit{e}}_y,\text{(7)}& {\displaystyle }{z}_{\mathrm{k}}& {\displaystyle }={\mathit{r}}^{\left(\mathrm{k}\right)}\cdot {\mathit{e}}_z.\end{array}$$

To complete the variable set for computing the potentials and the stream function, the azimuthal angle ϕ is introduced as

$$\begin{array}{}\text{(8)}& \mathit{\varphi }=\mathrm{atan}(z_{\mathrm{k}}/y_{\mathrm{k}}).\end{array}$$

2.2 Velocity potential

In the frame of potential field theory, the flow velocity U is obtained either from the velocity potential (a scalar potential) Φ^(vel) or from the stream function (also a scalar potential) Ψ as

$$\begin{array}{}\text{(9)}& \mathit{U}=-\mathrm{\nabla }{\mathrm{\Phi }}^{\left(\mathrm{vel}\right)}=-\mathrm{\nabla }\times \left(\mathrm{\Psi }{\mathit{e}}_{\mathit{\varphi }}\right).\end{array}$$

The symbol e_ϕ is the unit vector in the azimuthal direction around the symmetry axis (Sun-to-planet direction). Kobel and Flückiger (1994) and Guicking et al. (2012) obtained the analytic expression of the velocity potential Φ^(vel) using the shell variable v and the connector variable u.

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathrm{\Phi }}^{\left(\mathrm{vel}\right)}& =-U_x\left({\displaystyle \frac{v_{\mathrm{mp}}^{\mathrm{2}}{v}_{\mathrm{bs}}^{\mathrm{2}}}{v_{\mathrm{bs}}^{\mathrm{2}}-v_{\mathrm{mp}}^{\mathrm{2}}}}\right)\left({\displaystyle \frac{u^{\mathrm{2}}-v^{\mathrm{2}}}{\mathrm{2}{v}_{\mathrm{bs}}^{\mathrm{2}}}}+\ln v\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_x\left(u^{\mathrm{2}}-v^{\mathrm{2}}\right)+{\mathrm{\Phi }}_{\mathrm{0}}^{\left(\mathrm{vel}\right)},\end{array}\end{array}$$

where U_x is the upstream flow velocity, v_mp is the shell variable at the magnetopause, v_bs is the shell variable at the bow shock, v is the shell variable, u is the connector variable, and ${\mathrm{\Phi }}_{\mathrm{0}}^{\left(\mathrm{vel}\right)}$ is a free parameter (integration constant) which is set to zero without loss of generality. The boundary shell values v_mp and v_bs contain the information on the stand-off distances (R_mp and R_bs) in the subsolar region, and they are defined by Kobel and Flückiger (1994) as

$$\begin{array}{}\text{(11)}& {\displaystyle }{v}_{\mathrm{mp}}& {\displaystyle }=\sqrt{R_{\mathrm{mp}}},\text{(12)}& {\displaystyle }{v}_{\mathrm{bs}}& {\displaystyle }=\sqrt{\mathrm{2}{R}_{\mathrm{bs}}-R_{\mathrm{mp}}}.\end{array}$$

2.3 Stream function

Guicking et al. (2012) transformed the KF potential and obtained analytically the stream function Ψ as a function of the shell variable and the connector variable:

$$\begin{array}{}\text{(13)}& \mathrm{\Psi }=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_x\left({\displaystyle \frac{v_{\mathrm{mp}}^{\mathrm{2}}{v}_{\mathrm{bs}}^{\mathrm{2}}}{v_{\mathrm{bs}}^{\mathrm{2}}-v_{\mathrm{mp}}^{\mathrm{2}}}}\right){\displaystyle \frac{u}{v}}\left({\displaystyle \frac{v^{\mathrm{2}}}{v_{\mathrm{bs}}^{\mathrm{2}}}}-\mathrm{1}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{U}_{x}uv.\end{array}$$

Hereafter, one may set $U_x=-\mathrm{1}$ so that the velocity potential Φ^(vel) is normalized to the upstream velocity. Iso-contour lines of the stream function represent the streamline.

2.4 Magnetic scalar potential

The magnetic field in the magnetosheath is derived from the scalar potential in the same fashion as the flow velocity; that is,

$$\begin{array}{}\text{(14)}& \mathit{B}=-\mathrm{\nabla }{\mathrm{\Phi }}^{\left(\mathrm{mag}\right)}.\end{array}$$

The magnetic potential is a function of the shell variable v and the connector u (Kobel and Flückiger, 1994):

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathrm{\Phi }}^{\left(\mathrm{mag}\right)}& =-{\displaystyle \frac{v_{\mathrm{mp}}^{\mathrm{2}}{v}_{\mathrm{bs}}^{\mathrm{2}}}{v_{\mathrm{bs}}^{\mathrm{2}}-v_{\mathrm{mp}}^{\mathrm{2}}}}\times \\ & [\left(B_y^{\left(\mathrm{up}\right)}\cos \mathit{\varphi }+B_z^{\left(\mathrm{up}\right)}\sin \mathit{\varphi }\right)u\left({\displaystyle \frac{\mathrm{1}}{v}}+{\displaystyle \frac{v}{v_{\mathrm{bs}}^{\mathrm{2}}}}\right)\\ & +B_x^{\left(\mathrm{up}\right)}\left({\displaystyle \frac{u^{\mathrm{2}}-v^{\mathrm{2}}}{\mathrm{2}{v}_{\mathrm{bs}}^{\mathrm{2}}}}+\ln v\right)]\\ & -B_x^{\left(\mathrm{up}\right)}(-x)-B_y^{\left(\mathrm{up}\right)}y-B_z^{\left(\mathrm{up}\right)}z+{\mathrm{\Phi }}_{\mathrm{0}}^{\left(\mathrm{mag}\right)},\end{array}\end{array}$$

where $B_x^{\left(\mathrm{up}\right)}$ is the sunward component of the upstream magnetic field (corresponding to the GSE-X in near-Earth space), and $B_y^{\left(\mathrm{up}\right)}$ and $B_z^{\left(\mathrm{up}\right)}$ are two components of the upstream magnetic field perpendicular to the x direction. ϕ is the azimuthal angle of the position around the symmetry axis (the y direction is given by the angle ϕ=0). The integration constant is chosen as ${\mathrm{\Phi }}_{\mathrm{0}}^{\left(\mathrm{mag}\right)}=\mathrm{0}$. The magnetic potential cannot be further transformed into the form of a stream function since the magnetic field distribution is essentially three-dimensional in the magnetosheath.

3 Mapping method comparison

3.1 Mapping problem

Our task is to find the shell variable v and the connector u in the empirical magnetosheath by finding a suitable mapping of the position vector from the empirical magnetosheath (denoted by r) onto the KF parabolic system (denoted by r^(k)). The evaluated v and u are then readily used to obtain the scalar potentials and the stream function. The flow velocity and the magnetic field in the empirical magnetosheath are obtained by computing the gradient of the respective potential.

A practically useful mapping procedure to utilize the KF potential is proposed by Soucek and Escoubet (2012) by using the radial direction as a reference. While the radial mapping can reasonably (i.e., with a relatively high accuracy) transform the KF potential into the empirical magnetosheath domain on the dayside, the mapping quality becomes degraded in the flank region due to the conversion effect associated with the non-orthogonal grid construction. Our approach differs from the radial mapping by using the magnetopause-normal direction as a reference. The azimuthal coordinate ϕ is still orthogonal to the u and v coordinates. We briefly compare between the two mapping methods here.

3.2 Radial direction as reference

Soucek and Escoubet (2012) presented in their pioneering work an algorithm of radial mapping by transforming the KF magnetosheath model into a general, empirical magnetosheath shape by referring to the radial direction from the planet and scaling the radial position in the magnetosheath to the KF model. While the radial mapping can reasonably (i.e., with a relatively high accuracy) transform the KF potential into the empirical magnetosheath domain on the dayside, the mapping quality becomes degraded in the flank region due to the strongly non-orthogonal grids. Figure 1 displays a comparison of the radial grids between the KF magnetosheath model and the empirical magnetosheath model. The grids span the radial direction to the planet and transfinite interpolation between the bow shock and the magnetopause. The radial mapping has a drawback in a stronger grid non-orthogonality effect, which causes an artificial converging flow pattern in the flank region (velocity potential shown in Fig. 2) when the scalar potential is directly transformed. In the Soucek–Escoubet method, the problem of the flow conversion effect was avoided by solving the magnetohydrodynamics (MHD) Rankine–Hugoniot relation and tracking the streamline iteratively between the KF parabolic magnetosheath model and the empirical magnetosheath model.

Although the method introduced by Génot et al. (2011) and later adapted by Soucek and Escoubet (2012) is computationally less expensive than global magnetosheath simulations, the density and velocity fields from the bow shock to a given point in the magnetosheath still need to be computed along the streamline in an incremental way. Moreover, the Rankine–Hugoniot relations need to be solved before calculating iteratively the streamline, the flow velocity vector, to track the plasma density flow velocity along the streamline. Naturally, the uncertainty in this calculation depends on the step size (larger uncertainties for larger step sizes), and the errors accumulate along the streamline. The method introduced in this work overcomes this issue by constructing a magnetopause-normal grid system such that computational efforts are improved (no need to solve the Rankine–Hugoniot relations, and the error does not accumulate in the flank region).

Figure 1 Grid pattern generated by the radial mapping for the Kobel–Flückiger parabolic magnetosheath (a) and the non-parabolic, empirical magnetosheath (b).

Figure 2 Velocity potential in the Kobel–Flückiger model (a) and its radial mapping onto the non-parabolic empirical boundaries (b). Note that the same function is plotted here for different mapping methods. Contours represent the velocity potential normalized to the solar wind, ${\mathrm{\Phi }}^{\left(\mathrm{vel}\right)}/U_x$, which is in units of the planetary radii. The gray range is chosen for the visual demonstration purpose from 0.2R_E to 90R_E (a) and from 2R_E to 200R_E (b).

3.3 Magnetopause-normal direction as reference

Our mapping method differs from the radial mapping method in that the magnetopause-normal direction is used as a reference to the magnetopause. Our method guarantees the grid orthogonality around the magnetopause both on the dayside and in the flank region. The magnetopause-normal grids are shown in Fig. 3 for the KF magnetosheath model (with parabolic boundaries) and the empirical magnetosheath model (with non-parabolic boundaries). Even though the exact conformal mapping is not available, the magnetopause-normal mapping method retains the grid orthogonality around the magnetopause. This feature (orthogonality around the magnetopause) plays a crucial role in mapping the scalar potentials. An example of the scalar-potential mapping by referring to the magnetopause-normal direction (our final results) is shown in Sect. 4.8.

Figure 3 Mesh pattern used in the magnetopause-normal mapping in this work for the parabolic boundaries (a) and the non-parabolic, empirical boundaries (b).

4 Magnetopause-normal mapping

4.1 Overview of the procedure

The magnetopause-normal mapping is performed with two transformations. In the first transformation, the position vector is mapped from the empirical magnetosheath r onto the KF magnetosheath model r_k. This is based on the assumption that the distance to the magnetopause along the magnetopause-normal direction when scaled to the magnetosheath thickness (defined as the distance from the magnetopause to the bow shock along the magnetopause-normal direction) remains constant. The azimuthal angle ϕ is the same between the empirical magnetosheath and the KF model. The first transformation is divided into computing the distance to the magnetopause (step 1), the thickness of the empirical magnetosheath (step 2), the thickness of the KF magnetosheath (step 3), and the mapping of the position vector onto the KF model (step 4).

In the second transformation, the mapped position vector is used to compute the shell variable v and the connector variable u (step 5) and to obtain the potentials and the stream function in the empirical magnetosheath using Eqs. (10), (13), and (15) (step 6). Here again, the azimuthal angle ϕ is treated as the same.

Figure 4 illustrates the mapping procedure and graphically explains the variables that need to be determined to perform the mapping such as the zenith angle θ_mp associated with the minimum distance to the magnetopause the distance from the planet to the bow shock r_bs, the distance from the planet to the magnetopause r_mp, the relative distance to the magnetopause α_emp, and the magnetosheath thickness ${\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}$. The position vector r, the bow shock stand-off distance R_bs, the bow shock shape, the magnetopause stand-off distance R_mp, and the magnetopause shape are assumed to be known in our mapping.

Figure 4 Variables used in the magnetopause-normal mapping with the zenith angle θ_mp along the direction nearest to the magnetopause, the radial distance to the bow shock and magnetopause along the magnetosheath-normal direction (r_bs and r_mp, respectively), the distance from the magnetosheath to the magnetopause α_emp, and the magnetopause thickness ${\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}$. The position vector is denoted by r. The bow shock and magnetopause stand-off distances are denoted by R_bs and R_mp, respectively.

4.2 Setup

We begin with a position vector in the empirical magnetosheath domain, and we express the position vector as $\mathit{r}=x{\mathit{e}}_x+y{\mathit{e}}_y+z{\mathit{e}}_z$. Hereafter, we present the mapping procedure in the two-dimensional plane spanning the x and y directions for simplicity, but the computation in three dimensions is straightforward by representing the y component of position vector in the cylindrical fashion as ρcos ϕ and the z component as ρsin ϕ using the distance ρ to the x axis. The boundaries (bow shock and magnetopause) are specified by the users and do not need to be parabolic. In this paper, we use the following bow shock and magnetopause models.

• The empirical bow shock position expressed in GSE (Geocentric Solar Ecliptic) coordinates proposed and discussed by Farris et al. (1991) and Cairns et al. (1995):

  $$\begin{array}{}\text{(16)}& x=R_{\mathrm{bs}}-b_{\mathrm{emp}}{y}^{\mathrm{2}},\end{array}$$

  where R_bs is the bow shock stand-off distance, and b_emp is the empirical flaring parameter. We note here that the original Farris empirical bow shock model is not a paraboloid model; it is an ellipsoid model (with an eccentricity of 0.81), describing the bow shock on the dayside. It is not a proper representation of the far-flank bow shock. Also, the Cairns paraboloid bow shock model does not properly represent the far-flank bow shock. The distant bow shock shape approaches that of a hyperboloid.

• The empirical magnetopause position by Shue et al. (1997):

  $$\begin{array}{}\text{(17)}& x^{\mathrm{2}}+y^{\mathrm{2}}-{\displaystyle \frac{\mathrm{4}{R}_{\mathrm{mp}}^{\mathrm{4}}}{\mathrm{4}{R}_{\mathrm{mp}}^{\mathrm{2}}-y^{\mathrm{2}}}}=\mathrm{0},\end{array}$$

  in the Cartesian representation and

  $$\begin{array}{}\text{(18)}& r_{\mathrm{mp}}=R_{\mathrm{mp}}\sqrt{{\displaystyle \frac{\mathrm{2}}{\mathrm{1}+\cos \mathit{\theta }}}}\end{array}$$

  in the polar representation.

  We use a specific exponent for the Shue model (with an alpha exponent of 0.5) in an effort to show that the analytic model is “simple”. The solar wind conditions for which this exponent is applicable is not often encountered (e.g., interplanetary magnetic field has the B_z component larger than +8 nT, with specific values of solar wind dynamic pressure).

In our setup, the radial distance from the planet to the bow shock is expressed as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{r}_{\mathrm{bs}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{b}_{\mathrm{emp}}{\sin }^{\mathrm{2}}\mathit{\theta }}}(-\cos \mathit{\theta }+\\ & \sqrt{\mathrm{1}-(\mathrm{1}-\mathrm{4}{b}_{\mathrm{emp}}{R}_{\mathrm{bs}}){\sin }^{\mathrm{2}}\mathit{\theta }}).\end{array}\end{array}$$

By introducing the zenith angle θ and inserting x=r_bscos θ and y=r_bssin θ in Eq. (16), we obtain the equation for the radial distance to the empirical bow shock:

$$\begin{array}{}\text{(20)}& b_{\mathrm{emp}}{r}_{\mathrm{bs}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\theta }+r_{\mathrm{bs}}\cos \mathit{\theta }-R_{\mathrm{bs}}=\mathrm{0}.\end{array}$$

Equation (20) can be algebraically solved, and we take the positive value of the solution as presented in Eq. (19).

The radial distance to the magnetopause is given conveniently by Eq. (18). The Shue model reproduces the magnetopause stand-off distance R_mp in the subsolar direction (θ=0), and the cylindrical distance asymptotes to 2R_mp in the tail. It is worth noting here that one needs to compute the radial distance from the planet to the bow shock or magnetopause as a function of the zenith angle when using different shapes.

4.3 Step 1: measuring the distance to magnetopause

In the first step, the distance from the given position in the magnetosheath to the nearest magnetopause is computed (see Fig. 5). We express the position vector along the magnetopause-normal direction as

$$\begin{array}{}\text{(21)}& \mathit{r}={\mathit{r}}_{\mathrm{mp}}+{\mathit{\alpha }}_{\mathrm{emp}}{\mathit{e}}_{\mathrm{mp}},\end{array}$$

where r_mp is the magnetopause position nearest to the position vector, and e_mp is the unit vector in the magnetopause-normal direction. The unit vector points away from the planet and satisfies the condition

$$\begin{array}{}\text{(22)}& {\mathit{r}}_{\mathrm{mp}}\cdot {\mathit{e}}_{\mathrm{mp}}>\mathrm{0}.\end{array}$$

The symbol α_emp is the distance to the magnetopause along the magnetopause-normal direction e_mp in the empirical magnetosheath.

The nearest magnetopause position is obtained by searching for the zenith angle θ_mp for the minimum distance from the sample position to the magnetopause. The distance D is defined as

$$\begin{array}{}\text{(23)}& D=\sqrt{(r_x-r_{\mathrm{mp}}\cos {\mathit{\theta }}_{\mathrm{mp}}{)}^{\mathrm{2}}+(r_y-r_{\mathrm{mp}}\sin {\mathit{\theta }}_{\mathrm{mp}}{)}^{\mathrm{2}}}.\end{array}$$

The search for the minimum distance is implemented in a brute-force fashion as a function of μ_mp=cos θ_mp in our study.

Figure 5 Measuring the distance to the empirical magnetopause (step 1).

Using the minimum distance to the magnetopause r_mp and the zenith angle θ_mp, we are ready to compute the magnetopause-normal direction and the distance α_emp. To obtain the magnetopause-normal direction, we define the magnetopause shape function f_mp as

$$\begin{array}{}\text{(24)}& f_{\mathrm{mp}}=x^{\mathrm{2}}+y^{\mathrm{2}}-{\displaystyle \frac{\mathrm{4}{R}_{\mathrm{mp}}^{\mathrm{4}}}{\mathrm{4}{R}_{\mathrm{mp}}^{\mathrm{2}}-y^{\mathrm{2}}}},\end{array}$$

and we compute the normal direction by the gradient of f_mp as

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\partial f_{\mathrm{mp}}}{\partial x}}& {\displaystyle }=\mathrm{2}x,\text{(26)}& {\displaystyle \frac{\partial f_{\mathrm{mp}}}{\partial y}}& {\displaystyle }=\mathrm{2}y\left[\mathrm{1}-{\displaystyle \frac{\mathrm{4}{R}_{\mathrm{mp}}^{\mathrm{4}}}{(\mathrm{4}{R}_{\mathrm{mp}}^{\mathrm{2}}-y^{\mathrm{2}}{)}^{\mathrm{2}}}}\right].\end{array}$$

The magnetopause-normal direction is obtained by normalizing the gradient vector $({\partial }_x{f}_{\mathrm{mp}},{\partial }_y{f}_{\mathrm{mp}})$ and representing it with the basis vectors (e_x and e_y) as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\mathit{e}}_{\mathrm{mp}}& ={\displaystyle \frac{\mathrm{sgn}}{\sqrt{({\partial }_x{f}_{\mathrm{mp}}{)}^{\mathrm{2}}+({\partial }_y{f}_{\mathrm{mp}}{)}^{\mathrm{2}}}}}\times \\ & \left({\partial }_x{f}_{\mathrm{mp}}{\mathit{e}}_{\mathrm{x}}+{\partial }_y{f}_{\mathrm{mp}}{\mathit{e}}_{\mathrm{y}}\right),\end{array}\end{array}$$

evaluated at the magnetopause (x=r_mpcos θ_mp and y=r_mpsin θ_mp). The magnetopause-normal vector e_mp has a unit length, and the sign ($\mathrm{sgn}=\pm \mathrm{1}$) is chosen such that the normal vector is pointing outward (Eq. 22). The distance α_emp to the magnetopause along the normal direction is obtained from Eq. (21) as

$$\begin{array}{}\text{(28)}& {\mathit{\alpha }}_{\mathrm{emp}}={\displaystyle \frac{(x-r_{\mathrm{mp}}\cos {\mathit{\theta }}_{\mathrm{mp}})+(y-r_{\mathrm{mp}}\sin {\mathit{\theta }}_{\mathrm{mp}})}{{\mathit{e}}_{\mathrm{mp}}\cdot {\mathit{e}}_x+{\mathit{e}}_{\mathrm{mp}}\cdot {\mathit{e}}_y}}.\end{array}$$

Equation (28) is constructed to be robust against the singular behavior on the dayside (${\mathit{e}}_{\mathrm{mp}}\cdot {\mathit{e}}_y=\mathrm{0}$) and in distant tail (${\mathit{e}}_{\mathrm{mp}}\cdot {\mathit{e}}_x=\mathrm{0}$).

4.4 Step 2: computing the thickness of empirical magnetosheath

In the second step, the magnetosheath thickness is computed using the position vector and the magnetopause normal direction (Fig. 6). For our mapping purpose, the distance α_emp is normalized to the magnetosheath thickness ${\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}$ such that the relative distance ${\mathit{\alpha }}_{\mathrm{emp}}/{\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}$ serves as an invariant of the mapping from the empirical magnetosheath onto the KF magnetosheath. To achieve this, we combine Eq. (16) with Eq. (21) and analytically determine the thickness from the bow shock to the magnetopause in the empirical magnetosheath. That is, the thickness ${\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}$ is obtained by rewriting the bow shock quadratic equation (Eq. 16) for the position vector using the variable ${\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}$ (Eq. 21) extended to the bow shock location. The equation is again quadratic, and the solution is algebraically obtained as:

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{b}_{\mathrm{emp}}{e}_{\mathrm{mp},y}^{\mathrm{2}}}}\times \\ & \left[-(e_{\mathrm{mp},x}+\mathrm{2}{b}_{\mathrm{emp}}{y}_{\mathrm{mp}}{e}_{\mathrm{mp},y}{)}^{\mathrm{2}}+d_{\mathit{\alpha }}\right],\end{array}\end{array}$$

where d_α is an auxiliary variable defined as

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{d}_{\mathit{\alpha }}& =\left[\right(e_{\mathrm{mp},x}+\mathrm{2}{b}_{\mathrm{emp}}{y}_{\mathrm{mp}}{e}_{\mathrm{mp},y}{)}^{\mathrm{2}}-\\ & \mathrm{4}{b}_{\mathrm{emp}}{e}_{\mathrm{mp},y}^{\mathrm{2}}\times \\ & (x_{\mathrm{mp}}+b_{\mathrm{emp}}{y}_{\mathrm{mp}}^{\mathrm{2}}-R_{\mathrm{bs}}){]}^{\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

In the subsolar direction (y_mp=0), the thickness is simply given as

$$\begin{array}{}\text{(31)}& {\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}=R_{\mathrm{bs}}-R_{\mathrm{mp}}.\end{array}$$

Equation (29) becomes singular in the subsolar direction, and Eq. (31) needs to be set separately.

Figure 6 Computing the magnetosheath thickness in the empirical model (step 2).

4.5 Step 3: computing the magnetosheath thickness in the KF system

In the third step, the magnetosheath thickness is computed in the KF model (Fig. 7). We repeat the procedures of steps 1 and 2 for the KF system and determine the KF magnetosheath thickness as reference. We treat the zenith angle θ_mp and the relative distance ${\mathit{\alpha }}_{\mathrm{emp}}/{\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}$ as invariants of the mapping between the empirical magnetosheath and the KF system. The KF bow shock location is given as

$$\begin{array}{}\text{(32)}& x=R_{\mathrm{bs}}-b_{\mathrm{k}}{y}^{\mathrm{2}},\end{array}$$

where the KF bow-shock flaring parameter b_k is pre-fixed as (Kobel and Flückiger, 1994)

$$\begin{array}{}\text{(33)}& b_{\mathrm{k}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{R}_{\mathrm{bs}}-\mathrm{2}{R}_{\mathrm{mp}}}}.\end{array}$$

The radial distance from the planet to the KF bow shock is

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{r}_{\mathrm{bs}}^{\left(\mathrm{k}\right)}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{b}_{\mathrm{k}}{\sin }^{\mathrm{2}}\mathit{\theta }}}\times \\ & \left(-\cos \mathit{\theta }+\sqrt{\mathrm{1}+(\mathrm{4}{b}_{\mathrm{k}}{R}_{\mathrm{bs}}-\mathrm{1}){\sin }^{\mathrm{2}}\mathit{\theta }}\right).\end{array}\end{array}$$

The KF magnetopause is defined in Kobel and Flückiger (1994) as

$$\begin{array}{}\text{(35)}& x=R_{\mathrm{mp}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{R}_{\mathrm{mp}}}}{y}^{\mathrm{2}}.\end{array}$$

From Eq. (35), the radial distance from the planet to the KF magnetopause is computed as

$$\begin{array}{}\text{(36)}& r_{\mathrm{mp}}^{\left(\mathrm{k}\right)}={\displaystyle \frac{R_{\mathrm{mp}}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}\left(-\cos \mathit{\theta }+\sqrt{\mathrm{1}+{\sin }^{\mathrm{2}}\mathit{\theta }}\right).\end{array}$$

To obtain the magnetopause-normal direction in the KF system, we compute the gradient of the magnetopause shape function:

$$\begin{array}{}\text{(37)}& f_{\mathrm{mp}}^{\left(\mathrm{k}\right)}=x-R_{\mathrm{mp}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{R}_{\mathrm{mp}}}}{y}^{\mathrm{2}}.\end{array}$$

The gradient is analytically given as

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{\partial f_{\mathrm{mp}}^{\left(\mathrm{k}\right)}}{\partial x}}& {\displaystyle }=\mathrm{1}.\text{(39)}& {\displaystyle \frac{\partial f_{\mathrm{mp}}^{\left(\mathrm{k}\right)}}{\partial y}}& {\displaystyle }={\displaystyle \frac{y}{R_{\mathrm{mp}}}}.\end{array}$$

The magnetopause-normal direction ${\mathit{e}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}$ is then obtained by applying Eqs. (38) and (39) to Eq. (27), which reads as

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{\mathit{e}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}& ={\displaystyle \frac{\mathrm{sgn}}{\sqrt{({\partial }_x{f}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}{)}^{\mathrm{2}}+({\partial }_y{f}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}{)}^{\mathrm{2}}}}}\times \\ & \left({\partial }_x{f}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}{\mathit{e}}_{\mathrm{x}}+{\partial }_y{f}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}{\mathit{e}}_{\mathrm{y}}\right).\end{array}\end{array}$$

Figure 7 Computing the magnetosheath thickness in the KF model (step 3). The same zenith angle as that in step 2 is used.

The thickness in the KF system ${\mathit{\alpha }}_{\mathrm{k}}^{\left(\mathrm{bs}\right)}$ is determined by combining the bow shock shape (Eq. 32) with the position vector at the bow shock:

$$\begin{array}{}\text{(41)}& {\mathit{r}}_{\mathrm{bs}}^{\left(\mathrm{k}\right)}={\mathit{r}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}+{\mathit{\alpha }}_{\mathrm{k}}^{\left(\mathrm{bs}\right)}{\mathit{e}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}.\end{array}$$

Equation (32) becomes again a quadratic equation with respect to the thickness ${\mathit{\alpha }}_{\mathrm{k}}^{\left(\mathrm{bs}\right)}$, and the solution reads

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{k}}^{\left(\mathrm{bs}\right)}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{b}_{\mathrm{k}}(e_{\mathrm{mp},y}^{\left(\mathrm{k}\right)}{)}^{\mathrm{2}}}}\times \\ & [-(e_{\mathrm{mp},x}^{\left(\mathrm{k}\right)}+\mathrm{2}{b}_{\mathrm{k}}{y}_{\mathrm{mp}}{e}_{\mathrm{mp},y}^{\left(\mathrm{k}\right)})+d_{\mathit{\alpha }}^{\left(\mathrm{k}\right)}],\end{array}\end{array}$$

where the auxiliary variable $d_{\mathit{\alpha }}^{\left(\mathrm{k}\right)}$ is defined as

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{d}_{\mathit{\alpha }}^{\left(\mathrm{k}\right)}& =\left[\right(e_{\mathrm{mp},x}^{\left(\mathrm{k}\right)}+\mathrm{2}{b}_{\mathrm{k}}{y}_{\mathrm{mp}}{e}_{\mathrm{mp},y}^{\left(\mathrm{k}\right)}{)}^{\mathrm{2}}-\\ & -\mathrm{4}{b}_{\mathrm{k}}\left(e_{\mathrm{mp},y}^{\left(\mathrm{k}\right)}{)}^{\mathrm{2}}\right(x_{\mathrm{mp}}+b_{\mathrm{k}}{y}_{\mathrm{mp}}^{\mathrm{2}}-R_{\mathrm{bs}}){]}^{\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

4.6 Step 4: mapping the position vector onto the KF system

In the fourth step, the mapping of the position vector is performed from the empirical magnetosheath onto the KF system (Fig. 8). An assumption is made such that the relative distance to the magnetopause along the magnetopause-normal direction is the same between the two systems. The distance from the magnetosheath position vector to the magnetopause along the magnetopause-normal direction in the KF system α_k is then determined by the relative distance in the empirical magnetosheath α_emp, the thickness of the empirical magnetosheath ${\mathit{\alpha }}_{\mathrm{emp}}^{\mathrm{bs}}$, and magnetosheath thickness in the KF system ${\mathit{\alpha }}_{\mathrm{k}}^{\left(\mathrm{bs}\right)}$ as

$$\begin{array}{}\text{(44)}& {\mathit{\alpha }}_{\mathrm{k}}={\mathit{\alpha }}_{\mathrm{emp}}{\mathit{\alpha }}_{\mathrm{k}}^{\left(\mathrm{bs}\right)}/{\mathit{\alpha }}_{\mathrm{emp}}^{\left(\mathrm{bs}\right)}.\end{array}$$

The mapped position vector is then computed as

$$\begin{array}{}\text{(45)}& {\mathit{r}}^{\left(\mathrm{k}\right)}={\mathit{r}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}+{\mathit{\alpha }}_{\mathrm{k}}{\mathit{e}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)},\end{array}$$

using the nearest magnetopause position ${\mathit{r}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}$ (Eq. 36), the magnetosheath-to-magnetopause distance α_k (Eq. 44), and the magnetopause-normal direction ${\mathit{e}}_{\mathrm{mp}}^{\left(\mathrm{k}\right)}$ (Eq. 40).

Figure 8 Mapping the position vector onto the KF magnetosheath model (step 4).

4.7 Step 5: evaluating the shell and connector variables

In the fifth step, the shell variable v and the connector variable u are computed from the mapped position vector r^(k) using Eqs. (1) and (2), respectively. The variables v and u are the same as the parabolic coordinates used in the KF potential with a focus at $x_{\mathrm{0}}=R_{\mathrm{mp}}/\mathrm{2}$. In our algorithm, the focus is explicitly given in Eqs. (1), (2), and (3). The azimuthal angle around the symmetry axis ϕ is treated in the same way as in the KF model.

Figure 9 Evaluating the shell variable v and the connector variable u in the KF magnetosheath model (step 5).

Figure 10 compares the iso-contours of the shell v and the connector u represented in the KF system (left panel) and the empirical magnetosheath (right panel) for a bow shock stand-off distance of 12.8 R_E (Génot et al., 2011), a bow shock flaring of $\mathrm{0.0223}{R}_{\mathrm{E}}^{-\mathrm{1}}$ (Farris et al., 1991; Cairns et al., 1995), and a magnetopause stand-off distance 9.8 R_E (Génot et al., 2011). The shell variable v is characterized by the lines with the curvature center on the right side in the panel and contains the parabolic bow shock (at v=v_bs) and magnetopause (at v=v_mp) marked by thick lines. The connector variable u has the curvature center on the left side in the panel, and the iso-contour lines are orthogonal to the bow shock and magnetopause. The computation of the u and v variables and their gradient and curl is performed in Cartesian coordinates so that the connection represented by the Christoffel symbol vanishes in the computation. Computation in the Cartesian domain is also beneficial to the practical application because spacecraft trajectories are often represented in Cartesian coordinates.

Figure 10 Iso-contour lines with $u=\mathrm{const}.$ (center of curvature on the left side) and that with $v=\mathrm{const}.$ (center of curvature on the right side) in the KF magnetosheath model (a) and the empirical magnetosheath model (b). The bow shock stand-off distance is 12.8 Earth radii and the magnetopause stand-off distance is 9.8 Earth radii.

4.8 Step 6: computing the potentials and stream function

The scalar potentials (velocity potential and magnetic potential) and the stream function are obtained from the shell v and the connector u using Eqs. (10), (13), and (15). The velocity potential (normalized to the upstream flow) is displayed in Fig. 11 left panel, and the stream function is displayed in the right panel. The iso-contours of the velocity potential represent the lines of the same flow velocity. The iso-contours of the stream function represent the streamlines in the magnetosheath. The flow is deflected around the nose of the magnetopause (the subsolar point), and the streamlines are tangential to the magnetopause. Due to the grid orthogonality around the magnetopause, the streamlines are constructed as tangential to the magnetopause shape, which qualifies the magnetopause-normal mapping method as a useful tool for the magnetosheath model.

Figure 11 Velocity potential (a) and stream function (b) in the empirical magnetosheath domain obtained by mapping onto the shell variable v and the connector variable u. Note that two different functions are plotted here for the same grid and mapping method. The gray range is chosen for the visual demonstration from 5.5R_E to 314R_E (a) and from 0R_E to 15R_E (b).

The magnetic potential and the derived magnetic field are displayed in Fig. 12. The magnetic potential and the magnetic field (the gradient of the potential multiplied by the minus sign) depend on the upstream field. Figure 12 shows an example with an upstream field angle of 135° to the x axis (i.e., 45° to the upstream flow direction). The magnetic field is computed using the central difference scheme. Near the boundaries (bow shock or magnetopause), the mesh resolution is enhanced so that the mesh points do not cross the boundary when computing with the central difference scheme. The upstream field is deflected on the positive y side (right panel, lower half plane) and is draping the magnetopause on the negative y side (right panel, upper half plane).

Figure 12 Magnetic potential for the upstream magnetic field with an angle of 135° to the x axis (45° to the upstream flow direction (a) and sampled magnetic field vectors obtained by the negative gradient of the magnetic potential (b). The gray range of the magnetic potential is from −35R_E to 348R_E.

5 Discussion

5.1 Extension of the mapping approach

It is straightforward to extend our method to different shapes of the Shue magnetopause model. The following form of gradient can be used for a general value of the Shue exponent α_mp:

$$\begin{array}{}\text{(46)}& {\displaystyle \frac{\partial f_{\mathrm{mp}}}{\partial x}}=\cos \mathit{\theta }+{\displaystyle \frac{{\mathrm{2}}_{\mathrm{mp}}^{\mathit{\alpha }}{R}_{\mathrm{mp}}\mathit{\alpha }}{r}}{\left(\mathrm{1}+\cos \mathit{\theta }\right)}^{-({\mathit{\alpha }}_{\mathrm{mp}}+\mathrm{1})}{\sin }^{\mathrm{2}}\mathit{\theta },\text{(47)}& {\displaystyle \frac{\partial f_{\mathrm{mp}}}{\partial y}}=\sin \mathit{\theta }-{\displaystyle \frac{{\mathrm{2}}_{\mathrm{mp}}^{\mathit{\alpha }}{R}_{\mathrm{mp}}\mathit{\alpha }}{r}}{\left(\mathrm{1}+\cos \mathit{\theta }\right)}^{-({\mathit{\alpha }}_{\mathrm{mp}}+\mathrm{1})}\cos \mathit{\theta }\sin \mathit{\theta },\end{array}$$

where f_mp in Eqs. (46)–(47) is defined as

$$\begin{array}{}\text{(48)}& f_{\mathrm{mp}}=r-R_{\mathrm{mp}}{\left({\displaystyle \frac{\mathrm{2}}{\mathrm{1}+\cos \mathit{\theta }}}\right)}^{{\mathit{\alpha }}_{\mathrm{mp}}}.\end{array}$$

Mesh pattern applied to different values of the Shue exponent for an open-type magnetopause α_mp=0.6 (Fig. 13a) and a closed-type magnetopause α_mp=0.4 (Fig. 13b).

Figure 13 Mesh pattern applied to different values of the Shue exponent for an open-type magnetopause α_mp=0.6 (a) and a closed-type magnetopause α_mp=0.4 (b).

Also, our method can be extended to a three-dimensional, non-axisymmetric geometry of magnetosheath (e.g., Dimmock and Nykyri, 2013). To achieve this goal, a suitable set of the variables needs to be found for the mapping: the shell variable, the connector variable, and the azimuthal angle around the subsolar axis (solar wind direction intersecting the planetary magnetic dipole). Namely, one needs to give a non-axisymmetric bow shock shape and a non-axisymmetric magnetopause shape; compute the magnetopause-normal direction and construct grids in the magnetosheath; measure the distance to the magnetopause and the bow shock; scale and relate the distance to the KF model; and evaluate the scalar potentials through the u, v, and ϕ variables.

5.2 Other approaches

It is possible to obtain the steady-state magnetosheath potential in different ways.

• First, one may numerically solve the Laplace equation for a given set of boundary shapes (bow shock and magnetopause). Various numerical solvers are known for solving the Laplace equation such as the Jacobi method, the Gauss–Seidel method, and the successive over-relaxation (SOR) method. These Laplace solvers are numerically more expensive than the mapping method, but the computation in 3-D is feasible with the contemporary computational resources. On the other hand, the magnetosheath is not bounded but extends in the tail direction. The challenge here is thus to construct a properly bounded area for the Laplace equation.

• Second, one may expand the magnetosheath magnetic field in different orthogonal functions. The KF solution makes use of the Bessel functions (Kobel and Flückiger, 1994). For flexible magnetopause and bow shock boundary models, a magnetosheath magnetic field model is constructed by making use of Legendre polynomials (Romashets and Vandas, 2019).

• Third, one may introduce a suitable conformal mapping by limiting the magnetosheath modeling to a complex plane (two-dimensional domain). The harmonic functions such as the KF solution are transformed from the parabolic magnetosheath shape into the non-axisymmetric magnetosheath shape. The problem here is that finding the conformal mapping is not an easy task, because the magnetosheath is not a spatially-closed domain and one has to set the boundary in the magnetosheath to complete the domain bounded by the bow shock, the subsolar axis, and the magnetopause.

• Fourth, one may solve the Rankine–Hugoniot relations and track the streamline stepwise by referring to the KF solution, as is done in Soucek and Escoubet (2012). This method computes the magnetosheath flow and magnetic field along the streamline, and the computing for the entire magnetosheath domain is numerically expensive.

6 Concluding remarks

Our potential mapping method may be regarded as an updated version of the radial mapping method (Soucek and Escoubet, 2012) by retaining the orthogonality near the magnetopause in the flank to tail region and also by computing the field through the potential mapping. Velocity potential, stream function, and magnetic potential are evaluated in the empirical magnetosheath.

The advantages of our methods are as follows:

1. The method makes extensive use of the exact solution of the Laplace equation (the Kobel–Flückiger potential and the Guicking stream function). The plasma flow and magnetic field can be determined semi-analytically in a wide range of zenith angles in the magnetosheath when the solar wind conditions and the boundary shapes are given.

2. The method is applicable to a non-parabolic shape of magnetosheath domain, opening the door to develop a tool to assist numerical simulations and spacecraft observations of not only the Earth but also the planetary magnetosheath domain.

3. The method is computationally inexpensive. In particular, if the shapes of bow shock and magnetopause are analytically given, most of the computational steps in the potential mapping method have an analytic expression.

As stated in Sect. 1, one ideally needs to find a conformal mapping from the KF magnetosheath model onto the empirical magnetosheath. While the conformal mapping is known both for the empirical bow shock and the empirical magnetopause, the conformal mapping of the entire magnetosheath domain still remains a challenge. There are two problems with this approach. First, the closing boundary (the u contours) connecting between the bow shock and the magnetopause is not known, and moreover, the uniqueness of finding such a boundary is not guaranteed. Second, the gradients along u are not the same between the empirical bow shock and the empirical magnetopause such that a naive transfinite interpolation ends up with highly non-orthogonal grids in the magnetosheath.

Our method of computing the plasma flow and magnetic field should be compared against the observations and simulations. For example, THEMIS and ARTEMIS spacecraft (Angelopoulos, 2008) and MMS spacecraft (Burch et al., 2016) are providing a huge amount of data on both sides of the bow shock in the equatorial plane; Cluster spacecraft (Escoubet et al., 2001) are collecting data in polar orbit; ACE spacecraft data (Stone et al., 1998) may be used as an upstream monitor; and Earth Flyby data of planetary missions (such as Cassini, BepiColombo) cover the far-distance tail region. In reality, non-axisymmetric structure arises in the magnetosheath. There is no restriction regarding the choice of the magnetopause model. The magnetopause-normal direction needs to be computed either analytically using the gradient of the magnetopause function as ∇f_mp or numerically for a user-defined magnetopause shape.