# Magnetized Protostellar Bipolar Outflows

## 1 Introduction

General outflows and collimated jets are thought to be intimately related to infall and protostellar accretion. Recent observations (Donati et al. donal97 (1997), Guenther & Emerson gue96 (1996)) suggest the relevance of magnetic fields to star formation, first suggested by Mestel & Spitzer (ms56 (1956)). Tomisaka (tomisaka (1998)) has studied numerically the dynamical collapse of magnetized molecular cloud cores from the runaway cloud collapse phase to the central point mass accretion phase. He finds that the evolution of the cloud contracting under its self-gravity is well expressed by a self-similar solution. Moreover inflow-outflow circulation appears as a natural consequence of the initial configuration. His results support the magnetized self-similar models as recently presented by Fiege and Henriksen (fiege1 (1996) a,b, hereafter FH1, FH2) following Henriksen (henrik (1993), henrik2 (1994)), and Henriksen and Valls-Gabaud (HVG (1994), hereafter HV). In those models, the self-similar gravitationally driven convective circulation in a heated quadrupolar protostar envelope is solved rigorously. Moreover recent observations (Greaves & Holland GH (1998)) of the polarized dust emission from six star-forming cloud cores have revealed, for the first time, some very large twists in the magnetic field. This is consistent with FH1 and HV where circumstellar magnetic fields follow large-scale streamlines. This remains very nearly the case in the present model.

The self-similar models regard the molecular outflow as a natural consequence of the circulation established by the collapse of the pre-stellar cloud. The models describe an outflow velocity that increases toward the axis of rotation, a convective pattern of infall/outflow, self-consistent axial collimating magnetic fields and rotation, and “cored apple” type distributions of circum-protostellar gas (Andre et al. andre (1993)).

Such models may imply a natural connection between the fast ionized jets seen near the polar axis of the wind, and the slower and less-collimated molecular outflows that surround them. However velocities obtained by FH1 and HV are smaller than some observed velocities (e.g. EHV outflow ($`V100kms^1`$), Bachiller et al. bachiller (1990)) unless the solution is pushed very close to the central mass where the limiting speed is about the escape speed (and the jet speed) of a few hundred $`kms^1`$. But this is precisely the region where the inner boundary layer may be dominant and the assumptions of this type of model are liable to break down. Moreover the luminosity needed for radiative heating in FH1 is too large in order to drive high-velocity outflows if only dust opacity is used. In the present work we show that the Poynting flux increases both the velocity and collimation of the modeled outflows by helping to transport mass and energy from the equatorial regions. This is much as has been argued for some time by other authors (e.g. Pudritz and Norman PN (1986), Ouyed and Pudritz ouyedp (1997)), but our models are globally consistent in space. Self-similar models can not of course be globally consistent in time even if they are not stationary (we have found such models and will report them elsewhere), since they are ignorant of initial conditions. And in fact they must also be ‘intermediate’ in space since there are inevitably small regions near the equator, the axis and at the centre which are excluded from the domain of self-similarity. Thus, we derive globally self-consistent models of bipolar outflows and infall accretion, but specifically exclude the regions dominated by the disc, jet, and protostar.

Qualitatively we may understand the quadrupolar global circulation in the following fashion. As material falls towards the central object, it is gradually slowed down by increasing radial pressure gradients due to the rising temperature, density, and magnetic field encountered by material near the central object. This pressure barrier, with the help of centrifugal forces, deflects and accelerates much of the infalling matter into an outflow along the axis of symmetry. The outflow velocity is naturally the highest near the axis of rotation since the pressure gradients are strongest there. Mathematically the quadrupolar model can be seen as a form of instability of the spherical Bondi accretion flow, wherein rotation, magnetic fields and anisotropic heating transform a central nodal singularity into a saddle point singularity.

The pressure gradients also act to collimate the flow, forming consistently the now traditional de Laval nozzle (e.g. Blandford and Rees br (1974), Königl ko1 (1982)). We find however that the most convincing flows are everywhere super-Alfvénic, in agreement with FH1, but the nozzle may still be critical in terms of the fast magnetosonic speed as discussed in FH1 and below. We note that one can expect conical shocks near the rotational axis, which would seem to be necessary to explain the shocked molecular hydrogen emission often observed.

The paper is organized as follows: in §2 we introduce the model and its approximations, and discuss their consequences; numerical solutions are presented in §3 for virial-isothermal and radiative models describing outflows from low and high mass protostars. In §4 we study the ensuing synthetic radio maps; this section is followed by discussions of the results in §5, and conclusions are presented in §6.

## 2 The Analytical Model

The main difference between hydrodynamic and MHD outflows is the low asymptotic Mach number in the MHD case. For protostellar outflows, the magnetosonic speed is about $`100kms^1`$ at $`10^5`$ AU and therefore far exceeds the sonic speed. It is not a priori obvious that the internal Alfvénic Mach number of the jet or outflow is small. However most of our models are super-Alfvénic everywhere.

The plasma is supposed to be perfectly conducting and therefore the flow is governed by the basic equations of ideal MHD without taking into account resistivity or ambipolar diffusion. In order to make the system tractable we assume axisymmetric flow so that $`/\varphi =0`$ and all flow variables are functions only of $`r`$ and $`\theta `$.

Kudoh and Shibata (1997b ) have performed time dependent one-dimensional (1.5-dimensional) magnetohydrodynamic numerical simulations of astrophysical jets that are magnetically driven from Keplerian discs. They have found that the jets, which are ejected from the disc, have the same properties as the steady magnetically driven jets they had found before (Kudoh & Shibata 1997a ). Their numerical results suggest that a steady model, as we are assuming in the present work (i.e. $`/t=0`$), is a good first approximation to the time-dependent problem on time-scales short compared to that required to dissipate the circumstellar material. Time independence will be relaxed in future work.

Accretion discs near protostars are probably heavily convective and therefore prone to dynamo action (Brandenburg et al. brand (1995), FH1). Since the ensuing magnetic fields are loaded with disc plasma, currents are allowed to circulate between the disc surface and the wind region above the disc. For rapidly rotating discs, this dynamical system may evolve naturally into a quadrupolar structure (Camenzind cam (1990), Khanna and Camenzind kc (1996)). The strong differential rotation of accretion discs is responsible for the excitation of quadrupolar modes in discs. Consequently, in the present model, magnetic field and streamlines are required to be quadrupolar in the poloidal plane. A moderate amount of heating supplied near the turning point allows the outflow to possess finite velocities at infinity.

These considerations have led us to use the set of equations of steady, axisymmetric, ideal MHD and to seek solutions with a quadrupolar geometry. Even though the poloidal components of the magnetic and velocity fields have to be parallel, the toroidal components need not share the same constant of proportionality (Henriksen henrik3 (1996)). This permits a poloidal conservative electric field to exist in the inertial frame, and so admits steady Poynting flux driving. One needs then to introduce an electric potential function and an azimuthal magnetic field independent of the azimuthal velocity. This point is the major improvement of the model with respect to FH1. In the radiative case, the equations describing diffusive radiative transfer (to zero order in $`v/c`$) and a Kramer’s-type law for the Rosseland opacity are added as in FH1. We also neglect the self-gravity of the protostellar material by assuming that the gravitational field is dominated by the central mass that is assumed to be an external fixed parameter in our model. Henceforth the model only applies to a state where the protostellar core is already formed.

It is unlikely that this model can be extended to include the optical jets since these must originate very near to the central object where the model may not apply in its steady form. Nevertheless, the inclusion of a jet model in the axial region, such as given by lery et al.(lery1 (1998), 1999b ), could remedy this lack in order to make a more global model for protostellar objects.

### 2.1 Scales and Self-similar Laws

Our model is developed within the context of power law self-similar models as developed by Henriksen (henrik (1993)), and in HV and FH1. A model with such symmetry was already used by Bardeen & Berger (BB (1978)) for galactic winds, although the present development has been independent. Parallel developments have been largely occupied with the study of winds from an established accretion disc (Blandford & Payne bp (1982), Konigl ko (1989), Ferreira & Pelletier ferr1 (1993), Rosso & Pelletier RossoP (1996), Ferreira ferr2 (1997), Contopoulos & Lovelace conto (1994), Contopoulos contop (1995), Sauty & Tsinganos sautyt (1994), Tsinganos et al. tsinetal (1996), Vlahakis & Tsinganos VT (1998)). Our use of this model in order to study inflow and outflow as part of the global circulation dynamics around protostars seems unique. The form used corresponds to an example of ‘incomplete self-similarity’ in the classification of Barenblatt (baren (1996)), but fits into the general scheme of self-similarity advocated by Carter and Henriksen (CH (1991)) with $`r`$ as the direction of the Lie self-similar motion. We work in spherical polar coordinates $`r,\theta ,\varphi `$ centered on the point mass $`M`$ and having the polar axis directed along the mean angular momentum of the surrounding material.

