# An exact analytical solution for the weakly magnetized flow around an
axially symmetric paraboloid, with application to magnetosphere models

Jens Kleimann Theoretische Physik IV, Ruhr-Universität Bochum, 44780 Bochum,
Germany  Ruhr Astroparticle and Plasma Physics Center (RAPP Center), Germany
e-mail<EMAIL_ADDRESS>Christian Röken Department of Geometry and Topology,
Faculty of Science, University of Granada, 18071 Granada, Spain  Lichtenberg
Group for History and Philosophy of Physics, Institut für Philosophie,
Universität Bonn, 53115 Bonn, Germany

###### Abstract

Rotationally symmetric bodies with longitudinal cross sections of parabolic
shape are frequently used to model astrophysical objects, such as
magnetospheres and other blunt objects, immersed in interplanetary or
interstellar gas or plasma flows. We discuss a simple formula for the
potential flow of an incompressible fluid around an elliptic paraboloid whose
axis of symmetry coincides with the direction of incoming flow. Prescribing
this flow, we derive an exact analytical solution to the induction equation of
ideal magnetohydrodynamics for the case of an initially homogeneous magnetic
field of arbitrary orientation being passively advected in this flow. Our
solution procedure employs Euler potentials and Cauchy’s integral formalism
based on the flow’s stream function and isochrones. Furthermore, we use a
particular renormalization procedure that allows us to generate more general
analytic expressions modeling the deformations experienced by arbitrary scalar
or vector-valued fields embedded in the flow as they are advected first
towards and then past the parabolic obstacle. Finally, both the velocity field
and the magnetic field embedded therein are generalized from incompressible to
mildly compressible flow, where the associated density distribution is found
from Bernoulli’s principle.

##  1 Introduction

The study of fluid flow around solid obstacles is a classical topic of
hydrodynamics, with countless applications in both engineering and physics,
and especially astrophysics. Therein, the case of the fluid being magnetized
is of particular interest for several classes of objects, ranging from the
heliosphere (e.g. Röken et al., 2015; Kleimann et al., 2022), planets (e.g.
Petrinec & Russell, 1997; Kotova et al., 2020) and their magnetospheres (e.g.
Spreiter & Stahara, 1995; Richardson, 2002), comets and planetary satellites
(e.g. Combi et al., 2002, and references therein), up to clusters of galaxies
(e.g. Dursi & Pfrommer, 2008).

While quantitative applications are often faced with complications such as
turbulence, reactive flows, or high Mach numbers that typically can only be
investigated using numerical simulations, a limited set of analytical models
has been obtained by restriction to stationary flows that exhibit some form of
spatial symmetry. For instance, Nabert et al. (2013) employ certain axial and
point-like symmetries of the magnetohydrodynamic (MHD) quantities and a
stationary flow to derive the plasma properties of the Earth’s magnetosheath,
the region enclosed between the magnetopause and the bow shock. In a similar
vein, the present paper considers stationary flows $\vec{u}$ which satisfy the
constraints of being both irrotational

$\nabla\times\vec{u}=\vec{0}$ (1)

as well as solenoidal

$\nabla\cdot\vec{u}=0\,.$ (2)

The popularity of these two constraints is mainly due to the fact that they
allow $\vec{u}$ to be expressed through a scalar potential $\Phi$ via
$\vec{u}=-\nabla\Phi$, with $\Phi$ satisfying the Laplace equation

$\nabla^{2}\Phi=\nabla\cdot(\nabla\Phi)=\nabla\cdot(-\vec{u})=0\,.$ (3)

This implies in particular that the flow field is fully determined by the
shape of the obstacle and the flow conditions prevailing at the upstream
boundary. In the case of planar symmetry, in which the shapes and flows are
invariant along a Cartesian direction, classical examples of obstacles for
which such potential flows have been explicitly derived include infinite
planes, wedges, and the circular cylinder. By introducing a complex velocity
potential, it is possible to generate additional solutions from the existing
ones through the application of conformal mappings such as the Joukowski
transform (e.g. Batchelor, 2000, pp. 435).

This paper, however, is concerned with settings where both the solid object
and the flow around it feature an axial symmetry with respect to the inflow
direction, rather than a translational one. For these, the method of complex
potentials and conformal mappings is not available, and thus the set of known
examples is correspondingly smaller.

The set of known nontrivial solutions for magnetized generalizations of such
flows is smaller still, having until now essentially been restricted to the
sphere (e.g. Dursi & Pfrommer, 2008; Isenberg et al., 2015, and references
therein) and the Rankine-type heliosphere (Röken et al., 2015; Isenberg et
al., 2015). In the former case of a sphere, the exact magnetic field solution
is known only up to an integral which so far could not be expressed in terms
of analytic functions, whereas the latter Rankine-type case involves special
functions, namely incomplete elliptic integrals. In contrast to these, the MHD
solutions to be presented here are predominately of rather compact form and
expressible through basic functions without recourse to unevaluated integrals,
making them a viable alternative to numerical computations at least in simple
applications.

Moreover, besides the obvious puristic and educational motivations, this
geometrical setting is also of significant astrophysical interest. Starting
with the work of Kobel & Flückiger (1994) and its subsequent refinements (e.g.
Romashets & Vandas, 2019, and references therein), the upwind sides of
planetary magnetopauses continue to be frequently modeled as rotationally
symmetric paraboloids. Although the present approach cannot include
discontinuous surfaces like bow shocks because it computes its velocity
components from a smooth potential, it is certainly relevant at least in
situations where a bow shock can either be neglected or is absent altogether.

The paper is structured as follows. After this introduction, the analytic
expressions and properties of the potential flow field around an axially-
symmetric paraboloid are derived and discussed in Section 2. Thereafter,
Section 3 presents the analytic derivation and discussion of a magnetic field
$\vec{B}$ being passively advected in this flow by solving the stationary
induction equation of ideal MHD

$\nabla\times(\vec{u}\times\vec{B})=\vec{0}\,.$ (4)

This includes the extension to paraboloids being linearly scaled to a custom
aspect ratio. Section 4 presents further generalizations, most notably to the
advection of more general scalar and vectorial fields. Finally, Section 5
generalizes the ambient paraboloid flow field to subsonic compressible flow,
including the associated distribution of mass density, and derives the
correspondingly changed magnetic field, after which Section 6 concludes the
paper with a summary.

##  2 The inviscid flow field

In order to motivate the course of the present analysis, we start by
considering the example of inviscid flow around a sphere, whose well-known
flow potential and stream function (e.g. Batchelor, 2000, pp. 452) in terms of
cylindrical coordinates $(\rho,\varphi,z)\in\mathbb{R}_{\geq
0}\times[0,2\pi)\times\mathbb{R}$ read, respectively,

$\displaystyle\Phi_{0}$ $\displaystyle=$ $\displaystyle
z\left(1+\frac{1}{2\,r^{3}}\right)$ (5) $\displaystyle\Psi_{0}$
$\displaystyle=$
$\displaystyle-\frac{\rho^{2}}{2}\left(1-\frac{1}{r^{3}}\right),$ (6)

and result in the flow field

$\displaystyle\vec{u}_{0}$ $\displaystyle=$
$\displaystyle-\nabla\Phi_{0}=\nabla\times\left[(\Psi_{0}/\rho)\,\vec{e}_{\varphi}\right]$
(7) $\displaystyle=$
$\displaystyle\frac{3\,\rho\,z}{2\,r^{5}}\,\vec{e}_{\rho}+\left(\frac{z^{2}-\rho^{2}/2}{r^{5}}-1\right)\vec{e}_{z}\,,$

where we use the shorthand $r=\sqrt{\rho^{2}+z^{2}}$ wherever it allows for
more compact notation of formulas. Throughout the paper, quantities of lengths
are normalized to a characteristic length scale (such as the sphere’s radius
in the above case), while all other quantities, most notably the magnitudes of
the velocity and the magnetic field, are normalized to their respective values
at the upstream boundary.

Since streamlines are given by lines of constant $\Psi_{0}$, we may interpret
Equation (6) as a special case of the simpler stream function
$\Psi_{\mathrm{bg}}:=-\rho^{2}/2$ of the undisturbed background flow
$\vec{u}_{\mathrm{bg}}=-\vec{e}_{z}$ being modulated by a factor
$1-1/f(\rho,z)$ which characterizes the shape of the obstacle in such a way
that $f(\rho,z)=1$ represents its contour, in this case the sphere with $r=1$.
The particular streamline $\Psi_{0}=0$ will then include not only the entire
$z$-axis but split at two stagnation points (or one in cases where the
obstacle is of infinite extent in the downstream direction), from which the
emanating streamlines precisely trace the obstacle’s surface.

This observation invites the question whether the flow around other shapes
with axial symmetry might be obtained in an analogous manner. For instance,
one might attempt to generalize the sphere flow with its circular contour
$f(\rho,z)=r^{3}$ to the more general contour

$f(\rho,z)=g_{\varepsilon}(r+\varepsilon\,z)$ (8)

with $g_{0}(r)=r^{3}$ and the argument encoding the parametric equation

$r+\varepsilon\,z=1\quad\Leftrightarrow\quad
r=\frac{1}{1+\varepsilon\cos\vartheta}$ (9)

of a conic section with its focal point at the origin, where
$z=r\cos\vartheta$ and $\vartheta\in(0,\pi)$, covering ellipses
($0<\varepsilon<1$) and hyperbolas ($\varepsilon>1$) of eccentricity
$\varepsilon$, as well as the special cases of circle ($\varepsilon=0$) and
parabola ($\varepsilon=1$). But while the resulting flows around such
ellipsoidal or hyperboloidal bodies of revolution are indeed solenoidal by
construction (and might therefore be of actual interest for specialized
applications), direct calculations show that no $g_{\varepsilon}(\cdot)$ that
additionally satisfy Equation (1) exists for such bodies (see Appendix A). On
the other hand, similar considerations can be used to generalize the sphere’s
flow potential (5) to nonzero eccentricities, resulting in flow fields
satisfying Equation (1) but not (2). We note that the correct components for
flow around ellipsoids satisfying both constraints involve the solution of the
Laplace equation $\nabla^{2}\Phi_{\varepsilon}=0$ in ellipsoidal coordinates
(e.g. Milne-Thomson, 1968, pp. 475), and the rather complicated form of the
associated expressions renders the general case $0<\varepsilon<1$ less
attractive in the present context.

In the following, we will therefore focus on the parabolic case
$\varepsilon=1$ with $g_{1}(r+z)=r+z$ being the identity mapping. This leads
to stream function

$\Psi=-\frac{\rho^{2}}{2}\left(1-\frac{1}{r+z}\right)=-\frac{1}{2}\,(r-z)\,(r+z-1)\,,$
(10)

where the parabolic obstacle itself is given by

$r+z=1\quad\Leftrightarrow\quad 2\,z=1-\rho^{2}\ ,$ (11)

and the lower index ’1’, which was previously used to hint at $\varepsilon=1$,
will be suppressed from now on. Equation (10) reproduces the stream function
formula

$\psi=2\,c^{2}U\left(\xi^{2}-\xi_{0}^{2}\right)\eta^{2}$ (12)

from Milne-Thomson (1968, p. 479), expressed in elliptical coordinates
$(\xi,\eta)$ satisfying
$(\rho,z)=\left(2\,c\,\xi\,\eta,[\xi^{2}-\eta^{2}]\,c\right)$ when choosing
$c=U=1$ and $\xi_{0}=1/\sqrt{2}$, up to an overall minus sign (i.e.
$\psi=-\Psi$) that originates from an alternative definition of the stream
function used therein. Applying the relation

$\vec{u}=\nabla\times\left[(\Psi/\rho)\,\vec{e}_{\varphi}\right]$ (13)

to the stream function (10), we obtain flow components

$\displaystyle u_{\rho}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho}\frac{\partial\Psi}{\partial
z}=\frac{1}{2\,r}\,\sqrt{\frac{r-z}{r+z}}$ (14) $\displaystyle u_{z}$
$\displaystyle=$
$\displaystyle+\frac{1}{\rho}\frac{\partial\Psi}{\partial\rho}=\frac{1}{2\,r}-1$
(15)

which, by way of direct calculation, can indeed be shown to be irrotational
(see again Appendix A). Therefore, a flow potential for this field exists, and
is easily found to be of the simple form

$\Phi=z-\frac{\ln(r+z)}{2}.$ (16)

Figure 1 visualizes the structure of the flow field (14)–(15) and the
associated stream function (10) and potential (16). As can be seen from both
the plots and Equation (10), the parabolic surface comes about as the
separatrix between the incident background flow and that of a line source
extending along the entire negative $z$-axis, in contrast to the sphere case,
where the flow is created by a point dipole at the origin. As expected, the
stagnation point is located at $z=1/2$. Along the stagnation line
($\rho=0,z>1/2$), the flow decelerates rather gently according to
$u_{z}|_{\rho=0}=(2\,z)^{-1}-1$, unlike the flow ahead of a sphere, for which
Equation (7) yields a somewhat sharper profile, namely
$u_{0,z}|_{\rho=0}=z^{-3}-1$. This difference is likely caused by the
paraboloid’s infinite lateral extent, which requires a much larger body of
fluid to be displaced.

The flow field (14)–(15) is of course markedly simpler than the corresponding
viscous flow that has been investigated analytically by Miller (1969, 1971)
and numerically by, e.g. Davis & Werle (1972) because it avoids the intricate
modeling of boundary layers that is required when solving the Navier–Stokes
equation for the geometry in question.

