# Velocity asymmetries in young stellar object jets. Intrinsic and extrinsic mechanisms

T. Matsakos, N. Vlahakis, K. Tsinganos, K. Karampelas, C. Sauty, V.

Cayatte, S. P. Matt, S. Massaglia, E. Trussoni, A. Mignone

► To cite this version:

T. Matsakos, N. Vlahakis, K. Tsinganos, K. Karampelas, C. Sauty, et al. Velocity asymmetries in young stellar object jets. Intrinsic and extrinsic mechanisms. Astronomy and Astrophysics - A&A, 2012, 545, 10.1051/0004-6361/201219498. insu-03622054

HAL Id: insu-03622054

https://insu.hal.science/insu-03622054

Submitted on 28 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Velocity asymmetries in young stellar object jets

## Intrinsic and extrinsic mechanisms

T. Matsakos<fnref n="1" />,<fnref n="2" />, N. Vlahakis<fnref n="2" />, K. Tsinganos<fnref n="2" />, K. Karampelas<fnref n="2" />, C. Sauty<fnref n="3" />, V. Cayatte<fnref n="3" />, S. P. Matt<fnref n="4" />, S. Massaglia<fnref n="5" />,
E. Trussoni<fnref n="6" />, and A. Mignone<fnref n="5" />

Received 27 April 2012 / Accepted 30 July 2012

## ABSTRACT

**Context.** It is well established that some YSO jets (e.g. RW Aur) display different propagation speeds between their blue and red shifted parts, a feature possibly associated with the central engine or the environment in which the jet propagates.

**Aims.** To understand the origin of asymmetric YSO jet velocities, we investigate the efficiency of two candidate mechanisms, one based on the intrinsic properties of the system and the other on the role of the external medium. In particular, a parallel or anti-parallel configuration between the protostellar magnetosphere and the disk magnetic field is considered, and the resulting dynamics examined both in an ideal and in a resistive magneto-hydrodynamical (MHD) regime. Moreover, we explore the effects of a potential difference in the pressure of the environment, as a consequence of the nonuniform density distribution of molecular clouds.

**Methods.** Ideal and resistive axisymmetric numerical simulations were carried out for a variety of models, all of which are based on a combination of two analytical solutions, a disk wind and a stellar outflow. The initial two-component jet is modified by either inverting the orientation of its inner magnetic field or imposing a constant surrounding pressure. The velocity profiles are studied by assuming steady flows as well as after strong time variable ejection is incorporated.

**Results.** Discrepancies between the speeds of the two outflows in opposite directions can indeed occur both due to unaligned magnetic fields and different outer pressures. In the former case, the asymmetry appears only on the dependence of the velocity on the cylindrical distance, but the implied observed value is significantly altered when the density distribution is also taken into account. On the other hand, a nonuniform medium collimates the two jets unevenly, directly affecting their propagation speed. A further interesting feature of the pressure-confined outflow simulations is the formation of static knots whose spacing seems to be associated with the ambient pressure.

**Conclusions.** Jet velocity asymmetries are anticipated both when multipolar magnetic moments are present in the star-disk system and when nonuniform environments are considered. The latter is an external mechanism that can easily explain the large timescale of the phenomenon, whereas the former naturally relates it to the YSO intrinsic properties.

**Key words.** magnetohydrodynamics (MHD) – stars: pre-main sequence – stars: winds, outflows – ISM: jets and outflows

## 1. Introduction

Over the past few years, the two-component jet scenario emerges as a strong candidate for describing young stellar object (YSO) outflows. Observational data of classical T Tauri stars (CTTS; Edwards et al. 2006; Kwan et al. 2007) indicate the presence of two genres of winds: one ejected radially out of the central object (e.g. Sauty & Tsinganos 1994; Trussoni et al. 1997; Matt & Pudritz 2005) and the other launched at a constant angle with respect to the disk plane (e.g. Blandford & Payne 1982; Tzeferacos et al. 2009; Salmeron et al. 2011). Consequently, CTTS outflows may be associated with either a stellar or a disk origin, or else with both outflow components making comparable contributions. In addition, such a scenario is supported by theoretical arguments (Ferreira et al. 2006). An extended disk wind is required to explain the observed YSO mass loss rates, whereas a pressure-driven stellar outflow is expected to propagate in the

central region, possibly accounting for the protostellar spin down (Matt & Pudritz 2008; Sauty et al. 2011). Numerical simulations have been also employed recently to investigate the various aspects of two-component jets (Meliani et al. 2006; Fendt 2009; Matsakos et al. 2009, hereafter M09).

Detailed observations of some YSO with bipolar flows have shown peculiar velocity asymmetries between the blue and red shifted regions (e.g. Woitas et al. 2002; Coffey et al. 2004; Perrin et al. 2007). In particular, recent estimates of the RW Aur jet indicate that the speed of the approaching lobe is roughly 50% higher than the receding one (Hartigan & Hillenbrand 2009). Although the above studies suggest that the asymmetry originates close to the source, observations of the same object carried out by Melnikov et al. (2009) point out a similar mass-outflow rate for the two opposite jets and conclude that their different speeds are environmental effects. Thus, it is still an open question whether the discrepancy between

<footnote n="1">¹ CEA, IRAMIS, Service Photons, Atomes et Molécules, 91191 Gif-sur-Yvette, France
e-mail: titos.matsakos@cea.fr</footnote>

<footnote n="2">² IASA & Sect. of Astrophysics, Astronomy and Mechanics, Dept. of Physics, University of Athens, 15784 Zografos, Athens, Greece</footnote>

<footnote n="3">³ LUTH, Observatoire de Paris, UMR 8102 du CNRS, Université Paris Diderot, 92190 Meudon, France</footnote>

<footnote n="4">⁴ Laboratoire AIM Paris-Saclay, CEA/Irfu Université Paris-Diderot CNRS/INSU, 91191 Gif-Sur-Yvette, France</footnote>

<footnote n="5">⁵ Dipartimento di Fisica, Università degli Studi di Torino, via Pietro Giuria 1, 10125 Torino, Italy</footnote>

<footnote n="6">⁶ INAF/Osservatorio Astronomico di Torino, via Osservatorio 20, 10025 Pino Torinese, Italy</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Left panel: a quadrupolar disk field and a dipolar stellar magnetosphere. Right panel: a bipolar disk field and a quadrupolar stellar magnetosphere. The northern hemisphere of either case is referred to as a parallel configuration, whereas the southern as an anti-parallel. The dashed line indicates the surface that separates the disk wind from the stellar outflow. The simulations do not include the central region of the YSO, so we do not attempt to describe the star-disk interaction that takes place inside the gray ellipse. In addition, the requirement of $\nabla \cdot \mathbf{B} = 0$ means that field lines must return along the accretion disk.</figcaption></figure>

the two hemispheres relies on an intrinsic property, such as the magnetic field configuration, or an external factor, such as the physical conditions of the surrounding medium.

