Angular momentum in bipolar outflows: dynamical evolutionary model

TL;DR

This paper presents a dynamical model of bipolar molecular outflows driven by stellar winds interacting with rotating collapsing clouds, analyzing shell shapes, velocities, and angular momentum evolution.

Contribution

It introduces a time-dependent, numerical model of bipolar outflows considering wind anisotropy and rotation, providing new insights into shell dynamics and angular momentum transfer.

Findings

Shell elongation increases with wind anisotropy.

Rotation velocities are lower than observed in some sources.

A critical wind-to-accretion momentum ratio determines shell expansion.

Abstract

We model molecular outflows produced by the time dependent interaction between a stellar wind and a rotating cloud envelope in gravitational collapse, studied by Ulrich. We consider spherical and anisotropic stellar winds. We assume that the bipolar outflow is a thin shocked shell, with axial symmetry around the cloud rotation axis and obtain the mass and momentum fluxes into the shell. We solve numerically a set of partial differential equations in space and time, and obtain the shape of the shell, the mass surface density, the velocity field, and the angular momentum of the material in the shell. We find that there is a critical value of the ratio between the wind and the accretion flow momentum rates $\beta$ that allows the shell to expand. As expected, the elongation of the shells increase with the stellar wind anisotropy. In our models, the rotation velocity of the shell is the order to 0.1 - 0.2 km s$^{-1}$, a factor of 5-10 lower than the values measured in several sources. We compare our models with those of Wilkin and Stahler for early evolutionary times and find that our shells have the same sizes at the pole, although we use different boundary conditions at the equator.

Tables

Table 1. Table 1: Observational parameters of the molecular outflows with rotation.

Source	Observations	$M_{\ast }$	$\Delta v_{\phi }$	$\Delta r$	zcut
		(M⊙)	(km s-1)	(AU)	(AU)
CB26	HCO+(1-0) and 13CO (2-1)	0.5	1.1	100	560
Ori-S6	SO (65-54) and 12CO (2-1)	2.0	2.0	1000	1200
HH 797	12CO (2-1)	1.0	2.0	1000	7500
DG Tau B	12CO (2-1)	0.5	1.0	150	450
Orion Source I	Si18O and H2O	8.7	5.0	80	150
HH 30	12CO (2-1) and 13CO (2-1)	0.45	0.4	150	200
NGC 1333 IRAS 4C	CCH	0.18	0.4	470	700

Table 2. Table 2: Values of β crit subscript 𝛽 crit \beta_{\mathrm{crit}} for different values of α 𝛼 \alpha and initial radius r s ​ 0 ​ ( 0 ) subscript 𝑟 𝑠 0 0 r_{s0}(0) .

\diagbox[innerwidth=1cm,innerleftsep=-0.75cm,innerrightsep=-5pt]$\alpha$$r_{s0}\left(0\right)$	${10}^{-5}$	5$\times {10}^{-5}$	${10}^{-4}$
0.01	2.061	0.972	0.900
0.10	20.084	9.140	6.534
0.50	99.732	45.085	32.023

Equations

$$P_m=R_s^2\sin \theta \sigma ,$$

$$P_r$$

$$P_{\theta }$$

$$P_{\varphi }$$

$$\frac{\partial P_m}{\partial t}+\frac{\partial }{\partial \theta }\left(\frac{P_{\theta }}{R_s}\right)=$$

$$R_s^2\sin \theta \left[{\rho }_a\left(\frac{P_r}{P_m}-U_{ar}\right)-{\rho }_w\left(\frac{P_r}{P_m}-v_w\right)\right],$$

$$\frac{\partial P_r}{\partial t}+\frac{\partial }{\partial \theta }\left(\frac{P_r{P}_{\theta }}{R_s{P}_m}\right)-\frac{P_{\theta }^2+P_{\varphi }^2}{R_s{P}_m}+\frac{G{M}_{\ast }{P}_m}{R_s^2}=R_s^2\sin \theta$$

$$\frac{\partial P_{\theta }}{\partial t}$$

$$\frac{\partial P_{\varphi }}{\partial t}$$

$$\frac{\partial R_s}{\partial t}=\frac{P_r}{P_m}-\frac{1}{R_s}\frac{P_{\theta }}{P_m}\frac{\partial R_s}{\partial \theta },$$

$$\zeta \equiv \frac{R_{\mathrm{cen}}}{r},$$

$$U_{ar}=-v_0{\zeta }^{1/2}{\left(1+\frac{\cos \theta }{\cos {\theta }_0}\right)}^{1/2},$$

$$U_{a\theta }=v_0{\zeta }^{1/2}\left(\frac{\cos {\theta }_0-\cos \theta }{\sin \theta }\right){\left(1+\frac{\cos \theta }{\cos {\theta }_0}\right)}^{1/2},$$

$$U_{a\varphi }=-v_0{\zeta }^{1/2}\frac{\sin {\theta }_0}{\sin \theta }{\left(1-\frac{\cos \theta }{\cos {\theta }_0}\right)}^{1/2},$$

$${\rho }_a=-\frac{{\dot{M}}_a{\zeta }^2}{4\pi R_{cen}^2{U}_{ar}}[1+2\zeta P_2(\cos {\theta }_0){]}^{-1},$$