Figure 1:  Left: streamlines (lines of constant $\Psi$, blue) and
equipotentials (lines of constant $\Phi$, red) at 0.5 spacing, with the
paraboloid obstacle (green) given by its parametric form (11). Right: the
resulting vector field represented by arrows and contours of absolute speed
$\|\vec{u}\|$. The dashed black line marks the flow-generating, semi-infinite
line source. Arrow lengths have been scaled $\propto\|\vec{u}\|^{0.3}$ in
order to accommodate high speeds near the source, and contours above unity
have been suppressed since these occur exclusively inside the obstacle.

##  3 The frozen-in magnetic field

###  3.1 The longitudinal and transversal components

Next, we turn to the computation of the components of a magnetic field which
is embedded into the flow (14)–(15) and passively advected in it. As often in
such cases, we require the magnetic field at upstream infinity
($z\rightarrow+\infty$) to be of constant (though arbitrary) orientation and
magnitude. The boundary field is best specified through its constant Cartesian
components $(B_{x0},B_{y0},B_{z0})$, yielding the cylindrical representation

$\displaystyle B_{\rho 0}:=B_{\rho}|_{z=\infty}$ $\displaystyle=$
$\displaystyle+B_{x0}\cos\varphi+B_{y0}\sin\varphi$ (17) $\displaystyle
B_{\varphi 0}:=B_{\varphi}|_{z=\infty}$ $\displaystyle=$ $\displaystyle-
B_{x0}\sin\varphi+B_{y0}\cos\varphi$ (18) $\displaystyle B_{z}|_{z=\infty}$
$\displaystyle=$ $\displaystyle B_{z0}.$ (19)

To derive the full frozen-in magnetic field, we closely follow the procedure
used by Isenberg et al. (2015) for the Rankine half-body flow by employing the
decomposition

$\vec{B}=\vec{B}_{\mathrm{lo}}+\vec{B}_{\mathrm{tr}}$ (20)

of the magnetic field into a longitudinal (flow-parallel) component
$\vec{B}_{\mathrm{lo}}$ and a transversal component $\vec{B}_{\mathrm{tr}}$,
which at the upstream boundary (but generally not at other locations) is
perpendicular to the flow. When introducing the azimuthal angle $\varphi_{0}$
of the boundary field, i.e. the unique angle satisfying

$\frac{B_{x0}}{\cos\varphi_{0}}=\frac{B_{y0}}{\sin\varphi_{0}}=\sqrt{(B_{x0})^{2}+(B_{y0})^{2}}=B_{\mathrm{tr,0}},$
(21)

the two horizontal components (17) and (18) may alternatively be written as

$\displaystyle B_{\rho 0}$ $\displaystyle=$
$\displaystyle+B_{\mathrm{tr},0}\cos(\varphi-\varphi_{0})$ (22) $\displaystyle
B_{\varphi 0}$ $\displaystyle=$ $\displaystyle-
B_{\mathrm{tr},0}\sin(\varphi-\varphi_{0})\,.$ (23)

Now, being anti-parallel to $\vec{u}$, which at upstream infinity is just
$\vec{u}|_{z=\infty}=-\vec{e}_{z}$, the longitudinal component satisfying the
boundary condition (19) is trivially found to be

$\displaystyle\vec{B}_{\mathrm{lo}}$ $\displaystyle=$ $\displaystyle
B_{z0}\,[-\vec{u}]$ $\displaystyle=$ $\displaystyle
B_{z0}\left[-\frac{1}{2\,r}\,\sqrt{\frac{r-z}{r+z}}\,\vec{e}_{\rho}+\left(1-\frac{1}{2\,r}\right)\vec{e}_{z}\right].$

And just like any solenoidal magnetic field, the transverse component can be
written as

$\vec{B}_{\mathrm{tr}}=\nabla\alpha\times\nabla\beta$ (25)

using Euler potentials $\alpha$ and $\beta$ (e.g. Stern, 1966), which are
scalar functions of position. This ansatz implies in particular that field
lines may be identified as the lines at which surfaces of constant $\alpha$
intersect those of constant $\beta$. Furthermore, it guarantees that the
divergence constraint $\nabla\cdot\vec{B}_{\mathrm{tr}}=0$ is always
satisfied. Here, we choose the first Euler potential $\alpha$ as

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle
B_{\mathrm{tr,0}}\,\sqrt{-2\,\Psi}\,\sin(\varphi-\varphi_{0})$ (26)
$\displaystyle=$ $\displaystyle
B_{\mathrm{tr,0}}\,\sqrt{(r-z)\,(r+z-1)}\,\sin(\varphi-\varphi_{0})$

in complete analogy to Equation (7) in Isenberg et al. (2015). This ensures
that (i) a streamline cannot leave its surface of constant $\alpha$ and that
(ii) said surfaces thus form an infinite bundle of nonintersecting sheets that
are parallel to the boundary field and to each other at infinity, but then
drape around the parabolic obstacle like a vertical curtain.

###  3.2 Isochrones

The second Euler potential $\beta$ is chosen as the negative of the travel
time $t$ that a fluid particle spent moving from upstream infinity along a
trajectory $t\mapsto(\rho(t),z(t))$ to its current position, and lines (or, in
three dimensions, surfaces) of constant $\beta$ will therefore be referred to
as isochrones. The minus sign, which is absent from the formulas used by
Isenberg et al. (2015), has to be added to either $\alpha$ or $\beta$ to
ensure consistency with the desired boundary conditions (17)–(19), and we
chose $\beta$ for reasons that will become clear in Section 4.2.

From the definitions of the velocity components $u_{\rho}$ and $u_{z}$, which
then read

$\left(u_{\rho},u_{z}\right)=\left(-\frac{{\rm d}\rho(\beta)}{{\rm
d}\beta},-\frac{{\rm d}z(\beta)}{{\rm d}\beta}\right),$ (27)

we obtain the respective formulas for $\beta$ as

$-\int_{a}^{\rho}\frac{{\rm
d}\rho^{\prime}}{u_{\rho}(\rho^{\prime})|_{\Psi}}=\beta=-\int_{\infty}^{z}\frac{{\rm
d}z^{\prime}}{u_{z}(z^{\prime})|_{\Psi}}.$ (28)

The symbol $\Psi$ in the integrand indicates that integration has to be
performed along a streamline, i.e. at constant $\Psi$, and $a$ is the value
attained by the $\rho$ coordinate of the streamline in question at
$z\rightarrow\infty$, i.e. the value satisfying

$\displaystyle\Psi(\rho,z)$ $\displaystyle=$
$\displaystyle\Psi(a,\infty)=-a^{2}/2$ (29) $\displaystyle\Rightarrow\
a(\rho,z)$ $\displaystyle=$
$\displaystyle\sqrt{-2\,\Psi(\rho,z)}=\sqrt{(r-z)\,(r+z-1)}\,.$

This allows us to use $a$ rather than $\Psi$ to label streamlines.

A minor conceptual problem arises from the unavoidable fact that any travel
time computed via Equation (28) will necessarily be infinite. Practical
computations, such as that of the resulting magnetic field components via
Equation (25), therefore require the integral to start at a finite height
$z_{0}$ (and a corresponding axial distance $\rho_{0}$) to keep the travel
time finite, and only afterwards perform the limit
$(\rho_{0},z_{0})\rightarrow(a,\infty)$ along the streamline in question (but
see Section 3.3 for a remedy).

While both expressions in Equation (28) are fully equivalent from a
mathematical point of view, tentatively solving Equation (10) for $\rho$ and
substituting into $1/u_{z}$ in order to eliminate $\rho$ in favor of $a$
indicates that the first expression in Equation (28) should be preferred. We
may thus solve Equation (10) for $z$ to obtain

$\displaystyle 2\,z$ $\displaystyle=$
$\displaystyle\frac{\rho^{2}}{\rho^{2}+2\,\Psi}-\rho^{2}-2\,\Psi=\frac{\rho^{2}}{\rho^{2}-a^{2}}-\rho^{2}+a^{2}$
(30)

and, after a series of elementary simplifications, are lead to the indefinite
integral

$\displaystyle T(\rho,a)$ $\displaystyle:=$ $\displaystyle\int\frac{{\rm
d}\rho}{u_{\rho}(\rho,a)}=\int\left(1+\frac{\rho^{2}}{(\rho^{2}-a^{2})^{2}}\right)\rho\,{\rm
d}\rho$ (31) $\displaystyle=$
$\displaystyle\int\left(\frac{w+1}{2\,w}+\frac{a^{2}}{2\,w^{2}}\right){\rm
d}w=\frac{w+\ln(w)}{2}-\frac{a^{2}}{2\,w}$ $\displaystyle=$
$\displaystyle\frac{\rho^{2}-a^{2}}{2}+\frac{\ln\left(\rho^{2}-a^{2}\right)}{2}-\frac{a^{2}}{2\,(\rho^{2}-a^{2})}\,,$

where $w=\rho^{2}-a^{2}\geq 0$ has been substituted temporarily.

The travel time of a fluid element moving from $(\rho_{0},z_{0})$ to
$(\rho,z)$ along a streamline $a=a(\rho_{0},z_{0})=a(\rho,z)$ is therefore
given by the difference $T(\rho,a)-T(\rho_{0},a)$. In particular, we have

$\beta=T(a,a)-T(\rho,a)$ (32)

as the second Euler potential, wherein the $z$-dependence comes through
$a=a(\rho,z)$. It is both noteworthy and fortunate that the exact integral
(31) can be found easily and without recourse to special functions such as
incomplete elliptic integrals, as is the case for the Rankine heliosphere
(Röken et al., 2015).

Since the first term in Equation (32) tends to $-\infty$, we temporarily
evaluate it at a finite but undetermined position $(\rho_{0},z_{0})$ on the
same streamline, which results in

$\displaystyle\hat{\beta}(\rho,z;z_{0}):=T\big{(}\rho_{0},a(\rho,z)\big{)}-T\big{(}\rho,a(\rho,z)\big{)}$
$\displaystyle=\frac{1}{2}\ln\bigg{(}\frac{\sqrt{(r-z)\,(r+z-1)+(z_{0}-1/2)^{2}}}{r-z}-\bigg{.}$
$\displaystyle\quad\bigg{.}-\frac{z_{0}-1/2}{r-z}\bigg{)}+z-z_{0}\,,$ (33)

wherein the condition $a(\rho_{0},z_{0})=a(\rho,z)$ has been used to express
$\rho_{0}$ as a function of $\rho$, $z$, and $z_{0}$.

Next, we compute $\nabla\hat{\beta}$, and only then let
$z_{0}\rightarrow\infty$, obtaining

$\displaystyle\nabla\beta$ $\displaystyle=$
$\displaystyle\lim_{z_{0}\rightarrow\infty}\nabla\hat{\beta}=\frac{\rho}{2\,(r+z-1)\,r}\,\vec{e}_{\rho}+$
(34) $\displaystyle+\left(1+\frac{r+z}{2\,(r+z-1)\,r}\right)\vec{e}_{z}.$

With $\nabla\alpha$ known as well, $\vec{B}_{\mathrm{tr}}$ itself is finally
found via Equation (25). As a result of this procedure, we obtain

$\displaystyle\frac{(B_{\mathrm{tr}})_{\rho}}{B_{\mathrm{tr,0}}}$
$\displaystyle=$
$\displaystyle\cos(\varphi-\varphi_{0})\,\sqrt{\frac{r+z}{r+z-1}}\times$ (35)
$\displaystyle\times\left(1-\frac{r-z}{2\,(r+z)\,r}\right)$
$\displaystyle\frac{(B_{\mathrm{tr}})_{\varphi}}{B_{\mathrm{tr,0}}}$
$\displaystyle=$
$\displaystyle-\sin(\varphi-\varphi_{0})\,\sqrt{\frac{r+z}{r+z-1}}$ (36)
$\displaystyle\frac{(B_{\mathrm{tr}})_{z}}{B_{\mathrm{tr,0}}}$
$\displaystyle=$
$\displaystyle-\frac{\cos(\varphi-\varphi_{0})}{2\,r}\,\sqrt{\frac{r-z}{r+z-1}}\,.$
(37)

The total, solenoidal, most general magnetic field solution to Equation (4),
subject to constant boundary conditions, is then simply the sum of the
respective longitudinal and transversal components. Appendix B summarizes
explicit expression for this total field in several coordinate systems.

###  3.3 Isochrones revisited

For the applications to be discussed in Section 4, it would be convenient to
have at our disposal a set of isochrones which, while leading to the same
magnetic field, are finite-valued in the region of interest. We will now
introduce a procedure to arrive at such an isochrone field.

Our starting point is the “quasi-isochrone” (3.2), which only becomes a real
isochrone in the limit $z_{0}\rightarrow\infty$, passing through the point
$(\rho,z)$. The idea is to exploit the gauge freedom of the Euler potential by
adding to $\hat{\beta}$ the fixed travel time of a fluid element traveling
along the $z$-axis from $z=z_{0}$ to some other reference height
$z={z_{\mathrm{a}}}$. This reference travel time is conveniently calculated to
be

$\displaystyle\Delta\beta(z_{0},{z_{\mathrm{a}}})$ $\displaystyle:=$
$\displaystyle\int^{z_{\mathrm{a}}}_{z_{0}}\frac{{\rm
d}z}{u_{z}|_{\rho=0}}=\int^{z_{\mathrm{a}}}_{z_{0}}\frac{{\rm
d}z}{(2\,z)^{-1}-1}$ (38) $\displaystyle=$
$\displaystyle-\left[z+\frac{\ln\bigl{(}2\,z-1\bigr{)}}{2}\right]^{z_{\mathrm{a}}}_{z_{0}},$

wherein we have used the fact that ${z_{\mathrm{a}}},z_{0}>1/2$. The choice of
${z_{\mathrm{a}}}$ is arbitrary, and we choose it such that the square bracket
in Equation (38) vanishes at $z={z_{\mathrm{a}}}$, leading to the implicit
equation

