text:
# Non-uniform chiral phase in effective chiral quark models

We analyze the phase diagram in effective chiral quark models (the Nambu–Jona-Lasinio model, the $`\sigma `$-model with quarks) and show that at the mean-field level a phase with a periodically-modulated chiral fields separates the usual phases with broken and restored chiral symmetry. A possible signal of such a phase is the production of multipion jets travelling in opposite directions, with individual pions having momenta of the order of several hundred MeV. This signal can be interpreted in terms of disoriented chiral condensates.

H. Niewodniczański Institute of Nuclear Physics, PL-31342 Kraków, Poland

PACS: 12.39.Fe, 21.65.+f, 12.38.Mh, 64.70.-p

The problem of phase transitions in hadronic matter is one of the most important issues in strong-interaction physics. It has been studied theoretically for many years with a variety of methods, with the results gaining importance at the approach of the operation of RHIC. We now believe that in the ($`T`$-$`\mu `$) diagram, where $`T`$ denotes the temperature and $`\mu `$ the baryon chemical potential, one can recognize a quite complicated structure: phases with broken/restored chiral symmetry, confined/deconfined phases, color superconducting phase (in the case of two light flavors) , or color-flavor locked phase (in the case of three light flavors) . In certain limits detailed features of phase transitions are known from first-principle QCD calculations (color superconducting gaps at asymptotically-large baryon densities ), or general symmetry considerations (e.g. the universality class of the chiral phase transition  in the massless limit). Lattice calculations provide information for finite-temperature QCD at zero baryon density . However, these “exact” methods are inapplicable to the case of moderate baryon densities (up to a few times the nuclear saturation density), where out of necessity we have to rely on models . In fact, a lot of our experience comes from model calculations, just to mention the expectation that chiral symmetry is eventually restored as the baryon density is increased.

In this letter we consider the chiral phase structure for a class of low-energy chiral models: the Nambu –Jona-Lasinio model  and the linear $`\sigma `$-model with quarks . These models have been used extensively over the past years to describe low-energy properties of mesons and baryons. They have also been used to describe dense and hot system of the Fermi gas of quarks , in particular to study chiral restoration at high $`T`$ or $`\mu `$. What is usually claimed is that at the mean-field level (synonymous to one-quark-loop or leading-$`N_c`$) the model has two phases: with broken chiral symmetry, i.e. large dynamically-generated quark mass $`M`$, which appears at low $`T`$ and $`\mu `$, and with restored chiral symmetry, with massless (or almost massless) quarks, at higher values of $`T`$ or $`\mu `$. We argue that the situation is much more complicated if one allows for more general mean-field solutions of the model. We show that for a certain range of $`T`$ and $`\mu `$ there appears a non-uniform phase with broken rotational and isospin symmetry. In the $`T`$-$`\mu `$ diagram this phase separates the usual chiral-symmetry-broken phase from the chiral-symmetry-restored phase. The results for the linear $`\sigma `$-model with quarks at $`T=0`$ were obtained a long time ago by Kutschera, Kotlorz, and one of us (WB) . Here we generalize the calculation to finite temperatures, as well as apply it to the Nambu–Jona-Lasinio model.

The Lagrangian densities of the $`\sigma `$-model and of the partially-bosonized Nambu–Jona-Lasinio model have the form $`\overline{\psi }[i/g(\sigma +i\gamma _5\tau \pi )]\psi +L_\mathrm{m}(\sigma ,\pi )`$, where $`\psi `$ denotes the quark fields, $`g`$ is the quark-meson coupling constant, and the mesonic part $`L_\mathrm{m}(\sigma ,\pi )`$ is specific to the model. We consider spatially non-uniform, time-independent chiral fields:

$$\sigma (\stackrel{}{x})=\frac{M}{g}\mathrm{cos}(\stackrel{}{q}\stackrel{}{x}),\pi _0(\stackrel{}{x})=\frac{M}{g}\mathrm{sin}(\stackrel{}{q}\stackrel{}{x}),\pi _\pm (\stackrel{}{x})=0,$$
(1)
where $`M`$ and $`\stackrel{}{q}`$ are parameters, which are determined dynamically. Such an ansatz has been considered for the first time in Refs.  in context of the pion condensation in nuclear matter. In the case of vanishing $`\stackrel{}{q}`$ we recover $`\sigma =\frac{M}{g}`$, $`\stackrel{}{\pi }=0`$. The Dirac equation with fields (1) can be solved exactly , which facilitates greatly practical  and theoretical  applications. The quark spectrum has two branches:

$$E_\pm (\stackrel{}{p})=\sqrt{M^2+\stackrel{}{p}^2+\frac{1}{4}\stackrel{}{q}^2\pm \sqrt{M^2\stackrel{}{q}^2+(\stackrel{}{q}\stackrel{}{p})^2}},$$
(2)
where $`\stackrel{}{p}`$ labels the momentum of the quark. The basic quantity of our study is the grand thermodynamical potential density $`\mathrm{\Omega }(T,\mu ;M,\stackrel{}{q})`$. Its global minimum with respect to $`M`$ and $`\stackrel{}{q}`$ at fixed $`T`$ and $`\mu `$ determines the ground state of the system subject to conservation of baryon number. We stress that our calculation is self-consistent, i.e. both quark and meson fields are solutions to the equations of motion. At the mean-field (one-quark-loop) level we have explicitly

$`\mathrm{\Omega }(T,\mu ;M,\stackrel{}{q})=\mathrm{\Omega }_0(M,\stackrel{}{q})`$ (3)
$`2N_cT{\displaystyle \underset{i=\pm }{}}{\displaystyle \frac{d^3p}{(2\pi )^3}\mathrm{ln}\left[\left(1+e^{(\mu E_i(\stackrel{}{p}))/T}\right)\left(1+e^{(\mu +E_i(\stackrel{}{p}))/T}\right)\right]},`$
The medium part of $`\mathrm{\Omega }`$ describes the ideal gas of quarks with the spectrum (2). The vacuum part, $`\mathrm{\Omega }_0`$, depends on the model considered. In the $`\sigma `$-model

$$\mathrm{\Omega }_0^\sigma (M,\stackrel{}{q})=\frac{m_\sigma ^2F_\pi ^2}{8}(M^2/M_0^21)^2+\frac{F_\pi ^2}{2}(M^2/M_0^2)\stackrel{}{q}{}_{}{}^{2},$$
(4)
where the parameter $`M_0=gF_\pi `$ is the vacuum value of the constituent quark mass, $`F_\pi =93\mathrm{M}\mathrm{e}\mathrm{V}`$ is the pion decay constant, and $`m_\sigma `$ is the vacuum value of $`\sigma `$-meson mass, treated as a model parameter. The first term in (4) comes from the Mexican Hat potential, while the second term is the kinetic energy of the meson fields. In the Nambu–Jona-Lasinio model we have

$$\mathrm{\Omega }_0^{\mathrm{NJL}}(M,\stackrel{}{q})=(M^2M_0^2)/(4G)+V(M)+\frac{F_\pi ^2}{2}(M^2/M_0^2)\stackrel{}{q}{}_{}{}^{2}+𝒪(\stackrel{}{q}{}_{}{}^{4}),$$
(5)
where $`G`$ is the four-quark coupling constant, and $`V(M)`$ is the Dirac sea contribution to the energy for the case $`\stackrel{}{q}=0`$. Here we adopt the simple 3-momentum regulator, hence $`V(M)=2N_c/\pi ^2_0^\mathrm{\Lambda }𝑑pp^2\left(\sqrt{p^2+M^2}\sqrt{p^2+M_0^2}\right)`$ , where $`\mathrm{\Lambda }`$ is the cut-off parameter. In the calculation presented below we use the parameters of Ref. : $`G=5.01\mathrm{GeV}^2`$, $`\mathrm{\Lambda }=650\mathrm{M}\mathrm{e}\mathrm{V}`$. Equation (5) has been gradient-expanded in powers of $`\stackrel{}{q}`$. The coefficient of the quadratic term is model independent. Higher-order terms depend on the value of $`\mathrm{\Lambda }`$. However, their contribution is small for moderate values of $`\stackrel{}{q}`$  and we can safely drop them. With this simplification the two models become very similar. The same conclusion holds for other variants of the Nambu–Jona-Lasinio model, which differ by the regularization method. In fact, the basic difference between the models is the difference of values of $`\mathrm{\Omega }_0`$ at the minimum and at the point $`M=0`$, $`\stackrel{}{q}=0`$. The larger this difference, the more difficult it is to restore chiral symmetry.

The essence of the dynamics of our system can be understood as follows: as we increase the chemical potential, it becomes favorable for the system to develop a non-zero $`\stackrel{}{q}`$. This is because the quark energies on the $`E_{}`$ branch lower with $`\stackrel{}{q}`$, and the energy gain overcomes the repulsion of the meson kinetic term in $`\mathrm{\Omega }_0`$. Thus, for certain $`T`$ and $`\mu `$ the ground state has finite $`\stackrel{}{q}`$. At the further increase of $`\mu `$ the mass of the quark drops to $`0`$, and chiral symmetry is restored. The phase diagram for the Nambu–Jona-Lasinio model with parameters of Ref.  is shown in Fig. 1(a). We can see three phases: at low $`\mu `$ and $`T`$ we have the uniform chirally-broken phase, with $`M>0`$ and $`\stackrel{}{q}`$ $`=0`$, wrapped around it is the phase with $`M>0`$ and $`\stackrel{}{q}`$ $`0`$ , and at high $`\mu `$ or $`T`$ we find the chirally-restored phase, with $`M=0`$. The lower line in Fig. 1(a) describes the first-order phase transition between the uniform and non-uniform chirally broken phases. The discontinuities in the $`M`$ and $`\stackrel{}{q}`$ parameters can be read-off from Fig. 1(c). Correspondingly, in the baryon density – temperature diagram (Fig. 1(b)) we notice the region of equilibrium of the two phases (the region between the two lower lines). The upper curve in Fig. 1(a) shows the second-order phase transition between the non-uniform chirally broken phase to the chirally-restored phase. At this phase transition the order parameter $`M`$ vanishes at non-zero $`\stackrel{}{q}`$. Interestingly, the character of the transition between the uniform and non-uniform chiral phases depends on the parameters of the models. In the $`\sigma `$-model for $`m_\sigma =800`$MeV we find a very similar diagram as in Fig. 1. However for $`m_\sigma =1200`$MeV the phase transition is second-order. It changes character around $`m_\sigma =1040`$MeV.