Spectropolarimetric measurements of T Tauri stars suggest multipolar magnetospheres, with the dipolar component not always parallel to the rotation axis (Valenti & Johns-Krull 2004; Daou et al. 2006; Yang et al. 2007; Donati & Landstreet 2009). Nonequatorially symmetric field topologies may affect both the way that matter accretes on the protostar (Long et al. 2007, 2008) and the jet launching mechanisms. For jets, it is not clear whether the magnetic field asymmetry persists on comparable timescales to the jet propagation timescales (decades or centuries). In fact, a periodic stellar activity could be directly associated with jet variability. In a similar context, Lovelace et al. (2010) have simulated the central region of YSOs assuming complex magnetic fields. In one interesting case, where the magnetosphere is a combination of dipolar and quadrupolar moments, they find the extreme scenario of one-sided conical outflows being ejected from the star-disk interaction interface. On the other hand, the accretion disk could also have a quadrupolar magnetic field originating either in the dynamo mechanism or advection (Aburihan et al. 2001). In fact, it has been shown that nonbipolar disk field topologies can be a potential source of asymmetry in AGN jets (Wang et al. 1992; Chagelishvili et al. 1996).

The external medium may play a relevant role in the dynamical propagation of the outflow; we refer again to the detailed HST observations of RW Aur by Melnikov et al. (2009), whose findings can be consistent with the effects of inhomogeneities in the environmental conditions. Clumpy and filamentary structures are observed in molecular clouds on a wide range of lengthscales, implying the presence of an external medium that surrounds the jets. In addition, an environment affecting jet propagation may be present but remain almost undetectable in observations, as discussed recently by Teşileanu et al. (2012). Theoretical arguments suggest that density anisotropies in the vicinity of YSOs could collimate YSO outflows and even induce oscillations in the jet’s cross section (Königl 1982). Numerical simulations have investigated the effects of the surrounding gas, such as a collapsing environment whose ram pressure collimates a spherical wind (Delamarter et al. 2000) or an isothermal

medium with a toroidal density distribution that results in a cylindrically shaped outflow (Frank & Noriega-Crespo 1994; Frank & Mellema 1996). Another class of simulations has studied the collimating role of a vertical outer magnetic field, whose magnetic pressure effectively confines an isotropic stellar wind (Matt & Böhm 2003; Matt et al. 2003). We note that, although the external pressure is not thought to be responsible for the jet collimation observed on larger scales, it could still affect it and hence modify the propagation speed.

The goal of the present work is to study the feature of the velocity asymmetry within the context of the two-component jet models presented in M09. We take advantage of both theoretical and numerical approaches (Gracia et al. 2006), by setting a combination of two analytical YSO outflow solutions as initial conditions. Tabulated data of such MHD jets have been derived using the self-similarity assumption (Vlahakis & Tsinganos 1998) and are available for both disk winds and stellar jets. We employ the analytical disk outflow (ADO) solution from Vlahakis et al. (2000) and the analytical stellar outflow (ASO) model from Sauty et al. (2002).

Matsakos et al. (2008) have addressed the topological stability, as well as several physical and numerical properties of each class of solutions separately. In M09 the two complementary outflows were properly mixed inside the computational box, such that the stellar outflow dominated the inner regions and the disk wind the outer. The stability and co-evolution of several dual component jet cases were investigated as functions of the mixing parameters and the enforced time variability. Furthermore, the introduction of flow fluctuations generated shocks, whose large-scale structure demonstrated a strong resemblance to real YSO jet knots.

In this paper, two main scenarios are examined as possible mechanisms for producing velocity asymmetries. The first one refers to intrinsic YSO properties that introduce distinct magnetic field topologies on each side of the equator. The second one refers to external effects arising from differences in the ambient pressure.

In the first case, presented in Fig. 1, each hemisphere is considered to have one of the following two magnetic field configurations: the field lines coming out of the stellar surface
<figure><img src="image_3.png" /><figcaption>Fig. 2. The northern hemisphere assumes a lower pressure, hence a lower degree of collimation, as compared to the southern. The sketch is not in scale, the effects of this mechanism will hold at any distance wherein distinct environments are encountered. The magnetic field topology here is simple, because a dipolar stellar magnetosphere is surrounded by a bipolar disk field.</figcaption></figure>

are parallel (north) or anti-parallel (south) to the large-scale magnetic field. This implies a quadrupolar disk field and a dipolar magnetosphere (left) or, equivalently, a bipolar disk field surrounding a magnetic quadrupole (right).

In the second scenario, shown in Fig. 2, the properties of the medium along the propagation axis are assumed to affect the collimation of the jet. Specifically, the outer pressure could modify the diameter of the outflow and, in turn, would determine the wind speed. This mechanism does not depend on any YSO timescale, and it could appear at any height that the jet encounters heterogeneous environment. We note that, although a thermal pressure of a hot tenuous medium is used in the simulations, this could equivalently represent the pressure due to a cold dense environment, turbulence, or even a large-scale magnetic field. Stute et al. (2008) truncated ADO solutions by imposing a weaker outflow of the same type at outer radii. Here, we extend their study by assuming a constant surrounding pressure instead, also with different values in each jet direction.

All our models focus on large enough scales and avoid the complicated dynamics of the star-disk interaction (e.g. Bessolaz et al. 2008; Zanni et al. 2009; Lovelace et al. 2010). Essentially, these mechanisms are included merely as boundary conditions, an approach frequently adopted in the literature of jet propagation studies (e.g. Ouyed et al. 2003; Fendt 2006; Teşileanu et al. 2009).

Simulations are carried out to study the time evolution, as well as the potential final steady states, and to compare the results obtained between the two sides of the equatorial plane. We are especially looking for emerging asymmetries in the vertical velocity profiles. The resistive MHD regime is of particular interest (Čemeljić et al. 2008) since reconnection plays a key role, mainly at the locations where the magnetic field inverts sign. Furthermore, all cases are re-examined after assuming that the flow has a strong time-variable character, an effect that produces discontinuities along its axis. We do not attempt to model the observed properties of the RW Aur jet, but rather to investigate the intrinsic or extrinsic physical processes that could underlie the velocity asymmetry phenomenon, in a general manner.

The paper is structured as follows. Section 2 describes the jet models and provides information on the numerical setup.

Section 3 presents and discusses the obtained results, while Sect. 4 reports the conclusions of this work.

## 2. Numerical models and setup

### 2.1. MHD equations

The nondimensional resistive MHD equations written for the primitive variables are

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{V}) = 0, \qquad (1) $$

$$ \frac{\partial \mathbf{V}}{\partial t} + (\mathbf{V} \cdot \nabla) \mathbf{V} + \frac{1}{\rho} \mathbf{B} \times \mathbf{J} + \frac{1}{\rho} \nabla P = -\nabla \Phi, \qquad (2) $$

$$ \frac{\partial P}{\partial t} + \mathbf{V} \cdot \nabla P + \Gamma P \nabla \cdot \mathbf{V} = (\Gamma - 1) \eta \mathbf{J}^2 + \Lambda, \qquad (3) $$

