ASYMPTOTIC OPENING ANGLES FOR COLLIDING-WIND BOW SHOCKS: THE CHARACTERISTIC-ANGLE APPROXIMATION

When numerical results, rather than analytic forms, for the shape of a thin and steady shock are desired, the standard approach for determining the shock opening angle involves a detailed integration of the shock conditions all along its surface, starting from the stagnation point along the symmetry axis (when coriolis effects are neglected). The physical requirement for such an integration is that the shock cone has a "memory" of the history of the gas that it is carrying, encoded in the form of the mass, energy, and momentum fluxes along the shock cone, and the point where any windstream from one of the stars meets the working surface of the shock cone depends on that complete history—it is not a local attribute of that windstream. However, a simplification can be applied if it is only the asymptotic opening angle that is desired, not the detailed structure. Then the result should be expressible in a global form, breaking from the standard locally integrated numerical approach, and allowing the problem to be solved in a single step.

In other words, the opening angle should be solvable by applying purely global constraints, but so far that has only been possible for highly radiative thin shocks, as assumed in CRW. We suggest that a systematic, if idealized, approximation that allows for such global constraints follows from a simple extension of the thin-shock approximation to include adiabatic cooling, but which conserves different wind properties consistent with the presence of explosive gas pressure and the absence of significant radiative cooling. Then contrasts between the predicted opening angle in the adiabatic and radiative limits, for both mixed and unmixed winds, focus squarely on the effects of adiabaticity, and guide more complete hydrodynamic studies of these differences.

In the process, we also seek to clarify some potential misconceptions about the relative importance of the history of the mass loading of the shocked region, as opposed to the local constraints stemming from the collision of two individual windstreams without consideration of that history. These issues relate to the concept of centrifugal corrections and the "Dyson approximation" (Dyson 1975), so we begin our analysis there.

An approximation that has been applied in many contexts (e.g., Luo et al. 1990; Stevens et al. 1992; Antokhin et al. 2004) to highly supersonic bow shocks, is requiring a balance between the ram pressure perpendicular to the shock front in the local colliding windstreams (Dyson 1975). This amounts to neglecting the momentum requirements of turning the inertial flow already moving along the shock front, so is termed the neglect of "centrifugal" corrections. Explorations into the errors introduced by this approach (e.g., Mac Low et al. 1991) typically focus on the behavior near the head of the bow shock and on the generation of the hardest X-rays, but the approximation may lose validity farther downstream where the global shock morphology is determined.

It will be argued here that although this approximation is generally valid near the stagnation point, since there the history of the introduction of mass into the shocked region is of least importance, it is not a particularly useful approximation to apply to the asymptotic shape of the shock, as the latter is more related to globally integrated geometric factors that incorporate the history of the mass loading of the shock. As such, the Dyson approximation is of more value for analyzing high-temperature emissions from near the shock apex, than it is useful for understanding the more global shock-cone shape. On the other hand, the global character of our approach cannot produce local emission diagnostics near the shock apex, so is intended to be complementary to the use of the Dyson approximation, and applicable only to morphological observations.

The key global constraint is overall conservation of the components of vector momentum flux, tracking the wind inputs as well as any external exchange of momentum, in axial wedges conforming to the assumed azimuthal symmetry. The symmetry axis is taken to be the binary line of centers, so applies only on scales where the orbital effects which break that symmetry may be neglected. As will become more apparent below, tracking the global momentum flux is substantially different from asserting local perpendicular ram balance, because the shocked gas carries a momentum flux that must be deflected more so than stopped, thereby "soaking up" locally unbalanced ram pressure, as the shock angle turns from its initially perpendicular angle to its ultimate asymptotic angle. This effect represents a type of mass loading, wherein the stronger wind piles into previously shocked elements of itself, and that allows for both global and local violation of the perpendicular ram balance. This integrated memory of the interaction, not local physics, determines the asymptotic opening angle of the shock. As such, this paper explores the fundamental conceptual reason why the Dyson approximation is not suited to the asymptotic shock geometry, and why this is further exacerbated in the presence of adiabaticity.