$$v_0={\left(\frac{G{M}_{\ast }}{R_{\mathrm{cen}}}\right)}^{1/2},$$

$$\zeta =\frac{\cos {\theta }_0-\cos \theta }{{\sin }^2{\theta }_0\cos {\theta }_0}.$$

$${\rho }_w=\frac{{\dot{M}}_w}{4\pi r^2{v}_w}f(\theta ),$$

$$f(\theta )=\frac{A+B{\cos }^{2n}\theta }{A+B/(2n+1)}.$$

$$r_s=\frac{R_s}{R_{\mathrm{cen}}},$$

$$\tau =\frac{v_0}{R_{\mathrm{cen}}}t,$$

$$p_m=\frac{4\pi v_0}{{\dot{M}}_a{R}_{\mathrm{cen}}}{P}_m,$$

$$p_r=\frac{4\pi }{{\dot{M}}_a{R}_{\mathrm{cen}}}{P}_r,$$

$$p_{\theta }=\frac{4\pi }{{\dot{M}}_a{R}_{\mathrm{cen}}}{P}_{\theta },$$

$$p_{\varphi }=\frac{4\pi }{{\dot{M}}_a{R}_{\mathrm{cen}}}{P}_{\varphi }.$$

$$u_{ar}=\frac{U_{ar}}{v_0},$$

$$u_{a\theta }=\frac{U_{a\theta }}{v_0},$$

$$u_{a\varphi }=\frac{U_{a\varphi }}{v_0},$$

$$u_{wr}=\frac{v_w}{v_0}.$$

$$\alpha =\frac{{\dot{M}}_w}{{\dot{M}}_a},$$

$$\beta =\frac{{\dot{M}}_w{v}_w}{{\dot{M}}_a{v}_0}\equiv \alpha u_{wr}.$$

$${\rho }_a^{\prime }=\frac{4\pi R_{cen}^2{v}_0}{{\dot{M}}_a}{\rho }_a\equiv -\frac{{\zeta }^2}{u_{ar}}{\left[1+2\zeta P_2(\cos {\theta }_0)\right]}^{-1},$$

$${\rho }_w^{\prime }=\frac{4\pi R_{cen}^2{v}_0}{{\dot{M}}_a}{\rho }_w\equiv \frac{\alpha }{r_s^2{u}_{wr}}f(\theta ).$$

$$\frac{\partial p_m}{\partial \tau }+\frac{\partial }{\partial \theta }\left(\frac{p_{\theta }}{r_s}\right)=\sin \theta \times$$

Full text

accretion flow – angular momentum – hydrodynamics equations – molecular outflows – rotation velocity

\AuthorCallLimit

=1 \fullcollaborationNameThe Friends of AASTeX Collaboration

1 Introduction

The study of the molecular outflows and protostellar jets is fundamental to understand the star formation process. Molecular outflows probably limit the mass of the star-disk system (e.g., Shu et al. 1993) and can induce changes in the chemical composition of their host cloud (e.g., Bachiller 1996) since they are a mixture of entrained material from the cloud and the outflowing stellar wind (e.g., Snell et al. 1980).

The magneto-centrifugal mechanism (Blandford & Payne, 1982) is considered the principal candidate for producing the jets of young stars (see reviews by Königl & Pudritz 2000 and Shu et al. 2000). In this mechanism the magnetic field, anchored to the star-disk system, is responsable for accelerating the jet. However, it is still under debate where these magnetic fields are anchored to the disk: it could be at a narrow region at the truncation radius $R_{x}$ of the disk by a stellar magnetosphere (X-winds, e.g., Shu et al. 1994) or at a wider range of radii (disk winds, e.g., Pudritz & Norman 1983). Magnetohydrodynamic models predict that the material ejected from the disk has a toroidal angular momentum component related to the rotation at the disk foot point. Therefore, the observed rotational velocity of the jet can give information about its origin on the disk (Anderson et al., 2003). For example, Lee et al. (2009) and Lee et al. (2017) found that the protostellar jets HH 211 and HH 212, respectively, are ejected from very small radii, consistent with the X-wind model.

Molecular outflows have been explained as driven by fast stellar winds or as actual disk winds. In the former case, molecular outflows are produced when a fast stellar wind collides with the parent cloud accelerating and entraining cloud material (e.g., see reviews by Arce et al. 2007 and Bally 2016). In the latter case, the molecular outflow is ejected directly from the accretion disk (e.g., Pudritz & Norman 1986).

In recent years, rotation has been observed in a few molecular outflows, which are almost on the plane of the sky. 222As an alternative to the interpretation of a velocity difference as rotation, De Colle et al. (2016) showed that such differences can also be due to asymmetric shocks produced in the interaction of the stellar wind with the environment or by asymmetries in the ejection velocity of the disk-star system. These sources are: CB 26 (Launhardt et al., 2009), Ori-S6 (Zapata et al., 2010), HH 797 (Pech et al., 2012), DG Tau B (Zapata et al., 2015), Orion Source I (Hirota et al., 2017), HH 30 (Louvet et al., 2018), and NGC 1333 IRAS 4C (Zhang et al., 2018). Table 1 presents a summary of their characteristics.