There is always the question of applicability and reliability of a model calculation like ours. Firstly, we need to have sufficiently high baryon densities or temperatures to trust the calculation. Physical systems at low densities consist of nucleons, and we need to dissolve them into quarks for the model to be applicable. By geometrical arguments we may expect this to happen at densities of the order of a few nuclear-matter saturation densities. Note that in Fig. 1 this is the region of the non-uniform chirally-broken phase, thus we may hope that the model is valid there. Secondly, the size of the system should be sufficiently large such that our infinite-size calculation is meaningful. It should be larger than $`2\pi /|\stackrel{}{q}|`$ in order to accommodate at least one period of the chiral field. Using the values of $`|\stackrel{}{q}|200600`$ MeV this gives the minimum size in the range $`26`$ fm, easily accessible in heavy-ion collisions. In addition, a more elaborate calculation should incorporate meson fluctuations. At low temperatures pions are easily excited due to their small mass, and such effects should certainly be included. Note, however, that inclusion of meson loops in the study of the quark condensate in Ref.  modified, but did not invalidate the mean-field results. Finally, the effect of the explicit breaking of chiral symmetry by the finite current quark mass, $`m`$, should be incorporated. In studies with $`\stackrel{}{q}=0`$  this effect was not quantitatively small. Moreover, with finite $`m`$ chiral restoration is no longer a second-order phase transition, but becomes a smooth cross-over. Another question is whether within the mean-field approximation the solution with fields (1) is the ground state. What we have verified is that for a certain range of $`T`$ and $`\mu `$ we have a lower-energy state than the uniform phase. A priori, there may exist yet lower-energy, up-to-now unknown, mean-field solutions.

We conclude with a discussion of a possible signal of the non-uniform chirally broken phase. Suppose that sufficiently large domains of the non-uniform chirally-broken phase are created in the cooling process of the plasma formed in a heavy-ion collision. Such a domain can be treated as a classical source of pions, in full analogy to disoriented chiral condensates (see e.g. ). The coherent decay of such a source would lead to a characteristic signal in the form of two multi-pion “jets”, with pion momenta peaked around $`\pm \stackrel{}{q}`$ (in the rest frame of the source), where $`\stackrel{}{q}`$ is taken at the phase-transition point, and is typically $`200300\mathrm{M}\mathrm{e}\mathrm{V}`$. Indeed, according to Bjorken’s description, pion quanta $`\delta \varphi ^a(x)`$ satisfy the equation $`\left(\mathrm{}+m_\pi ^2\right)\delta \varphi ^a(x)=j^a(\stackrel{}{x})`$, where $`j^a(\stackrel{}{x})`$ is the source, in our case proportional to the classical pion field, i.e. to $`\mathrm{sin}(\stackrel{}{q}\stackrel{}{x})`$. Thus the Fourier transform of $`\delta \varphi ^a`$ picks up the components at $`\pm \stackrel{}{q}`$. Finite size of the domain spreads out the distribution, but if the domain has the size of, say, $`3`$ wavelengths $`2\pi /|\stackrel{}{q}|`$, the half-width of peaks is only $`20\%`$ of $`|\stackrel{}{q}|`$. Of course, larger domains lead to sharper peaks. If multiple domains contribute, then the signal should be folded with the momentum distribution of the domain and averaged over directions. In Eq. (1) we have used the neutral classical pion field. This case leads to production of neutral pion quanta. However, we can rotate ansatz (1) in isospin, which leads to a degenerate solution. Such a configuration will decay into jets of charged pions. In addition, the quanta of the $`\sigma `$ field decay subsequently into pairs on neutral as well as charged pions. A more detailed calculation of the spectrum of pions requires additional assumptions for the dynamics of the phase transitions and is left for a further study. Also, the analysis of the formation and stability of domains is outside of the scope of this paper.

We thank Marek Kutschera, Jacek Dziarmaga, and Wojciech Florkowski for many useful conversations. One of us (MS) acknowledges the support of the Polish State Committee for Scientific Research, grant 2P03B 086 14.


format: text

text:
# Modeling Non-Axisymmetric Bow Shocks: Solution Method and Exact Analytic Solutions

## 1 Introduction

Supersonic stellar winds shock the surrounding gas and drive expanding bubbles into the interstellar medium. These shocks provide an opportunity to probe the properties of both the driving stellar wind and the ambient medium. If the star is moving with respect to the interstellar gas, the bubble will be distorted into a cometary shape. When the stellar motion is supersonic, we refer to these as stellar wind bow shocks (Baranov, Krasnobaev & Kulikovskii 1971; Dyson 1975). Since the discovery of such bow shocks around young B stars (Van Buren & McCray 1988), bow shocks have been found associated with many classes of objects, such as pulsars (Kulkarni et al. 1992) and cataclysmic variables (e.g. Vela X-1: Kaper et al. 1997); examples include well-known naked-eye stars (e.g. Betelgeuse: Noriega-Crespo et al. 1997). Bow shocks have been proposed as an explanation for cometary, ultracompact HII regions (Van Buren et al. 1990; Mac Low et al. 1991) and as a means of explaining the lifetimes of ultracompact HII regions. In a recent survey of the IRAS database using HiRes processing, Van Buren, Noriega-Crespo, & Dgani (1995) found 58 candidate bow shocks.

Non-axisymmetric stellar wind bow shocks occur when a star with an anisotropic wind moves supersonically with respect to the local medium, or if the star has an isotropic wind but moves in an ambient medium containing a transverse density gradient. Models of non-axisymmetric bow shocks are relevant to cometary ultracompact HII regions due to wind-blowing O stars moving supersonically with respect to the surrounding molecular cloud, when the ambient material does not have a constant density. A non-axisymmetric bow shock has also been invoked to explain the morphology of Kepler’s supernova remnant, where the supernova ejecta collide with a non-axisymmetric bow shock generated by the pre-supernova wind (Bandiera 1987; Borkowski, Blondin & Sarazin 1992). Another example is the bow shock due to the head of a jet propagating into a region with a density gradient. Non-axisymmetric, ram-pressure balance models of the collision between a stellar wind and the photoevaporating flow from an externally illuminated circumstellar disk have been given by Henney et al. (1996). A formulation for steady-state non-axisymmetric bow shocks and colliding winds was given by Bandiera (1993). However, Bandiera’s numerical method is sufficiently complicated that a simpler, analytic method is desirable.

In this contribution, I present a method for solving the problem of steady-state, momentum-conserving, non-axisymmetric, thin-shell bow shocks and colliding winds. This is an extension of the previous analytic solution method of Wilkin (1996, hereafter Paper I) and of Cantó, Raga, & Wilkin (1996, hereafter Paper II) to non-axisymmetric problems (see also Wilkin 1997a).

An outline of the paper is as follows. In § 2, we formulate the problem of the steady-state collision of two winds, and in § 3 we treat the problem of a bow shock resulting from an isotropic wind interacting with a plane-parallel flow containing a transverse density gradient. In § 4, we allow for non-isotropic winds, especially an axisymmetric wind with random orientation of the symmetry axis with respect to the direction of stellar motion. The rate at which kinetic energy is thermalized for the bow shock is treated in § 5. Results and future directions of this research are summarized in § 6.

## 2 Mathematical Formulation and Solution Method

### 2.1 Description of the Collision Surface

The hypersonic collision of two winds will in general result in a system of two shocks. A specific example, that of a stellar wind bow shock, is shown schematically in Figure Modeling Non-Axisymmetric Bow Shocks: Solution Method and Exact Analytic Solutions. Because this paper is concerned with steady-state solutions, the colliding winds are assumed to be unchanging in time. The stellar wind bow shock arises when a wind-blowing star moves supersonically with respect to the intersellar gas. In this case, we formulate the problem in the reference frame of the star, so the ambient medium is described as a wind of parallel streamlines impinging on the bow shock. For the collision of two winds in a binary star, we will neglect orbital motion in order to consider a steady-state problem in an inertial reference frame. In steady-state, the amount of mass and momentum within a given volume does not increase, and a flow pattern exists between the two shocks that carries away the mass and momentum deposited by the colliding winds. If the shocks are radiative, post-shock cooling lowers the temperature of the gas and leads to a large compression. We make the strong thin-shell assumption that cooling is so efficient that the shocked shell collapses to an infinitesimally thin layer, with a finite surface density $`\sigma `$ of matter. For this thin-shell assumption to apply, it is also necessary that the post-shock gas not be supported by magnetic fields, which, if well-coupled to the gas, could maintain a finite thickness even in the presence of efficient cooling. We assume there are no other forces such as those due to radiation or magnetic fields. The two shocked winds may in principle be separated by a contact discontinuity. However, there would then be supersonic shear across the interface that is expected to be unstable, leading to a mixing of the two fluids. A detailed treatment of the mixing is beyond the scope of this work. Instead, the mixing is assumed to be instantaneous, so the shell will have a unique velocity $`𝐕_t`$ at any location within it. This velocity represents an average of the turbulent fluctuations that would be present in a more detailed treatment. In steady-state, the geometric shape of the shell is unchanging in time, so the velocity of matter within the infinitely thin layer must be purely tangent to the shell. There will, however, be acceleration normal to the shell, because the fluid typically does not follow a straight path. The assumption that the incident streams are hypersonic, combined with the perfect cooling assumption, means that the flow within the shell will also be hypersonic, and pressure forces may be neglected in describing the motions along the shell.111Clearly this assumption must break down for real, finite temperature systems near the stagnation point, because the tangential flow velocity in the shell vanishes at that point. However, pressure forces depend upon the gradient of the pressure, which also vanishes at the stagnation point, so the pressure forces may still be small compared to momentum deposition in the stagnation region (Wilkin 1997b). In any case, beyond a small region near the stagnation point, the tangential flow in the shell will be supersonic, provided the incident flows are supersonic and the shocked fluid cools efficiently. Defining a spherical coordinate system $`(r,\theta ,\varphi )`$, and denoting the radius of the surface by $`r=R(\theta ,\varphi )`$, a complete description of this idealized shell is given by specifying the quantities $`R,\sigma ,𝐕_t`$ as a function of position $`(\theta ,\varphi )`$ within the shell. There are in fact only four independent quantities, because the condition that the motion within the shell be tangential implies that one need only solve for two velocity components.

