# Bipolar molecular outflows driven by hydromagnetic protostellar winds

## 1. Introduction

As they form, young stars emit powerful winds whose mechanical luminosity amounts to a fair fraction of the stars’ binding energy. These winds are often observed as jets, and it is generally believed that a toroidal magnetic field is responsible for their collimation (Benford 1978, Blandford & Payne 1982). When a protostellar wind strikes the ambient medium, a bipolar molecular outflow is produced. These outflows have been credited with the support (Norman & Silk (1980), McKee (1989), Bertoldi & McKee (1996)) as well as the disruption (e.g., Bally et al. (1999)) of molecular clouds and the dense clumps within them that produce star clusters. We seek a model of protostellar outflows that is sufficiently detailed to permit a quantitative study of both their supportive and destructive roles in the lives of their parent clouds.

Bipolar molecular outflows display a number of common features that must be reproduced in any viable model. As summarized by Lada & Fich (1996), these include a nearly linear position-velocity relation (a “Hubble law”) and a mass-velocity relation $`dM/dvv^\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }1.8`$. It is debated whether these features result primarily from turbulent mixing in the region affected by the jet (e.g., Raga et al. (1993), Stahler (1994)), or from the dynamics of a shocked shell bounding the wind cocoon (e.g., Shu et al. (1991), Masson & Chernin 1992, 1993). In this paper we investigate shells of ambient material set into motion by hydromagnetic protostellar winds, and demonstrate that these naturally produce the primary characteristics of observed outflows. (We do not address traces of high-velocity molecular gas sometimes found within these shells.)

Our analysis of the wind force distribution (§2) follows Shu et al. (1995) and the suggestion made by Ostriker (1997); our investigation of shell motion (§3) follows Shu et al. (1991), Masson & Chernin (1993), and Li & Shu (1996). We generalize these results and show (§4), contrary to the conclusions of Masson & Chernin, that a combined model can reproduce the observed mass-velocity relationship as well as the position-velocity law.

## 2. Force distributions from disk winds and X-winds

We seek the distribution of wind momentum flux on scales of the dense clumps in molecular clouds, far larger than any scale associated with an accretion disk whose wind produces the outflow. The angular distribution of the wind momentum injection rate $`\dot{p}_w`$ at time $`t`$ can be written as

$$\frac{d\dot{p}_w(t)}{d\mathrm{\Omega }}=r^2\rho _wv_w^2\frac{\dot{p}_w(t)}{4\pi }P(\mu ),$$
(1)   
where $`\mu =\mathrm{cos}\theta `$ labels directions from the outflow axis. We wish to determine $`P(\mu )`$, the normalized force distribution, which we assume is constant in time.

Shu et al (1995) have shown that $`\rho _w1/(r\mathrm{sin}\theta )^2`$ is a good approximation for X-winds. Since the wind velocity $`v_w`$ is approximately the same on different streamlines in this model, it follows that $`\rho _wv_w^21/(r\mathrm{sin}\theta )^2`$. In fact, this should be a reasonably good approximation for any radial hydromagnetic wind that has expanded to a large distance, as the following heuristic argument indicates. Such winds expand more rapidly than the fast magnetosonic velocity $`c_fB/(4\pi \rho _w)^{1/2}`$, where we have assumed that thermal pressure is negligible. At large distances, the field wraps into a spiral with $`BB_\varphi `$. First, consider the flow along streamlines. Flux conservation gives $`2\pi r\mathrm{\Delta }rB_\varphi =`$ const in radial flow. Since $`v_w`$ is about constant at large distances, it follows that $`\mathrm{\Delta }r`$ const, so that $`B_\varphi 1/r`$. Since $`\rho _w1/r^2`$ in a constant velocity, radial wind, it follows that $`v_w/c_f`$ is approximately constant along a streamline at large distances. The flow approximates an isothermal wind, for which $`(v_w/c_f)^2`$ increases logarithmically with distance from unity near the source. This implies that at large distances $`(v_w/c_f)^2`$ is approximately a constant across streamlines as well, so that $`\rho _wv_w^2B_\varphi ^2`$. But note that the field must become approximately force free at large distances: balancing the tension $`B_\varphi ^2/(4\pi \varpi )`$, where $`\varpi r\mathrm{sin}\theta `$ is the cylindrical radius, against the pressure gradient $`(1/8\pi )B_\varphi ^2/\varpi `$ gives $`B_\varphi 1/\varpi =1/r\mathrm{sin}\theta `$. We conclude that $`\rho _wv_w^2B_\varphi ^21/(r\mathrm{sin}\theta )^2`$ is a general characteristic of radial hydromagnetic winds.

