An Improved Analytical Model of the Local Interstellar Magnetic Field: The Extension to Compressibility

TL;DR

This paper extends an existing analytical model of the local interstellar magnetic field near the heliopause from incompressible to compressible flow, incorporating isothermal and adiabatic conditions, and validates it against MHD simulations.

Contribution

It introduces a compressible extension of a previous magnetic field model, providing exact expressions and approximations for density and flow velocity, and compares results with MHD simulations.

Findings

The extended model accurately describes magnetic draping in compressible flows.

Analytical solutions align well with self-consistent MHD simulation results.

The model enhances understanding of interstellar magnetic field behavior near the heliopause.

Abstract

A previously published analytical magnetohydrodynamic model for the local interstellar magnetic field in the vicinity of the heliopause (R\"oken et al. 2015) is extended from incompressible to compressible, yet predominantly subsonic flow, considering both isothermal and adiabatic equations of state. Exact expressions and suitable approximations for the density and the flow velocity are derived and discussed. In addition to the stationary induction equation, these expressions also satisfy the momentum balance equation along stream lines. The practical usefulness of the corresponding, still exact analytical magnetic field solution is assessed by comparing it quantitatively to results from a fully self-consistent magnetohydrodynamic simulation of the interstellar magnetic field draping around the heliopause.

Equations

$$\Phi (\mathbf{r}):=u_0\left(\frac{q}{r}+z\right)$$

$${\mathbf{u}}_{\mathrm{R}}(\mathbf{r}):=-\nabla \Phi (\mathbf{r})=u_0\left[\frac{q\rho }{r^3}{\mathbf{e}}_{\rho }+\left(\frac{qz}{r^3}-1\right){\mathbf{e}}_z\right],$$

$$H(\rho ,z):=2q-{\rho }^2-z\sqrt{4q-{\rho }^2}=0.$$

$$\nabla \times [{\mathbf{u}}_{\mathrm{R}}\times \mathbf{B}]=0$$

$$\nabla \cdot \mathbf{B}=0,$$

$$\mathbf{s}=-\nabla \Phi (\mathbf{r}).$$

$$m_{\mathrm{p}}n(\mathbf{u}\cdot \nabla )\mathbf{u}=-\nabla P$$

$$P(n)=C{n}^{\gamma },$$

$$-\mathbf{s}\cdot \nabla P$$

$$\nabla P=\left(\frac{\mathrm{d}P}{\mathrm{d}n}\right)\nabla n=C\gamma n^{\gamma -1}\nabla n$$

$$2\left[\frac{{\mathbf{s}}^2}{n}-\frac{C\gamma }{m_{\mathrm{p}}}{n}^{\gamma }\right]\mathbf{s}\cdot \nabla n=\mathbf{s}\cdot \nabla \left({\mathbf{s}}^2\right).$$

$$\mathbf{s}\cdot \nabla =s_{\rho }{\partial }_{\rho }+s_z{\partial }_z=p{\partial }_{\alpha }=p\left(\frac{\mathrm{d}\rho }{\mathrm{d}\alpha }{\partial }_{\rho }+\frac{\mathrm{d}z}{\mathrm{d}\alpha }{\partial }_z\right)$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}\alpha }=\frac{s_{\rho }}{p}\text{mbox}and\frac{\mathrm{d}z}{\mathrm{d}\alpha }=\frac{s_z}{p}.$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}\alpha }=1\iff \rho =\alpha +\mathcal{F}(\beta )$$

$$\frac{\mathrm{d}z}{\mathrm{d}\alpha }=\frac{s_z}{s_{\rho }},$$

$$\left[-\frac{m_{\mathrm{p}}}{C\gamma n_0^{\gamma +1}}{\mathbf{s}}^2+{\left(\frac{n}{n_0}\right)}^{\gamma +1}\right]\frac{2}{n}\frac{\partial n}{\partial \alpha }=-\frac{m_{\mathrm{p}}}{C\gamma n_0^{\gamma +1}}\frac{\partial {\mathbf{s}}^2}{\partial \alpha }.$$

$$g:=-\frac{m_{\mathrm{p}}}{C\gamma n_0^{\gamma +1}}{\mathbf{s}}^2,$$

$$\left[g+{\left(\frac{n}{n_0}\right)}^{\gamma +1}\right]\frac{2}{n}\frac{\partial n}{\partial \alpha }=\frac{\partial g}{\partial \alpha }.$$

$$\left[g+{\bar{n}}^{\gamma +1}\right]\frac{2}{\bar{n}}\frac{\mathrm{d}\bar{n}}{\mathrm{d}g}=1$$

$$c:={\left(\frac{\mathrm{d}P}{\mathrm{d}(m_{\mathrm{p}}n)}\right)}^{1/2}={\left(\frac{C\gamma }{m_{\mathrm{p}}}{n}^{\gamma -1}\right)}^{1/2},$$

$$C\gamma n_0^{\gamma +1}=m_{\mathrm{p}}{n}_0^2{c}_0^2.$$

$$g=-m^2\frac{{\mathbf{s}}^2}{{\mathbf{s}}_0^2}=-m^2{\bar{\mathbf{s}}}^2,$$

