© ESO, 2014

1. Introduction

Whenever a massive object passes through a rarefied medium, it draws surrounding matter toward it. As a result, this material creates an overdense wake behind the object that exerts its own gravitational pull, retarding the original motion. Such dynamical friction arises whether the medium consists of non-interacting point particles, e.g., a stellar cluster, or a continuum fluid, e.g., an interstellar cloud. Chandrasekhar (1943) provided the essential theory when the background is collisionless, and his solution has been extensively used in studies of both star clusters and galaxies. Interaction of a gravitating object with gas also occurs in a wide variety of situations. A partial list of topics and references includes: the interaction of planets and gaseous disks (Teyssandier et al. 2013); the orbital decay of common-envelope binaries (Ricker & Taam 2008); the settling of massive stars in dense molecular clouds (Chavarría et al. 2010); the coalescence of massive black holes in both isolated galactic nuclei (Narayan 2000) and colliding galaxies (Armitage & Natarajan 2005); and the heating of intracluster gas by infalling galaxies (El-Zant et al. 2004).

Despite the widespread occurrence of gaseous dynamical friction, there is still no generally accepted derivation of the force, even after 70 years of effort. The flow in the vicinity of the gravitating mass is complex both temporally and spatially, as many simulations have shown (see, e.g. Matsuda et al. 1987). Theorists seeking a fully analytic expression for the time-averaged force have employed various strategems to circumvent a detailed description of this region. In the course of their classic studies of stellar accretion, Bondi & Hoyle (1944) took the star to be traversing a zero-temperature gas. In their model, fluid elements follow hyperbolic orbits in the star’s reference frame and land in an infinitely thin, dense spindle behind the object. Bondi & Hoyle analyzed the transfer of linear momentum from the spindle to the star, and hence obtained the force. Dodd & McCrae (1952) similarly investigated this hypersonic limit, as did, much more recently, Cantó et al. (2011). Dokuchaev (1964) first treated a finite-temperature gas. He determined the force by integrating the total power emitted by the object in acoustic waves (see also Rephaeli & Salpeter 1980). Ruderman & Spiegel (1971) used an impulse approximation in the reference frame of the background gas. Finally, Ostriker (1999) calculated the force by integrating directly over the wake, whose density she obtained through a linear perturbation analysis.

Fig. 1 Sketch of mass and momentum flow. Surrounding the central gravitating body of mass M is a large, imaginary sphere. The thin curves represent streamlines of background gas. After entering the sphere, most gas simply exits again, while a small portion joins directly onto the mass. In addition, some gas temporarily forms an overdense wake, also sketched here. The wake tugs gravitationally on the mass, thereby imparting momentum to it (broad arrow). All gas entering the wake also leaves it, either joining the mass or exiting the sphere.

These researchers focused principally on the supersonic case, which is often the most interesting one astrophysically. That is, they took V > c_s , where V is the speed of the gravitating mass relative to the distant background, and c_s the sound speed in that gas. (In the hypersonic calculations of Bondi, Hoyle, and their successors, c_s was implicitly set to zero.) While their answers differed in detail, all agreed that the friction force in this regime varies as V^-2, with a coefficient that includes a Coulomb logarithm. This latter term also appears in the force derived by Chandrasekhar (1943), and arises from an integration in radius away from the mass. In all derivations, at least one of the integration limits is rather ill-defined.

The previous studies made two important, simplifying assumptions. First, they neglected any accretion of background gas by the moving object. Quantitatively, the assumption was that R ≫ r_acc , where R is the object’s physical radius and the accretion radius r_acc ≡ 2   G   M/V²is the distance from the mass M within which its gravity qualitatively alters the background flow. They further ignored R compared to the scale for spatial variations in the surrounding flow. If L denotes the latter scale, then these assumptions may be summarized as L ≫ R ≫ r_acc . Unfortunately, the inequality R ≫ r_accis often not satisfied. For example, a low-mass star moving through a cluster-producing molecular cloud has R ~ 10¹¹  cmand r_acc ~ 10¹⁶  cm . An analysis that covers this regime should treat the object as being point-like in all respects, allowing the possibility of mass accretion through infall.

This infall cannot occur via direct impact, since the geometrical cross section of the body is negligible by assumption. Instead, some fluid elements that initially miss the object are pulled back into it. Mass accretion is, in fact, closely related to dynamical friction. In the reference frame attached to the mass, the background gas flows by with a speed that approaches V far away. The steady-state accretion rate onto the object is simply the net influx of mass through any closed surface surrounding it. Similarly, the friction force is the net influx of linear momentum. Much closer to the mass, this momentum influx manifests itself as two distinct force components. One is the gravitational pull from the wake, as described earlier. A second component is the direct advection of linear momentum from any background gas that falls into the object. Any determination of the force by an asymptotic surface integration cannot tease apart these two contributions.

To clarify these ideas, Fig. 1 shows graphically the mass and momentum flow within the extended gas cloud surrounding the central, gravitating object. Gas enters the region via the dotted sphere shown in the sketch. This gas then follows one of three paths. Most of it leaves the sphere downstream of the gravitating mass, as shown by the two outermost streamlines (thin solid curves). Only a small fraction of the gas makes its way to the deep interior. Some of it accretes onto the mass, first missing it and then looping back. Other gas that misses more widely is temporarily slowed by the gravitational pull of the mass and forms an overdense wake, also sketched in the figure. In steady state, the wake cannot gain net mass, so all of the gas entering it also leaves, either by joining onto the mass or exiting the sphere directly. Meanwhile, the wake itself tugs on the mass gravitationally. The broad arrow in the figure depicts this second form of momentum input to the mass. It is possible to determine the force analytically by integrating the net momentum flux over the bounding sphere. Numerical simulations can obtain this same force by calculating the advective and gravitational components acting directly on the mass.

In a previous paper (Lee & Stahler 2011, hereafter Paper I), we did the analytic surface integration. We focused on the subsonic case, V < c_s , which traditionally has been less explored¹. Working in the reference frame whose origin is attached to the mass, we developed a perturbative method to analyze the small deviations of the background gas from a uniform, constant-density, flow. By integrating the perturbed variables over a large sphere to obtain the net momentum influx, we arrived at a surprisingly simple result for the force: F = Ṁ   V . Here, Ṁ is the mass accretion rate onto the moving object. To evaluate this quantity from first principles, one would have to follow the trajectories of fluid elements as they accrete onto the central mass, and ensure that, following turnaround, they smoothly cross the sonic transition. This transition occurs well inside the region of validity for our calculation. Lacking a fundamental theory for Ṁ, we adopted a variant of the interpolation formula of Bondi (1952) and thus found the force explicitly as a function of velocity. In this subsonic regime, our expression agrees reasonably well with past simulations (e.g., Ruffert 1996).

Here we apply the same technique to study the supersonic case. We quickly encounter the technical difficulty that some perturbation variables diverge, when expressed as functions of the angle from the background velocity vector. The true flow around the gravitating mass is, of course, well-defined everywhere, and the divergences simply indicate that the adopted series expansions fail at certain locations. Despite this mathematical inconvenience, we are able once again to derive the dynamical friction force through spatial integration of the linear momentum influx. The force has exactly the same form as previously: F = Ṁ   V . In the high-speed limit, we recover the V^-2 behavior found by others, but not the Coulomb logarithm.

In Sect. 2 below, we describe our solution strategy, and then formulate the problem using a convenient, non-dimensional scheme. Whenever material repeats that of Paper I, we abbreviate its presentation as much as possible. Section 3 analyzes the perturbed flow to first order only. In this approximation, we find there is no mass accretion or friction force, just as in the subsonic case. We extend the analysis to second order in Sect. 4, thus accounting for mass accretion. Here we also describe our method for calculating the flow numerically, and show sample results. Section 5 presents the analytic derivation of the force itself, and compares our expression to past simulations. Those found a greater force than we derive in the supersonic regime; we indicate possible causes for this discrepancy. Assuming our result to be correct, we present several representative applications. Finally, Sect. 6 compares our derivation to previous ones, and indicates directions for future work.

