text:
# A Global Jet/Circulation Model for Young Stars

## 1 Introduction

Star formation occurs in molecular clouds. Observations have shown that the accretion phase, during which the central object builds up its mass, is very often accompanied by the powerful ejection of prominent bipolar outflows, that are comprised of fast jets (Ray et al.  (1996); Eislöffel et al. 2000) ; Reipurth et al.  (2000)) and slower and less collimated molecular outflows (Bachiller (1996)). A topic of great theoretical and observational interest is the dynamical connection between the jets and the molecular outflows. Strong arguments, such as the large masses observed in the CO flows, suggest that most of the molecular outflow material consists of swept up ambient gas (Bachiller & Tafalla (1999)). The only possible source for this gas is the surrounding ambient cloud. This implies that a driving agent accelerates molecular cloud material and gives rise to the observed CO lobes. An obvious candidate for such an agent is the central fast jet. The dynamical relationship between the optical jet and the surrounding bipolar molecular outflows, however, is controversial (Snell et al.  (1980); Rodriguez et al.  (1980); Bally & Lada (1983); Welch et al.  (1983)). At the heart of the controversy are the following questions. Has the central jet or wind sufficient thrust to drive the molecular outflow? If yes, how is the momentum transferred to the CO lobes?

There seems to be a consensus that jets have a magneto-centrifugal origin, and they are either launched from the accretion disk (disk wind) (Blandford & Payne (1982); Pelletier & Pudritz (1992); Wardle & Königl (1993); Sauty & Tsinganos (1994); Paatz & Camenzind (1996); Ferreira (1997); Königl & Pudritz (2000); Contopoulos & Sauty (2001)) or from the location of the interaction of the protostar’s magnetosphere with the disk (X-wind) (Shu et al.  (1994); Ferreira et al.  (1999); Shu et al.  (2000)). In the disk wind model, the engine consists of a Keplerian disk threaded by a magnetic field that is either generated in situ, or advected in from larger scales, while X-winds are magnetised stellar winds where the interaction of a protostar’s magnetosphere with the surrounding disk results in the opening of some magneto-spheric field lines. On the other hand, the precise mechanism generating molecular outflows is poorly understood. The latter are generally believed to be driven by a jet (e.g. Padman et al. 1997, Chernin & Masson (1991)) or by a wide-angle wind (Barral & Cantó (1981); Smith (1986); Shu et al.  (1991)), but all the proposed possibilities present difficulties (Cabrit & Bertout (1992)). For example, in the case of a wide-angle wind with a steady wind-blown cavity, the shocked wind converging at the tip of the cavity should form a jet which is not observed. However, in the case of the wide-angle X-wind model, Shang, Shu & Glassgold (1998), for example, have shown that a Herbig-Haro jet appears along the rotation axis, i.e., , in the central part of the cavity. On the other hand, for such models, the ambient medium has a density that decreases like $`r^2`$ away from the core, which cancels the decrease in wind ram pressure due to the radial expansion, ensuring a constant shell velocity. This, however, does not seem to be consistent with the shape observed for the largest flows (Cabrit & Bertout (1992)).

In the case of jet-driven molecular outflow models, or “prompt entrainment” mechanism, the bulk of the molecular outflow is accelerated ambient gas near the head of the jet or more precisely along the wings of its associated bow shock (e.g. Chernin & Masson 1995). This idea is supported by some recent simulations (Micono et al. 2000), where it is shown that a large fraction of the initial jet momentum can be transferred to the ambient medium. However, in some cases, the jet may fail to provide the necessary momentum, as in L1551 IRS 5 (Fridlund & Liseau (1998)). Moreover, observations frequently suggest that jets decelerate near their outer edges, probably due to entrainment of the surrounding material, but show little or no deceleration near the flow axis (Giovanardi et al.  (2000)). Also, in some cases, the axis of the cavity differs from the axis of the jet, e.g. by about 20 for HL Tau (Close et al. 1997b ), or by about 40 in the case of NGC 2261 (Close et al. 1997b ; Cantó et al.  (1981)). Another major problem appears to be that steady jets produce molecular lobes with large aspect ratios (Chernin & Masson (1991)). This would mean that the jet should be changing its position over time to accelerate different parcels of molecular gas (“wandering jet”), or varying its flow velocity resulting in internal shocks which push material sideways. Finally, based on numerical simulations, the efficiency of prompt entrainment in accelerating molecules is found to be poor in the case of YSO jets with densities comparable to their ambient environment (Downes & Ray (1999)).

In addition to these problems, many bipolar outflows from massive protostars transport masses largely exceeding those of their associated stars. The bipolar outflow masses appear to range from about 10 to 4800 solar masses, with a mean value around 130 (Churchwell 1997). This probably indicates that much of the outflowing molecular gas could not have been accelerated close to the central YSO, since the protostar only represents a small fraction of the outflowing mass. Most of the models, discussed previously, face difficulties in explaining such massive outflows. Indeed, in this case, it seems likely that deflection of infalling matter into bipolar outflows may be a crucial mechanism, since entrainment and swept-up mass do not appear to be able to account for the very large observed outflow masses. This is the mechanism that we propose in the present article.

Finally, we point out that a correlation exists between the bolometric luminosity of the central exciting source and the rate of momentum injection into the outflow (Bally & Lada (1983); Cabrit et al.  (1997)). These quantities are roughly consistent with a single power law (with slope of about 0.7-0.8) across the full range of source luminosities extending over five orders of magnitude (Cabrit et al.  (1997)). This continuous trend suggests that similar mechanisms may be responsible for the production of molecular outflows from both low and high mass stars. As we will show, the present model can accommodate this observational constraint.

We start by presenting the global model in section 2, and by writing down our set of self-similar MHD equations in section 3. Then, in section 4, we show the properties of the solutions along the streamlines, while, in section 5, we discuss the main quantities that characterize the flow such as the fluxes of mass, momentum and energy, the plasma beta or the electric current. In section 6 we present some physical consequences of the model and we discuss several observational implications. Finally, we give our conclusions in section 7.

## 2 The Model

We propose a global model in which the central jet is driven by a magneto-centrifugal accretion-ejection engine, while the molecular outflow is powered by the infalling matter through a quadrupolar circulation pattern around the central object. The situation is schematically represented in Fig. 1. The molecular outflow may still be affected by entrainment from the fast jet, but this would be limited to the polar regions and it would not be the dominant factor for its acceleration.

The global model combines a jet model (I in Fig. 1), a circulation model (II) (Henriksen & Valls-Gabaud (1994); Fiege & Henriksen (1996); Lery et al.  (1999); Aburihan, Fiege, Henriksen & Lery (2001)) and an infalling envelope model (III). We suggest that the gas in all the three regions I, II, and III can be described by the self-similar heated, quadrupolar and axisymmetric magnetohydrodynamic model (Henriksen & Valls-Gabaud (1994); Fiege & Henriksen (1996); Lery et al.  (1999)), described here in section 3. In fact, it is possible to incorporate the jet model inside the self-similar circulation model, and to include the infalling envelope model near the equator. The torus and the jet were not part of the previously published models of this type (Fiege & Henriksen (1996); Lery et al.  (1999)), and their inclusion is a new result reported here. Therefore a description of the flows around YSOs can be obtained from the axis to the equator with a single set of equations written in the framework of a self-similar hypothesis.

Note that the model does not directly treat the thin accretion disk near the equator nor the accretion-ejection engine that provides the jet in proximity of the central source, since the physical processes governing these regions are substantially different from the large-scale infall/outflow region.

### 2.1 The Circulation Model

The circulation model (II in Fig. 1) has already been studied in Lery et al.  (1999), and we only recall its main features here. In the infall phase, the flow is slowed down by the increasing radial pressure gradients due to heating by the central source, and by an increase with proximity to the central object of the strength of the barrier due to the rotation and the magnetic field. This pressure ‘barrier’ deflects and accelerates, by means of its poloidal gradient, much of the infalling matter into an axial outflow. The Poynting flux included in the model increases both the velocity and collimation of the outflows by helping to transport mass and energy from the equatorial to the axial regions (see Lery et al.  (1999) for details). The outflow velocity is naturally more pronounced as the axis of rotation is approached since the pressure gradients are strongest there and the density is smallest. These gradients also act to collimate the flow, forming a “throat” in the axial region.

The main attraction for the circulation model is that large quantities of neutral material can be accelerated directly into the outflow without the need for jet entrainment, whose plausibility has long been debated. Moreover, the model can generally explain the observed opening angles and velocity structure seen in high mass systems (Fiege & Henriksen (1996); Lery et al.  (1999)).

Self-similar models cannot be globally consistent in time even if they are non-stationary, since they are, in general, ignorant of the initial conditions. The self-similar circulation model is intermediate in space between a region around the symmetry axis containing the jet, and the infalling envelope close to the disk. These regions are described in the following. Note that the axis itself and the equatorial plane are strictly excluded from the domain of self-similarity although they may be approached asymptotically.

### 2.2 The Jet Model

In most jet models, the matter is thought to be transferred from an accretion disk and/or a thick torus surrounding the central forming star, to the jets. In the most widely accepted scenario the outflow acceleration is mediated by magnetic fields that thread Keplerian type disks. In the present model, no detail is given about the origin of the jet itself but our description is compatible with most models mentioned in the introduction. We consider that the system consisting of the protostar and the thin accretion disk acts as a reservoir for the jet and only study the flow inside the jet (I in Fig. 1). The same set of equations used for region II can be applied to the jet. As for the circulation solutions, the velocity is larger closer to the axis but the dynamics is dominated by the variations of the various quantities in the radial direction, i.e., along the jet. The jet model applies only to a limited range of finite radii, where the self-similarity applies. Specifically, our model cannot address the launching mechanism of the jet near the star. The structure of a supersonic jet far from the source is disconnected from its origin. Thus, it is not essential to include a detailed description of the central engine in our model.

### 2.3 The Self-similar Infalling Envelope Model