This argument can be made more precise using a result due to Ostriker (1997). Suppose the wind arises from a Keplerian disk, where the wind density varies with initial radius $`\varpi _0`$ as $`\rho _0\varpi _0^q`$ and the Alfvén velocity varies with orbital velocity at the disk. The value $`q=3/2`$ corresponds to the solution of Blandford & Payne (1982), whereas values $`0.5<q<1`$ were considered by Ostriker (1997). Assume that each streamline expands to $`\varpi \varpi _0`$, so that the wind is significantly super-Alfvénic. The conservation of specific energy, angular frequency, and mass flux along streamlines, along with the “isorotation” relation between $`𝐁`$ and $`𝐯`$, give $`\rho _wv_w^2C(r)\varpi _0^{(1q)/2}\varpi ^2`$, where $`C(r)`$ is a slowly varying function of $`r`$ (Ostriker 1997). We see that the heuristic argument above is valid provided that the disk radius $`\varpi _0`$, which varies from $`\varpi _{0,\mathrm{in}}`$ to $`\varpi _{0,\mathrm{out}}`$, say, has a smaller range of variation than $`\varpi `$, which varies from an innermost radius $`r\theta _{\mathrm{core}}`$ to $`r`$. More precisely, if there is a power law relation between $`\varpi _0`$ and $`\mathrm{sin}\theta `$, then $`P(\mu )r^2\rho _wv_w^2`$ gives

$$P(\mu )(\mathrm{sin}\theta )^{2(1+ϵ)};ϵ\left(\frac{q1}{4}\right)\frac{\mathrm{ln}(\varpi _{0,\mathrm{in}}/\varpi _{0,\mathrm{out}})}{\mathrm{ln}(\theta _{\mathrm{core}})},$$
(2)   
so we recover our earlier result if $`\theta _{\mathrm{core}}\varpi _{0,\mathrm{in}}/\varpi _{0,\mathrm{out}}`$, which allows the flow to be quasi-radial. For an X-wind, $`\varpi _{0,\mathrm{in}}=\varpi _{0,\mathrm{out}}`$ and the heuristic result should be quite accurate. Conditions at the axis set $`\theta _{\mathrm{core}}`$. Although the inner boundary is generally not considered in disk wind models, the Shu et al. (1995) theory posits a core of open field lines from the pole of the accreting star. The balance of magnetic pressure and tension in the fiducial X-wind model gives $`\varpi _{\mathrm{core}}=2.5(1+0.18\mathrm{log}_{10}r_{\mathrm{pc}})\mathrm{AU}`$, or $`\theta _{\mathrm{core}}=1.2\times 10^5r_{\mathrm{pc}}^1(1+0.18\mathrm{log}_{10}r_{\mathrm{pc}})`$, at a distance of $`r_{\mathrm{pc}}`$ parsecs.

We have evaluated the accuracy of equation (2) by calculating the wind force distribution analytically using the method outlined by Shu et al. (1995), as Ostriker (1997) suggested. We present this calculation in Matzner (1999), where we find that equation (2) is a good approximation at all angles more than fifteen degrees from the equator, and best matches the actual solution between one and ten degrees from the outflow axis. The fiducial X-wind model (Shu et al. 1995) has $`P(\mu )\mathrm{sin}^2\theta `$ greater by $`40\%`$ toward the equator, for reasons not considered in equation (2); this corresponds to $`ϵ1/50`$ within $`10^{}`$ of the axis at a distance of $`0.1\mathrm{pc}`$.