$2\,{z_{\mathrm{a}}}+\ln(2\,{z_{\mathrm{a}}}-1)=0\,.$ (39)

The actual single-valued solution to this equation, which could be written as

${z_{\mathrm{a}}}=\frac{W\left(e^{-1}\right)+1}{2}\approx 0.639$ (40)

using Lambert’s $W$ function, is immaterial for our purpose. What does matter
is that it causes our renormalized Euler potential to assume the particularly
simple form

$\displaystyle{\beta_{\mathrm{ren}}}(\rho,z)$ $\displaystyle:=$
$\displaystyle\lim_{z_{0}\rightarrow\infty}\big{[}\hat{\beta}(\rho,z;z_{0})+\Delta\beta({z_{\mathrm{a}}},z_{0})\big{]}$
(41) $\displaystyle=$ $\displaystyle z+\frac{\ln(r+z-1)}{2},$

a delightfully compact expression satisfying ${\beta_{\mathrm{ren}}}(0,1)=1$
as a result of our gauge (39). Using Equation (34), it is straightforward to
verify that $\nabla\beta=\nabla{\beta_{\mathrm{ren}}}$, which ensures that
replacing $\beta$ by ${\beta_{\mathrm{ren}}}$ in Equation (25) does not change
the resulting magnetic field.

###  3.4 Different aspect ratios

The fact that our parabolic surface has a fixed aspect ratio may pose a
problem for applications which require the obstacle’s shape to be more
elongated or compressed in the direction of its axis of symmetry. Suppose that
one wishes to apply a linear scaling transformation

$(\rho,z)\rightarrow(\ell_{\rho}\,\rho,\ell_{z}\,z)$ (42)

with $\ell_{\rho},\ell_{z}\in\mathbb{R}_{>0}$, thus keeping the focal point
fixed at the origin but moving the stagnation point from $z=1/2$ to
$z=\ell_{z}/2$ and the cylindrical radius at which the paraboloid intersects
the $(z=0)$ plane from $\rho=1$ to $\rho=\ell_{\rho}$. In an application, we
might for instance imagine $\ell_{\rho}$ and $\ell_{z}$ to represent the
object’s actual spatial scale, expressed in problem-specific physical units of
length. Then, we may construct from the original $\vec{u}$ and $\vec{B}$ the
scaled fields

$\displaystyle\vec{u}^{\,\mathrm{sc}}(\vec{r}\,)$ $\displaystyle:=$
$\displaystyle(\ell_{\rho}/\ell_{z})\,u_{\rho}(\vec{r}_{\ell})\,\vec{e}_{\rho}+u_{z}(\vec{r}_{\ell})\,\vec{e}_{z}$
(43) $\displaystyle\vec{B}_{\mathrm{lo}}^{\mathrm{sc}}(\vec{r}\,)$
$\displaystyle:=$ $\displaystyle-B_{z0}\,\vec{u}^{\,\mathrm{sc}}(\vec{r})$
(44) $\displaystyle\vec{B}_{\mathrm{tr}}^{\mathrm{sc}}(\vec{r}\,)$
$\displaystyle:=$ $\displaystyle
B_{\mathrm{tr},\rho}(\vec{r}_{\ell})\,\vec{e}_{\rho}+B_{\mathrm{tr},\varphi}(\vec{r}_{\ell})\,\vec{e}_{\varphi}+$
(45)
$\displaystyle+(\ell_{z}/\ell_{\rho})\,B_{\mathrm{tr},z}(\vec{r}_{\ell})\,\vec{e}_{z},$

in which the shorthands

$\displaystyle\vec{r}$ $\displaystyle:=$
$\displaystyle\rho\,\vec{e}_{\rho}+z\,\vec{e}_{z}$ (46)
$\displaystyle\vec{r}_{\ell}$ $\displaystyle:=$
$\displaystyle(\rho/\ell_{\rho})\,\vec{e}_{\rho}+(z/\ell_{z})\,\vec{e}_{z}$
(47)

have been used as arguments. These fields $\vec{u}^{\,\mathrm{sc}}$ and
$\vec{B}^{\mathrm{sc}}=\vec{B}_{\mathrm{lo}}^{\mathrm{sc}}+\vec{B}_{\mathrm{tr}}^{\mathrm{sc}}$
can be shown to be (i) divergence-free, (ii) tangential to the scaled
paraboloid, and (iii) an exact solution to the induction equation (4), subject
to the same upstream boundary conditions as the original fields $\vec{u}$ and
$\vec{B}$ prior to scaling. The last statement also follows from the more
general proof given by Kleimann et al. (2016) in the context of so-called
distortion flows, for which Equation (42) is a simple example.

Furthermore, given that said proof does not require any spatial symmetries, we
may even take this idea one step further and replace the cylindrically
symmetric scaling of Equation (42) with a still more general

$(x,y,z)\rightarrow(\ell_{x}\,x,\ell_{y}\,y,\ell_{z}\,z)\,.$ (48)

This leads to what could be called a “triaxial paraboloid,” which intersects
planes of constant $z$ in ellipses rather than circles, with the ellipse’s
semi-major axes at $z=0$ given by $\ell_{x}$ and $\ell_{y}$. The corresponding
(nonpotential) Cartesian flow field is then

$\vec{u}^{\,\mathrm{sc}}=\frac{1}{2\,r_{\ell}}\left[\frac{x\,\vec{e}_{x}+y\,\vec{e}_{y}}{r_{\ell}\,\ell_{z}+z}+\vec{e}_{z}\right]-\vec{e}_{z}\,,$
(49)

wherein

$r_{\ell}=\sqrt{(x/\ell_{x})^{2}+(y/\ell_{y})^{2}+(z/\ell_{z})^{2}}$ (50)

is the triaxial analog to the absolute value of Equation (47). Explicit
cylindrical and Cartesian components of $\vec{B}^{\mathrm{sc}}$ for the
respective scaling laws (42) and (48) are provided in Appendix B alongside the
unscaled ones.

It seems worth noting that, while Equations (43)–(45), as well as Equation
(49) together with its associated magnetic field (see Equations (B11)–(B13) in
Appendix B) do constitute valid MHD solutions, the scaled velocity fields (43)
and (49) are irrotational only for $\ell_{\rho}=\ell_{z}$ and
$\ell_{x}=\ell_{y}=\ell_{z}$, respectively, and therefore in general do not
have scalar potentials. The exact potential flows around scaled paraboloids
could probably only be found be solving the Laplace equation (3) for these
geometries, which is beyond the scope of this paper.

##  4 Discussion and application

###  4.1 General field structure

Figure 6 shows renderings of selected field lines from several different
vantage points. As expected, the field is almost homogeneous at larger
distances from the obstacle, only to smoothly drape around it in the vicinity
of the latter. Along the $z$-axis, the respective magnitudes of the
longitudinal and transversal components are

$\displaystyle\|\vec{B}_{\mathrm{lo}}|_{\rho=0}\|$ $\displaystyle=$
$\displaystyle\left(1-\frac{1}{2\,z}\right)|B_{z0}|$ (51)
$\displaystyle\|\vec{B}_{\mathrm{tr}}|_{\rho=0}\|$ $\displaystyle=$
$\displaystyle\left(1-\frac{1}{2\,z}\right)^{-1/2}\sqrt{(B_{x0})^{2}+(B_{y0})^{2}}$
(52)

implying that in the limit $z\rightarrow 1/2$ (i.e. at the stagnation point,
located where the $z$-axis intersects the obstacle’s surface), the
longitudinal component vanishes while the transversal component tends to
infinity. This behavior is analogous to that of the Rankine heliosphere, where
the corresponding term in round brackets reads $(1-1/z^{2})$ and likewise
diverges as $z\rightarrow 1$ (see Appendix C of Röken et al., 2015). Indeed,
it can be verified that this behavior is not restricted to the stagnation
point but extends to the entire surface of the paraboloid, also in exact
analogy to the heliosphere case. The reason for this divergent field strength
is found in a mismatch between the conflicting assumptions of stationarity and
idealness, which enforces an infinite pile-up of magnetic flux than cannot be
dissipated away. In a more realistic nonideal setting, the rate of incoming
advected flux would be balanced by resistive dissipation, allowing for an
equilibrium situation with finite field strength. This, however, is beyond the
scope of the present approach, which crucially relies on the frozen-in
condition that is ensured by the ideal induction equation (4). A viable option
to mitigate this problem for the heliopause, which would likely also work in
the present case, has been presented by Florinski et al. (2019).

Figure 2:  Left: four perspective renderings of field lines for boundary
conditions $(B_{x0},B_{y0},B_{z0})=(1,0,1)$ draping around the obstacle
(yellow) and colored according to $\|\vec{B}\|$. Since the latter tends to
infinity near the obstacle’s surface, contour values have been clamped within
the interval $[0.5,4.0]$. Right: slices along the $y=0$ plane showing field
lines of the longitudinal (top) and transversal (bottom) parts of the magnetic
field, with the background again color-coded to represent field strength in
this plane. The longitudinal field (whose field lines coincide with
streamlines) vanishes at the stagnation point, while the transversal field
(whose field lines coincide with isochrones) diverges on the paraboloid’s
surface, but otherwise remains close to its upstream boundary value almost
everywhere else.

###  4.2 Arbitrary base vector fields

The choice of homogeneous boundary conditions (17)–(19) was chiefly motivated
both as the simplest nontrivial choice, and also because typical astrophysical
applications usually consider the incoming medium indeed to be largely
homogeneous. We may, however, generalize the procedure employed here to almost
arbitrary fields and study how they are deformed and modified by the presence
of the solid obstacle which they have to smoothly avoid.

The underlying idea, which was already employed by Röken et al. (2015), is to
note that the streamlines and (renormalized) isochrones as given by the
respective Equations (3.2) and (41) can be interpreted as establishing a non-
orthogonal coordinate system filling the entire volume exterior to the
obstacle, and that, according to Cauchy’s integral (Cauchy, 1816), the
components of a field advected in the flow, when expressed in these
coordinates, will remain constant throughout the entire transport. The
transformation expressed in Equations (3.2) and (41) then serves to map
coordinates $(\rho,z)$ of the distorted space, at which the solution is well-
defined only exterior to the obstacle, to
$(a,{\beta_{\mathrm{ren}}})\in\mathbb{R}_{\geq 0}\times\mathbb{R}$, the
coordinates of undistorted space.

A rather illustrative way of applying this mapping to vector-valued quantities
is to think of a vector $\vec{B}(\rho,\varphi,z)$ as an actual “arrow” of
infinitesimal length pointing from coordinates $(\rho,\varphi,z)$ to
$(\rho+\delta\rho,\varphi+\delta\varphi,z+\delta z)$, and to apply the above
transformation from undistorted $(a,\varphi,{\beta_{\mathrm{ren}}})$ to
distorted $(\rho,\varphi,z)$ coordinates separately to the “base” and the
“tip” of the arrow. The respective vector components are then

$\displaystyle\left(B_{\rho},B_{\varphi},B_{z}\right)$ $\displaystyle=$
$\displaystyle\left(\delta\rho,\rho\,\delta\varphi,\delta z\right)$ (53)
$\displaystyle\left(B_{0,a},B_{0,\varphi},B_{0,{{\beta_{\mathrm{ren}}}}}\right)$
$\displaystyle=$ $\displaystyle\left(\delta
a,a\,\delta\varphi,\delta{\beta_{\mathrm{ren}}}\right)$ (54)

wherein $\vec{B}_{0}(\rho,\varphi,z)$ is an “undistorted” magnetic field,
which no longer needs to be restricted to constant components from the
upstream boundary, but may vary across and fill the entire space (which we may
think of as being obstacle-free). Without the obstacle, both coordinate
systems coincide exactly. To achieve this very desirable property was the
motivation behind amending the second, rather than the first, Euler potential
with a minus sign in Equation (28), since this ensures that both
$\nabla{\beta_{\mathrm{ren}}}|_{\rho=0}$ and $\nabla z=\vec{e}_{z}$ are
pointing into the same direction.

When the obstacle is present, however, arguments $(\rho,\varphi,z)$ of
$\vec{B}_{0}$ are to be replaced by $(a,\varphi,{\beta_{\mathrm{ren}}})$,
which are themselves functions of $\rho$ and $z$. Comparing the second
($\varphi$-)component of Equations (53) and (54), we see that the toroidal
component is therefore simply

$B_{\varphi}(\rho,\varphi,z)=\frac{\rho}{a(\rho,z)}B_{0,\varphi}\big{(}a(\rho,z),\varphi,{\beta_{\mathrm{ren}}}(\rho,z)\big{)}\,.$
(55)

For the poloidal components, however, we have

$\displaystyle\left(B_{a},B_{{\beta_{\mathrm{ren}}}}\right)$ $\displaystyle=$
$\displaystyle\left(\delta a,\delta{\beta_{\mathrm{ren}}}\right)$
$\displaystyle=$ $\displaystyle\big{(}a(\rho+\delta\rho,z+\delta
z)-a(\rho,z),$ $\displaystyle\hskip
7.11317pt{\beta_{\mathrm{ren}}}(\rho+\delta\rho,z+\delta
z)-{\beta_{\mathrm{ren}}}(\rho,z)\big{)}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial a}{\partial\rho}\delta\rho+\frac{\partial
a}{\partial z}\delta
z,\frac{\partial{\beta_{\mathrm{ren}}}}{\partial\rho}\delta\rho+\frac{\partial{\beta_{\mathrm{ren}}}}{\partial
z}\delta z\right).$

The derivatives of $a(\rho,z)$ can be easily evaluated using the general
relations (13) and (29) as