During the formation of the star, it is expected that an accreting disk forms and feeds the central protostar from infalling ambient cloud core material; these disks may have radii from 20 AU up to several hundred AU, and masses containing a fraction of the central protostar mass. Hence, the last element of the model is the region where material is falling from the cloud onto the central object. It consists of a standard flattened dense accretion disk, not treated by the present model (i.e., the equatorial-most region in Fig. 1), and of an extended and less dense torus of infalling gas and dust (III in Fig. 1). The latter can be as large as a few thousand AU and is also described by the set of self-similar equations of the circulation model. The material, for both the circulation and the infalling envelope models, originate from the same region in the molecular cloud. We thus suggest that the same set of equations may apply in both cases. This statement would not be applicable to the thin disk where resistivity, ionization fraction and other parameters can substantially differ from the outer regions.

## 3 The Equations of the Model

### 3.1 The Variables

The solutions are developed within the context of $`r`$-self-similarity wherein a power of $`r`$ multiplies an unknown function of $`\theta `$; the spherical coordinates $`r`$, $`\theta `$ and $`\varphi `$ being used. The only physical scales that enter into our calculation are the gravitational constant $`G`$, the fixed central mass $`M`$, and a fiducial radius $`r_o`$ (Fiege & Henriksen (1996); Lery et al.  (1999)). Note that Fiege et al. expressed $`r_o`$ in terms of $`G`$, $`M`$, and the central luminosity $`L_{}`$, while Lery et al. offered a different choice. In the present work, we leave $`r_o`$ as a scale that can be freely specified for any observed outflow source. We also introduce two other parameters that are the limiting angles between the jet, the infall and the circulation zones. The first angle, $`\theta _{min}`$, goes from the axis of rotation to the boundary of the jet (angle between zone I and II in Fig. 1). The second angle, $`\theta _{max}`$, corresponds to the limit of the circulation region (angle between zone II and III in Fig. 1). The power laws of the self-similar system 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 the 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 self-similar index $`\alpha `$ is a free parameter of the solution, but must lie in the range $`1/2<\alpha 1/4`$, for simultaneous infall and outflow to occur. The radiation field can help us to understand the physical meaning of $`\alpha `$. Indeed, in this context, the self-similar behaviour that we impose on the radiation field takes the form $`F_{rad}=\left(GM/r_o\right)^{3/2}M/r_o^3\left(r/r_o\right)^{\alpha _f2}f(\theta )`$ (Fiege & Henriksen (1996); Lery et al.  (1999)). In this equation, the index $`\alpha _f`$ is a measure of the radiation energy balance as a function of radial distance. When $`\alpha _f`$ is negative, it corresponds to a loss in energy due to the net radiation leaving the region, while, for a positive $`\alpha _f`$, the system has an energy input due to the absorption of radiation. If one supposes that the opacity is predominantly due to dust, the index $`\alpha _f`$ is related to $`\alpha `$ by $`\alpha _f2(1/4\alpha )`$. Consequently, the case $`\alpha _f=0`$, where there is no net flux of radiation through the domain of our solution, corresponds to $`\alpha =1/4`$. For larger (respectively smaller) values of $`\alpha `$, there is a net loss (gain) of energy by radiation (Fiege & Henriksen (1996)).

Hence, if we assume that the gravitational potential is dominated by the central mass, i.e., self-gravitation is negligible, the equations of radiative MHD admit the following radial scaling relations for the variables

$$𝐯=\sqrt{\frac{GM}{r_o}}\left(\frac{r}{r_o}\right)^{1/2}𝐮(\theta ),$$
(1)
$$\rho =\frac{M}{r_o^3}\left(\frac{r}{r_o}\right)^{2\alpha 1/2}\mu (\theta ),$$
(2)
$$B_{\varphi ,p}=\sqrt{\frac{GM^2}{r_o^4}}\left(\frac{r}{r_o}\right)^{\alpha 3/4}\frac{u_{\varphi ,p}(\theta )}{y_{\varphi ,p}(\theta )},$$
(3)
$$p=\frac{GM^2}{r_o^4}\left(\frac{r}{r_o}\right)^{2\alpha 3/2}P(\theta ),$$
(4)
$$\frac{kT}{m_\mu m_H}=\frac{GM}{r_o}\left(\frac{r}{r_o}\right)^1\mathrm{\Theta }(\theta ),$$
(5)
$$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 ),$$
(6)
where $`𝐯`$, $`\rho `$, $`𝐁`$, $`p`$, $`T`$ and $`𝐅_{rad}`$ respectively correspond to the velocity, the density, the magnetic field, the pressure, the temperature and the radiative flux. In these equations the microscopic constants are represented by $`k`$ for Boltzmann’s constant, $`m_\mu `$ for the mean atomic weight, $`m_H`$ for the mass of the hydrogen atom. In the last equation, the index $`\alpha _f`$ is a measure of the loss (if negative) or gain (if positive) in energy by radiation as a function of radial distance. The self-similar variable directly related to magnetic field $`y`$ can be divided into poloidal and toroidal components. In the present model the two components are not equal. They are respectively defined by $`y_{p,\varphi }(\theta )=M_{ap,a\varphi }/\sqrt{4\pi \mu (\theta )}`$. Consequently the system also 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 )v_{p,\varphi }^2/\left(B_{p,\varphi }^2/4\pi \rho \right)`$. We will refer to $`\mathrm{\Theta }_0`$ as the value of the self-similar temperature $`\mathrm{\Theta }(\theta )`$ on the axis.

### 3.2 The Equations

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 `$. We further restrict ourselves to steady models (i.e., $`/t=0`$). Magnetic field and streamlines are required to be quadrupolar in the poloidal plane for the circulation model (see Lery et al.  (1999) for the specification of the boundary conditions). Under these assumptions the self-similar equations are:

1. Mass flux conservation

$$\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,$$
(7)
2. 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,$$
(8)
3. 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,`$ (9)
4. $`\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 }^21\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,`$ (10)
5. 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,`$ (11)
6. 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.$$
(12)

This set of equations can produce either circulating flows, pure infall, or pure outflow. Note that the boundary conditions guarantee mass conservation in the circulating region (II), while in the two other cases, mass fluxes are not conserved. Indeed, the source acts as a reservoir for the pure outflow (region I) while the accretion disk would act as sink for the pure infall (region III). Therefore, any discussion on net mass infall/outflow can only apply to region I and III.

## 4 Example of Solution

From a numerical point of view, we integrate the set of equations from a given solution close to the axis up to a given angle $`\theta _{min}`$. This solution near $`\theta =0`$ has finite values for $`V_r`$, $`B_r`$ and the density. The other components of the velocity and magnetic field are vanishing, otherwise they would show jump discontinuities at the boundaries. The last point from the jet model solution can be used as initial condition for the integration of the same set of equations of the circulation model. However, we have restricted ourselves to solutions in which the boundary conditions for the circulation region are zero for $`B_\varphi `$, $`B_\theta `$, $`V_\varphi `$, $`V_\theta `$, a good behaviour for the circulation solution. We note that under these hypotheses no material crosses the boundaries defined by the cones delimiting the three regions, and that the solutions are in pressure balance.

Once the circulation solution is obtained, we integrate then from this solution down to a given angle ($`\theta _{max}`$) in the equatorial region, excluding a small zone corresponding to the thin accretion disk. We also ensure that $`B_\varphi `$, $`B_\theta `$, $`V_\varphi `$, $`V_\theta `$ vanish at the interface between the circulation solution and the infalling envelope. The main difference between the two solutions is the change of signs of the $`\theta `$-components of the velocity and the magnetic field. This corresponds to a solution that goes towards the axis in the case of the circulation solution, and towards the equator for the infalling envelope solution.

The three components of the velocity and of the magnetic field for a jet solution are shown in Fig.2 as a function of the angle between the axis and a limiting angle ($`\theta _{min}=0.1`$ in this case), and at a distance of $`10^4`$ AU, the mass of the central protostar being set to one solar mass. Radial velocities of several hundred of km s<sup>-1</sup> are obtained near the axis. There, the radial component of the magnetic field also dominates the other component. The rapid increase in velocity towards the axis is easily explained as follows. From Bernoulli’s law, a decrease in density along a streamline, and hence in pressure, corresponds to an increase in velocity. This is what we observe towards the axis, where material encounters a low pressure region. Bernoulli’s law strictly applies only in the absence of a Poynting flux, in which energy is directly transferred between the magnetic field and the kinetic energy. However, it is clear from our solutions that the Poynting flux is not dominant in this region. Note that the flatness of $`b_r`$ in the innermost part of the jet is not a general behavior of the self-similar solutions but corresponds to the present set of parameters. It is noteworthy that the strength of the radial component of the field could help to stabilize the jet w.r.t. internal instabilities (Lery et al.  (2000)). In most of the jet, the angular momentum and the magnetic pressure are negligible with respect to the kinetic energy, whereas, close to the interface with region II, the hoop stress becomes stronger, due to the increase in $`B_\varphi `$. The change of slope in the variables at $`\theta 0.07`$ corresponds to the location of the steep gradient in $`B_r`$ and of the maximum amplitude of $`B_\varphi `$ in the jet. If this angle is taken as a measure of the broadening of the jet during the propagation at large scales, we may derive a full opening angle of the flow of about 8, which is in the range of values derived from the observations (Eislöffel et al. 2000) ).

The solutions in the circulation and the infalling envelope torus region are reported in Fig.3. Following the flow, the radial velocity is negative in the infalling region, $`\theta _{tp}>\pi /4`$, it changes signs at the turning point (the closest point to the source) located at $`\theta _{tp}\pi /4`$ in the present case. Finally it increases and passes from a few km s<sup>-1</sup> to more than 10 km s<sup>-1</sup> in the axial region. Note that the radial velocity vanishes in the equatorial region where it connects to the thin accretion disk, for this solution. This is neither a general result, nor a boundary condition that we impose. Solutions with negative values of the radial velocity can also be obtained. The rotation is maximum at the turning point but it is still dynamically important closer to the disk. The $`\theta `$-component of the velocity shows that the gas is moving towards the axial region in the circulation model while it is directed towards the equator for the infalling envelope. Finally, the maximum of $`V_\theta `$ corresponds to a maximum in the toroidal component of the magnetic field $`B_\varphi `$ which produces a radial “hoop-stress” in the outflow that helps the collimation (Lery & Frank (2000)).