Although the approximation $`P(\mu )(\mathrm{sin}\theta )^{2(1+ϵ)}`$ is valid for ideal axisymmetric winds, we expect the force distribution to flatten within some angle $`\theta _0>\theta _{\mathrm{core}}`$ in reality. This flattening could be produced by any of the mechanisms thought to create Herbig-Haro objects, e.g., jet precession, internal shocks from a fluctuating wind velocity, or the magnetic kink instability. Assuming that $`ϵ`$ is negligible and that $`\theta _01`$, we can therefore approximate the force distribution of a magnetized protostellar wind as

$$P(\mu )\frac{1}{\mathrm{ln}(2/\theta _0)\left(1+\theta _0^2\mu ^2\right)},$$
(3)   
where the prefactor assures $`_0^1P(\mu )𝑑\mu =1`$. Outflows may require a larger value of $`\theta _0`$ than appropriate for winds themselves, if mixing between sectors (neglected here) dilutes the momentum on axis; the current theory will still apply, with this larger $`\theta _0`$. However, so long as $`\theta _01`$, the wind force is tightly concentrated along the axis: The formation of jets is an inevitable consequence of a hydromagnetic wind.

## 3. Shells driven by protostellar winds

We shall now explore the structure and motion of a shell of ambient material struck by a wind with the force distribution given by equation (3). Following Shu et al. (1991), Masson & Chernin (1992), and Li & Shu (1996), we idealize the swept-up shell as thin and momentum-conserving. This is justified because both shocks are radiative for protostellar wind velocities (Koo & McKee 1992a; however, see §5 for a consideration of magnetic pressure). We will also adopt the assumption that the flow is entirely radial, so that mass and momentum are conserved in each angular sector and there is no relative motion of the shocked fluids.

A shell is driven by a “heavy” wind if less ambient material than wind material has been swept up; such shells travel at nearly the wind velocity, and the crossing time of the wind is therefore comparable to the outflow age. Alternatively, a shell driven by a “light” wind is one that has swept up more ambient gas, and has decelerated significantly. In the limit of a very light wind, both the wind’s mass and its flight time can be neglected; this limit is approached rapidly once a comparable mass has been swept up (Koo & McKee 1992b). Molecular outflows expand five to twenty times more slowly than their driving winds, and are comparably more massive. We may therefore neglect the wind’s mass and flight time, and integrate the equation of momentum conservation in each direction $`\mu `$, $`dp_w/d\mathrm{\Omega }=v_sM_a(R_s,\mu )/\mathrm{\Omega }`$, where $`M_a`$ is the ambient mass inside the shell radius $`R_s`$. We assume that the ambient gas has a density $`\rho _a=\rho _{01}r^{k_\rho }Q(\mu )`$, where $`\rho _{01}`$ is a constant and the angular factor $`Q(\mu )`$ is normalized so that $`_0^1Q(\mu )𝑑\mu =1`$. We find that the shell radius is

$$R_s^{4k_\rho }=\frac{(4k_\rho )(3k_\rho )P(\mu )}{4\pi \rho _{01}Q(\mu )}_0^tp_w(t^{})𝑑t^{}.$$
(4)   
If the wind momentum is a power law in time, $`p_wt^{\eta _{\mathrm{in}}}`$, then the shell velocity is given by $`v(\mu )=\eta R_s/t`$, where

$$\eta \mathrm{ln}R_s(\mu ,t)/\mathrm{ln}t=(\eta _{\mathrm{in}}+1)/(4k_\rho ).$$
(5)   
Equation (4) shows that the shell expands self-similarly, regardless of the wind history: its radial and velocity structures are fixed, while its scale expands as $`p_w𝑑t`$ increases. Self-similarity is expected, because we have chosen a scale-free medium; other radial scales, such as the scale of the light-heavy wind transition, the wind collimation scale, and the cooling length, are all small compared to a typical outflow.