2. Method of solution

2.1. Physical assumptions

In the reference frame whose origin coincides with the mass, the background gas has speed V far from the object and a spatially uniform density, which we denote as ρ₀. We take the gas to be isothermal, with associated sound speed c_s. As depicted in Fig. 2, we will be working in a spherical coordinate system (r,θ), and will assume that the gas flow is axisymmetric about the polar (z-) axis, which is parallel to the asymptotic fluid velocity.

Fig. 2 Mathematical treatment of the flow. We erect a spherical coordinate system centered on the gravitating body. The gas is isothermal, and its velocity far upstream is β   cs  (β > 1). Indicated are the Mach angle θM and its supplement, ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\mathrm{\equiv }\mathit{\pi }\mathrm{-}{\mathit{\theta }}_{\mathit{M}}$ .

We neglect the self-gravity of the gas, and will be analyzing perturbations to the flow relatively far from the central mass. Specifically, our expansions are valid for r ≫ r_s , where ${\mathit{r}}_{\mathrm{s}}\mathrm{\equiv }\mathit{G}\hspace{0.17em}\mathit{M}\mathit{/}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$is the sonic radius. We assume that in the far-field region of interest, the flow is steady-state. Of course, a steady flow cannot be established over arbitrarily large distances, and one must judge, in each astrophysical situation, whether the assumption is justified. In numerical simulations, which we later discuss by way of comparison, the flow occurs within some fixed computational volume, at the center of which lies the mass. For the non-relativistic case of relevance, the force of gravity propagates instantaneously, but alterations in the flow density and velocity take time. In practice, a steady flow is indeed approached in many such experiments (e.g., Pogorelov et al. 2000). The observed decay of transients is instructive, but again is only broadly suggestive of what may occur in Nature.

Figure 2 singles out two special angles, both of which only occur in supersonic flow. The first is the familiar Mach angle, defined through the relation $\sin \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}\mathrm{\equiv }\mathrm{1}\mathit{/}\mathit{\beta ,}$(1)where β ≡ V/c_s . The second is the supplement of the first: ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\mathrm{\equiv }\mathit{\pi }\mathrm{-}{\mathit{\theta }}_{\mathit{M}}$. The figure of revolution swept out by the radius lying along θ = θ_Mis the Mach cone, while we dub the analogous figure swept out by $\mathit{\theta }\mathrm{=}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$the “anti-Mach cone”. As Fig. 3 illustrates, we denote as the “upstream” region that portion of the flow bounded by the anti-Mach cone, while the “downstream” region lies within the Mach cone. Finally, the “intermediate” region has ${\mathit{\theta }}_{\mathit{M}}\mathit{<}\mathit{\theta }\mathit{<}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ .

Fig. 3 Regions of the flow. The upstream and downstream regions lie within the anti-Mach and Mach cones, respectively. Between the two cones is the intermediate region. In the far-field flow, there are no physical barriers between the intermediate region and its neighbors, although mathematical divergences appear at the two cones.

A significant feature of the flow surrounding a gravitating mass is the presence of an accretion bowshock. Any fluid element that joins onto the mass penetrates to r ≲ r_acc , and encounters a shock at θ ~ θ_M . At much larger r, such bowshocks weaken and degenerate to acoustic pulses, before fading away entirely (Zel’dovich & Raizer 1968, Chap. 1). No such shock arises near ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$, even at relatively small distances. In summary, for r ≫ r_s ≳ r_acc , we do not expect any discontinuities at θ_M or ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ in steady-state flow, either in the fluid variables themselves or their derivatives. We need to keep this key point in mind as we encounter apparent discontinuities at both angles.

2.2. Mathematical formulation

As in Paper I, we describe the flow through two dependent variables, the mass density ρ(r,θ) and the stream function ψ(r,θ). Individual velocity components may be recovered from these through the relations $\begin{array}{ccc}{\mathit{u}}_{\mathit{r}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\rho }\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }}\hspace{0.17em}\hspace{0.17em}\frac{\mathit{\partial \psi }}{\mathit{\partial \theta }}\\ {\mathit{u}}_{\mathit{\theta }}& \mathrm{=}& \frac{-1}{\mathit{\rho }\hspace{0.17em}\mathit{r}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }}\hspace{0.17em}\hspace{0.17em}\frac{\mathit{\partial \psi }}{\mathit{\partial r}}\mathit{,}\end{array}$which automatically ensure mass continuity. In the extreme far-field limit, the stream function approaches that for the background, uniform flow: $\begin{array}{ccc}\underset{\mathit{r}\hspace{0.17em}\mathrm{\to }\hspace{0.17em}\mathrm{\infty }}{\lim }\mathit{\psi }& \mathrm{=}& \frac{{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}\mathit{V}\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}{\mathrm{2}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}\frac{\mathit{\beta }\hspace{0.17em}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}{\mathrm{2}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)}^{\mathrm{2}}\mathrm{·}\end{array}$Our second form of the stream function’s limit suggests how to expand ψ(r,θ) in a perturbation series valid for finite, but large, r: $\begin{array}{ccc}\mathit{\psi }& \mathrm{=}& {\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}[{\mathit{f}}_{\mathrm{2}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)}^{\mathrm{2}}\mathrm{+}{\mathit{f}}_{\mathrm{1}}\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}\mathrm{+}{\mathit{f}}_{-1}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)}^{-1}\mathrm{+}\mathit{...}\hspace{0.17em}]\mathit{.}\end{array}$(4)Here, f₂    ≡    β   sin²θ/2 , while f₁, f₀, f_-1, etc., are as yet unknown, non-dimensional functions of β and θ. We write a similar perturbation expansion for the density: $\begin{array}{ccc}\mathit{\rho }& \mathrm{=}& {\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}[\mathrm{1}\mathrm{+}{\mathit{g}}_{-1}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)}^{-1}\mathrm{+}{\mathit{g}}_{-2}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)}^{-2}\hspace{0.17em}\mathrm{+}{\mathit{g}}_{-3}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{s}}}\right)}^{-3}\mathrm{+}\mathit{...}\hspace{0.17em}]\mathit{.}\end{array}$(5)The quantities g_-1, g_-2, and g_-3 are again non-dimensional functions of β and θ, all of them unknown at this stage².

The dynamical equations will be simplified once we recast all variables into non-dimensional form. We let the fiducial radius, density, and speed be r_s, ρ₀, and c_s, respectively, while we normalize the stream function to ${\mathit{\rho }}_{\mathrm{0}}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}$. We will not change notation, but alert the reader whenever we revert to dimensional varables. Our fully non-dimensional perturbation series are $\begin{array}{ccc}\mathit{\psi }& \mathrm{=}& {\mathit{f}}_{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{f}}_{\mathrm{1}}\mathit{r}\mathrm{+}{\mathit{f}}_{\mathrm{0}}\mathrm{+}{\mathit{f}}_{-1}{\mathit{r}}^{-1}\mathrm{+}\mathit{...}\\ \mathit{\rho }& \mathrm{=}& \mathrm{1}\mathrm{+}{\mathit{g}}_{-1}{\mathit{r}}^{-1}\mathrm{+}{\mathit{g}}_{-2}{\mathit{r}}^{-2}\mathrm{+}{\mathit{g}}_{-3}{\mathit{r}}^{-3}\mathrm{+}\mathit{...}\end{array}$The task will be to solve for the various functions f_i(θ) and g_i(θ) appearing in these two series. (Henceforth, we will suppress the β-dependence for simplicity.) As explained in Sect. 2.3 of Paper I, the appropriate boundary conditions are ${\mathit{f}}_{\mathit{i}}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}$for i = 1, 0, −1, −2 , etc., and f_i(0) = 0for i = 1, −1, −2 , etc. These conditions ensure regularity of both u_r and u_θ, as given by Eqs. (2) and (3), respectively, on both the upstream and downstream axes. Since sinθ → 0 on both axes, the r- and θ-derivatives of ψ must also vanish. Using the expansion of ψ in Eq. (6), we derive the aforementioned conditions. There is no associated restriction on f₀(0), which is tied to the mass accretion rate onto the central object, as shown in Sect. 4 below.