We also present the density iso-contours, with superimposed poloidal velocity vectors, in Fig.4 for the entire model in the poloidal plane. The density in the jet is lower than the density in the molecular outflow by at least a decade in the present example, but this is not a general result of the model. In the molecular outflow the density increases from the outflow to the infalling region, decreases at the interface and reaches again a maximum near the equator. Other examples of density profiles as well as variations of the main quantities along a streamline can be found in Lery et al. (1999) for the circulation model.

We finally present the streamlines obtained with the solution given above in Fig. 5. The streamlines are projected in the poloidal plane and are shown with arrows that indicate the direction of the flow. The outermost region is dominated by the circulation solution except close to the axis and in the equatorial region. On the other hand, the accretion-ejection is dominated by the dynamics near the protostar.

## 5 General properties of the Global Model

### 5.1 The Fluxes

The empirical relationships between the mass, momentum, and energy fluxes provide important observational constraints on physical models of bipolar outflows. Indeed, in our model, the total mass flux $`\dot{M}`$, momentum flux $`\dot{P}`$ and kinetic energy flux $`\dot{E}`$ are respectively given by

$$\dot{M}\left(r/r_o\right)^{1+2\alpha },\dot{P}\left(r/r_o\right)^{1/2+2\alpha },\dot{E}\left(r/r_o\right)^{2\alpha }.$$
(13)
It follows that models with $`\alpha <1/4`$ and $`\alpha >1/4`$, i.e., with loss in energy by radiation through the boundaries, have mass and momentum fluxes dominated by the inner and outer regions respectively. Moreover, the energy flux is dominated by the interior region if $`\alpha <0`$ (and by the exterior for $`\alpha >0`$).

Let us now consider in more detail the mass fluxes in this solution. The total mass flux coming from the molecular cloud is given by $`\dot{M}_{in}=\dot{M}_{infall}+\dot{M}_{acc}`$, where $`\dot{M}_{infall}`$ corresponds to the infalling portion of the circulation pattern, and $`\dot{M}_{acc}`$ to the mass flux in the envelope which is subsequently channeled toward the star through the accretion disk. The total mass loss rate reads $`\dot{M}_{out}=\dot{M}_{MO}+\dot{M}_{jet}`$, $`\dot{M}_{MO}`$ and $`\dot{M}_{jet}`$ being the mass loss rate of the molecular outflow and of the jet, respectively Since the infalling gas is conserved along the path in the circulation model, $`\dot{M}_{infall}=\dot{M}_{MO}`$. Thus, the global mass flux conservation reads

$$\dot{M}_{jet}=\dot{M}_{acc}\dot{M}_{}=f\dot{M}_{acc},$$
(14)
where $`\dot{M}_{}`$ is the accretion rate onto the star, and $`f`$ is the fraction of the accretion flow that goes into the jet.

The mass fluxes are calculated as $`\dot{M}=_{\theta _{in}}^{\theta _{out}}\rho (\theta )V_r(\theta )2\pi R^2𝑑\theta `$, where $`\theta _{in}`$ and $`\theta _{out}`$ are the poloidal angles delimiting each region of the system, and $`R`$ a fiducial distance in the radial direction. We thus consider fluxes through sections of a spherical surface centered on the star, the distance of which from the source is arbitrary as long as the material keep flowing in conical volumes (far enough from the disk). By using this method, the values that we obtain can be tunable within a few orders of magnitude of $`G^{1/2}M^{3/2}r_o^{3/2}`$. For the solution presented here, the mass fluxes in the jet, the molecular outflow and the accreting torus are found to be, respectively:

$$\dot{M}_{jet}6\times 10^7M_{}yr^1,$$
(15)
$$\dot{M}_{infall}=\dot{M}_{MO}8.5\times 10^6M_{}yr^1,$$
(16)
$$\dot{M}_{acc}6\times 10^6M_{}yr^1.$$
(17)
The value of the mass loss rate in the jet is of the same order of magnitude of the ones found in several observational studies, e.g. $`\dot{M}_{jet}\mathrm{3.7\; 10}^7`$ for HH34 (Bacciotti et al.  (1999)), $`\dot{M}_{jet}\mathrm{2.4\; 10}^7`$ for DG Tau (Bacciotti et al.  (2001)). Our value for $`\dot{M}_{MO}`$ is also similar to the ones estimated observationally for typical molecular flows (Bachiller (1996)). We also get $`f=\dot{M}_{jet}/\dot{M}_{acc}0.1`$, for this particular model, again in good agreement with observations and predicted by other MHD models specifically treating the launch of the jet. For example, $`f0.03`$ for self-similar disk-wind models (e.g. Ferreira 1997), and $`f0.3`$ for X-wind models (e.g. Shu et al. 1994).

Our model only applies to the YSO environment from $`r_o10^210^3`$ AU to $`10^410^5`$ AU, which is outside the region where the jet is launched. However, we speculate that a quadrupolar jet could form on a smaller scale. In order to understand the type of acceleration in the jet, it can be useful to look at $`\nu =4\pi \eta v_{\varphi o}/B_o`$, where $`\eta =\rho v_r/B_r`$ is the mass to magnetic flux ratio. In the example presented here, we have $`\nu =0.031`$. This shows that the flow is not accelerated magnetically but centrifugally or thermally. This could be done through series of giant flares close to the central object (Feigelson & Montmerle (1999)), but this is beyond the scope of the present work.

Finally we point out that in our model the poloidal momentum flux, calculated as $`\dot{P}=_{\theta _{in}}^{\theta _{out}}\rho (\theta )V_r^2(\theta )2\pi R^2𝑑\theta `$, is about the same in the jet and in the molecular outflow. This is also in agreement with observations (Bacciotti et al.  (1999)). In our model, however, there is no strong dynamical relationship between the jet and the molecular outflow. Thus it seems that the observational similarity between the momentum rates in such outflows may have no relationship with the mechanism accelerating the molecular lobes. In a future paper, we plan to constrain the global model with observations of specific sources for which the opening angle as well as the mass fluxes can be estimated.

### 5.2 The Energetics of the Flow along Streamlines

The acceleration of a magneto-hydrodynamic wind is usually described in terms of the centrifugal force acting along the field lines. Alternatively, the development of the flow can be described in terms of energies. At the closest point to the source, part of the energy is in the form of magnetic energy, which gets converted (in part) into kinetic energy, as we will see. In this section and in the next ones we will describe and study the variations of the energies, of the plasma beta and of the electric current along streamlines in order to illustrate the physical aspects of the model. We will concentrate, in particular, on the circulation model that is the corner stone of the global model.

In analogy to Bernoulli’s theorem, it is useful for our discussion to define a function $`Be`$ that corresponds to the sum of the specific kinetic energy, the specific potential energy (or gravitational energy), the specific enthalpy of the gas, and the specific magnetic energy, all divided by the local specific gravitational energy to make $`Be`$ dimensionless. In the present case, the gravitational energy per unit mass is approximately the square of the local free-fall velocity. One should note that $`Be`$ is not a streamline integral due to the presence of the magnetic field. Nevertheless, it remains useful for discussing the flow of energy between the various components. This function of $`\theta `$ is given by

$`Be(\theta )`$ $``$ $`(E_k+E_G+E_H+E_B)/(GM/r)`$ (18)
$`=`$ $`(u_r^2+u_\theta ^2+u_\varphi ^2)/21+\mathrm{\Theta }+u^2/(8\pi \mu y^2),`$
where $`u`$, $`\mathrm{\Theta }`$ and $`\mu `$ respectively represent the dimensionless velocity, temperature and density. The variables $`E_k`$, $`E_G`$, $`E_B`$ respectively correspond to the kinetic, gravitational, and magnetic energies, while $`E_H`$ is the enthalpy. The variable $`y`$ is related to the magnetic field by $`Bu/y`$. The parameter $`\alpha `$ is the self-similar exponent, $`M`$ is the mass of the central object and $`r_o`$ is a fiducial radius (See paper I for details). Note that the full energy equation that has $`Be`$ as a first integral has been given previously (See equation 15 in paper I or Henriksen 1997). In the flow, the positivity of $`Be`$ suggests that the gas should reach infinity with a net positive energy, although this is not strict since the energy is not conserved. On the other hand, when $`Be`$ is negative, the gas should not spontaneously escape to infinity. Note that the positivity of $`Be`$ does not imply a lack of conservation of energy, but the excess $`Be`$ represents an energy transfer from hotter regions at smaller radii to motion, i.e., mainly from the enthalpy to the kinetic energy. Of course when proper boundary conditions are applied and the full global problem is solved, the total energy would be conserved. The net gain in $`Be`$ from the infall region to the outflow region in its various components (See Fig. 6) represents the energy transferred to the gas when it is closest to the star.

In Fig. 6, the solutions are for $`\alpha =0.12`$,$`\mathrm{\Theta }_0=0.5`$, $`r_o100`$AU, $`M1M_{}`$. Note that $`Be`$ is smaller than the kinetic energy in this figure because the gravitational energy, not represented here, is negative. Far from the central object, the dominating energy is gravitational at the beginning of the infall. There, repulsing forces such as the pressure barrier, the centrifugal force or the magnetic force are weaker than gravity. Therefore, $`Be`$ starts with a negative value. When the gas gets closer, both the kinetic and magnetic energies increase. Then $`Be`$ becomes positive, which allows the gas to escape in the form of an outflow. It seems most natural that more powerful outflows preferentially form along the rotational axis since there $`Be`$ is most positive in all cases. Thus a bipolar morphology for outflows arises quite naturally in the present scenario. It is worth comparing this result with the zero-pressure limit (Lery et al.  (1999)) where $`Be`$ is conserved on streamlines. In this case, the model does not produce velocities greater than the escape velocity anywhere. This 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 in the outflow phase.

The most rapidly outflowing gas is always near the symmetry axis because these streamlines pass closest to the star, deeper into the gravitational potential well. Also, the material on these streamlines is heated the most vigorously. We note that we do not include radiation pressure in our model. In the present context, the role of the radiation produced by the source is to transport outward thermal energy. The radiation is eventually absorbed contributing to the gas heating.