Shu et al. (1991), Masson & Chernin (1992) and Li & Shu (1996) have all considered a steady wind ($`\eta _{\mathrm{in}}=1`$) and an ambient distribution appropriate for pre-stellar cores at the point of collapse: $`k_\rho =2`$, and $`Q(\mu )`$ larger toward the equator, because of magnetic or rotational flattening. However, this is appropriate only in the region whose gravity is dominated by the pre-stellar core, or where it prescribes a density higher than the ambient density: $`r0.07(\sigma _{\mathrm{th}}/0.2\mathrm{km}\mathrm{s}^1)(n_\mathrm{H}/10^4\mathrm{cm}^3)^{1/2}\mathrm{pc}`$, using the theory of Shu (1977). <sup>1</sup><sup>1</sup>1 Note that if the wind, star and core masses are each within about a factor of ten of the last, then the outflow must have emerged from its core if the wind is to be light near its axis. This follows from the overwhelming factor ($`10^4`$) by which the axial force is enhanced; it is only exacerbated by any flattening of the core. The core mass may continue to affect the low-velocity, equatorial flow, however.

Outside this radius, anisotropies in the ambient gas are unlikely to correlate with the outflow direction, so we may assume $`Q(\mu )1`$. We also then expect $`k_\rho 0`$, $`1`$ or $`2`$, if the lobe in question is smaller than, comparable to, or emerging from its parent star-forming “clump”. Assuming the ambient medium is isotropic and the outflow has not escaped its clump, and that $`ϵ1`$ so that its effect can be neglected in the outflow shape, our model reduces to:

$$\frac{R_s}{R_{\mathrm{head}}}=\left[1+(1\mu ^2)\theta _0^2\right]^{1/(4k_\rho )},$$
(6)   
where the radius of the lobe head, $`R_{\mathrm{head}}`$, expands according to equations (3) and (4) with $`\mu =1`$.

## 4. Comparison with observations

The shell described above is a Hubble flow in the sense that $`𝐯_s(\mu ,t)=\eta 𝐑_s(\mu ,t)/t`$. Since the relative line-of-sight velocity is related to the line-of-sight distance by $`v_{\mathrm{obs}}=\eta z_{\mathrm{los}}/t`$, the position-velocity (PV) diagram along the extent of an optically thin outflow is the same as its image to an observer situated in the plane of the sky along the short axis of the outflow. Whereas inclination causes a foreshortening of the outflow in a sky map, it causes the PV diagram to rotate. An elongated outflow ($`\theta _01`$) will display a nearly linear PV diagram from the greatest to the least values of velocity and position. The agreement of self-similar outflows with the observed Hubble law was pointed out by Shu et al (1991) for the particular case $`k_\rho =2`$ and $`\eta _{\mathrm{in}}=1`$; here we see that it is a quite general property.

The mass-velocity relationship, $`dM/dv_{\mathrm{obs}}`$ is a projection of the PV diagram onto the velocity axis. Typical outflows show $`dM(v_{\mathrm{obs}})/dv_{\mathrm{obs}}v_{\mathrm{obs}}^\mathrm{\Gamma }`$, with $`\mathrm{\Gamma }1.8`$. In an elongated outflow of inclination $`i`$, $`v_{\mathrm{obs}}v\mathrm{cos}i`$, where $`v=|𝐯|`$, because all but the lowest velocities are achieved at small angles from the outflow axis. Therefore, $`\mathrm{\Gamma }=d\mathrm{ln}(dM/dv)/d\mathrm{ln}(v)`$ except at low velocities. Because $`dM/dv=(M/\mu )/(v/\mu )v^{3k_\rho }/(\mu v^{(5+ϵk_\rho )/(1+ϵ)})`$ from equations (3) and (4) with $`Q(\mu )=1`$, we find