### 2.2 Previous Treatment of the Problem

The conservation laws of mass and momentum may be used to derive the properties of the shell in terms of those of the two incident winds. Since the momentum conservation law has three components, there will be four equations in four unknowns, or including the condition of zero normal velocity, five equations in five unknowns. Because the assumption of vanishing thickness eliminates one spatial variable, these equations will be partial differential equations (PDEs) in two spatial coordinates. Suitable boundary conditions must also be specified. In practice, this means identifying the stagnation (or standoff) point, where the two winds collide head-on, and where for steady-state, the ram pressures of the two fluids balance. The equations may then be integrated away from the standoff point, following the motion of a fluid element. In order to begin the numerical integrations, one generally finds that an expansion about the conditions at the stagnation point is needed.

This approach was taken by Bandiera (1993), who derived a set of PDEs in curvilinear coordinates matched to the shape of the shell. These equations were then solved under the assumption of radial, constant velocity, isotropic winds from two point sources. The moving stellar wind bow shock problem is formally obtained by taking the limit of one source placed infinitely far away, while allowing its mass loss rate to be infinite so as to produce a finite density near the other star (at the bow shock).

Bandiera noted that for the specific problems cited, the motion within the shell would be along planes. One may readily see this for the bow shock shown in Figure 1. Consider a point on the bow shock surface, and the plane of constant azimuthal angle containing it and the stellar trajectory. The radial wind striking the shell at this point has momentum lying in this plane. Similarly, the ambient medium striking this point has momentum lying in the plane. If the two incident streams mix instantaneously on impact, the resultant momentum must also lie in the same plane. As the fluid flows along the shell, it continues to incorporate momentum contributions lying in the same plane. The collision of radial winds emanating from two point sources gives the same result. Thus, given the perfect mixing asumption, the flow within the shell will lie along planes, a situation we refer to as meridional flow. If the shocked fluid does not mix, the fluid trajectories in the shell are not confined to a plane, and further work is needed to solve this more complicated problem.222In the case of axisymmetric problems with shear, one may show that for divergent flow fields such as those encountered in bow shocks and colliding winds, the flow geometry and the fluxes of mass, momentum and angular momentum are the same as for the mixed case (Wilkin, Cantó and Raga 2000).

The ensuing geometric simplification allowed Bandiera to construct a two-dimensional grid, consisting of neighboring trajectories within the shell at a set of azimuthal angles. His computational method was to integrate the PDEs for mass, normal momentum, and tangential momentum, by marching along the trajectories within the shell as if the equations were ODEs. By differencing neighboring solutions, numerical values of the cross-streamline derivatives were obtained and supplied to the ODE integrator. Provided the spacing of the two-dimensional grid is sufficiently fine, his method will yield accurate solutions. Bandiera postulated that the dependence of the equations on derivatives of quantities across streamlines was “fictitious,” although he did not succeed in eliminating it from the equations.

In this contribution, the solution method will be simplified to an integral approach that avoids the need for PDEs, solving purely algebraic equations. Additionally, the solution may conveniently be obtained in terms of ordinary Cartesian, cylindrical polar, or spherical polar coordinates, and does not require a coordinate system matching the shape of the shell.

### 2.3 The Solution Method

The solution method is based upon the observation that thin shells driven by hypersonic winds are momentum-conserving in the vector sense (Paper I). In order to conserve momentum in steady-state, the momentum flux in the shell must be the vector sum of the two incident momentum fluxes, integrated from the standoff point to the point of interest. Such an integration is performed over the area of the shell between two planes of constant azimuthal angle $`\varphi `$, in the limit that their separation $`\mathrm{\Delta }\varphi `$ is infinitesimal (Fig. 1). For axisymmetric flow with no rotational motion about the symmetry axis, no mass or momentum crosses these bounding planes, so the mass and momentum flowing within such a wedge must exactly equal the sum of the mass or momentum fluxes onto the external faces of the bounding surfaces (shocks). Given the known vector momentum flux, one may determine the shape of the shell, since the fluid must move in the direction of its momentum. While this description is complete, the mathematics was subsequently simplified with the inclusion of an additional, angular momentum flux integral (Paper II). In Papers I and II, only constant velocity, isotropic winds were considered. For such winds, the vector momentum flux incident onto the shocked shell is independent of the detailed shape of the shell, so it is possible to specify the flux of momentum onto the shell from each side analytically. The methods of Papers I and II will now be extended to non-axisymmetric bow shocks and colliding winds, provided the flow is meridional. For meridional flow, the bounding surfaces are again planes of constant $`\varphi `$, allowing us to integrate the external fluxes and determine the internal fluxes within the shell. For anisotropic, radial, constant velocity winds, such as those considered by Bandiera, the momentum flux will depend only upon the number of streamlines intersected, and we will be able to obtain solutions analytically.

Proceeding with the solution method, we note the location of the stagnation point, where the two streams collide head-on. For the two-wind collision, this point lies between the two stars, along the line connecting their centers. For the bow shock due to a wind-blowing star moving with respect to the ambient medium, it will lie on the stellar trajectory, in the direction of stellar motion. For either problem, the z-coordinate axis is chosen to contain the stagnation point at $`\theta =0`$, with the coordinate origin at the wind source (Fig. 1). In what follows, the wind from the coordinate origin will be referred to as “the wind”, while the “second wind” may be either a radial wind from a second point source, or the ambient medium in the frame of the moving star, where the ambient velocity is $`𝐕_a=V_a\widehat{𝐳}`$. However, we will show explicit formulas only for the bow shock, reserving detailed treatment of the binary star colliding winds for a future paper. One should bear in mind that the methods shown below will be applicable as well for that problem. Returning to the description of a thin slice of the shell bounded by planes of constant $`\varphi `$, the mass, momentum, and angular momentum flux functions per unit azimuthal angle are defined by

$$\mathrm{\Phi }_m=\varpi \sigma V_t\mathrm{sec}\eta ,$$
$`(1)`$
$$𝚽=\mathrm{\Phi }_mV_t\widehat{𝐭},$$
$`(2)`$
$$𝚽_J=𝐑\times 𝚽.$$
$`(3)`$
Here $`\varpi =R\mathrm{sin}\theta `$ is the cylindrical radius, and the unit vector at constant $`\varphi `$ tangent to the shell is $`\widehat{𝐭}=\widehat{\varphi }\times \widehat{𝐧}/|\widehat{\varphi }\times \widehat{𝐧}|`$, where $`\widehat{𝐧}`$ is the unit outward normal to the shell. Noting that the arc length element $`ds`$ traced along the shell is given by $`(ds)^2=(dR)^2+R^2(d\theta )^2+R^2\mathrm{sin}^2\theta (d\varphi )^2,`$ the arc length element traced along the shell at constant $`\theta `$ is $`\varpi \mathrm{sec}\eta d\varphi `$, where the angle $`\eta `$ is given by

$$\mathrm{tan}\eta =\frac{R}{\varphi }/R\mathrm{sin}\theta .$$
$`(4)`$
For example, the flux of mass crossing a surface of constant polar angle $`\theta `$, between the azimuthal angles $`\varphi `$ and $`\varphi +d\varphi `$, is $`\mathrm{\Phi }_md\varphi `$. If $`V_\theta `$ is the $`\theta `$-component of velocity within the shell, and we write the cylindrical components of the momentum flux function as $`\mathrm{\Phi }_\varpi `$ and $`\mathrm{\Phi }_z`$, then these components are related to the angular momentum flux by

$$𝚽_J=\mathrm{\Phi }_mRV_\theta \widehat{\varphi },$$
$`(5)`$
$$𝚽_J=R(\mathrm{\Phi }_\varpi \mathrm{cos}\theta \mathrm{\Phi }_z\mathrm{sin}\theta )\widehat{\varphi }.$$
$`(6)`$
For meridional flow, the angular momentum within the shell will be in the $`\widehat{\varphi }`$-direction, so for the remainder of the discussion we will only need the magnitude of the angular momentum flux $`\mathrm{\Phi }_J=𝚽_J\widehat{\varphi }`$.

The mass, momentum, and angular momentum flux functions due to the incident winds are similarly defined. For the first wind, located at the coordinate origin, these are given in terms of the wind density $`\rho _w`$ and velocity $`𝐕_w`$ as

$$\mathrm{\Phi }_{m,w}\mathrm{\Delta }\varphi =\rho _w(𝐕_w\widehat{𝐧})𝑑A,$$
$`(7)`$
$$𝚽_w\mathrm{\Delta }\varphi =\rho _w(𝐕_w\widehat{𝐧})𝐕_w𝑑A,$$
$`(8)`$
$$𝚽_{J,w}\mathrm{\Delta }\varphi =\rho _w(𝐕_w\widehat{𝐧})(𝐑\times 𝐕_w)𝑑A.$$
$`(9)`$
The area of integration is that between two planes of constant $`\varphi `$, from the standoff point to polar angle $`\theta `$, following the shape of the collision surface, which is determined along the course of the integrations. The sign of the normal direction is chosen so as to point away from the origin. The flux functions for the ambient medium are strictly analogous, except the unit vector normal to the collision surface would be in the reverse direction:

$$\mathrm{\Phi }_{m,a}\mathrm{\Delta }\varphi =\rho _a(𝐕_a\widehat{𝐧})𝑑A,$$
$`(10)`$
$$𝚽_a\mathrm{\Delta }\varphi =\rho _a(𝐕_a\widehat{𝐧})𝐕_a𝑑A,$$
$`(11)`$
$$𝚽_{J,a}\mathrm{\Delta }\varphi =\rho _a(𝐕_a\widehat{𝐧})(𝐑\times 𝐕_a)𝑑A.$$
$`(12)`$
If $`R(\theta ,\varphi )`$ is the radius of the shell, and its partial derivatives are $`R_\theta `$ and $`R_\varphi `$, the surface area element is

$$dA=\sqrt{R^2+R_\theta ^2+R_\varphi ^2\mathrm{csc}^2\theta }R\mathrm{sin}\theta d\theta d\varphi .$$
$`(13)`$
The unit outward normal $`\widehat{𝐧}`$ to the surface is given by

$$\widehat{𝐧}=\frac{\left(R\widehat{𝐫}R_\theta \widehat{\theta }R_\varphi \mathrm{csc}\theta \widehat{\varphi }\right)}{\sqrt{R^2+R_\theta ^2+R_\varphi ^2\mathrm{csc}^2\theta }}.$$
$`(14)`$
Combining these,