As the gas gets closer to the source, it rotates faster. The infalling plasma therefore has a larger electric current driven by the rotational motion. This increases the magnetic energy, which is eventually converted into kinetic energy as the gas is redirected outwards. The magnetic field acts to collimate and accelerate the gas towards the polar regions. There the flow presents a strong poloidal velocity and a very low magnetic energy.

In Fig. 7, we show the ratio of the final to initial magnitudes of the kinetic energy in the flow as a function of the initial angle of infall. The solution has been obtained for the same set of parameters as in Fig. 6. Since the gas falls from the equatorial regions and ends up in the axial zone, the ratio gives the gain in energy acquired during the circulation. The gain in kinetic energy varies drastically. It goes from less than unity in the region far above the disk ($`\theta _{in}45^{}`$) to more than 10 for streamlines that start close to the equatorial region ($`\theta _{in}10^{}`$). In the latter case, if the gas starts with an infalling velocity around 5 km s<sup>-1</sup>, it will acquire an escaping velocity of at least 50 km s<sup>-1</sup> (Bachiller & Tafalla (1999)). Again, the circulation is thus very efficient in accelerating the flow for the streamlines that start the closest to the equator, and then travel the closest to the central protostar, going deeper in the potential well.

Therefore, in the axial region, the velocity will appear to decrease as a function of the distance from the axis. Such properties has been observed for collimated molecular outflows (e.g. Guilloteau et al. 1997), and more recently, also for the central optical jet (Bacciotti et al.  (2000)). On the right panel of Fig. 7, the kinetic energy is represented at the turning point as a function of the initial angle of infall. The energy is maximum for an intermediate range of angles, between $`25^{}`$ and $`35^{}`$.

Summarizing, the material ejected in the outflow acquires a positive energy due to a combination of heating by the central protostar, pressure gradients developed by the infalling material, rotation and magnetic forces due to the non-zero Poynting flux. The study of the energetics along a streamline helps to describe how the local infalling flow is channeled into a bipolar outflow.

### 5.3 The Electric Current and the Plasma $`\beta `$ Parameter

The electric current is a major indicator of the collimation of outflows. Indeed, it has been shown that the asymptotic shape of steady axisymmetric magnetized outflows is either paraboloidal or cylindrical according to whether the electric current carried at infinity vanishes or not (See Lery, Heyvaerts, Appl & Norman 1999 and reference therein). In the present model, the axial electric current carried by the flow is given by

$$I\left(r/r_o\right)^{\alpha +1/4},$$
(19)
with $`1/4\alpha >1/2`$. If $`\alpha <1/4`$ (i.e., solutions with a gain in energy by radiation), the current is only large close to the protostar and will vanish at infinity. Thus the outflow should have a paraboloidal or conical shape. On the other hand, for the special case $`\alpha =1/4`$, i.e., solutions with radiative balance, the current is not a function of the distance to the central object but only a function of the angle. The current will have a finite value at infinity and therefore the flow should focus cylindrically. A similar conclusion is obtained for $`\alpha >1/4`$ where the current increases with the distance. Consequently, for this range, the outflow should collimate cylindrically far from the source, at least in the axial regions.

It is also interesting to study the plasma $`\beta `$ parameter that corresponds to the gas to magnetic pressure ratio. For our model, $`\beta `$ is given by $`\beta =y_\varphi ^2P/u_\varphi ^2`$. This quantity varies drastically with the angle in the circulation model as shown in Fig. 8. $`\beta `$ starts from 10 in the infall region, reduces to $`10^2`$ at the turning point and finally increases up to several hundreds in the jet region.

### 5.4 Down to the Thin Disk

The last point that we shall discuss in this section concerns the boundary conditions at the interface with the thin disk, that in the model, is regarded as an open boundary. Physically, this corresponds to the entry of the gas into this region, while it is still falling onto the central object. The gas is not effectively falling onto the thin disk since self-gravity is not taken into account in the model. We have represented the shape of the poloidal component of the magnetic field above the disk in Fig. 9 for several solutions, that differ only by the initial conditions for the the self-similar magnetic variable. We find that the angle of the trajectory of the infalling gas when it goes out of the boundary directly depends on the initial magnetic field. For example, nearly perpendicular field lines can be obtained for a large magnitude of the field. When the magnetic field decreases, the lines bend and present a larger angle w.r.t. the perpendicular to the disk as shown in Fig. 9. In this case, the gas entering the disk has a non-vanishing radial velocity, and the magnetic field has an helical structure.

## 6 Discussion

### 6.1 The Molecular Outflow Acceleration

From an observational point of view, molecular outflows generally consist of a pair of oppositely-directed and poorly collimated lobes, which are symmetrically located about an embedded young stellar object, and detected in broad mm-wave emission lines, especially of carbon monoxide (CO) (Bachiller & Tafalla (1999)). They are usually moderately collimated (with initial opening angles $`\theta _i20^{}90^{}`$), but highly collimated molecular outflows ($`\theta _i<10^{}`$) are also observed (e.g. André et al. 1990). The most luminous YSOs, which generally correspond to the most massive ones, have outflows with relatively wide-opening angles (Churchwell (1997)). Their outflow masses are estimated to be up to several hundred solar masses (e.g. $`200M_{}`$ for Mon R2 (Tafalla et al.  (1997))), and in general they can be more than 5 times larger than the final mass of the central star.

In the present model, most of the properties of molecular outflows can be understood in terms of the jet and the outflow having separate origins, although they may interact at their interfaces. Observationally, the spatial coincidence of shocked molecular hydrogen bows with peaks in the CO outflow emission suggest that the jet and outflow interact (e.g. Davis & Eislöffel 1995). In our model, the jet comes from the inner-most part of the central object, i.e., the inner accretion disk and the forming star, while the molecular outflow is produced by the circulation and is tied to the physics of the molecular cloud. We self-consistently model the infall, outflow, and jet regions using the same set of equations around the same central object. Even though the three different parts of the model are formally disconnected, they are related through strict boundary conditions at their interface surfaces.

This has interesting consequences concerning the subsequent interaction of the flow with the coaxial jet. Since in our case, the difference in velocities between the jet and the molecular outflow material are reduced from the start, the shocks in the zone of acceleration due to their interaction should be less strong. The central fast jet has still the largest part of the total momentum per unit area, and the molecular outflow could undergo a prompt entrainment from the head of the jet. The effects of this entrainment, however, are reduced because of the motion induced in the CO outflow by the circulation around the protostar.

But the most interesting feature of the circulation model is probably that it could produce solutions where the mass of the molecular outflow is larger than the final mass of the forming star. This would particularly true if self-gravity were included in the model. One might then understand how bipolar outflows from massive protostars are observed to transport masses largely exceeding those of the associated stars. For high luminosity objects, solutions can be obtained that show well-collimated jets, even though the CO flow may appear poorly collimated (Bachiller & Tafalla (1999)).

Other important features of the model are that the distribution of the flow mass with velocity approximately follows a power law, and the flow collimation increases with velocity. The model also shows solutions in which the flow velocity increases with distance along a line parallel to the outflow axis. Also, the synthetic CO spectra corresponding to these solutions nicely reproduce most of the observational features of molecular outflows (Fiege & Henriksen (1996); Lery et al.  (1999)). The acceleration is relatively gentle in this type of model, so that molecules are advected with the flow with relatively little destruction. This is in contrast to the more severe acceleration that occurs when gas is entrained by a fast jet, which tends to destroy the molecules (Tafalla et al.  (1997); Downes & Ray (1999)). In our model, several tens of solar masses are moving at supersonic velocities while remaining neutral and in molecular form. The calculations produce solutions where the outflow can have large opening angles, and where the most massive protostars produce the fastest and the less collimated outflows (Lery et al.  (1999)).

In conclusion, we suggest that CO outflows are dominated by the global circulation of material around the protostar, except for in a thin layer surrounding the jet, where the dynamics is governed by entrainment. We stress that, in the present model, the two flows (atomic and molecular) are not strongly linked dynamically, and, hence, there is no need to transfer large momentum from the jet to the molecular outflow through the entrainment processes.

### 6.2 Molecular Cavities

At CO mm-wavelengths, e.g. L1551 (Moriarty-Schieven & Snell (1988); Hartigan et al.  (2000)), NGC 2071 (Moriarty-Schieven et al.  (1989)), or Mon R2 (Tafalla et al.  (1997)) and also in the infrared band (Close et al. 1997a ), cavities inside molecular outflows are observed. The edges of such cavities, located between the jet and the molecular outflow, are delineated by $`H_2`$ emission. Indeed, close to the source, “throats” are observed 100 to 200 AU above the accretion disk (Close et al. 1997b ). The association of outflows with cavities is not limited to low-luminosity sources (Campbell et al.  (1986); Tafalla et al.  (1997)). For example, the blue lobe of AFGL 490 (Campbell et al.  (1986)) coincides with a conical nebula of scattered light, suggesting that the outflow interacts with a shell of dense gas surrounding a partly evacuated channel.

The cavities can have large opening angle in the range 90-125 (e.g. Velusamy & Langer 1998), and their nature is unclear (Chandler et al.  (1996); Gomez et al.  (1997); Gueth et al.  (1998)). These observations are in good agreement with the opening angles in our model and in wide-angle wind models, which produce wider cavities than jet models. It has also been proposed that the cavities could be produced by the wandering of a precessing jet. However, it appears that the dynamical time-scale of the jet is too short for this explanation to work in some observed sources (Chandler et al.  (1996)).

In the present model, the molecular cavities may be identified with the substantial decrease in density in the intermediate region between the jet and the molecular outflow (See Fig.4 in the axial region). Note that solutions with wider opening angles than the one presented here can easily be obtained. Moreover, the central jet region may be comprised mainly of atomic gas, which is not observable in molecular lines. This is likely because the gas will probably be dissociated and ionized near the protostar. As the gas escapes and cools, atomic gas will form rapidly. However, the timescale for the production of molecules is much longer. Thus, the molecular outflows appear as a hollow conical structure (Lery et al.  (1999)). The model indicates a wide range of opening angles from 40 to 140, depending on the input parameters, in agreement with the observed values. However, this angle may vary with time but the timescale for the widening is expected to be much longer than the typical dynamical timescale. Thus our model may represent a quasi-steady flow pattern at some instant in a slowly changing outflow. Anyhow, the cavity is a result of the circulation itself and we do not need to assume that the jet is precessing.

