Two-jet structure of the flow produced by magnetized hypersonic spherical source into the steady unmagnetized medium

In the present paper we consider the steady flow produced by the hypersonic magnetized spherical source into the steady unmagnetized medium. The magnetic field of the source is considered purely azimuthal. The basic dimensionless parameter of the problem is Alfvén number at the inflow boundary at the source equatorial plane, which we denote as ε. We first review the numerical results, that was already presented in our recent paper (Golikov et al, 2017). According to our numerical results there should be a tangential discontinuity dividing the flow and the ambient medium. In the present paper we analytically derive the asymptotic behaviour of the tangential discontinuity radius as epsilon approaches zero.


Introduction
An interest for the flow the present paper is devoted to, was inspired by the problem of stellar wind/interstellar medium interaction -an astrophysical problem, that rose up about half a century ago. We would like to briefly introduce the reader to the physical background of the problem first.
It was realized in the middle of the last century that the Sun constantly produces a flow of fully ionized hydrogen plasma (see historical review in Chapter 1 of [1]). First direct measurements showed that this flow is highly supersonic at the Earth orbit [2]. This flow, nowadays known as the solar wind, proceeds up to approximately 90 astronomical units, where it undergoes a shock transition, because of the ambient interstellar gas pressure. The region of the solar wind/interstellar medium interaction is commonly referred as the heliospheric interface. The structure of this region is strongly affected by a number of physical factors; for details, see, for example, reviews and papers [3], [4].
It should be noted now, that not only the Sun has wind : a number of images of the structures similar to the heliospheric interface in the vicinity of other stars where observed in recent years. These structures are called astrospheres. The shape of an astrosphere may greatly differ from that of the heliosphere.
In the simplest case a star can be considered as a steady hypersonic spherical source of the fully ionized hydrogen plasma -the stellar wind. Following [5], [6], we treat the flow hydrodynamically in our study. The flow is spherically symmetric and is governed by the system of Euler equations together with the polytrope law, where we can omit the pressure gradient due to the hypersonic character of the flow. In this case, the solution is trivial: where R stands for the distance from the origin, γ is the polytrope index (usually equals to 5/3, that corresponds to monoatomic gases). However, the flow cannot proceed freely up to infinity due to existence of ambient interstellar gas. If we assume this gas to be at rest relatively to the star, then it causes the flow to overcome shock transition at certain distance R T S . Index "TS" stands for termination shock ; this term is traditionally used to distinguish this shock from the bow shock, that appears in ambient gas flow, if we assume it to be relatively moving with a supersonic velocity. The post-shock flow is subsonic, and if we treat it as incompressible, than the solution in the post-shock region is simple: where p ∞ is an ambient gas pressure. Let us now discuss the stellar magnetic field. The global magnetic field structure on the stellar surface corresponds to the magnetic field of a dipole. However, this magnetic field subsists not in vacuum, but in hydrogen plasma flow the star emits -the stellar wind. The stellar wind stretches the magnetic field lines, making their topology farther from the initial dipole field on the surface, but closer to the "split monopole" field (see Chapter 10 of [1]): such a field corresponds to the magnetic field of a (hypothetical) monopole, but with opposite polarities on opposite hemispheres. In our study the polarity of the field doesn't matter, since the magnetic field appears quadratically in equations of motion (see (7) and (8) below). Hence, we can assume the stellar magnetic field to be produced by a single "magnetic monopole".
There are no electric field in the stellar wind reference frame, due to its high electric conductivity. Since the magnetic dissipation is very low, the stellar magnetic field is completely "frozen" in the flow: It means that the magnetic flux through any fluid surface is constant. For purely radial field on the stellar surface this implies: where W mag stands for "magnetic power" of the source. Stars usually rotate, and this rotation affects their magnetic field. Let z-axis be the stellar rotation axis, and Ω be the stellar rotation frequency. Let us for now assume, that the magnetic field does not affect the flow (kinematic approximation). Applying the frozen-in condition (1) to the pre-shock flow, we get: Here θ is the angle counted from z-axis, ϕ is the azimuthal angle. This magnetic field configuration is nowadays referred as Parker spiral, after Parker (see [5]); see figure 1, left panel.
Kinematic approximation is not necessarily valid, and we generally have to add the magnetic force into the equation of gas motion: Magnitudes of magnetic force F mag and of gas inertia ρ(V · ∇)V are related to each other as squares of alfvén velocity and bulk velocity of the gas: The square root of this quantity is usually referred as Alfvén number : So, one can say that the Alfvén number characterizes the dynamical impact of the frozen-in magnetic field on the flow. In the pre-shock region for kinemtically-treated magnetic field, we get: Hence, as long as this quantity is small at the stellar surface, it remains so in the whole pre-shock region, and we can neglect the dynamical impact of the magnetic field there. However, it is not so in post-shock region; if we treat the magnetic field kinematically there, then: From this follows, that Al ∼ R 3 there, the magnetic force becomes very significant, and the kinematic approach is not valid in the post-shock region.
In the present paper we restrict ourselves to the case for which the "magnetic power" W mag is vanishing, whereas the rotation frequency Ω is huge, so that the azimuthal magnetic field is finite: In this case the azimuthal component of the magnetic field is much larger, than the radial one (see 2), and the magnetic force (3) is directed towards the stellar rotation axis. This force collimates the flow in front of the axis, forcing it to form a jet (see schematic picture of the flow: figure 1, right panel). The radius of the jet appears to be finite, in other words, the tangential discontinuity (TD) dividing the stellar wind and the interstellar medium is formed.
Since the magnetic field is purely azimuthal in our consideration, it is convenient to denote the proportionality factor of the magnetic field B and sin θ/R in the pre-shock region with a single symbol: We have to note, that a stellar wind flow with magnetically-driven jets was proposed earlier by Drake (see [7]). However, the results presented in this paper, both analytical and numerical, are questionable. We would like to elaborate the main points of this paper in a tidier manner.