Our core assumption, for the purposes of achieving simple closed-form results, is that the flow of the shocked gas is governed by a single characteristic angle, θ_s, which forms the separatrix of the two winds downstream of the interaction region. The basic assumed geometry, and the meaning of θ_s in any azimuthal wedge dϕ, is indicated in Figure 1. There, θ_s is interpreted as the opening half-angle of the asymptotic bow shock, to within the limitations of the approximation.

Zoom In Zoom Out Reset image size

This is also the formal limit considered by CRW, which they justify by assuming rapid radiative cooling and lateral confinement of the shocked gas between the two unshocked winds. In our case, this assumption is not intended to be taken as literally as in the thin-shock approximation, because when gas pressure is significant, as with adiabatic shocks, such lateral confinement is not physically realizable. Nevertheless, the return of thermal energy to bulk flow energy, endemic to adiabatic cooling, is assumed to result in the bulk of the interacting gas exiting the system more or less along a single most prominent or characteristic direction, and it is that characteristic direction we seek to understand, even though the shock is not formally thin. Future hydrodynamic simulations are planned to explore the applicability of this approximation, but at this stage we view it primarily as an asymptotic benchmark, to enable systematic comparisons between radiative and adiabatic shocks.

We now lay out our fundamental treatment of the quasi-thin shock approximation, in terms of conserved fluxes of vector momentum components in the z and ρ directions. The key addition to the CRW treatment is the inclusion of gas pressure, capable of generating a new source of momentum flux away from the binary line of centers, as occurs in the adiabatic limit. In axial symmetry, gas pressure does not alter the integrated z component of momentum flux within each axial wedge under consideration, because the z direction is globally constant in cylindrical coordinates, and there are no external forces in the z direction. However, gas pressure in the shocked gas can and does increase the flux of the component of momentum that points perpendicular to the axis, which we denote as the ρ direction in cylindrical coordinates. It is this accounting of the ρ component of the momentum flux that allows us to track the influences of adiabaticity.

We consider a sphere at very large radius centered on the binary system, and track the flux in the colliding winds of various scalar and momentum quantities through that sphere. Owing to the axial symmetry, we may restrict to a thin wedge, of axial angular width dϕ and with its vertex all along the line of centers of the binary, shaped like the section of an orange. Considering such identical wedges avoids cancellation when summing the flux of vectors with ρ components. We may then apply appropriate global conservation laws, resulting from the given flux sources in the spherical winds of each star, and the absence of any external forces (for simplicity we neglect gravity and radiative forces and assume the winds appear at their terminal speeds), to constrain the nature of the fluxes through the sphere, independently of any wind/wind interaction. Hence, we obtain quantities that must be the same in the actual case as they are in an imaginary case where the winds do not interact at all, and these quantities supply us with suitable constraints for deriving the desired expression for the characteristic shock angle θ_s.

As mentioned above, we define the characteristic angle along which the shocked gas flows to be θ_s, and let P_s be the scalar momentum flux along that angle (a sum over both shocked winds). Purely for convenience, we normalize the scalar momentum flux of star A, which is the mass-loss rate ${\dot{M}}$ times the terminal speed v_A, to be 8π. As all momentum fluxes in this paper are expressed per angular width dϕ of the azimuthal wedge under consideration, this implies that they may all be converted to real momentum flux units by multiplying them by ${\dot{M}}_A {v_A}d\phi /8\pi$.

The essential device we use is to consider separately the flux of the component of momentum in the z direction (i.e., the component along the binary line of centers), the flux of the component of momentum in the ρ direction (along the radius of the wedge), and where possible, the scalar momentum flux (the flux of the magnitude of the momentum being advected). Due to the absence of any external forces on the wedge in the z direction, we may assert equality of the flux of the z component of momentum, between the actual wedge, and an imaginary wedge where no wind interactions of any kind occur. This conservation of z-component of momentum is written

$$\begin{equation} P_{s}\cos {\theta _s}\ = \ P_{A,z}({\theta _s}) \ + \ P_{B,z}({\theta _s}), \, \end{equation} \tag{ 1 }$$