### 6.3 Asymmetric Outflows

Another important observational issue is the large number of bipolar optical outflows showing properties that differ substantially in the blue-shifted and red-shifted parts (Hirth et al.  (1994); Bacciotti et al.  (1999)). In the axisymmetric quadrupolar circulation model, solutions are not limited to an hemisphere, and consequently can differ from side to side of the disk (Fiege & Henriksen (1996); Lery et al.  (1999)). . We present an example of such a solution in Fig. 10.

In this figure, the streamlines are represented for all the three different parts of the model, i.e., the circulation model, the jet model and the infall envelope. The solutions have the same set of parameters on both sides of the equator, i.e., $`\alpha =0.12`$,$`\mathrm{\Theta }_0=0.35`$, $`r_o100`$AU, $`M1M_{}`$, but have different opening angles, that we prescribe as input parameters ($`\theta _{max}=2\pi /3>\pi /2`$ for the upper part, and $`\theta _{max}=\pi /3<\pi /2`$ for the lower part). These solutions may be generated by different conditions on both sides of the disk, for example in a region of a cloud that had a strong pre-existing density gradient, or at the edge of a molecular cloud. The latter case is similar to the standard interpretation of the asymmetric objects as an outflow breaking out of the cloud on one side. Note, however, that contrary to the optical jet case, for the molecular outflow the weak lobe will be pointing out of the cloud as observed for many unipolar CO outflows known, such as NGC2024-FIR5 (Richer et al.  (1992)), or the almost-unipolar HH46-47 system (Chernin & Masson (1991)). Our model also predicts that the dominant outflow lobe will be on the side pointing into the cloud because that is the side where material can be gathered. Thus, the various asymmetries observed on opposite sides of many bipolar outflows, could find a natural explanation in the quadrupolar geometry. In any event, detailed comparisons of solutions for a specific asymmetric source could be a good test of our hypothesis. It is also noteworthy that purely unipolar models had already been found in the context of the circulation model when treated without Poynting flux (Fiege (1996)).

### 6.4 Collimation

Several mechanisms, from purely hydrodynamic (Mellema & Frank (1997)) to magnetic ones, have been proposed to account for the collimation of jets. Magnetic processes seem very promising and they are generally accepted as the primary option. In this case, the outflow is collimated by the pinch effect of the toroidal component of the magnetic field that arises as a consequence of the rotation of the disk (Lery et al.  (1999); Breitmoser & Camenzind (2000)). In cylindrical jet asymptotics, the outflowing plasma reaches an equilibrium wherein it is confined by magnetic forces or gas pressure gradients, while it is supported by centrifugal forces or gas pressure gradients. The idea that the external medium has, in fact, a non-negligible effect has been raised recently (Okamoto (1999); Lery & Frank (2000)). Indeed, magnetized winds may not cylindrically collimate without an external help, such as the channeling effects of thick accretion discs and/or confinement resulting from the ambient medium. Moreover, in spite of the success of magnetic models, recent numerical studies have shown that pure hydrodynamic collimation could be effective at producing jets (Icke et al.  (1992); Mellema & Frank (1997)), thanks to the “shock focusing” mechanism due to the interaction of the outflow with the external environment. In the present model, the collimation of the axial jet arises to a large extent from the hoop stress of the toroidal field component but it is also reinforced by the density gradients in the surrounding molecular outflow. Thus, both hydrodynamic and magnetic effects play a role. Note that the present model predicts that the collimation does tend to increase with the velocity of the gas regardless of the mass of the central object (Fiege & Henriksen (1996); Lery et al.  (1999)).

The presence of low density polar regions may also give a preferential direction and help the alignment of jets coming from a binary system, as in the case of L1551. Here, HST imaging and spectroscopy from the ground (Fridlund & Liseau (1998); Reipurth et al.  (2000)) indicates the presence of two individual jets almost aligned (angle between the jets less than 20) and lying inside a CO cavity.

### 6.5 Possible origin of the quadrupolar topology

Most of current outflow models are based on a dipolar geometry where the magnetic field lines are either advected via accretion towards the central object in the accretion disk or are created in situ. Aburihan et al. (2001) proposed a simple, qualitative mechanism to generate a quadrupolar field topology, which we briefly summarize here. The basic idea is that advection deforms the field into a predominantly quadrupolar topology as the magnetized gas collapses. During the earliest stages of formation, the magnetic field near the protostar would be mainly dipolar as the gas starts to collapse. Since the magnetic field is frozen into the gas, the infalling material advects the magnetic field toward the accreting core until it is strongly pinched toward the protostar in the equatorial region. At this point, the field is mainly radial. Poloidal pressure gradients, with both radial and angular components, develop as the protostar forms. Eventually these pressure gradients become strong enough to deflect some of the gas toward the polar regions (Lery et al.  (1999)). The magnetic field near the protostar is advected with this outflowing gas, and eventually forms magnetic arches near the polar axis. A predominantly quadrupolar field topology remains as these magnetic arches are advected to large radii. Other mechanisms might exist to convert an initially dipolar field into a quadrupolar geometry, but these might rely on reconnection or dynamos.

## 7 Conclusion

We have proposed a new global MHD model for flows around YSOs. The global model combines a jet model, a circulation model and an infalling envelope model. Instead of the usual mechanisms invoked for the origin of molecular outflows, the outflow is powered by the infalling matter through a heated quadrupolar circulation pattern around the central object. The molecular outflow may still be affected by entrainment from the fast jet, but this would be limited to the polar regions and it would not be the dominant factor behind its acceleration.

The solutions are developed within the context of $`r`$-self-similarity. We use the set of steady axisymmetric MHD equations for all the three parts of the model. The only physical scales that enter into our calculation are the gravitational constant $`G`$, the fixed central mass $`M`$, and a fiducial radius $`r_o`$. The parameters of the model are the indexes $`\alpha `$ and $`\alpha _f`$ of the self-similar system, and $`\theta _{min}`$ and $`\theta _{max}`$, the limiting angles between the three regions of the model. Magnetic field and streamlines are required to be quadrupolar in the poloidal plane for the circulation model.

We have explored the physical behaviour of the model by studying the flow through its various energy components along streamlines. We have reported here solutions for the three different parts of this self-similar model. An example of solution is presented for a one solar mass object, even though this class of model can apply to both low and high mass protostars. The solutions show dynamically significant density gradients in the axial region, precisely where the radial velocity and collimation are the largest. Rotation is maximum at the turning point but it is still dynamically important close to the disk.

Our model can reproduce mass loss rates that are in agreement with values derived from observations. For the specific case that is presented here, we obtain $`\dot{M}_{jet}/\dot{M}_{acc}0.1`$, as observed in several systems, even if our model does not deal with the precise jet launching mechanism. Other values, typically ranging from $`0.04`$ to $`0.4`$ are possible, which may impose additional observational constraints on the parameter space of our model.

From an observational point of view, we clarify the nature of the molecular outflow acceleration and its relation with the fast jet, by providing a global picture of the jet/outflow system, which does not primarily rely on entrainment (prompt or turbulent). In this scenario, the observed molecular cavities result from the circulation pattern itself, and there is no need for jet wandering. Our model helps to explain the presence of high mass molecular outflows around massive objects, being the molecular material simply circulating around the central object from the nebula to the bipolar outflows. The various asymmetries observed on opposite sides of many bipolar outflows could find a natural explanation in the quadrupolar geometry; our model predicts that the dominant outflow lobe will be on the side pointing into the cloud. The present study suggest that radiative heating and the Poynting flux may ultimately be the main energy sources driving molecular outflow. Finally, although the details of the jet mechanism may be peculiar to individual objects we believe the infall-outflow circulation to arise naturally given accretion, and thus could also be present in AGN’s and around suitably placed objects such as SS433.

###### Acknowledgements.

TL and AF acknowledge support from NSF Grant AST-0978765 and by the University of Rochester’s Laboratory for Laser Energetics. JDF acknowledges support from CITA and NSERC fellowships.


format: text

text:
# Multifragmentation in central collisions at intermediate energies: a fast process?

## 1 Introduction

In central nuclear collisions at intermediate energies, multifragmentation, defined as the emission in a short time scale of several species of atomic number larger than two, is an important process as compared to other decay mechanisms such as the formation of heavy residues or fission (see for instance ). Such reactions are thus believed to be the ideal tool to study the transition from a liquid-like state (nuclei at normal density) towards a gas-like state associated with the vaporisation of the system . In such a picture, nuclear multifragmentation is associated with the excursion of the excited system at low densities inside the coexistence region of the nuclear phase diagram . There, a volume instability, the so-called spinodale decomposition, leads the system to fragmentation. Such a scenario is supported by microscopic transport calculations based on the nuclear Boltzmann equation and its stochastic extensions . These models predict a compression-expansion phase driving the system at low density.

Recent experimental analyses based on nuclear calorimetry have claimed evidence for a liquid-gas phase transition through the study of various signals . Some of these analyses make extensive use of the thermal multifragmentation statistical models  to prove the existence of thermal equilibrium . In such a framework, charge distributions and fragment multiplicities are fitted with help of several parameters characterising the thermalised source. More involved studies are based on ’back-tracing’ techniques and they allow to obtain not only the mean values but also the distributions of the parameters of the source . It is worth noting that, in all cases, not all data are fitted. Indeed, it is very often assumed that a strong pre-equilibrium phase is present and this is why kinematical cuts are applied to the data before comparison with the predictions of the models. In such thermal models, the only source of fluctuations regarding the kinetic energy distributions of the emitted species is the thermal motion inside the thermalised source associated with a possible (usually assumed self-similar) collective motion . Recently, deformation has been introduced to account for the observed anisotropies in the data .