Mathematical formulation of the problem
As we have already noted in the Introduction, we shall omit the radial component of the magnetic field at the source surface. Then the magnetic field becomes purely azimuthal in the whole domain of the flow. From this follows, that there is no magnetic force acting in the azimuthal direction, and the flow becomes axisymmetric with no twist around the symmetry axis (z-axis).
The system of ideal MHD equations in cylindrical coordinates (z, r, ϕ) is then follows (see [8]): Here ρ, V z V r , p are the density, two components of velocity and pressure, respectively; B is the ϕ-component of the magnetic field (index ϕ is omitted since the two other components are assumed to be zero); γ is the ratio of specific heats. We shall also use spherical coordinates (R, θ, ϕ) when it is convenient; spherical and cylindrical coordinates are connected as follows: In order to set the supersonic inflow conditions at certain distance R = R in (index in stands for "inflow") we should specify the inflow velocity V in , density ρ in , pressure p in and the azimuthal component of the magnetic field that depends on angle θ counted from the z-axis as B = B in sin θ, where B in is a constant.
However, instead of specifying density ρ in , velocity V in , pressure p in and magnetic field B in sin θ at certain distance R in , it is more convenient to specify the mass loss rateṀ = in , the terminal velocity V term (i.e. the velocity at infinite distance from the origin if we propose free supersonic unmagnetized spherical outflow into vacuum), the entropy s = log(p/ρ γ ) and the magnetic field multiplied by the distance from the origin at the equatorial plane F B = B in R in (see (5) with θ = π/2) there. It is so, because as long as we neglect the dynamical impact of the magnetic field on the flow, these quantities do not change with distance from the origin, hence, the flow does not depend on the inner boundary distance R in . Whereas we do not neglect the magnetic field impact in the present consideration, the second set of quantities is still more convenient.
In the present paper we consider hypersonic outflow (M 1) at the inflow boundary R = R in , so the gas velocity equals to its terminal velocity V in = V term .
The outer boundary is a tangential discontinuity with an unknown shape and position. The total pressure p + B 2 /8π at this boundary equals to the ambient pressure p ∞ , and the normal component of the velocity at the discontinuity is zero: V n = 0. Notice that the purely azimuthal magnetic field is tangential to the outer boundary and the boundary condition B n = 0 is satisfied automatically. (B n is the normal component of magnetic field to the discontinuity.) To finish the formulation of the problem we have to set the boundary conditions far from the source in the jets (at infinity, z → ±∞). In the numerical model we assume the so-called soft boundary conditions there meaning that ∂/∂z = 0 for all quantities.