$$\begin{array}{rl}{\bar{\mathbf{s}}}^2=& \frac{1}{{\mathbf{s}}_0^2}{\left(-\nabla \Phi \right)}^2={\left[\left(\frac{q\rho }{r^3}\right){\mathbf{e}}_{\rho }+\left(\frac{qz}{r^3}-1\right){\mathbf{e}}_z\right]}^2\\ =& 1-\frac{2qz}{r^3}+\frac{q^2}{r^4}=:\mathcal{A}.\end{array}$$

$$\begin{array}{rl}g+{\bar{n}}^{\gamma +1}=& -m^2\frac{{\mathbf{s}}^2}{{\mathbf{s}}_0^2}+{\left(\frac{n}{n_0}\right)}^{\gamma +1}\\ =& -\frac{{\mathbf{u}}_0^2}{c_0^2}\frac{(n\mathbf{u}{)}^2}{(n_0{\mathbf{u}}_0{)}^2}+{\left(\frac{n}{n_0}\right)}^2{\left(\frac{c}{c_0}\right)}^2\\ =& {\left(\frac{n}{n_0}\right)}^2\frac{c^2-{\mathbf{u}}^2}{c_0^2},\end{array}$$

$$\frac{\mathrm{d}g}{\mathrm{d}\bar{n}}-\frac{2g}{\bar{n}}=2{\bar{n}}^{\gamma }.$$

$$\frac{\mathrm{d}{g}_{\mathrm{hom}}}{\mathrm{d}\bar{n}}-\frac{2g_{\mathrm{hom}}}{\bar{n}}=0,$$

$$g_{\mathrm{hom}}=c_1{\bar{n}}^2,$$

$$\frac{\mathrm{d}h}{\mathrm{d}\bar{n}}=2{\bar{n}}^{\gamma -2}$$