Our aim in the present work is to explore the possibilities of a scenario in which the main hypothesis underlying the statistical models is abandoned. This means that we wish to compare the data (with as less kinematical cuts as possible) with a model taking into account explicitely the entrance channel characteristics of the reaction as well as the initial correlations of the nucleons in momentum space. This is, of course, not the first attempt of this kind. For instance, molecular dynamics models have been used to understand fragmentation in such a scenario . Indeed, all molecular dynamics approaches, whether they are ’classical’ or ’quantal’, conclude that fragmentation is a fast process for which there is no time for the system to build compression or reach equilibrium. Here, we develop, in this spirit, a phenomenological and transparent approach allowing a detailed comparison with the data. Since the statistical models are able to correctly reproduce the ’static’ properties of nuclear fragmentation, we focus our study on the ’kinematical’ properties for both fragments and light particles obtained in central collisions in the Fermi energy range. The data used for the study have been collected by the INDRA collaboration near the GANIL facility . We consider central collisions for Xe+Sn at 32 and 50 Mev/u . This contribution is organised as follows: the model (hereafter called the ELIE event generator) is described in section 2 and a detailed comparison with experimental data is developped in section 3. Conclusions and perspectives are drawn at the end of the paper.

## 2 The ELIE event generator

The main idea of the model is to consider the extreme hypothesis for which the initial momentum distribution of the nucleons has no time to relax before the fragmentation of the system and for which a (small) fraction of the nucleons suffers hard nucleon-nucleon collisions leading them to escape the system as free particles. For the (large) remaining part of the system, only a rearrangement of the nucleons in the final state inside the various emitted products is considered.

As mentioned in the introduction, we are essentially interested in the kinematical properties of multifragmentation. Therefore, the definition of the partition (that is the ’chemical’ composition in terms of complex particles and fragments present in the final state of the reaction) is not determined from ab initio calculations but we rather use the information given by the experimental data. In other words, the mean multiplicity of each species is put ’by hand’. However, since it is assumed that the fragments are excited, the total multiplicity of detected particles is the sum of two components. The first one is associated with the ’fast’ production during the reaction time and the other corresponds to longer time scales and is associated with the evaporation of particles by the fragments. An iteration on the magnitude of these two components is performed until a reasonable agreement with the data is reached. As far as fragments are concerned, only the mean and the second moment of the multiplicity distribution are needed. Indeed, it appears that the charge distribution can be reproduced rather accurately by randomly choosing sequentially the size of the fragments (see Fig.5). Once the partition is fixed, the energetics of the process is considered. The energy balance writes:

$$E_0=E_{ko}+E_{pot}+Q+E_K+E^{}$$
(1)
where $`E_0`$ is the available centre-of-mass energy, $`E_{ko}`$ is the energy corresponding to the few knocked-out nucleons, $`Q`$ is the mass energy balance between the entrance channel and the considered partition, $`E^{}`$ is the excitation in the fragments while $`E_K`$ is the total kinetic energy carried away by the fragments and the light particles. Last, $`E_{pot}`$ is the potential energy associated with the long range Coulomb force acting between the charged species of the partition. As far as light complex particles and fragments are considered, their composition is obtained by randomly chosing nucleons in the two initial Fermi momentum spheres separated by the relative momentum corresponding to the beam energy. The kinetic energy of each species is thus obtained by summing the momenta of all nucleons belonging to the considered complex product. The total kinetic energy obtained by this procedure is such that there is little chance to fullfill energy conservation. Thus, a Metropolis procedure is developped. An exchange of a nucleon (a Metropolis move) between two species chosen randomly is performed and the kinetic energies recalculated. The procedure is iterated until the value $`E_K^{(i)}`$ at step $`i`$ matches the value of $`E_K`$ needed to fulfill Eq.1 with a given accuracy of typically less than a percent. Then, the partition is ready for propagation in order to ’burn’ the excitation energy of the fragments and to transform the potential Coulomb energy into kinetic energy. The initial momenta are given by the procedure described above. The initial positions are obtained by randomly placing the species inside a volume sufficiently large so that nuclei and particles do not overlap: this gives the value of $`E_{pot}`$. The volume in question is not a critical parameter since the key quantity to be compared with the data is not the separate values of $`E_{pot}`$ and $`E_K`$ but rather their sum (which is constrained by Eq.1) (see for instance  for a detailed discussion on this point). Then, the partition is propagated in a similar way as in the SIMON event generator . In particular, the decay in flight of excited species is considered in order to preserve space-time correlations. The generated event is then filtered through the INDRA software filter and stored for analysis.

Apart from the fact that the partition is to a large extent put ’by hand’, there remains essentially two parameters: The amount of fast knocked-out nucleons and the mean excitation energy in each fragment. These two parameters are chosen in order to get the best agreement with the experimental data (see next section). Let us discuss first the question of the excitation energy, $`E^{}`$. Fig.1 shows the kinematical characteristics of the nucleons inside the fragments. The angular distributions are more and more forward-backward peaked as the beam energy increases meaning that the fragments in this approach are more and more deformed and thus more and more excited. Since, in the model, the short range interaction among the particles is only taken into account by means of the Q-values, nucleons do not interact individually and therefore the velocity distributions shown in Fig.1 are not physical. They however show little evolution and thus, it is expected that the excitation energy resulting from the deformation increases slowly with the beam energy.

Anticipating the results of the next section, we have found a good agreement with the data for the following values of $`E^{}`$: 3.5 and 5. MeV/u for resp. 32 and 50 MeV/u. These values are compatible with those found in . The other parameters of the simulation, that is the amount of first chance knocked-out nucleons and the initial multiplicities of light particles and fragments are adjusted to reproduce simultaneously the observed multiplicities and the center-of-mass kinetic energy distribution of the proton as shown in Fig.2. The percentage of knocked-out nucleons increases with the beam energy as it should be since the phase space for in-medium nucleon-nucleon collisions opens gradually. In terms of the fraction of the total number of nucleons, one finds: 2. and 4. $`\%`$ for resp. 32 and 50 MeV/u. The mean and the width of the fragment multiplicity distribution are 6.5, 3 and 9, 3 resp. for 32 and 50 Mev/u.

## 3 Comparison with the data

The selection of the data has been performed using first completeness criteria, requiring that at least 80 $`\%`$ of the total charge and total linear momentum (for charged particles) be detected. This is a necessary condition to perform an event by event analysis. Then, for each event, a sorting is applied by means of two variables. First, a momentum tensor analysis is developped. The diagonalisation of the tensor gives three eigen-values with which several sorting variables may be defined. Here, we use the so-called $`\theta _{flow}`$ angle. This is the angle between the main axis of the tensor and the beam axis. For more details, see . Second, a variable using the characteristics of the light charged particles is built: This is called $`E_{trans12}`$  defined as the sum of the transverse energies of all detected charged particles with charge lower or equal to 2.

Fig.3 shows the correlation between $`\theta _{flow}`$ and $`E_{trans12}`$ for an incident energy of 50 MeV/u. A comparison between the model and the experimental data shows that the model can reproduce the largest values of the $`\theta _{flow}`$ distribution while it misses the low energy part. This is tentatively interpreted by assuming that the low values of the flow angle distribution are associated with non-central collisions for which the binary character of the reaction remains until fragmentation. Indeed, in our model, this possibility is not taken into account due to the lack of geometry of the simulation. In other words, our model in its present shape is only able to treat reactions corresponding to a total overlap of the two partners and is thus unable to treat deep inelatic reactions. From now, we consider only events for which the $`\theta _{flow}`$ distributions of the model and the data coincide. This corresponds to flow angles larger than 30 degrees for 50 MeV/u and 60 degrees for 32 MeV/u.

All forthcoming figures have been normalised to the total number of events fulfilling this condition. Fig.4 shows the atomic number-parallel velocity plot for central Xe+Sn collisions at 50 MeV/u. The general trend of the experimental correlation is reproduced by the simulation as shown in the lower part of the figure. Fragments are essentially emitted at mid-rapidity. However, a closer look at the figures shows that the model underestimates the kinetic energy of the fragments around Z=20. This point will become clearer in the following.

A tentitative ’substraction’ procedure has been applied to the data and is displayed in the upper part of Fig. 4. It consists in subtracting from the experimental data the contribution predicted by the model and associated with central collisions with the normalisation discussed above. It reveals the presence of a rather ’pure’ fast component associated with the remnants of a projectile-like and a target-like in mid-central collisions. However, for fragments with Z lower than 10, a so-called neck emission corresponding to fragments located around mid-velocity is also evidenced. This mechanism has already been studied for the same system but for a different selection of the data . A detailed analysis of neck emission with such a technique is in progress. We now consider some properties of the fragments in Fig. 5 and 6. Due to problems with the software filter in the CM backward direction, we show only kinematical characteristics of particles emitted in the forward part of the angular range in the CM frame. The mean number and kinetic energy of the fragments as a function of the CM angle are rather well reproduced. The multiplicity is somehow overestimated. Here, we would like to point out that the comparison between the model and the experiment is made without kinematical cuts: In other words, the whole event is considered. This induces strong constraints on the model since any departure between the data and the model for some observables has strong feed-back effects on another set of observables. A general ’fitting’ procedure such as a back-tracing analysis should be envisaged in this context. The kinetic energy distribution is nicely reproduced until large values showing the important role played by the internal motion of the nucleons inside the two partners of the reaction. Same conclusions can be drawn concerning the results at 32 MeV/u (Fig.6). However, here, a sizeable deviation in the angular distribution is observed at small angles indicating the presence of a contribution which is not accounted for in the model. We will come back to that point later.

We now discuss the characteristics of light charged particles. We recall that most of them are produced by picking randomly nucleons from the initial two Fermi spheres. However, some of them are evaporated by the excited fragments on longer time scales. Both contributions lead to anisotropic angular distributions as shown in Fig.7. We stress here that even if the final angular distribution of the fragments may be isotropic (as is the case if one considers only events with $`\theta _{flow}`$ larger than 60 degrees), those particles which have been evaporated will show an anisotropic angular distribution. Although the proton distribution is underestimated at very forward angles, the general trend for other particles is well reproduced. A comparison of the mean kinetic energy as a function of the CM emission angle (Fig.8) exhibits a sizeable deviation, in particular for tritons and alphas, between the model and the data. This means that the model is far from perfect and some aspects of the emission process for light particles is missing in our formalism. It is however surprising that the general trends of the angular and kinetic energy distributions of light particles are rather correctly reproduced by the model using our crude ’coalescence’ algorithm.