$\displaystyle\frac{\partial a}{\partial\rho}$ $\displaystyle=$
$\displaystyle\frac{\partial~{}}{\partial\rho}\sqrt{-2\,\Psi}=-\frac{1}{\sqrt{-2\,\Psi}}\frac{\partial\Psi}{\partial\rho}=-\frac{1}{a}\,(\rho\,u_{z})$
(57) $\displaystyle\frac{\partial a}{\partial z}$ $\displaystyle=$
$\displaystyle\frac{\partial~{}}{\partial
z}\sqrt{-2\,\Psi}=-\frac{1}{\sqrt{-2\,\Psi}}\frac{\partial\Psi}{\partial
z}=\frac{1}{a}\,(\rho\,u_{\rho})$ (58)

independently of the specific form of $a$ or $\Psi$ (and hence independently
of the shape of the obstacle). Inserting these into Equation (4.2) and solving
for the desired poloidal field components (53) provides the latter as

$\displaystyle\left(\begin{array}[]{c}B_{\rho}\\\
B_{z}\end{array}\right)=\frac{1}{D}\left(\begin{array}[]{rc}-(a/\rho)(\partial_{z}{\beta_{\mathrm{ren}}})&\
\,\,\,u_{\rho}\\\ (a/\rho)(\partial_{\rho}{\beta_{\mathrm{ren}}})&\
\,\,\,u_{z}\end{array}\right)\left(\begin{array}[]{c}B_{a}\\\
B_{\beta_{\mathrm{ren}}}\end{array}\right)$ (65)

with the determinant

$\displaystyle D$ $\displaystyle:=$
$\displaystyle\vec{u}\cdot\nabla{\beta_{\mathrm{ren}}}=\frac{\partial\rho}{\partial
t}\,\frac{\partial{\beta_{\mathrm{ren}}}}{\partial\rho}+\frac{\partial
z}{\partial t}\,\frac{\partial{\beta_{\mathrm{ren}}}}{\partial z}$ (66)
$\displaystyle=$ $\displaystyle\frac{{\rm d}{\beta_{\mathrm{ren}}}}{{\rm
d}t}=-1\,.$

The shape of the obstacle only enters through the partial derivatives of
${\beta_{\mathrm{ren}}}$, which are already known from Equation (34) because
$\beta$ and ${\beta_{\mathrm{ren}}}$ share the same gradient by construction.
Inserting these into Equation (65) together with the components (14) and (15)
of $\vec{u}$ and performing a series of algebraic simplifications, the
explicit form of the two poloidal components finally becomes

$\displaystyle B_{\rho}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r+z}{r+z-1}}\left(1-\frac{r-z}{2\,(r+z)\,r}\right)B_{a}$
$\displaystyle-\frac{1}{2\,r}\,\sqrt{\frac{r-z}{r+z}}\,B_{\beta_{\mathrm{ren}}}$
$\displaystyle B_{z}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\,r}\,\sqrt{\frac{r-z}{r+z-1}}\,B_{a}+\left(1-\frac{1}{2\,r}\right)B_{\beta_{\mathrm{ren}}}\,.$
(68)

Note that this is consistent with the formulas (B1) and (B3) found earlier for
constant boundary conditions.

###  4.3 A simple example

To illustrate the procedure outlined above using an astrophysically relevant
example, let us consider the force-free magnetic field

$\vec{B}=\frac{\rho\,\vec{e}_{\varphi}+\vec{e}_{z}}{1+\rho^{2}}$ (69)

whose field lines form helical screw-like spirals of constant axial radius and
constant thread lead $2\pi$. To make this example more interesting, we reduce
its symmetry by rotating it by $90^{\circ}$ around the $y$-axis such that the
$x$-axis becomes the new axis of symmetry, displace this axis of symmetry from
$(y,z)=(0,0)$ to $(y_{\mathrm{c}},z_{\mathrm{c}})$, and then re-express it in
cylindrical coordinates. The result reads

$\displaystyle\vec{B}$ $\displaystyle=$
$\displaystyle\bigl{(}\left[\cos\varphi-(z-z_{\mathrm{c}})\sin\varphi\right]\vec{e}_{\rho}\bigr{.}+$
(70)
$\displaystyle\left[-\sin\varphi-(z-z_{\mathrm{c}})\cos\varphi\right]\vec{e}_{\varphi}+\big{.}\left[\rho\sin\varphi-
y_{\mathrm{c}}\right]\vec{e}_{z}\bigr{)}\times$
$\displaystyle\bigl{[}(\rho\sin\varphi-
y_{\mathrm{c}})^{2}+(z-z_{\mathrm{c}})^{2}+1\bigr{]}^{-1}.$

The application of Equations (55), (4.2), and (68) to this base field then
leads to the explicit form of the distorted field (not shown here for
brevity). Figure 3 presents a comparison of the original and the distorted
field. In an astrophysical application, we might imagine the paraboloid to
represent, for instance, a planetary magnetopause, while the sample field (70)
could model a solar/stellar magnetic flux tube that is incident upon the
magnetosphere.

In such a situation, the travel-time interpretation of
$-{\beta_{\mathrm{ren}}}$ allows the deformed magnetic field to be interpreted
in two distinct ways:

  1. 1.

According to the interpretation used in the previous derivation, one may
simply compare an undisturbed field to a version of itself disturbed by the
presence of the (in this case) parabolic obstacle.

  2. 2.

Alternatively, since these fields are frozen into a stationary flow field, the
situation can be viewed as a time-dependent advection in this flow. Therein,
the undistorted field, which is advected in the homogeneous flow
$\vec{u}_{\mathrm{bg}}=-\vec{e}_{z}$, satisfies

$\vec{B}(\rho,\varphi,z,t)=\vec{B}(\rho,\varphi,z+{\beta_{\mathrm{ren}}},t-{\beta_{\mathrm{ren}}})$
(71)

and can in principle be traced to a time-dependent boundary condition at
$z\rightarrow+\infty$.

This property carries over to the distorted version of said fields, an idea
which is also illustrated in the example of Figure 3 showing three instances
of the same flux rope at consecutive times (or, alternatively, a sequence of
distinct, vertically displaced flux ropes at a single moment of time).

Figure 3:  Left: perspective rendering of several instances of the sample flux
tube field (70) with its axis centered on $y_{\mathrm{c}}=\pm 2$ and
$z_{\mathrm{c}}\in\\{-3,0,3\\}$ in its original form ($y<0$) vs. its shape
when distorted by the presence of the paraboloid obstacle ($y>0$).
Additionally, a semi-transparent cylinder of radius $0.8$ marks the outer
layer of field lines in each case. Color of field lines encodes total field
strength, increasing from blue to red. These sets of field lines may be
interpreted as successive instances, viewed at uniformly spaced time
intervals, of the flux tube being advected in the flow. Middle: parallel
projection of the same situation, viewed from the $+x$ direction. Right: again
the same situation viewed from $+y$. Note that, since this is also a parallel
projection, the deformed flux tube segments in the front are indeed physically
longer than those on the far side, rather than just appearing so because of
their closer proximity to the observer.

###  4.4 Advecting arbitrary scalar fields

The method outlined above may also be applied to the simpler case of scalar,
rather than vector-valued quantities. Given such a scalar quantity
$Q(\rho,\varphi,z)$, its distorted counterpart is just
$Q\big{(}a(\rho,z),\varphi,{\beta_{\mathrm{ren}}}(\rho,z)\big{)}$, with the
dependence of $(a,{\beta_{\mathrm{ren}}})$ on $(\rho,z)$ again given by
Equations (29) and (41). A sample visualization using spheres, cylinders
aligned with either the $x$\- and $y$-direction, and cubes centered on
Cartesian coordinates $(x_{\mathrm{c}},y_{\mathrm{c}},z_{\mathrm{c}})$ is
shown in Figure 4. Note that, due to the incompressibility of the underlying
flow field, the total volume of all shapes depicted here is conserved at any
time. (This would no longer be the case for the compressible flow discussed in
Section 5, which has $\nabla\cdot\vec{u}^{\mathrm{C}}\neq 0$.) In the
aforementioned application, the initially spherical blobs might, for instance,
stand in for dust clouds or regions of enhanced or decreased temperature in
the stellar wind.

Figure 4:  The deformation of basic geometric shapes induced by the presence
of the parabolic obstacle. (A): sequences of red spheres of radius $0.8$ at
$(x_{\mathrm{c}},y_{\mathrm{c}},z_{\mathrm{c}})=(-1,-1,t_{z})$ for
$t_{z}\in\\{-4,-2,0,2\\}$, together with a second sequence of blue cubes of
side length $1.2$ at
$(x_{\mathrm{c}},y_{\mathrm{c}},z_{\mathrm{c}})=(1,1,t_{z})$ arranged parallel
to, but not intersecting, the inflow ($z$-)axis. (B): The same shapes are
severely distorted when passively advected in the flow, visualizing the
pancake-like compression and flattening of these shapes as they are advected
towards the obstacle, eventually sliding off it. (C): The same situation as
(A) but for a sequence of cylinders of height $1.6$ and radius $0.8$
alternating in color and orientation and placed right on the $z$-axis. (D):
The cylinders are shown being compressed into (eventually) infinitesimally
thin shells coating the paraboloid’s surface.

We call attention to the fact that the extension to arbitrary base fields only
became possible through the renormalization of the second Euler potential,
which is now directly associated with finite travel times.

###  4.5 The inverse transform

To see more clearly how the presence of the parabolic obstacle displaces a
fluid element, it would be very convenient to know the direct inverse
$(a,{\beta_{\mathrm{ren}}})\mapsto(\rho,z)$ of the transformation given by
Equations (29) and (41). Unfortunately, the implicit nature of the second
equation forestalls the determination of an analytic formula to this effect,
so here at last, a numerical procedure becomes mandatory. We will now present
and discuss such a procedure.

The first transformation formula (29) may be easily solved for $\rho$ to yield

$\rho=\sqrt{a^{2}-(z-1/2)+w_{a}(z)}$ (72)

with the shorthand

$w_{a}(z):=\sqrt{a^{2}+(z-1/2)^{2}}$ (73)

that will be of further use shortly. We then insert Equation (72) into the
second transformation (41) to obtain

${\beta_{\mathrm{ren}}}=z+\frac{\ln\big{(}w_{a}(z)+z-1/2\big{)}}{2}=:\beta_{a}(z)$
(74)

as an equation for the desired second coordinate $z$, which, while still being
implicit in nature, for given $(a,{\beta_{\mathrm{ren}}})$ only depends on a
single argument. Equation (74) may be interpreted as giving the isochrone
label of a fluid element located at height $z$ while sliding along streamline
$a$. It is clear from this interpretation, and also easy to verify via

$\beta_{a}^{\prime}(z)\equiv\frac{{\rm d}\beta_{a}(z)}{{\rm
d}z}=1+\frac{1}{2\,w_{a}(z)}>0,$ (75)

that $z\mapsto\beta_{a}(z)$ is strictly monotonous. Being additionally also
strictly concave ($\beta_{a}^{\prime\prime}(z)<0$ for $z>1/2$), our implicit
Equation (74) is well suited for an iterative solver based on Newton’s method,
which in this case reads

$\displaystyle z_{k+1}$ $\displaystyle=$ $\displaystyle
z_{k}-\frac{\beta_{a}(z_{k})-{\beta_{\mathrm{ren}}}}{\beta_{a}^{\prime}(z_{k})}$
$\displaystyle=$
$\displaystyle\frac{z_{k}+w_{a}(z_{k})\left[2\,{\beta_{\mathrm{ren}}}-\ln\big{(}w_{a}(z_{k})+z_{k}-1/2\big{)}\right]}{1+2\,w_{a}(z_{k})}$

and for a given starting value, which we place at $z_{1}=1$, quickly converges
to the desired solution $z=\lim_{k\rightarrow\infty}z_{k}$ of Equation (74).

For the special case $a=0$, we need to observe that, since this is the only
streamline that does not extend infinitely into the $-z$ direction but rather
ends at the stagnation point, it may happen that $z_{k+1}\leq 1/2$ for some
$k$, indicating a position on or inside the paraboloid. A possible strategy,
which we shall adopt here without evaluating further options for its
optimization, is then to replace $z_{k+1}$ by $(1/2+z_{k})/2$, i.e. the
midpoint between the stagnation point and the previous estimate. This ensures
that $z_{k+1}$ is both a valid external position and closer to the desired
result than $z_{k}$ had been.

Since Newton’s method converges quadratically, satisfactory results are
typically obtained within just a few iterations. Tentative computations for
target values on a grid of $(a,{\beta_{\mathrm{ren}}})\in[0,4]\times[-8,4]$
with a spacing of $0.1$ in both directions show that it takes between $N=2$
and $7$ iterations ($3.34$ on average) to ensure
$|\beta_{a}(z_{1+N})-{\beta_{\mathrm{ren}}}|<10^{-9}$ if $a>0$. On the $a=0$
axis, $\sim$13 iterations are typically needed (except for
${\beta_{\mathrm{ren}}}=1$, where $z=1$ by construction).

##  5 Extension to Compressible flow

Kleimann et al. (2017) presented a generalization of the Röken et al. (2015)
MHD model of the outer heliosheath to compressible flow, thus allowing number
density $n$ to vary in a physically more realistic way. We will now derive the
corresponding case for the flow around the paraboloid.

###  5.1 Flow and density field