$$ \frac{\partial \mathbf{B}}{\partial t} + \nabla \times (\mathbf{B} \times \mathbf{V}) = -\eta \nabla \times \mathbf{J}. \qquad (4) $$

The quantities $\rho$, $P$, and $\mathbf{V}$ are the density, pressure, and velocity, respectively. The magnetic field $\mathbf{B}$, which includes the factor $(4\pi)^{-1/2}$, satisfies the condition $\nabla \cdot \mathbf{B} = 0$. Moreover, the electric current is defined as $\mathbf{J} = \nabla \times \mathbf{B}$, and has also absorbed a factor $(4\pi)^{-1/2}$. The gravitational potential is given by $\Phi = -GM/R$, where $G$ is the gravitational constant, $M$ the mass of the central object and $R$ the spherical radius. The magnetic resistivity tensor has been assumed to be diagonal, with all its elements equal to the scalar quantity $\eta$. Finally, $\Lambda$ represents the volumetric energy gain/loss terms and $\Gamma = 5/3$ is the ratio of the specific heats.

We adopt the axisymmetric cylindrical coordinates $(r, z)$ and we express all quantities in code units. A code variable $U$ is converted to its physical units $U'$ with the help of $U_0$, i.e. $U' = U_0 U$. The constants $U_0$ are specified such that the values of the quantities in the final configuration of the simulations match the typical values of real jets, such as observed velocities and densities. Specifically, $r_0 = z_0 = 1.0$ AU, $\rho_0 = 2.1 \times 10^{-18}$ g cm$^{-3}$, $V_0 = 4.0$ km s$^{-1}$, $P_0 = 3.3 \times 10^{-7}$ dyne cm$^{-2}$, $B_0 = 2.0 \times 10^{-3}$ G, and $t_0 = 1.2$ yr.

### 2.2. Analytical solutions

The employed ADO solution, denoted with subscript D, is appropriate for describing disk winds. It is a steady-state, radially self-similar MHD outflow; namely, each quantity has a certain scaling along the disk radius based on the Keplerian rotation profile. Therefore, if all flow variables are known along a particular fieldline, the spatial distribution of the entire solution can be reconstructed. This type of models have conical critical surfaces and consequently diverge on the axis. However, we substitute the inner region with the ASO solution, denoted with S, which is a meridionally self-similar jet. In this case, the critical surfaces have a spherical shape that effectively models stellar outflows. Some kind of energy input is required at the base of the wind to accelerate the plasma, but we do not treat this region because we simulate high altitudes that correspond to its super Alfvénic domain. Since the ADO solution is described by the polytropic index $\Gamma = 1.05$, we specify the same value in the simulations which implies an energy source term $\Lambda = (5/3 - \Gamma)P\nabla \cdot \mathbf{V}$. We have verified that this does not significantly affect the dynamics (see also Matsakos et al. 2008). In other words, although thermal effects are thought to be crucial for the mass loading close to the disk surface (Ferreira & Casse 2004), we can neglect them above it where the magneto-centrifugal mechanism
takes over. The adopted analytical solutions and their parameters are presented in detail in M09.

## 2.3. Normalization and mixing

The two-component jet model parameters can be classified into two categories. The first one describes the relative normalization of the analytical solutions with the help of the following ratios:

$$ \ell_L = \frac{R_*}{r_*}, \quad \ell_V = \frac{V_{S*}}{V_{D*}}, \quad \ell_B = \frac{B_{S*}}{B_{D*}}, \qquad (5) $$

where $R_*$ is the Alfvénic spherical radius of the ASO model, $r_*$ the cylindrical radius of the Alfvénic surface on a specific field-line of the ADO solution and the star denotes the characteristic values at these locations. The subscripts L, V and B stand for length, velocity and magnetic field, respectively. We specify $\ell_V = 5.96$, such that the two solutions correspond to the same protostellar gravitational field. Following the guidelines of M09 we set $\ell_L = 0.1$. Owing to the nonexistence of a characteristic lengthscale in the ADO solution, the value of this parameter does not play any role in the combination of the two wind components. We also define $\ell_B = 0.01$, in order to provide comparable magnetic fields for the two outflows in the initial conditions of the numerical setup. Higher values would result in a dominant stellar outflow and a negligible disk wind, whereas lower ones in a weak stellar outflow and a strong disk wind.

The second class of parameters controls the combination of the two analytical models. To achieve a smoothly varying magnetic field and ensure the dominance of the protostellar magnetosphere in the central region, we follow an improved mixing process rather than the one presented in M09. The merging depends on the magnetic flux functions $A_D$ and $A_S$, which label the field lines of each initial one-component solution. Once a total magnetic flux $A$ is defined, then the poloidal component of the magnetic field is given by

$$ \mathbf{B}_p = (\nabla A \times \hat{\phi})/r, \qquad (6) $$

where $\hat{\phi}$ is the unit toroidal vector. However, a steep transition from one profile to the other could create artificial gradients in the flux, which would be reflected in the magnetic field. As a result, $B$ could end up being excessively strong or even inverted across the matching surface.

The footpoint of the fieldline that marks the transition, $A_{\text{mix}}$, is therefore chosen at a location where both $A_D$ and $A_S$ have roughly the same slope. Then, a first approximation of the total magnetic flux is defined as $A_1 = A_D + A_S + A_C$, where $A_C$ is a constant that helps normalize $A_D$ and $A_S$. We require that the disk and stellar fields are exponentially damped in the central and outer regions respectively, which leads to an improved approximation for the magnetic flux:

$$ A_2 = \left\{ 1 - \exp\left[-\left(\frac{A_1}{A_{\text{mix}}}\right)^2\right]\right\} A_D + \exp\left[-\left(\frac{A_1}{A_{\text{mix}}}\right)^2\right] A_S. \qquad (7) $$

Finally, with the help of $A_2$, all two-component physical variables, $U$, are initialized using the following Gaussian-type mixing function:

$$ U = \left\{ 1 - \exp\left[-\left(\frac{A_2}{A_{\text{mix}}}\right)^2\right]\right\} U_D + \exp\left[-\left(\frac{A_2}{A_{\text{mix}}}\right)^2\right] U_S. \qquad (8) $$

This relation is also applied to generate the total initial two-component magnetic flux, $A$. The poloidal field is then derived from Eq. (6) and it is divergence free by definition.

Moreover, since accretion and protostellar variability are expected to introduce fluctuations in the ejection, the velocity prescribed on the bottom boundary is multiplied by the following function:

$$ f_S(r, t) = 1 + p \sin\left(\frac{2\pi t}{T_{\text{var}}}\right) \exp\left[-\left(\frac{r}{r_{\text{var}}}\right)^2\right], \qquad (9) $$

where $T_{\text{var}} = 2.5$ is the period of the pulsation and $r_{\text{var}} = 5$ is roughly the cylindrical radius at which the matching surface intersects the lower boundary of the computational box. The fractional variability is set to $p = 0.5$ in the time variable models, in order to check the stability in the presence of strong perturbations and enhance reconnection when resistivity is included. Essentially, Eq. (9) provides a Gaussian spatial distribution of a periodic velocity variability of ~3 yr, which produces shocks and knot-like structures.