Dimensionless formulation
In order to formulate the problem in dimensionless form, we choose the distance to the termination shock for the purely gas-dynamic case (B = 0) as a characteristic distance: As the characteristic density and velocity we choose their values downstream from the termination shock for the purely gas-dynamic case: 1/2 , so that the form of the governing equations system (6)-(10) remains unchanged after rescaling. The pressure balance condition at tangential discontinuity then becomes: is obtained from the solution of purely gas-dynamic problem (zero magnetic field).
As we have noticed previously, the inflow boundary conditions are formulated in terms of stellar mass-loss rateṀ , stellar wind terminal velocity V term , and the constant characterizing the azimuthal magnetic field F B . In the hypersonic outflow case we do not have to specify the entropy s, since the pressure at the inflow boundary do not affect the flow. In dimensionless form these quantities are given as follows: where epsilon is a dimensionless parameter of the problem: The second equality shows that epsilon is an Alfvén number at the inflow boundary R = R in at the equatorial plane z = 0. From here and then we omit "hats" assuming all of the quantities to be dimensionless.

First integrals
First let us introduce a streamline function ψ, which is defined in axisymmetric case as follows: At the inner boundary sphere (R = R in ) ψ can be expressed as following: The problem formulated above has three first integrals along a streamline: The first integral (equation (13)) is the Bernoulli integral generalized for the case of MHD (see, e.g. [8]). Equation (14) is the integral that follows from eqs. (10) and (6). Equation (15) is the adiabatic condition.
The right-hand parts of these equations are constants along streamlines. The expressions for these constants have been easily deduced from the inner boundary conditions in their dimensionless form. Note, however, that the entropic integral (15) is not conserved at the TS. S(ψ; ε) could be explicitly derived downstream the TS only if the pre-shock solution is known, and the exact shape of the TS is given a priori. The constants in the other two integrals (eqs. (13), (14)) do not depend on the TS, though.
Note that if the flow at the inflow boundary is subalfvénic (i.e. if ε > 1), than the termination shock is absent, and the flow differs qualitatively from the one for which it is superalfvénic. We consider only the superalfvénic case in the present paper: ε < 1. As for the case of the Sun, ε 0.08.