The key idea behind said generalization is to reinterpret the potential $\Phi$
as one for momentum density $n\,\vec{u}$ rather than $\vec{u}$, thereby
ensuring that the flow’s streamline structure, and in particular the shape of
the obstacle, remains unchanged while permitting $n$ to vary along
streamlines. This approach has the appealing feature that Equation (3)
automatically ensures $\nabla\cdot\left(n\,\vec{u}\right)=0$, implying that
the compressible extension of our model continues to conserve mass by
construction. The newly introduced degree of freedom is handled by
additionally requiring the momentum balance equation

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

to be satisfied when projected upon streamlines, with $m_{\mathrm{p}}$
denoting the particle (proton) mass and $P$ the gas pressure. The latter is in
turn given by the polytropic equation of state

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

where $\gamma$ is the adiabatic index and the constant $C$ is to be determined
from a boundary condition at infinity.

Fortunately, most of the derivation presented in Kleimann et al. (2017) does
not make reference to a particular flow field, and is thus also applicable to
the present situation without any modifications. In particular, the desired
density field $n$ may be obtained from the same respective implicit equations,
namely

$m^{2}{\cal A}=n^{2}\left[m^{2}-\frac{2\,(n^{\gamma-1}-1)}{\gamma-1}\right]$
(79)

if $\gamma\neq 1$, or

$m^{2}{\cal A}=n^{2}\left[m^{2}-2\,\ln(n)\right]$ (80)

for the isothermal case $\gamma=1$. Here, the parameter $m\in[0,1]$ is the
sonic Mach number at upstream infinity and ${\cal A}$ is a scalar function of
meridional position $(\rho,z)$, which, for the present paraboloid case,
evaluates to

$\displaystyle{\cal A}:=(\nabla\Phi)^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{2\,r}\,\sqrt{\frac{r-z}{r+z}}\,\right)^{2}+\left(\frac{1}{2\,r}-1\right)^{2}$
(81) $\displaystyle=$ $\displaystyle 1-\frac{1}{r}+\frac{1}{2\,(r+z)\,r}\,.$

This function maps the entire region exterior to the obstacle into the
interval $[0,1]$, with the limiting cases of $0$ and $1$ being exclusively
obtained at the stagnation point and at upstream infinity, respectively. Note
that both Equation (79) and (80) are consistent with the uniform density case
$n=1$ being recovered in the limit $m\rightarrow 0$.

Continuing the analogy to Kleimann et al. (2017), the isothermal case (80) can
be “solved” using Lambert’s W function as

$n({\cal
A})=\exp\left(\frac{m^{2}}{2}+\frac{1}{2}\,W\left(-m^{2}\exp(-m^{2})\,{\cal
A}\right)\right),$ (82)

indicating that the peak density at the stagnation point amounts to
$n(0)=\exp(m^{2}/2)\leq\exp(1/2)\approx 1.649$.

For the mono-atomic ideal gas case $\gamma=5/3$, upon which we shall
specialize in the following besides $\gamma=1$, solving Equation (79) is
equivalent to finding the only real-valued root of

$\xi^{4}-(1+m^{2}/3)\,\xi^{3}+(m^{2}/3){\cal A}\,,$ (83)

a fourth-order polynomial in $\xi:=n^{2/3}$.

Although the quartic (83) can in principle be solved using Ferrari’s formula,
a more convenient expression would be desirable, and the same is certainly
true for Equation (82) involving Lambert’s $W$ function. For this reason, we
repeat the same steps already applied by Kleimann et al. (2017), namely to
replace the both Equation (82) and the exact (but implicit) expression (83)
for $n({\cal A})$ by a second-order polynomial $n_{2}({\cal A})$ whose
coefficients are adjusted such that it agrees with $n({\cal A})$ at ${\cal
A}=0$ (the stagnation point), at ${\cal A}=1$ (infinity) and also has the same
derivative at ${\cal A}=0$. A simpler, less accurate alternative is the first-
order polynomial $n_{1}({\cal A})$ which only fulfills the first two
conditions. Reproducing Equations (44)–(47) from Kleimann et al. (2017), we
get

$\displaystyle n_{1}({\cal A})$ $\displaystyle=$ $\displaystyle
n_{\mathrm{sp}}+(1-n_{\mathrm{sp}})\,{\cal A}$ (84) $\displaystyle n_{2}({\cal
A})$ $\displaystyle=$ $\displaystyle n_{\mathrm{sp}}+D_{0}\,{\cal
A}+(1-n_{\mathrm{sp}}-D_{0})\,{\cal A}^{2}$ (85)

with