3. First-order flow

3.1. Upstream and downstream regions

The r- and θ-components of Euler’s equation, in non-dimensional form, are

$\begin{array}{ccc}& & {\mathit{u}}_{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{r}}}{\mathit{\partial r}}\mathrm{+}\frac{{\mathit{u}}_{\mathit{\theta }}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{r}}}{\mathit{\partial \theta }}\mathrm{-}\frac{{\mathit{u}}_{\mathit{\theta }}^{\mathrm{2}}}{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\hspace{0.17em}\frac{\mathit{\partial \rho }}{\mathit{\partial r}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\\ & & {\mathit{u}}_{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{\theta }}}{\mathit{\partial r}}\mathrm{+}\frac{{\mathit{u}}_{\mathit{\theta }}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{\theta }}}{\mathit{\partial \theta }}\mathrm{+}\frac{{\mathit{u}}_{\mathit{r}}\hspace{0.17em}{\mathit{u}}_{\mathit{\theta }}}{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }\hspace{0.17em}\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial \rho }}{\mathit{\partial \theta }}\mathrm{·}\end{array}$Our procedure is first to express u_r and u_θ in terms of ψ and ρ, using Eqs. (2) and (3). We then expand ψ and ρ themselves in their respective perturbation series and equate the coeffcients of various powers of r. To avoid cumbersome division by the series for ρ, we first multiply Eqs. (8) and (9) through by ρ³.

Equating the coefficients of the highest power of r, which is r^-1, we find that they are identically equal. Equating coefficients of r^-2 yields the first-order equations. From Eq. (8), we find $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}& \mathrm{+}& {\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\\ & \mathrm{+}& \mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\cos }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{+}\mathrm{1}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(10)while Eq. (9) yields $\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{=}\mathrm{0.}$(11)These equations are identical in form to Eqs. (18) and (19), respectively, of Paper I, hereafter designated Eqs. (I.18) and (I.19). However, their solutions may differ. The factor (1 − β²   sin²   θ)in Eq. (11), which was positive in the subsonic case for any angle θ, can now be positive, negative, or zero.

Consider first the upstream region. For $\mathit{\theta }\mathit{>}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ , the term (1 − β²   sin²   θ)is indeed positive, and we may recast Eq. (11) as $\begin{array}{ccc}\frac{\mathrm{d}}{\mathrm{d}\mathit{\theta }}\left[{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{g}}_{-1}\right]\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$whose solution is ${\mathit{g}}_{-1}\mathrm{=}\frac{\mathit{C}}{{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}$(12)Here, C is independent of θ, but is possibly a function of β. The analogous result appeared in our subsonic analysis as Eq. (I.20). In that case, we ultimately found C to be unity.

Returning to the supersonic flow, we see that, for any non-zero C, the function g_-1 diverges at $\mathit{\theta }\mathrm{=}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ . This fact does not mean that the density itself diverges at that location; there is no reason for it to do so. Rather, it is our series expansion of ρ   (r,θ) that fails in this region. The issue here is mathematical, rather than physical. Although the function g_-1(θ) diverges, the associated term in the perturbation expansion, Eq. (7), is multiplied by r^-1. At any finite angular separation from ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$, the first-order correction to the background density can be made arbitrarily small by considering a sufficiently large r-value. This observation suggests that the divergences encountered here and elsewhere in our analysis will vanish if we change independent variables from r and θ to another set that mixes the two. In any case, we do not explore that possibility in the present paper.

Substitution of Eq. (12) for g_-1 into (10) yields the governing equation for f₁: ${\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{+}{\mathit{f}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\beta }}\mathrm{-}\frac{\mathit{C}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\mathrm{)}}{\mathit{\beta }{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}$(13)This equation, being identical to (I.21), has the same general solution: $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\beta }}\mathrm{-}\frac{\mathit{C}\hspace{0.17em}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{\beta }}\mathrm{+}\mathit{D}\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\mathrm{+}\mathit{E}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta ,}\end{array}$(14)where D and E are additional constants. The vanishing of f₁   (π) tells us that D = (1 − C)/β , while ${\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}\mathrm{=}\mathrm{0}$implies that E = 0 . Both relations also held in the subsonic case.

We next turn to the downstream region. Here again, the critical term (1 − β²   sin²   θ)is positive. The solutions for g_-1 and f₁ are thus identical to the upstream solutions in Eqs. (12) and (14), if we replace C, D, and E by new constants C′, D′, and E′. Application of the appropriate boundary conditions then tells us that E′ = 0and D′ = (C′ − 1)/β .

Referring again to Fig. 3, the supersonic flow we are now analyzing differs from the subsonic one by the presence of the intermediate region, lying between the Mach and anti-Mach cones. As β decreases to unity from above, this region narrows symmetrically about θ = π/2and ultimately vanishes, as does the physical distinction between the supersonic and subsonic flows. Continuity strongly suggests that, while (C,D)and (C′,D′)could, in principle, be β-dependent, they are actually not. Instead, they have the truly constant values (1,0)found in the subsonic case.

In summary, the first-order flow in both the upstream and downstream regions is described by

$\begin{array}{ccc}& & {\mathit{f}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{\beta }}\\ & & {\mathit{g}}_{-1}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}\end{array}$After substituting these coefficients into the series expansions for ψ and ρ, one may calculate u_r and u_θ to linear order from Eqs. (2) and (3), respectively. At any fixed, r-value, however, the series representation for ρ fails sufficiently close to either Mach cone, and the velocity components cannot be obtained in this manner.

Comparing our results thus far with the previous literature, Dokuchaev (1964); Ruderman & Spiegel (1971); and Ostriker (1999) all obtained Eq. (12) for the downstream, first-order density distribution, but with C = 2 . They also found that g_-1 vanishes for θ > θ_M . In the subsonic analysis of Ostriker (1999), g_-1 is everywhere finite. Thus, the density distributions in her supersonic and subsonic flows do not approach one another in the β = 1limit. Ruderman & Spiegel (1971) further argued that the apparent divergence in the downstream density perturbation at θ = θ_Msignifies the presence of a bowshock. While a bowshock certainly arises relatively close to the gravitating mass, it does not persist into the far field, the only regime where an analysis based on small perturbations of the background gas is justified.

3.2. Intermediate region