where again P_s is the scalar momentum flux summed over both interacting winds and assumed to be flowing characteristically along θ_s. Here,P_A,z(θ_s) and P_B,z(θ_s) are

$$\begin{equation} P_{A,z}({\theta _s}) \ = \ 2 \int _0^{\theta _s}d\theta \ \sin \theta \cos \theta \ = \ \sin ^2 {\theta _s} \end{equation} \tag{ 2 }$$

and

$$\begin{equation} P_{B,z}({\theta _s}) \ = \ 2 \eta \int _{\theta _s}^\pi d\theta \ \sin \theta \cos \theta \ = \ - \eta \sin ^2 {\theta _s}, \, \end{equation} \tag{ 3 }$$

which are respectively the flux of the z component of momentum that is embroiled in the shock from wind A and from wind B, for η the ratio of the total wind momentum fluxes in the weak wind (B) to the strong wind (A). In the latter expressions, the sin θ is from the azimuthal extent of the solid angle contributing along θ, and cos θ gives the desired projection onto the z-component of the momentum.

It is worth special mention of the reasons why gas pressure in the shocked winds cannot alter the z-component of the momentum flux balance. Although the forces induced by gas pressure are mediated by pressure gradients in the shocked gas, we need not consider such gradients in detail, opting instead to simply track the conserved momentum-flux quantity. Because the system as a whole is closed, and may be translated without consequence in the z direction, the combined shocked and unshocked gas must exhibit zero total z-component momentum flux. Thus, the flux of the z-component of momentum summed over the shocked gas much equal the negative of the z-component of momentum summed over the unshocked gas. Then the fact that the shocked and unshocked gases cannot exchange any momentum between them (because of the supersonic nature of their interaction), the flux of the z-component of momentum in both the unshocked gas, and the shocked gas independently, must be the same as they would have been had there been no wind interaction at all. In particular, the presence of gas pressure within the shocked component cannot alter this total z-component momentum flux, regardless of the detailed nature of the pressure gradients it produces, and this simplification motivates considering the z-component, along the symmetry axis, separately from the radial momentum flux, away from the symmetry axis.

In contrast with the simple net outcome of the z-momentum balance, the situation in the ρ direction is complicated whenever gas pressure in the shocked winds is included in the analysis (so whenever the shocked gas does not rapidly radiatively cool). The fundamental difference stems from the symmetry of the azimuthal wedges we consider, which support no ignorable radial translation, and instead exhibit pressure forces on their boundaries that may be treated as an equal squeezing force on opposite faces of the wedge. These squeezing forces are perpendicular to the faces of the wedge, so when added vectorially, they must yield a net force, and a net increase in the flux of momentum, entirely in the ρ direction. Let this source of ρ momentum, per angular width dϕ, be denoted P_e. Then in a similar spirit to Equation (1), we have for the flux of the ρ component of momentum

$$\begin{equation} P_{s}\sin {\theta _s}\ = \ P_{A,\rho }({\theta _s}) \ + \ P_{B,\rho }({\theta _s}) \ + \ P_e, \, \end{equation} \tag{ 4 }$$

where

$$\begin{equation} P_{A,\rho }({\theta _s}) \ = \ 2 \int _0^{\theta _s}d\theta \sin ^2 \theta \ = \ {\theta _s}\ - \ \sin {\theta _s}\cos {\theta _s} \end{equation} \tag{ 5 }$$

and

$$\begin{equation} P_{B,\rho }({\theta _s}) \ = \ 2 \eta \int _{\theta _s}^\pi d\theta \sin ^2 \theta \ = \ \eta \left(\pi \ - \ {\theta _s}\ + \ \sin {\theta _s}\cos {\theta _s}\right) \end{equation} \tag{ 6 }$$

give the sources of ρ momentum from the two winds respectively, again integrated over the solid angle of embroiled gas appropriate for each wind.

Again, although it may at first glance seem that the P_e momentum flux in the ρ direction appears somewhat magically, it can be physically traced to the curved axial symmetry, because the azimuthal transport of vector momentum flux endemic to the high gas pressure in an adiabatic shock will inevitably lead to an enhancement of momentum flux in the ρ direction. This is related to the fact that if we use cylindrical coordinates to treat an arbitrary distribution of particles that are individually moving in a straight line, the particles will always experience an increase in their radial component of momentum, at the expense of their azimuthal momentum, with no change in their axial momentum.