The self-similar symmetry is identical to that used in FH1 except that two quantities $`y_\varphi `$ and $`y_p`$ defined such that $`\sqrt{4\pi \mu }y_\beta v_\beta /(B_\beta /\sqrt{4\pi \rho })`$ (where $`\beta `$ indicates the appropriate component) replace $`y`$ in FH1. The power laws of the self-similar symmetry are determined, up to a single parameter $`\alpha `$, if we assume that the local gravitational field is dominated by a fixed central mass. In terms of a fiducial radial distance, $`r_o`$, the self-similar symmetry is sought as a function of two scale invariants, $`r/r_o`$ and $`\theta `$, in a separated power-law form. The equations of radiative MHD require the following radial scaling relations for the variables:

$$𝐯=\sqrt{\frac{GM}{r_o}}\left(\frac{r}{r_o}\right)^{1/2}𝐮(\theta ),$$
(1)   
$$B_\varphi =\sqrt{\frac{GM^2}{r_o^4}}\left(\frac{r}{r_o}\right)^{\alpha 3/4}\frac{u_\varphi (\theta )}{y_\varphi (\theta )},$$
(2)   
$$B_p=\sqrt{\frac{GM^2}{r_o^4}}\left(\frac{r}{r_o}\right)^{\alpha 3/4}\frac{u_p(\theta )}{y_p(\theta )},$$
(3)   
$$\rho =\frac{M}{r_o^3}\left(\frac{r}{r_o}\right)^{2\alpha 1/2}\mu (\theta ),$$
(4)   
$$p=\frac{GM^2}{r_o^4}\left(\frac{r}{r_o}\right)^{2\alpha 3/2}P(\theta ),$$
(5)   
$$\frac{kT}{\overline{\mu }m_H}=\frac{GM}{r_o}\left(\frac{r}{r_o}\right)^1\mathrm{\Theta }(\theta ),$$
(6)   
$$F_{rad}=\left(\frac{GM}{r_o}\right)^{3/2}\frac{M}{r_o^3}\left(\frac{r}{r_o}\right)^{\alpha _f2}f(\theta ).$$
(7)   
In these equations the microscopic constants are represented by $`k`$ for Boltzmann’s constant, $`\overline{\mu }`$ for the mean atomic weight, $`m_H`$ for the mass of the hydrogen atom. The self-similar index $`\alpha `$ is imposed as a parameter of the solution, but as in FH1 it is constrained to lie in the range $`1/4\alpha >1/2`$. In the last equation, the index $`\alpha _f`$ is a measure of the loss (if negative) or gain in radiation energy as a function of radial distance.

Whenever we take account of radiative diffusion explicitly, as well as for the definition of the luminosity, we will use the same formulation as in FH1. The simplest approach however is to assume that the temperature is some fraction of the ‘virialized’ value. That is, according to relation (6), that $`\mathrm{\Theta }`$ takes some constant value $`1`$. By using the first law and assuming an ideal gas (so that with the above relations for $`p`$ and $`\rho `$ we have $`P(\theta )=\mu \mathrm{\Theta }`$ ) this assumption implies that on a sphere the specific entropy $`s`$ varies according to

$$Tds=\frac{kT}{\overline{\mu }m_H}d\mathrm{ln}\mu .$$
(8)   
Consequently we are explicitly adding or subtracting heat from the system on a spherical surface according to the sign of $`d\mathrm{ln}\mu `$. In our models this normally increases towards the rotational axis. Thus, at least implicitly, we always employ some form of heat transfer, presumably radiative, to the gas over the poles relative to the equatorial material.

The same procedure for variations in the radial direction yields

$$Tds=\frac{GM}{r_o}\mathrm{\Theta }\left[\frac{1+2(\gamma 1)(\alpha 1/4)}{\gamma 1}\right]\left(\frac{r_o}{r}\right)^2\frac{dr}{r_o}.$$
(9)   
Energy is therefore lost to the material with increasing radius provided that the quantity in the square bracket is positive. This requires $`\alpha 1/4>\frac{1}{2(\gamma 1)}`$, which is true for reasonable ratios of specific heats throughout our permitted range of $`\alpha `$.

In order for radiative heating to be adequate we must have $`\tau F_{rad}/c`$ sufficiently large, where $`F_{rad}`$ is the net radiative flux through a mass element. We note that because of the dependence on $`F_{rad}`$, this is not the same as requiring an optical depth of order unity. Nevertheless in order for the radiative transfer by diffusion to be plausible, this latter condition must also be satisfied.

In either case the issue depends on the nature of the opacity $`\kappa `$. We follow FH1 in supposing that this is predominantly a dust opacity (at least at some distance from the ‘star’) which is taken in Kramer’s form

$$\kappa =\kappa _o\rho ^aT^b.$$
(10)   
A fit to a dust model in FH1 yields $`a=0`$, and $`b2`$ and $`\kappa _o=6.9\times 10^6cgs`$. We can calculate the mean optical depth between an inner radius $`r_{}`$ and the fiducial radius $`r_o`$, namely $`<\tau >\frac{1}{4\pi }_r_{}^{r_o}\kappa \rho 𝑑r𝑑\mathrm{\Omega }`$, using this opacity and the self-similar scaling relations. We find

$$<\tau >=\frac{<\mu \mathrm{\Theta }^2>}{|2\alpha 3/2|}\left[\left(r_{}/r_o\right)^{(2\alpha 3/2)}1\right],$$
(11)   
where we have defined the fiducial radius as

$$r_o\left[\kappa _o\left(\frac{\overline{\mu }m_H}{k}\right)^2G^2\right]^{1/4}M^{3/4},$$
(12)   
which is numerically

$$r_o1400(M/M_{})^{3/4}AU.$$
(13)   
Thus we obtain a fiducial radius that is characteristic of the bipolar outflow sources and is a convenient unit for our discussions. A similar scale (to within a factor 5) was previously deduced in Henriksen (henrik (1993)) based on the source luminosity and fundamental constants (including the radiation constant). We note that the external opacity from $`\mathrm{}`$ to $`r_o`$ is given simply by $`<\tau >_{ext}=<\mu \mathrm{\Theta }^2>/|2\alpha 3/2|`$ and the ‘photosphere’ of the cloud is therefore found where $`\tau _{ext}2/3`$.

The corresponding temperature at the fiducial radius is given by

$$T164\mathrm{\Theta }(M/M_{})^{1/4}K.$$
(14)   
We observe that if we identify our radiative fiducial scale $`r_o`$ with an empirical scale $`r_{em}`$ from $`<\mu >/r_{em}^2=2\times 10^{9\pm 1}(AU)^2`$ (Richer et al. PPIV (1998)), we can infer a relation between the mean density parameter $`<\mu >`$ and the central mass in the form

$$<\mu >=3.38\times 10^{3\pm 1}\left(M/M_{}\right)^{3/2}.$$
(15)   
This is a useful way to relate the key model parameter $`\mu `$ to the physical mass and luminosity. We note however that the optical depth is not likely to approach unity for reasonable temperatures, except for exceptionally massive objects. But this does not prevent radiative heating from playing an important dynamical role near the protostar.

In fact using $`r_o`$, as above, and the definition of radiation field $`F_{rad}`$ (Eq. 7), the local energy flow can be expressed by

$$r^2F_{rad}=0.015\left(r/r_o\right)^{\alpha _f}\left(M/M_{}\right)^{5/8}fL_{}.$$
(16)   
The radiative force per unit mass is given by

$$\kappa F_{rad}/c10^{13}\left(M/M_{}\right)^{1/8}\left(r/r_o\right)^{(\alpha _f4)}fcms^2,$$
(17)   
which must be less than the Eddington limit of $`GM/r^2=3.5\times 10^7(M/M_{})\left(r_o/r\right)^2cms^2`$, for consistency. These are helpful considerations for our radiative models below that always verify this condition a posteriori.

### 2.2 The Equations

We use the self-similar forms in the usual set of ideal MHD equations together with the radiative diffusion equation when applicable. The corresponding system of equations is reported in Appendix A, and the ensuing first integrals in Appendix B. As noted above, in the present model, there is no constraint requiring parallelism between velocity and magnetic field as used in FH1. This allows us to discern clearly the effects of a non-zero Poynting flux. The self-similar variable directly related to magnetic field $`y`$ has two components that are not equal in the present model. These quantities are either projected in the poloidal plane or correspond to the toroidal components and are respectively defined by

$$y_{p,\varphi }(\theta )=\frac{M_{ap,a\varphi }}{\sqrt{4\pi \mu (\theta )}}.$$
(18)   
Consequently the system deals with two different components of the Alfvénic Mach number $`M_{ap}`$ and $`M_{a\varphi }`$ defined by:

$$M_{ap,a\varphi }^2(\theta )=\frac{v_{p,\varphi }^2}{B_{p,\varphi }^2/4\pi \rho }.$$
(19)   
In this model, since $`r`$ and $`\varphi `$ respectively correspond to the directions of self-similarity and axisymmetry, only waves that propagate along $`\theta `$ in the poloidal plane can preserve these two symmetries. Therefore the relevant Alfvén mode propagates in the direction $`\theta `$ in the poloidal plane with a phase velocity:

$$V_{a\theta }=\frac{B_\theta }{\sqrt{4\pi \rho }}.$$
(20)   
Moreover we define an Alfvénic point where the total poloidal flow velocity is equal to this value in magnitude and direction (i.e. $`u_r=0`$), which by (18) is equivalent to:

$$M_{ap}=1=4\pi \mu y_p^2.$$
(21)   
One can then define a critical scaled density corresponding to the Alfvénic scaled density

$$\mu _a=\frac{1}{4\pi y_p^2}.$$
(22)   
The flow is super-Alfvénic if the density $`\mu `$ is larger than the critical Alfvén density $`\mu _a`$. Moreover, the compressible slow and fast MHD waves propagate in a poloidal direction $`\theta `$ with phase velocities that satisfy the quartic (Sauty and Tsinganos sautyt (1994))

$$\left(V_{s/f}\right)_\theta ^4\left(V_{s/f}\right)_\theta ^2\left(V_a^2+C_s^2\right)+C_s^2\left(V_a\right)_\theta ^2=0,$$
(23)   
where $`C_s^2`$ is the isothermal sound speed. The condition that $`u_\theta `$ be equal to one of these wave speeds (where we expect a singular point in the family of solutions) can be expressed in terms of self-similar variables as:

$$\left(u_\theta ^2\right)_{crit}=\mathrm{\Theta }+\frac{u_r^2+u_\varphi ^2y_p^2y_\varphi ^2}{4\pi \mu y_p^21},$$
(24)   
and in physical units (cf FH1):

$$\left(v_\theta ^2\right)_{crit}=\frac{P}{\rho }\left(1M_{ap}^2\right)+\frac{B^2}{4\pi \rho }.$$
(25)   
Equation (23) shows that whenever $`V_a>>C_s`$ the slow and fast speeds become $`V_{a\theta }C_s/V_a`$ and $`\sqrt{V_a^2V_{a\theta }^2C_s^2/V_a^2}`$ respectively. Thus a super-Alfvénic flow as defined above can only encounter the fast speed and then only where $`V_aV_{a\theta }`$. In a region where the sound speed dominates, the fast and slow speeds become $`\sqrt{C_s^2V_{a\theta }^2}`$ and $`V_{a\theta }`$ respectively. Either the sound speed or the Alfvén speed is dominant according as $`\mathrm{\Theta }`$ is larger or smaller than $`u_p^2/M_{ap}^2+u_\varphi ^2/M_{a\varphi }^2`$.

As in Henriksen (henrik3 (1996)) we can consider the energy balance in these models. A convenient form for the energy equation is

$$.\left(\rho v(\frac{v^2}{2}+\mathrm{\Phi })+S\right)=v.p,$$
(26)   
where the Poynting flux vector is

$$S=\frac{B^2}{4\pi }v\left(\frac{B.v}{4\pi }\right)B,$$
(27)   
and $`\mathrm{\Phi }=GM/r`$. In self-similar variables the Poynting flux becomes

$$S\left(\frac{GM}{r_o}\right)^{\frac{3}{2}}\left(\frac{M}{r_o^3}\right)\left(\frac{r}{r_o}\right)^{2\alpha 2}\sigma (\theta )$$
(28)   
where the self-similar Poynting flux $`\sigma (\theta )`$ can be defined in terms of the other variables as:

$$\sigma (\theta )=\left(\frac{1}{y_p}\frac{1}{y_\varphi }\right)\left(u_\varphi \frac{u_p^2}{y_p}u_p\frac{u_\varphi ^2}{y_\varphi }\right),$$
(29)   
The energy equation becomes

$`2\alpha \left[\mu u_r({\displaystyle \frac{u^2}{2}}1)+\sigma _r\right]+(2\alpha {\displaystyle \frac{3}{2}})u_rP`$   
$`={\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{d}{d\theta }}\left[\mathrm{sin}\theta (\mu u_\theta ({\displaystyle \frac{u^2}{2}}1)+\sigma _\theta )\right]{\displaystyle \frac{u_\theta }{\mathrm{sin}\theta }}{\displaystyle \frac{dP}{d\theta }}.`$ (30)   
This equation is generally true for our models but it does not provide an independent integral for non-zero pressure. This we treat in the next sub-section.

### 2.3 The Zero Pressure Limit

In the zero pressure limit one can use the self-similar energy equation together with $`\sigma _r=(u_r/u_\theta )\sigma _\theta `$ and $`u_r/u_\theta =d\mathrm{ln}r/d\theta `$ on a stream line to write a Bernoulli integral

$$r^{2\alpha }\mathrm{sin}\theta \left[\mu u_\theta (u^2/21)+\sigma _\theta \right]=B,$$
(31)   
where $`B`$ is the Bernoulli constant on each stream line. In view of Eq. 40 we can write this as

$$\frac{B}{\eta _1}=\frac{1}{r(\theta )}\left[u^2/21+\frac{u_\varphi ^2}{4\pi \mu y_\varphi }\left(\frac{1}{y_\varphi }\frac{1}{y_p}\right)\right],$$
(32)   
where the second term evidently describes the Poynting driving (magnetic energy term) and the first term gives the kinetic and gravitational energies. Hence if $`B0`$, we require the expression contained in brackets on the RHS of this last equation to tend towards infinity both at the equator and on the polar axis since $`r(\theta )\mathrm{}`$ in these limits. This is an improbable outer boundary condition although physically realizable at some large but finite radius (which might be a function of angle and so trace out an excluded region for the self-similar solution). It does permit the exchange between say a large driving Poynting term near the equator and a large kinetic energy term near the axis, but the source of the ‘free’ magnetic energy is left undetermined.

We turn to the special case $`B=0`$, which requires this same factor to be zero everywhere on a streamline. In order that there be an exchange between the kinetic term and the Poynting term in the sum we must have either $`u^2<2`$ and $`y_\varphi <y_p`$ or the reverse of both inequalities. In the first case we do not produce velocities greater than the escape velocity anywhere, while in the second case we require $`M_{ap}<M_{a\varphi }`$ everywhere on the streamline which also seems improbable. We are left then with the first possibility, which includes the case where each term vanishes separately on the boundaries. This suggests that starting from a physically likely free-fall boundary condition, $`v`$ is less than or equal to the escape velocity on every stream line in the absence of pressure, becoming equal again to the escape velocity at large distances.