For ${\mathit{\theta }}_{\mathit{M}}\mathit{<}\mathit{\theta }\mathit{<}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ , the term (1 − β²   sin²   θ)in Eq. (11) is negative. We therefore rewrite this equation as $\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{=}\mathrm{0.}$(17)This is equivalent to $\begin{array}{ccc}\frac{\mathrm{d}}{\mathrm{d}\mathit{\theta }}\left[{\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{g}}_{-1}\right]\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$which in turn implies that ${\mathit{g}}_{-1}\mathrm{=}\frac{\mathit{C}}{{\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}$(18)Here we are reverting to our original notation for the integration constants, as the previous ones have all been evaluated.

Substitution of the new expression for g_-1 into Eq. (10) yields ${\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{+}{\mathit{f}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\beta }}\mathrm{-}\frac{\mathit{C}\hspace{0.17em}\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}{\mathit{\beta }{\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}$(19)We may solve this equation using the method of variation of parameters. Adding the two homogeneous solutions yields the general result $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\beta }}\mathrm{+}\frac{\mathit{C}\hspace{0.17em}{\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{\beta }}\mathrm{+}\mathit{D}\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\mathrm{+}\mathit{E}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\mathit{.}\end{array}$(20)As in the upstream and downstream regions, the series representation of ψ is well-behaved at either Mach cone, at least to linear order. Approaching either cone from the outside, Eq. (15) tells us that f₁ → 1/β . Equation (20) above implies that f₁ has the additional term D   cos   θ + E   sin   θwhen we approach the cones from the intermediate region.

To proceed, we consider the physical interpretation of the stream function. Equation (2) tells us that ψ(r,θ) is the θ-integral of the mass flux ρ   u_r over a surface of radius r. More generally, the stream function is the net rate of mass transport into any surface of revolution extending from θ = 0to the angle of interest. A discontinuity in ψ at an angle θ thus represents a thin sheet of mass being injected or ejected along the corresponding cone. If we are to reject such a solution as unphysical, then we must demand continuity of ψ. That is, f₁ must again approach 1/β at the two Mach cones, and D = E = 0in Eq. (20).

Our mathematical description of the first-order flow in the intermediate region is now

$\begin{array}{ccc}& & {\mathit{f}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{C}\hspace{0.17em}{\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{\beta }}\\ & & {\mathit{g}}_{-1}\mathrm{=}\frac{\mathit{C}}{{\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}\end{array}$Since the intermediate region does not reach either the upstream or downstream axis, we cannot appeal to our usual boundary conditions in order to evaluate C. We defer this issue until our analysis of the second-order flow in Sect. 4, when we will show that the requirement of mass continuity again settles the matter. For any C-value, the divergence of g_-1 at the Mach and anti-Mach cones indicates a breakdown of the series representation for the density.

3.3. Vorticity

A key, simplifying property of the flow is that it is irrotational. To reprise the argument from Paper I, we first write Euler’s equation in the steady state as ${u}\mathrm{\times }{\omega }\mathrm{=}{\nabla }\mathit{B,}$(23)where ω    ≡ ∇ × uis the vorticity and the Bernoulli function on the right-hand side is $\mathit{B}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}{\mathit{u}}^{\mathrm{2}}\mathrm{+}\ln \hspace{0.17em}\mathit{\rho }\mathrm{-}\frac{\mathrm{1}}{\mathit{r}}\mathrm{·}$(24)Dotting both sides of Eq. (23) with u, we find that $\begin{array}{ccc}\mathrm{(}{u}\hspace{0.17em}{·}\hspace{0.17em}{{\nabla }}^{\mathrm{)}}\hspace{0.17em}\mathit{B}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$which is the familiar statement that B is constant along streamlines. Moreover, B approaches β²/2 at large r, both upstream and downstream. (Recall that the non-dimensional density ρ becomes unity in this limit.) In these regions, therefore, B has the same value on every streamline, and ∇B = 0 . Since ω is orthogonal to u in our poloidal flow, Eq. (23) imples that ω = 0 , i.e., the flow is irrotational outside the Mach cones.

If our derived intermediate solution is correct, then irrotationality must hold there as well, since there is no physical barrier between the intermediate flow and those upstream and downstream, at least in the laminar, far field. The only non-zero component of the vorticity is ω_φ, so the condition of irrotationality becomes $\frac{\mathit{\partial }{\mathit{u}}_{\mathit{r}}}{\mathit{\partial \theta }}\mathrm{=}\frac{\mathit{\partial }\mathrm{\left(}\mathit{r}{\mathit{u}}_{\mathit{\theta }}\mathrm{\right)}}{\mathit{\partial r}}\mathrm{·}$(25)To test the validity of this relation in the intermediate region, we use Eqs. (2) and (3) for u_r and u_θ, respectively, and then substitute in the series expansions for ψ and ρ, Eqs. (6) and (7).

Equating the coefficient of the highest power of r, which is r⁰, is equivalent to testing for irrotationality in the uniform, background flow. Specifically, we require ${\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime \prime }}\mathrm{-}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime }}\hspace{0.17em}\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\mathrm{=}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{f}}_{\mathrm{2}}\mathit{.}$(26)Using f₂ = β   sin²   θ/2 , we verify that the above equation does hold. This result is to be expected, as a uniform flow is manifestly irrotational.

We next equate the coefficients of r^-1, effectively testing irrotationality in the first-order flow. The required condition is now ${\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{-}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime \prime }}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{-}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime }}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\mathrm{=}\left({\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{-}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime }}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\right)\cot \hspace{0.17em}\mathit{\theta ,}$(27)where f₁ and g_-1 are given by Eqs. (21) and (22), respectively. Using these functional forms, we find that this last equation is satisfied for any value of C. Thus, the first-order flow in the intermediate region is indeed irrotational.

4. Second-order flow

4.1. Dynamical equations

Having established the first-order flow, at least up to the constant C, we consider the next higher approximation. We return to the perturbative expansion of Euler’s equation, as described at the beginning of Sect. 3.1. By equating coefficients of r^-3, we derive the second-order equations, which govern the variables f₀ and g_-2. Here the source terms involve f₁ and g_-1. Prior to substituting in the explicit solutions for f₁ and g_-1, the equations are identical to those in the subsonic problem; we display them again for convenient reference.