Equations (1) and (4) may be viewed as a matrix equation in (η, P_s) in terms of the unknown P_e and the given θ_s, which solves to

$$\begin{equation} \eta \ = \ {\sin {\theta _s}\ {-} \ (P_e+ {\theta _s}) \cos {\theta _s}\over \sin {\theta _s}\ {+} \ (\pi - {\theta _s}) \cos {\theta _s}} \end{equation} \tag{ 7 }$$

and

$$\begin{equation} P_{s}\ = \ {(\pi \ {+}\ P_e) \sin ^2 {\theta _s}\over \sin {\theta _s}\ {+} \ (\pi - {\theta _s}) \cos {\theta _s}}. \end{equation} \tag{ 8 }$$

If we imagine that θ_s is known from observation, then the above represents two equations in the three unknowns η, P_s, and P_e, so solving them will require an auxiliary assumption about the energy transport. We will address that assumption in the form of limiting cases that deal with the degree of adiabaticity and the degree of wind mixing.

We first consider the simplest case of purely radiative shocks, for which gas pressure plays no role and we may set P_e = 0. This situation is handled in CRW; we include it here only for completeness. The solution from Equation (7) is immediately

$$\begin{equation} \eta \ = \ {\tan {\theta _s}\ {-} \ {\theta _s}\over \tan {\theta _s}\ {-} \ {\theta _s}\ + \pi }. \end{equation} \tag{ 9 }$$

Note that it is not necessary to make any assumptions about the degree of mixing in the two shocked winds, nor to what extent P_s is locally carried by a single speed or a spread in speeds within the shock, because our result for θ_s is a global attribute of the momentum balance independently of how that momentum is partitioned among shock components. The resulting η(θ_s) is depicted in Figure 2.

Zoom In Zoom Out Reset image size

In the opposite limit of negligible radiative cooling, we have energy conservation. As such, once the shock gas is essentially completely adiabatically cooled, its bulk-flow kinetic energy must return to its pre-shock value. This does not explicitly violate the need for entropy to rise, because the expanding gas will experience a dropping density, such that by the time the temperature is low, the density is also quite low, and the specific entropy can rise continuously, as required by the irreversible nature of the shock and the expansion.

Thus in the adiabatic limit, we assume that quasi-thin adiabatic shocks will ultimately convert all the heat thermalized in the shock back into bulk kinetic energy flowing approximately along θ_s, but the impact on P_s will depend on whether the two winds mix and reach a single combined velocity, or if they retain their individual character and return only to their original terminal speed, or some combination thereof. First we treat the case where each wind returns adiabatically to its original terminal speed without mixing. This means that the scalar momentum flux missing from the two original winds will all show up in the asymptotic flow along the shock angle, i.e., that scalar momentum flux will all end up being P_s. Hence instead of substituting for P_e directly, we may replace Equation (4) with the constraint

$$\begin{eqnarray} P_{s} &=& 2\int _0^{\theta _s}d\theta \, \sin \theta \ {+} \ 2 \eta \int _{\theta _s}^\pi d\theta \, \sin \theta \nonumber\\ &=& 2 (1-\cos {\theta _s}) \ {+} \ 2 \eta (1+\cos {\theta _s}). \end{eqnarray} \tag{ 10 }$$

This results in the solutions to Equations (7) and (8) becoming

$$\begin{equation} \eta \ = \ \tan ^4 \left({{\theta _s}\over 2} \right) \end{equation} \tag{ 11 }$$

and

$$\begin{equation} P_{s}\ = \ 2 (1-\cos {\theta _s}) \ + \ 2 \tan ^4 \left({{\theta _s}\over 2} \right) (1+\cos {\theta _s}) \, \end{equation} \tag{ 12 }$$

subject to