From the observed rotation velocities, assuming that the outflows are disk winds, different authors obtained a disk launching radii between 10-50 AU (e.g., Launhardt et al. 2009; Pech et al. 2012). Nevertheless, Zapata et al. (2015), showed that magneto-centrifugal and photoevaporated disk winds do not have enough linear or angular momentum to account for the observed linear momentum and angular momentum rates in the molecular outflow of DG Tau B. They found that the observed rates are larger by a factor of 100, because the disk winds are not very massive. They pointed out that to account for the large masses of the observed molecular outflows they must be mainly entrained material from the parent cloud.

Several authors have modeled the molecular outflow as a wind-driven shell formed by the interaction between a radial stellar wind and the ambient cloud (e.g., Shu et al. 1991; Matzner & McKee 1999; Cantó et al. 2006). The ambient cloud can also be an accreting envelope. For example, Mendoza et al. (2004) described the hydrodynamical interaction between a rotating accretion flow and a spherically symmetric stellar wind. However, they did not consider neither the gravitational pull from the central star nor the centrifugal terms in the momentum equation. Later, Wilkin & Stahler (2003) took into account these effect but considered only the early evolution of the outflow.

Here we present a time dependent model of the interaction between a rotating accretion flow and a fast radial stellar wind. The molecular outflow is a thin shell driven by the fast stellar wind, that gains mass from both the stellar wind and the accretion flow. In this model we consider the gravitational pull of the central star and the centrifugal terms in the momentum equations. Also, we follow the evolution of the shell from the stellar surface up to large distances from the central star. We consider the molecular outflows produced by both isotropic and axisymmetric stellar winds with a polar angle dependance.

This paper is organized in the following way: in section 2 we show the equations of the dynamic evolution of the shell. Section 3 presents the description of the accretion flow and the stellar wind. The method of solution is presented in section 4, where we show the non dimensional equations and boundary conditions. In this section we also find semi analytic solutions for expansions around both the pole and the equator. Section 5 presents the results for different stellar winds. In section 6 we discuss our results. Finally, the conclusions are presented in section 7.

2 General formulation

We assume that the molecular outflow is formed by the supersonic collision between a stellar wind and an accretion flow. This collision leads to the formation of both an inner and an outer shock front. We assume that the cooling behind these shocks is relatively efficient because the shock velocities are expected to be less than $<$ 100 km s*-1* (Hartigan et al., 1987). Thus, the region between the shocks is described by a cold and thin shell. Within the shell, two fluids with a different density and velocity come into contact producing internal shearing layers which are subject to the Kelvin-Helmholtz instability, quickly leading to a turbulent mixing. Here we assume that the mixing is so efficient that one may describe the shell as a single fluid (e.g., Wilkin & Stahler 2003).

The evolution of the shell is governed by the fluxes of mass and momentum from the stellar wind and the accretion flow, by the gravitational influence of the central star, and by the centrifugal effects.

2.1 Shell equations

To derive the shell equations, we use spherical coordinates $r$, $\theta$, and $\phi$ for the radial, the polar, and the azimuthal coordinates, respectively. The coordinate system is centered on the star and we assume axial symmetry. The shell has a radius $R_{s}$, a mass surface density $\sigma$, and velocity components $U_{r}$, $U_{\theta}$, and $U_{\phi}$. All of these functions depend on $\theta$ and $t$, although, for simplicity, we will omit these dependences.

We assume that the accretion and wind flows are axisymmetric, and that they vary in a timescale much longer than the shell evolution time. Therefore, their properties depend only on the coordinates $r$ and $\theta$. The accretion flow has a mass volume density $\rho_{a}$ and velocity components $U_{ar}$, $U_{a\theta}$, and $U_{a\phi}$. The stellar wind has a mass volume density $\rho_{w}$ and velocity components $U_{wr}$, $U_{w\theta}$, and $U_{w\phi}$. Figure 1 shows the outflow model where a thin shell is formed by the interaction of the stellar wind and the accretion flow.

In Appendix A, we show the derivation of the equations of the shell evolution in a general form. In order to write these equations in more compact form, we define the mass flux

[TABLE]

and the momentum fluxes

[TABLE]

Also, we consider that the stellar wind has only a radial velocity component ($v_{w}=U_{wr}$).

Then, the continuity equation (eq. [A3]) can be written in terms of the mass and momentum fluxes as

[TABLE]

where the RHS shows the contribution to the shell mass from the stellar wind and the accretion flow.

The equation of the momentum in radial direction (eq. [A7]) is given by

[TABLE]

where, $G$ is the gravitational constant and $M_{*}$ is the stellar mass. The third term in the LHS comes from the centrifugal effect, and the last term is due to the weight of the shell. The RHS has the contribution from the stellar wind and the accretion flow.

The momenta in the $\theta$ and the azimuthal directions (eqns. [A8] and [A9]) in terms of the mass and momentum fluxes, respectively, can be written as

[TABLE]

[TABLE]

In these two equations, the last terms on the LHS are due to the centrifugal effect on the shell, while in the RHS only the accretion flow contributes to the momentum fluxes in these directions.

Finally, the evolution of the shell radius can be written as

[TABLE]

where the first term in the RHS corresponds to the radial velocity and the second term is the contribution of the tangential motion along the shell.

To solve these equations for the evolution of the shell radius $R_{s}(\theta,t)$, one needs to specify the properties of the accretion flow and the stellar wind. This model allows an accretion flow and a stellar wind with a general velocity field. The only constrain is that they have to be axisymmetric. We note that our formulation of the equations is different from Wilkin & Stahler (2003). In our model, the vectors are expressed in spherical coordinates, while in the model of Wilkin & Stahler (2003) the vectors are decomposed in directions orthogonal and parallel to the shell.