From the r-component of Euler’s equation, we derive Eq. (I.27), which is $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime \prime }}& \mathrm{-}& \mathit{\beta }\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\\ & & \mathrm{+}\hspace{0.17em}\mathrm{(}\mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\cos }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{2}}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{=}{\mathrm{𝒜}}_{\mathrm{1}}\mathrm{+}{\mathrm{𝒜}}_{\mathrm{2}}\mathrm{+}{\mathrm{𝒜}}_{\mathrm{3}}\mathit{.}\end{array}$(28)The three right-hand terms are $\begin{array}{ccc}{\mathrm{𝒜}}_{\mathrm{1}}& \mathrm{\equiv }& \frac{{\mathit{f}}_{\mathrm{1}}^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}}\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }}{{\sin }^{\mathrm{3}}\hspace{0.17em}\mathit{\theta }}\mathrm{+}\frac{{\mathrm{(}{{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}}^{\mathrm{)}}}^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\mathrm{+}\frac{{\mathit{f}}_{\mathrm{1}}\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime \prime }}}{{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\\ {\mathrm{𝒜}}_{\mathrm{2}}& \mathrm{\equiv }& \mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{-}\mathrm{2}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\\ & & \mathrm{+}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\\ {\mathrm{𝒜}}_{\mathrm{3}}& \mathrm{\equiv }& \mathrm{2}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{2}}\mathrm{-}\mathrm{3}\hspace{0.17em}{\mathit{g}}_{-1}\mathit{.}\end{array}$The θ-component of Euler’s equation yields Eq. (I.31): $\begin{array}{ccc}& \mathrm{-}& \mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}\mathrm{𝒟}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{=}{\mathrm{ℬ}}_{\mathrm{1}}\mathrm{+}{\mathrm{ℬ}}_{\mathrm{2}}\mathrm{+}{\mathrm{ℬ}}_{\mathrm{3}}\mathit{.}\end{array}$(29)Here, $\mathrm{𝒟}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta ,}$(30)and $\begin{array}{ccc}{\mathrm{ℬ}}_{\mathrm{1}}& \mathrm{\equiv }& {\mathit{f}}_{\mathrm{1}}^{\mathrm{2}}\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}}\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}}{{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\\ {\mathrm{ℬ}}_{\mathrm{2}}& \mathrm{\equiv }& \mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\mathrm{+}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\\ {\mathrm{ℬ}}_{\mathrm{3}}& \mathrm{\equiv }& \mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{g}}_{-1}\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\mathit{.}\end{array}$In the upstream and downstream regions, there are no unknown constants in the source terms, i.e., both f₁ and g_-1 are identical to their subsonic counterparts. Hence, the explicit form of the second-order equations is also the same. Referring to Eqs. (I.35) and (I.36), we have $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime \prime }}& \mathrm{-}& \mathit{\beta }\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\\ & \mathrm{+}& \mathrm{(}\mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\cos }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{2}}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{=}\frac{\mathrm{1}}{\mathrm{𝒟}}\mathrm{-}\frac{\mathrm{3}}{\sqrt{\mathrm{𝒟}}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{1}\mathrm{+}\sqrt{\mathrm{𝒟}}}\mathit{,}\end{array}$(31)and $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}& \mathrm{+}& \mathrm{𝒟}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-2}\\ & \mathrm{=}& {\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }[\mathrm{-}\frac{\mathrm{2}}{{\mathrm{𝒟}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}}{{\mathrm{𝒟}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\hspace{0.17em}\mathrm{-}\frac{\mathrm{1}}{\mathrm{𝒟}}\mathrm{+}\frac{\mathrm{1}}{{\mathrm{(}\mathrm{1}\mathrm{+}{\sqrt{\mathrm{𝒟}}}^{\mathrm{)}}}^{\mathrm{2}}}]\mathrm{·}\end{array}$(32)In the intermediate region, however, f₁ and g_-1 are given by the new expressions in Eqs. (21) and (22), respectively. After lengthy manipulation, we find the new source terms, and hence the explicit dynamical equations in this region. These equations, which still contain the unknown constant C, are

$\begin{array}{ccc}& & \mathrm{-}\hspace{0.17em}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime \prime }}\mathrm{-}\mathit{\beta }\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\mathrm{+}\mathrm{(}\mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\cos }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}{\mathrm{2}}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{=}\\ & & \frac{\mathrm{2}\hspace{0.17em}\mathit{C}\hspace{0.17em}{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\mathrm{-}\frac{\mathrm{3}\hspace{0.17em}\mathit{C}}{{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{+}\frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{ℰ}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}{\mathit{C}}^{\mathrm{2}}}{{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\mathit{,}\end{array}$(33)and $\begin{array}{ccc}& & \mathrm{-}\hspace{0.17em}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{-}\mathrm{ℰ}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{=}\hspace{0.17em}\\ & & \frac{\mathrm{2}\hspace{0.17em}{\mathit{C}}^{\mathrm{2}}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }}{{\mathrm{ℰ}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{2}\hspace{0.17em}\mathit{C}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }}{{\mathrm{ℰ}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\\ & & \mathrm{-}\hspace{0.17em}\hspace{0.17em}\frac{{\mathit{C}}^{\mathrm{2}}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }}{\mathrm{ℰ}}\mathrm{+}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{C}}^{\mathrm{2}}\mathrm{)}\cot \hspace{0.17em}\mathit{\theta }}{{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\\ & & \mathrm{+}\hspace{0.17em}\hspace{0.17em}{\mathit{C}}^{\mathrm{2}}\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\mathrm{+}\hspace{0.17em}\frac{\mathrm{2}\hspace{0.17em}\mathit{C}\hspace{0.17em}{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }}{{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\mathrm{·}\end{array}$(34)Here, we have defined $\mathrm{ℰ}\hspace{0.17em}\mathrm{\equiv }\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\mathrm{-}\mathrm{1}\mathit{,}$(35)which is positive throughout this region.

4.2. Near-cone divergence

Our plan is to integrate numerically two, coupled second-order equations in order to determine the variables ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 throughout the flow. With ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ in hand, another integration will yield f₀ itself. Before we embark on this program, let us consider more carefully the behavior of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2. In the upstream and downstream regions, the term vanishes as one approaches either Mach cone. Thus, the source terms in Eqs. (31) and (32) diverge in that limit. Similarly, the source terms in Eqs. (33) and (34) diverge because of the vanishing of ℰ. It is possible, then, that both ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 also diverge at the cones.

This is indeed the case. In this section, we first establish the divergences of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 in the downstream region only. That is, we assume that θ approaches θ_M from below. Later in the section, we outline the derivation for the remaining regions of the flow. We then use the full result to shed light on the unknown parameter C.

Since θ    <    θ_M, we define the positive angle α    ≡    θ_M − θ, and assume that ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 take the following asymptotic forms as α diminishes:

$\begin{array}{ccc}& & {\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{=}\mathit{\zeta }\hspace{0.17em}\hspace{0.17em}{\mathit{\alpha }}^{\mathrm{-}\mathit{m}}\\ & & {\mathit{g}}_{-2}\mathrm{=}\mathit{\gamma }\hspace{0.17em}\hspace{0.17em}{\mathit{\alpha }}^{\mathrm{-}\mathit{n}}\mathit{,}\end{array}$where ζ, γ, m, and n are all functions of β only. If both ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 truly diverge, then m and n are positive. Within Eqs. (31) and (32), we replace sin   θ and cos   θ by their limiting values, 1/β and $\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\mathit{/}\mathit{\beta }$, respectively. We further note that $\mathrm{𝒟}\mathrm{\to }\mathrm{2}\mathit{\alpha }\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}$. After differentiating the asymptotic forms of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 with respect to θ, we derive simplified, limiting forms of Eqs. (31) and (32). Expressed using α as the independent variable, these are $\begin{array}{ccc}& & \mathrm{-}\hspace{0.17em}\hspace{0.17em}\mathit{m}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}\mathit{\zeta }\hspace{0.17em}{\mathit{\alpha }}^{\mathrm{-}\mathit{m}\mathrm{-}\mathrm{1}}\mathrm{+}\mathit{n}\hspace{0.17em}\mathit{\gamma }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}{\mathit{\alpha }}^{\mathrm{-}\mathit{n}\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{-}\frac{{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{4}\hspace{0.17em}\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\hspace{0.17em}\hspace{0.17em}{\mathit{\alpha }}^{-1}\\ & & \mathit{\beta }\hspace{0.17em}\mathit{\zeta }\hspace{0.17em}{\mathit{\alpha }}^{\mathrm{-}\mathit{m}}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{n}\mathrm{)}\hspace{0.17em}\mathit{\gamma }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\hspace{0.17em}\hspace{0.17em}{\mathit{\alpha }}^{\mathrm{-}\mathit{n}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}}\hspace{0.17em}\hspace{0.17em}{\mathit{\alpha }}^{-2}\mathit{.}\end{array}$In deriving these equations, we have retained only the most rapidly diverging terms on the right-hand sides of Eqs. (31) and (32). We have also dropped terms proportional to α^−m and α⁻ⁿ in the left-hand side of Eq. (31), since both terms are dominated, in the small-α limit, by the two we have kept.

We next consider the β-dependence of m and n. The right-hand side of Eq. (39) is proportional to α^-2. For this equation to balance as α diminishes, one left-hand term could also diverge as α^-2, and the other more slowly. However, if this were the case, i.e., if either m = 2 > nor n = 2 > m , then Eq. (38) would be unbalanced. It is also posssible that the two left-hand terms in (39) diverge more rapidly than α^-2, balancing each other. We would then have m = n > 2 . Finally, all three terms in Eq. (39) could diverge at the same rate: m = n = 2 .

Assuming provisionally that m = n > 2 , then we may ignore the right-hand terms of both Eqs. (38) and (39) in the asymptotic limit. After dividing out all terms containing α, the two equations reduce to

$\begin{array}{ccc}& & \mathrm{-}\mathit{\beta }\hspace{0.17em}\mathit{\zeta }\mathrm{+}\mathit{\gamma }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{0}\\ & & \mathit{\beta }\hspace{0.17em}\mathit{\zeta }\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{n}\mathrm{)}\hspace{0.17em}\mathit{\gamma }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{0.}\end{array}$Adding these two yields $\begin{array}{ccc}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathit{n}\mathrm{)}\hspace{0.17em}\mathit{\gamma }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$which implies that n = 3/2 . Since we assumed n > 2at the outset, our original hypothesis, m = n > 2 , was incorrect.

We have established that m = n = 2 . Equations (38) and (39) now reduce asymptotically to $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}\mathit{\zeta }\mathrm{+}\mathit{\gamma }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}& \mathrm{=}& \mathrm{0}\\ \mathit{\beta }\hspace{0.17em}\mathit{\zeta }\mathrm{-}\mathrm{2}\hspace{0.17em}\mathit{\gamma }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}}\mathit{,}\end{array}$which have the unique solution $\begin{array}{ccc}\mathit{\zeta }& \mathrm{=}& \frac{-1}{\mathrm{2}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}}\\ \mathit{\gamma }& \mathrm{=}& \frac{-1}{\mathrm{2}\hspace{0.17em}\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{·}\end{array}$Consider next the region just upstream from the Mach cone. Here, we may assume the same asymptotic forms for ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 as in Eqs. (36) and (37), provided the independent variable α remains positive: α    ≡    θ − θ_M. Inserting these functional forms into the second-order Eqs. (33) and (34), we note that dα/dθ changes sign from −1 to +1, and that ℰ now approaches $\mathrm{2}\hspace{0.17em}\mathit{\alpha }\hspace{0.17em}\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}$ near the cone. We then derive equations analogous to (38) and (39), which may again be solved in the near-cone limit. After applying similar reasoning to the two regions surrounding the anti-Mach cone, we arrive at the following asymptotic forms for ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 near θ_M and ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$: $\begin{array}{ccc}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}& \mathrm{=}& \frac{-1}{\mathrm{2}\mathit{\beta }\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\hspace{0.17em}{\mathrm{(}{\mathit{\theta }}_{\mathit{M}}\mathrm{-}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\lesssim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}\\ & \mathrm{=}& \frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{2}\mathit{\beta }\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\hspace{0.17em}{\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathit{M}}\mathrm{)}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\gtrsim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}\\ & \mathrm{=}& \frac{\mathrm{-}{\mathit{C}}^{\mathrm{2}}}{\mathrm{2}\mathit{\beta }\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\hspace{0.17em}{\mathrm{(}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\mathrm{-}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\lesssim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}\mathit{\beta }\sqrt{{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\hspace{0.17em}{\mathrm{(}\mathit{\theta }\mathrm{-}{{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}}^{\mathrm{)}}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\gtrsim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\mathit{,}\end{array}$(40)and $\begin{array}{ccc}{\mathit{g}}_{-2}& \mathrm{=}& \frac{-1}{\mathrm{2}\hspace{0.17em}\left({\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\right)\hspace{0.17em}{\mathrm{(}{\mathit{\theta }}_{\mathit{M}}\mathrm{-}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\lesssim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}\\ & \mathrm{=}& \frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\left({\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\right)\hspace{0.17em}{\mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathit{M}}\mathrm{)}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\gtrsim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}\\ & \mathrm{=}& \frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\left({\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\right)\hspace{0.17em}{\mathrm{(}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\mathrm{-}{\mathit{\theta }}^{\mathrm{)}}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\lesssim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\\ & \mathrm{=}& \frac{-1}{\mathrm{2}\hspace{0.17em}\left({\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\right)\hspace{0.17em}{\mathrm{(}\mathit{\theta }\mathrm{-}{{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}}^{\mathrm{)}}}^{\mathrm{2}}} \mathit{\theta }\hspace{0.17em}\gtrsim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}\mathit{.}\end{array}$(41)Obtaining f₀ requires that we integrate ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ over θ. Starting at the upstream axis, with f₀(π) = 0 , the integrated f₀(θ) diverges to positive infinity as θ approaches ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$. However, the stream function has a direct physical meaning as the mass transfer rate into a surface of revolution. Since this rate is finite for any θ, the function ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$ must be integrable across the anti-Mach cone. That is, the upward divergence of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ for $\mathit{\theta }\hspace{0.17em}\gtrsim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$must be cancelled by a matching downward divergence for $\mathit{\theta }\hspace{0.17em}\lesssim \hspace{0.17em}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ . A similar antisymmetry of the divergences must occur at the Mach cone, θ = θ_M(see Fig. 4). Inspection of Eq. (40) shows that this requirement forces C² to be unity. The parameter C itself is either +1 or −1, independent of β. We will demonstrate presently that the positive solution is the physically relevant one.

4.3. Vorticity

Fig. 4 Asymmetric divergences at the Mach angle, for β = 1.1 . The solid and dashed curves trace ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g-2, respectively, close to the Mach angle θM, marked here by the vertical, dotted line. The functions diverge antisymmetrically once we choose C2 = 1 . Both curves are taken from the numerical integration outlined in Sect. 4.5. The shaded region is that which we remove before enforcing continuity of f0 and g-2 across the cones.

The second-order flow must be irrotational, as we argued in Sect. 3.3. Starting with Eq. (25), we again replace the velocity components by ψ and ρ, and develop the latter in our perturbation series. Equating coefficients of r^-2 tests irrotationality in the second-order flow. Specifically, we require

$\begin{array}{ccc}& & \mathrm{(}\mathrm{-}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{2}}\mathrm{+}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{g}}_{-2}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{-}{{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}}^{\mathrm{)}}\cot \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}\\ & & \mathrm{+}\hspace{0.17em}\hspace{0.17em}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime \prime }}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{2}}\mathrm{+}\mathrm{2}\hspace{0.17em}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{g}}_{-1}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\mathrm{-}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime \prime }}\hspace{0.17em}{\mathit{g}}_{-2}\\ & & \mathrm{-}\hspace{0.17em}\hspace{0.17em}{\mathit{f}}_{\mathrm{2}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\mathrm{-}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\hspace{0.17em}{\mathit{g}}_{-1}\mathrm{-}{\mathit{f}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{g}}_{-1}^{\mathrm{\prime }}\mathrm{+}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime \prime }}\hspace{0.17em}\\ & & \end{array}$(42)Within the upstream and downstream regions, Eqs. (15) and (16) give us f₁ and g_-1, respectively. Substituting these functions into Eq. (42) then yields the condition of irrotationality: $\begin{array}{ccc}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime \prime }}\mathrm{-}\mathit{\beta }\hspace{0.17em}\cot \mathit{\theta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}& \mathrm{-}& {\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\hspace{0.17em}\\ & \mathrm{+}& \hspace{0.17em}\mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{=}\frac{\mathrm{1}}{\mathrm{𝒟}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{\mathrm{𝒟}}}\mathrm{·}\end{array}$(43)Finally, we may combine this last relation with Eq. (31), the r-component of Euler’s equation, to obtain $\begin{array}{ccc}\mathrm{-}\mathrm{2}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}{\mathit{g}}_{-2}\mathrm{=}\frac{\mathrm{2}}{\mathrm{𝒟}}\mathrm{-}\frac{\mathrm{4}}{\sqrt{\mathrm{𝒟}}}\mathrm{+}\frac{\mathrm{2}}{\mathrm{1}\mathrm{+}\sqrt{\mathrm{𝒟}}}\mathrm{·}& & \end{array}$(44)This is identical to Eq. (I.64), and will prove useful when we evaluate the force.