$$\begin{eqnarray} P_e\ &=& \ {1 \over 8} {\rm sec} ^4 \left({{\theta _s}\over 2 } \right) [ 8 \sin {\theta _s}\ {-} \ \pi \cos (2{\theta _s}) \nonumber\\ && + 4 (\pi -2{\theta _s}) \cos {\theta _s}\ - \ 3 \pi ]. \end{eqnarray} \tag{ 13 }$$

The resulting η(θ_s) is depicted in Figure 2.

If we make the opposite assumption that when the thermalized energy is adiabatically returned to the flow along θ_s, the winds are completely mixed and reach a single joint characteristic velocity, then we get the maximal explosive support of the opening angle θ_s. The additional support comes from the fact that not only is the ram pressure perpendicular to the shock thermalized, but even some of the ram pressure initially along the shock is thermalized, when the winds mix and reach a common speed. This is essentially an increase in gas pressure due to frictional heating, and so we expect the smallest η(θ_s) in this case. Here our global constraint on P_s comes from the total scalar kinetic energy flux and the total scalar mass flux that enter the shock zone, which must in turn flow out at a single characteristic speed, approximately along θ_s, in a manner consistent with P_s.

Following this logic, the scalar mass flux per wedge thickness dϕ, again in units where the total scalar momentum flux from star A is 8π, is

$$\begin{equation} M_s \ = \ {2 \over {v_A}} (1-\cos {\theta _s}) \ + \ {2 \eta \over {v_B}} (1+\cos {\theta _s}) \, \end{equation} \tag{ 14 }$$

and the scalar kinetic energy flux in those units is

$$\begin{equation} K_s \ {=} \ { P_{s}^2 \over 2 M_s} \ {=} \ (1 - \cos {\theta _s}) {v_A}\ {+} \ \eta (1+\cos {\theta _s}) {v_B}, \end{equation} \tag{ 15 }$$

which results in the constraint

$${\fontsize{8.5}{9.5}\selectfont{\begin{eqnarray} &&P_{s}\ {=} \ 2 \sqrt{ [1 \ - \ \cos {\theta _s}\ + \ (1+\cos {\theta _s}) \eta u] \left[1 \ - \ \cos {\theta _s}\ + \ (1+\cos {\theta _s}) {\eta \over u}\right] },\nonumber\\ \end{eqnarray}}} \tag{ 16 }$$

where we have defined u = v_A/v_B. Applying Equations (7) and (8) then gives

$${\fontsize{8}{9.5}\selectfont{\begin{eqnarray} &&\lefteqn{ \eta \ = \ {\cos ^2 {\theta _s}\over 8 u (3 \cos {\theta _s}\ - \ 1)} {\rm sec}^6 \left({{\theta _s}\over 2} \right) [2(1+u^2)\cos ^2 {\theta _s}}\nonumber\\ &&- 2(1+u^2)-u \sin ^2 {\theta _s}\tan ^2 {\theta _s}\nonumber \\ &&+ \sqrt{2}\sqrt{1+u+4u^2+u^3+u^4+(1-u)^2(1+u+u^2)\cos (2{\theta _s})} \sin {\theta _s}\tan {\theta _s}],\nonumber\\ \end{eqnarray}}} \tag{ 17 }$$

and the expressions for P_s and P_e may also be written in closed form but they are quite long and involved. Although it is not immediately obvious, the result in Equation (17) does indeed reduce to Equation (11) when u = 1, since then the presence or absence of mixing is irrelevant to the global dynamics. This result for η(θ_s) is also included in Figure 2.

The above results simplify in the limit θ_s ≪ 1, which occurs when η ≪ 1, so it is informative to consider what physical insights may be conferred in that simple limit. Note these limits apply asymptotically only when the weak-wind star has a negligibly small radius; in real situations, when η is small enough for asymptotic expressions to apply, the possibility must be considered separately that the strong wind may crash directly into the photosphere of the companion, invalidating our assumptions. Nevertheless, the low-η limits do convey general insights into the reasons that different degrees of adiabaticity require different amounts of weak-wind momentum to support the shock cone at a given opening angle θ_s, and it is one of the primary advantages of closed-form expressions that they submit to this type of scaling analysis.