3 Accretion flow and stellar wind

In this section we will apply our model to flows with specific properties in order to solve for the shell evolution.

3.1 The accretion flow

The accretion flow is given by the gravitational collapse of a rotating cloud described by Ulrich (1976). In this model, the fluid particles have a uniform rotation rate at large distances and fall to the center conserving the specific angular momentum $j$. The collapse is assumed to be pressureless and thus, the orbits are ballistic. The latter assumption is true when the flow is supersonic, and heating by radiation and viscosity effects are negligible. The collapse of the gas reaches a centrifugal barrier at $R_{\mathrm{cen}}=j^{2}/GM_{*}$.

In terms of the non dimensional radial variable

[TABLE]

and the polar angle $\theta$, the velocity field and the density profile of the accretion flow are given by

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where $\theta_{0}$ is the initial polar angle of the orbit of the fluid element at the beginning of the collapse towards the center, $v_{0}$ is the free fall velocity

[TABLE]

$\dot{M}_{a}$ is the mass accretion rate, and the Legendre polynomial is $P_{2}(\cos\theta_{0})=\frac{1}{2}\left(3\cos\theta_{0}^{2}-1\right)$. The angle $\theta_{0}$ is given implicitly in terms of variables $\theta$ and $\zeta$ by

[TABLE]

Eq. (13) of Mendoza et al. (2004) gives an explicit solution of this equation.

Note that eq. (11) for the azimuthal velocity differes in sign with that given by Ulrich (1976). This only means that we consider the accretion flow to rotate with a negative angular momentum, as in Figure 1 of Mendoza et al. (2004).

3.2 The stellar wind

We assume an anisotropic stellar wind with a mass loss rate $\dot{M}_{w}$ and only a radial velocity component $U_{wr}=v_{w}$, assumed to be constant. The density is given by

[TABLE]

where $f(\theta)$ is anisotropy function given by

[TABLE]

The constants are $A\geq 0$ and $B\geq 0$ and $n$ is an integer. For $B=0$ or $n=0$ one recovers an isotropic stellar wind, while for $B>0$ and $n>0$, the density profile is anisotropic. This function is normalized such that the integral of the mass flux around the star recovers the total mass loss rate, $\dot{M}_{w}=2\pi\int^{\pi}_{0}\rho_{w}v_{w}r^{2}\sin\theta d\theta$.

4 Solution of the equations

4.1 Non dimensional equations

To solve the equations we define the following non dimensional variables: the non dimensional radius

[TABLE]

the non dimensional time

[TABLE]

the non dimensional mass flux

[TABLE]

and the non dimensional momentum fluxes

[TABLE]

The non dimensional velocities of the accretion flow and the stellar wind are

[TABLE]

and

[TABLE]

Also, we define the ratio of the wind mass loss rate and the mass accretion rate

[TABLE]

and the ratio between the stellar wind and the accretion flow momentum rates

[TABLE]

Finally, the non dimensional densities of the accretion flow and the stellar wind are given by

[TABLE]

and

[TABLE]

In terms of the new variables, eqns. (3) - (7) can be written as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $f(\theta)$ is defined in eq. (16).

These equations need initial conditions to advance in time, and boundary conditions (BCs) at the pole ($\theta=0$) and at the equator ($\theta=\pi/2$).

4.2 Boundary conditions

In the next section we expand the variables in powers of $\theta$ and obtain equations for their time evolution at the pole and at the equator. These solutions provide BCs for the partial differential equations (31) - (35).

4.2.1 Expansions around the pole

We expand in power series the mass and momentum fluxes as well as the radius of the shell. The equations are expanded to second order in $\theta$ for $\theta\ll 1$, such that the variables are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where the coefficients $b_{m1}$, $b_{r1}$, $b_{\theta 2}$, $b_{\phi 2}$, and $r_{s0}$ are functions of the non dimensional time ($\tau$) given by the solution of the set of differential equations described in Appendix B.

4.2.2 Expansions around the equator

In the equator the density of the accretion flow at centrifugal radius diverges. If the shell evolves in this direction, eventually, it is going to find a barrier of infinite density. At that point the shell will stagnate at $R_{cen}$ or collapse back to stellar surface.

We expand the variables around the equator as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where $\Theta=\left(\frac{\pi}{2}-\eta\right)-\theta\ll 1$, and the angle $\eta$ defines a physical boundary in the equatorial region (e.g., a disk).

The coefficients $q_{m0}$, $q_{m1}$, $q_{r0}$, $q_{r1}$, $q_{\theta 0}$, $q_{\theta 1}$, $q_{\phi 0}$, $q_{\phi 1}$, $q_{rs0}$, and $q_{rs1}$ are functions of the non dimensional time ($\tau$) given by the solution of a set of differential equations described in Appendix C.

5 Results