$n_{\rm sp}:=n(0)=\left\\{\begin{array}[]{lcl}\exp\left(m^{2}/2\right)&\
:&\gamma=1\\\\[5.69046pt] \left(1+m^{2}/3\right)^{3/2}&\
:&\gamma=5/3\end{array}\right.$ (86)

the stagnation point density and

$D_{0}:=-\frac{m^{2}}{2}\times\left\\{\begin{array}[]{lcl}\exp\left(-m^{2}/2\right)&\
:&\gamma=1\\\\[5.69046pt] \displaystyle\left(1+m^{2}/3\right)^{-5/2}&\
:&\gamma=5/3\end{array}\right.$ (87)

the derivative $n^{\prime}(0)$ of $n({\cal A})$ at that point.

Figure 5 shows contours of density for selected values of $m$ relevant for the
isothermal case, in each case comparing the exact formula (82) with its
approximation (85). As can be seen, both sides show a very similar behavior,
and, by construction, agree exactly at both the stagnation point and towards
infinity. As can also be seen from the figure, the density at any exterior
point is never below the density of incompressible flow. This constitutes a
noteworthy departure from the heliosphere case (Kleimann et al., 2017)
featuring an extended region in the flanks in which the flow speed is higher,
and the density correspondingly lower, than in the incompressible situation.
In general, the density of the isothermal case is very similar to that of the
adiabatic case, which is not shown for exactly this reason.

It is worth noting that, while the polynomial approximation can be evaluated
at any position $(\rho,z)$, the “exact” solution exhibits a region in which it
is undefined (indicated in white in Figure 5). In the heliosphere case, this
region grows as $m$ increases, eventually engulfing the entire downwind half-
space as $m\rightarrow 1$. This, however, is apparently different for the
paraboloid, where it remains confined to the interior of the obstacle, thereby
allowing the exact solution (82) to be used for any $m\leq 1$ (and technically
even beyond that, though of course no bow shock can be expected) if so
desired.

Figure 5:  Density contours for three different values of upstream Mach number
$m$, in each case plotted using both the exact profile (Equation (82), left
side) vs. the polynomial approximation (Equation (85), right side). The red
curve marks the parabolic obstacle. Values below 0.7 are not shown, and no
exact solution exists in the white region.

###  5.2 The modified magnetic field

Because the magnetic field previously derived from Equation (4) is frozen into
$\vec{u}$ rather than $n\,\vec{u}$, the question arises as to how the
transition from incompressible to compressible flow changes the associated
magnetic field. This question was already investigated by Kleimann et al.
(2017) for the Rankine-type heliosphere, resulting again in fully analytic,
albeit more involved, formulas for the modified components of $\vec{B}$.
Rather than applying the same procedure to the present parabolic case, we
approach the problem using the framework of Euler potentials, as was done for
the incompressible case in Section 3.

We begin by deriving the longitudinal component, which is parallel to
$\vec{u}^{\mathrm{C}}$, and may therefore be written as
$\vec{B}_{\mathrm{lo}}^{\mathrm{C}}=\omega\,\vec{u}^{\mathrm{C}}$ with
$\omega$ being a scalar function of position. From

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\nabla\cdot\vec{B}_{\mathrm{lo}}^{\mathrm{C}}=\nabla\cdot[\omega\,\vec{u}^{\mathrm{I}}/n_{1,2}]$
(88) $\displaystyle=$
$\displaystyle(\omega/n_{1,2})\,\underbrace{\nabla\cdot\vec{u}^{\mathrm{I}}}_{=0}+\left(\vec{u}^{\mathrm{I}}\cdot\nabla\right)(\omega/n_{1,2})$

we can conclude that $\omega/n_{1,2}$ is constant along any streamline, though
it may still vary across streamlines. However, consistency with the upstream
boundary condition can only be satisfied if $\omega/n_{1,2}$ equals the global
constant $-B_{z0}$. Therefore,
$\vec{B}_{\mathrm{lo}}^{\mathrm{C}}=-B_{z0}\,\vec{u}^{\mathrm{I}}=\vec{B}_{\mathrm{lo}}^{\mathrm{I}}$,
implying that the longitudinal component remains completely unaffected by the
transition to compressibility.

For the transversal component, we use the same first Euler potential
$\alpha=B_{\mathrm{tr,0}}\,a\,\sin(\varphi-\varphi_{0})=-a\,B_{\varphi 0}$ as
in the incompressible case, whose gradient evaluates to

$\displaystyle\nabla\alpha=\frac{\rho}{a}\,B_{\varphi
0}\,\bigl{(}u_{z}^{\mathrm{I}}\,\vec{e}_{\rho}-u_{\rho}^{\mathrm{I}}\,\vec{e}_{z}\bigr{)}+\frac{a}{\rho}\,B_{\rho
0}\,\vec{e}_{\varphi}$ (89)

according to Equations (57) and (58). Similarly to the associated
incompressible expression in Section 3.2, the gradient of the second Euler
potential then becomes

$\displaystyle\nabla\beta^{\mathrm{C}}$ $\displaystyle=$
$\displaystyle\lim_{z_{0}\rightarrow\infty}\nabla\hat{\beta}^{\mathrm{C}}=-\lim_{z_{0}\rightarrow\infty}\nabla\int_{\rho_{0}(\rho,z;z_{0})}^{\rho}\frac{{\rm
d}\rho^{\prime}}{u_{\rho}^{\mathrm{C}}}$ (90) $\displaystyle=$
$\displaystyle-\lim_{z_{0}\rightarrow\infty}\nabla\int_{\rho_{0}(\rho,z;z_{0})}^{\rho}\frac{n_{\mathrm{sp}}+\nu_{1}{\cal
A}+\nu_{2}{\cal A}^{2}}{u_{\rho}^{\mathrm{I}}}\,{\rm d}\rho^{\prime}$
$\displaystyle=$ $\displaystyle
n_{\mathrm{sp}}\,\nabla\beta^{\mathrm{I}}-\nu_{1}\,\vec{K}_{1}(\rho,z)-\nu_{2}\,\vec{K}_{2}(\rho,z)$

with the functions

$\vec{K}_{k}(\rho,z):=\lim_{z_{0}\rightarrow\infty}\nabla\int_{\rho_{0}(\rho,z;z_{0})}^{\rho}\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\,{\rm d}\rho^{\prime}\,,$ (91)

$k\in\\{1,2\\}$, explicitly specified in Appendix C and the constants
$\nu_{1,2}$ to be read off Equation (84) or (85).

With the above, we may express the magnetic field

$\vec{B}^{\mathrm{C}}=\vec{B}_{\mathrm{tr}}^{\mathrm{C}}+\vec{B}_{\mathrm{lo}}^{\mathrm{C}}=\nabla\alpha\times\nabla\beta^{\mathrm{C}}-B_{z0}\,\vec{u}^{\mathrm{I}}$
(92)

for compressible flow in terms of the magnetic field $\vec{B}^{\mathrm{I}}$
for incompressible flow given in Equations (4.2) and (68). Fortunately, the
toroidal part $B_{\varphi}^{\mathrm{C}}$ of the compressible magnetic field
can be related to the respective incompressible component
$B_{\varphi}^{\mathrm{I}}$ in a way that is independent of the specific forms
of the poloidal flow field and the gradient of the compressible second Euler
potential. This is evident from the short computation

$\displaystyle B_{\varphi}^{\mathrm{C}}$ $\displaystyle=$ $\displaystyle
B_{\varphi,\mathrm{tr}}^{\mathrm{C}}+B_{\varphi,\mathrm{lo}}^{\mathrm{C}}=\left[\nabla\alpha\times\nabla\beta^{\mathrm{C}}\right]_{\varphi}+0$
(93) $\displaystyle=$ $\displaystyle-\underbrace{(\rho/a)B_{\varphi
0}}_{=B_{\varphi}^{\mathrm{I}}}\,n_{1,2}\,\underbrace{\vec{u}^{\mathrm{C}}\cdot\nabla\beta^{\mathrm{C}}}_{=-1}=n_{1,2}\,B_{\varphi}^{\mathrm{I}}\,,$

which uses the compressible version of the determinantial relation (66). For
the poloidal parts $B_{\rho}^{\mathrm{C}}$ and $B_{z}^{\mathrm{C}}$, however,
the explicit form of the flow field as well as of the gradient of the
compressible second Euler potential become relevant. We thus obtain

$\displaystyle B_{\rho}^{\mathrm{C}}$ $\displaystyle=$ $\displaystyle
n_{\mathrm{sp}}\,B_{\rho}^{\mathrm{I}}-\mathcal{Z}_{z}+(n_{\mathrm{sp}}-1)\,u_{\rho}\,B_{z0}$
(94) $\displaystyle B_{z}^{\mathrm{C}}$ $\displaystyle=$ $\displaystyle
n_{\mathrm{sp}}\,B_{z}^{\mathrm{I}}+\mathcal{Z}_{\rho}+(n_{\mathrm{sp}}-1)\,u_{z}\,B_{z0}\,,$
(95)

where

$\displaystyle\vec{\mathcal{Z}}:=\frac{a}{\rho}\,\bigl{(}\nu_{1}\,\vec{K}_{1}+\nu_{2}\,\vec{K}_{2}\bigr{)}\,B_{\rho
0}\,.$ (96)

Despite being well-defined at any position outside the paraboloid’s surface,
the vector-valued function $\vec{K}_{2}$ is somewhat cumbersome to numerically
evaluate in the vicinity of the $(a=0)$ streamline (see Equation (C15) in
Appendix C). Therefore, we now compute expressions for the compressible field
components on the $z$-axis that are both well-defined and easy to evaluate
directly via the induction equation, yielding

$\displaystyle\frac{B_{x}^{\mathrm{C}}|_{\rho=0}}{B_{x0}}$ $\displaystyle=$
$\displaystyle
n_{1,2}|_{\rho=0}\,\biggl{(}1-\frac{1}{2\,z}\biggr{)}^{-1/2}=\frac{B_{y}^{\mathrm{C}}|_{\rho=0}}{B_{y0}}$
(97) $\displaystyle\frac{B_{z}^{\mathrm{C}}|_{\rho=0}}{B_{z0}}$
$\displaystyle=$ $\displaystyle 1-\frac{1}{2\,z}$ (98)

in analogy to Kleimann et al. (2017). We point out that these expressions are
of course compatible with the magnetic field components (93)–(95) properly
evaluated at $\rho=0$. Also, taking the limit $z\rightarrow\infty$ of the
latter components gives rise to the boundary field $\vec{B}_{0}$ at upstream
infinity. Furthermore, it seems worth noting that these results are
“approximative” only in the sense that $n_{1,2}$ do not exactly fulfill the
momentum balance equation (77). The magnetic field thus derived, however, is
again an exact solution to the induction equation (4) for
$u^{\mathrm{C}}=n_{1,2}\,u^{\mathrm{I}}$.

Figure 6 compares the compressible field to its incompressible counterpart.
While both drape around the obstacle in a qualitatively similar way, it can be
seen that the compressible field tends to be less sharply curved, and also
starts to deviate from straight lines already at larger upwind distances,
which is consistent with the general behavior found in the compressible
heliospheric case. Regarding the total field strength, we note in particular
the occurrence of magnetic traps of $\vec{B}^{\mathrm{I,C}}$, regions of
markedly reduced field strength in which charged particles may gyrate for
extended periods of time before they can escape again. (See Florinski et al.
(2019) for an investigation of such traps for the heliospheric case based on
an extension of the Röken et al. (2015) field model.)

While the compressible magnetic field itself is fully described by its
vectorial components, the compressible analog to the scalar advection
processes discussed in Section 4.4 requires a correspondingly modified
renormalized second Euler potential $\beta_{\mathrm{ren}}^{\mathrm{C}}$. A
derivation of the explicit form of the latter is presented in Appendix D,
allowing applications such as those presented in Figure 4 to be obtained also
for compressible flow with relative ease. As a further example of its
continuing usefulness, we note that, for constant Cartesian components
$(B_{x0},0,B_{z0})$ of a planar homogeneous field, the general field line
equation

$\frac{{\rm d}x}{B_{x}}=\frac{{\rm d}z}{B_{z}}$ (99)

can trivially be integrated to the condition

$B_{z0}\,x-B_{x0}\,z=\mathrm{const.}$ (100)

describing the associated straight field lines. However, since the
transformation
$(\rho,z)\leftrightarrow\left(a(\rho,z),{\beta_{\mathrm{ren}}}(\rho,z)\right)$
discussed in Section 4.2 continues to be valid for the compressible case,
field lines of the latter can easily be drawn as contours of

$B_{z0}\ a\Big{(}\frac{x}{\cos\varphi},z\Big{)}\,\cos\varphi-B_{x0}\
\beta^{\mathrm{C}}_{\mathrm{ren}}\Big{(}\frac{x}{\cos\varphi},z\Big{)}$ (101)

where $x=\rho\,\cos\varphi$, $\varphi\in\\{0,\pi\\}$ is used rather than
$\rho$ to be able to access negative $x$ coordinates. The usefulness of this
method, which allows us to draw field lines without the otherwise mandatory
recourse to a numerical field line tracer, is exemplified in Figure 6.
Therein, contours of expression (101) are drawn using the compressible second
Euler potential of Equation (D9), evaluated for the first-order density
approximation (84), i.e. with $\nu_{1}=1-n_{\mathrm{sp}}$ and $\nu_{2}=0$.
With the corresponding field lines (derived by traditional numerical
integration) overplotted, and both adjusted to pass through the same set of
seed points, the excellent agreement of the results from both methods is
evident.

Although in this plot the first-order approximation is preferred over its
second-order counterpart for technical reasons (detailed in Appendix C), we
see that even this simpler first-order variant provides a realistic and
physically plausible model of a magnetic field draping around the obstacle
that is immersed in a compressible flow.

Figure 6:  Field lines of $\vec{B}^{\mathrm{C}}$ (compressible, Equations
(94)–(95), $m=1$, $\gamma=1$, solid, red) for the first-order density
approximation (84), constrasted with those of $\vec{B}^{\mathrm{I}}$
(incompressible, Equations (B1)–(B3), $m=0$, dashed, blue) in the $y=0$ plane.
The former are amended with the corresponding contours of
$\beta^{\mathrm{C}}_{\mathrm{ren}}$ (thick gray lines), and both clearly
coincide. Additionally, thin contour lines mark the regions at which the
respective field strengths of $\vec{B}^{\mathrm{C}}$ (red, dash-dotted) and
$\vec{B}^{\mathrm{I}}$ (blue, dotted) have decreased to $[100,80,60,40]$% of
their common boundary value of
$\|\vec{B}_{0}^{\mathrm{C}}\|=\|\vec{B}_{0}^{\mathrm{I}}\|=\sqrt{2}$.

##  6 Summary and conclusions

In this work, a (magneto)hydrodynamical model has been presented for laminar,
inviscid flow ahead of and around a solid obstacle of paraboloid form whose
axis of symmetry coincides with that of the incident flow. For the case of the
fluid being both incompressible and irrotational, compact formulas are derived
for the scalar flow potential and the stream function, and hence for the
velocity. Next, this velocity field was used to solve the induction equation
of ideal magnetohydrodynamics for the magnetic field, subject to a homogeneous
upstream boundary condition of arbitrary inclination. Again, the resulting
field components can be written in a compact and explicit form, and involve
neither integrals nor special functions. This model therefore constitutes the
first (and therefore so far also the only) known nontrivial axially symmetric
MHD solution for obstacle-induced flow that uses no integrals whatsoever.

Furthermore, a renormalization procedure was used to address the problem of
the infinite travel time spent by a fluid element moving from upstream
infinity to the vicinity of the obstacle. As a result, a simple one-to-one
mapping from undistorted (obstacle-free) to distorted coordinate space is
established that not only eases visualization, but also allows to describe the
profound deformation experienced by arbitrarily prescribed scalar or vector-
valued fields through the presence of the obstacle around which they are
advected. In one direction, this mapping is available as an explicit formula,
while in the opposite direction, recourse to a (quickly converging) numerical
algorithm becomes mandatory. When interpreting the central obstacle as the
magnetopause of a planet in a subsonic stellar wind, typical astrophysical
applications include the modeling of transients such as magnetic clouds or
local density enhancements embedded in the flow that is incident on the
magnetopause.

In a further step, the above model has been generalized first to arbitrary
aspect ratios and then from incompressible to mildly compressible (subsonic)
flow. The resulting spatial variations in density, which exhibit the expected
pile-up ahead of the stagnation point and a return to the undisturbed density
at large lateral distances, are consistent with conservation of mass and
approximately consistent with conservation of linear momentum along
streamlines. Even for this compressible extension, an exact solution of the
associated magnetic field frozen into the flow could be found at the expense
of more involved, albeit still explicitly defined and fully analytical
expressions. Already in the simpler, incompressible case, this property is
unique among known analytical solutions for magnetized flow around solid
obstacles, and therefore makes the present model a valuable tool not only in
the magnetospheric context, but also for qualitative investigations of similar
geometries such as the sphere/ellipsoid or the heliopause.

##  Appendix A Modulation factor for obstacle shapes with conic sections

We show that exact solutions of the irrotationality constraint (1) for the
particular ansatz (8) only exist if $\varepsilon\in\\{0,1\\}$. To this end, we
substitute the stream function

$\Psi_{f}=-\frac{\rho^{2}}{2}\biggl{(}1-\frac{1}{f(\rho,z)}\biggr{)}$ (A1)

into the equation

$\nabla\times\left[\nabla\times\left(\frac{\Psi_{f}}{\rho}\,\vec{e}_{\varphi}\right)\right]=\vec{0}\,,$
(A2)

which yields

$\partial_{\rho\rho}f-\frac{2}{f}\,(\partial_{\rho}f)^{2}+\frac{3}{\rho}\,\partial_{\rho}f=-\partial_{zz}f+\frac{2}{f}\,(\partial_{z}f)^{2}\,.$
(A3)

Then, using Equation (8) and $\eta=r+\varepsilon\,z$, we immediately obtain

$\left(2\,\eta+[\varepsilon^{2}-1]\,r\right)\biggl{[}\partial_{\eta\eta}g_{\varepsilon}(\eta)-\frac{2\,[\partial_{\eta}g_{\varepsilon}(\eta)]^{2}}{g_{\varepsilon}(\eta)}\biggr{]}+4\,\partial_{\eta}g_{\varepsilon}(\eta)=0\,.$
(A4)

This nonlinear equation for $g_{\varepsilon}$ can be reduced to a simpler
linear equation for $h_{\varepsilon}=1/g_{\varepsilon}$, namely

$\left(2\,\eta+[\varepsilon^{2}-1]\,r\right)\partial_{\eta\eta}h_{\varepsilon}(\eta)+4\,\partial_{\eta}h_{\varepsilon}(\eta)=0\,.$
(A5)

By analyzing the prefactor $2\,\eta+[\varepsilon^{2}-1]\,r$ of the second-
order contribution, we find that the parameter value $\varepsilon=0$, for
which $\eta=r$, results in the ordinary differential equation

$r\,\partial_{rr}h_{0}(r)+4\,\partial_{r}h_{0}(r)=0$ (A6)

with the general solution

$h_{0}(r)=b+\frac{c}{r^{3}}$ (A7)

for $b,c\in\mathbb{R}$, reproducing Equation (6) if $b=0$ and $c=1$. The
choice $\varepsilon=1$, on the other hand, leads to an ordinary differential
equation of the form

$\eta\,\partial_{\eta\eta}h_{1}(\eta)+2\,\partial_{\eta}h_{1}(\eta)=0$ (A8)

with the general solution

$h_{1}(r+z)=b+\frac{c}{r+z}\,,$ (A9)

which gives rise to Equation (10) for the same choices of $b$ and $c$. Note
that, since these constants of integration can be absorbed into the
normalizations of length and velocity, different choices for $b$ and $c$ will
not lead to qualitatively different solutions.

For the remaining parameter ranges $0<\varepsilon<1$ and $\varepsilon>1$, it
is obvious that Equation (A5) does not have any solutions, as the prefactor,
in addition to $\eta$, now contains the independent variable $r$. Accordingly,
general solutions $h_{\varepsilon}\in C^{2}(\mathbb{R},\mathbb{R})$ of the
irrotationality constraint only exist for the cases $\varepsilon=0$ and
$\varepsilon=1$. Employing a more general ansatz than the simple modulation
factor ansatz (A1) together with the conic section condition (8) may, however,
resolve this problem.

##  Appendix B Magnetic field components in other coordinate systems

To ease the application of our model, we summarize here, without derivation,
the components of our magnetic field solution (for incompressible flow and
homogeneous upstream boundary conditions) in cylindrical, Cartesian, and
spherical coordinates.

###  B.1 Cylindrical coordinates

This is simply the sum of Equations (3.1) and (35)–(37), written in terms of
Cartesian boundary components according to Equation (21):

$\displaystyle B_{\rho}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r+z}{r+z-1}}\left(1-\frac{r-z}{2\,(r+z)\,r}\right)\left[B_{x0}\cos\varphi+B_{y0}\sin\varphi\right]-\frac{1}{2\,r}\sqrt{\frac{r-z}{r+z}}\,B_{z0}$
(B1) $\displaystyle B_{\varphi}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r+z}{r+z-1}}\left[-B_{x0}\sin\varphi+B_{y0}\cos\varphi\right]$
(B2) $\displaystyle B_{z}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\,r}\,\sqrt{\frac{r-z}{r+z-1}}\left[B_{x0}\cos\varphi+B_{y0}\sin\varphi\right]+\left(1-\frac{1}{2\,r}\right)B_{z0}\,.$
(B3)

For an obstacle shape that has been linearly scaled using Equation (42), these
components are generalized to

$\displaystyle B^{\mathrm{sc}}_{\rho}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r_{\ell}\,\ell_{z}+z}{(r_{\ell}-1)\,\ell_{z}+z}}\left(1-\frac{r_{\ell}\,\ell_{z}-z}{2\,(r_{\ell}\,\ell_{z}+z)\,r_{\ell}}\right)\left[B_{x0}\cos\varphi+B_{y0}\sin\varphi\right]-\frac{\ell_{\rho}}{2\,r_{\ell}\,\ell_{z}}\sqrt{\frac{r_{\ell}\,\ell_{z}-z}{r_{\ell}\,\ell_{z}+z}}\,B_{z0}$
(B4) $\displaystyle B^{\mathrm{sc}}_{\varphi}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r_{\ell}\,\ell_{z}+z}{(r_{\ell}-1)\,\ell_{z}+z}}\left[-B_{x0}\sin\varphi+B_{y0}\cos\varphi\right]$
(B5) $\displaystyle B^{\mathrm{sc}}_{z}$ $\displaystyle=$
$\displaystyle-\frac{\ell_{z}}{2\,r_{\ell}\,\ell_{\rho}}\,\sqrt{\frac{r_{\ell}\,\ell_{z}-z}{(r_{\ell}-1)\,\ell_{z}+z}}\left[B_{x0}\cos\varphi+B_{y0}\sin\varphi\right]+\left(1-\frac{1}{2\,r_{\ell}}\right)B_{z0}$
(B6)