We begin by noting that if one adopts the local approximation of equating the perpendicular momentum fluxes across the shock front (e.g., Dyson 1975; Luo et al. 1990; Stevens et al. 1992; Antokhin et al. 2004; Falceta-Goncalves et al. 2008) and extrapolates it globally to the asymptotic opening angle, one should expect the required η to scale like θ_s² when θ_s ≪ 1. This is because the strong wind effectively has no inertia after it is shocked, so the weak wind bears the full burden at every point along the shock of maintaining that shock against the strong wind ram pressure. The weaker wind can maintain a perpendicular ram balance against at most a solid-angle fraction of order η of the strong wind, and θ_s² determines the solid-angle fraction subtended by the interaction zone, so η ∝ θ_s².

However, when a more accurate global accounting of the momentum requirements is undertaken, in the limit of radiative shocks we find from Equation (9) that

$$\begin{equation} \eta \ \cong \ {{\theta _s}^3 \over 3 \pi }, \, \end{equation} \tag{ 18 }$$

as seen in the exact result of CRW and the heuristic fit of Eichler & Usov (1993). The extra power of θ_s may be interpreted as being due to the fact that the global requirement for the weak wind is not to stop the strong wind, but merely to deflect it through an angle θ_s. So the η momentum flux must deflect through an angle θ_s a fraction ∼θ_s² of the strong wind, requiring η ∼ θ_s³. Indeed, we can be even more quantitative and note that the strong-wind momentum flux along the z direction is $2 \int _0^{\theta _s}d\theta \ \sin \theta \cos \theta \cong {\theta _s}^2$, and the deflecting momentum flux in the ρ direction from the weak wind is $2\eta \int _{\theta _s}^\pi d\theta \sin ^2 \theta \cong \pi \eta$, along with a contribution in the ρ direction from the strong wind itself of $2\int _0^{\theta _s}d\theta \ \sin ^2 \theta \cong 2{\theta _s}^3/3$. Adding the momentum flux in the ρ direction to that in the z direction bends the total shocked momentum flux an angle θ_s, where from the above estimates we have

$$\begin{equation} {\theta _s}\ \cong \ {\pi \eta \over {\theta _s}^2} \ + \ {2 \over 3} {{\theta _s}^3 \over {\theta _s}^2}, \, \end{equation} \tag{ 19 }$$

which results directly in η ≅ θ_s³/3π.

For adiabatic shocks with no mixing, our approximation in Equation (11) yields when θ_s ≪ 1

$$\begin{equation} \eta \ \cong \ {{\theta _s}^4 \over 16}. \end{equation} \tag{ 20 }$$

Here, we find yet another added power of θ_s, this time because the explosive heating of the shocked winds, upon re-expansion away from the axis, provides substantial momentum support for the bending of the shock angle. Indeed, the combination of the pre-shocked momentum flux in the ρ direction, and the explosive gas pressure contribution in that direction, produce almost enough momentum to support a small deflection θ_s by themselves, without help from the weak wind. Hence, the weak wind needs to provide only a small "coaxing" on top of these two momentum fluxes along ρ (each which scales ∼θ_s³), so this coaxing appears at an even higher order of θ_s, at order θ_s⁴. So this analysis elucidates the physical reasons why radiative shocks require less weak-wind momentum flux to maintain a given narrow shock angle than would a putative perpendicular ram balance, and adiabatic shocks require less still.

Considering the case of adiabatic shocks with complete mixing, we find from Equation (17) that when θ_s ≪ 1,

$$\begin{equation} \eta \ {\cong} \ {{\theta _s}^4 u \over 8 (1+u^2)}, \, \end{equation} \tag{ 21 }$$

where again u is the ratio of the terminal speeds in the two winds, and it does not matter which wind is in the numerator, only the contrast expressed by u. Here we see the now-familiar θ⁴ scaling of adiabatic shocks, but we also find that when there is a strong contrast in the wind speeds, adiabatic mixing allows additional thermalization and additional explosive expansion away from the axis, further supporting the deflection of the strong wind and allowing for an even smaller η to suffice, in light of the identity u/8(1 + u²) < 1/16.