Eqns. (31) - (35) describe the evolution of the shell. The non dimensional equations are solved numerically. We assume that initially, the shell is spherical and massless, with a radius close to the stellar surface. We assume an initial non dimensional radius of $r_{s0}(0)\simeq R_{*}/R_{cen}\simeq 10^{-4}$. We assume that the ratio of the wind mass loss rate and the mass accretion rate is $\alpha=0.1$, a typical value for molecular outflows (e.g., see figure 14 of Ellerbroek et al. 2013). Also, we assume $\beta=21$333The stellar wind and the free fall velocities correspond to parameters of the central star of the molecular outflow CB 26 (Launhardt et al., 2009). These parameters are a stellar mass $\mathrm{M}_{*}=0.5\mathrm{M}_{\odot}$ and a centrifugal radius of $R_{cen}=200$ AU (Launhardt & Sargent, 2001). We also assumed a radius $R_{*}=2$ R⊙. The stellar wind velocity is taken as the escape velocity of the central star. With these assumptions, the stellar wind velocity of $309$ km s*-1* and free fall velocity of $v_{0}=1.5$ km s*-1*. . The integration is done from $t=0$ to $t=1000$ yr.

In this section we study the shell evolution for different stellar wind models. As an example, consider the molecular outflow produced by an anisotropic stellar wind with $A=1$, $B=20$, and $n=2$. Figure 2 shows the shape of the shell $R_{s}(\theta,t)$ for different times from $t=250$ yr to $t=1000$ yr. The shells are elongated along the cloud rotational axis. We define the shell collimation as the ratio

[TABLE]

where $R_{s}(0,t)$ is the shell radius at the pole, and $\varpi_{\rm max}$ is the maximum width of the shell. This ratio measures the shell elongation. During the shell evolution, this model has a collimation $C\sim 2.5$, similar to the observed outflows CB 26 (Launhardt et al., 2009) and DG Tau B (Zapata et al., 2015).

Figure 3 shows the radial, the $\theta$, and the azimuthal velocities as functions of $\theta$ for this model, for different times. The velocities of the shell are obtained from eqns. (2). The left panel shows the radial velocity of the shell. This velocity decreases with time, i.e., the shell is slowing down. It also decreases with the angle such that at the equator, $\theta=\pi/2$ this velocity tends zero, because the shell finds a barrier at $R_{\rm cen}$ where the density is infinite. The middle panel shows the $\theta$-velocity of the shell. This velocity decreases with the time, but increases with the polar angle $\theta$, due to the material that slides from the pole to the equator; this material feeds the accretion disk $U_{\theta}>0$. The right panel shows the azimuthal velocity. This rotation velocity decreases with the time and increases with angle: at the equator the rotation velocity is maximum because the orbits of the accretion flow that lands at this point have the largest angular momentum with respect to the pole.

The mass surface density of the shell along the radial direction, obtained from eq. (1), is plotted in Figure 4. In this figure, we observe that, for angles close to the pole, the surface density decreases with time, while for angles close to the equator, the surface density increases with time.

The total mass of the shell is given by

[TABLE]

where $dA=2\pi R^{2}\sin\theta d\theta$, and $P_{m}$ defined in eq. (1). Figure 5 shows that the mass increases with time.

The specific angular momentum of the shell in $z-$direction is

[TABLE]

and the total angular momentum is

[TABLE]

Figure 6 shows the total angular momentum, which, increases with time.

In order to compare the outflow model with observations, we projected the velocity field of the shell along the line of sight for an inclination angle $i=5^{\circ}$ with respect to the plane of the sky. The velocity of the line of sight $v_{los}$ is shown in Figure 7. The left and right panels show the velocity at 250 yr and 1000 yr, respectively. For an inclination angle larger than $0^{\circ}$, the velocity along the line of sight is a combination of the radial and the $\theta$ velocities. Figure 8 shows cuts at different heights $z_{cut}$ of the map at 250 yr. The left panel shows a position-velocity diagram for $z_{cut}=-560$ AU. The middle panel and right panel have cuts at $z_{cut}=0$ and $z_{cut}=420$, respectively.

Now we consider the effect of degree of anisotropy of the stellar wind on the shape of the shell. Figure 9 shows the shape of the shell for parameters $\alpha=0.1$, $\beta=21$, $r_{s0}(0)=10^{-4}$, and a stellar wind model with $A=1$, a time of $t=1000$ yr, for different values of the anisotropy parameters $B$ and $n$. The left panel shows the shape of the shells for $n=2$, one can see that as $B$ increases, the shell becomes more elongated. The right panel shows the shape shell for $B=20$, as the exponent $n$ increases, the shell are more elongated.

Figure 10 shows the collimation of the shell for parameters $\alpha=0.1$, $\beta=21$, $r_{s0}(0)=10^{-4}$, and different values of the anisotropy parameter $B$ and different values of the exponent $n$ from an isotropic stellar wind ($B=0$), to very anisotropic stellar winds ($B=20$ and $n=4$). As expected, the collimation increases with both $B$ and $n$.

6 Discussion

We find that at the pole the shell collapses for isotropic and anisotropic stellar winds if the value of the ratio between the stellar wind and the accretion flow momentum rates is less than a critical value, $\beta<\beta_{\rm crit}$ (see appendix B). This happens because the shell does not have enough momentum to escape. At the equator, for a given value of $\beta$, the shell will always stagnate near the centrifugal radius (see appendix C). This happens because the density diverges at $R_{\rm cen}$.