## 2.4. Numerical models

The anti-parallel configuration is initialized by assuming a sudden inversion of the magnetic field beyond the fieldline rooted at $r_{\text{inv}} = 5$. In nonresistive MHD, the two systems are entirely equivalent thanks to the invariance of the axisymmetric MHD equations under the following transformations: i) $\mathbf{B} \rightarrow -\mathbf{B}$, ii) $B_{\phi} \rightarrow -B_{\phi}$ and $V_{\phi} \rightarrow -V_{\phi}$. The finite values of the electric conductivity and its dependence on quantities such as $J$ is an unknown factor in YSO outflows. Therefore, we follow a simple approach by applying a constant physical resistivity, with the value $\eta = 1.0$, an approximation adequate to provide a rough estimate of how the system evolves in the presence of nonideal effects. As explained in Sect. 3, this value can be considered low because it ensures the dominance of the advection terms over the resistive ones by several orders of magnitude throughout most parts of the computational box. However, this does not hold true for reconnection separatrices, so its impact cannot be neglected there. Finally, we note that even in the ideal MHD cases, the finite size of the grid cells inevitably introduces artificial reconnection. We minimize its effects by increasing the resolution, and we estimate its influence by comparing results of different mesh refinements.

On the other hand, the truncation of the jet for the low-pressure hemisphere is initialized with the following process. We calculate the total pressure of the outflow, thermal plus magnetic, at the location $r_p = 30$ on the bottom boundary<fnref n="1" />. Then, a constant pressure $P_L$ of the same magnitude is applied on the domain beyond the limiting fieldline rooted at $r_p$. In addition, the density is arbitrarily reduced by two orders of magnitude, whereas the velocity and magnetic fields are set to zero. Although pressure equilibrium holds at the bottom boundary, the weaker field and jet pressure at higher altitudes cannot compensate for the push by the external medium. As a result, a radial collapse is expected during the first steps of the simulation until the horizontal forces are balanced. For the jet propagating in the opposite direction, we assume another simulation case where the truncation is still at $r_p$, but a higher outer pressure is applied with the value $P_H = 2P_L$. Finally, we include a case where the medium has $P_{\text{VH}} = 6P_L$ in order to discuss the morphology of the final steady state configuration.

Table 1 lists the fourteen two-component jet models considered, along with a brief description of their setup. In particular,

<footnote n="1">¹ The plasma-$\beta$ (thermal over magnetic pressure) is approximately equal to 0.03 there.</footnote>
<table><caption>**Table 1.** Numerical models, the name of which follows from the initials of their description.</caption><thead><tr><th>Name</th><th>Variability</th><th>MHD</th><th><i>B</i> orientation</th><th>Outer <i>P</i></th></tr></thead><tbody><tr><td>NIP</td><td>Nonvariable</td><td>Ideal</td><td>Parallel</td><td>—</td></tr><tr><td>NIA</td><td>Nonvariable</td><td>Ideal</td><td>Anti-parallel</td><td>—</td></tr><tr><td>NRP</td><td>Nonvariable</td><td>Resistive</td><td>Parallel</td><td>—</td></tr><tr><td>NRA</td><td>Nonvariable</td><td>Resistive</td><td>Anti-parallel</td><td>—</td></tr><tr><td>VIP</td><td>Variable</td><td>Ideal</td><td>Parallel</td><td>—</td></tr><tr><td>VIA</td><td>Variable</td><td>Ideal</td><td>Anti-parallel</td><td>—</td></tr><tr><td>VRP</td><td>Variable</td><td>Resistive</td><td>Parallel</td><td>—</td></tr><tr><td>VRA</td><td>Variable</td><td>Resistive</td><td>Anti-parallel</td><td>—</td></tr><tr><td>NIPL</td><td>Nonvariable</td><td>Ideal</td><td>Parallel</td><td>Low</td></tr><tr><td>NIPH</td><td>Nonvariable</td><td>Ideal</td><td>Parallel</td><td>High</td></tr><tr><td>NIPVH</td><td>Nonvariable</td><td>Ideal</td><td>Parallel</td><td>Very High</td></tr><tr><td>VIPL</td><td>Variable</td><td>Ideal</td><td>Parallel</td><td>Low</td></tr><tr><td>VIPH</td><td>Variable</td><td>Ideal</td><td>Parallel</td><td>High</td></tr><tr><td>VIPVH</td><td>Variable</td><td>Ideal</td><td>Parallel</td><td>Very High</td></tr></tbody></table>

models NIP and NIA focus on the effects of the parallel and anti-parallel magnetic field configuration in the ideal MHD formulation. Cases NRP and NRA include resistivity in order to understand the role of magnetic reconnection and Ohmic heating, especially in the regions where the magnetic field inverts sign. On the other hand, different values of an external pressure are imposed in models NIPL, NIPH, and NIPVH to investigate the effects of distinct environments. We also consider the time-variable version of all previous cases to examine their behavior in the presence of flow fluctuations, i.e. models VIP, VIA, VRP, VRA, VIPL, VIPH, and VIPVH.

## 2.5. Numerical setup

The simulations are performed with PLUTO<fnref n="2" />, a modular shock-capturing numerical code for computational astrophysics (Mignone et al. 2007, 2012). Cylindrical coordinates are chosen in 2.5D, and second-order accuracy is specified in both time and space. All simulations adopt the TVDLF solver (Lax-Friedrichs scheme), a choice that is essential for retaining stability but comes at the cost of increased numerical diffusivity. Outflow conditions are specified at the top and righthand boundaries, axisymmetric on the left, and all variables are kept fixed to their initial values on the bottom boundary. The size of our domain is $0 \le r \le 100$ in the radial direction and $10 \le z \le 310$ in the vertical. We are interested in radial distances below 50, but a larger box ensures that any boundary effects coming from the imposed zero derivative of $B$ will not influence our results (Zanni et al. 2007). The detailed launching mechanisms of each component are not taken into account since the computational box starts from a finite height. The grid resolution is set to $500 \times 1500$, and the simulations are carried out until $t_{stop} = 30$. We point out that the nonvariable and non-pressure-confined simulations reach an exact steady state before $t = 10$. In M09, we studied the steady state of such two-component jet models on larger timescales and found that the configuration is time-independent to very high accuracy. The pressure-confined models require an order of magnitude longer evolution time to attain a quasi steady outcome, so we specify $t_{stop} = 300$ and a resolution of $300 \times 900$ for efficiency.

## 3. Results

As an overall remark, all nontruncated models reach steady states in an equivalent physical time of a few years. The final

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Logarithmic density contours along with magnetic field lines (left) and vertical velocity contours together with poloidal current lines (right) for the initial conditions of model NIP. The dashed line indicates the surface $A_{mix}$ around which the mixing is assumed, and also marks the polarity reversal for the anti-parallel configurations. The thick solid line shows the separatrix beyond which a constant pressure is initially applied on the pressure-confined models. The initial setup of NIA is identical except for the opposite sign of $B$ and $J$ in the central regions.</figcaption></figure>