$$\mathrm{\Gamma }=2+ϵ\frac{4k_\rho }{1+ϵ},$$
(7)   
for $`v_{\mathrm{obs}}2v_{\mathrm{min}}`$ where $`v_{\mathrm{min}}\theta _0^{2/(4k_\rho )}v_{\mathrm{head}}`$ is the minimum space velocity (see eq. 6). Essentially all of the (non-equatorial) flow shows a value of $`\mathrm{\Gamma }`$ very close to $`2`$ (if $`ϵ1`$), in excellent agreement with a typical observed value of $`1.8`$. The slope is shallower for $`v_{\mathrm{obs}}v_{\mathrm{min}}`$: because $`dM/dv_{\mathrm{obs}}`$ is a symmetric function of $`v_{\mathrm{obs}}`$, $`\mathrm{\Gamma }=0`$ when $`v_{\mathrm{obs}}=0`$. The velocity at which $`dM/dvv^2`$ fails is greater for more inclined outflows, or for larger $`\theta _0`$; this could in principle constrain the inclination. The ability of our model to fit observational mass-velocity curves is demonstrated in Figure 1 for L1551, NGC2071, and NGC2264G.

The qualitative agreement of X-winds in isothermal toroids with observed mass-velocity curves has previously been shown by Li & Shu (1996). Again, we see that this is a general feature of momentum–conserving shells driven by hydromagnetic winds, and is not linked to particular models for the source of the wind nor the ambient medium, so long as this medium has a power law density distribution.

The broadening angle $`\theta _0`$ can be constrained using the range of velocities for which $`\mathrm{\Gamma }2`$: this law holds for at least a factor of 4 in the outflows shown in Figure 1. This must be less than about $`v_{\mathrm{head}}/(2v_{\mathrm{min}})=\theta _0^{2/(4k_\rho )}`$, so $`\theta _010^{1.8(1k_\rho /4)}`$. However, the velocity factor is also reduced by inclination at a given $`\theta _0`$; in Figure 1 we show that $`\theta _0=10^2`$ allows an inclination of $`45^{}`$ to fit the data; $`\theta _0=10^{1.5}`$ gives too small a velocity range unless $`i=0`$.

If $`p_{\mathrm{obs},\mathrm{lobe}}`$ is the net momentum of an outflow lobe along the line of sight, another constraint on $`\theta _0`$ comes from the ratio,

$`{\displaystyle \frac{p_{\mathrm{obs},\mathrm{lobe}}}{(v_{\mathrm{obs}}^2dM/dv_{\mathrm{obs}})_{\mathrm{head}}}}={\displaystyle \frac{2\mathrm{ln}(\theta _0^1)}{4k_\rho }},`$ (8)   
valid if $`ϵ,\theta _01`$. The result was obtained by integrating the momentum along the line of sight for an uninclined ($`i=0`$) outflow; we did not assume that $`\mu 1`$, but instead used the exact expression $`v_{\mathrm{obs}}=\mu v`$ appropriate for each shell. It is interesting to note that this agrees exactly with the expression $`\mathrm{ln}(v_{\mathrm{head}}/v_{\mathrm{min}})`$ one would estimate from $`dm/dvv^2`$ with $`\mu =1`$. Although the ratio was derived for zero inclination, it is actually independent of $`i`$, because the numerator and denominator scale together. Observations place a lower limit on the ratio (an upper limit on $`\theta _0`$), since the net momentum may be underestimated. We find that $`\theta _0<10^{1.5}`$ for L1551 (Moriarty-Schieven & Snell (1988)) and $`\theta _0<10^{1.3}`$ for NGC2071 (Moriarty-Schieven, Hughes & Snell (1989)), assuming $`k_\rho =0`$. Again, note that $`\theta _010^2`$ is consistent with the data.

The current model reproduces the typical extents and velocities of observed outflows. For instance, equations (3) and (4) imply that a wind of momentum $`20M_{}\mathrm{km}\mathrm{s}^1`$ with $`\theta _0=10^2`$, blowing steadily into a uniform density $`10^4\mathrm{cm}^3`$ for $`10^5`$ years, drives lobes whose heads decelerate to $`7.4\mathrm{km}\mathrm{s}^1`$ and expand to $`1.5\mathrm{pc}`$ (each) in this time.