Thus without pressure effects one is unable to produce net outflows at infinity without appealing to a Poynting flux from an excluded equatorial region. However one does thus deflect and collimate the material (magnetically and centrifugally) in a lossless way, and so all energy added near the star may appear in the outflow at infinity rather than be lost to work against gravity. The fact that there is no net gain of kinetic energy from the magnetic field energy on a stream line is peculiar to the case with $`B=0`$. It is due to the fact that energy that is first gained by the magnetic field at the expense of gravity and rotation during infall is subsequently returned as gravitational potential energy and kinetic energy during the outflow. With $`B0`$, there may be a source of Poynting flux in the excluded region (see e.g. also Henriksen (henrik3 (1996)).

## 3 Numerical Analysis

### 3.1 The Numerical Method

The numerical solutions in this work are obtained by initiating the integrations of the system of six equations from a chosen angle between the axis ($`\theta =0`$) and the equator ($`\theta =\frac{\pi }{2}`$). Starting values are varied until the boundary conditions are met. Solutions can not strictly reach the axis and therefore are limited by a minimum angle $`\theta _{min}`$. Moreover $`\frac{\pi }{2}`$ bounds the solution only if we demand exact reflection symmetry about the equator. Solutions are thus generally defined on a wedge given by $`0<\theta _{min}\theta \theta _{max}`$. In order to be able to start at will from either a sub- or super-Alfvénic point, a systematic search for singular points of the system has been conducted. This allows us to restrict a first guess for the typical values to use as input for the system.

Several general criteria have been employed to further constrain the solutions. First the radial velocity $`u_r`$ must vanish once, to have inflow-outflow, and only once, to avoid a higher order convection pattern. We also require $`u_\theta `$ and $`u_\varphi `$ to vanish at both boundaries, since there must be no mass flux across them, and since only zero angular momentum material can reach the axis (which all material lying strictly in the equator does). In practice, the first of these two conditions is sufficient to impose the second one. The same constraints apply to the magnetic field components. We take $`u_\theta `$ to be negative in order to get a circulation pattern rising over the poles, and $`u_\varphi `$ to be positive as a convention. We require also that $`u_r`$ is either decreasing or constant with increasing angle and that the density and pressure attain their maxima close to $`\theta _{max}`$.

To find solutions satisfying all of the previous requirements a Monte Carlo shooting method has been used starting from values corresponding to either sub- or super-Alfvénic starting points. Solutions have been found using as starting points successively $`\theta _{min}`$, $`\theta _{max}`$ and the stream-line turning point where the radial velocity vanishes. It has been found that the most convenient starting point is $`\theta _{min}`$. Results of these investigations and representative cases are shown below.

### 3.2 Low Mass Protostars

Most of the present models for star formation deal only with low mass stars. It is then natural to begin with this type of object. The solution shown in Fig.1 is representative of this class in the virial-isothermal approximation. It shows the variations with angle of the six self-similar variables (related to velocity ($`u_r`$, $`u_\theta `$, $`u_\varphi `$), density ($`\mu `$ and $`\mu _a`$), magnetic field ($`y_p`$, $`y_\varphi `$)) together with poloidal and toroidal Mach numbers ($`M_{ap}`$, $`M_{a\varphi }`$) and poloidal component of the Poynting flux ($`S_p`$). In this solution we have taken $`\alpha =0.3`$ and $`\mathrm{\Theta }=0.41`$ in order to compare with solutions given in FH1. As discussed in the previous paragraph, various requirements have been met; $`u_\varphi `$ and $`u_\theta `$ both vanish at the polar and equatorial boundaries, as do the magnetic field components ($`y`$ is diverging). The maxima for rotation velocity and Poynting flux, corresponding to the minimum in $`u_\theta `$, occur where the radial velocity changes signs. This point defines the total opening angle of the outflow which in this model is of the order of $`60^{}`$. The radial component of the Poynting flux $`S_r`$ changes signs as $`u_r`$ does, but its poloidal component defined by $`S_p=\sqrt{S_r^2+S_\theta ^2}`$ remains positive by definition and shows only a slight dip near maximum.

The high velocity outflow is of course much more narrowly collimated than the general outflow ($`u_r>0`$). The density $`\mu `$ is always superior to the Alfvénic critical density and both densities are almost equal where $`u_r`$ changes sign. Thus the flow is always super-Alfvénic as confirmed by the Mach number components in Fig.1. Since differences between $`y_p`$ and $`y_\varphi `$ are not evident to the eye in this figure, we have plotted the poloidal component of Poynting flux as given by Eq. 29. It shows the peak energy carried by the magnetic field to be at the turning point of the flow.

Thus this calculation produces wide angle relatively slow outflows surrounding a fast component that is identifiable with the EHV outflows or ‘molecular jets’. We shall see below that inclusion of the Poynting flux allows more collimated and faster ‘jets’ surrounded by a slower wide angle outflow than is the case without it. Therefore the magnetic field plays a crucial role in the YSO outflow despite its conservative nature.

For embedded sources (class 0/I) most of the circumstellar matter is distributed in an envelope with a typical size of about $`10^4`$ AU (Adams et al. adams (1993), Terebey et al. tere (1998), André & Montmerle andreM (1994)). The size that we will use to present our solutions will be of this order. For example, two poloidal projections of the streamlines are represented in Fig.2 out to a radius of 4000 AU. On the left panel (case A), the previous solution is shown from (nearly) $`\theta =0`$ to (nearly) $`\frac{\pi }{2}`$. In the panel on the right (case B) we show the result of varying the magnetic field and density from the previous values. The solution shown stops at $`\theta =1.33`$ and the empty region is outside of the domain of the solution. If the pressure and tangential velocity matching were done properly at the boundaries, a rotating Bondi accretion flow (Henriksen and Heaton HenriksenHeaton (1975)) might describe naturally a true accretion disc.

The self-similar structure of the streamlines does not in fact strictly match an arbitrary boundary condition at infinity, but the calculated symmetry seems to be a natural small scale limit of a more general self-gravitating circulation flow. Moreover the maximal meridional velocity close to the axial region (see Fig.1) shows the tendency for the flow to collimate cylindrically along the axis of rotation. This is in agreement with Lery et al. (lery1 (1998), 1999b ) where it is shown that cylindrical collimation is a general asymptotic behavior of magnetized outflows surrounded by an external pressure, due in the present case to the outer part of the molecular cloud or the interstellar medium.

Fig.3 presents a zoom of streamlines in the axial region view from above . Streamlines start in the equatorial region and end in the axial region with small angles ($`\theta =0.005`$ to 0.15). In addition to the meridional behavior shown previously in Fig.2, each streamline makes a spiraling approach to the axis and then emerges in the form of an helix wrapped about the axis of symmetry. The closer in angle to the equator the streamline starts, the larger the angle of rotation it makes around the axis, the closer it gets to the central mass, and the nearer to the axis it asymptotically emerges. On the other hand, for streamlines beginning well above the equator, rotation is very nearly negligible and the path almost remains in a poloidal plane.

Indirect measurements of the magnetic field in protostars are presently available (inferred from the radio flux from gyro-synchrotron emission by Ray et al. Ray (1997), and Hughes hughes (1999)). Such measurements can provide a critical test of models. In Fig.4 the angular dependences of the components of the magnetic field are shown in Gauss at a distance of $`5\times 10^4`$ AU for the same illustrative example of a low mass protostar. Magnetic fields of the order of a milligauss or less are obtained with a peak around 1 mGauss that coincides in angle with the radial turning point of the flow. As required by the quadrupolar geometry $`B_\varphi `$ and $`B_\theta `$ vanish both on the equator and on the axis, while $`B_r`$ vanishes only on the axis and at the opening angle, where the largest field intensity exists.

One remarkable prediction of these models is that the magnetic field at a given radius varies dramatically with angle (and hence probably with inclination to the line of sight). The values range from 10 microGauss to a milliGauss at $`10^4`$ AU, from $`10^2`$ to a few Gauss at 20 to 40 AU in agreement with some observations (Ray et al. Ray (1997)). Peak field strengths may reach values as high as 1 to 100 Gauss at 1 AU (e.g. Hughes hughes (1999)), as shown in Fig.5. There is also a weaker dependence on the self-similar index as is also shown in Fig.5. These predicted magnetic field strengths are surprisingly large. This is nevertheless consistent with existing (rare) observations and reinforce the idea that magnetic fields play a major role around young stellar objects (Donati et al. donal97 (1997), Guenther & Emerson gue96 (1996)).

### 3.3 Massive Protostars

To adapt the model for high mass protostars, we must find solutions with higher self-similar density (according to Eq. 15) using the method outlined in Sect. 3.1 and adjust the other variables, mainly the magnetic field in fact, until the boundary conditions are satisfied. For radiative high mass protostars, we may speculate that we may not need a magnetic field as strong as in the low mass case to either launch or collimate the outflows since strong pressure gradients are produced by radiative heating (see also Shepherd & Churchwell shch (1996), Churchwell church (1997)). However the virial-isothermal model that we give below has a stronger magnetic field.

#### 3.3.1 The Virial-Isothermal Case

Fig.6 shows the self-similar variables for a virial-isothermal massive object. This solution is also completely super-Alfvénic, just as is the previous low mass example that is overplotted with dashed lines. Although the overall opening angle of the flow (as measured by the angle at which the radial velocity changes sign) is the same in both cases, we observe that massive protostars produce faster (scaled) radial velocities throughout the positive range. There is also a more negative $`u_\theta `$ throughout this range which leads to an enhanced ‘collimation’, corresponding to an increased axial density of material. In addition, the peak rotation at the ‘turning point’ and the infall velocity just before this peak are both larger for the massive star. The smaller Alfvénic Mach numbers also show that the magnetic field is more important in this flow. This seems consistent with the ability to confine the higher physical temperature.

Our high mass virial-isothermal protostar yields larger scaled velocities than does the low mass object. This is mainly due to the larger luminosity in the high mass case. Sufficiently close to the axis both the low ($`1M_{}`$) and high ($`30M_{}`$) mass objects can produce $`300kms^1`$. However as noted previously it is not clear that streamlines at such small angles are part of the present flow.

#### 3.3.2 The Radiative Case

Massive molecular outflows present large bolometric luminosity, and therefore radiation probably plays an important role in the dynamics. Consequently radiative heating has to be investigated, taking into account consistently the diffusion equation. The radiation field and the temperature are no longer constant with angle and the radiation flux is not purely radial. The index $`\alpha _f`$ in Eq. 7 is chosen to be negative ($`0.1`$ as in FH1) and so simulates a slight radiation loss. The numerical method is similar to that used above. The self-similar radiation flux $`f(\theta )`$ is fixed at $`\theta _{min}`$ to be zero , and the non-zero value at $`\theta _{max}`$ measures the radiation ‘loss’ from the self-similar region.

Fig.7 shows a solution for the radiative case where $`\alpha =0.2`$, $`a=0`$ and $`b=2`$. The temperature $`\mathrm{\Theta }`$ is at its maximum (0.5) in the jet region and decreases towards the equator. The components of the radiation flux are much smaller than in FH1. Thus the required heating is less important than in FH1.

Fig.7 presents two illustrative radiative solutions for low ($`M0.3M_{}`$) and high mass ($`M7M_{}`$) cases. Streamlines in the poloidal plane are overplotted on contours of the hydrogen number density in the upper part of the figure, while the intensity of the radiation field is represented at the bottom. Fig.7 clearly shows the saddle type singularity of the central point. We see that a Keplerian disc naturally appears as part of the global solution as does the collimated outflow.

\[\[Caption of the figure that is given in gif format because of its size\] Streamlines and density contours of the hydrogen number density (in logarithmic scale) for low and high mass radiative solutions with $`\alpha =0.2`$. The density levels shown set between $`10^3`$ to $`10^5`$ cm<sup>-3</sup> for the low mass case, and between $`10^5`$ to $`10^7`$ cm<sup>-3</sup>in the other case. Length scales are in AU. The lower panel presents the total intensity of the radiation field (on a logarithmic scale).\]

The radiation field is smaller in the massive case in Fig.7. This is remarkable since one should expect significant radiative heating surrounding massive protostars. However, we have chosen to have the same components of the self-similar velocity for both cases, as well as the same self-similar temperature and radiation flux as input parameters. Therefore Fig.7 shows that the necessary radiation field has to be larger in low mass outflows in order to get the same scaled velocity as in the massive case. Particularly in the latter case, the global shape of the radiation field is rather similar to the solution obtained by Madej et al. (MLH (1987)), where an anisotropic radiation field was produced by scattering in thick accretion discs. It was shown in their article that collimation was clearly displayed by the radiation field in the polar region. We also find, in the massive case, that the radiation field increases away from the equator. On the other hand the low mass case almost shows an almost spherical symmetry, except on the axis.

Radiative heating increases the pressure gradient pushing outwards and therefore peak and asymptotic velocities along streamlines are larger as represented in Fig.8. This figure shows variation of total velocity with integrated distance along a streamline for both radiative and virial-isothermal cases for the same input parameters. We begin the integration in the equatorial region at about $`10^4`$ AU from the protostar and integrate out to $`10^6`$ AU in the axial region, where we attain a final angle of 0.1 radian (see FH1 for comparison). Velocity is maximum at the turning point ($`u_r=0`$), which is the closest point to the source, and is almost completely toroidal there. The inclusion of radiative heating provides a net acceleration of material which clearly increases the asymptotic axial velocity. Moreover the asymmetry in this graph indicates that the angular variation dominates the radial decrease as dictated by the self-similar forms. In fact this asymmetry is the clearest indication of the production of high velocity outflow from the relatively slow (but massive) infall.

Global differences in the solutions between virial isothermal and radiative cases remain as described in FH1. The temperature is maximum in the jet region and decreases towards the equator. The radial component of the radiation flux decreases monotically in the same direction and its $`\theta `$-component increases, being always positive.

In order to show the influence of the mass on solutions, we present variations of the total velocity with the integrated distance along streamlines as defined in the previous plot for different masses in Fig.9. For all the curves, the input parameters are kept the same ($`\alpha =0.2`$, $`\mathrm{\Theta }=0.5`$) except for the mass of the central object and therefore the self-similar density $`\mu `$. Even the initial components of velocity and magnetic field remain the same. Only the starting distance of integration from the central object varies in order to differentiate the various curves. The larger luminosities for massive objects explain the larger velocities. In the solutions shown in Fig.9, the total velocity is multiplied by approximately a factor 5 as material flows from the equatorial region to the axial region. The efficiency of heating as a driving mechanism grows with the input values of the temperature $`\mathrm{\Theta }`$ and of the components of the self-similar radiation flux $`f`$. Indeed, a factor of 10 or more can be obtained just by increasing these parameters. Moreover, the figure shows that the highest velocity components ($`100kms^1`$) are most tightly collimated along the axis ($`\theta =0.02`$). Thus the most massive protostars produce the fastest flows with maximum radial velocities in the axial region.

We conclude that radiative heating combined with Poynting flux driving provides an efficient mechanism for producing high velocity outflows when the opacity is dominated by dust.

### 3.4 Parameter Study

\[\[Caption of the figure that is given in gif format because of its size\] Monte Carlo shooting: Plots of maxima of $`u_\varphi `$ at the turning point ($`u_r=0`$) with respect to the maxima of $`u_\theta `$ for different values of the temperature parameter ($`\mathrm{\Theta }`$=0.01;0.16;0.26;0.36). ($`\alpha =0.2`$)\]

We use a Monte Carlo exploration of our parameter space in order to study general trends in the model. The self-similar index and the temperature are fixed and the six input parameters (velocity and magnetic field scaled components and the density) are randomly chosen. The conditions presented in Sect. 3.1 are required for any solution to be considered ‘good’. One should note that there exists a broad range of solutions that possess favorable characteristics. Only the most significant results are reported here.

We plot the maxima of $`u_\varphi `$ as functions of maxima of $`u_\theta `$ for various values of the self-similar temperature ($`\mathrm{\Theta }`$=0.01;0.16;0.26;0.36), and for a given index ($`\alpha =0.2`$). These maxima occur at the turning point where the radial velocity vanishes. Each point shown here represents a ‘good’ solution. Rotation is directly measured by $`u_\varphi `$, and the tendency for material to go towards the axis, i.e. the collimation, is related to $`u_\theta `$. We find, as a global trend, that when $`(u_\theta )_{max}`$ decreases, $`(u_\varphi )_{max}`$ remains almost constant ($`(u_\varphi )_{max}0.25(u_\theta )_{max}+constant`$) until $`(u_\theta )_{max}`$ reaches a threshold value. Then $`(u_\varphi )_{max}`$ decreases faster than $`(u_\theta )_{max}`$, which is what was anticipated analytically for super-Alfvénic flow when $`\alpha =1/4`$ (see appendix B). In addition, there is a limiting $`(u_\varphi )_{max}`$ for each $`(u_\theta )_{max}`$ at a given temperature, below which solutions are not found. When $`\mathrm{\Theta }`$ increases, which physically corresponds to larger pressure gradients, $`(u_\varphi )_{max}`$ decreases while $`(u_\theta )_{max}`$ becomes less negative. If the magnitude of the latter quantity is taken as a measure of the collimation, we see that rotation and collimation decrease together with increasing temperature and therefore probably with the bolometric luminosity of the central object.

It is found that solutions are less influenced by density (and hence mass) than by temperature. Nevertheless the distribution of opening angle is not found to change its form with temperature. In particular the most probable value is always found essentially at the angle shown in the two sample solutions. The variation is such that the larger opening angles are found at the lower left of the band for each temperature, that is with maximal poloidal collimation flow and the smallest possible maximum rotation.

It is also found that larger opening angles are associated with smaller magnetic fields. The corresponding results are not reported here, but in fact the variation of the opening angle is rather large for a relatively small magnetic field variation (a reduction of 25$`\%`$ in the magnetic field magnitude can widen the opening angle from $`35^{}`$ to $`60^{}`$). If the magnetic field varies secularly as the evolution proceeds, then a sequence of our models could be regarded as a series of ‘snapshots’ of the protostellar evolution. This evolutionary sequence would show that the opening angle increases with time, as the magnetic field becomes less important. This would be consistent with the notion that the youngest outflows are generally the most collimated. Indeed, this has been observationally verified by Velusamy and Langer (VeluLan (1998)), who provide evidence that outflows widen in time. Moreover when the magnetic field reaches a lower threshold in our model, the solution becomes a pure outflow. Therefore, as obtained with our model, if the opening angle broadens sufficiently, it may ultimately cut off the accretion supplying new material for the outflow.

The angular dependence of the self-similar variables are represented in Fig.10. The only parameter that varies is the self-similar temperature. The figure (left panel) illustrates the fact that collimation decreases with temperature as already mentioned. Moreover variations of the self-similar temperature give rise to different density profiles in the infalling ‘disc’ region, as seen on the right hand panel of Fig.10. The density decreases close to the equator for the largest values of the temperature with a maximum well away from the equator. Such solutions accrete so rapidly in the equatorial plane that they have a reduced density relative to their rotationally supported surroundings. These solutions resemble the special case $`\alpha =1/4`$ as discussed in appendix B. On the other hand the temperature has little effect on the density profile in the axial region, although it broadens the central ‘funnel’ substantially.

We note that solutions for cold flows ($`\mathrm{\Theta }=0`$) have also been obtained. In such a case, the circulation is powered only by Poynting flux since our Bernoulli integral (Eq. 32) applies. We do not expect asymmetric velocities at the two infinities on the stream lines. For the same set of parameters, the total velocity in the ‘interacting’ region is smaller for cold flows and the fast rotating region (where $`u_\varphi `$ is maximum) is narrower relative to finite pressure flows. Moreover the Alfvénic Mach numbers almost reach unity at the peak rotation point, showing that the flow is only slightly super-Alfvénic there, and that in fact the two terms in Eq. 32 are probably comparable.

Dependences of the self-similar index have also been investigated. The main results are illustrated by Fig.11 where the radial and longitudinal components of velocity are plotted for three values of the self-similar index but with the same self-similar temperature. For more negative self-similar index one gets dramatically smaller $`\theta `$\- and $`\varphi `$\- velocities (recall that $`\alpha =1/2`$ yields pure radial accretion) and a smaller change in the radial velocity.

## 4 Observational Consequences of the Model: Synthetic $`{}_{}{}^{13}CO`$ Spectral Lines and Maps.

In this Section, we compute $`{}_{}{}^{13}CO`$ (J=$`10`$) emission lines for several of the outflow solutions discussed in Sect. 3.1. Specifically, we provide examples of spectra, channel maps, position-velocity diagrams, and intensity-velocity diagrams for one example of each of the following: the low mass solution of Sect. 3.2, the radiative high mass solution of Sect. 3.3.2, and the virial-isothermal high mass solution of Sect. 3.3.1. The importance of this analysis is that it allows us to directly compare the observable features of our models with real observational data. However, we cannot fully explore the observational consequences of our models in the present paper. Our purpose in this Section will be mainly to outline the qualitative features of our models, leaving a complete exploration of the parameter space for the second paper in our series.

The method that we follow is essentially the same as that outlined in FH2, except that the code used there has been substantially improved in a number of respects. Firstly, we are now able to generate results with much higher spectral and spatial resolution. Secondly, although we compute spectra on a grid of pencil beams through the outflow source, we convolve our spectra with a Gaussian telescope beam to more accurately simulate observational results. FH2, on the other hand, did not perform this convolution, which resulted in many spuriously sharp features in their maps. Thirdly, we do not embed our solutions in a background of molecular gas to simulate the molecular cloud, which presumably surrounds the outflow. FH2 demonstrated that the presence of background material has no significant effect on the spectra or maps, except for in the lowest, and hence least interesting, velocity channels. Since we are primarily interested in emission at relatively high velocities, corresponding to the outflow, emission or absorption by relatively slow moving background gas is of no consequence. As in FH2, the primary limitation of our analysis is that we assume local thermodynamic equilibrium. Although this is probably not strictly true, we do expect the $`{}_{}{}^{13}CO`$ level populations to be approximately in equilibrium with the local temperature provided that the optical depth is at least moderately high. Therefore, we hope to capture at least the character of the emission, if not the precise intensity.

### 4.1 Low Mass Solutions

We have assumed a mass of $`1M_{}`$ for the low mass solution, which implies a fiducial radius of $`r_o=1400AU`$ using Eq. 13. We have also chosen a inclination angle of $`45^{}`$ between the symmetry axis of the outflow and the plane of the sky. With the mass and $`r_o`$ determined, the density, temperature, and velocity are determined at any position by the dimensionless parameters $`\mu (\theta )`$ and $`\mathrm{\Theta }(\theta )`$, and the radial scalings given by Eqs. 1-7. Therefore, by constructing a line of sight through the cloud, we can compute the emission (expressed as a radiation temperature) as a function of velocity parallel to the line of sight; examples of these spectra are shown in Fig.12. We have computed spectra on a $`61\times 120`$ grid of map positions across the outflow.

We find that the spectra shown in Fig.12 vary substantially with map position across the outflow. Spectra computed near the map projection of the outflow axis typically show relatively high velocity wings that may extend up to several $`kms^1`$, but the emission generally becomes too weak to show past $`2kms^1`$. We note that the strongest emission always occurs at low velocities, with relatively weak emission in the wings, as is the case for real outflows. Spectra may contain either a single wing or both red and blue-shifted wings depending on whether a particular line of sight crosses only a single lobe of the outflow or both. There are also many cases where line wings are absent, but well-defined “shoulders” are present on the low velocity dominated spectra.

In the equatorial regions of the maps, double-peaked spectra are found with relatively low velocity peaks. These spectra are mainly due to the radial infall that dominates in the equatorial regions. It is unlikely that the double-peaked profiles are due to rotation, since radial infall velocities exceed the rotational velocities at all angles $`\theta `$ except for very near the radial turning point of the outflow. Furthermore, many of the spectra have a slight asymmetry in which the blue-shifted peak shows stronger emission. This is suggestive of the well-known outflow asymmetry, which is due to the more efficient self-absorption of blue-shifted emission by overlying material moving towards the observer (Shepherd & Churchwell shch (1996), Gregersen et al. greg (1997), Zhou et al. zhou (1996)). The effect is slight in this case since the low mass solution, with a $`1M_{}`$ central object, produces a peak optical depth of only $`0.2`$ along most lines of sight. The massive solutions, discussed in Sect. 3.3.2 and 3.3.1, produces much higher optical depths and leads to a much more pronounced asymmetry.

Fig.13 shows a set of channel maps for the low mass solution. To make these synthetic maps, we have integrated the emission over $`1kms^1`$ wide channels from $`6.5`$ to $`6.5kms^1`$. We have deliberately excluded the emission with $`|v|<1kms^1`$, since we find that the lowest velocity emission is essentially featureless. We have convolved the maps with a Gaussian telescope beam with a FWHM of $`800AU`$; for example, this corresponds to a $`5arcsec`$ beam at a distance of $`160pc`$.

The outflow is apparent at all velocities shown in our maps. We obtain a relatively wide outflow “cone” at the lowest velocities ($`0.5`$ to $`1.5kms^1`$) and very well collimated jet-like emission apparent at velocities between $`1.5`$ and $`3.5kms^1`$. At higher velocities, the jet-like feature breaks up into outflow “spots” that move away from the central object as the velocity increases. The most important feature of these channel maps is that the opening angle of the outflow gradually decreases as the magnitude of the velocity increases. As discussed in FH2, this effect is entirely due to the angular velocity sorting of the outflow solutions. Our models always have the property that the outflow velocity increases towards the axis of symmetry. Thus, our models naturally produce outflows in which most of the material is poorly collimated and moves at relatively low velocities, but the fastest jet-like components are very well collimated towards the axis of symmetry. This property has, in fact, been observed in a number of outflow sources (Bachiller et al. bachiller (1990), Guilloteau et al. Guilloteau (1992), Gueth et al Gueth (1996)).

Several authors have noted that molecular outflows seem to be characterized by a power-law dependence of total integrated intensity as a function of velocity (Rodriguez et al. Rodriguez (1982), Masson & Chernin Masson (1992), Chandler et al. Chandler (1996)). In Fig.14, we have summed all of the spectra in order to show how the intensity is distributed with velocity. The solid line in the figure represents the blue-shifted emission, while the dashed line represents the red-shifted emission. We note that there are no important differences due to the somewhat low optical depth of our solution. We find that the intensity $`I`$ is well fit by a power law $`Iv^{3.5}`$ over all velocities greater than approximately $`0.2kms^1`$. This power-law behavior has, in fact, been found for several real outflows , and the index is in reasonable agreement with the available data (Cabrit et al. cab (1996), and Richer et al. PPIV (1998)). It remains to be seen, however, how sensitive the power-law index is to the parameters of our model. This important issue will be resolved in the second paper in this series, where we more completely explore the line emission of our models.

An important feature of Fig.14 is that the emission falls very rapidly with increasing velocity. For example, the intensity at even a rather modest velocity of $`3kms^1`$ is only $`0.01\%`$ of the peak intensity. At some velocity, the emission will fall below the noise threshold that is always present in real observations. Past this velocity, the emission shown in Fig.13 would be undetectable.

In Fig.15, we show a velocity-position diagram for a cut along the outflow axis. We have actually computed spectra along several cuts parallel to the outflow axis and convolved them with a Gaussian beam to produce a single velocity-position diagram. We find that there is substantial emission at low velocities along the entire cut. The highest velocities are found near the map centre, as must be the case since $`vr^{1/2}`$ by Eq. 1. We note that the emission contours in this figure are logarithmically spaced; thus, we find that the most intense emission occurs at low velocities at all map positions.

The outflow is apparent in Fig.15 by the emission extending out to approximately $`\pm 8kms^1`$ on either side of the outflow. Moreover, the general appearance of the position-velocity diagram is suggestive of a globally accelerated flow since higher velocities generally appear at greater distances from the central source. How can this be reconciled with the fact that velocity decreases with radius in our model? It was shown in FH2 that the apparent acceleration is due to the unique velocity sorting that is predicted by our model. The radial velocity always increases with decreasing polar angle near the outflow axis. Our cut crosses progressively smaller angles with increasing distance from the central object, and therefore encounters higher velocity material.

The ragged appearance of the highest velocity contours is entirely due to the limitations of our numerical procedures at high velocities. In our model, the highest velocity components are extremely localized near the central star and the outflow axis. Thus, whether or not a pencil beam passes through such a component is a matter of chance when the spatial extent falls below the grid spacing. This apparently happens when $`|v|5kms^1`$.

### 4.2 Massive Solutions

#### 4.2.1 The Radiative Case

It is useful to compare the $`{}_{}{}^{13}CO`$ emission of the radiative solution discussed in Sect. 3.3.2, with the low mass solution discussed above. We have assumed a mass of $`10M_{}`$, which implies a fiducial radius of $`r_o=7900AU`$ according to Eq. 13. We find, heuristically, that the peak optical depth in the line centre decreases with increasing central mass. Since the mass is a free parameter, we have tuned it to achieve realistic optical depths of order unity. The mass is larger than for the low mass solution discussed above, but the effects of radiative heating are more likely to become important in outflows surrounding more massive protostars, in any case. The only other parameter that we have changed is the FWHM of the Gaussian beam that we convolve with our solution; here, we can use a slightly larger and more realistic 10 arcsec beam without smearing out the details of the solution. As in Sect. 4.1, we have computed spectra on a $`61\times 120`$ grid of map positions assumed a $`45^{}`$ inclination angle.

Qualitatively, we find a great deal of similarity with the low mass solution. The line profiles shown in Fig.16 have similar structure, and the outflow shown in Fig.17 is collimated to a similar degree. The main difference is the velocity of the flow. In the radiative case, we find that there is significant emission out to at least $`10kms^1`$, and we find some low intensity emission out to $`15kms^1`$; this is most clearly illustrated by the position velocity diagram shown in Fig.19. Examining the intensity-velocity diagram shown in Fig.18, we find that the intensity $`Iv^{3.6}`$ for velocities greater than about $`0.5kms^1`$. It is striking that nearly the same power law index should be obtained for both the low mass and radiative solutions. Admittedly, we do not fully understand the reasons for this similarity in behavior at present. We also doubt that this sort of behavior is a generic feature of our models. Nevertheless, we are encouraged that at least some of our models are capable of reproducing such realistic observed properties.

#### 4.2.2 The Virial-Isothermal Case

We have also computed the $`{}_{}{}^{13}CO`$ emission of the high mass solution discussed previously in Sect. 3.3.1. We find that the optical depth is extremely high unless the mass of the central protostar is chosen to be rather large; therefore, we have used a mass of $`100M_{}`$ and a corresponding fiducial radius of $`r_o=44000AU`$ to compute the spectra. Even using such a large central mass, we find that the peak optical depth in the line centre is very high, with values exceeding 5 in many positions. Such high optical depths cause real difficulty for the present version of our code, since any given line of sight encounters the “photosphere”, where the optical depth at some velocity is approximately unity, rather abruptly. Thus, a great deal of care was taken to ensure that the solution was well-sampled along each line of sight. We note that the corresponding increase in computational time required us to decrease the spatial resolution of our maps. Here, we employ a $`41\times 80`$ grid, as opposed to the $`61\times 120`$ grid used elsewhere. The inclination angle is $`45^{}`$, as in the previous two cases.

The spectra shown in Fig.20 are much more complex than the spectra of the low mass or radiative solutions. Many show multiple sharp peaks and very extended wing emission out to well past $`20kms^1`$ (as in Shepherd & Churchwell shch (1996), Gregersen et al. greg (1997) for CO maps of molecular outflows, and Doeleman et al. doeleman (1999) for SiO masers). We are quite certain that the multiple peaks are real, and not numerical artifacts, since we have verified that their location and appearance is independent of the spectral resolution and the degree to which a solution is sampled along a line of sight. It remains unclear to us in detail why such jagged peaks should be present, although it must depend on the variation of optical depth, density and velocity field within the ‘beam’. Thus behavior should be studied as a function of inclination to the line of sight.

Despite the unusual appearance of the spectra, the channel maps shown in Fig.21 are quite reasonable. A well-collimated outflow is apparent at velocities greater than about $`10kms^1`$. The position velocity diagram shown in Fig.23 shows significant emission at most positions out to approximately $`20kms^1`$. We obtain emission at higher velocities (up to approximately $`40kms^1`$, but the emission is restricted to the central map positions nearest the protostar. The intensity-velocity diagram (Fig.22) again shows a power law behavior, but the slope is different than that obtained in the previous two cases; here we find that $`Iv^{4.1}`$, which is still inside the range allowed by the available observations (Cabrit et al. cab (1996), Richer et al. PPIV (1998)).

## 5 Discussion

### 5.1 The Alfvénic Point

All the solutions presented previously were completely super-Alfvénic. Here we report a trans-Alfvénic solution in order to show why they are generally not preferred. In such a solution, the outflow size is reduced and the ultimate velocities are smaller. Therefore this type of solution is less efficient in producing high velocity outflows. Entirely sub-Alfvénic solutions can also be found but they admit only pure accretion. Fig.24 shows the density profile for two super-Alfvénic solutions (the two upper most plots A and B) and one trans-Alfvénic solution (C). Case (A) presents a model where the density is well above the critical Alfvénic density. In this case, the Alfvénic Mach number is larger than unity. The second case (B) has almost the critical density at the turning point, and consequently the Alfvénic Mach number almost reaches unity at this point. However the solution always remains super-Alfvénic. Finally the lower panel (case C) displays a trans-Alfvénic solution showing that the size of the outflow has been reduced drastically. Therefore the most interesting solutions for circulation flow appear to be super-Alfvénic solutions.

### 5.2 Comparison with FH1

Fiege & Henriksen (fiege1 (1996)a, FH1) have shown that $`r`$-self-similar quadrupolar circulation could be a relatively realistic model for bipolar outflow. They did not include Poynting flux that could assist the outflow as in the present paper. The inclusion of this effect allows us to obtain faster and more collimated outflows. The path followed by material also differs from FH1 as shown in Fig.25, where streamlines from FH1 and the present models are plotted. All the solutions start from approximately the same initial position in the equatorial region with the same set of parameters. Streamlines from the present model get closer to the source than the FH1 solution, for both radiative and virial-isothermal cases. Thus, in the present model, the flow passes relatively close to the central mass. The limiting speed is about the escape speed from the central mass (and the jet speed) of a few hundred $`kms^1`$. We note that the most collimated of the the three solutions is the radiative case.

### 5.3 Stability Analysis

The steady configurations presented in this article should develop strong shocks in the axial region due to the strong collimation created by magnetic and pressure forces. Moreover current-carrying jets, as in the present model, are liable to Kelvin-Helmholtz (KH), Pressure Driven (PD) and magnetic instabilities driven by the electrical current. These so called current driven (CD) instabilities have been studied only recently in the context of astrophysical jets (Appl appl (1996), Appl et al. LBA (1999)). We use the results of the latter linear stability analysis in order to get a rough estimate of CDI and KHI in our magnetic configuration. Growth rates $`\mathrm{\Gamma }`$ of the different modes show that the CD helical mode ($`m=1`$) dominates the CD pinching mode ($`m=0`$) and that CD and KH instabilities become comparable when the Mach number approaches unity. This is the case in the present self-similar model. The CD instability (with a typical wavelength $`\lambda _{CD}\frac{2}{5}R_{outflow}`$, $`R_{outflow}`$ being the characteristic size of the outflow) should be located on a current sheet where reconnection and particle acceleration might occur. KH pinching modes ($`\lambda _{KH}3\times R_{outflow}`$) should dominate KH helical modes, and consequently give rise to internal shocks in the outflow. All these results require numerical simulations that we plan to study in another paper.

## 6 Conclusions

In this paper we have presented a model based on the $`r`$-self-similarity assumption applied to the basic equations of ideal axisymmetric and stationary MHD, including Poynting flux. Detailed comparisons with the observations have been computed and characteristic scales of the problem have been given as functions of source properties. The luminosity needed for radiative heating is smaller even while driving faster outflows than in the previous model (FH1). The model geometry implies a natural connection between the fast ionized jets seen near the polar axis of the wind, and the slower and less-collimated molecular outflows that surround them, although the jets may be partly due to central activity not included in our model.

The most massive protostars produce the fastest flows with maximum radial velocities in the axial region. Radiative heating produces faster outflows compared to the virial-isothermal case, for both low and high mass objects. Larger opening angles are associated with smaller magnetic fields. Consequently, a gradual evolutionary loss of magnetic flux may result in outflows that widen as they age.

Synthetic spectral lines from $`{}_{}{}^{13}CO`$ (J=$`10`$) allow direct comparison with observational results via channel maps, maps of total emission, position-velocity and intensity-velocity diagrams. The model reproduces well observational features. Due to internal instabilities in the most collimated ‘jet’ part of the flow, the time evolution of the steady model should give rise to regularly spaced knots, and possibly excitation, ionization and non-thermal particle acceleration due to field annihilation and shock dissipation. Thus the model in its current state of development shows that radiative heating combined with Poynting flux driving is efficient in producing high velocity outflows when the characteristic luminosity of the forming star (as deduced from observations) is used. More extensive searches in parameter space, the detailed fit to specific sources, and time dependent modeling all remain to be done. Although we do not discuss the regions excluded from the self-similar region of the flow, the current correspondence to observational features is sufficiently exact that we expect simultaneous infall/outflow (‘circulation’) to be an essential part of any realistic model.

###### Acknowledgements.

We would like to thank Sylvie Cabrit for constructive and useful discussions.    

## Appendix A The Equations

The mass flux conservation equation remains the same as in FH1:

$$\left(1+2\alpha \right)\mu u_r+\frac{1}{\mathrm{sin}\theta }\frac{d}{d\theta }\left(\mu u_\theta \mathrm{sin}\theta \right)=0,$$
(33)   
while the self-similar variable $`y`$ is simply replaced by its poloidal component in the following set of equations:

Magnetic Flux conservation:

$$\frac{\left(\alpha +5/4\right)u_r}{y_p}+\frac{1}{\mathrm{sin}\theta }\frac{d}{d\theta }\left(\frac{u_\theta \mathrm{sin}\theta }{y_p}\right)=0$$
(34)   
Radial component of momentum equation:

$`u_\theta {\displaystyle \frac{du_r}{d\theta }}\left(1M_{ap}^{}{}_{}{}^{2}\right)\left(u_{\theta }^{}{}_{}{}^{2}+u_{\varphi }^{}{}_{}{}^{2}\right)\left(1{\displaystyle \frac{\alpha +1/4}{M_{ap}^2}}\right)`$   
$`{\displaystyle \frac{u_r^2}{2}}\left({\displaystyle \frac{3}{2}}2\alpha \right)\mathrm{\Theta }+1+{\displaystyle \frac{u_ru_\theta }{M_{ap}^2y_p}}{\displaystyle \frac{dy_p}{d\theta }}=0`$ (35)   
Other equations of the system contain terms with mixed toroidal and poloidal components; Angular Momentum conservation:

$`{\displaystyle \frac{1}{u_\varphi }}{\displaystyle \frac{du_\varphi }{d\theta }}\left(1{\displaystyle \frac{1}{M_{ap}M_{a\varphi }}}\right)+{\displaystyle \frac{u_r}{u_\theta }}\left(1/2{\displaystyle \frac{\alpha +1/4}{M_{ap}M_{a\varphi }}}\right)`$   
$`+\mathrm{cot}\theta \left(1{\displaystyle \frac{1}{M_{ap}M_{a\varphi }}}\right)+{\displaystyle \frac{1}{y_\varphi M_{ap}M_{a\varphi }}}{\displaystyle \frac{dy_\varphi }{d\theta }}=0`$ (36)   
Faraday’s law plus zero comoving electric field:

$$\frac{d(u_\varphi u_\theta )}{d\theta }+\left[\alpha \frac{1}{4}\right]u_\varphi u_r+u_\varphi u_\theta \frac{d}{d\theta }\mathrm{ln}\left[\frac{1}{y_p}\frac{1}{y_\varphi }\right]=0$$
(37)   
$`\theta `$-component of momentum equation:

$`{\displaystyle \frac{u_ru_\theta }{2}}\left(1{\displaystyle \frac{2\alpha +1/2}{M_{ap}^{}{}_{}{}^{2}}}\right)+u_{\varphi }^{}{}_{}{}^{2}\mathrm{cot}\theta \left(M_{a\varphi }^{}{}_{}{}^{2}1\right)`$   
$`+{\displaystyle \frac{u_r}{M_{ap}^{}{}_{}{}^{2}}}{\displaystyle \frac{du_r}{d\theta }}+u_\theta {\displaystyle \frac{du_\theta }{d\theta }}+{\displaystyle \frac{u_\varphi }{M_{ap}^2}}{\displaystyle \frac{du_\varphi }{d\theta }}{\displaystyle \frac{u_{r}^{}{}_{}{}^{2}}{y_pM_{ap}^{}{}_{}{}^{2}}}{\displaystyle \frac{dy_p}{d\theta }}`$   
$`{\displaystyle \frac{u_{\varphi }^{}{}_{}{}^{2}}{y_\varphi M_{a\varphi }^{}{}_{}{}^{2}}}{\displaystyle \frac{dy_\varphi }{d\theta }}+{\displaystyle \frac{d\mathrm{\Theta }}{d\theta }}+{\displaystyle \frac{\mathrm{\Theta }}{\mu }}{\displaystyle \frac{d\mu }{d\theta }}=0`$ (38)   
To treat the radiative heating we either hold $`\mathrm{\Theta }`$ constant (virial-isothermal case) as discussed above, or we include the equations of radiative diffusion (radiative case) as in FH1. The equations (17), (18) and (22) of FH1 continue to apply provided that the optical depth given above is $`2/3`$.

## Appendix B The First integrals

Equations (33) and (34) together yield the first integral

$$\mu (u_\theta \mathrm{sin}\theta )^{\frac{1/4\alpha }{5/4+\alpha }}y_p^{\frac{1+2\alpha }{5/4+\alpha }}q_1/4\pi ,$$
(39)   
where $`y_p`$ replaces $`y`$ in FH1. Moreover directly from Eq. 33 and the poloidal stream-line equation $`dr/rd\theta =u_r/u_\theta `$ we have the stream-line integral

$$r^{(1+2\alpha )}\mu u_\theta \mathrm{sin}\theta \eta _1.$$
(40)   
From this we see immediately that $`\alpha >1/2`$ if we are to have quadrupolar stream-lines for which $`r\mathrm{}`$ as $`u_\theta 0`$. A second stream-line integral follows from Eq. 34 and the stream-line equation, but it is not independent of Eqs. 39 and 40 since it follows by eliminating $`\mu `$ between these two. The limit $`\alpha <1/4`$ follows immediately from Eq. 39 if we require finite $`\mu `$, $`u_\theta =0`$ and odd symmetry in the magnetic field at the equator. The integral then implies implies that $`y_p\mathrm{}`$ and so the field passes through zero in the equatorial plane. This conclusion is avoided on the polar axis since the radial velocity is free to become very large there. This choice has the additional merit that the equator is super Alfvénic as well as the axis.

Using Eq. 34 together with Faraday’s law (37) we obtain the integral

$$u_\varphi \left[\frac{1}{y_p}\frac{1}{y_\varphi }\right]\left(\frac{\mathrm{sin}\theta }{y_p}\right)^{\frac{1/4\alpha }{5/4+\alpha }}u_\theta ^{\frac{3}{2(5/4+\alpha )}}=q_2.$$
(41)   
Note that the constant $`q_2`$ is very interesting because it can be used to measure the strength of the electric field which is in fact

$$E=\left[\frac{GM}{cr_o}\right]\sqrt{\frac{M}{r_o^3}}\left(\frac{r}{r_o}\right)^{\alpha \frac{5}{4}}\left[\frac{1}{y_p}\frac{1}{y_\varphi }\right]\left(u_\varphi \times u_p\right).$$
(42)   
The integral is not present in FH1 since the electric field is everywhere zero there.

In general we can not find an angular momentum integral explicitly (this is related to the difficulty of expressing scale invariance in an action principle; B. Gaffet, private communication) but in the special case of similarity index $`\alpha =1/4`$ such an explicit integral exists in the form (cf FH1):

$$\frac{u_\theta \mathrm{sin}\theta }{M_{ap}^2}q_3\left(u_\varphi \mathrm{sin}\theta \right)^3\left(1\frac{1}{M_{ap}M_{a\varphi }}\right)^3.$$
(43)   
Although this is a very special case not normally related to our numerical solutions (the mass density is constant in radius so that eventually the dominance of the central mass is broken as $`r`$ increases), it is instructive to consider it further. Eq. 39 becomes simply $`\mu =(q_1/4\pi )y_p^1`$ so that $`M_{ap}^2=q_1y_p`$ and $`M_{a\varphi }^2=q_1y_\varphi ^2/y_p`$. Then Eqs. 41 and 43 can be solved together for $`y_p`$ and $`y_\varphi `$ to give

$$\frac{1}{y_p}=\frac{q_2q_3\mathrm{sin}^2\theta +q_1q_3u_\varphi u_\theta \mathrm{sin}^2\theta }{q_3u_\varphi u_\theta \mathrm{sin}^2\theta +u_\theta ^2/u_\varphi ^2},$$
(44)   
and

$$\frac{1}{y_\varphi }=\frac{q_1q_3u_\varphi ^4\mathrm{sin}^2\theta q_2}{q_3u_\varphi ^4\mathrm{sin}^2\theta +u_\varphi u_\theta }.$$
(45)   
We therefore observe that near the equator $`y_p`$ can be infinite only if $`u_\varphi 0`$ faster than $`u_\theta `$. This feature is shared with our numerical solutions, but of course the equator is strictly outside the domain of the solution. From the second equation we see that $`y_\varphi `$ tends to zero at the equator under the same conditions. If one were to insist that $`u_\varphi 0`$ at the equator, then Eq. 43 shows that we need $`M_{a\varphi }M_{ap}=1`$ there, or by the above Mach number definitions $`y_\varphi =1/q_1`$. Our solutions for $`y_p`$ and $`y_\varphi `$ show that this is only possible if $`q_2=0`$, that is with zero Poynting flux. Thus it seems from our analysis of this special case (however we always find $`u_\varphi =0`$ at the equator in our solutions) that the presence of a non-zero Poynting flux is not compatible with a Keplerian equatorial disc. The material there is falling radially towards the star. This is a possible non-linear end state for an instability that couples disc rotation to Alfvén waves propagating out of the plane (Shu et al. shu94 (1994), Dendy et al. dentagger (1998)).