configuration of the disk wind remains close to the initial setup, demonstrating the stability of the ADO solution despite its significant modification. On the other hand, the hot stellar component gets compressed around the axis, which also results in an increase in the flow speed there. Even in the cases where the flow is strongly fluctuating, the two-component jet on average retains the properties of the unperturbed model. A steady shock that forms along a somewhat diagonal line at the lower part of the box is found to causally disconnect the jet launching regions from the propagation domain. Although a different mixing has been assumed, all the above features are common to the two-component jet models presented in M09, where a detailed discussion can be found.

### 3.1. Parallel/anti-parallel configurations in ideal MHD

The lefthand panel of Fig. 3 presents the density and the magnetic field lines of the starting setup of NIP, whereas the right-hand panel shows the distribution of the vertical velocity and the poloidal current contours. The initial configuration of NIA has an identical morphology, apart from the reversal of $B$ and $J$ that takes place inside the region marked with a dashed fieldline. Essentially, this surface separates the hot stellar outflow from the cold magneto-centrifugal disk wind, which is supposed to carry most of the mass and angular momentum extracted from the system. Moreover, the pressure-driven flow is faster than the outer parts of the two-component model, a property associated with the normalization and the acceleration efficiency of the specific analytical solutions employed. In relation to the discussion on the mixing function, Fig. 3 demonstrates that there is a restricted choice for the location of the matching surface due to the requirement that both solutions ought to have there a magnetic field of

<footnote n="2">² Freely available at http://plutocode.ph.unito.it</footnote>
<figure><img src="image_5.png" /></figure>

**Fig. 4.** In the left panel, the image shows the logarithmic density and the magnetic field lines for the final steady states reached for models NIP (left hand side; parallel configuration) and NIA (right hand side; anti-parallel configuration). In the right panel, the vertical velocity distribution and electric current contours are shown for the same models. Each panel displays side by side both hemispheres of the YSO system, described in Fig. 1, to compare the steady-states reached by the two models.

a similar shape and magnitude. In other words, the roughly vertical field of the stellar component cannot be matched smoothly with the bent lines coming out of the disk at outer radii. Finally, the thick solid fieldline represents the truncation surface that is imposed in the simulations that study the effects of an external pressure.

Each one of the two panels of Fig. 4 shows side by side the jet evolution of the parallel and anti-parallel configuration described in Fig. 1. In particular, the left pair shows the density and magnetic field lines of the final steady states reached by NIP (left) and NIA (right). The flipping of the field is evident in the anti-parallel case, although the general structure is maintained. From the theoretical point of view of ideal MHD, the magnetic field reversal ($B \rightarrow -B$) of a flux tube would not change the dynamics. However, the thickness of the current sheet that would form on its surface cannot be assumed zero in a physical situation.

The simulation of NIA gives a decreased value of density at inner radii and a peak along the field inversion surface. The latter feature forms simultaneously at all heights during the first steps of the simulation, and it appears more prominent in lower resolutions. In fact, numerical reconnection inevitably appears when a curved current sheet is dragged through a square cell grid. The field and, in particular, the strong toroidal component that dominates the magnetic pressure are destroyed at the interface, giving rise to a perpendicular force. Matter accumulates along that surface, and eventually the thermal pressure compensates the lack of the magnetic one. More accurate solvers and higher resolution would treat the separatrix of the inversion of $B$ in a more consistent way; however, on the one hand, the applied physical resistivity of the simulations presented below is well above the numerical diffusivity, and on the other, “smoothing out” effects are an important stability factor for the initial two-component magnetic field.

<figure><img src="image_6.png" /></figure>

**Fig. 5.** The logarithmic distribution of the magnetic Reynolds number $R_m = rV_z/\eta$ (left panel) and $R_\beta = (P/B^2)R_m$ (right panel) for both YSO hemispheres, i.e. NRP (on the left of each pair) and NRA (on the right of each pair).

The righthand pair of Fig. 4 plots the final vertical velocity distribution and the poloidal currents for the corresponding models shown on the left. A wider radial profile of $V_z$ is observed in NIA than in NIP. Apart from the current sheet that forms in both models along the weak diagonal shock, another one also appears on the surface of the magnetic field reversal, as expected from the previous discussion.

## 3.2. Parallel/anti-parallel configurations in resistive MHD

Based on the analysis of Čemeljić et al. (2008) we consider the magnetic Reynolds number $R_m = rV_z/\eta$ and the quantity $R_\beta = (P/B^2)R_m$, which measure the contribution of the resistive terms in the induction and energy equations, respectively. Figure 5 displays $\log(R_m)$ (left pair) and $\log(R_\beta)$ (right pair) for the final states of NRP (left) and NRA (right). Clearly, the magnetic Reynolds number is large in most parts of the computational domain, i.e. $R_m \gg 1$, indicating that the diffusion of the magnetic field is not of particular importance. However, this definition of $R_m$ does not take into account the magnetic field reversal whose characteristic length is much smaller than $r$. The above argument thus excludes this region. Therefore, approximating resistivity with a low constant value is a useful approach to allow reconnection in our models without significantly affecting the validity of the ideal MHD results of the inner and outer regions.

Moreover, we observe that the condition $10 < R_\beta < 100$ holds in most parts of our domain, which indicates that Ohmic heating is not important either. The value of $R_\beta$ is also ill-defined in the anti-parallel case, since it is based on $R_m$. Nevertheless, $R_\beta$ is large at the axis where the plasma is hot and along the matching surface that the magnetic field reconnects, thanks to the high values of the plasma-$\beta$ there.

Figure 6 shows the distribution of $\log(\rho)$ and the magnetic field lines (left pair) of the final steady states of NRP (left) and NRA (right). The right pair displays the streamlines over the contours of $V_z$ for the same models. Initially,
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Side-by-side logarithmic density contours and magnetic field lines (**left panel**) for the two opposite propagating jets, i.e. NRP (**on the left**) and NRA (**on the right**). Contours of $V_z$ and streamlines (**right panel**) for the same models.</figcaption></figure>

the asymmetric magnetic topologies lead to distinct current distributions. Reconnection is manifested at the regions where the magnetic field inverts sign, driving a different temporal evolution between the two hemispheres. Evidently, the presence of resistive effects in the anti-parallel case results in a spatial readjustment of the velocity profile in the inner regions.

Another observed feature is the reversal of $B_z$ that takes place at the upper and outer domains of both NRP and NRA. This is a steady state configuration that comes from the resistive effects at the base of the outflow and a shear in the velocity profile of the flow. In particular, although the two vector fields, **B** and **V**, are parallel at the bottom boundary, the low $R_m$ values suggest that diffusion dominates advection at the zones right above it. As a consequence, a small angle forms between the poloidal components of **B** and **V** before the magnetic Reynolds number of the flow acquires high values. Shearing effects increase this angle resulting in a negative $B_z$ component.