$$\widehat{𝐧}dA=\left(R\widehat{𝐫}R_\theta \widehat{\theta }R_\varphi \mathrm{csc}\theta \widehat{\varphi }\right)R\mathrm{sin}\theta d\theta d\varphi .$$
$`(15)`$
For a radial wind, the partial derivatives with respect to $`\theta `$ and $`\varphi `$ do not enter the expression for $`𝐕_w\widehat{𝐧}dA`$, which is simply

$$𝐕_w\widehat{𝐧}dA=V_wR^2\mathrm{sin}\theta d\theta d\varphi .$$
$`(16)`$
Similarly, for the ambient medium $`𝐕_a=V_a\widehat{𝐳}=V_a(\widehat{𝐫}\mathrm{cos}\theta \widehat{\theta }\mathrm{sin}\theta )`$, so noting that $`d\varpi /d\theta =R\mathrm{cos}\theta +R_\theta \mathrm{sin}\theta `$, we have

$$𝐕_a\widehat{𝐧}dA=V_a\varpi d\varpi d\varphi .$$
$`(17)`$
This eliminates some of the partial derivatives that complicated Bandiera’s very general treatment of the geometry. The rate at which conserved quantities are advected onto the shell does not depend upon the detailed description of the shell, provided we ensure the correct number of streamlines is counted. Using eq.(16), the resulting forms of eqs.(7-9) for the radial stellar wind are

$$\mathrm{\Phi }_{m,w}=_0^\theta R^2\rho _wV_w\mathrm{sin}\theta ^{}d\theta ^{},$$
$`(18)`$
$$𝚽_w=_0^\theta R^2\rho _wV_w^2[\widehat{\varpi }\mathrm{sin}\theta ^{}+\widehat{𝐳}\mathrm{cos}\theta ^{}]\mathrm{sin}\theta ^{}d\theta ^{},$$
$`(19)`$
$$\mathrm{\Phi }_{J,w}=0.$$
$`(20)`$
Here the radial unit vector has been written in terms of its cylindrical polar components, and the integrations are performed at constant $`\varphi `$. Because the stellar wind is radial, it imparts no angular momentum to the shell, and the wind angular momentum flux vanishes. Using eq.(17) to simplify eqs.(10-12) for the ambient medium,

$$\mathrm{\Phi }_{m,a}=V_a_0^\varpi \rho _a\varpi ^{}𝑑\varpi ^{},$$
$`(21)`$
$$𝚽_a=V_a^2\widehat{𝐳}_0^\varpi \rho _a\varpi ^{}𝑑\varpi ^{},$$
$`(22)`$
$$𝚽_{J,a}=V_a^2\widehat{\varphi }_0^\varpi \rho _a\varpi _{}^{}{}_{}{}^{2}𝑑\varpi ^{}.$$
$`(23)`$
The mass conservation law for the pre-shock radial wind in steady-state is

$$\frac{}{r}\left(r^2\rho _wV_w\right)=0.$$
$`(24)`$
Consequently, if the wind is of mass loss rate $`\dot{M}_w`$ and streamline-average speed $`\overline{V}_w`$, its properties are given by

$$\rho _wV_w=\frac{\dot{M}_w}{4\pi r^2}f_w(\theta ,\varphi ),$$
$`(25)`$
$$\rho _wV_w^2=\frac{\dot{M}_w\overline{V}_w}{4\pi r^2}g_w(\theta ,\varphi ).$$
$`(26)`$
Here the dimensionless functions $`f_w`$ and $`g_w`$ are normalized to have unit streamline-average over $`4\pi `$ steradians. Note that eq.(26) assumes the stellar wind to be coasting. Given this assumption, $`g_w`$ is independent of the radius. The wind mass and momentum fluxes onto the shell may now be written

$$\mathrm{\Phi }_{m,w}=\frac{\dot{M}_w}{4\pi }_0^\theta f_w(\theta ^{},\varphi )\mathrm{sin}\theta ^{}d\theta ^{},$$
$`(27)`$
$$𝚽_w=\frac{\dot{M}_w\overline{V}_w}{4\pi }_0^\theta g_w(\theta ^{},\varphi )[\widehat{\varpi }\mathrm{sin}\theta ^{}+\widehat{𝐳}\mathrm{cos}\theta ^{}]\mathrm{sin}\theta ^{}d\theta ^{}.$$
$`(28)`$
In order to obtain steady-state solutions, the ambient density must be independent of the z-coordinate, although its dependence on the remaining $`\varpi ,\varphi `$ coordinates may be arbitrary. Thus, we assume an ambient density of the form

$$\rho _a=\rho _0f_a(\varpi ,\varphi ),$$
$`(29)`$
where $`\rho _0`$ is the value of the ambient density along the stellar trajectory and $`f(0,\varphi )=1`$. The final forms of the ambient flux functions are thus

$$\mathrm{\Phi }_{m,a}=\rho _0V_a_0^\varpi f_a\varpi ^{}𝑑\varpi ^{}.$$
$`(30)`$
$$𝚽_a=V_a\widehat{𝐳}\mathrm{\Phi }_{m,a}.$$
$`(31)`$
$$\mathrm{\Phi }_{J,a}=\rho _0V_a^2_0^\varpi f_a\varpi ^2𝑑\varpi ^{}.$$
$`(32)`$
The mass, momentum, and angular momentum flux functions for both the wind and ambient medium are clearly streamline integrals that do not depend upon the detailed shape of the shell - the integration is essentially one over the solid angle of the $`\mathrm{\Delta }\varphi `$ wedge, as seen by the origin. This is a consequence of the fact that the pre-shock media conserve these quantities in steady-state.

Suppose we have performed the integrations for the mass, momentum, and angular momentum onto the narrow slice of the shell, for both incident winds, according to eqs.(20,27,28,30-32). In steady-state, these quantities do not accumulate at any location within the shell, but are carried away by the flow within the layer. The conservation laws of mass, momentum, and angular momentum take the form

$$\mathrm{\Phi }_m=\mathrm{\Phi }_{m,w}+\mathrm{\Phi }_{m,a},$$
$`(33)`$
$$𝚽=𝚽_w+𝚽_a,$$
$`(34)`$
$$\mathrm{\Phi }_J=\mathrm{\Phi }_{J,w}+\mathrm{\Phi }_{J,a}.$$
$`(35)`$
Now we may describe the properties of the shell explicitly in terms of those of the external winds. Equation (6) yields the shell radius

$$R=\frac{\mathrm{\Phi }_J}{\mathrm{\Phi }_\varpi \mathrm{cos}\theta \mathrm{\Phi }_z\mathrm{sin}\theta }.$$
$`(36)`$
This equation combines the formulation of Paper II, which was in terms of azimuthally integrated fluxes, with the original treatment of Paper I, using fluxes in an infinitesimally thin wedge. In eq.(36), each momentum flux function within the shell is to be evaluated as the sum of appropriate source terms, using the conservation laws (eqs.). It is to be stressed that if the flow is not meridional, eq.(5) does not hold, and the solution method is considerably more complicated. However, eqs.(6,36) are still valid provided the angular momentum flux is replaced by its $`\widehat{\varphi }`$ component. From eqs.(1,2), the tangential velocity of matter in the shell is

$$V_t\widehat{𝐭}=𝚽/\mathrm{\Phi }_m,$$
$`(37)`$
while the mass surface density in the shell is given by

$$\sigma =(\mathrm{\Phi }_m^2/\varpi |𝚽|)\mathrm{cos}\eta .$$
$`(38)`$
### 2.4 The Applied Torque Method and an Example

To provide a specific example of the solution method, as well as a reference solution to be compared to in the following, we consider the simple problem of a bow shock from an isotropic wind in a uniform ambient medium of density $`\rho _a`$. In this case, $`f_w=g_w=f_a=1`$, and the flux functions are

$$\mathrm{\Phi }_{m,w}=\frac{\dot{M}_w}{4\pi }(1\mathrm{cos}\theta ),$$
$`(39)`$
$$𝚽_w=\frac{\dot{M}_wV_w}{8\pi }[\widehat{\varpi }(\theta \mathrm{sin}\theta \mathrm{cos}\theta )+\widehat{𝐳}\mathrm{sin}^2\theta ],$$
$`(40)`$
$$\mathrm{\Phi }_{J,w}=0,$$
$`(20)`$
$$\mathrm{\Phi }_{m,a}=\frac{1}{2}\varpi ^2\rho _aV_a,$$
$`(41)`$
$$𝚽_a=\frac{1}{2}\varpi ^2\rho _aV_a^2\widehat{𝐳},$$
$`(42),`$
$$\mathrm{\Phi }_{J,a}=\frac{1}{3}\varpi ^3\rho _aV_a^2.$$
$`(43)`$
Because the stellar wind has uniform speed, we have replaced $`\overline{V}_w`$ with $`V_w`$. One could immediately substitute these results into eq.(36) to determine the shell’s shape. However, before proceeding to apply the formalism, note that the equation for the shape of the shell depends upon a specific combination of three flux functions in eq.(36). One may divide the equation into two parts, one part for each wind, according to

$$𝒯_k=\varpi [\mathrm{\Phi }_{\varpi ,k}\mathrm{cot}\theta \mathrm{\Phi }_{z,k}]\mathrm{\Phi }_{\mathrm{J},k},$$
$`(44)`$
and $`k=w`$ or $`a`$ for the wind and ambient (or second wind) sources. Physically, $`𝒯`$ represents the applied torque necessary to compress the fan of streamlines between the stagnation point and position $`R(\theta ,\varphi )`$ into a unidirectional beam possessing the same momentum flux (we shall call $`𝒯`$ the required torque). This is clearly seen by considering the radial wind from the coordinate origin. It has no angular momentum, but an equivalent beam containing the same momentum flux, located at a specific point $`R(\theta ,\varphi )`$ does indeed have angular momentum about the origin, because the beam will not be radial. The required torque for the isotropic wind follows from eqs. (20,40,44):

$$𝒯_w=\frac{\dot{M}_wV_w\varpi }{8\pi }(1\theta \mathrm{cot}\theta ).$$
$`(45)`$
Using eqs.(41-44), the torque necessary to compress the ambient streamlines to a unidirectional beam is