We have also compared the results of our model with those of Wilkin & Stahler (2003). For an isotropic wind and the same values of the parameters $\alpha$ and $\beta$, we find that our shells have the same sizes at the pole. Also, the collimation factor is the same as in their model $C\sim 1.6$ (see appendix D). The only differences are due to the assumption of different BCs at the disk surface, close to the equator.

The collimation of the shells depends on the anisotropy of the stellar wind and the accretion flow. In the models with an anisotropic stellar wind and the Ulrich accretion flow, it is difficult to get collimation factors $C$ much larger than 3 (figure 10), while the collimation factors of observed sources have values $\sim$ 3 - 10, Bontemps et al. (1996).

We also compare our model with the molecular outflow CB 26 (Launhardt et al., 2009), at the time when they both have the same size in the polar direction. The result of this comparison is:

The dynamical time of the model is $t=$ 250 yr which is half of the kinematic age calculated with the observed size and current velocity. This discrepancy is due to the fact that the shell is decelerating. 2. 2.

The radial velocity of the model, is of the order to 10 km s*-1*, consistent with the observed expansion velocity (e.g. Lee et al. 2018). 3. 3.

The collimation factor is similar to the observed value. 4. 4.

The shell mass in the model is 2$\times$10*-3* M⊙, twice the observed value. 5. 5.

The rotation velocity of the model is lower by an order of magnitude than the observed value (see Table 1). 6. 6.

The total angular momentum is also lower than the observed value.

The low rotation of the model may be resolved, if that the stellar wind has angular momentum, or the parent cloud has more angular momentum than the Ulrich’s flow, or with a combination of both mechanisms. Part of the problem is that the accreting envelope does not have large rotation velocities. In addition, the model shell has more mass than the observed shell, with slowly rotating material.

Finally, we note that in our thin shell model, we assume that the pressure effects are negligible. Pressure gradients will not affect the thin shell approximation which depends on an efficient cooling of the shocked gas. On the other hand, pressure gradients inside the shell could change the gas tangential dynamics, accelerating or decelerating the flow along the shell. This effect is not expected to be important when the flow is supersonic. Thus, to evaluate the effect of the pressure, one has to calculate the temperature of the shell, which is out of the scope of this paper.

7 Conclusions

To understand the evolution and properties of molecular outflows we developed a model of the interaction between a stellar wind and an accretion flow that follows the evolution of a thin shell that is pushed by and entrains material from both flows. We have formulated the problem in such a way that we can consider an accretion flow and a stellar wind with general velocity fields (collimated or not collimated and with or without rotation) provided that they have axial symmetry.

In the present paper we have considered isotropic and anisotropic stellar winds with only radial velocity. The accretion flow is given by the collapse of a slowly rotating molecular cloud from Ulrich (1976). The evolution of the outflow was followed from its origin, close to the stellar surface, to large distances from the central star.

The shell evolution has a strong dependence on the ratio between the wind and the accretion flow mass and momentum rates, $\alpha$ and $\beta$ (eqns. [27] and [28]). In order for the molecular outflow shells to expand, it is necessary that $\beta\geq\beta_{crit}$ for a given value of $\alpha$. If $\beta<\beta_{crit}$, the whole shell collapses back to stellar surface.

The interacting flows considered in this work, produce moderate outflow collimation ($C\sim 3$) and low rotation ($v_{\phi}\sim 0.1$ km s*-1*). These values are lower than observed in the sources in Table 1. It is left as future work to explore other physical collapsing envelopes and stellar winds to determine the outflow characteristics.

Acknowledgements. J. A. López-Vázquez and S. Lizano acknowledge support from PAPIIT-UNAM IN101418 and CONACyT 23863. J. Cantó acknowledges support from PAPIIT-UNAM-IG 100218. We thank an anonymous referee for useful suggestions that improved the presentation of this paper.

Appendix A Derivation of the equations

In spherical coordinates the continuity equation is given by

[TABLE]

where $\rho$ is the mass volume density and $v_{r}$, $v_{\theta}$, and $v_{\phi}$ are the fluid velocities. We assume axisimmetry with respect to $\phi$ direction and multiply the above equation by $r^{2}\sin\theta$. Then, the continuity equation can be written as

[TABLE]

Integrating the latter equation in the radial direction from $R$ to $R+\Delta R$ for fixed $\theta$ (see Figure 1), the continuity equation of the shell can be written as

[TABLE]

In this equation the mass surface density of the shell in the radial direction is $\sigma=\rho_{s}\delta$, and $U_{r}$ and $U_{\theta}$ are the radial and $\theta$ velocity components of the shell material, respectively.

The equation of the fluid momentum in the radial direction is given by

[TABLE]

where $F_{g}=-GM_{*}\rho/r^{2}$ is the gravitational force per unit volume.

The equations of the fluid momentum in the $\theta$ and the azimuthal directions are given by

[TABLE]

and

[TABLE]

One multiplies eq. (A1) by $v_{r}$, $v_{\theta}$, and $v_{\phi}$, and add the result to eqns. (A4) - (A6), respectively. Then, one integrates each equation in the radial direction for a thin shell, considering the axial symmetry, and obtains the following equations for the momenta of the shell in the radial, the $\theta$ and the azimuthal directions

[TABLE]

[TABLE]

and

[TABLE]

where $U_{\phi}$ is the azimuthal velocity of the shell material.

To obtain the evolution of the shell radius, one can write

[TABLE]