Results of numerical solution: parametric study
In this section we review results of numerical solution of the full two-dimensional problem (6)-(10) formulated in Section 2, that were presented previously in our paper [9]. In order to get the solution we have used the two-dimensional shock-fitting Godunov [10] scheme with Chakravarthy-Osher TVD limiter [11]. The numerical scheme is identical to that employed in [4].
We have used flow-adapting computational grid. The grid was adapted to the tube-like structure of the considered problem; it has (60 × 200) cells in the region between the termination shock and the tangential discontinuity. Note that the grid has been constructed in such a way that cells are nearly squares everywhere, except the region far from the star, where the flow does not change much in the z-direction; hence there the cells are stretched along the z-axis.
The convergence of the result with respect to the spatial grid was verified via the cell bisection procedure. Also, in order to verify the quality of the numerical results, we have checked that the total mass flux through the outer boundary in the jet equals to the mass flux at the inner inflow boundary. We have also checked the conservation of the three first integrals (13) Figure 4 shows one-dimensional distributions of the flow parameters in the pre-shock region (top row), and on the r-and z-axis of the post-shock region (middle and bottom rows respectively), for ε = 0.05, 0.1, 0.25 , and 0.5. Note that on the latter figure the distance is normalized to the actual (i.e. ε-dependent) termination shock distance.
The distances to the TS and TD decrease with increasing parameter ε, meaning that the jet becomes more collimated. For ε = 0.25 the distance to the termination shock is 0.655 and the distance to the tangential discontinuity in the equatorial plane is 1.29. The jet radius is 0.851. For ε = 0.5 the distances are even smaller (bottom row of figure 3). It is interesting to note that for ε = 0.25 the termination shock is still very close to the spherically symmetric one (although, a slight deflection from a sphere could be seen in the figure). For ε = 0.5 the termination shock becomes distorted.
One can see from the plots of density that the flow between the TS and the TD is nearly incompressible in the case of small epsilon and the magnetic field grows linearly with the distance from the symmetry axis as it should follow from the integral (14) (see plots (a) and (d) on the top row of figure 3). At the same time for large epsilon the flow is sufficiently compressible (see plots (a) on middle and bottom rows of figure 3, and figure 4). Thus the behaviour of the magnetic field is not linear with respect to the distance from the symmetry axis (see plots (d) on middle and bottom rows of figure 3 and plots (b) of figures 4).
To understand the plasma and magnetic field distribution in figure 3 in detail it is very useful to plot (see plots (f) in figure 3) the r-projection of the sum of pressure and magnetic forces acting on the flow in equation (8): For all values of ε the maximum of the force (directed toward the axis of symmetry) in the subsonic region is in equatorial plane just after the termination shock because the magnetic field jumps significantly (by approximately the factor of four for γ = 5/3). Therefore the last term increases by the factor of 16. The magnetic force is partially compensated by the pressure gradient when the distance increases. F r remains negative in the entire region. The flow is decelerated in the r-direction and the jet oriented toward the z-axis is formed.
While for ε = 0.01 the role of the magnetic force is restricted to the subsonic region, for ε = 0.25 and ε = 0.5 the magnetic force is also significant in the supersonic region. Similar to the subsonic region, the magnetic force acts toward the axis of symmetry (z-axis). The streamlines are slightly deflected towards the axis. This leads to a relative increase of the number density in the pole direction and a slight relative decrease of it in the equatorial direction as compared to the weak magnetic field case (e.g. ε = 0.01). The described effect is seen in the isolines of density in plots (a) on middle and bottom row of figure 3 and especially in plot (a) on top row of figure  4. Therefore for the cases of ε = 0.25 and ε = 0.5 the pre-shock density has the maximum at the pole (r = 0). The maximum remains beyond the shock for the nearly spherical (as for ε = 0.25) or ellipsoidal (as for ε = 0.5) shock. Further from the shock the density still slightly increases due to negative F r force. Further small decrease of the density at the z-axis in the jet may be connected with the fact that the pressure gradient slightly overcomes the magnetic pressure. The slight pressure overcome may result in a small increase of the jet radius with z. Although the effect is small it can be recognized in bottom row of figure 3. It is also interesting to note that the pressure maximum at the z-axis (see plots (c) in bottom row of figure 3) corresponds to the minimum of the plasma velocity in this region, as it follows from the Bernoulli integral (see plots (b) in bottom row of figure 3 and the plot (c) in bottom row of figure 4). This local deceleration at the axis results in the larger velocities at the tangential discontinuity bounding the jet.
Since the velocity is larger at the TD, the Mach number is larger there. Plots (e) in figure 3 show the fast-magnetosonic Mach number M ms = V /a ms . It is also shown on one-dimensional plots (d) in figure 4. Note that fast magnetosonic wave speed is the only meaningful disturbance propagation speed in our case of axial symmetry and purely azimuthal magnetic field: For ε = 0.5 the flow in the jet becomes transsonic (plot (e) on bottom row of figure 3 and plot (d) on bottom row of figure 4) in a sense that the fast-magnetosonic Mach number in the jet becomes very close to unity near the z-axis and even greater near the tangential discontinuity. For even larger values of ε (say, for ε = 0.8) the flow stops at some point on the symmetry axis downwind from the TS and turns around forming an eddy that causes sufficient numerical problems. The flow in the jet becomes completely supersonic. We do not show the results for this case because of possible computational uncertainties.
3.1. Radii of the tangential discontinuity at z = 0 and at z = ∞ Points on the figure 5 present the radii of the tangential discontinuity at z = 0 (crosses) and at z = ∞ (squares) with respect to parameter ε obtained from the numerical solution. We notice that in the large interval of the parameter (ε 0.1) both dependencies are very well fitted with  Figure 5. Radii of the tangential discontinuity at z = 0 (circles) and at z = ∞ (squares) with respect to parameter ε obtained from the numerical solution. Dashed lines correspond to the power-law fits for the numerical points (see (17)).
following power-laws: Note that exponents of both fits are the same and equal to −1/3. In the following section we explain it for the case of small epsilon without addressing our numerical results.