with
$r_{\ell}=\|\vec{r}_{\ell}\|=\sqrt{(\rho/\rho_{\ell})^{2}+(z/z_{\ell})^{2}}$.

###  B.2 Cartesian coordinates

These are found from the cylindrical components (B1)–(B3) in the usual way
through

$(B_{x},B_{y},B_{z})=(B_{\rho}\cos\varphi-
B_{\varphi}\sin\varphi,B_{\rho}\sin\varphi+B_{\varphi}\cos\varphi,B_{z})$ (B7)

and the substitutions $\cos\varphi=x/\rho$, $\sin\varphi=y/\rho$, and
$\rho=\sqrt{x^{2}+y^{2}}=\sqrt{r^{2}-z^{2}}$ as

$\displaystyle B_{x}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r+z}{r+z-1}}\left(B_{x0}-\frac{(x\,B_{x0}+y\,B_{y0})\,x}{2\,(r+z)^{2}\,r}\right)-\frac{x\,B_{z0}}{2\,(r+z)\,r}$
(B8) $\displaystyle B_{y}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r+z}{r+z-1}}\left(B_{y0}-\frac{(x\,B_{x0}+y\,B_{y0})\,y}{2\,(r+z)^{2}\,r}\right)-\frac{y\,B_{z0}}{2\,(r+z)\,r}$
(B9) $\displaystyle B_{z}$ $\displaystyle=$
$\displaystyle-\frac{x\,B_{x0}+y\,B_{y0}}{2\,r\,\sqrt{(r+z)(r+z-1)}}+\left(1-\frac{1}{2\,r}\right)B_{z0}$
(B10)

after a series of elementary algebraic manipulations. Furthermore, when
subjected to the triaxial scaling law of Equation (48), they become

$\displaystyle B_{x}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r_{\ell}\,\ell_{z}+z}{(r_{\ell}-1)\,\ell_{z}+z}}\left(B_{x0}-\frac{\left[\left(x/\ell_{x}^{2}\right)B_{x0}+\left(y/\ell_{y}^{2}\right)B_{y0}\right]x}{2\,(r_{\ell}+z/\ell_{z})^{2}\,r_{\ell}}\right)-\frac{x\,B_{z0}}{2\,(r_{\ell}\,\ell_{z}+z)\,r_{\ell}}$
(B11) $\displaystyle B_{y}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r_{\ell}\,\ell_{z}+z}{(r_{\ell}-1)\,\ell_{z}+z}}\left(B_{y0}-\frac{\left[\left(x/\ell_{x}^{2}\right)B_{x0}+\left(y/\ell_{y}^{2}\right)B_{y0}\right]y}{2\,(r_{\ell}+z/\ell_{z})^{2}\,r_{\ell}}\right)-\frac{y\,B_{z0}}{2\,(r_{\ell}\,\ell_{z}+z)\,r_{\ell}}$
(B12) $\displaystyle B_{z}$ $\displaystyle=$
$\displaystyle-\frac{\left[\left(x/\ell_{x}^{2}\right)B_{x0}+\left(y/\ell_{y}^{2}\right)B_{y0}\right]\ell_{z}^{2}}{2\,r_{\ell}\,\sqrt{\left(r_{\ell}\,\ell_{z}+z\right)\left([r_{\ell}-1]\ell_{z}+z\right)}}+\left(1-\frac{1}{2\,r_{\ell}}\right)B_{z0}$
(B13)

with $r_{\ell}$ given by Equation (50). Choosing
$\ell_{x}=\ell_{y}=\ell_{\rho}$ reproduces the symmetrically scaled case
(B4)–(B6), which in turn degenerates into the standard case (B8)–(B10) of
potential flow if $\ell_{\rho}=\ell_{z}$.

###  B.3 Spherical coordinates

For the sake of brevity, we content ourselves with the potential flow case for
these coordinates, which is obtained from Equations (B1)–(B3) through the
standard substitutions
$(\rho,\varphi,z)=(r\sin\vartheta,\varphi,r\cos\vartheta)$, resulting in

$\displaystyle B_{r}$ $\displaystyle=$
$\displaystyle\frac{[2\,(1+\cos\vartheta)\,r-1]\,\sin\vartheta}{2\,\sqrt{(1+\cos\vartheta)\,(r+r\cos\vartheta-1)\,r}}\,[B_{x0}\cos\varphi+B_{y0}\sin\varphi]+\left(\cos\vartheta-\frac{1}{2\,r}\right)B_{z0}$
(B14) $\displaystyle B_{\vartheta}$ $\displaystyle=$
$\displaystyle\frac{[2\,(1+\cos\vartheta)\,r-1]\,\cos\vartheta+1}{2\,\sqrt{(1+\cos\vartheta)\,(r+r\cos\vartheta-1)\,r}}\,[B_{x0}\cos\varphi+B_{y0}\sin\varphi]-\sin\vartheta\left(1-\frac{1}{2\,(1+\cos\vartheta)\,r}\right)B_{z0}$
(B15) $\displaystyle B_{\varphi}$ $\displaystyle=$
$\displaystyle\left(1-\frac{1}{(1+\cos\vartheta)\,r}\right)^{-1/2}[-B_{x0}\sin\varphi+B_{y0}\cos\varphi]$
(B16)

through the corresponding transformation

$(B_{r},B_{\vartheta},B_{\varphi})=(B_{\rho}\sin\vartheta+B_{z}\cos\vartheta,B_{\rho}\cos\vartheta-
B_{z}\sin\vartheta,B_{\varphi})\,.$ (B17)

##  Appendix C Determination of the functions $\vec{K}_{1}$ and $\vec{K}_{2}$

We present explicit analytic expressions for the functions $\vec{K}_{1}$ and
$\vec{K}_{2}$ defined in Equation (91). For many applications, it will likely
be sufficient to content oneself with the first-order density approximation
(84), whose magnetic field merely requires the simpler $\vec{K}_{1}$ but not
the much more involved $\vec{K}_{2}$ (see below). Besides its more complex
formulas, a further potential complication with the second-order case may
arise in numerical applications because several of the terms in $\vec{K}_{2}$
individually diverge in the vicinity of the $(a=0)$ streamline, and only
produce values of order unity in their total sum through mutual cancellation.
This is the main reason why the second-order magnetic field formulas, which we
only provide here for the sake of completeness, have not been used in Figure
6.

Suppressing the limit for the time being and using the Leibniz integral rule,
we find that Equation (91) may be written in the form

$\displaystyle\nabla\int_{\rho_{0}}^{\rho}\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\,{\rm
d}\rho^{\prime}=\int_{\rho_{0}}^{\rho}\nabla\biggl{(}\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\biggr{)}\,{\rm d}\rho^{\prime}+\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\,\nabla\rho-\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\bigg{|}_{\rho^{\prime}=\rho_{0}}\,\nabla\rho_{0}\,,$
(C1)

where

$\displaystyle\rho_{0}=\rho_{0}(\rho,z;z_{0})=\sqrt{\Bigl{[}\sqrt{(r-z)\,(r+z-1)+(z_{0}-1/2)^{2}}+1/2\Bigr{]}^{2}-z_{0}^{2}}\,.$
(C2)

We point out that the latter relation is a direct result of the condition
$a(\rho_{0},z_{0})=a(\rho,z)$. Evaluating the gradients on the right hand side
of (C1), we obtain

$\displaystyle\nabla\int_{\rho_{0}}^{\rho}\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\,{\rm d}\rho^{\prime}=$
$\displaystyle\int_{\rho_{0}}^{\rho}\partial_{a}\biggl{(}\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\biggr{)}\,{\rm
d}\rho^{\prime}\times\biggl{[}\frac{\partial
a}{\partial\rho}\,\vec{e}_{\rho}+\frac{\partial a}{\partial
z}\,\vec{e}_{z}\biggr{]}+\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\,\vec{e}_{\rho}-\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\bigg{|}_{\rho^{\prime}=\rho_{0}}\,\biggl{[}\frac{\partial\rho_{0}}{\partial\rho}\,\vec{e}_{\rho}+\frac{\partial\rho_{0}}{\partial
z}\,\vec{e}_{z}\biggr{]}\,.$ (C3)

With

$\displaystyle\frac{{\cal
A}}{u_{\rho}^{\mathrm{I}}}=\rho+\frac{a^{4}}{(\rho^{2}-a^{2})^{2}\,\rho}$ (C4)

and therefore

$\displaystyle\partial_{a}\biggl{(}\frac{{\cal
A}}{u_{\rho}^{\mathrm{I}}}\biggr{)}=\frac{4\,a^{3}\rho}{(\rho^{2}-a^{2})^{3}}\,,$
(C5)

the remaining integral on the right hand side for the case $k=1$ readily
yields

$\displaystyle\int_{\rho_{0}}^{\rho}\partial_{a}\biggl{(}\frac{{\cal
A}}{u_{\rho}^{\mathrm{I}}}\biggr{)}\,{\rm
d}\rho^{\prime}=\frac{(\rho^{2}-\rho_{0}^{2})\,(\rho^{2}+\rho_{0}^{2}-2\,a^{2})\,a^{3}}{(\rho^{2}-a^{2})^{2}\,(\rho_{0}^{2}-a^{2})^{2}}\,.$
(C6)

Substituting this integral as well as Equation (29) into Equation (C3) for
$k=1$, the expression for the function $\vec{K}_{1}$ results, after having
performed the limit, in

$\displaystyle\vec{K}_{1}=\frac{1}{2\,r}\,\sqrt{\frac{r-z}{r+z}}\,\,\vec{e}_{\rho}+\biggl{(}\frac{1}{2\,r}-1\biggr{)}\,\vec{e}_{z}\,.$
(C7)

For the first-order approximation, this is all that is needed to evaluate the
corresponding compressible magnetic field. To obtain the expression for
$\vec{K}_{2}$ that is required for the second-order case, we use

$\displaystyle\frac{{\cal
A}^{2}}{u_{\rho}^{\mathrm{I}}}=\frac{\bigl{[}(\rho^{2}-a^{2})^{2}\,\rho^{2}+a^{4}\bigr{]}^{2}}{\bigl{[}\rho^{2}+(\rho^{2}-a^{2})^{2}\bigr{]}\,(\rho^{2}-a^{2})^{2}\,\rho^{3}}$
(C8)

and

$\displaystyle\partial_{a}\biggl{(}\frac{{\cal
A}^{2}}{u_{\rho}^{\mathrm{I}}}\biggr{)}=\frac{4\,\bigl{[}(\rho^{2}-a^{2})^{2}\,\rho^{2}+a^{4}\bigr{]}\,\bigl{[}(2\,\rho^{2}-a^{2})\,a^{2}-(\rho^{2}-a^{2})^{2}\,(\rho^{2}-2\,a^{2})\bigr{]}\,a}{\bigl{[}\rho^{2}+(\rho^{2}-a^{2})^{2}\bigr{]}^{2}\,(\rho^{2}-a^{2})^{3}\,\rho}\,,$
(C9)

such that the integral on the right hand side of Equation (C3) for $k=2$
becomes

$\displaystyle\int_{\rho_{0}}^{\rho}\partial_{a}\biggl{(}\frac{{\cal
A}^{2}}{u_{\rho}^{\mathrm{I}}}\biggr{)}\,{\rm d}\rho^{\prime}$
$\displaystyle=$
$\displaystyle-\frac{2}{a}\,\biggl{(}\frac{1}{a^{4}}-\frac{4}{a^{2}}+2\biggr{)}\,Q-\frac{2\,(\rho^{2}-\rho_{0}^{2})\,\bigl{[}\rho^{2}\,\rho_{0}^{2}+(1-3\,a^{2})\,(\rho^{2}+\rho_{0}^{2})-5\,(1-a^{2})\,a^{2}+1\bigr{]}}{\bigl{[}\rho^{2}+(\rho^{2}-a^{2})^{2}\bigr{]}\,\bigl{[}\rho_{0}^{2}+(\rho_{0}^{2}-a^{2})^{2}\bigr{]}\,a}-$
$\displaystyle-\frac{(\rho^{2}-\rho_{0}^{2})\,\bigl{[}2\,\rho^{2}\,\rho_{0}^{2}-3\,(\rho^{2}+\rho_{0}^{2})\,a^{2}+4\,a^{4}\bigr{]}\,a}{(\rho^{2}-a^{2})^{2}\,(\rho_{0}^{2}-a^{2})^{2}}+\frac{2\,a^{2}-1}{a^{5}}\,\ln{\Biggl{(}\frac{\bigl{[}\rho^{2}+(\rho^{2}-a^{2})^{2}\bigr{]}\,\rho_{0}^{4}}{\bigl{[}\rho_{0}^{2}+(\rho_{0}^{2}-a^{2})^{2}\bigr{]}\,\rho^{4}}\Biggr{)}}\,,$

where