$$𝒯_a=\frac{1}{6}\varpi ^3\rho _aV_a^2.$$
$`(46)`$
For the ambient medium, the required torque expression simplifies because the ambient medium has no cylindrical radial momentum, and its z-momentum flux is related to its mass flux by eq.(31), so in general,

$$𝒯_a=\varpi V_a\mathrm{\Phi }_{m,a}\mathrm{\Phi }_{J,a}.$$
$`(47)`$
Similarly, the radial wind’s required torque is simplified due to its lack of angular momentum about the origin.

The solution surface is now defined by

$$𝒯_w+𝒯_a=0.$$
$`(48)`$
The interpretation of eq.(48) is that each wind supplies the torque necessary to compress the other wind. These torques are equal and opposite only if the correct value of the radius (or cylindrical radius $`\varpi `$) is used, conserving angular momentum as well as linear momentum. The utility of this approach is that if we vary the properties of one of the winds, holding the other constant, the shape of the shell is affected only by the change in the applied torque from the second wind.

Now the applied torque expressions for the wind (eq.) and the ambient medium (eq.), when substituted into eq.(48), quickly recover the solution of Paper I,

$$R(\theta )=R_0\mathrm{csc}\theta \sqrt{3(1\theta \mathrm{cot}\theta )},$$
$`(49)`$
where $`R_0`$ is the standoff radius:

$$R_0=\sqrt{\frac{\dot{M}_wV_w}{4\pi \rho _aV_a^2}}.$$
$`(50)`$
## 3 Solutions for Bow Shock with Ambient Density Gradient

In this section we assume the asymmetry of the bow shock to be due to a density gradient in the ambient medium, while the the stellar wind is assumed to be isotropic. The next section will treat bow shocks from anisotropic winds.

The mass and momentum contributions of the stellar wind to the shell are given by eqs.(39,40,20), and its required torque is given by eq.(45). The ambient medium is assumed to have a plane-parallel density stratification, so in eq.(29), we replace $`f_a(\varpi ,\varphi )`$ simply by $`f_a(x)`$, with $`f_a`$ is normalized so that $`f_a(0)=1`$, and where

$$x=\varpi \mathrm{cos}\varphi .$$
$`(51)`$
The integrated mass and angular momentum flux onto the wedge from the ambient medium are

$$\mathrm{\Phi }_{m,a}=\rho _0V_a_0^\varpi f_a(\varpi ^{}\mathrm{cos}\varphi )\varpi ^{}𝑑\varpi ^{}.$$
$`(52)`$
$$\mathrm{\Phi }_{J,a}=\rho _0V_a^2_0^\varpi f_a(\varpi ^{}\mathrm{cos}\varphi )\varpi ^2𝑑\varpi ^{}.$$
$`(53)`$
A complete description of the bow shock is now obtained by adding the contributions of both the wind and ambient medium using the conservation laws (eqs.\[33-35\]).

### 3.1 Solution for Exponential Density Law

We first obtain the solution for an exponential distribution of ambient density. Let the density scale height be $`H`$. Define for brevity the coordinate $`y=x/H`$, so the density law is

$$f_a(x)=\mathrm{exp}(y),$$
$`(54)`$
and the mass and angular momentum flux integrals are

$$\mathrm{\Phi }_{m,a}=H^2\rho _0V_a\mathrm{sec}^2\varphi [1\mathrm{exp}(y)(1+y)],$$
$`(55)`$
$$\mathrm{\Phi }_{J,a}=H^3\rho _0V_a^2\mathrm{sec}^3\varphi [2\mathrm{exp}(y)(y^2+2y+2)].$$
$`(56)`$
The ambient momentum flux follows from eq.(31), while the required torque follows from eq.(47), which yields

$$𝒯_a=H^3\rho _0V_a^2\mathrm{sec}^3\varphi \left[y2+\mathrm{exp}(y)(y+2)\right].$$
$`(57)`$
Using this torque formula and that for the stellar wind given by eq.(45), the shell’s shape is given by eq.(48), which yields upon simplification

$$l^2\mathrm{sec}^2\varphi y^1[y2+\mathrm{exp}(y)(y+2)]=\frac{1}{2}(1\theta \mathrm{cot}\theta ).$$
$`(58)`$
Here $`l=H/R_0`$ is the density scale height in units of the standoff radius. Letting $`y_{ax}`$ be the value of $`y=x/H`$ appropriate for the axisymmetric solution of eq.(49),

$$y_{ax}=\frac{R_0}{H}\mathrm{cos}\varphi \sqrt{3(1\theta \mathrm{cot}\theta )},$$
$`(59)`$
eq.(58) takes the simple form

$$\frac{1}{y}\left[y2+\mathrm{exp}(y)(y+2)\right]=\frac{1}{6}y_{ax}^2.$$
$`(60)`$
This formula may be solved numerically for $`y(\theta ,\varphi )`$, so letting $`y_s`$ denote the solution for $`y`$ to eq.(60), the shell’s shape is given by

$$R(\theta ,\varphi )=H\mathrm{sec}\varphi \mathrm{csc}\theta y_s(\theta ,\varphi ).$$
$`(61)`$
The result is a family of bow shock solutions distinguished by the value of the nondimensional parameter $`l`$. By examining eq.(60) we may deduce several properties of the solution. The left-hand side of the equation is monotonically increasing as a function of $`y`$. The signs of $`y`$ and $`y_{ax}`$ must be the same, since it is due to the $`\mathrm{cos}\varphi `$ factor, so it follows that $`y`$ increases monotonically, although in a very non-linear fashion, as a function of $`y_{ax}`$. Because the parameter $`l`$ does not enter eq.(60), we see that there is one universal solution for the problem, although this solution is scaled by the value of $`R_0`$ and stretched (distorted) depending upon the value of $`H/R_0`$. The solution for the radius at all angles $`\theta ,\varphi `$ and for all values of $`l`$ follows from Fig. Modeling Non-Axisymmetric Bow Shocks: Solution Method and Exact Analytic Solutions. In particular, as $`y\mathrm{}`$, the left hand side of eq.(60) approaches unity, implying a breakout angle for the shell when $`y_{ax}=\sqrt{6}`$, which yields

$$1\theta _{\mathrm{}}\mathrm{cot}\theta _{\mathrm{}}=2l^2\mathrm{sec}^2\varphi .$$
$`(62)`$
The opening angle $`\theta _{\mathrm{}}`$ depends upon $`\mathrm{cos}\varphi `$, and applies only to $`\mathrm{cos}\varphi 0`$. In the tail of the bow shock, for $`\mathrm{cos}\varphi <0`$, we have $`y\mathrm{}`$. In this case, $`y`$ depends logarithmically upon $`y_{ax}`$. The most interesting part of the solution, for $`|y|1`$ may be described by an expansion in terms of small $`y_{ax}`$. Noting that for $`\mathrm{cos}\varphi =0`$ the solution is identical to that of eq.(49), we anticipate that an expansion for small $`\mathrm{cos}\varphi `$ corresponds to an expansion in terms of small $`y_{ax}`$ about the standard bow shock solution. Letting $`R_{ax}`$ be the value of $`R(\theta )`$ for that solution (given by eq.), we obtain

$$RR_{ax}[1+\frac{1}{4}y_{ax}+\frac{13}{160}y_{ax}^2+\frac{7}{240}y_{ax}^3+\frac{11843}{1075000}y_{ax}^4+\mathrm{}.].$$
$`(63)`$
The behavior of the solution near the stagnation point is given by

$$R=R_0\left[1+\frac{\mathrm{cos}\varphi }{4l}\theta +\left(\frac{1}{5}+\frac{13\mathrm{cos}^2\varphi }{160l^2}\right)\theta ^2+\mathrm{}\right].$$
$`(64)`$
Unlike the axisymmetric bow shock, there is a linear term in the behavior near the standoff point, so it is not describable as parabolic with the z coordinate axis as the axis of a parabola, except for the special angles $`\varphi =\pi /2\mathrm{or}3\pi /2`$. As a consistency check, note that as the scale height increases relative to the standoff radius $`(l=H/R_0\mathrm{})`$, the solution reduces to the standard axisymmetric bow shock (Paper I). The lowest order effect near the standoff point is a tilt of the bow shock head, described by the linear term in $`\theta `$. This means that although the wind and ambient streamlines meet head-on at the standoff point, they are not normal to the shell at this location, quite different from the axisymmetric case. This effect is due to the instantaneous mixing we have assumed. If shear is present in the shell, the two colliding flows could have separate stagnation points, which would be located where the incident stream is indeed normal to the shell. These bow shock solutions differ from the standard model (Baranov, Krasnobaev & Kulikovskii 1971) in that the exponential ambient density distribution implies a finite total mass on one side of the bow shock. The solution then more closely resembles those of two-wind collisions (Paper II) where there is a finite opening angle for the bow shock tail. Similar problems include blast waves in finite mass media (Koo & McKee 1990) and wind breakout from a stratified medium (Cantó 1980; Borkowski, Blondin & Harrington 1997). For the high density region, $`\mathrm{cos}\varphi <0`$, the shell’s asymptotic shape in the tail region is

$$z\frac{2l^2}{\pi }\varpi \mathrm{sec}^2\varphi \mathrm{exp}\left(\frac{\varpi }{H}\mathrm{cos}\varphi \right).$$
$`(65)`$
The shell continues to be more distorted as it expands due to the $`\mathrm{exp}(y)`$ factor. Because the shell has a finite opening angle on one side, not all wind streamlines intersect the shell, but those with $`\theta >\theta _{\mathrm{}}(\varphi )`$ freely expand to infinity. Further consequences of the finite opening angle of the shell are discussed in §5.

### 3.2 Solution for Linear or Polynomial Density Law

Consider a stratification of the ambient medium according to

$$f_a(x)=1+a_1x+a_2x^2+\mathrm{}.$$
$`(66)`$
Of course, one must ensure that this expression is non-negative over the domain of interest. The fluxes of mass and angular momentum onto the shell are

$$\mathrm{\Phi }_{m,a}=\varpi ^2\rho _0V_a\left[\frac{1}{2}+\frac{a_1}{3}x+\frac{a_2}{4}x^2+\mathrm{}\right],$$
$`(67)`$
$$\mathrm{\Phi }_{J,a}=\varpi ^3\rho _0V_a^2\left[\frac{1}{3}+\frac{a_1}{4}x+\frac{a_2}{5}x^2+\mathrm{}\right].$$
$`(68)`$
The ambient momentum flux follows from eq.(31), while the required torque follows from eq.(47), giving the result