Figure 7 plots $V_z$ (middle) and the logarithms of $\rho$ (top), and $\rho V_z$ (bottom) along the radial direction at $z = 160$ for models NIP (**left**), NIA (**left**), NRP (**right**) and NRA (**right**). From the parallel setups we deduce that the prescribed value of resistivity has a minor effect on dynamics. In contrast, non-ideal MHD effects amplify the features of the anti-parallel configuration, namely the decrease in the density at inner radii accompanied by an increase in the speed. In addition, the peak observed in $\rho$ along the current sheet is more extended in NRA due to the physical reconnection applied. The mass flux profiles also reflect the asymmetry, despite the fact that $\rho$ and $V_z$ change in opposite directions. Notably, beyond the matching surface all four models demonstrate similar behaviors.

### 3.3. Time variable flows in parallel/anti-parallel configurations

This section investigates the previously presented models when a time variable velocity is applied on the bottom boundary, localized at the stellar component. The simulations are performed until the final time of the corresponding unperturbed cases is

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Radial profile of $\log(\rho)$, $V_z$ and $\log(\rho V_z)$ at $z = 160$ for models NIP (**left; dotted line**), NIA (**left; dashed dotted**), NRP (**right; dotted**), and NRA (**right; dashed dotted**).</figcaption></figure>

reached, during which several shocks propagate throughout the computational domain.

Figure 8 compares the logarithmic density contours and the magnetic field between models VRP (top) and VRA (bottom). The shock front is slightly faster in the anti-parallel configuration and seems to have a larger radius, too. High matter concentrations along the jet axis are observed at substantially different heights; i.e., $z = 160$ for VRP and $z = 210$ for VRA. Furthermore, higher density is also found beyond the matching fieldline of VRA, suggesting that the disk wind is affected as well. In fact, the matter trapped along the current sheet is pushed radially out by the inner shocks that travel across the jet flow.

Figure 9 plots the logarithm of $\rho$ along the $z$ axis for the ideal (bottom) and resistive (top) cases. Density has lower overall values in the anti-parallel topology whereas flow fluctuations alter substantially its vertical profile. In particular, the high concentration regions are located at different heights than in the parallel setup, but without demonstrating any systematic lagging. Both features seem to become stronger when physical resistivity is included. On the one hand, this implies the non-negligible role of numerical diffusivity (required for stability), and on the other, that higher $\eta$ values could lead to significantly different results between the two hemispheres. Since emissivity is proportional to the density squared, it is expected that the above discrepancies could have a direct impact on observations. The pressure approximately follows the behavior of $\rho$, hence is not shown. As a result, the temperature has an almost flat profile for the ideal cases, whereas Ohmic heating slightly increases the average temperature of the resistive models, although without any distinguishable effects among them.

In Fig. 10, which follows the notation of Fig. 9, the vertical velocity $V_z$ (top) is plotted versus the $z$ direction. It was already anticipated from Fig. 7 that the parallel/anti-parallel configurations would show similar speeds. The small displacement of the shock front between the variable resistive models is evident here as well. However, having found discrepancies in the radial profile of $V_z$, as well as in the density distribution, we introduce
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Logarithmic density and field lines for models VRP (top) and VRA (bottom). The two jets correspond to the two hemispheres of a YSO and are shown side by side for comparison. The shock front and matter condensations along the axis are found at different locations between the two cases.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Logarithmic density profiles along the axis for the ideal MHD models (bottom), i.e. NIP (dotted line), NIA (dashed dotted), VIP (solid), VIA (dashed), and the resistive MHD models (top), NRP (dotted), NRA (dashed dotted), VRP (solid), VRA (dashed).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Vertical velocity and weighted average of $V_z$, Eq. (10), along the z axis. The line style of each model follows from Fig. 9. Although the flow speed is the same among all cases, the quantity $\langle V_z \rangle$ does demonstrate asymmetries when nonideal MHD effects are included.</figcaption></figure>

the following weighted average flow velocity, that incorporates emissivity in an approximate way:

$$ \langle V_z \rangle = \frac{\int_{0}^{r_{\text{max}}} \rho^2 V_z 2\pi r dr}{\int_{0}^{r_{\text{max}}} \rho^2 2\pi r dr} \quad (10) $$

This quantity provides a better estimate for the observational implications of the simulations than $V_z$ alone. The integration takes place within the cylinder of radius $r_{max} = 30$, where $\rho$ has its highest values. For altitudes $z > 200$ the density peak along the current sheet is not inside the integrated domain.

The bottom panels of Fig. 10 display the vertical dependence of Eq. (10) for the corresponding models. The ideal MHD regime does not show any significant discrepancies beyond $z = 200$. However, resistive cases do demonstrate velocity asymmetries both in the variable and nonvariable models. The difference is of the order of 30% and 20%, respectively, but since it depends on the specific $\eta$ value assumed, it could be

much higher or even negligible in a real astrophysical situation. In addition, including the current sheet within the integration domain results in a slower mean speed for the anti-parallel model. This can be seen by comparing the regions $z < 200$ and $z > 200$. In other words, the choice of $r_{max}$ also influences the absolute value of the estimated discrepancy, so the resulting percentage should not be taken as a robust number<fnref n="3" />.

Although a parametric study of $\eta$ would clarify the asymmetry dependence on finite conductivity, such a task is beyond the scope of this work. Besides, the actual resistivity value is an unknown factor in YSOs. It suffices that even in the case of $\eta = 1.0$, a low value that does not have significant effects on the parallel

<footnote n="3">³ On the contrary, our aim is qualitative and not quantitative, i.e. to show that velocity asymmetries can occur under some general considerations. A more complex average velocity could have been assumed, but we deliberately choose a simple function to constrain the number of unknown factors that could influence our results.</footnote>
<figure><img src="image_12.png" /></figure>

**Fig. 11.** Logarithmic density and field lines for models NIPL, VIPL, NIPH, VIPH, NIPVH, and VIPVH (from left to right). The *left and middle panels* correspond to the northern and southern hemispheres of Fig. 2, respectively. The *right panel* displays the scenario of a very high external pressure in order to visualize the oscillation of the jet's cross section. Each panel shows the unperturbed model on the left and the time variable counterpart on the right. The outer structure and the condensations along the jet axis are static and unaffected by flow fluctuations.

case, an asymmetry of ~20% to ~30% is found. Nevertheless, we have performed simulations with higher and lower resistivity values that seem to agree with the above conclusions. We also note that the vertical flow velocity is constant along *z*, which indicates that there is no need to simulate larger spatial scales.

### 3.4. External pressure