Let us conclude this brief comparison between ELIE and the experimental data by coming back to the kinematics of the fragments. In Fig.9 and 10, the mean and the width of the CM kinetic energy of the fragments as a function of their atomic number are displayed. At 50 MeV/u, the second moment is very correctly reproduced while the mean kinetic energy is underestimated as far as the charge exceeds values of 15-20. This could have been anticipated looking at the results of Fig.4: The memory of the entrance channel is stronger for the data than for the simulation. In other words, the model (by means of the Metropolis moves described above) allows too much mixing between the nucleons of the target and of the projectile as far as medium and large mass fragments are concerned. This means that an improvment of the model is needed to account for ’geometrical’ effects. In its present version, our approach can only deal with those collisions for which an almost complete overlap of the projectile and the target is achieved. This corresponds to very small cross-sections and obviously does not cover the whole impact parameter range associated with multifragmentation. By considering only the transverse degrees of freedom (lower part of the figures), one gets rid more or less of such difficulties and the model then gives better results. An analysis (not shown here) of the data at larger impact parameters shows that the transverse energies of the clusters are almost independent of the centrality of the collision except for ’controlable’ final state interaction. Our model should be able to reproduce such trends. Work on this point is in progress.

At 32 MeV/u, a large deviation between the data and the model is observed for the widths although a rather correct agreement is obtained for the mean values. We think the deviation is due to an additional mechanism in the data that is not considered in the model: This is the remaining binary character of the collision despite the already severe selection on the impact parameter. Deep inelastic collisions lead to a strong increase of the fluctuations in the kinetic energy distributions of the fragments, a fact that cannot be accounted for in the present version of our model. It is worth noting that no better agreement is achieved when comparing with tranverse quantities. This suggests that the disagreement is not only due to geometrical effects but to collective effects such as, for,instance, the influence of the mean field on the dynamics of the reaction.

## 4 Conclusions and perspectives

We have studied the kinematical characteristics of nuclear products emitted in central highly fragmented nuclear collisions at intermediate energies. Data have been compared with a phenomenological event generator, ELIE, in which it is assumed that the initial momentum distribution of the nucleons has no time to relax before the disassembly of the system. This distribution is thus used to sample the kinetic energies of the various species produced in the reaction under the constraint of conservation laws. It is worth noting that this technique is flexible enough to test rapidly other (possibly thermalised) momentum distributions when the system breaks. As far as partitions are concerned, they are produced ’by hand’ except for the charge distributions which are reproduced by a combinatorial random algorithm, a puzzling result which remains to be understood. A rather good agreement with the data has been obtained. The kinetic energy and angular distributions of both fragments and light charged particles are well reproduced.

In particular, in our approach, both the collective energy and the deformation often claimed as being necessary to reproduce the data in the framework of thermal statistical models are ’naturally’ obtained.

We thus believe that our model is a valuable alternative to thermal statistical approaches based on shape and/or momentum equilibration in the matter at low density. In our approach, most light particles as well as fragments are produced rapidly at nearly normal density. As such, they keep a strong memory of the entrance channel. In particular, fragments emerge from the reaction highly deformed. Then, they decay sequentially to reach their ground state by evaporation on longer time scales by a rearrangement and emission of nucleons and possibly other complex particles. The strong memory of the entrance channel and the importance of the internal motion of the nucleons inside the two partners of the reaction are such that, in our picture, the origin of kinetic energy fluctuations is to a large extent non-thermal. We would like to point out that this picture is very similar to the one proposed in .

Deviations between the model and the experimental data are expected as far as the incident energy increases. Work along this line is in progress. Indeed, nuclear transparency is expected to decrease as the phase space for in-medium nucleon-nucleon collisions increases. Therefore, the two Fermi distributions should be progressively more and more relaxed on shorter and shorter time scales as the incident energy increases. In the near future, we also plan to extend the model to asymetric systems and to non-central collisions.


text:
# Construction of Orthonormal Bases from the Vectors of the Problem in Minkowski SpaceContribution to the Proceedings of the VIth International School-seminar “Actual Problems of High Energy Physics”. Gomel, Belarus, August 7-16, 2001

## 1 Introduction

Nowadays, methods of direct calculation of the amplitudes of processes are intensively developed (see e.g.  and references there). In this connection the task of construction of orthonormal bases from the vectors of the problem (in particular, polarization bases for vector bosons) becomes actual.

In the present work, we propose to use the modified Gram – Schmidt orthonormalization process in Minkowski space for these purposes.

We use Feynman metrics:
$`\mu =0,1,2,3,a^\mu =(a_0,\stackrel{}{a}),a_\mu =(a_0,\stackrel{}{a}),ab=a_\mu b^\mu =a_0b_0\stackrel{}{a}\stackrel{}{b}`$;
the sign of the Levi-Civita tensor is determined as $`\epsilon _{0123}=+1`$.

Moreover, the generalized Gram determinants are used (see also ). For example:

$$\begin{array}{c}G\left(\begin{array}{cc}_\mu & _\nu \\ _\alpha & _\beta \end{array}\right)=\left|\begin{array}{cc}g_{\mu \alpha }& g_{\mu \beta }\\ g_{\nu \alpha }& g_{\nu \beta }\end{array}\right|=g_{\mu \alpha }g_{\nu \beta }g_{\mu \beta }g_{\nu \alpha },\hfill \end{array}g_{\mu \alpha }=\{\begin{array}{cc}\hfill 1,& \mu =\alpha =0\hfill \\ \hfill 1,& \mu =\alpha =1,2,3\hfill \\ \hfill 0,& \mu \alpha \hfill \end{array}$$

$`\begin{array}{c}G\left(\begin{array}{cc}a& b\\ c& _\beta \end{array}\right)=G\left(\begin{array}{cc}_\mu & _\nu \\ _\alpha & _\beta \end{array}\right)a^\mu b^\nu c^\alpha =\left|\begin{array}{cc}(ac)& a_\beta \\ (bc)& b_\beta \end{array}\right|=(ac)b_\beta (bc)a_\beta \hfill \end{array}`$, and so on.

## 2 Gram – Schmidt orthonormalization process in the $`n`$-dimensional space

We consider briefly the Gram – Schmidt orthonormalization process, with the help of which one can construct an orthonormal basis from the vectors of the problem (see e.g. ).

Let a system of linearly independent vectors $`\{x_1,x_2,\mathrm{},x_n\}`$ is in the $`n`$-dimensional space. Then the following system of vectors is orthogonal:

$$y_1=x_1,$$

$$y_2=x_2\frac{(x_2y_1)}{(y_1y_1)}y_1,$$

$$y_3=x_3\frac{(x_3y_1)}{(y_1y_1)}y_1\frac{(x_3y_2)}{(y_2y_2)}y_2,$$

$$$$

$$y_n=x_n\frac{(x_ny_1)}{(y_1y_1)}y_1\frac{(x_ny_2)}{(y_2y_2)}y_2\mathrm{}\frac{(x_ny_{n1})}{(y_{n1}y_{n1})}y_{n1}.$$

Note that111Here the Greek indices denote the components of vectors, $`\rho =1,2,\mathrm{},n`$

$$\begin{array}{c}(y_k)_\rho =\frac{G\left(\begin{array}{ccccc}x_1& x_2& \mathrm{}& x_{k1}& x_k\\ x_1& x_2& \mathrm{}& x_{k1}& _\rho \end{array}\right)}{G\left(\begin{array}{cccc}x_l& x_2& \mathrm{}& x_{k1}\\ x_1& x_2& \mathrm{}& x_{k1}\end{array}\right)},k=1,2,\mathrm{},n.\hfill \end{array}$$

The orthonormal system is formed by the vectors:

$$\begin{array}{c}(z_k)_\rho =\frac{(y_k)_\rho }{\sqrt{(y_ky_k)}}=\frac{G\left(\begin{array}{ccccc}x_1& x_2& \mathrm{}& x_{k1}& x_k\\ x_1& x_2& \mathrm{}& x_{k1}& _\rho \end{array}\right)}{\left[G\left(\begin{array}{cccc}x_l& x_2& \mathrm{}& x_{k1}\\ x_1& x_2& \mathrm{}& x_{k1}\end{array}\right)G\left(\begin{array}{ccccc}x_l& x_2& \mathrm{}& x_{k1}& x_k\\ x_1& x_2& \mathrm{}& x_{k1}& x_k\end{array}\right)\right]^{1/2}}.\hfill \end{array}$$
(1)
## 3 Modified orthonormalization process in Minkowski space

In the Minkowski space, the Gram – Schmidt orthonormalization process has certain features.

The first, one of the vectors of the basis should be time-like. Therefore, the 4-momentum of a massive particle can be chosen as the first vector (further we also consider the situation where only the massless particles participate in the reaction).

The second, because of some symmetric Gram determinants entering into the denominator (1) are negative (see Appendix), their signs should be taken into account.

Let $`p`$ is an arbitrary 4-momentum, such that $`p^2=m^20`$, $`a`$, $`b`$ and $`c`$ are arbitrary vectors. Using (1), we obtain that four vectors $`l_0`$, $`l_1`$, $`l_2`$, $`l_3`$ form an orthonormal basis:

$$l_0=\frac{p}{m},$$
(2)
$$\begin{array}{c}(l_1)_\rho =\frac{G\left(\begin{array}{cc}p& a\\ p& _\rho \end{array}\right)}{m\left[G\left(\begin{array}{cc}p& a\\ p& a\end{array}\right)\right]^{1/2}}=\frac{m^2a_\rho (pa)p_\rho }{m\sqrt{(pa)^2m^2a^2}},\hfill \end{array}$$
(3)
$$\begin{array}{c}(l_2)_\rho =\frac{G\left(\begin{array}{ccc}p& a& b\\ p& a& _\rho \end{array}\right)}{\left[G\left(\begin{array}{cc}p& a\\ p& a\end{array}\right)G\left(\begin{array}{ccc}p& a& b\\ p& a& b\end{array}\right)\right]^{1/2}}\hfill \\ =\frac{\left[a^2(pb)(pa)(ab)\right]p_\rho +\left[m^2(ab)(pa)(pb)\right]a_\rho +\left[(pa)^2m^2a^2\right]b_\rho }{\sqrt{(pa)^2m^2a^2}\sqrt{2(pa)(pb)(ab)+m^2a^2b^2m^2(ab)^2a^2(pb)^2b^2(pa)^2}},\hfill \end{array}$$
(4)
$$\begin{array}{c}(l_3)_\rho =\frac{G\left(\begin{array}{cccc}p& a& b& c\\ p& a& b& _\rho \end{array}\right)}{\left[G\left(\begin{array}{ccc}p& a& b\\ p& a& b\end{array}\right)G\left(\begin{array}{cccc}p& a& b& c\\ p& a& b& c\end{array}\right)\right]^{1/2}}.\hfill \end{array}$$

Finally we take into account that in the Minkowski space

$$\begin{array}{c}G\left(\begin{array}{cccc}p& a& b& c\\ p& a& b& c\end{array}\right)=\epsilon _{\mu \nu \sigma \tau }p^\mu a^\nu b^\sigma c^\tau \epsilon _{\alpha \beta \lambda \kappa }p^\alpha a^\beta b^\lambda c^\kappa ,\hfill \end{array}$$

$$\begin{array}{c}G\left(\begin{array}{cccc}p& a& b& c\\ p& a& b& _\rho \end{array}\right)=\epsilon _{\mu \nu \sigma \tau }p^\mu a^\nu b^\sigma c^\tau \epsilon _{\rho \alpha \beta \lambda }p^\alpha a^\beta b^\lambda .\hfill \end{array}$$

As result we have

$$\begin{array}{c}(l_3)_\rho =\frac{\epsilon _{\rho \alpha \beta \lambda }p^\alpha a^\beta b^\lambda }{\left[G\left(\begin{array}{ccc}p& a& b\\ p& a& b\end{array}\right)\right]^{1/2}}\hfill \\ =\frac{\epsilon _{\rho \alpha \beta \lambda }p^\alpha a^\beta b^\lambda }{\sqrt{2(pa)(pb)(ab)+m^2a^2b^2m^2(ab)^2a^2(pb)^2b^2(pa)^2}}.\hfill \end{array}$$
(5)
Note that in the final formulae the vector $`c`$ is absent.

Thus, an orthonormal basis in the Minkowski space, using three vectors of the problem (one of which is a 4-momentum of a massive particle) and a total antisymmetric Levi-Civita tensor, can be always constructed.

Construction of a basis for the reaction with massless particles can be performed in a following way. Let $`p_1`$ and $`p_2`$ are $`4`$-momentums of the problem, such that $`p_1^2=p_2^2=0`$. We consider a particular form of the basis (2) – (5) at $`p=p_1+p_2`$, $`m=\sqrt{2(p_1p_2)}`$, $`a=p_2`$:

$$l_0=\frac{p_1+p_2}{\sqrt{2(p_1p_2)}},$$
(6)
$$l_1=\frac{p_1+p_2}{\sqrt{2(p_1p_2)}},$$
(7)
$$(l_2)_\rho =\frac{(p_2b)(p_1)_\rho (p_1b)(p_2)_\rho +(p_1p_2)b_\rho }{\sqrt{2(p_1p_2)(p_1b)(p_2b)b^2(p_1p_2)^2}},$$
(8)
$$(l_3)_\rho =\frac{\epsilon _{\rho \alpha \beta \lambda }p_1^\alpha p_2^\beta b^\lambda }{\sqrt{2(p_1p_2)(p_1b)(p_2b)b^2(p_1p_2)^2}}.$$
(9)
In particular when $`p_1`$ is the 4-momentum of a photon, vectors

$$e_\rho ^\pm (p_1,p_2,b)=\frac{(l_2)_\rho \pm i(l_3)_\rho }{\sqrt{2}}=\frac{Tr\left[(1\pm \gamma _5)\gamma _\rho \widehat{p}_1\widehat{b}\widehat{p}_2\right]}{4\sqrt{2(p_1p_2)}\sqrt{2(p_1b)(p_2b)b^2(p_1p_2)}}$$
(10)
can be used as the polarization basis for this photon. Using the algebra of $`\gamma `$-matrices, it is easy to show that

$$\widehat{e}^\pm (p_1,p_2,b)=\frac{(1\gamma _5)\widehat{p}_1\widehat{b}\widehat{p}_2+(1\pm \gamma _5)\widehat{p}_2\widehat{b}\widehat{p}_1}{2\sqrt{2(p_1p_2)}\sqrt{2(p_1b)(p_2b)b^2(p_1p_2)}}.$$
(11)
Note that replacement of the vector of the problem $`b`$ in (10) by another vector $`c`$ leads to the appearance of the phase factor. Really, using the identity (see e.g. )

$$(1\pm \gamma _5)\widehat{q}Q(1\pm \gamma _5)\widehat{q}=Tr\left[(1\pm \gamma _5)\widehat{q}Q\right](1\pm \gamma _5)\widehat{q},$$

where $`q`$ is an arbitrary massless vector, $`Q`$ is an arbitrary $`4\times 4`$-matrix, we have

$$e_\rho ^\pm (p_1,p_2,b)=e_\rho ^\pm (p_1,p_2,c)e^{\pm i\phi (p_1,p_2,b,c)},$$
(12)
$$e^{\pm i\phi (p_1,p_2,b,c)}=\frac{Tr\left[(1\gamma _5)\widehat{p}_1\widehat{b}\widehat{p}_2\widehat{c}\right]}{4\sqrt{2(p_1b)(p_2b)b^2(p_1p_2)}\sqrt{2(p_1c)(p_2c)c^2(p_1p_2)}}.$$
(13)
Some details of calculations within the considered bases may be found in .

Appendix

Consider the signs of symmetric Gram determinants in denominators of the formulae (2) – (5).

1. $`\begin{array}{c}G\left(\begin{array}{cc}p& a\\ p& a\end{array}\right)=p^2a^2(pa)^2\hfill \end{array}`$, $`p^2=p_{0}^{}{}_{}{}^{2}\stackrel{}{p}^2=m^2`$, $`(p_0>|\stackrel{}{p}|)`$.
When $`a^20`$, it is evident that

$$\begin{array}{c}G\left(\begin{array}{cc}p& a\\ p& a\end{array}\right)0.\hfill \end{array}$$
(14)
Now we consider the situation when $`a^2=a_{0}^{}{}_{}{}^{2}\stackrel{}{a}^2=m_{}^{}{}_{}{}^{2}`$, $`(a_0>|\stackrel{}{a}|)`$:

$$\begin{array}{c}p^2a^2(pa)^2=(p_{0}^{}{}_{}{}^{2}\stackrel{}{p}^2)(a_{0}^{}{}_{}{}^{2}\stackrel{}{a}^2)(p_0a_0|\stackrel{}{p}||\stackrel{}{a}|\mathrm{cos}\varphi )^2\hfill \\ =\stackrel{}{p}^2\stackrel{}{a}^2(1\mathrm{cos}^2\varphi )p_{0}^{}{}_{}{}^{2}\stackrel{}{a}^2a_{0}^{}{}_{}{}^{2}\stackrel{}{p}^2+2p_0a_0|\stackrel{}{p}||\stackrel{}{a}|\mathrm{cos}\varphi \hfill \\ =\stackrel{}{p}^2\stackrel{}{a}^2(1\mathrm{cos}^2\varphi )\left(p_0|\stackrel{}{a}|a_0|\stackrel{}{p}|\right)^22p_0a_0|\stackrel{}{p}||\stackrel{}{a}|(1\mathrm{cos}\varphi )\hfill \\ \stackrel{}{p}^2\stackrel{}{a}^2(1\mathrm{cos}^2\varphi )\left(p_0|\stackrel{}{a}|a_0|\stackrel{}{p}|\right)^22\stackrel{}{p}^2\stackrel{}{a}^2(1\mathrm{cos}\varphi )\hfill \\ =\stackrel{}{p}^2\stackrel{}{a}^2(1\mathrm{cos}\varphi )^2\left(p_0|\stackrel{}{a}|a_0|\stackrel{}{p}|\right)^20.\hfill \end{array}$$

Thus, for any vector $`a`$, (14) is fulfill always.
2. By the direct calculation it can be shown that

$$\begin{array}{c}G\left(\begin{array}{ccc}p& a& b\\ p& a& b\end{array}\right)\hfill \\ =\left|\begin{array}{ccc}p_0& p_x& p_y\\ a_0& a_x& a_y& \\ b_0& b_x& b_y\end{array}\right|^2+\left|\begin{array}{ccc}p_0& p_x& p_z\\ a_0& a_x& a_z& \\ b_0& b_x& b_z\end{array}\right|^2+\left|\begin{array}{ccc}p_0& p_y& p_z\\ a_0& a_y& a_z& \\ b_0& b_y& b_z\end{array}\right|^2\left|\begin{array}{ccc}p_x& p_y& p_z\\ a_x& a_y& a_z& \\ b_x& b_y& b_z\end{array}\right|^2.\hfill \end{array}$$

Under a Lorentz transformation of $`p^\mu =(p_0,p_x,p_y,p_z)`$ into $`(m,0,0,0)`$, the last term reduced to zero. Due to Gram determinant is Lorentz invariant, we have

$$\begin{array}{c}G\left(\begin{array}{ccc}p& a& b\\ p& a& b\end{array}\right)0.\hfill \end{array}$$
(15)
3. $$\begin{array}{c}G\left(\begin{array}{cccc}p& a& b& c\\ p& a& b& c\end{array}\right)=\left|\begin{array}{cccc}p_0& p_x& p_y& p_z\\ a_0& a_x& a_y& a_z& \\ b_0& b_x& b_y& b_z\\ c_0& c_x& c_y& c_z\end{array}\right|^20.\hfill \end{array}$$
(16)

Anyway, as it is well known, the equals signs in (14) – (16) can be only when the vectors entering into the Gram determinants are linearly dependent.