$$𝒯_a=\frac{\varpi ^3}{6}\rho _0V_a^2\left[1+\frac{a_1}{2}x+\frac{3}{10}a_2x^2+\mathrm{}+\frac{6(n+1)!}{(n+3)!}a_nx^n+\mathrm{}\right].$$
$`(69)`$
Substitution into eq.(48), we obtain upon simplification

$$\frac{\varpi ^2}{R_0^2}\left[1+\frac{a_1}{2}\varpi \mathrm{cos}\varphi +\frac{3a_2}{10}\varpi ^2\mathrm{cos}^2\varphi +\mathrm{}\right]=3(1\theta \mathrm{cot}\theta ).$$
$`(70)`$
Restricting the treatment to a linear density gradient, with only $`a_1`$ non-zero, this equation is cubic in the variable $`\varpi `$. One may solve analytically or numerically for the function $`\varpi (\theta ,\varphi )`$, so the solution surface is then given by $`R(\theta ,\varphi )=\varpi (\theta ,\varphi )\mathrm{csc}\theta `$. Inclusion of higher order terms in the density law simply increases the degree of the polynomial, but the solution is obtained in the same way.

The behavior of the solution near the stagnation point is given by

$$\frac{R}{R_0}=1\frac{a_1R_0\mathrm{cos}\varphi }{4}\theta +\left(\frac{1}{5}+(\frac{5}{32}a_1^2\frac{3}{20}a_2)R_0^2\mathrm{cos}^2\varphi \right)\theta ^2+\mathrm{}.$$
$`(71)`$
As a consistency check of the solutions, for the case of a uniform ambient medium $`(a_1,a_2,\mathrm{}=0)`$, eqs.(70,71) reduce to the standard solution of Paper I. As a further check, note that this result is consistent with the solution for an exponential mass distribution, if we choose the coefficients corresponding to the exponential distribution of the previous subsection, $`a_1=1/H,a_2=1/2H^2,\mathrm{}`$, recovering eq.(64).

## 4 Treatment of Anisotropic Winds

### 4.1 The Axisymmetric Wind

We now relax the assumption that the wind is isotropic to permit an axisymmetric wind, where the axis of symmetry is misaligned with the direction of stellar motion, or in the case of a two-wind collision in a binary star, where the axis does not point to the other star. Treatment of the accelerated wind will be deferred to a future contribution; in this paper the wind is assumed to be coasting. The coordinate axes are chosen so that the z-direction points in the direction of stellar motion, or towards the second star. In terms of a spherical coordinate system, where the azimuthal angle $`\varphi `$ is measured about the z-axis, motion of the shocked fluid in the shell is along planes of constant $`\varphi `$. We also define starred coordinates so that the $`z_{}`$-axis is the symmetry axis of the stellar wind. The two coordinate systems are related by

$`\mathrm{sin}\theta _{}\mathrm{cos}\varphi _{}`$ $`=`$ $`\mathrm{sin}\theta \mathrm{cos}\varphi ,`$ (72)
$`\mathrm{sin}\theta _{}\mathrm{sin}\varphi _{}`$ $`=`$ $`\mathrm{sin}\theta \mathrm{sin}\varphi \mathrm{cos}\lambda \mathrm{cos}\theta \mathrm{sin}\lambda ,`$ (73)
$`\mathrm{cos}\theta _{}`$ $`=`$ $`\mathrm{sin}\theta \mathrm{sin}\varphi \mathrm{sin}\lambda +\mathrm{cos}\theta \mathrm{cos}\lambda .`$ (74)
The wind mass and momentum flux densities depend on the polar angle $`\theta _{}`$:

$$\rho _wV_w=\frac{\dot{M}_w}{4\pi r^2}f_w(\theta _{}),$$
$`(75)`$
$$\rho _wV_w^2=\frac{\dot{M}_w\overline{V}_w}{4\pi r^2}g_w(\theta _{}).$$
$`(76)`$
The nondimensional functions $`f_w`$ and $`g_w`$ are normalized to have unit average value over $`4\pi `$ steradians.

The incident fluxes of mass, momentum, and angular momentum from the wind onto the shell are

$`\mathrm{\Phi }_{m,w}(\theta ,\varphi )`$ $`=`$ $`{\displaystyle \frac{\dot{M}_w}{4\pi }}F_w(\theta ,\varphi ),`$ (77)
$`𝚽_w(\theta ,\varphi )`$ $`=`$ $`{\displaystyle \frac{\dot{M}_w\overline{V}_w}{4\pi }}𝐆_w(\theta ,\varphi ),`$ (78)
$`\mathrm{\Phi }_{\mathrm{J},w}(\theta ,\varphi )`$ $`=`$ $`0.`$ (20)
where the nondimensional functions $`F_w`$ and $`𝐆_w=G_{w,\varpi }\widehat{\varpi }+G_{w,z}\widehat{𝐳}`$ are given by

$`F_w`$ $`=`$ $`{\displaystyle _0^\theta }f_w\mathrm{sin}\theta ^{}\mathrm{d}\theta ^{},`$ (79)
$`𝐆_w`$ $`=`$ $`{\displaystyle _0^\theta }g_w[\widehat{\varpi }\mathrm{sin}\theta ^{}+\widehat{𝐳}\mathrm{cos}\theta ^{}]\mathrm{sin}\theta ^{}\mathrm{d}\theta ^{}.`$ (80)
Here $`f_w`$ and $`g_w`$ have argument $`\theta _{}(\theta ^{},\varphi )`$, so the integrations are performed using the transformation eqs.(72-74) to evaluate $`\theta _{}(\theta ,\varphi )`$. We also define a nondimensional function $`T_w`$ associated with the required torque $`𝒯_w`$ to compress the wind streamlines according to

$$𝒯_w=\frac{\dot{M}_w\overline{V}_w}{4\pi }\varpi T_w.$$
$`(81)`$
Because the wind’s angular momentum flux vanishes, eq.(44) implies that

$$T_w=G_{w,\varpi }\mathrm{cot}\theta G_{w,z}.$$
$`(82)`$
The functions $`f_w(\theta _{}),g_w(\theta _{})`$ may now be expanded in terms of powers of $`\mathrm{cos}\theta _{}`$,

$$f_w(\theta _{})=\underset{i=0}{\overset{\mathrm{}}{}}b_i\mathrm{cos}^i\theta _{},$$
$`(83)`$
$$g_w(\theta _{})=\underset{i=0}{\overset{\mathrm{}}{}}c_i\mathrm{cos}^i\theta _{}.$$
$`(84)`$
Defining for brevity

$$p=\mathrm{sin}\varphi \mathrm{sin}\lambda ,$$
$`(85)`$
$$q=\mathrm{cos}\lambda ,$$
$`(86)`$
the transformation given by eq.(74) is $`\mathrm{cos}\theta _{}=p\mathrm{sin}\theta +q\mathrm{cos}\theta `$. Denoting the trigonometric integrals by

$$\mathrm{I}_{j,k}=_0^\theta \mathrm{sin}^j\theta ^{}\mathrm{cos}^k\theta ^{}d\theta ^{},$$
$`(87)`$
we evaluate $`F_w`$ and $`𝐆_w`$ using the above series, and bringing $`p`$ and $`q`$ outside of the integrals, we finally have

$$F_w=\underset{i=0}{\overset{\mathrm{}}{}}b_i\underset{j=0}{\overset{i}{}}p^{ij}q^j\left(\genfrac{}{}{0pt}{}{i}{j}\right)I_{1+ij,j},$$
$`(88)`$
$$𝐆_w=\underset{i=0}{\overset{\mathrm{}}{}}c_i\underset{j=0}{\overset{i}{}}p^{ij}q^j\left(\genfrac{}{}{0pt}{}{i}{j}\right)\left[\widehat{\varpi }I_{2+ij,j}+\widehat{𝐳}I_{1+ij,1+j}\right].$$
$`(89)`$
Here we have used the binomial coefficients to write the sums. The applied torque necessary to compress the wind streamlines to a thin shell, in nondimensional form, is now

$$T_w=\underset{i=0}{\overset{\mathrm{}}{}}c_iT_w^{(i)},$$
$`(90)`$
where the responses due to the individual $`cos^i\theta _{}`$ terms are given by

$$T_w^{(i)}=\underset{j=0}{\overset{i}{}}p^{ij}q^j\left(\genfrac{}{}{0pt}{}{i}{j}\right)\left[\mathrm{cot}\theta \mathrm{I}_{2+ij,j}\mathrm{I}_{1+ij,1+j}\right].$$
$`(91)`$
### 4.2 General Solution for Bow Shock Driven by an Axisymmetric, Misaligned Wind

The wind description of the previous subsection may now be applied to the problem of a bow shock driven by an anisotropic wind. The wind is assumed to be driven by a star moving at speed $`V_a`$ in a medium of uniform density $`\rho _a`$. We describe the bow shock’s properties in the frame of the star.

The standoff radius, defined as the shell radius at $`\theta =0,`$ which corresponds to $`\theta _{}=\lambda `$, is given by

$$R_\lambda =\sqrt{\frac{\dot{M}_w\overline{V}_wg_w(\lambda )}{4\pi \rho _aV_a^2}}.$$
$`(92)`$
We also define $`R_0`$ to be the standoff radius for the equivalent isotropic wind

$$R_0=\sqrt{\frac{\dot{M}_w\overline{V}_w}{4\pi \rho _aV_a^2}}.$$
$`(93)`$
It is important to recall for the results below that the standoff radius will be $`R_0\sqrt{g_w(\lambda )}`$, where

$$g_w(\lambda )=\underset{i=0}{\overset{\mathrm{}}{}}c_i\mathrm{cos}^i\lambda .$$
$`(94)`$
By eqs.(46,48,82), the solution for the shell’s radius is

$$R(\theta ,\varphi )=R_0\mathrm{csc}\theta \sqrt{6T_w},$$
$`(95)`$
where $`T_w`$ is given by eqs.(90,91). The total mass and momentum fluxes, including the contribution from the ambient medium, are