All simulations of pressure-confined two-component jets reach quasi steady-state configurations within a much longer simulation time than for the models presented above. This is because the initial conditions are farther away from the equilibrium described by the ADO and ASO solutions. Moreover, all final configurations possess static knots along the flow axis as well as a wave-like cylindrical outer surface that recollimates the flow at multiple locations. The separation of these high-density regions seems to increase the lower the ambient pressure applied. Notably, the variability of the inner flow does not seem to disrupt this structure, the shape, or position of its surface. The jet models become asymptotically cylindrical, after having passed from a stage of oscillations in their radius, speed, density, and other physical parameters. Vlahakis & Tsinganos (1997) have shown that, under rather general assumptions, this is a common behavior of magnetized outflows that start noncylindrically before they reach collimation.

The logarithmic density along with the field lines are shown in Fig. 11 for models NIPL, VIPL, NIPH, VIPH, NIPVH, and VIPVH (from left to right). All of them have the same truncation radius, indicated in Fig. 3, but a different external pressure is imposed. Despite the pressure equilibrium holding right above the bottom boundary, the radial ram pressure of the flow expands the jet for small *z*. Eventually it becomes compensated

<figure><img src="image_13.png" /></figure>

**Fig. 12.** Vertical velocity contours for the models shown in Fig. 11, i.e. NIPL and VIPL (**left panel**), NIPH and VIPH (**middle panel**), NIPVH and VIPVH (**right panel**). The velocity maps follow the density distributions.

and then recollimation occurs, compressing the flow at almost its initial diameter. At higher altitudes, the jet radius increases again and the process is repeated creating a sequence of static knots, a configuration that is stable and does not propagate with time. These features smooth out with *z* and the jet acquires a cylindrical shape far away, as seen in model NIPVH. We have verified the stability of the static knots with simulations carried out up to *t* = 500, though in a lower resolution. During time evolution, Kelvin-Helmholtz instabilities appear on the interface with the external medium, but they are gradually suppressed by the time the quasi steady state is reached.

The jet structure of NIPH seems to be a smaller copy of NIPL, and the same holds true for NIPVH with respect to NIPH. This is due to the intrinsic self-similarity property of the ADO solution; namely, all field lines have the same shape but on a different scale. Consequently, a similar outer structure is recovered in all pressure-confined cases, depending on the value of the external pressure imposed. Finally, the mixing of the two analytical outflows introduces a lengthscale in the system, so that the above argument will not hold for arbitrarily small radii.

Figure 12 displays the velocity distribution for models NIPL (**left pair; left**), VIPL (**left pair; right**), NIPH (**middle pair; left**), VIPH (**middle pair; right**), NIPVH (**right pair; left**) and VIPVH (**right pair; right**). As expected, the velocity profile reflects the static oscillations of the jet's cross section since the mass flux is conserved. In addition, its radial dependence consists of layers of decreasing speed. Enforced time variability barely modifies the total distribution demonstrating the stability of the configuration. The low velocity observed in the outer medium is a product of the initial strong transient state in which the outflow is squeezed as the truncation surface (Fig. 3) moves upwards.

The logarithmic density, the vertical velocity and $\langle V_z \rangle$ are plotted along the axis in Fig. 13 for the models NIPH, NIPL, VIPH, and VIPL. We notice that the discrepancies are stronger in the weighted velocity, a quantity possibly relevant for observations. Moreover, $\langle V_z \rangle$ demonstrates a difference between the cases NIPH and VIPH which is relatively smaller than the
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** The quantities log(ρ), Vz, and ⟨Vz⟩ (Eq. (10)) along the z axis, for models NIPH (dotted), NIPL (dashed dotted), VIPH (solid), VIPL (dashed).</figcaption></figure>

variations between the models NIPH and NIPL. This implies that the recollimation structure dominates, producing density peaks much higher than those of the enforced time variability.

Velocity asymmetries are evident in this class of simulations even without accounting emissivity. In particular, the middle plot of Fig. 13 demonstrates a ~10% deviation in Vz between the two cases. The bottom panel suggests that the observed speed of the one jet could be twice the value of the other when weighted with the square of the density, at least locally. In these cases, the integration is calculated for rmax = 20. Therefore, if a YSO outflow finds itself inside a nonuniform environment in which the pressure ratio between the two hemispheres is as low as 2, a significant velocity asymmetry could be measured.

### 3.5. Mass loading

A third possible scenario is that of an increased mass loading at the base of the outflow, such as a long-lived asymmetric accretion or a strong radiation source heating unevenly the disk surface. In that case, the wind density can be different above and below the equator, and so will the acceleration. However, on the one hand, our numerical setup is not fully appropriate to investigating this mechanism in a consistent way, and on the other, the results of such simulations did not provide adequate evidence of any formation of velocity asymmetries. Therefore, we briefly discuss this case without presenting an extensive analysis.

We have carried out simulations assuming the initial setup of NIP, but setting a two times higher density on the bottom boundary. These models allow us to examine the readjustment of the velocity profiles within a higher mass loading regime, naively assuming that all other jet properties remain the same at its base. Despite the significant modification of the boundary conditions, no considerable asymmetries have been found. Regarding the stellar component, we attribute this negative result to the fact that the outflow is already super-Alfvénic. A more consistent way to check how a high plasma mass will affect the velocity is to take the acceleration regions into account. However, this is not a straightforward task since it is not clear

how the energy input required to drive the jet will change in a high-mass loading regime. In addition, the much shorter characteristic time and lengths involved cannot be easily included in our comparably larger scale simulations. On the other hand, although the disk wind does include sub-Alfvénic regions, they appear at large radii of low density and speed that barely contribute to the velocity profile.

To sum up, the large scale two-component jet models presented in this paper, cannot fully capture the case of an asymmetric mass loading. The implicit assumption that the flow on the bottom boundary of each hemisphere has different mass but exactly the same velocity might not be valid if the acceleration mechanisms vary strongly with respect to the ejected mass. Nevertheless, we have verified that in the cases that this holds true, no substantial asymmetries are expected.

## 4. Summary and conclusions

In this paper we address the velocity asymmetries of YSO outflows by investigating two classes of candidate mechanisms that could possibly generate this feature.

The first class depends on the intrinsic properties of the YSO and assumes a parallel and an anti-parallel magnetic field configuration, one for each hemisphere. The application of physical reconnection results in different density and velocity distributions between the two sides, leading to observable speed discrepancies. This scenario can provide asymmetric jets coming from intrinsic physical conditions and processes, even without considering the complex dynamics of the star-disk interaction. The relative questions that arise are whether multipolar fields in the star-disk system can exist and survive on a timescale comparable to that of jet propagation and what the actual value of resistivity in YSO outflows is that could amplify or even suppress the asymmetry phenomenon.

The second class of mechanisms is based on external effects, namely when the YSO resides in a nonuniform environment. Imposing distinct outer pressures at some boundary line of the jet is found to directly affect the degree of collimation of each flow, which in turn results in significantly modified propagation speeds. This mechanism is appropriate for explaining an outflow asymmetry on large timescales. In addition, static knots are found to manifest in the jet structure that are due to multiple recollimation locations. Their stability is robust despite the enforced flow perturbations, whereas their separation is found to be closely associated with the ambient pressure value. There is some evidence for stationary shocks in some systems (Matt & Böhm 2003; Bonito et al. 2011; Schneider et al. 2011), which are possibly associated with the collimation of the flow. However, shocks in protostellar jets are generally observed to propagate with the flow (Reipurth & Bally 2001, and references therein). That standing recollimation shocks, as predicted in non-variable simulations with a confining pressure, are not generally observed in YSO jets implies that either the dynamical evolution of the outflow is not strongly affected by the environment or that recollimation shocks are present but undetected.