Turning to the intermediate region, we need to use Eqs. (21) and (22) for f₁, g_-1, and their derivatives. In both equations, the unknown parameter C appears. Substitution into Eq. (42) above yields the condition for irrotationality: $\begin{array}{ccc}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime \prime }}& \mathrm{-}& \mathit{\beta }\hspace{0.17em}\cot \mathit{\theta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}\hspace{0.17em}\\ & \mathrm{+}& \mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{=}\frac{{\mathit{C}}^{\mathrm{2}}}{\mathrm{ℰ}}\mathrm{-}\frac{\mathit{C}}{\sqrt{\mathrm{ℰ}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{ℰ}}\mathrm{-}\frac{\mathit{C}}{\sqrt{\mathrm{ℰ}}}\mathrm{·}\end{array}$(45)In the last form of this equation, we have used the fact that C² = 1 . Combination with the Euler Eq. (31) gives $\begin{array}{ccc}\mathrm{-}\mathrm{2}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}\cot \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}& \mathrm{+}& \mathrm{2}\hspace{0.17em}\mathrm{(}{\mathit{\beta }}^{\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}{\mathit{g}}_{-2}\\ & \mathrm{=}& \hspace{0.17em}\frac{\mathrm{2}\hspace{0.17em}{\mathit{C}}^{\mathrm{2}}}{\mathrm{ℰ}}\mathrm{-}\frac{\mathrm{4}\hspace{0.17em}\mathit{C}}{{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}{\mathit{C}}^{\mathrm{2}}}{{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\mathrm{+}\frac{\mathrm{2}\hspace{0.17em}\mathit{C}\hspace{0.17em}{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\\ & \mathrm{=}& \frac{\mathrm{2}}{\mathrm{ℰ}}\mathrm{-}\frac{\mathrm{4}\hspace{0.17em}\mathit{C}}{{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}\hspace{0.17em}\mathit{C}\hspace{0.17em}{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}\mathit{\theta }}\mathit{,}\end{array}$(46)which will again aid in the force evaluation.

Fig. 5 Enforcing irrotationality. In our steady-state, isothermal flow, the Bernoulli function B is constant along each streamline. In principle, B could vary from one streamline to the next. However, if we also make B constant along any cone, then the function is a universal constant and the flow is irrotational.

We emphasize a key difference between Eqs. (44) and (46). Because the upstream and downstream regions join smoothly onto the uniform background, the flow in both is guaranteed to be irrotational to all orders. In our numerical determination of the flow, to be described in Sect. 4.5 below, we verified that Eq. (44) indeed holds numerically both upstream and downstream. In contrast, we need to impose Eq. (46), or an equivalent condition, to obtain a physically acceptable solution in the intermediate region.

Enforcing irrotationality in the intermediate region requires only that we impose this condition along any interior cone, i.e., at a θ-value such that ${\mathit{\theta }}_{\mathit{M}}\mathit{<}\mathit{\theta }\mathit{<}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ . To see why, return to Euler’s equation, in the form given by Eq. (23). Projecting this vector equation along the r-direction gives $\begin{array}{ccc}{\mathit{u}}_{\mathit{\theta }}{\mathit{\omega }}_{\mathit{\phi }}\mathrm{=}\frac{\mathit{\partial B}}{\mathit{\partial r}}\mathrm{·}& & \end{array}$This relation holds at fixed θ and φ. Since the flow is axisymmetric, the equation also holds at all φ-values, i.e., along a cone. If ω_φ vanishes along such a cone, then B is constant on that surface: ∂B/∂r = 0 . But we already know that B does not vary along any streamline. As Fig. 5 illustrates schematically, the additional constraint tells us that the B-values characterizing each streamline are identical, i.e., B is a true spatial constant. Thus, ∇B = 0and, from Eq. (23), the flow is irrotational.

In practice, we apply ω_φ = 0at θ = π/2 , the equatorial plane of the flow. The irrotationality condition, expressed as Eq. (46), implies that there is a unique (but β-dependent) value of g_-2 at that angle: $\begin{array}{ccc}{\mathit{g}}_{-2}\hspace{0.17em}\mathrm{(}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{)}& \mathrm{=}& \frac{\mathrm{2}}{({\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}{)}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{2}\hspace{0.17em}\mathit{C}}{({\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}{)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{+}\hspace{0.17em}\frac{\mathit{C}}{\mathit{\beta }\hspace{0.17em}({\mathit{\beta }}^{\mathrm{2}}\mathrm{-}\mathrm{1}{)}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}\end{array}$(47)When we find the flow numerically in the intermediate region, we will impose Eq. (47) as an initial condition. The resulting flow is then irrotational. That is, all solutions of Euler’s Eqs. (33) and (34) also obey Eq. (46).

4.4. Mass accretion

It is only in the second-order flow that mass accretion onto the central object appears. Moreover, the function f₀(θ), evaluated on the downstream axis, gives the actual accretion rate. Higher-order approximations to the flow contribute no additional information in this regard. We now review the argument, first advanced in Paper I, and apply it to the present, supersonic, case.

Consider the net transfer rate of mass into a sphere of radius r. In a steady-state flow, this rate is the same through any closed surface; we simply use a sphere for convenience. The dimensional result is $\begin{array}{ccc}\mathit{Ṁ}& \mathrm{=}& \mathrm{-}\mathrm{2}\hspace{0.17em}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\mathit{\rho }\hspace{0.17em}{\mathit{u}}_{\mathit{r}}\hspace{0.17em}{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\hspace{0.17em}\mathrm{d}\mathit{\theta }\mathrm{=}\mathrm{2}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\psi }\mathrm{\left(}\mathit{r,}\mathrm{0}\mathrm{\right)}\mathit{.}\end{array}$(48)Here we have substituted Eq. (2) for u_r and utilized the normalization for the stream function, ψ(r,π) = 0 ³. We have also implicitly assumed that ρ and u_r are smooth functions, despite the divergences arising in their perturbation expansions.

If we identify $\mathrm{2}\mathit{\pi }{\mathit{\rho }}_{\mathrm{0}}{\mathit{c}}_{\mathrm{s}}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}$as our fiducial mass accretion rate, then the non-dimensional counterpart of the last equation becomes $\begin{array}{ccc}\mathit{Ṁ}\mathrm{=}\mathit{\psi }\mathrm{\left(}\mathit{r,}\mathrm{0}\mathrm{\right)}\mathit{.}& & \end{array}$Substitution of the series expansion for ψ(r,0) in Eq. (6) and application of the boundary condition that f_i(0) = 0for i = 1, −1, −2, etc. leads to $\mathit{Ṁ}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(49)It is noteworthy that Ṁ depends only on a coefficient from the second-order expansion; this fact calls for a more physical explanation. One reason is that the streamlines of the first-order flow are symmetric, so that accretion does not occur to this approximation. Secondly, the mass flux ρ   u_r, when expanded using approximations higher than second-order, generates terms that fall off faster than r^-2. After integrating these terms over a large bounding sphere and taking the large-r limit, they do not contribute to Ṁ. In any event, Eq. (49) underscores the fact that the function ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$must be integrable across both Mach cones, despite its divergent behavior at the cones themselves.

We may now use our second-order equations to obtain a useful relation between f₀(0), and hence the mass accretion rate, and the flow density. Consider first the upstream and downstream regions. The left-hand side of Eq. (32) is a perfect derivative: $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{+}\mathrm{𝒟}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}& \mathrm{-}& \mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{g}}_{-2}\\ & \mathrm{=}& \frac{\mathrm{d}}{\mathrm{d}\mathit{\theta }}\mathrm{[}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{]}\mathit{.}\end{array}$In the right-hand side of the same equation, we note that is an even function, in the sense that . Since cos   θ is an odd function, cos   θ =  −cos   (π − θ) , the entire right-hand side of Eq. (32) is odd.

Within the intermediate region, the left-hand side of Eq. (34) is $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{-}\mathrm{ℰ}\hspace{0.17em}{\mathit{g}}_{-2}^{\mathrm{\prime }}& \mathrm{-}& \mathrm{2}\hspace{0.17em}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}\sin \hspace{0.17em}\mathit{\theta }\hspace{0.17em}\cos \hspace{0.17em}\mathit{\theta }\hspace{0.17em}{\mathit{g}}_{-2}\hspace{0.17em}\\ & \mathrm{=}& \frac{\mathrm{d}}{\mathrm{d}\mathit{\theta }}\mathrm{[}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}\hspace{0.17em}{\sin }^{\mathrm{2}}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{)}}\hspace{0.17em}{\mathit{g}}_{-2}\mathrm{]}\mathit{.}\end{array}$Again, the right-hand side of Eq. (34) is odd. Now the integral of an odd function from the downstream to the upstream axes vanishes. Thus, if we integrate Eq. (32) from θ = 0to θ_M, Eq. (34) from θ_M to ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$, and Eq. (32) from ${\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ to π, we find that $\begin{array}{ccc}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{+}{\mathit{g}}_{-2}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{\beta }\hspace{0.17em}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}\mathrm{+}{\mathit{g}}_{-2}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}\mathit{.}& & \end{array}$Using f₀(π) = 0 , we have the desired result: ${\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{{\mathit{g}}_{-2}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{-}{\mathit{g}}_{-2}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}}{\mathit{\beta }}\mathrm{·}$(50)The same relation appeared in the subsonic study as Eq. (I.46).