$$\mathrm{\Phi }_m=\frac{\dot{M}_w}{4\pi }\left[F_w+\frac{1}{2\alpha }\stackrel{~}{\varpi }^2\right],$$
$`(96)`$
$$𝚽=\frac{\dot{M}_w\overline{V}_w}{4\pi }\left[G_{w,\varpi }\widehat{\varpi }+(G_{w,z}\frac{1}{2}\stackrel{~}{\varpi }^2)\widehat{𝐳}\right],$$
$`(97)`$
where $`\stackrel{~}{\varpi }=\varpi /R_0`$, and $`\alpha =V_a/\overline{V}_w`$.

To obtain a complete description of the bow shock’s properties, we must specify the velocity of material in the shell and the mass surface density of matter. The velocity is given by the ratio of the mass and momentum fluxes:

$$V_t\widehat{𝐭}=V_a\frac{\left[2G_{w,\varpi }\widehat{\varpi }+(2G_{w,z}\stackrel{~}{\varpi }^2)\widehat{𝐳}\right]}{\left[2\alpha F_w+\stackrel{~}{\varpi }^2\right]}.$$
$`(98)`$
The surface density is given by eq.(38), which yields

$$\sigma =\rho _aR_0\frac{\left[2\alpha F_w+\stackrel{~}{\varpi }^2\right]^2}{2\stackrel{~}{\varpi }\sqrt{\left[4G_{w,\varpi }^2+(2G_{w,z}\stackrel{~}{\varpi }^2)^2\right]}}\mathrm{cos}\eta .$$
$`(99)`$
In order to evaluate $`\sigma `$, we need to know the angle $`\eta `$. Using eqs.(4,95), we have

$$\mathrm{tan}\eta =\frac{1}{2}\mathrm{csc}\theta \frac{1}{T_w}\frac{T_w}{\varphi },$$
$`(100)`$
where differentiation of eqs.(90,91) with respect to $`\varphi `$, which enters only in the variable $`p`$, yields

$$\frac{T_w}{\varphi }=\underset{i=1}{\overset{\mathrm{}}{}}c_i\frac{T_w^{(i)}}{\varphi },$$
$`(101)`$
$$\frac{T_w^{(i)}}{\varphi }=\mathrm{cot}\varphi \underset{j=0}{\overset{i1}{}}(ij)p^{ij}q^j\left(\genfrac{}{}{0pt}{}{i}{j}\right)\left[\mathrm{cot}\theta \mathrm{I}_{2+ij,j}\mathrm{I}_{1+ij,1+j}\right].$$
$`(102)`$
Note that the summation limits are restricted so that the terms with $`j=i`$ vanish, including the $`c_0`$ term.

In the region of the standoff point, an expansion to second order in $`\theta `$ gives

$$RR_\lambda \{1+\theta \frac{pg_w^{}(\lambda )}{4g_w(\lambda )}+\theta ^2[\frac{1}{5}\frac{p^2g_w^2(\lambda )}{32g_w^2(\lambda )}$$

$$+\frac{3(p^2g_w^{\prime \prime }(\lambda )qg_w^{}(\lambda ))}{40g_w(\lambda )}]\}.$$
$`(103)`$
Primes on the function $`g_w`$ represent differentiation with respect to $`q=\mathrm{cos}\lambda `$. As was the case of asymmetry resulting from an ambient density gradient, there is a $`\theta ^1`$ term indicating a tilt of the standoff region relative to the direction of the stellar motion. The mass and momentum flux functions near the standoff point are

$$\mathrm{\Phi }_m\frac{\dot{M}_w}{4\pi }\left\{\frac{1}{2}f_w(\lambda )+\frac{1}{2\alpha }g_w(\lambda )\right\}\theta ^2.$$
$`(104)`$
$$𝚽\frac{\dot{M}_w\overline{V}_w}{4\pi }\left\{\frac{1}{3}g_w(\lambda )\widehat{\varpi }+\frac{p}{3}g_w^{}(\lambda )\widehat{𝐳}\right\}\theta ^3.$$
$`(105)`$
Equations (4,28,29) now give the tangential velocity in the shell and the mass per unit area as

$$V_t\widehat{𝐭}\frac{2}{3}V_a\theta \frac{g_w(\lambda )\widehat{\varpi }+pg_w^{}(\lambda )\widehat{𝐳}}{\alpha f_w(\lambda )+g_w(\lambda )},$$
$`(106)`$
and