In the case of the RW Aur jet, it is not clear enough from the available observations whether the asymmetry is associated with the central engine (Woitas et al. 2002; Hartigan & Hillenbrand 2009) or with environmental effects (Melnikov et al. 2009). Both scenarios could in principle be feasible. We have arbitrarily selected the values of our parameters to explore both candidate mechanisms in a general manner without attempting to model this particular object. To understand the applicability of a particular model to interpreting the velocity asymmetry, more detailed
observations are required, as well as more advanced numerical models, including 3D geometry and radiative processes.

**Acknowledgements.** We would like to thank an anonymous referee for helpful suggestions and the A&A editor M. Walmsley for useful comments that lead to a better presentation of this work. We would also like to thank P. Tzeferacos, O. Teşileanu, and E. M. de Gouveia Dal Pino for fruitful discussions, as well as J.-P. Chièze, F. Thais, C. Stehlé, and E. Audit for their support during the preparation of this paper. This research was supported by a Marie Curie European Reintegration Grant within the 7th European Community Framework Programme, ("TJ-CompTON", PERG05-GA-2009-249164).

## References

Aburihan, M., Fiege, J., Henriksen, R., & Lery, T. 2001, MNRAS, 326, 1217

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Bessolaz, N., Zanni, C., Ferreira, J., Keppens, R., & Bouvier, J. 2008, A&A, 478, 155

Bonito, R., Orlando, S., Miceli, M., et al. 2011, ApJ, 737, 54

Čemeljić, M., Gracia, J., Vlahakis, N., & Tsinganos, K. 2008, MNRAS, 389, 1022

Chagelishvili, G. D., Bodo, G., & Trussoni, E. 1996, A&A, 306, 329

Coffey, D., Bacciotti, F., Woitas, J., Ray, T. P., & Eislöffel, J. 2004, ApJ, 604, 758

Daou, A. G., Johns-Krull, C. M., & Valenti, J. A. 2006, AJ, 131, 520

Delamarter, G., Frank, A., & Hartmann, L. 2000, ApJ, 530, 923

Donati, J.-F., & Landstreet, J. D. 2009, ARA&A, 47, 333

Edwards, S., Fischer, W., Hillenbrand, L., & Kwan, J. 2006, ApJ, 646, 319

Fendt, C. 2006, ApJ, 651, 272

Fendt, C. 2009, ApJ, 692, 346

Ferreira, J., & Casse, F. 2004, ApJ, 601, 139

Ferreira, J., Dougados, C., & Cabrit, S. 2006, A&A, 453, 785

Frank, A., & Mellema, G. 1996, ApJ, 472, 684

Frank, A., & Noriega-Crespo, G. 1994, A&A, 290, 643

Gracia, J., Vlahakis, N., & Tsinganos, K. 2006, MNRAS, 367, 201

Hartigan, P., & Hillenbrand, L. 2009, ApJ, 705, 1388

Königl, A. 1982, ApJ, 261, 115

Kwan, J., Edwards, S., & Fischer, W. 2007, ApJ, 657, 897

Long, M., Romanova, M. M., & Lovelace, R. V. E. 2007, MNRAS, 374, 436

Long, M., Romanova, M. M., & Lovelace, R. V. E. 2008, MNRAS, 386, 1274

Lovelace, R. V. E., Romanova, M. M., Ustyugova, G. V., & Koldoba, A. V. 2010, MNRAS, 408, 2083

Matsakos, T., Tsinganos, K., Vlahakis, N., et al. 2008, A&A, 477, 521

Matsakos, T., Massaglia, S., Trussoni, E., et al. 2009, A&A, 502, 217 (M09)

Matt, S., & Böhm K.-H. 2003, PASP, 115, 334

Matt, S., & Pudritz, R. 2005, ApJ, 632, 135

Matt, S., & Pudritz, R. 2008, ApJ, 678, 1109

Matt, S., Winglee, R., & Böhm, K.-H. 2003, MNRAS, 345, 660

Meliani, Z., Casse, F., & Sauty, C. 2006, A&A, 460, 1

Melnikov, S. Yu, Eislöffel, J., Bacciotti, F., Woitas, J., & Ray, T. P. 2009, A&A, 506, 763

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7

Ouyed, R., Clarke, D., & Pudritz, R. E. 2003, ApJ, 582, 292

Perin, M. D., & Graham, J. R. 2007, ApJ, 670, 499

Reipurth, B., & Bally, J. 2001, ARA&A, 39, 403

Salmeron, R., Königl, A., & Wardle, M. 2011, MNRAS, 412, 1162

Sauty, C., & Tsinganos, K. 1994, A&A, 287, 893

Sauty, C., Trussoni, E., & Tsinganos, K. 2002, A&A, 389, 1068

Sauty, C., Meliani, Z., Lima, J. J. G., et al. 2011, A&A, 533, 46

Schneider, P. C., Günther, H. M., & Schmitt, J. H. M. M. 2011, A&A, 530, 123

Stute, M., Tsinganos, K., Vlahakis, N., Matsakos, T., & Gracia, J. 2008, A&A, 491, 339

Teşileanu, O., Massaglia, S., Mignone, A., Bodo, G., & Bacciotti, F. 2009, A&A, 507, 581

Teşileanu, O., Mignone, A., Massaglia, S., & Bacciotti, F. 2012, ApJ, 746, 96

Trussoni, E., Tsinganos, K., & Sauty, C. 1997, A&A, 325, 1099

Tzeferacos, P., Ferrari, A., Mignone, A., et al. 2009, MNRAS, 400, 820

Valenti, J. A., & Johns-Krull, C. M. 2004, Ap&SS, 292, 619

Vlahakis, N., & Tsinganos, K. 1997, MNRAS, 292, 591

Vlahakis, N., & Tsinganos, K. 1998, MNRAS, 298, 777

Vlahakis, N., Tsinganos, K., Sauty, C., & Trussoni, E. 2000, MNRAS, 318, 417

Wang, J. C. L., Sulkanen, M. E., & Lovelace R. V. E. 1992, ApJ, 390, 46

Woitas, J., Ray, T., Bacciotti, F., Davis, C., & Eislöffel, J. 2002, ApJ, 580, 336

Yang, H., Johns-Krull, C. M., & Valenti, J. A. 2007, AJ, 133, 73

Zanni, C., & Ferreira, J. 2009, A&A, 508, 1117

Zanni, C., Ferrari, A., Rosner, R., Bodo, G., & Massaglia, S. 2007, A&A, 469, 811