Fig. 6 First-order streamlines. Shown are contours of constant f2   r2 + f1   r, for the case β = 1.1 . The dotted, diagonal lines trace the Mach and anti-Mach cones, while the central, dotted circle is r/rs = 1 . Within the intermediate region, the dashed streamlines were constructed assuming C =  −1 . The solid ones correspond to C =  +1and are more realistic.

4.5. Numerical solution

4.5.1. Evaluation of C

Our description of the flow is clearly incomplete until we identify the still unknown parameter C, and not just its absolute value. We again employ physical reasoning. The value of C affects the shape of the streamlines, even in the first-order flow. Only one shape is reasonable dynamically. Figure 6 shows streamlines of the first-order flow for the case β = 1.1 . That is, we plot contours of constant f₂   r² + f₁   r . In the upstream and donwstream regions, f₁ is taken from Eq. (15). For the intermediate region, we use Eq. (21), with both C =  +1(solid curves) and C =  −1(dashed curves). For either choice of the parameter, we see that all first-order streamlines are symmetric about the equatorial plane, as we found in the subsonic flow (see, e.g., Fig. I.2). The kinks at both Mach cones result from the crudeness of this first-order approximation, and would disappear in successively more accurate treatments.

The solid curves in Fig. 6 bow inward, toward the central mass, both upstream of the anti-Mach cone and inside it. Such inward turning results naturally from the gravitational pull of the mass. Eventually, the concurrent rise in density and pressure pushes the gas outward. If we were instead to choose C =  −1 , the streamlines in the intermediate region would bow outward, which is not expected.