$\displaystyle
Q=Q(\rho,a;\rho_{0}):=\left\\{\begin{array}[]{lcl}\displaystyle\displaystyle\frac{1}{\sqrt{1-4\,a^{2}}}\,\mathrm{arccoth}\Biggl{(}\frac{\rho^{2}+\rho_{0}^{2}+2\,(\rho^{2}-a^{2})\,(\rho_{0}^{2}-a^{2})}{\sqrt{1-4\,a^{2}}\,(\rho^{2}-\rho_{0}^{2})}\Biggr{)}&\
:&\ a<1/2\\\\[18.49411pt]
\displaystyle\frac{8\,(\rho^{2}-\rho_{0}^{2})}{(1+4\,\rho^{2})\,(1+4\,\rho_{0}^{2})}&\
:&\ a=1/2\\\\[14.22636pt]
\displaystyle\frac{1}{\sqrt{4\,a^{2}-1}}\,\mathrm{arccot}\Biggl{(}\frac{\rho^{2}+\rho_{0}^{2}+2\,(\rho^{2}-a^{2})\,(\rho_{0}^{2}-a^{2})}{\sqrt{4\,a^{2}-1}\,(\rho^{2}-\rho_{0}^{2})}\Biggr{)}&:&\
a>1/2\,.\end{array}\right.$ (C14)

Again substituting the integral as well as Equations (29) and (C2) into
Equation (C3) for $k=2$ and performing the limit, the function $\vec{K}_{2}$
computes to

$\displaystyle\vec{K}_{2}=\frac{1}{2\,r}\,\sqrt{\frac{r+z}{r-z}}\,\biggl{[}2\,(r-1)+\frac{1}{r+z}\biggr{]}^{2}\,\vec{e}_{\rho}+\biggl{[}\biggl{(}1-\frac{1}{2\,r}\biggr{)}\,\sqrt{\frac{r+z}{r-z}}\
\vec{e}_{\rho}+\frac{1}{2\,r}\,\vec{e}_{z}\biggr{]}\times\biggl{[}2\,\biggl{(}1-\frac{\rho^{2}}{a^{2}}\biggr{)}\,\biggl{(}\frac{1}{a^{4}}-\frac{4}{a^{2}}+2\biggr{)}\times\biggr{.}$
(C15)
$\displaystyle\biggl{.}\times\,Q(\rho,a;a)+3-2\,r-\frac{r^{2}+2\,z^{2}+3\,r\,z-4\,r-5\,z+4+1/(r-z)}{(r+z-1)^{2}\,r}+\frac{(2\,a^{2}-1)\,(r-z)}{a^{6}}\,\ln{\biggl{(}\frac{2\,(r+z-1)\,r}{(r+z)^{2}}\biggr{)}}\biggr{]}\,.\hskip
34.14322pt$

##  Appendix D The renormalized second Euler potential for the compressible
case

For the determination of the second renormalized Euler potential

$\beta^{\mathrm{C}}_{\mathrm{ren}}=\lim_{z_{0}\rightarrow\infty}\left[-\int_{\rho_{0}}^{\rho}\frac{{\rm
d}\rho^{\prime}}{u_{\rho}^{\mathrm{C}}|_{a}}+\int_{z_{0}}^{z_{\mathrm{a}}}\frac{{\rm
d}z}{u_{z}^{\mathrm{C}}|_{\rho=0}}\right]=\sum_{k=0}^{2}\nu_{k}\lim_{z_{0}\rightarrow\infty}\big{(}\hat{\beta}_{k}+\Delta\beta_{k}\big{)}\,,$
(D1)

in which

$\hat{\beta}_{k}:=-\int_{\rho_{0}}^{\rho}\left.\frac{{\cal
A}^{k}}{u_{\rho}^{\mathrm{I}}}\right|_{a}{\rm
d}\rho^{\prime}\quad\mbox{and}\quad\Delta\beta_{k}:=\int_{z_{0}}^{z_{\mathrm{a}}}\left.\frac{{\cal
A}^{k}}{u_{z}^{\mathrm{I}}}\right|_{\rho=0}{\rm d}z$ (D2)

with $\rho_{0}$ given by Equation (C2), we first compute explicit expressions
for the quantities defined in Equation (D2), yielding

$\displaystyle\hat{\beta}_{1}$ $\displaystyle=$
$\displaystyle-\int_{\rho_{0}}^{\rho}\left(\rho^{\prime}+\frac{a^{4}}{({\rho^{\prime}}^{2}-a^{2})^{2}\,\rho^{\prime}}\right){\rm
d}\rho^{\prime}=\frac{1}{2}\Biggl{[}-{\rho^{\prime}}^{2}+\frac{a^{2}}{{\rho^{\prime}}^{2}-a^{2}}-\ln\biggl{(}\frac{{\rho^{\prime}}^{2}}{{\rho^{\prime}}^{2}-a^{2}}\biggr{)}\Biggr{]}\Bigg{|}_{\rho_{0}}^{\rho}$
(D3) $\displaystyle\hat{\beta}_{2}$ $\displaystyle=$
$\displaystyle-\int_{\rho_{0}}^{\rho}\left(\frac{\big{[}a^{4}+({\rho^{\prime}}^{2}-a^{2})^{2}\,{\rho^{\prime}}^{2}\big{]}^{2}}{{\big{[}{\rho^{\prime}}^{2}+({\rho^{\prime}}^{2}-a^{2})^{2}\big{]}\,({\rho^{\prime}}^{2}-a^{2})^{2}\,\rho^{\prime}}^{3}}\right){\rm
d}\rho^{\prime}=\biggl{(}\frac{1}{2\,a^{4}}-\frac{3}{a^{2}}+4\biggr{)}\,Q(\rho,a;\rho_{0})+$
(D4)
$\displaystyle+\frac{1}{2}\left[\left.-{\rho^{\prime}}^{2}+\frac{1}{{\rho^{\prime}}^{2}}+\frac{a^{2}}{{\rho^{\prime}}^{2}-a^{2}}+\ln\biggl{(}\frac{({\rho^{\prime}}^{2}-a^{2})^{3}}{{\rho^{\prime}}^{4}}\biggr{)}-\frac{1}{a^{2}}\,\biggl{(}\frac{1}{2\,a^{2}}-2\biggr{)}\,\ln\biggl{(}\frac{{\rho^{\prime}}^{2}+({\rho^{\prime}}^{2}-a^{2})^{2}}{{\rho^{\prime}}^{4}}\biggr{)}\right]\right|_{\rho_{0}}^{\rho}\,,$

where the function $Q$ is defined in (C14), and

$\displaystyle\Delta\beta_{1}$ $\displaystyle=$
$\displaystyle\int_{z_{0}}^{z_{\mathrm{a}}}\left(\frac{1}{2\,z}-1\right){\rm
d}z=\left.\left[\frac{\ln(z)}{2}-z\right]\right|_{z_{0}}^{z_{\mathrm{a}}}$
(D5) $\displaystyle\Delta\beta_{2}$ $\displaystyle=$
$\displaystyle\int_{z_{0}}^{z_{\mathrm{a}}}\left(\frac{1}{2\,z}-1\right)^{3}{\rm
d}z=\left.\left[\frac{3\,\ln(z)}{2}-z+\frac{3}{4\,z}-\frac{1}{16\,z^{2}}\right]\right|_{z_{0}}^{z_{\mathrm{a}}}\,.$
(D6)

The quantities $\hat{\beta}_{0}$ and $\Delta\beta_{0}$ can be found in
Sections 3.2 and 3.3. Then, the limits
$L_{k}:=\lim_{z_{0}\rightarrow\infty}\big{(}\hat{\beta}_{k}+\Delta\beta_{k}\big{)}$
read

$\displaystyle L_{1}$ $\displaystyle=$ $\displaystyle
z-\frac{1}{2}\ln\left(\frac{r+z}{2}\right)$ (D7) $\displaystyle L_{2}$
$\displaystyle=$ $\displaystyle
z-\frac{1}{2}\,\ln\biggl{(}\frac{(r+z)^{2}\,(r+z-1)}{8}\biggr{)}-\biggl{(}\frac{1}{2\,a^{4}}-\frac{3}{a^{2}}+4\biggr{)}\,Q(\rho,a;a)-\frac{1}{2\,(r+z-1)\,\rho^{2}}+$
(D8)
$\displaystyle+\frac{(r-1)\,r-(z-1/2)^{2}}{(r-z)^{2}\,(r+z-1)^{2}}\,\ln\left(\frac{2\,(r+z-1)\,r}{(r+z)^{2}}\right)+\frac{11}{16}\,,$

with each ${z_{\mathrm{a}}}$ adjusted such that $L_{k}|_{\rho=0,z=1}=1$ (as
done previously for $L_{0}$ in Section 3.3). Substituting these expressions
together with $\nu_{0}=n_{\mathrm{sp}}$ and
$L_{0}=\beta^{\mathrm{I}}_{\mathrm{ren}}$ (cf. Equation (41)) into Equation
(D1), we are finally lead to

$\displaystyle\beta^{\mathrm{C}}_{\mathrm{ren}}$ $\displaystyle=$
$\displaystyle\underbrace{(n_{\mathrm{sp}}+\nu_{1}+\nu_{2})}_{=1}\,z+\frac{n_{\mathrm{sp}}\ln(r+z-1)}{2}-\frac{\nu_{1}}{2}\,\ln\left(\frac{r+z}{2}\right)+\nu_{2}\,\Biggl{[}-\frac{1}{2}\,\ln\biggl{(}\frac{(r+z)^{2}\,(r+z-1)}{8}\biggr{)}-\Biggr{.}$
(D9)
$\displaystyle\Biggl{.}-\biggl{(}\frac{1}{2\,a^{4}}-\frac{3}{a^{2}}+4\biggr{)}\,Q(\rho,a;a)-\frac{1}{2\,(r+z-1)\,\rho^{2}}+\frac{(r-1)\,r-(z-1/2)^{2}}{(r-z)^{2}\,(r+z-1)^{2}}\,\ln\left(\frac{2\,(r+z-1)\,r}{(r+z)^{2}}\right)+\frac{11}{16}\Biggr{]}\
.$

Just as with the $\vec{K}_{1,2}$ terms for the magnetic field components, we
again see that the first-order case ($\nu_{2}=0$) is quite compact and void of
apparent practical difficulties, while the second-order case ($\nu_{2}\neq 0$)
is not only more involved but has several terms which individually diverge
near $a=0$ and only cancel when written as a sum. This is most clearly seen in
the final two terms in the square bracket of Equation (D9) diverging to
$\pm\infty$ as $a\rightarrow 0$.

## References

  * Batchelor (2000) Batchelor, G. 2000, An Introduction to Fluid Dynamics, Cambridge Mathematical Library (Cambridge University Press)
  * Cauchy (1816) Cauchy, A.-L. 1816, Théorie de la propagation des ondes à la surface d’un fluide pesant d’une profondeur indéfinie (Académie Royale des Sciences)
  * Combi et al. (2002) Combi, M. R., Gombosi, T. I., & Kabin, K. 2002, Geophys. Monogr. Ser. 130, 151, doi: 10.1029/130GM10
  * Davis & Werle (1972) Davis, R. T., & Werle, M. J. 1972, AIAA Journal, 10, 1224, doi: 10.2514/3.50354
  * Dursi & Pfrommer (2008) Dursi, L. J., & Pfrommer, C. 2008, Astrophys. J., 677, 993, doi: 10.1086/529371
  * Florinski et al. (2019) Florinski, V. A., Kleimann, J., & Ghanbari, K. 2019, in AGU Fall Meeting Abstracts, Vol. 2019, SH51B–05
  * Isenberg et al. (2015) Isenberg, P. A., Forbes, T. G., & Möbius, E. 2015, Astrophys. J., 805, 153, doi: 10.1088/0004-637X/805/2/153
  * Kleimann et al. (2017) Kleimann, J., Röken, C., & Fichtner, H. 2017, Astrophys. J., 838, 75, doi: 10.3847/1538-4357/aa5f15
  * Kleimann et al. (2016) Kleimann, J., Röken, C., Fichtner, H., & Heerikhuisen, J. 2016, ApJ, 816, 29, doi: 10.3847/0004-637X/816/1/29
  * Kleimann et al. (2022) Kleimann, J., Dialynas, K., Fraternale, F., et al. 2022, Space Sci. Rev, 218, 36, doi: 10.1007/s11214-022-00902-6
  * Kobel & Flückiger (1994) Kobel, E., & Flückiger, E. O. 1994, J. Geophys. Res., 99, 23617, doi: 10.1029/94JA01778
  * Kotova et al. (2020) Kotova, G., Verigin, M., Gomboshi, T., & Kabin, K. 2020, Solar-Terrestrial Physics, 6, 44, doi: 10.12737/stp-64202007
  * Miller (1969) Miller, D. R. 1969, Proc. Camb. Philos. Soc., 65, 285, doi: 10.1017/S0305004100044248
  * Miller (1971) —. 1971, Proc. Camb. Philos. Soc., 70, 123, doi: 10.1017/S0305004100049732
  * Milne-Thomson (1968) Milne-Thomson, L. M. 1968, Theoretical Hydrodynamics, 5th edn. (MacMillan Press Ltd, London)
  * Nabert et al. (2013) Nabert, C., Glassmeier, K. H., & Plaschke, F. 2013, Ann. Geophys., 31, 419, doi: 10.5194/angeo-31-419-2013
  * Petrinec & Russell (1997) Petrinec, S. M., & Russell, C. T. 1997, Space Sci. Rev., 79, 757, doi: 10.1023/A:1004938724300
  * Richardson (2002) Richardson, J. D. 2002, Planet. Space Sci., 50, 503, doi: 10.1016/S0032-0633(02)00029-6
  * Röken et al. (2015) Röken, C., Kleimann, J., & Fichtner, H. 2015, Astrophys. J., 805, 173, doi: 10.1088/0004-637X/805/2/173
  * Romashets & Vandas (2019) Romashets, E. P., & Vandas, M. 2019, J. Geophys. Res. (Space Physics), 124, 2697, doi: 10.1029/2018JA026006
  * Spreiter & Stahara (1995) Spreiter, J. R., & Stahara, S. S. 1995, Adv. Space Res., 15, 433, doi: 10.1016/0273-1177(94)00128-N
  * Stern (1966) Stern, D. P. 1966, Space Sci. Rev., 6, 147, doi: 10.1007/BF00222592