$$g_{\rm inhom}={\bar{n}}^{2}\times\left\{\begin{array}[]{ccc}\ln\left({\bar{n}}^{2}\right)+c_{2}&:&\gamma=1\\ &&\\ \displaystyle\frac{2\,{\bar{n}}^{\gamma-1}}{\gamma-1}+c_{2}&:&\gamma\neq 1\end{array}\right.$$

$${\bar{n}}(g)=\exp\left[\frac{F\big{(}g\,\exp(c_{2})\big{)}}{2}\right]$$

Full text

1. Introduction

Recently, the problem of an exact analytical magnetohydrodynamical (MHD) solution for an idealized structure of the local interstellar magnetic field draping around the heliopause was solved by Röken et al. (2015), hereafter referred to as Paper I. (See Isenberg et al. (2015) for an approximate, singularity-hampered solution.) It was obtained under the assumptions that (i) the field is frozen into an (ii) axisymmetric and (iii) incompressible interstellar plasma flow. The assumption of axisymmetry has been dropped in the second paper of this series (Kleimann et al., 2016) by employing so-called distortion flows that allow for more realistic cross sections of the heliospheric tail flattened by the interstellar magnetic field. In that paper, the assumption of incompressibility was addressed briefly by pointing out that the use of a solenoidal distortion flow does not imply any constraints on the compressibility of the interstellar plasma flows. In the present, third paper of the series, we generalize the solution presented in Paper I to the case of a compressible plasma flow.

The interest in such solutions lies in the fact that, on the one hand, the interstellar plasma is likely to be supersonic and super-fast-magnetosonic (Ben-Jaffel et al., 2013; Scherer & Fichtner, 2014) and, thus, is compressed at an interstellar bow shock. (For an analytical treatment of bow shocks in other astrophysical scenarios, see the recent paper by Keshet & Naor (2016).) The resulting subsonic low-Mach number flow in the outer heliosheath between the heliopause and the bow shock is limited to a relatively narrow region bounded by the so-called sonic lines (e.g., Scherer et al., 2016), at which the flow again becomes supersonic toward the flanks of the heliosphere. The ensuing flow compressibility in this region has been discussed in the context of the stability of the heliopause by, e.g., Caillol & Ruderman (2007) and Belov & Ruderman (2010). On the other hand, even for the case that such a bow shock would not exist (McComas et al., 2012), the plasma flow in the region of the so-called bow wave (Zank et al., 2013) would, at least partly, be characterized by Mach numbers below but close to unity (Gayley et al., 1997), and the usual incompressibility assumption for subsonic flows with low Mach numbers ($M\lesssim 0.3$) would not hold. So, in any case, there are reasons to drop the strict incompressibility assumption for the local interstellar medium.

As is demonstrated here, physically meaningful solutions for compressible, yet predominantly subsonic flows can be worked out analytically. These density structures and their approximations are derived and discussed in Section 2. In Section 3, we explicitly compute the resulting improved magnetic field frozen into such an approximative flow by exploiting the exact result of Paper I, taking into account modifications induced by the non-constant density. The usefulness of this improved analytical field solution is assessed in Section 4 by comparing it to a fully self-consistent numerical MHD model. A summary and conclusions are given in Section 5.

2. Density Structure

In this section, we determine the physical (number) density structure of the compressible model and discuss suitable approximations.

2.1. The Rankine-type Heliosphere Model

Since its introduction as a simple model for the interstellar flow in and around the heliopause by Parker (1961), the incompressible Rankine half-body flow continues to be popular in the heliophysics community (Yu, 1974; Nerney & Suess, 1995; Fahr et al., 2014, 2016; Isenberg et al., 2015; Sylla & Fichtner, 2015; Zirnstein & McComas, 2015). At its core lies a point-like mass source of strength $(4\pi u_{0})q$ located at the origin, superimposed on an undisturbed flow ${\bf u}_{0}=-u_{0}\,{\bf e}_{z}$ that is incident from the $+z$ direction. The Rankine velocity field thus derives from a flow potential

[TABLE]

as

[TABLE]

where $(\rho,\varphi,z)$ are cylindrical coordinates, and $r:=\|{\bf r}\|$. It can easily be shown that the heliopause, defined as the set of all flow lines emanating from the stagnation point $(\rho,z)_{\rm sp}:=(0,\sqrt{q})$, is given by the surface

[TABLE]

The problem of an interstellar magnetic field being passively advected in the Rankine flow (2) was recently considered in Paper I, which provided the derivation of the exact analytical magnetic field solution of the steady-state induction equation

[TABLE]

of ideal MHD and the magnetic divergence constraint

[TABLE]

subject to boundary conditions consisting of an arbitrarily inclined, homogeneous magnetic field at upstream infinity.

2.2. The Extension to Compressibility

While this idealized magnetic field was shown to yield a reasonable approximation to corresponding results from a self-consistent numerical model, one weakness that was identified in the comparison is that the pile-up of magnetic flux ahead of the heliopause is restricted to a relatively narrow layer, whereas this layer appears much broader in the numerical simulations, at least if parameters are chosen such that no bow shock forms. The cause of this shortcoming can be traced back to the fact that the underlying flow field (2) is incompressible and stems from an ad-hoc choice which, however reasonable, does not honor any conservation laws except that for mass.

In the improved model, which we present in this work, we relax the condition of incompressibility in favor of a physically more realistic description, while at the same time retaining as much as possible of the original flow structure. The key idea here is that while we continue to employ the flow potential (1), it is now re-interpreted as a potential for the particle flux density ${\bf s}:=n\,{\bf u}$, rather than for the velocity ${\bf u}$, where the number density $n$ is no longer constant, but may vary in space, such that Eq. (2) is replaced by

[TABLE]

Consequently, we observe that

the flow line structure, and in particular the shape of the heliopause as given by Eq. (3), remains unchanged, 2. 2.

the solenoidality of $\nabla\Phi$ now implies $\nabla\cdot(n\,{\bf u})=0$, i.e., mass continues to be conserved, and 3. 3.

the flow ${\bf u}$ itself is no longer incompressible (or irrotational), and the density may vary along stream lines.

Being derived from a global potential, the particle flux density is bound to be differentiable in the entire domain of interest. This implies that our model cannot accommodate oblique shocks, and thus, in particular, no bow shock ahead of which ${\bf s}$ could be genuinely undisturbed (i.e., constant); it is a model for predominantly subsonic flow. However, it is worth noting that the use of a flow potential as such does not preclude supersonic or even transsonic flow as long as the physical quantity that is computed from the potential’s gradient can rightfully be taken as irrotational and void of boundary layers (see, e.g., Curle & Davies, 1971; Caughey, 1982). Indeed, the model presented here does allow for smooth transitions to moderately supersonic velocities in the heliotail’s flanks, as will be shown in Section 2.6. Moreover, the absence of shocks should not be viewed as too severe a restriction in its applicability to the real heliosphere, the outer flow of which is clearly supersonic but may or may not be super-Alfvénic. The existence of a bow shock depends on the flow being faster than the fastest-propagating signal, which in this case are fast magnetosonic waves. Therefore, the upstream Mach number $m$ that will be introduced later in this section should not so much be identified with the actual sonic Mach number of the incoming interstellar flow, but rather be viewed as a parameter to be chosen close to the fast magnetosonic Mach number, which may well be around or even below unity.

Since, in this new framework, we have so far only fixed the product of $n$ and ${\bf u}$, an additional equation is needed. One of the most straightforward choices would the momentum balance equation

[TABLE]

with $m_{\rm p}$ being the proton mass and $P$ the gas pressure. (We refrain from introducing the otherwise canonical symbol $\rho$ for the density $m_{\rm p}n$ in order to reserve it for the cylindrical radius.) This yields two differential equations for $n$, as the azimuthal component vanishes identically due to rotational symmetry. For our choice of ${\bf s}$ there is no solution for $n$ satisfying both equations simultaneously. Rather than arbitrarily picking one of them, we instead consider the projection of the momentum equation onto stream lines. This is essentially identical to the use of Bernoulli’s equation for compressible, adiabatic flow (e.g., Fahr & Neutsch, 1983). The system of equations is then to be closed with the polytropic equation of state

[TABLE]

in which $\gamma$ is the adiabatic index, and the constant $C$ is determined from a boundary condition at infinity, see Eq. (21).

The projection of Eq. (7) onto stream lines reads

[TABLE]

Substituting the gradient of Eq. (8)

[TABLE]

into the left-hand side of Eq. (2.2) yields

[TABLE]

Since ${\bf s}$ is prescribed, Eq. (11) constitutes a nonlinear first-order partial differential equation for $n$ that will now be solved using the method of characteristics. To this end, we introduce new coordinates $\alpha$, $\beta$ such that

[TABLE]

with an arbitrary function $p=p(\alpha,\beta)$. This condition results in the coupled system of ordinary differential equations (ODEs)

[TABLE]

Choosing $p=s_{\rho}$ yields

[TABLE]

and in particular

[TABLE]

which implies that the characteristics and the stream lines of ${\bf s}$ coincide. We may thus replace the operator ${\bf s}\cdot\nabla$ in Eq. (11) by the derivative $s_{\rho}\,\partial_{\alpha}$ along the stream lines. After multiplication with -m_{\rm p}/\big{(}C\gamma n_{0}^{\gamma+1}s_{\rho}\big{)}, where $n_{0}$ is the number density at upstream infinity, we find

[TABLE]

In terms of the function

[TABLE]

Eq. (16) can be rewritten as

[TABLE]

Using $g$ — instead of $\alpha$ — as a new coordinate, we apply the chain rule to arrive at the ODE

[TABLE]

for the normalized density $\bar{n}:=n/n_{0}$. (For the remainder of this paper, a subscript [math] marks quantities taken at upstream infinity, and a bar denotes normalization with respect to this boundary value, i.e., $\bar{X}=X/X_{0}$ for any quantity $X$.) Introducing the sound speed

[TABLE]

and evaluating this expression at upstream infinity, the denominator of $g$ in (17) becomes

[TABLE]

Furthermore, since ${\bf s}_{0}^{2}=(n_{0}\,{\bf u}_{0})^{2}$, we may write (17) as

[TABLE]

in which the parameter $m:=\|{\bf u}_{0}\|/c_{0}$ denotes the hydrodynamic Mach number at upstream infinity, and the factor $\bar{\bf s}^{2}$ evaluates to

[TABLE]

This quantity attains its minimum value of zero at the stagnation point $(\rho,z)=(0,\sqrt{q})$ and tends to unity both for $r\rightarrow\infty$ and on the surface $2rz=q$, which passes through $(\rho,z)=(0,\sqrt{q/2})$ and approaches the $z=0$ plane for large $r$. The largest value outside the heliopause is reached on this surface at $(\rho,z)=(2\sqrt{6}/3,\,-\sqrt{3}/3)\,\sqrt{q}$ and amounts to $4/3$. In other words, ${\cal A}$ maps the entire region exterior to the heliopause onto the interval $[0,4/3]$.

2.3. Analytical Solutions

In order to solve Eq. (19), we rewrite it as an ODE for the function $g({\bar{n}})$. To this end, we need to establish that ${\bar{n}}(g)$ is invertible. Since, for ${\bar{n}}$ being continuous, this is equivalent to $({\rm d}{\bar{n}}/{\rm d}g)$ never changing sign, we see from Eq. (19) that the square bracket must not pass through zero. Since it evaluates to

[TABLE]

such a change of sign could only occur at a sonic transition. This implies that our model’s momentum-conserving flow solution is restricted to the purely subsonic case ${\bf u}^{2}<c^{2}\ \Rightarrow\ m<1$.

Having established that the inverse ${\bar{n}}^{-1}=g$ exists, we multiply Eq. (19) by $({\rm d}g/{\rm d}{\bar{n}})$ and again apply the chain rule, which yields

[TABLE]

The homogeneous solution, i.e., the solution of the ODE

[TABLE]

is simply obtained by integration with respect to ${\bar{n}}$ as

[TABLE]

where $c_{1}\in\mathbb{R}$ is a constant. Accordingly, the inhomogeneous solution can be derived using the ansatz $g_{\rm inhom}=h({\bar{n}})\,{\bar{n}}^{2}$, leading to the ODE

[TABLE]

for $h$, which can also be solved by simple integration with respect to ${\bar{n}}$. The inhomogeneous solution $g_{\rm inhom}$ is then

[TABLE]

with a constant $c_{2}\in\mathbb{R}$. Since the full solution is a superposition of the homogeneous and the inhomogeneous solution, and the constant $c_{1}$ can be absorbed into the constant $c_{2}$, it follows that $g$ itself is of the form (29). Finally, the density can be determined by “solving” the transcendental equations in (29) for ${\bar{n}}$. This is done in the following two subsections, considering separately the isothermal ($\gamma=1$) and the adiabatic ($\gamma=5/3$) cases.

2.4. Isothermal Flow

Inserting the ansatz

[TABLE]

into Eq. (29) for $\gamma=1$ results in

[TABLE]

Using the abbreviations

[TABLE]

this can also be written as

[TABLE]

which is the defining equation for Lambert’s W function. Thus, we obtain

[TABLE]

for the density. The range of density values to be covered stipulates that the principal branch of $W$ be employed. Furthermore, the constant $c_{2}$ is fixed by the boundary values

[TABLE]

according to Eqs. (22) and (29). It then directly follows that $c_{2}=-m^{2}$. In cylindrical coordinates, the density finally becomes

[TABLE]

with ${\cal A}$ defined in (23).

2.5. Mono-atomic Ideal Gas Flow

For $\gamma\in\mathbb{N}\backslash\{1\}$, Eq. (29) gives the zeros of a polynomial in ${\bar{n}}$ with powers of 0, 2, and $\gamma+1$. Analytical inversion is therefore only possible for $\gamma\in\{2,3\}$, and none of these cases is particularly meaningful in the heliospheric context. Moreover, since we are mainly interested in the case of $\gamma=5/3$ describing mono-atomic ideal gases, we substitute this value into the adiabatic flow equation (cf. the second case of Eq. (29)), obtaining

[TABLE]

In analogy to (36), the constant $c_{2}$ is found from

[TABLE]

The density is then implicitly given by

[TABLE]

with ${\bar{n}}$ being the solution to

[TABLE]

In a cylindrical representation, the density is again obtained by replacing $g$ with $-m^{2}\,{\cal A}$.

2.6. Density Approximations

The formulas for the density profiles

[TABLE]

for $\gamma=1$ (cf. Eq. (37)) and

[TABLE]

for $\gamma=5/3$ (cf. Eq. (41)) are somewhat cumbersome to handle due to the involvement of Lambert’s W function, which is only implicitly defined by the transcendental equation (34), and the fact that Eq. (41) leads to a fourth-order polynomial in ${\bar{n}}^{2}$. For this reason, we consider in the following suitable analytical approximations for ${\bar{n}}({\cal A})$ and the corresponding approximate flow fields. These open up the possibility to find exact analytical magnetic field solutions of the steady-state induction equation (4) and the magnetic divergence constraint (5) with respect to these approximative flows.

Both density profiles (42) and (43) can be reasonably well approximated by linear functions ${\bar{n}}^{(1)}({\cal A})$ that pass through ${\bar{n}}(1)=1$ (thus ensuring $n\rightarrow n_{0}$ at infinity, where we prescribe the boundary conditions) and the maximum value ${\bar{n}}(0)=:{\bar{n}}_{\rm sp}$ reached at the stagnation point. Although in principle any other values of ${\cal A}\in[0,4/3]$ could be used, values close to $4/3$ are only approached for large, finite distances in $(\rho,z)$ space and, therefore, of minor relevance for our intended application. Using second-order polynomials ${\bar{n}}^{(2)}({\cal A})$, the additional degree of freedom may be fixed by requiring the correct derivative $D_{{\cal A}_{\rm C}}:=\partial_{{\cal A}}{\bar{n}}|_{{\cal A}_{\rm C}}$ at either point ${\cal A}_{\rm C}\in\{0,1\}$ (or, alternatively, the correct density at some intermediate value of ${\cal A}$). We settle for a correct derivative at the stagnation point ${\cal A}=0$. In this case, suitable approximations for ${\bar{n}}({\cal A})$ to first and second order are

[TABLE]

with

[TABLE]

and

[TABLE]

for which we have used the fact that Eq. (19) at ${\cal A}=0$ can be transformed into

[TABLE]

and evaluated separately for both values of $\gamma$. For instance, when choosing the specific Mach number $m=0.6$, we obtain

[TABLE]

as suitable approximations for the respective density profiles.

As can be deduced from (46), the peak density ${\bar{n}}_{\rm sp}$ is always slightly higher in the isothermal case (by a factor of about $(1+m^{4}/14)$, i.e., at most $\sim{}7$%). This is plausible, since in this case the pressure gradient that decelerates the incoming flow ahead of the stagnation point is not caused by the combined gradients of density and temperature, but has to come from density alone.

Fig. 1 compares the exact isothermal solution for various values of $m$ against its first- and second-order approximations. As expected, the agreement is most favorable for small $m$, but continues to be useful also for larger values. The corresponding plot for $\gamma=5/3$ is not shown since it would look very similar to its isothermal counterpart.

It should be noted that, since Lambert’s W function can, by definition, only accommodate real-valued arguments larger than or equal to $-1/e$ (corresponding to subsonic flow), the density (42) ceases to be well-defined in regions for which

[TABLE]

However, since ${\cal A}\leq 4/3$ holds everywhere outside the heliopause, such regions can only arise for Mach numbers larger than $m_{\rm crit,isoth}:=\sqrt{-W(-3/(4\,e))}\approx 0.648$ (see, for instance, the solid curves for $m=0.75$ and $m=1.0$). For $\gamma=5/3$, a corresponding critical Mach number is found at $m_{\rm crit,adiab}\approx 0.619$.

As can prominently be seen from the right plot, the approximation is not at all restricted to subsonic flow, but does indeed feature a smooth sonic transition, leading to moderately high local Mach numbers of up to about 1.83 before eventually slowing down again toward its original, subsonic speed.

Fig. 2 shows the density contours in $(\rho,z)$ space that result from the exact solution and its second-order approximation, confirming the expectation of mass piling up in front of the heliopause and being washed down the flanks. The chosen value of $m=0.6$ is close to $m_{\rm crit,isoth}$, beyond which the exact solution is no longer well-defined everywhere outside the heliopause. However, we remark that Mach numbers up to $m=1$ may be used for the approximation in the entire space and the exact solution restricted to the upwind half-space near the inflow axis. As expected, the density is highest at the stagnation point. Inside the white, approximately circular region around the origin, no solution is available. This region grows with larger $m$, eventually engulfing not only the heliopause but the entire downwind half-space $z\leq 0$ as $m\rightarrow 1$. Since the approximative solution is a simple second-order polynomial, it can be defined in the entire space for any $m$, and thus exhibits no such white region.

It may at first sight seem questionable to rely on an “approximation” in a region where the exact solution is not at all available. What happens is that the approximation is first derived from the exact solution where it exists, and subsequently continued into the “white” region in which the exact solution is not defined, exploiting the fact that the former, due to its simple polynomial form, is void of any domain restrictions. As a consequence, we have no measure of the flow’s departure from the exact momentum-conserving solution in the continuation region. What we do know, however, is that it proceeds along the same Rankine-type stream lines, conserves mass exactly, and connects smoothly into the original solution’s domain. Given that momentum is not exactly conserved in either region anyway, we feel that these characteristics are nevertheless sufficient for the flow field to clearly qualify as physically meaningful for its intended purpose of deriving an improved formula for the heliospheric interstellar magnetic field. In the same vein, we note that none of the existing models (Whang, 2010; Schwadron et al., 2014; Isenberg et al., 2015, etc.) come anywhere near this degree of physical realism, despite their unimpeachably acknowledged usefulness for the heliophysics community. In the following, we will continue to refer to the thus extended approximation simply as “approximation,” irrespectively of the position at which it is evaluated.

The inflow axis ($\rho=0$) is of particular interest, especially since the full momentum equation is satisfied there. Therefore, Fig. 3 depicts the density and flow velocity along this axis as functions of the upstream heliocentric distance $z$ for both the isothermal ($\gamma=1$) and the adiabatic ($\gamma=5/3$) cases. The difference between them is obviously negligible except for large $m$. A pile-up of mass is apparent, which is of course absent from the incompressible case ($m=0$). The flow also decelerates much stronger at high Mach numbers. Additionally shown are the corresponding profiles obtained from a fully self-consistent hydrodynamical simulation, details of which are described in Section 4.1. The normalized numerical peak density of 1.338 thus obtained is in excellent agreement with its theoretical value of 1.337. Moreover, the value of $\bar{n}=1.078$ at $z=5\sqrt{q}$ is only slightly above 1.056, which is the normalized density predicted by Eq. (43).

3. Derivation of the Local Interstellar

Magnetic Field

3.1. Cylindrical Field Components

The components of the advected magnetic field ${\bf B}$ are derived using the method of line conservation (see Elsasser, 1956; Naor & Keshet, 2015), which applies the idea behind Cauchy’s integral (first formulated by Cauchy (1816) in the context of fluid mechanics) to ideal MHD. We thus exploit the fact that in ideal MHD, the quantity ${\bf B}/n$ satisfies the same equation of motion as the line element ${\bf\delta x}$ connecting two neighboring particles that are passively advected in the flow field ${\bf u}$.

While at upstream infinity, $\left({\bf B}/n\right)|_{\infty}={\bf B}_{0}/n_{0}$ is spanned by the basis vectors $\{{\bf e}_{\rho},{\bf e}_{\varphi},{\bf e}_{z}\}$, the corresponding basis at finite position is $\{{\bf c},(\rho/a)\,{\bf e}_{\varphi},-\bar{\bf u}\}$, where ${\bf c}$ is a vector connecting two fluid elements that start at the same “height” $z_{0}$ on adjacent stream lines labelled by $a$ and $a+\delta a$, and $a$ is a known function of position (cf. Fig. 2 and Eq. (49) in Paper I). Since the basis vectors are co-moving with the flow, the coefficients of ${\bf B}/n$ with respect to both bases are the same. We are therefore led to

[TABLE]

and thereby arrive at the relation

[TABLE]

that generalizes Eq. (67) of Paper I to the compressible case. Furthermore, it was shown in Paper I that the condition of equal travel times

[TABLE]

can be used to derive the components of ${\bf c}$ as

[TABLE]

The function $z_{a}(\rho)$ in Eq. (55) designates the $z$ coordinate of the flow line labelled by $a$ and parameterized by $\rho$, which passes through $(\rho_{a},z_{a}(\rho_{a}))$. Its derivative evaluates to

[TABLE]

according to Eq. (B.8) in Paper I.

We now relate the magnetic field of the incompressible case (as derived in and known from Paper I) to the compressible case (with $m>0$). For this purpose, we label the respective quantities with superscripts “I” and “C,” such that

[TABLE]

where ${\bar{n}}^{\rm C}={\bar{n}}$, ${\bar{n}}^{\rm I}=1$, and

[TABLE]

In both cases, the momentum density is prescribed through the same potential as ${\bf s}^{\rm I,C}=-\nabla\Phi$, implying

[TABLE]

Since ${\bar{n}}$ (and thus $\bar{\bf u}^{\rm C}$) is given only implicitly, we may not hope to obtain analytical expressions for the integral $K_{0}^{\rm C}$ in Eqs. (57) and (3.1). However, when approximating the normalized density as

[TABLE]

with the values for $\nu_{1,2}$ read off from Eq. (45), $K_{0}^{\rm C}$ can be related to $K_{0}^{\rm I}$ via

[TABLE]

in which the definition

[TABLE]

generalizes that in (59). This allows us to express the desired magnetic field solution for the compressible case in terms of the known one for the incompressible case as

[TABLE]

where

[TABLE]

The incompressible case is evidently recovered in the limit $m\rightarrow 0$, as in that case ${\bar{n}}_{\rm sp}\rightarrow 1$, ${\bar{n}}\rightarrow 1$, ${\cal C}\rightarrow 0$, and ${\cal Z}\rightarrow 0$.

We now determine an explicit expression for the correction term ${\cal C}$. In order to evaluate the integrals $K_{1,2}$, we first consider the normalized flow potential $\bar{\Phi}=\bar{\Phi}(\rho,z_{a}(\rho))$ along the streamline in question. The total $\rho$ derivative of this quantity becomes

[TABLE]

For $k=1$, we thus simply obtain

[TABLE]

where

[TABLE]

The integral for $k=2$ requires explicit evaluation. Using Eqs. (68) and (3.1) yields, after lengthy but straightforward computation,

[TABLE]

At this point, all prerequisites needed for the explicit construction of the ’compressible’ magnetic field components $B^{\rm C}_{\rho,\varphi,z}$ from the known ’incompressible’ ones $B^{\rm I}_{\rho,\varphi,z}$ are in place, and this construction proceeds as follows.

Choose an upstream Mach number $m<1$ (preferably near unity). The fast magnetosonic Mach number might serve as an educated guess (Spreiter & Stahara, 1995). 2. 2.

Compute the constant coefficients ${\bar{n}}_{\rm sp}$ and $\nu_{1,2}$ using Eqs. (46) and either (44) or (45). 3. 3.

At each desired position $(\rho,\varphi,z)$, evaluate

• (a)

functions $K_{1,2}$ from Eqs. (69) and (71),

• (b)

functions ${\cal C}$ and ${\cal Z}$ from Eqs. (3.1) and (67),

• (c)

function $\partial_{a}z_{a}(\rho)$ from Eq. (56). 4. 4.

Finally, substitute all quantities into Eqs. (64)–(66) to arrive at the desired field components, bearing in mind that ${\bf u}^{\rm I}=-\nabla\Phi$ according to Eq. (2).

In principle, the approximation (61) could be continued to even higher orders in ${\cal A}$. However, tentative computations indicate that the expressions for $K_{k}^{\rm I}$ become much more involved for $k>2$, and that comparatively little could be gained by going to order ${\cal A}^{3}$ or higher. We concede that the presented procedure is already more involved than the original field derived in Paper I, but also note that the result is clearly more realistic, not only because it satisfies more physical constraints, but also because it does indeed perform better in direct comparison to a numerical model. This will be demonstrated in Section 4.

3.2. Cartesian Field Components on the Inflow Axis

The cylindrical components (64)–(66) of the magnetic field can easily be converted into a Cartesian representation. The only region for which this conversion is not straightforward is the $z$ axis, where both $B_{\rho,\varphi}$ and the transformation factors $\sin\varphi$ and $\cos\varphi$ are ill-defined, and the limit $\rho\rightarrow 0$ is somewhat cumbersome to evaluate. But since the conditions along the inflow axis are of particular physical interest (especially given that they satisfy the momentum conservation equation (7) exactly), it is fortunate that an easier avenue is available to treat this special case. For this purpose, we revert to the induction equation (4), now written in the form

[TABLE]

which is equivalent to

[TABLE]

due to $\nabla\cdot{\bf s}=0=\nabla\cdot{\bf B}$. On the $z$ axis, where the only non-zero components of both ${\bf s}$ and $\nabla n$ are those in $z$ direction, the left-hand side of Eq. (73) simplifies to

[TABLE]

The Cartesian components of the right-hand side evaluate to

[TABLE]

In total, we obtain the differential equations

[TABLE]

which are readily solved to yield the remarkably simple expressions

[TABLE]

as a generalization of the corresponding formulas from Appendix C of Paper I.

It is interesting to note that, while $B_{x}$ and $B_{y}$ each contribute a factor ${\bar{n}}$ in the compressible case, $B_{z}$ is identical to its incompressible analog. This observation may at first sight come as a surprise, but is indeed consistent with the notion of line conservation: For the distance $\delta z$ between two neighboring particles both travelling along the inflow axis, we have $B_{z}/n\propto\delta z\propto u_{z}$ and hence $B_{z}\propto n\,u_{z}=s_{z}=-\partial_{z}\Phi$, which is the same in both cases. On the other hand, the separation $\delta x$ between two particles travelling on adjacent flow lines at the same height $z$ only depends on the geometry of the flow lines at that height, which is determined by $s_{z}$, hence $B_{x}/n\propto\delta x\propto s_{x}$, so $B_{x}\propto n\,s_{x}$.

4. Comparison to (M)HD Simulations

Whenever the need for a prescription of an interstellar magnetic field arises in the heliospheric context, it is obtained either analytically or by way of numerical simulations. Therefore, it seems reasonable to asses the usefulness of our new model by comparing it to a field that would typically be generated using a self-consistent MHD code. In the same vein, it is of interest to investigate to what extent the hydrodynamical aspects of this model are compatible with corresponding simulation data.

4.1. Pure HD for Flow Structure

We first provide details of the simulation underlying the gray plot profiles in the panel of Fig. 3. For these, the grid-based MHD code Cronos was employed to simulate an unmagnetized, axisymmetric heliosphere. The poloidal $(\rho,z)$ domain of size $[0,2000]$ AU$\,\times\,[-1000,1000]$ AU was covered with cells of size $\Delta\rho=1$ AU for $\rho<1000$ AU and 5 AU for $\rho\in[1000,2000]$ AU, and $\Delta z=1$ AU throughout. The outer layer of coarser cells served to shield the inner part from reflections possibly emanating from this boundary. Parameters were those from the plasma-only version of the Müller et al. (2008) benchmark (cf. Table 1 in that paper), except that the ISM temperature was raised from 6530 K to 39553 K in order to achieve the desired upstream Mach number of $m=0.8$. Such high temperature is consistent with numerical findings (see, e.g., Pauls & Zank, 1997; Fahr et al., 2000) that take explicitly into account a shock transition and charge exchange processes. The displayed profiles were extracted at simulation time $t=1500$ (in units of 1 AU$/c_{\rm s}$, corresponding to about 280 years in physical units), at which the configuration was deemed sufficiently stationary by visual inspection.

4.2. Full MHD for Magnetic Structure

The simulation that we employ to assess the quality of our improved magnetic field solution, as well as the resulting data, are identical to those used in Paper I. In particular, the undisturbed interstellar magnetic field of strength 0.3 nT is oriented in the $x$–$z$ plane and makes an angle of $50^{\circ}$ with the inflow axis. For further details of the numerical setup, the reader is again referred to that paper.

Fig. 4 compares all magnetic field components of the numerical solution to its isothermal, second-order, analytical counterpart along straight lines parallel to the Cartesian $x$, $y$, and $z$ axes and passing close to the stagnation point. The superior performance of our improved compressible model over its incompressible predecessor is evident from the observation that in all cases, the red curves yield much better approximations to the numerical data (black) than the blue ones. Together with the field line plot in the lower left quadrant, this panel illustrates and confirms the expected differences between the incompressible and the compressible model: In the compressible case, the heliosphere’s region of influence extends much further out in the upstream direction, the field lines are less sharply curved, and the field strength amplification resulting from magnetic pile-up ahead of the stagnation point is more pronounced. All of these characteristics contribute to a markedly improved degree of realism and physical usefulness.

5. Summary and Conclusions

In this work, the exact analytical MHD model by Röken et al. (2015) for the interstellar magnetic field in the outer vicinity of the heliopause has been significantly extended from incompressible to compressible, yet predominantly subsonic flow, with the upstream sonic Mach number $m<1$ as a free parameter. Although the use of flow potentials precludes the occurrence of a bow shock, the model’s qualitative applicability may well be extended to supersonic but sub-alfvénic flows, since it should be of minor relevance whether the field lines are distorted due to thermal or magnetic pressure.

This improved model transcends its predecessor not only conceptually by approximately satisfying conservation of linear momentum along stream lines, but also yields results that are much closer to those from a self-consistent, grid-based MHD simulation.

While the derived spatial structures of flow lines and density profiles are exact for both isothermal and adiabatic settings, the corresponding solution for a magnetic field being passively frozen into this flow required the exact density solution to be replaced by a suitable polynomial approximation to warrant analytical tractability. As an additional benefit, this approximation allows itself to be continued into spatial regions in which the exact solution is not defined. And unlike the exact flow solution, its continued approximation is not restricted to subsonic flows, but features a smooth sonic transition (reminiscent of “sonic lines,” e.g., Scherer et al., 2016) in the heliotail’s flanks.

Furthermore, for the velocity field associated with this approximative density distribution, the derived magnetic field is again an exact, solenoidal solution to the induction equation of ideal MHD. In particular, the magnetic field’s orientation at the upstream boundary can still be chosen freely. But even disregarding the entire magnetic part of the model, its hydrodynamic configuration can still offer a marked improvement over models in which incompressible or even compressible but isothermal flows are prescribed (Parker, 1961; Nerney & Suess, 1995). We note in passing that, although this option was not exploited in the present work, a physically meaningful temperature distribution of the local, heliosphere-dominated interstellar medium becomes readily accessible from the ideal gas law and the density field, which our model can provide.

We acknowledge financial support via the projects FI 706/15-1 and FI 706/21-1 funded by the Deutsche Forschungsgemeinschaft (DFG), as well as through the RAPP (Ruhr Astroparticle and Plasma Physics) Center, funded as MERCUR project St-2014-040. Furthermore, we also thank the anonymous referee for helpful and constructive comments.