We may also consider the matter more quantitatively by examining the velocity component u_R. Here R is the cylindrical radius: R    ≡    r   sin   θ . In the far-field limit, u_R tends to zero, since the background flow is in the z-direction. Closer to the equatorial plane, but still upstream, u_R should be slightly negative. Downstream, this velocity component should reverse sign as each streamline rejoins the background flow. Starting with the expreessions for f₁ and g_-1 in the intermediate region, we first calculate u_r and u_θ from Eqs. (2) and (3), and then obtain u_R    ≡    u_r   sin   θ +    u_θ   cos   θ . We find that $\begin{array}{ccc}{\mathit{u}}_{\mathit{R}}\mathrm{=}\left[\frac{\mathit{C}\mathrm{-}{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{r}\hspace{0.17em}\mathit{\beta }\hspace{0.17em}{\mathrm{ℰ}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\right]\cot \hspace{0.17em}\mathit{\theta }\mathit{.}& & \end{array}$Upstream from the mass ${\mathrm{(}}^{}\mathit{\theta }\hspace{0.17em}\mathit{>}\hspace{0.17em}\mathit{\pi }\mathit{/}{\mathrm{2}}^{\mathrm{)}}$ , we have cot   θ    >    0 . Thus, if C =  −1 , then u_R    >    0in this direction. Conversely, u_R    <    0downstream from the mass. This behavior is contrary to the physically reasonable one, so we conclude that C =  +1 .

4.5.2. Integration scheme

Determining the second-order flow requires that we numerically integrate Eqs. (31) and (32) in the upstream and downstream regions, and Eqs. (33) and (34), with C =  +1 , in the intermediate region. In all cases, we need to specify appropriate boundary values of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2. We then perform an additional integration to obtain f₀ and thereby the streamlines.

One of our boundary conditions is the value of f₀(0). As we have seen, this quantity is also the mass accretion rate Ṁ. There is still no fundamental theory to supply this rate, except at β = 0(Bondi 1952), and in the hypersonic (β    ≫    1)limit (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944). Extending the work of Bondi (1952), Moeckel & Throop (2009) suggested an interpolation formula that both respects the analytic limits and agrees reasonably well with numerical simulations. As we did in Paper I, we adopt this prescription, which is $\begin{array}{ccc}\mathit{Ṁ}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}{\mathrm{(}{\mathit{\lambda }}^{\mathrm{2}}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\left(\mathrm{1}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}\right)}^{\mathrm{2}}}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}& & \end{array}$(51)Here, 2   λ = e^3/2/2 = 2.24is the analytic value of Ṁ(0) in the isothermal case (Bondi 1952).

We are now in a position to outline our numerical procedure. Starting at the upstream axis, θ = π , we set ${\mathit{f}}_{\mathrm{0}}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\pi }\mathrm{\right)}\mathrm{=}\mathrm{0}$ . We provisionally treat g_-2(π) as a free parameter, to be determined later in order to ensure smoothness of the flow across the Mach cones, as we describe in more detail below. For any selected g_-2(π), we may then integrate Eqs. (31) and (32) to the anti-Mach cone, $\mathit{\theta }\mathrm{=}{\mathit{\theta }}_{\mathit{M}}^{\mathrm{\prime }}$ . Simultaneous integration of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ yields f₀ itself.

But knowledge of both g_-2(π) and f₀(0) also gives us g_-2(0), according to Eq. (50). Since ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}$ , we may again integrate the second-order Eqs. (31) and (32), along with ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$, from θ = 0to θ = θ_M . At this point, we have established a one-parameter family of flows covering both the upstream and downstream regions. For any value of g_-2(π), both flows are irrotational, as required physically.

Turning to the intermediate region, we begin at the equatorial plane, θ = π/2 . We use Eq. (47) for g_-2(π/2) , again setting C =  +1 . Our free parameter is now ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{(}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{)}$. For any value of this quantity, we integrate Eqs. (33) and (34) away from the plane in both directions until we come to the Mach cones. We have then established the run of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 throughout the intermediate region. Again, the flow is irrotational for any value of ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{(}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{)}$⁴.

4.5.3. Enforcing smoothness

Within the intermediate region, we only know f₀ up to a constant of integration, which has yet to be fixed. In addition, we have not yet determined our two free parameters, ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{(}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{)}$ and g_-2(π). All three quantities are specified by requiring that the flow be smoothly varying. If ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 were well-defined everywhere, this task would be straightforward. For example, we could tune one parameter until g_-2, as calculated in the downstream region, matched the intermediate g_-2 at the Mach cone. However, both ${\mathit{f}}_{\mathrm{0}}^{\mathrm{\prime }}$ and g_-2 diverge at the cones. Hence, we can only require smoothness outside some finite region surrounding each cone.