Then, the evolution of the radius of the shell is given by

[TABLE]

Eqns. (A3), (A7), (A8), (A9), and (A11) describe the evolution of the shell formed by the shock between a stellar wind and an accretion flow with a given velocity field. These equations can be written the more compact form shown in the eqns. (3) - (7) in terms of the mass and momentum fluxes (eqns. [1] and [2]).

Appendix B Radius and velocities around the pole

Around the pole, the mass and the momentum fluxes, and the radius can be expanded as power series of $\theta$ to second order as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In the above equations, we can note that the mass and the momentum fluxes do not have component of order zero, since these momenta are zero at $\theta=0$.

Substituting the eq. (40) in eq. (14) and expanding in Taylor’s series for $\theta\ll 1$ and $\theta_{0}\ll 1$, one finds the relation between the radius at the pole $r_{s0}$, the angle $\theta$, and the angle $\theta_{0}$,

[TABLE]

In addition, the velocities and the density of the accretion flow are expanded in Taylor’s series to first order in $\theta$. Using the latter equation for $\theta_{0}$ one obtains

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Also, the expansion for the wind density gives

[TABLE]

Substituting the eqns. (B1) - (B11) in eqns. (31) - (35), one finds that the coefficients $b_{m2}(\tau)=b_{r2}(\tau)=b_{\theta 1}(\tau)=b_{\phi 1}(\tau)=r_{s1}(\tau)=r_{s2}(\tau)=0$. In addition, one obtains a set of ordinary differential equations for the functions $b_{m1}(\tau)$, $b_{r1}(\tau)$, $b_{\theta 2}(\tau)$, $b_{\phi 2}(\tau)$, and $r_{s0}(\tau)$. These equations are,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

The coefficients $b_{m1}(\tau)$, $b_{r1}(\tau)$, $b_{\theta 2}(\tau)$, $b_{\phi 2}(\tau)$, and $r_{s0}(\tau)$ are functions of the time ($\tau$). The shell starts at an initial radius $r_{s0}(0)$, close to the stellar radius. Because initially the shell is massless, $b_{m1}(0)=b_{r1}(0)=b_{\theta 2}(0)=b_{\phi 2}(0)=0$. Thus, to find the ratios of $b_{r1}/b_{m1}$, $b_{\theta 2}/b_{m1}$, and $b_{\phi 2}/b_{m1}$, we expand the coefficients to first order in $\tau$ for $\tau\ll 1$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Substituting these equations in eqns. (B12) - (B16), one obtains $c_{m}$, $c_{r}$, $c_{\theta}$, $c_{\phi}$, and $c_{rs}$ as function of the initial radius $r_{s0}(0)$ and the ratio

[TABLE]

Then,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $Q_{1}$ is given by,

[TABLE]

One also obtains a quadratic equation for $\lambda$ substituting eqns. (B23) and (B24) in eq. (B22),

[TABLE]

Eqns. (B12) - (B16) are solved numerically for the coefficients $b_{m1}$, $b_{r1}$, $b_{\theta 2}$, $b_{\phi 2}$, and $r_{s0}$. The initial conditions are obtained from eqns. (B17) - (B21). The latter are evaluated at a very small non dimensional time, $\tau=10^{-9}$. From the solution of these equations one finds BCs at the pole for the mass and momentum fluxes, and the shell radius as a function of time.

For the case of the isotropic stellar wind ($B=0$ in eq. [16]), we explore the critical value of $\beta$ required for the expansion of the shell. Table 2 shows the critical value $\beta_{crit}$, for different values of $\alpha$ and initial radius $r_{s0}(0)$. Figure 11 shows the shell radius at the pole for early times for a model with $\alpha=0.1$, initial radius $r_{s0}(0)=10^{-4}$, and two cases: $\beta=\beta_{crit}$ and $\beta<\beta_{crit}$. In the former case, the shell always expands; in the latter case the shell collapses back onto the stellar surface.

To obtain $\beta_{crit}$, one takes into account the weight of the shell, the change of the radial momentum in the $\theta$ direction, and the momentum added to the shell by both the stellar wind and the accretion flow. If one neglects the weight of the shell for the models in Table 2, one obtains that a smaller value $\beta>0.5$ is enough for the shell to expand.

Assuming $\alpha=0.1$ and $\beta=21$, Figure 12 shows the evolution of the shell radius for models with parameters $A=1$ and $n=2$, and different values of the anisotropy parameter $B$ (left panel), and models with parameters $A=1$ and $B=20$, and different values of the exponent $n$ (right panel). This figure shows that the shell radius increases with the anisotropy parameter $B$ and the exponent $n$. For the same models, Figure 13 shows the radial velocity at the pole (the shell expansion velocity). The radial velocity also increases with the anisotropy parameter $B$ and the exponent $n$, and it tends to a constant value at large times. The $\theta$ and the azimuthal velocities at the pole are zero because $v_{\theta}=(b_{\theta 2}/b_{m1})\theta$ and $v_{\phi}=(b_{\phi 2}/b_{m1})\theta$ are linear functions of $\theta$.

Appendix C Radius and velocities around the equator

Around the equator the mass and the momentum fluxes, and the radius are expanded as power series around the angle $\left(\frac{\pi}{2}-\eta\right)-\theta$. These expansions are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where $\left(\frac{\pi}{2}-\eta\right)-\theta\ll 1$.