Analytical consideration
As we have already mentioned in Section 2.2, the entropic integral constant S is not conserved at the termination shock. In order to derive the expression for S downwind the TS we have to know the shape of the shock R T S (ψ) and the flow in the pre-shock region itself.
However, we are primarily interested in asymptotic behaviour of S as epsilon approaches zero. In the purely gas-dynamic case (ε = 0) the flow is known, and according to our problem rescaling (see Section 2.1) S in the post-shock region is expressed as follows: The expression for S(ψ; ε) for a non-zero epsilon should approach this one as epsilon approaches zero, i.e: At the moment we do not know what is the order of the first non-constant term in the series expansion. To denote it somehow, we shall introduce a new quantity -kappa: For now, kappa can be any positive number (if it is negative or zero, the O-term does not approach zero).

4.1.
Radius of the tangential discontinuity at the equatorial plane (z = 0) Let quantities with index "TD" and argument z denote the corresponding quantities at the tangential discontinuity at the height z, i.e.: The radius of the tangential discontinuity at the plane z = 0 (i.e. at the equatorial plane) can be obtained if we consider the streamline ψ = 0 (θ = π/2) that passes through the critical point at the TD.
At the critical point V = 0, therefore, the three first integrals (13), (14), (15) together with the boundary condition (11) provide four relations connecting five quantities. The five quantities are: 1) the dimensionless parameter ε, 2) the distance to the critical point r T D (0), and 3-5) the values of magnetic field, pressure and density B T D (0), p T D (0), ρ T D (0) at the critical point.
Using the four relations we can establish how r T D (0) depends on the parameter epsilon. The solution of the following system gives the dependency in the explicit form: where p ∞ denotes the right-hand side of equation (11). Using the expression (18) for entropy in the post-shock region, we find from the first equation of system (19), that: where is obtained from the solution of purely gas-dynamic problem (zero magnetic field case).
We do not know which of the O-terms of expression (20) is dominant as long as we do not know whether kappa is below unity or not. For convenience we shall introduce another quantity -bounded kappa -as follows:κ = κ, 0 < κ < 1, 1, κ ≥ 1 Than, the expression (20) becomes: Substituting (21) We immediately notice that ifκ = 1 (or, equivalently, κ ≥ 1, see the definition of bounded kappa), then the radius of the TD at z = 0 does not approach infinity, which is impossible, since the solution with vanishing epsilon must approach a solution for zero magnetic field; this solution does not have a tangential discontinuity. Consequently, κ =κ < 1 4.2. Radius of the tangential discontinuity in the jet (z → ∞) Reasonable assumptions for the flow far from the source, i.e. for z → ∞, are the following: 1) all parameters are independent of z coordinate, i.e. ∂/∂z = 0, and 2) V r component is negligible meaning that the inertial terms in equation (8) can be neglected. Therefore this equation can be written as d dr and the streamline function satisfies the following ODE: Two differential equations (23) and (24) together with the first integrals (13), (14), (15) form closed system of equations for ρ, V z , p, B, ψ as functions of r in the jet. This system can be reduced to the following system for ρ, V z and r 2 as functions of ψ: System (25) is the second-order system of ordinary differential equations (ODE) for the functions ρ(ψ), V z (ψ) and r 2 (ψ). The boundary conditions are posed as r 2 = 0 at z-axis (ψ = 1) and p + B 2 /8π = p ∞ at the tangential discontinuity (ψ = 0), where p and B are expressed in terms of r and ρ: p(ψ) = S(ψ; ε)ρ(ψ) γ , Using the expression (18) for entropy in the post-shock region, we shall find from the second equation of system (25), that: V 2 z (ψ; ε) = O(ε 2κ ); Then, it follows from the first equation of system (25), that: The radius of the TD in the jet is defined as r T D (∞) = r| ψ=0 ; hence,