$$\sigma \frac{3R_0\rho _a}{4}\frac{(\alpha f_w(\lambda )+g_w(\lambda ))^2}{\sqrt{g_w(\lambda )(g_w^2(\lambda )+p^2(g_w^{}(\lambda ))^2}}.$$
$`(107)`$
One may readily see that these solutions reduce to the results of Paper I for an isotropic wind ($`f_w=g_w=1`$).

### 4.3 Simple Solutions: Quadratic Dependence on $`cos\theta _{}`$

The simplest non-trivial solution is for a linear dependence of mass and momentum fluxes upon $`\mathrm{cos}\theta _{}`$. However, in astrophysical applications one is frequently concerned with stellar winds that have symmetry with respect to the equatorial plane, such as when the asymmetry arises due to rotation of the star. For winds that are symmetric with respect to $`\theta =\pi /2`$, the functions $`f_w(\theta _{}),g_w(\theta _{})`$ will have expansions in terms of even powers of $`\mathrm{cos}\theta _{}`$ only. Thus, we shall give the solution for a wind that depends upon $`\mathrm{cos}^2\theta _{}`$. This is sufficiently general to include as a subset previous models of non-axisymmetric bow shocks (Bandiera 1993; Chen & Huang 1997). We assume the mass and momentum fluxes are described by

$$f_w=b_0+b_1\mathrm{cos}\theta _{}+b_2\mathrm{cos}^2\theta _{},$$
$`(108)`$
$$g_w=c_0+c_1\mathrm{cos}\theta _{}+c_2\mathrm{cos}^2\theta _{},$$
$`(109)`$
where the normalization requires $`b_0=(1b_2/3)`$ and $`c_0=(1c_2/3)`$.

The mass flux integral $`F_w`$ is

$$F_w=\{b_0^{}(1\mu )+\frac{b_1}{2}[p(\theta \mathrm{sin}\theta \mathrm{cos}\theta )+q\mathrm{sin}^2\theta ]$$

$$+\frac{b_2}{3}\mathrm{sin}^2\theta [(q^2p^2)\mathrm{cos}\theta +2pq\mathrm{sin}\theta ]\},$$
$`(110)`$
where $`b_0^{}=b_0+b_2(2p^2+q^2)/3`$. The components of $`𝐆_w`$ are

$$G_{w,\varpi }=\frac{1}{4}\{c_0^{}(2\theta \mathrm{sin}2\theta )+\frac{c_1}{3}[p(89\mathrm{cos}\theta +\mathrm{cos}3\theta )+4q\mathrm{sin}^3\theta ]$$

$$+c_2\mathrm{sin}^3\theta [2pq\mathrm{sin}\theta +(q^2p^2)\mathrm{cos}\theta ]\},$$
$`(111)`$
$$G_{w,z}=\frac{1}{4}\{c_0^{}(1\mathrm{cos}2\theta )+\frac{c_1}{3}[4p\mathrm{sin}^3\theta +q(43\mathrm{cos}\theta \mathrm{cos}3\theta )]$$

$$+c_2[\frac{pq}{4}(4\theta \mathrm{sin}4\theta )+\frac{(q^2p^2)}{2}(2+\mathrm{cos}2\theta )\mathrm{sin}^2\theta ]\}.$$
$`(112)`$
Here the coefficient $`c_0^{}=c_0+c_2(3p^2+q^2)/4`$ depends on $`\varphi `$ and $`\lambda `$. Using eqs.(82,111,112) we obtain

$$T_w=\{\frac{c_0^{}}{2}(1\theta \mathrm{cot}\theta )\frac{c_1}{3}(1\mathrm{cos}\theta )[q+p\mathrm{tan}(\frac{\theta }{2})]$$

$$+\frac{c_2}{8}[(p^2q^2)\mathrm{sin}^2\theta pq(2\theta \mathrm{sin}2\theta )]\}.$$
$`(113)`$
Using eq.(95), with $`T_w`$ given by eq.(113), we obtain the shell’s shape

$$R=R_0\mathrm{csc}\theta \{3(1\theta \mathrm{cot}\theta )(c_0+\frac{c_2}{4}(3p^2+q^2))$$

$$+2c_1(1\mathrm{cos}\theta )[q+p\mathrm{tan}(\frac{\theta }{2}\left)\right]$$

$$+\frac{3c_2}{4}[(q^2p^2)\mathrm{sin}^2\theta +pq(2\theta \mathrm{sin}2\theta )]\}^{1/2},$$
$`(114)`$
where the explicit dependence upon stellar inclination and azimuthal angle is obtained using $`p=\mathrm{sin}\varphi \mathrm{sin}\lambda `$ and $`q=\mathrm{cos}\lambda `$. Examples are shown in Figures 4-6 for even parity winds with $`c_2=3,1.5`$, and $`1.5`$. Noting that Bandiera’s notation for this problem was $`\mathrm{\Delta }_Y=2c_2/3`$, these correspond to his cases of $`\mathrm{\Delta }_Y=2,1`$ and $`+1`$.

In the neighborhood of the standoff point, the shell’s shape is given by

$$R^2/R_0^2[(c_0+c_1q+c_2q^2)+\frac{p}{2}(c_1+2c_2q)\theta $$

$$+(8c_0+5c_1q+2c_2(3p^2+q^2))\theta ^2/20].$$
$`(115)`$
### 4.4 The Non-Axisymmetric Wind

The previous formulation is not restricted to an axisymmetric wind, so one may obtain solutions for bow shocks driven by radial, non-axisymmetric winds, by defining the momentum flux from the wind in terms of spherical harmonics. Because the method is clear from the previous discussion of an axisymmetric wind, we will not give further details here. We note only that the response from each spherical harmonic contribution to the wind momentum flux may be obtained from the required torque formula. Then the shape of the shell will depend upon the square root of the sum of these responses.

## 5 Energy Thermalization Rates

The energy thermalization rate of axisymmetric bow shocks and colliding winds has been considered by Wilkin, Cantó & Raga (2000). The thermalization occurs in the shocks, post-shock (radiative) relaxation, and mixing. The maximum possible rate of thermalization is given by the total incident kinetic energy flux in the center of mass frame. Here we are concerned with the total rate of energy thermalized, over the entire shell, rather than the amount per unit area. By center of mass, we mean the center of mass of the matter deposited onto the shell per unit time. This maximum thermalization rate is for the case of complete mixing. For the bow shock in a uniform ambient medium, the center of mass frame corresponds to the ambient frame, because a formally infinite amount of ambient mass per unit time strikes the shell, while the stellar wind has finite mass loss rate. For the bow shock driven by a star with a radial wind, the total energy thermalization rate is given by

$$\dot{E}_{\mathrm{therm}}=_0^{2\pi }_0^\pi \frac{1}{2}R^2\rho _wV_w^2(𝐕_w^{}𝐕_{\mathrm{shell}}^{})\widehat{𝐫}\mathrm{sin}\theta d\theta d\varphi .$$
$`(116)`$
Here primes indicate velocities in the ambient frame, so $`𝐕_w^{}=V_w\widehat{𝐫}+V_a\widehat{𝐳}`$, giving

$$V_w^2=[V_w^2+2V_wV_a\mathrm{cos}\theta +V_{}^2].$$
$`(117)`$
The velocity of the shell is its pattern speed and corresponds to the stellar velocity $`𝐕_{}=V_a\widehat{𝐳}`$. Because we are assuming a non-accelerating wind, the integral depends only upon the streamlines and not the detailed shape of the shell, so we may perform the integration over a spherical surface.

For the case of an isotropic wind from a star moving through a plane-parallel stratified ambient medium, the rate of energy thermalization is precisely the same as for a uniform medium (Wilkin, Cantó, & Raga 1997), since it is the kinetic energy in the ambient frame that is thermalized:

$$\dot{E}_{\mathrm{therm}}=\frac{1}{2}\dot{M}_w(V_w^2+V_{}^2).$$
$`(118)`$
This rate of energy thermalization is appropriate for an arbitrary ambient distribution of matter, provided only that it is truly infinite in mass, and that the bow shock is steady-state in the star’s frame. The second requirement demands that the ambient mass distribution be independent of the z-coordinate, although its $`(\varpi ,\varphi )`$ distribution is arbitrary. An exception, for example, is the case of an exponential stratification, for which case the total ambient mass flux may be finite (on one side of the bow shock). For this case, the center of mass frame is not equal to the ambient frame, and it is important to separately define the center of mass frame for each $`\varphi `$ slice, as it will depend upon azimuthal angle for the portion of the bow shock that contains finite opening angle. Because of the resulting opening angle, not all of the stellar wind kinetic energy (in the ambient frame) will be thermalized, for two reasons: (1) some streamlines miss the shell, and (2) there is a net flow of momentum along the shell even at infinite distance, whereas for the standard bow shock $`V_t0`$ in the ambient frame, far in the bow shock tail.

We now consider the case of an anisotropic wind, to determine whether the energy thermalization rate depends on the orientation of the wind. However, the thermalization rate remains independent of the stratification of the ambient medium, subject to the caveat about finite-mass systems. For the remainder of this section, we allow the wind to be non-axisymmetric, so in equations (75,76) we use $`f_w=f_w(\theta _{},\varphi _{}),g_w=g_w(\theta _{},\varphi _{})`$. When we refer to an axisymmetric wind, its axis of symmetry will be $`\theta _{}=0`$ as before.

The total kinetic energy loss rate of the wind is

$$\dot{E}_w=\frac{1}{8\pi }\dot{M}_w\overline{V}_w^2_0^{2\pi }_0^\pi \frac{g_w^2}{f_w}\mathrm{sin}\theta _{}d\theta _{}d\varphi _{}.$$
$`(119)`$
The functions $`f_w`$ and $`g_w`$ are assumed to have unit average over $`4\pi `$ steradians. In terms of this energy loss rate, we define the mean square wind speed as

$$<V_w^2>=\mathrm{\hspace{0.17em}2}\dot{E}_w/\dot{M}_w.$$
$`(120)`$
We also define a total vector momentum flux for the wind, a quantity that is non-vanishing only if the wind is not point-symmetric with respect to the origin

$$𝐏_w=\frac{\dot{M}_w\overline{V}_w}{4\pi }_0^{2\pi }_0^\pi g_w\widehat{r}\mathrm{sin}\theta d\theta d\varphi ,$$
$`(121)`$
$$=_0^{2\pi }𝚽_w(\pi ,\varphi )𝑑\varphi .$$
$`(122)`$
We now wish to calculate the kinetic energy deposition rate to the shell in the ambient frame. Letting $`\alpha =V_{}/\overline{V}_w`$, the wind velocity (squared) in the ambient frame becomes

$$V_w^2=\overline{V}_w^2\left[\frac{g_w^2}{f_w^2}+2\alpha \frac{g_w}{f_w}\mathrm{cos}\theta +\alpha ^2\right].$$
$`(123)`$
The wind density is given by

$$4\pi R^2\rho _w=\frac{\dot{M}_w}{\overline{V}_w}\frac{f_w^2}{g_w}.$$
$`(124)`$
The total kinetic energy flux onto the shell in the ambient frame is

$$\dot{E}_{\mathrm{therm}}=\frac{\dot{M}_w\overline{V}_w^2}{8\pi }_0^{2\pi }_0^\pi \left[\frac{g_w^2}{f_w}+2\alpha g_w\mathrm{cos}\theta +\alpha ^2f_w\right]\mathrm{sin}\theta _{}d\theta _{}d\varphi _{}.$$
$`(125)`$
The normalization condition for $`f`$ permits the last term to be easily integrated. The second term is $`V_{}`$ times the $`z`$ component of the wind momentum flux. In terms of the total wind kinetic energy flux and momentum flux, we have

$$\dot{E}_{therm}=\frac{1}{2}\dot{M}_w(<V_w^2>+V_{}^2)+𝐏_w𝐕_{}.$$
$`(126)`$
We see that a sufficient condition for the total thermalization rate to be independent of the wind orientation is the vanishing of $`P_{z,w}`$. Because we do not wish to invoke a special alignment of the wind orientation with the star’s direction of motion, more typically this requires a point-symmetric wind such that $`𝐏_w=0`$. As a special case, this can be more explicitly confirmed for the axisymmetric wind, assuming it is symmetric with respect to its midplane $`\theta _{}=\pi /2`$. Because we are considering the total thermalization rate, we may perform the solid angle integration in starred coordinates. In this case, we need only the transformation equation for $`\mathrm{cos}\theta `$ which is given by eq.(74). We assume the wind to be axisymmetric and symmetric with respect to $`\theta _{}=\pi /2`$, so it has an expansion in $`\mathrm{cos}\theta _{}`$ with only even powers. This implies that the $`2\alpha \mathrm{cos}\theta `$ term in the integral will give no contribution, by parity, because it yields only a $`\mathrm{sin}\varphi _{}`$ term, which vanishes in the azimuthal integration and terms with odd powers of $`\mathrm{cos}\theta _{}`$, which vanish in the polar angle integration. The remaining contribution is then

$$\dot{E}_{\mathrm{therm}}=\frac{\dot{M}_w\overline{V}_w}{4}_0^\pi \left[\frac{g_w^2(\theta _{})}{f_w(\theta _{})}+\alpha ^2f_w(\theta _{})\right]\mathrm{sin}\theta _{}d\theta _{},$$
$`(127)`$
$$=\frac{1}{2}\dot{M}_w(<V_w^2>+V_{}^2).$$
$`(128)`$
The total energy thermalization rate is independent of the orientation of the midplane-symmetric, axisymmetric wind. For the case of a non-point-symmetric wind, including an axisymmetric wind that doesn’t possess mirror symmetry with respect to the equator, the total thermalization rate will depend upon the stellar orientation because of the contribution from the $`2\alpha \mathrm{cos}\theta `$ term. This is the case of the example solution given in Section 4.3, for the wind with quadratic dependence on $`\mathrm{cos}\theta _{}`$. For that wind, the total vector momentum loss rate is

$$𝐏_w=\frac{\dot{M}_w\overline{V}_w}{3}c_1\widehat{𝐳}_{},$$
$`(129)`$
where $`\widehat{𝐳}_{}`$ points in the direction of the wind symmetry axis $`z_{}`$. The thermalization rate for the bow shock driven by this wind therefor has a fluctuation depending on orientation of the wind symmetry axis with respect to the stellar velocity of magnitude

$$\mathrm{\Delta }\dot{E}_{therm}=\frac{\dot{M}_w\overline{V}_wV_{}}{3}c_1\mathrm{cos}\lambda .$$
$`(130)`$
If the quadratic term in the wind momentum vanishes ($`c_2=0`$), then the absolute value of $`c_1`$ may not exceed unity. The amount this changes the total thermalization rate then depends upon the parameter $`\alpha =V_{}/V_w`$.

We must note that changing orientation of the point-symmetric wind did not change the total thermalization rate, but the spectrum of shock emission will clearly be different. For example the peak post-shock temperature will depend upon the incident normal component of velocity. The sum of thermalization by shock, relaxation, and mixing is independent of orientation, but the individual contributions will vary. This highlights the fact that the properties of bow shocks due to non-axisymmetric winds may be determined by measuring the bow shock shape, mass surface density (or column density), kinematics, and radiated energy. Although detailed shock calculations are beyond the scope of this work, a substantial amount of information should be obtainable from these global quantities and may be sufficient for sources where the data are sparse.

## 6 Summary

I have shown how to solve the problem of non-axisymmetric bow shocks and wind collisions with a simple formalism that requires only algebraic equations. Most often one considers constant wind speed for such problems, and in this case the method leads to exact, analytic (although possibly implicit) solutions. The availablility of simple analytic solutions makes it much easier to model observed sources and derive the properties of the driving winds. Among the principal applications of these solutions would be to determine the cause of the asymmetry in observed bow shocks, whether it be due to the ambient medium or an anisotropic wind, or both.

Future improvements are needed to include non-radial and accelerating winds and shearing motions in the shell, in which case the fluid elements are not restricted to a plane. Also, non-axisymmetric bow shocks due to colliding winds in binary systems, including the orbital motion, require a formulation in a non-inertial frame. A forthcoming paper will describe how to solve the general problem of colliding winds from two stars, including both anisotropy and acceleration in the winds.

Portions of this work were done as an NRC Research Associate at NASA Ames Research Center, and at IPAC with support by Long Term Space Astrophysics grant to D.Van Buren. IPAC is operated by JPL and the California Institute of Technology under contract with the National Aeronautics and Space Administration. I am greatful for helpful discussions with J.Anderson and C.McKee.