The expansion in Taylor’s series to first order around the angle $\left(\frac{\pi}{2}-\eta\right)-\theta$ of the velocity field and the density of the accretion flow gives

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where $\theta_{0}$ is obtained of the eq. (14).

The wind density is

[TABLE]

Substituting eqns. (C1) - (C5), the velocity field and the density of the accretion flow, and the density of the stellar wind given in eqns. (C6) - (C10), in eqns. (31) - (35), one obtains a set of differential equations for the coefficients

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where the functions $f_{ai,j}$ are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In these equations the functions $\Gamma_{i}$ are defined as,

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In the above equation, the partial derivatives of $\cos\theta_{0}$ are obtained writing eq. (14) as

[TABLE]

where $\zeta=1/r_{s}$, and are given by

[TABLE]

[TABLE]

The functions $P_{i}$ are given by

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Finally, the functions $T_{i}$ can be written as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The coefficients $q_{m0}(\tau)$, $q_{m1}(\tau)$, $q_{r0}(\tau)$, $q_{r1}(\tau)$, $q_{\theta 0}(\tau)$, $q_{\theta 1}(\tau)$, $q_{\phi 0}(\tau)$, $q_{\phi 1}(\tau)$, $q_{rs0}(\tau)$, and $q_{rs1}(\tau)$ are functions of the time ($\tau$). The shell starts at an equatorial radius $q_{rs0}(0)=r_{s0}(0)$ close to the star. Initially the shell is massless, so $q_{m0}(0)=q_{m1}(0)=q_{r0}(0)=q_{r1}(0)=q_{\theta 0}(0)=q_{\theta 1}(0)=q_{\phi 0}(0)=q_{\phi 1}(0)=q_{rs1}(0)=0$. Thus, to obtain the ratios between the coefficients $q_{ai}$, one expands them for early times $\tau\ll 1$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Substituting these equations in eqns. (C11) - (C20), one obtains

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

The functions $Q_{1}$, $Q_{2}$, and $\gamma_{i}$ are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Substituting eqns. (C66) and (C67) in eq. (C76), one obtains a quadratic function for $\Lambda$

[TABLE]

Solving numerically eqns. (C11) - (C20) with initial conditions given by eqns. (C56) -(C65) for small $\tau$ ($\tau=10^{-9}$), one finds the values of the coefficients $q_{m0}$, $q_{m1}$, $q_{r0}$, $q_{r1}$, $q_{\theta 0}$, $q_{\theta 1}$, $q_{\phi 0}$, $q_{\phi 1}$, $q_{rs0}$, and $q_{rs1}$. Thus, one obtains BCs at the equator for the mass and the momentum fluxes, and the radius as a function of time.

Assuming $\alpha=0.1$, $\beta=21$, and $r_{s0}(0)=10^{-4}$, Figure 14 shows the evolution of the normalized shell radius for models with parameters $A=1$ and $n=2$, and different values of the anisotropy parameter $B$ (left panel), and models with parameters $A=1$ and $B=20$, and different values of exponent $n$ (right panel). This figure shows that the normalized shell radius grows to a stagnation point at the centrifugal radius. For the same models, Figure 15 shows from left to right, the radial, the $\theta$, and the azimuthal velocities at the equator. With small oscillations, the radial and $\theta$ velocities tend to a constant value, while the azimuthal velocity decreases with time.

Appendix D Models of Wilkin $\&$ Stahler

We compare our models with the outflow models of Wilkin & Stahler (2003) for an isotropic stellar wind ($B=0$). The parameters of the model in their Figure 5 are: a ratio of the wind mass loss rate and the accretion rate $\alpha=1/3$, a wind velocity $v_{w}=159{\rm\,km\,s^{-1}}$, a stellar radius $R_{*}=3\,R_{\odot}$, an angular speed of the rotating envelope $\Omega=2\times 10^{-14}{\rm\,s^{-1}}$, and a sound speed $a_{0}=0.2{\rm\,km\,s^{-1}}$.

We recover the parameters we use in our model in the following way: from their eq. [1], the mass accretion rate is $\dot{M}_{a}=1.85\times 10^{-6}\,M_{\odot}{\rm\,yr^{-1}}$. Thus, the stellar mass, given by their eq. [9], is $M_{*}=0.07\,M_{\odot}$. From their eq. [2], the angular and sound speeds give a centrifugal radius $R_{cen}=0.055$ AU. The ratio of the wind and the accretion momentum rates $\beta$ is given by $\beta=\alpha v_{w}/v_{0}=1.57$, where $v_{0}=(GM_{*}/R_{\rm cen})^{1/2}=33.8{\rm\,km\,s^{-1}}$. Finally, the non dimensional time (eq. [18]) is $\Delta\tau=2.03(\Delta t/0.016{\rm\,yr})$. We also assume a disk as a boundary condition in the equatorial region with an angle $\eta=5^{\circ}$.

The results for this integration are plotted in figure 16. The shell sizes at the pole are the same as the models in Figure (5) of Wilkin & Stahler. Nevertheless, the boundary conditions of both models at the disk surface are different. In our model, the BC is that the shell radius on the disk surface cannot be larger than the centrifugal radius $R_{\rm cen}$. In their model, the shell expands continuously beyond the centrifugal radius.