## 5. Conclusions

Perhaps the most remarkable properties of bipolar molecular outflows, apart from their intensity and frequency in regions of active star formation, are their high degrees of collimation and the commonality of the relation $`dM/dvv^{1.8}`$. We have shown that hydromagnetic winds are naturally collimated so that the force distribution, $`\rho _wv_w^21/\mathrm{sin}^2\theta `$, approximately leads to $`dM/dvv^2`$ in any power-law ambient medium, provided the interaction is momentum conserving. Our results show that that the conclusions reached by Shu et al (1991) and (1995) and Li & Shu (1996) for steady X-winds in media with $`1/r^2`$ density distributions are far more general, and in addition are in good agreement with observations of protostellar outflows.

The fact that outflows are often observed to have $`dM/dv`$ slightly shallower than $`v^2`$ indicates that the theory is only approximate. A number of effects we have not considered could lead to a deviation from the -2 slope: CO self-absorption at lower velocities, mixing of radial momentum between angles, or the generation of lateral momentum when the wind impacts the shell (Masson and Chernin 1993). In the current model, $`\mathrm{\Gamma }`$ can differ from $`2`$ either because of the wind force does not exactly follow $`\rho _wv_w^2(\mathrm{sin}\theta )^2`$ ($`ϵ0`$), or because the ambient medium is not a single power law ($`k_\rho `$ varies). The model also predicts a flattening of $`dM/dv`$ at low velocities, which could raise the estimate of $`\mathrm{\Gamma }`$.

Let us consider the possibility that $`ϵ`$ is to blame for $`\mathrm{\Gamma }2`$. From equation (7), $`\mathrm{\Gamma }1.8`$ requires $`ϵ^15(4k_\rho )`$: the wind is slightly more concentrated toward the axis. Equation (2) implies that the number of decades of disk radius required to give this value of $`ϵ`$ is approximately $`\mathrm{log}_{10}(1/\theta _{\mathrm{core}})/[5(q1)(1k_\rho /4)]`$. The Blandford & Payne (1982) model has $`q=3/2`$; it would therefore require about 2-4 decades of disk radius to give $`\mathrm{\Gamma }=1.8`$ for $`\theta _{\mathrm{core}}10^5`$ and $`2k_\rho 0`$. A disk with $`q<1`$ (initial wind density increasing outward, e.g., Ostriker’s  models), has its wind force weighted toward the equator relative to $`\rho _wv_w^2(\mathrm{sin}\theta )^2`$, and produces outflows steeper than $`dM/dvv^2`$. X-wind models share this behavior, as they predict $`ϵ1/50`$.

Our model assumes that the shocked wind and ambient gas form a thin shell. Although this is appropriate for unmagnetized gas (Koo & McKee 1992a ), the fact that winds are collimated magnetically raises the possibility that outflows might become inflated with a cocoon of magnetically-supported shocked wind, before or after the wind shuts off. This depends on the wind’s terminal Alfvén Mach number (inversely proportional to Poynting flux), and also on whether the field remains ordered or becomes tangled. Because the Poynting flux decreases as the wind collimates ($`\theta _{\mathrm{core}}`$ decreases), and also as the kink instability develops (Choudhuri & Königl 1986), it is reasonable to ignore the magnetic pressure of the shocked wind.

###### Acknowledgements.

E. Ostriker informs us of similar, unpublished results she obtained independently. We are grateful to F. Shu for thoughtful suggestions, and to R. Plambeck and M. Hogerheijde for sharing and discussing their data. CDM appreciates comments from J. Monnier. The research of both CDM and CFM is supported in part by the National Science Foundation through NSF grant AST 95-30480, in part by a NASA grant to the Center for Star Formation Studies, and, for CFM, in part by a Guggenheim Fellowship. CFM gratefully acknowledges the hospitality of John Bahcall of the Institute for Advanced Study; his visit there was supported in part by a grant from the Alfred P. Sloan Foundation.