Twocomponent jet simulations
II. Combining analytical disk and stellar MHD outflow solutions
T. Matsakos^{1}  S. Massaglia^{1}  E. Trussoni^{2}  K. Tsinganos^{3}  N. Vlahakis^{3}  C. Sauty^{4}  A. Mignone^{1}
1 
Dipartimento di Fisica Generale, Università degli Studi di Torino,
via Pietro Giuria 1, 10125 Torino, Italy
2 
INAF/Osservatorio Astronomico di Torino, via Osservatorio 20,
10025 Pino Torinese, Italy
3 
IASA and Section of Astrophysics, Astronomy and Mechanics,
Department of Physics, University of Athens,
Panepistimiopolis, 15784 Zografos, Athens, Greece
4 
Observatoire de Paris, L.U.Th., 92190 Meudon, France
Received 27 September 2008 / Accepted 22 April 2009
Abstract
Context. Theoretical arguments along with observational data of YSO jets suggest the presence of two steady components: a disk wind type outflow needed to explain the observed high mass loss rates and a stellar wind type outflow probably accounting for the observed stellar spin down. Each component's contribution depends on the intrinsic physical properties of the YSOdisk system and its evolutionary stage.
Aims. The main goal of this paper is to understand some of the basic features of the evolution, interaction and coexistence of the two jet components over a parameter space and when time variability is enforced.
Methods. Having studied separately the numerical evolution of each type of the complementary disk and stellar analytical wind solutions in Paper I of this series, we proceed here to mix together the two models inside the computational box. The evolution in time is performed with the PLUTO code, investigating the dynamics of the twocomponent jets, the modifications each solution undergoes and the potential steady state reached.
Results. The coevolution of the two components, indeed, results in final steady state configurations with the disk wind effectively collimating the inner stellar component. The final outcome stays close to the initial solutions, supporting the validity of the analytical studies. Moreover, a weak shock forms, disconnecting the launching region of both outflows with the propagation domain of the twocomponent jet. On the other hand, several cases are being investigated to identify the role of each twocomponent jet parameter. Time variability is not found to considerably affect the dynamics, thus making all the conclusions robust. However, the flow fluctuations generate shocks, whose large scale structures have a strong resemblance to observed YSO jet knots.
Conclusions. Analytical disk and stellar solutions, even sub modified fast ones, provide a solid foundation to construct twocomponent jet models. Tuning their physical properties along with the twocomponent jet parameters allows a broad class of realistic scenarios to be addressed. The applied flow variability provides very promising perspectives for the comparison of the models with observations.
Key words: ISM: jets and outflows  stars: formation  stars: premain sequence  magnetohydrodynamics (MHD)  methods: numerical  stars: winds, outflows
1 Introduction
Jets are supersonic and highly collimated plasma outflows emanating from a plethora of astrophysical objects. In particular, those associated with Young Stellar Objects (YSO) have been found to be accretion powered (Cabrit et al.1990; Hartigan et al. 1995), to have narrow opening angles and to propagate for several hundreds of AU (Dougados et al. 2000; Hartigan et al. 2004). Although their large scale properties are rather well known, the conditions at the launch regions are still unclear. The new generation of high angular resolution instrumentation is expected to adequately resolve the central regions of YSOs and hence constrain the various theoretical models that currently exist.
A promising scenario supported by both observational data and theoretical arguments is that of a twocomponent jet, wherein a pressure driven stellar outflow is surrounded by a disk wind. In particular, He I profiles of classical T Tauri stars (CTTS) indicate the presence of two genres of wind (Edwards et al. 2006; Kwan et al. 2007). One is ejected radially with respect to the central object and the other is launched at a constant angle with respect to the equatorial plane. As a result, CTTS may be classified according to their outflow properties. Some of them seem to be associated with a stellar origin, others with a disk origin and the rest with both components having roughly equivalent contributions. Therefore, it is suggested that both types of winds participate, with the dominance being dictated by the intrinsic physical factors of the specific YSO.
Such a scenario (e.g. Sauty & Tsinganos 1994; Shu et al. 1994) is supported by theoretical arguments as well. Ferreira et al. (2006) conclude that YSO jets consist of two types of steady winds plus a sporadic outflow. An extended disk wind, which is required for the explanation of the high mass fluxes observed in optical jets and an inner pressure driven outflow of stellar origin (Bogovalov & Tsinganos 2001) collimated by the disk wind. A third component is expected to be launched due to the variable conditions of the thin layer between the protostellar magnetosphere and the disk's magnetic field. Their interaction may drive weak sporadic mass ejections probably associated with jet variability.
In favor of the twocomponent jet scenario, there is also the yet unresolved question of the protostellar spin down. Matt & Pudritz (2005, 2008a, b) have shown that the disklocking mechanism, which was believed to slow down the rotation of the central object, is not in good agreement with observations. On the contrary, they propose that the stellar wind is capable of and most likely responsible for the spin down of the protostar. A wide parameter space has been investigated to support such a conclusion, whereas it is argued that the physical mechanisms which drive the actual launching are less important, hence allowing all sorts of stellar wind models.
A plethora of studies exists in the literature concerning numerical simulations performed to investigate the launching and propagation of jets. Two approaches are adopted: in the one, the disk is treated as a boundary (e.g. Pudritz et al. 2006; Fendt 2006, 2009 and references therein), while in the other the disk is included in the computational box, hence studying its dynamics simultaneously and self consistently with those of the jet (first studied in Casse & Keppens 2002, 2004). More recently, Meliani et al. (2006) effectively incorporated a stellar type outflow accelerated by turbulent heating and in Meliani & Keppens (2007), the transverse stability of relativistic twocomponent jets was examined. Furthermore, adopting a different initial setup, Zanni et al. (2007) studied the effects of resistivity on the dynamics of the diskjet system and Tzeferacos et al. (2009) performed an interesting parameter study on disk magnetization.
Despite the complexity of the nonlinear MHD equations, the derivation of analytical steady state outflow solutions has proved successful in the context of selfsimilarity (Vlahakis & Tsinganos 1998). Each family of these solutions (radially or meridionally selfsimilar) manages to capture the physical mechanisms involved in either disk winds (Blandford & Payne 1982; Ferreira 1997; Vlahakis et al. 2000, hereafter VTST00), or stellar outflows (Sauty & Tsinganos 1994; Trussoni et al. 1997; Sauty et al. 2002, hereafter STT02). The geometrical properties of these two classes of solutions are complementary. Although radially selfsimilar models become singular at small polar angles, the meridionally selfsimilar ones are by definition appropriate for modeling of the outflow at the axis.
In the first paper of this series (Matsakos et al. 2008, Paper I), we addressed the topological stability, as well as several physical and numerical properties, separately for typical radially and meridionally selfsimilar solutions. Such analytically derived wind models were defined as ADO (Analytical Disk Outflow) and ASO (Analytical Stellar Outflow), respectively.
Concerning the ADO model, its main feature is the formation of a shock in the super fast magnetosonic region. Upstream of this shock, the analytical and the asymptotic numerical solutions are basically coincident, while the downstream flow converges to a consistent physical solution, overcoming the singularity of the analytical model at the symmetry axis (first achieved in Gracia et al. 2006). This shock corresponds to the numerically modified fast magnetosonic separatrix surface (FMSS, Tsinganos et al. 1996) that causally disconnects the downstream flow from its launching region. This property is quite robust to variation of the physical parameters, and has been recently confirmed also in Stute et al. (2008), where an outer radial truncation of the disk wind was imposed in the simulations. Moreover, a particular model was initialized by specifying a sub modified fast solution both at the initial conditions and at the boundaries (i.e. a flow causally connected throughout the whole computational domain). Over time, the shock was still found, with its position marking the FMSS that causally separates the upstream and downstream regions.
On the contrary, the ASO model does not show singularities at its boundaries and therefore, the evolution of its super Alfvénic region does not show any readjustments. However, since energy input is a vital constituent of the model's acceleration, the modifications of the energy source terms in the sub Alfvénic domain were demonstrated to strongly affect the outcome of the flow. In particular, we verified that an adiabatic evolution resulted in a collapse of the jet to an almost static atmosphere, whereas specifying a polytropic index to mimic almost isothermal conditions produces a weak collimated turbulent wind.
The goal of the present work is to study the twocomponent jet scenario, taking advantage of both analytical and numerical approaches. Specifically, we construct models by properly defining the initial conditions with a mixture of two analytical (ADO & ASO) solutions connected through a transition region. The introduction of a few normalization and mixing parameters, along with enforced time variability applied to the stellar component or at the matching surface, allows the examination of several interesting cases.
The paper is structured as follows. Section 2 revises a few basic properties of the analytical solutions, Sect. 3 describes the mixing procedure followed to set the initial conditions. In the same section, the different cases investigated are presented along with the numerical setup. In Sect. 4 we discuss the results of the simulations performed. Section 5 summarizes and reports the conclusions of this work.
2 MHD equations and the analytical solutions
Our starting point is the ideal MHD equations for the conservation of mass, momentum, energy and magnetic flux together with the fluxfreezing condition:
where , P, , are the density, pressure, velocity and magnetic field (over ), respectively. The gravitational potential, is equal to with , and R denoting the gravitational constant, the mass of the central object and the spherical radius, respectively. represents the volumetric energy gain/loss terms ( , with the energy source terms per unit mass), and is the ratio of the specific heats.
Assuming steady state conditions and axisymmetry, several conserved quantities exist along the fieldlines (e.g. Tsinganos 1982). These are the mass to magnetic flux ratio , the angular velocity of the footpoints of the fieldlines and the total angular momentum flux to mass flux ratio^{} L. If , where is the polytropic index (see Sect. 3.4), also the total energy flux to mass flux E and the specific entropy Q are conserved along the streamlines.
In the paper we adopt the following notation: subscripts D and S are used to refer to the ADO and ASO solutions, respectively, while and are the cylindrical and spherical coordinates. Note that in Paper I the subscript r was used for the ADO model and for the ASO solution. The subscript * denotes a constant of the order of unity which is used for the relative normalization of the two solutions in order to correspond, for instance, to a solution of the same protostellar mass, as will be explained in Sect. 3. The values of the starred quantities correspond to the nondimensional physical variables at the Alfvénic surfaces of each model at the reference fieldline (see below). Finally, subscript zero in a quantity U_{0} is used to introduce dimensions in the code units U, i.e. U' = U_{0}U, where U' is the physical value of a variable given in cgs.
2.1 The analytical models
We employ the ADO solution which is described in VTST00 and implemented in Paper I, that successfully crosses all three critical surfaces. The ASO model we adopt is a solution similar to the one presented in the first article of this series, taken from STT02, but with different parameter values: higher mass loss rate, larger magnetic lever arm and a non spherically symmetric gas pressure. Here we only provide a few aspects of the analytical solutions, whereas the model parameters are reported in Table 1^{}^{,}^{} and the explicit formulae of the physical variables are provided in Appendix A. Further technical information on the solutions can be found in Paper I and references therein.
Table 1: Parameters characterizing the adopted analytical solutions.
Recalling a few useful expressions, the starred quantities for each analytical
model are related in the following manner:
where r_{*} and R_{*} correspond to the nondimensional distances of the Alfvénic surfaces of the ADO and ASO solutions, respectively, and g = 4 is the constant of the gravitational force in code units.
The magnetic field of each solution is given by the following formula:
labels the isosurfaces that enclose constant poloidal magnetic flux, i.e. the magnetic fieldlines. In particular, for the ADO solution, A is given by:
(8) 
Similarly, for the ASO model:
(9) 
The values of and are provided by the analytical solutions (see Paper I and references therein for more details).
We provide here the measure of the magnetic lever arm, braking the disk or star
for each solution, as defined in Ferreira et al. (2006).
This is the same for all fieldlines, and is given by the relation
:
(10) 
(11) 
where is the cylindrical distance of the footpoint of a particular fieldline and is the cylindrical distance of its Alfvénic point. and correspond to the values of the analytical solutions at the equatorial plane and at the base of the stellar wind, , respectively.
3 The numerical models
In order to choose physical scales, we set the length, and density code units equivalent to r_{0} = 1 AU and g cm^{3}. In addition we assume the protostar to be of one solar mass, . Then, since the MHD Eqs. (1)(4) are written in non dimensional form, it can be easily derived that: km s^{1}, dyne cm^{2} and G. Hence, the time unit corresponds to t_{0} = 0.32 y.
3.1 Normalization
Now, we normalize the solutions to each other by defining the three ratios,
which are parameters of the twocomponent jet models:
(12) 
where the subscripts l, V and B stand for length, velocity and magnetic field, respectively. As it will be seen, only can be chosen freely, while the other two are fixed by physical arguments and the properties of the analytical solutions. More precisely, observations indicate that the launching region of disk winds lies in the range 0.2  3 AU (Bacciotti et al. 2002; Anderson et al. 2003; Coffey et al. 2004). Therefore, demanding that the reference fieldline is rooted at 0.16 AU on the equatorial plane, we find r_{*} = 1. Moreover, assuming that the region where the stellar wind is being launched is roughly at 0.01 AU or at in code units, we derive R_{*} = 0.1and hence . It follows from relations (5) and (6) that , , and thus . Finally, we arbitrarily set and the choice of will control .
3.2 The mixing function
Since the mixing will depend on the magnetic fieldlines, we define a trial
magnetic flux function by the simple sum:
.
We point out that this quantity will help only in the mixing procedure and will
not be used to generate the magnetic field present in the initial conditions.
We further define the mixing function:
with the weights and given by:
(14) 
In the latter expressions, is a constant corresponding to the matching surface rooted at 0.16 r_{*}, i.e. at 0.16 AU on the equatorial plane, q is a parameter that effectively moves this surface closer to the protostar and d sets the steepness of the transition from the inner ASO to the outer ADO solution.
The initial values of the physical variables ,
P, V_{z}, ,
are set up using relation (13).
Moreover, with the help of the same expression, we initialize the twocomponent
magnetic flux function^{} A, from which the poloidal component of the magnetic field is generated
using Eq. (7).
Finally, V_{r} is initialized following the ideal MHD condition, i.e. demanding
that the poloidal magnetic field is parallel to the poloidal velocity:
(15) 
Essentially, such a mixing function provides an exponential damping of each solution around a particular fieldline of the combined magnetic field. Therefore, close to the axis, the ASO model dominates, whereas the ADO becomes the main contributor at the outer regions.
The twocomponent jet numerical models can be constructed by specifying the three normalization parameters, , , , and the three mixing parameters, , q, d. As it has already been explained, , and are given a fixed value, leaving , q and d free to examine a variety of twocomponent scenarios. The latter three parameters control the respective dominance, the location of the matching surface in between the protostardisk region and the steepness of the transition region. Different values in this parameter space may address the various T Tauri outflow types and their evolutionary stage. One would expect that in many cases the efficiency of disk winds would manifest the early phases of the YSOdisk system, whereas stellar winds would eventually dominate, especially after the disk has accreted and during the arrival of the star on the main sequence.
3.3 Time variability
Accretion, which controls the conditions at the base of stellar winds, is not
steady in time but rather varies over different time scales ranging from hours,
to days, months, even years (Alencar & Batalha 2002; Stempels &
Piskunov 2002; Johns & Basri 1995).
On the other hand, the protostar is expected to show some sort of variability as
well, for instance the phenomenon of the 11 yr solar cycle.
Therefore, the introduction of time variability in the inner stellar component
will allow us to study the stability issues of more general and realistic
scenarios.
In order to achieve this we prescribe the following function:
where is the period of the pulsation and is roughly the cylindrical radius at which the matching separatrix intersects the lower boundary z = 10 of the computational box. We enforce a sinusoidal time variability depending on , by multiplying a physical quantity of the lower boundary with f(t). The exponential in Eq. (16) helps to contain the perturbation only at the inner regions, i.e. the stellar component. Note that flow fluctuations induce the formation of knotlike structures. Therefore, also including radiation cooling during the evolution (Tesileanu et al. 2008) would allow a direct comparison with observational data. However, this is left to a future article of this series.
Since it is believed that a sporadic outflow is driven by the stardisk magnetic
interaction (Ferreira et al. 2000; Matt et al. 2002), we examine
such cases as well.
In this case we adopt a similar function:
3.4 Energetics
We set
in
Eq. (3) with
.
This assumption, originally made for the derivation of the ADO solution, is
equivalent to a polytropic relation
along each
fieldline.
Essentially, it represents the adiabatic evolution of a gas with a ratio of
specific heats ,
which corresponds to the following energy conservation
law that is solved over time:
(18) 
Recall that in Paper I, simulations were carried out both for the ADO and the super Alfvenic regions of the ASO solution to test the effects of such an energetic assumption ( ), as well as an isothermal ( ) or an adiabatic one ( ). For each model, these different cases produced almost identical results, thus allowing us to safely adopt this simplification of the energy equation.
3.5 The numerical twocomponent jet models
Table 2: A short description and the parameters of the unperturbed numerical models.
Table 3: The time variable numerical models.
Table 2 lists the unperturbed numerical twocomponent jet models along with their parameters. Table 3 presents those constructed to investigate the stability and structure when time variability is applied at the stellar wind or at the matching surface, effectively mimicking an Xtype wind. In particular, the second column of Table 3 reports the ratio of the periodicity of the enforced fluctuation, , over the Keplerian rotation period ( days) of the protostellar radius, located roughly at 0.01 AU. This means that we address phenomena with time scales associated with accretion and the physical conditions present at the stardisk region. Note that in our models, the Keplerian period of the equatorial footpoint of the matching surface is of the order of 10 days. In the third and fourth columns of Table 2 we indicate the physical quantity that is varied and where it is perturbed, respectively (i.e. adopting Eq. (16) or (17)).
3.6 PLUTO code and the numerical setup
The simulations are performed with PLUTO^{} (Mignone et al. 2007), a versatile shockcapturing numerical code suitable for the solution of highmach number flows. The grid is set up in axisymmetric cylindrical coordinates (2.5D), leaving the study of azimuthal stability for a future work. Second order accuracy is applied in both space and time, and the LaxFriedrichs solver is adopted. However, the choice of a particular solver was not found to influence the results. The condition is ensured with the 8wave formulation.
Figure 1: Logarithms of the density at different times for a typical model (5q02). The time unit is 0.32 y (or 10 Keplerian rotations of the footpoint of the matching surface). In the rightmost plot, a weak shock is observed along the diagonal. In addition, the initial matching surface, approximately at r = 5 on the lower boundary, is still evident. To check this, compare it with the leftmost plot which describes a configuration very close to the initial one. 

Open with DEXTER 
The length code unit is equivalent to 1 AU, and therefore the correspondence to real physical distances is straightforward. We consider a computational box with and for the unperturbed models, and with and for the time variable ones, omitting the acceleration region of the ASO solution. There are two reasons for doing so. The first argument concerns the complexities appearing when the ADO solution is initialized in a computational domain that approaches the origin in cylindrical coordinate systems (as demonstrated in Paper I). Second, the ASO solution provides a time independent energy source term, which if included, will artificially constrain time evolution, as shown in Paper I. In addition, the complicated processes of the ejection and acceleration of stellar winds are as yet unresolved and hence it is better to first address the simpler dynamics of twocomponent jets with the stellar outflow already being super Alfvénic. Besides, the launching of each component takes place at different and extended locations of the YSOdisk system and therefore the interaction happens at higher altitudes. Moreover, the low frequency models, 4aS10^{3}, 4bS10^{3} and 5S10^{5}, are obviously associated with larger length scales and therefore the vertical direction is chosen for the former and for the latter. Essentially, we address the length scales of a few tens AU radially to a thousand AU vertically.
All models have a uniform resolution of 256 zones for every 100 AU. However, we have evolved a typical model, 5q02, also in a finer grid of to investigate the properties common to all models, such as time evolution features, potential steady states, deviations from analytical solutions and shock formation. Nevertheless, the cell size was not found to affect the outcome of the numerical evolution, a feature of the selfsimilar models that is also supported by Paper I. Furthermore, the unperturbed simulations have been carried out up to a final time of 80 y, equivalent to , i.e. 80 000 Keplerian rotations of the protostar, or for up to t = 250 in code units. The time unit in the code corresponds to 0.32 y. Due to the greater length and time scales involved in models 4aS10^{3}, 4bS10^{3} and 5S10^{4} the simulations were run up to a final time of 160 y.
At the lower boundary, we keep all variables fixed to their initial values after the mixing, in agreement with the conclusions of Paper I (wherein a detailed discussion of the correct treatment of the boundary conditions can be found). Outflow or extrapolated boundary conditions in the region where the flow enters the computational domain might artificially influence the long term simulations. At the axis we apply axisymmetric boundary conditions, whereas at the upper and right boundaries, we apply outflow conditions. Note that setting the derivative of equal to zero at the right boundary could cause artificial collimation. However, the ADO solution dominates at this boundary both in the initial conditions and over time (as will be seen in the next section). Therefore, recalling from Paper I that the ADO model maintains its exact equilibrium in the rightmost regions despite the specification of outflow conditions, we argue that the configurations studied here are not prone to such a numerical enforcement.
4 Results
We outline first the results obtained that are common to all twocomponent jet simulations and then we discuss the effects of the mixing parameters and the time variability.
4.1 Time evolution and steady state
The logarithm of the density is plotted in Fig. 1 for different evolutionary stages of a typical twocomponent jet model (5q02). The initial conditions correspond to equilibrium in the regions where each analytical solution dominates. However, around the matching surface, the models are modified and hence a strong perturbation is generated during the first timesteps of the simulation. An MHD wave propagates through the ADO solution without leaving behind any significant rearrangements. On the contrary, the equilibrium of the ASO model is substantially restructured, with its density dropping roughly by an order of magnitude. In only Keplerian rotations of the footpoint of the matching surface, the stellar component has already been totally and self consistently modified in the presence of the ADO solution.
From the rightmost plot of Fig. 1, notice that the initial matching surface is still evident. Indeed, this is expected, due to the fixed boundary conditions at the lower boundary. At t = 50, the initial perturbation has almost left the domain, with the twocomponent jet having reached a steady state. Notice the formation of a weak steady shock, which can be seen almost along the diagonal direction of the computational domain.
Figure 2: The density fluctuations as a function of time calculated at the points (r_{0} ,z_{0}) = (5 ,50), (15 ,50) and (30 ,50). The first is located in the ASO dominated region, the second close to the matching surface upstream of the shock and the third in the ADO dominated region (model 5q02). 

Open with DEXTER 
In order to establish the conclusion that the two components can coexist in a steady state, we plot in Fig. 2 the density fluctuations for different time scales of three specific points. One is located inside the stellar component and the other two upstream of the shock, at the matching surface and in the disk wind, respectively. Evidently, for t > 2.5 the solution remains almost unchanged up to t = 250, i.e. a time longer by two orders of magnitude. The disk wind reaches the final exact equilibrium a bit later (at ), due to the slower wave velocities of this region. Therefore, the rightmost plot of Fig. 1 represents a very well preserved steady state of the two intrinsically different jet components.
Both the physical and geometrical properties of the two winds are by definition considerably diverse, since their selfsimilar symmetries are orthogonal to each other. In turn, the same holds true for their respective poloidal critical surfaces. Therefore, the steady state of such a twocomponent jet was not a straightforward expectation. Nevertheless, it is clearly shown in Figs. 1 and 2 that the two complementary winds manage to coexist. Also taking into account the artificial boundary effects present in long term simulations, investigated in Paper I, the above results are adequate to argue that the twocomponent jet models reach a well defined steady state.
Despite the fact that the plots concern a particular model, the same conclusions are valid for all the unperturbed scenarios presented in Table 2.
4.2 Deviations from the analytical solutions
Another crucial question that arises concerns how close the final outcome of the simulations is to the initial analytical solutions. In particular, the smaller the deviations are found to be, the more valid and robust are the analytical studies on the selfsimilar MHD outflows. This also implies the easy and appropriate extension of their conclusions to the twocomponent jet scenario, especially for the analytically derived disk winds.
Therefore, in Fig. 3 we plot the critical surfaces of the ADO solution (red crosses) and those of the final numerical twocomponent jet (thick blue lines), along with the logarithmic density contours and the magnetic fieldlines (red lines). The propagation of the perturbation described in Sect. 4.1 results in the slight modification of the fast magnetosonic and Alfvénic critical surface, as can be seen from their almost perfect match in Fig. 3. On the other hand, the slow magnetosonic critical surface seems to have collapsed towards the lower boundary. The slow waves generated initially at the matching surface by definition cannot pass to the sub slow domain. Consequently, if the surface stayed at its initial location, the waves would not have been able to leave the lower right region of the computational box. Such a phenomenon is observed in SC3 and SC5 runs of Stute et al. (2008), where matter accumulates downstream of the slow critical surface. However, in our case, the separatrix is being bent downward, tangentially to the lower boundary, hence allowing the initial perturbation to exit the simulated box. As a result, the lower right region shows a significantly higher degree of deviation from the initial conditions than the rest of the domain. Nonetheless, it does reach a steady state asymptotically. Note that the critical surface cannot be dragged away, due to the fixed boundary conditions, which describe a sub slow flow at z = 10. Similar features are also observed in Gracia et al. (2006) and in most runs of Stute et al. (2008).
The fact that the system finds an equilibrium so close to the analytical solution is due to the topological stability of the ADO model discussed in Paper I. No matter if the central part of the disk wind is substituted as a whole with a physically and geometrically different kind of outflow, the solution maintains all its properties, proving its stability.
Figure 3: Logarithmic density contours (thin blue lines) for model 5q02 at t = 50. The magnetic poloidal fieldlines are overplotted with red lines. In the lower right part, going clockwise, the fast magnetosonic, the Alfvénic and the slow magnetosonic critical surfaces are plotted with red crosses for the ADO solution and with thick blue lines for the final numerical twocomponent flow. 

Open with DEXTER 
The poloidal critical surfaces plotted in Fig. 3 also provide other insights into the twocomponent scenario. Close to the axis, they have an elliptical shape, as can be seen in the region very close to (0, 10), and eventually become conical, after the matching surface. Intuitively, this makes sense due to the difference in the symmetry of the accelerating mechanisms of the two winds. In Kwan et al. (2007), two types of outflows are observed, one emanating radially out of the protostar and the other being ejected at a constant angle with respect to the disk midplane. This implies a geometry of the poloidal critical surfaces similar to Fig. 3.
Figure 4: The physical variables plotted at z = 50 for the ADO model alone (solid line), the ASO model alone (dashed line), the initial setup of 5q02 (diamonds) and its final configuration (crosses). The quantities displayed from left to right are: in the top row , and , in the middle row V_{r}, V_{z} and , in the bottom row B_{r}, B_{z} and . 

Open with DEXTER 
Figure 5: The normalized integrals of motion plotted along three fieldlines of model 5q02, rooted at the positions (from bottom to top) (3, 10), (6 , 10), and (9 , 10), corresponding to the ASO, mixing and ADO dominated regions, respectively. The distance from the lower boundary is parametrized by s. 

Open with DEXTER 
In Fig. 4, the physical variables are plotted at the constant height z = 50, for the initial setup (diamonds) and final configuration (crosses) of model 5q02. In addition, the initial ADO (solid lines) and ASO (dashed lines) solutions are also shown before their combination. All plots present the effect of the mixing function. Close to the axis, the ASO model dominates, whereas approaching the right boundary, the ADO becomes the main contributor. A jump can be observed in most quantities, which represents the weak shock discussed in Sect. 4.3. Apart from the density and the poloidal component of the magnetic field, the initial and final configurations converge at large distances, showing the stability of the ADO solution. However, this happens far from the slow magnetosonic critical surface. The modifications the initial ASO solution undergoes can be seen from the final equilibrium reached close to the axis. Note that the temperature plot can be used as a guide when looking for twocomponent jet parameters appropriate to address observed jets.
Finally, the normalized integrals of motion (Eqs. (5)(9) of Paper I) are plotted in Fig. 5 along three selected fieldlines rooted at the points (3, 10), (6, 10), and (9, 10). One is in the ASO dominated part, one is in the ADO domain and the other is almost along the matching surface crossing the shock. In all cases the integrals are conserved with high accuracy, varying only within a few . At large distances from the shock, they tend to become constant, which indicates that the system reaches a steady state in all three regions. For the two inner fieldlines, at s = 5 and s = 100, respectively, the observed jumps are related to the crossing of the shock. In particular, the larger deviation from constancy occurs for the specific entropy integral Q, as expected.
4.3 Shock formation
Figure 6: Normalized discontinuities (from left to right) of the density, total pressure (thermal plus magnetic) and the entropy Q) across the shock, close to the point (5 , 30). Notice that s is increasing in the inverse direction of r. 

Open with DEXTER 
Figure 7: A family of the characteristics (thin solid lines) of the fast magnetosonic waves in a zoomed super fast magnetosonic region around the shock (thick solid line) for model 5q02. The thick dashed line is the initial matching surface. 

Open with DEXTER 
In Fig. 6, we plot the normalized density, total pressure (P + B^{2}/2) and entropy Q across the shock (direction right to left) around the position (5 , 30), very close to the one assumed in Paper I. This point is located inside the domain where the stellar outflow dominates. Apparently, the density seems to increase by a factor of 2, whereas the pressure increases by a factor of 4. On the contrary, the jump seen in the entropy is very weak, being an order of magnitude smaller: this is not surprising, recalling that , i.e. conditions very close to isothermal, wherein entropy remains unchanged across shocks. Also there is no heating/cooling present in any simulation of this paper, thus making the above analysis simpler.
In Fig. 7 we plot one of the two families of the characteristics of the fast magnetosonic waves (thin lines) for model 5q02, along with the initial matching surface (dashed line). It is evident that the shock (thick solid line) is not crossed by the downstream characteristics. This shows the causal disconnection of the two domains, upstream and downstream of the shock. In other words, the shock represents the horizon for the propagation of all MHD waves, coinciding with the numerical FMSS.
This feature is closely related to the ADO solution and was studied in detail in Paper I. However, the twocomponent case we present here is especially interesting for the following reason. The shock manifests even in the central area, where the contribution of the ASO model is total. This implies that it is not associated with the lower boundary, but on the contrary, it forms above it, intersecting the symmetry axis. Taking also into account the results of Paper I, the shock seems to be an intrinsic feature of the ADO solution. Consequently, the presence of the disk wind model in the twocomponent jet scenario has the remarkable characteristic of producing outflows that are causally disconnected to their launching region, despite the fact that the initial conditions causally connect the whole computational box.
4.4 Parameter study
In this section, we present the behavior of the twocomponent jets when we change the model parameters.
Figure 8: Logarithm of the density for the final numerical solutions of models 1q01, 3q05, 4q01, 6q05 (left to right). The two leftmost cases have a weaker stellar component compared to the two rightmost. 

Open with DEXTER 
Figure 8 shows the logarithm of the final density of the simulations carried out for models 1q01, 3q05, 4q01, 6q05 (left to right). When the position of the matching surface is rooted closer to the disk rather than the star, the shock seems to bend towards the midplane, confining the unmodified ADO solution in a smaller domain. This result indicates that as the spatial domination of the ASO solution becomes larger, the ADO model controls a smaller portion of the box, thus forming the shock closer to the disk.
Figure 9: Logarithmic poloidal velocity and streamlines (dashed lines) for the unmixed ADO (top left) and ASO (top right) models separately. In the lower panel, models 7B05 ( ; left) and 8B5 ( ;right) are shown. The maximum values of the poloidal velocity of both of these twocomponent cases are km s^{1}, despite the misleading colorbar, which was accordingly chosen to match that of the top panel. 

Open with DEXTER 
Recalling that models 1q01 and 3q05 have a weaker ASO contribution ( ), compared to 4q01 and 6q05 ( ), Fig. 8 also suggests that the relative strength of the ASO model does not seem to considerably affect the position of the shock, although larger deviations are seen around the slow critical surface of the initial ADO model. The same result is derived from models 7B05, 8B5 (bottom of Fig. 9) and 9B10, where the location of the shock is found farther from the axis, the larger the value of . Nevertheless, this might be related to the previous result since a large also spatially reduces the contribution coming from the ADO model.
Furthermore, Fig. 9 presents the logarithm of the poloidal velocity and the streamlines (dashed lines) for the ADO and ASO solutions separately, as well as for models 7B05 ( ) and 8B5 ( ). The left plot of the bottom panel suggests that for disk wind dominated jets, the ADO solution is effectively collimating the central component. However, we know that polytropically evolved ASO solutions become more collimated and less dense than the non polytropic initial ASO models (Paper I). So, it is rather difficult to disentangle the collimation due to the disk wind and that due to the change in energetics.
Moreover, increasing by one order of magnitude the contribution of the ASO model, the streamlines take an almost radial geometry (lower panel, right plot of Fig. 9). A similar result was obtained by Meliani et al. (2006) when the mass loss rate of the inner stellar wind becomes comparable to the disk mass loss rate. Although this might contradict the parallel flow structure seen in the right plot of the top panel of Fig. 9, where the ASO solution is plotted alone, we note that such a strong collimation comes from the linear increase of (Fig. 4, dashed line). However, the twocomponent jet presents a more realistic distribution of current, with a decreasing toroidal field at large distances (Fig. 4, crosses) and hence the hoop stress is not capable of collimating the flow.
Figure 10: Logarithmic pressure contours for model 10d1 on the left and 11d4 on the right. 

Open with DEXTER 
Finally, we examine how the third free parameter, d, which defines the steepness of the transition region, influences the final steady state reached by the twocomponent jets.
The pressure contours of models 10d1 and 11d4 shown in Fig. 10 suggest that no matter how smoothly the variables change from one solution to the other, the matching surface reaches the same sort of structure at the end of the simulations. On the other hand, the shock is affected more dramatically. In the d = 1 case, it has a shape very similar to the one forming without the presence of the ASO solution (see Paper I), with a polar angle of as calculated close to the origin. On the contrary, the shock intersects the axis with a wider polar angle in the d = 4 case. Note that although the value of the parameter may change inside the computational box, it is kept fixed at the lower boundary and hence influences the evolution.
4.5 Time variable stellar or Xtype winds
This last section is dedicated to the stability issues raised by a potential time variability in the YSO's outflow. We apply time dependency (Eq. (16) or (17)) either at the stellar wind's base or around the Xpoint located at the interface between the stellar magnetosphere and the disk. The twocomponent jet parameters adopted are identical to model 5q02.
Figure 11: Magnetic fieldlines (red) and logarithmic density contours at t = 250 for models 1S1 (left) and 3S100 (right). 

Open with DEXTER 
High frequency velocity (or density) variations, associated with the Keplerian rotation at roughly a stellar radius, seem to fade away on larger scales, as shown on the left of Fig. 11. The structure remains very close to the unperturbed model. Two orders of magnitude lower frequency fluctuations, result in stronger gradients along the flow, as seen in the right panel of Fig. 11. Considering that the velocity varies by of its initial value, it is surprising how well the twocomponent jet structure is retained. Despite the ``wiggling'' of the magnetic field, the same flow features are found as in the unperturbed cases.
Figure 12: Logarithmic density contours (thin blue lines) and the magnetic field (red lines) for model 3S10^{2} at t = 50. The red crosses and the thick blue lines denote the critical poloidal velocity surfaces of the ADO solution and those of the final equilibrium, respectively. This plot is equivalent to Fig. 3. 

Open with DEXTER 
Figure 12 displays a plot equivalent to Fig. 3 in order to understand how the shock and critical surfaces change in the timevariable stellar wind case. The picture is very similar, apart from the perturbations seen in the density throughout the computational box. The poloidal critical surfaces show the same behavior as in the unperturbed models and the weak steady shock is still present, being slightly curved locally as the fluctuations propagate.
Figure 13: Logarithm of the density along with magnetic fieldlines (red lines) at t = 50 for models 6X1 and 8X10^{2}. 

Open with DEXTER 
Analogous results are derived by the models where time variability is enforced at the density and velocity of the outflow around the Xpoint. The momentum changes periodically by an order of magnitude. However, at the ASO and ADO dominated regions, the wind characteristics do not seem to be affected, especially in the high frequency variability case (model 6X1, plot of left panel of Fig. 13). On the other hand, more evident structures are produced in the 100 times lower frequency fluctuations, still without destroying the basic pattern (model 8X10^{2}, plot of right panel Fig. 13). This behaviour is similar to the stellar wind variability, but with a lesser degree of collimation.
In order to see how the low frequency variability can affect the jet far away from the launching region, we plot the quantity for models 4aS10^{3} (top left), 4bS10^{3} (bottom left) and 5S10^{4} (right) in Fig. 14. Close to the base, the numerical solutions remain close to the initial ADO and ASO models. However, at higher altitudes, the fluctuations create knotlike structures. It is evident that such models can be associated with some jet variability. We have checked that both stellar and Xwind type pulsations produce very similar structures far away from the diskstar system. This was expected, since kAU scales cannot distinguish the ejections coming from within 1 AU. The regular knot spacing observed in the jet of HH30 ( AU, Bacciotti et al. 1999) can be reasonably compared with our models 4aS10^{3} and 4bS10^{3}, with a structure periodicity of yr. Model 5S10^{4}, with a periodicity of 10 yrs, could be associated with the knots detected in the jet of HH34 where the condensation spacing is AU (Cohen & Jones 1987). Nevertheless, in this case there is a gap between the blobs and the star, suggesting a contribution of other processes to the knot formation. Note that the time scales of such fluctuations also correspond to typical stellar variabilities (e.g. the 11 y period of the solar cycle).
Figure 14: The quantity , which is roughly related with emissivity, is plotted for the low frequency variability models examined in large scales. In the left, models 4aS10^{3} (top) and 4bS10^{3} (bottom) are shown, whereas in the right model 5S10^{4} is displayed. Note that for model 4aS10^{3}, for model 4bS10^{3} and for model 5S10^{4}. However, the colorbars in the online version use a lower maximum value in order to enhance the displayed features. The length code unit corresponds to 1 AU. 

Open with DEXTER 
Figure 15: Density, pressure and sonic Mach number plotted along the jet axis for model 4aS10^{3}. 

Open with DEXTER 
Finally, Fig. 15 provides the proof that these knotlike structures are in fact shocks. The top panel displays the periodic density and pressure jumps along the jet axis, with the change being approximately an order of magnitude for both. Note that close to the base the discontinuities are not yet well developed. We also remark that these shocks are stronger than that associated with the FMSS (see Fig. 6). The lower panel reports the sonic Mach number as a function of z. Its mean value of the background flow is , in good accordance with YSO jet observations. The shocks propagate faster by , as expected in agreement with the inflow time variability. Although Fig. 15 suggests that the flow values converge to a similar periodic structure, a larger computational box is needed to verify such an argument. In a future study, we plan to apply radiation cooling effects to these time variable models, effectively producing realistic emission maps to be compared to real data.
5 Summary and conclusions
In this work, we have constructed twocomponent numerical jet models by properly combining two well studied analytical solutions, each one describing separately a disk wind (ADO) and a stellar outflow (ASO). We have investigated the features of the time evolution and the characteristics of the final outcome of the simulations as a function of the twocomponent jet parameters and the enforced time variability. Although the detailed launching mechanisms of each component are not taken into account, the twocomponent jet models presented here seem able to capture the dynamics and describe a variety of interesting scenarios.
The main conclusions of this work are the following:
 
 The twocomponent jet models show remarkable stability and always reach a well defined steady state. This result is robust despite the fact that the two solutions have orthogonal symmetries, different geometry and different physics (i.e. launching mechanisms). In addition, the conclusion holds true independently of the choice of the parameters and even in the cases where time variability is enforced at the stellar wind's base or around the Xpoint. Therefore, the analytical solutions provide solid foundations for realistic twocomponent jet scenarios.
 
 The system remains close to the initial analytical solutions. In particular, the disk wind dominated regions are barely changed in the presence of the stellar outflow, with the exception of the slow magnetosonic regions. On the other hand, the central component is self consistently modified due to the assumption of a polytropic equation of state and because of the effective collimation caused by the surrounding disk wind. This implies that specific YSO systems can be addressed more accurately by constructing analytical outflow solutions with the desirable characteristics, before merging them into a twocomponent regime.
 
 A shock manifests during the time evolution, preventing any information from the downstream domain from reaching the base of the outflow. This separatrix causally disconnects the twocomponent jet from its launching regions, although initially there is no such ``horizon'' present in the computational box. The initial ASO solution does not exhibit any modified fast separatrix (Sauty et al. 2002), whereas despite the existence of the FMSS in the initial ADO model (at small polar angles), it is effectively replaced by the stellar wind in the initial setup. Nonetheless, the final equilibria reached by the numerical models show the formation of a weak shock corresponding to such a surface, causally disconnecting the acceleration regions from the jet propagation physics and subsequent interaction with the outer medium.
 
 We may address various twocomponent jet scenarios, by means of two parameters controlling the relative contribution of each component, , and the time variability function, f(t). With the former, we can smoothly switch the physics from a totally magnetocentrifugal wind to a pressure driven jet. With the latter, flow fluctuations are introduced, producing knotlike structures on large scales that are quantitatively similar to HH30 and HH34 observations.
Acknowledgements
We acknowledge V. Cayatte and the rest of the group in LUTh for fruitful discussions, and an anonymous referee for helpful comments and suggestions that resulted in a better presentation of this work. We would also like to thank Capt. D. Kalogeras whose support during the revision of this paper prevented a delay of several months. The present work was supported in part by the European Community's Marie Curie Actions  Human Resource and Mobility within the JETSET (Jet Simulations, Experiments and Theory) network under contract MRTNCT2004 005592 and in part by the HPCEUROPA++ project (project number: 211437), with the support of the European Community  Research Infrastructure Action of the FP7 ``Coordination and support action'' Programme.
Appendix A: The selfsimilar outflow formulation
Axisymmetry, steady state and selfsimilarity assumptions simplify the ideal MHD equations to a set of coupled ODEs in spherical coordinates. These equations are solved numerically, providing the values of some key functions for each model.
For the ADO solution (radially selfsimilar), the physical variables are provided in terms of the tabulated key functions , and :
where p denotes the poloidal component.
The ASO solution (meridionally selfsimilar) is described with the help of the key functions , , and :
References
 Alencar, S. H. P., & Batalha, C. 2002, ApJ, 571, 378 [NASA ADS] [CrossRef] (In the text)
 Anderson, J. M., Li, Z.Y., Krasnopolsky, R., & Blandford, R. D. 2003, ApJ, 590, L107 [NASA ADS] [CrossRef] (In the text)
 Bacciotti, F., Eislöffel, J., & Ray, T. P. 1999, A&A, 350, 917 [NASA ADS] (In the text)
 Bacciotti, F., Ray, T. P., Mundt, R., Eislöffel, J., & Solf, J. 2002, ApJ, 576, 222 [NASA ADS] [CrossRef] (In the text)
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] (In the text)
 Bogovalov, S., & Tsinganos, K. 2001, MNRAS, 325, 249 [NASA ADS] [CrossRef] (In the text)
 Cabrit, S., Edwards, S., Strom, S. E., & Strom, K. M. 1990, ApJ, 354, 687 [NASA ADS] [CrossRef] (In the text)
 Casse, F., & Keppens, R. 2002, ApJ, 581, 988 [NASA ADS] [CrossRef] (In the text)
 Casse, F., & Keppens, R. 2004, ApJ, 601, 90 [NASA ADS] [CrossRef] (In the text)
 Coffey, D., Bacciotti, F., Woitas, J., Ray, T. P., & Eislöffel, J. 2004, ApJ, 604, 758 [NASA ADS] [CrossRef] (In the text)
 Cohen, M., & Jones, B. F. 1987, ApJ, 321, 846 [NASA ADS] [CrossRef] (In the text)
 Dougados, C., Cabrit, S., Lavalley, C., & Ménard, F. 2000, A&A, 357, L61 [NASA ADS] (In the text)
 Edwards, S., Fischer, W., Hillenbrand, L., & Kwan, J. 2006, ApJ, 646, 319 [NASA ADS] [CrossRef] (In the text)
 Fendt, C. 2006, ApJ, 651, 272 [NASA ADS] [CrossRef] (In the text)
 Fendt, C. 2009, ApJ, 692, 346 [NASA ADS] [CrossRef] (In the text)
 Ferreira, J. 1997, A&A, 319, 340 [NASA ADS] (In the text)
 Ferreira, J., Pelletier, G., & Appl, S. 2000, MNRAS, 312, 387 [NASA ADS] [CrossRef] (In the text)
 Ferreira, J., Dougados, C., & Cabrit, S. 2006, A&A, 453, 785 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gracia, J., Vlahakis, N., & Tsinganos, K. 2006, MNRAS, 367, 201 [NASA ADS] [CrossRef] (In the text)
 Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736 [NASA ADS] [CrossRef] (In the text)
 Hartigan, P., Edwards, S., & Pierson, R. 2004, ApJ, 609, 261 [NASA ADS] [CrossRef] (In the text)
 Johns, C. M., & Basri, G. 1995, ApJ, 449, 341 [NASA ADS] [CrossRef] (In the text)
 Kwan, J., Edwards, S., & Fischer, W. 2007, ApJ, 657, 897 [NASA ADS] [CrossRef] (In the text)
 Matt, S., & Pudritz, R. 2005, ApJ, 632, L135 [NASA ADS] [CrossRef] (In the text)
 Matt, S., & Pudritz, R. 2008a, ApJ, 678, 1109 [NASA ADS] [CrossRef] (In the text)
 Matt, S., & Pudritz, R. 2008b, ApJ, 681, 391 [NASA ADS] [CrossRef]
 Matt, S., Goodson, A. P., Winglee, R. M., & Bohm, K.H. 2002, ApJ, 574, 232 [NASA ADS] [CrossRef] (In the text)
 Matsakos, T., Tsinganos, K., Vlahakis, N., et al. 2008, A&A, 477, 521 [NASA ADS] [CrossRef] [EDP Sciences] (Paper I) (In the text)
 Meliani, Z., & Keppens, R. 2007, A&A, 475, 785 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Meliani, Z., Casse, F., & Sauty, C. 2006, A&A, 460, 1 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] (In the text)
 Pudritz, R., Rogers, C., & Ouyed, R. 2006, MNRAS, 365, 1131 [NASA ADS] [CrossRef] (In the text)
 Sauty, C., & Tsinganos, K. 1994, A&A, 287, 893 [NASA ADS] (In the text)
 Sauty, C., Trussoni, E., & Tsinganos, K. 2002, A&A, 389, 1068 [NASA ADS] [CrossRef] [EDP Sciences] (STT02) (In the text)
 Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781 [NASA ADS] [CrossRef] (In the text)
 Stempels, H. C., & Piskunov, N. 2002, A&A, 391, 595 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Stute, M., Tsinganos, K., Vlahakis, N., Matsakos, T., & Gracia, J. 2008, A&A, 491, 339 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tesileanu, O., Mignone, A., & Massaglia, S. 2008, A&A, 488, 429 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Trussoni, E., Tsinganos, K., & Sauty, C. 1997, A&A, 325, 1099 [NASA ADS] (In the text)
 Tsinganos, K. C. 1982, ApJ, 252, 775 [NASA ADS] [CrossRef] (In the text)
 Tsinganos, K., Sauty, C., Surlantzis, G., Trussoni, E., & Contopoulos, J. 1996, MNRAS, 283, 811 [NASA ADS] (In the text)
 Tzeferacos, P., Ferrari, A., Mignone, A. et al. 2009, MNRAS, in press (In the text)
 Vlahakis, N., & Tsinganos, K. 1998, MNRAS, 298, 777 [NASA ADS] [CrossRef] (In the text)
 Vlahakis, N., Tsinganos, K., Sauty, C., & Trussoni, E. 2000, MNRAS, 318, 417 [NASA ADS] [CrossRef] (VTST00) (In the text)
 Zanni, C., Ferrari, A., Rosner, R., Bodo, G., & Massaglia, S. 2007, A&A, 469, 811 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
Footnotes
 ... ratio^{}
 In Paper I, VTST00, STT02 and all previous studies on selfsimilar outflows, this integral was defined as ``specific angular momentum''.
 ...
^{}  The value of the parameter x = 0.75, is related to the ejection index of Ferreira (1997) and corresponds to zero ejection according to its expression. However, from Figs. 5 and 6 of VTST00 it is evident that the solution with x = 0.7575, i.e. , is almost identical to the one with x = 0.75 for . Therefore, we argue that the ADO solution employed here should not contradict the theoretical arguments presented in Ferreira (1997).
 ...^{}
 and are related to the rotational velocity, and to the longitudinal profile of the pressure and density, respectively, to the energetic balance across the poloidal fieldlines and to the energy input. The constants and measure the gravitational potential for each solution, whereas is associated to the relative magnitudes of magnetic and thermal pressure.
 ... function^{}
 Note that should not be confused with A. Although the former is a simple sum of and , the latter is computed from the mixing function, as the rest of the variables.
 ... PLUTO^{}
 Publicly available at http://plutocode.to.astro.it
All Tables
Table 1: Parameters characterizing the adopted analytical solutions.
Table 2: A short description and the parameters of the unperturbed numerical models.
Table 3: The time variable numerical models.
All Figures
Figure 1: Logarithms of the density at different times for a typical model (5q02). The time unit is 0.32 y (or 10 Keplerian rotations of the footpoint of the matching surface). In the rightmost plot, a weak shock is observed along the diagonal. In addition, the initial matching surface, approximately at r = 5 on the lower boundary, is still evident. To check this, compare it with the leftmost plot which describes a configuration very close to the initial one. 

Open with DEXTER  
In the text 
Figure 2: The density fluctuations as a function of time calculated at the points (r_{0} ,z_{0}) = (5 ,50), (15 ,50) and (30 ,50). The first is located in the ASO dominated region, the second close to the matching surface upstream of the shock and the third in the ADO dominated region (model 5q02). 

Open with DEXTER  
In the text 
Figure 3: Logarithmic density contours (thin blue lines) for model 5q02 at t = 50. The magnetic poloidal fieldlines are overplotted with red lines. In the lower right part, going clockwise, the fast magnetosonic, the Alfvénic and the slow magnetosonic critical surfaces are plotted with red crosses for the ADO solution and with thick blue lines for the final numerical twocomponent flow. 

Open with DEXTER  
In the text 
Figure 4: The physical variables plotted at z = 50 for the ADO model alone (solid line), the ASO model alone (dashed line), the initial setup of 5q02 (diamonds) and its final configuration (crosses). The quantities displayed from left to right are: in the top row , and , in the middle row V_{r}, V_{z} and , in the bottom row B_{r}, B_{z} and . 

Open with DEXTER  
In the text 
Figure 5: The normalized integrals of motion plotted along three fieldlines of model 5q02, rooted at the positions (from bottom to top) (3, 10), (6 , 10), and (9 , 10), corresponding to the ASO, mixing and ADO dominated regions, respectively. The distance from the lower boundary is parametrized by s. 

Open with DEXTER  
In the text 
Figure 6: Normalized discontinuities (from left to right) of the density, total pressure (thermal plus magnetic) and the entropy Q) across the shock, close to the point (5 , 30). Notice that s is increasing in the inverse direction of r. 

Open with DEXTER  
In the text 
Figure 7: A family of the characteristics (thin solid lines) of the fast magnetosonic waves in a zoomed super fast magnetosonic region around the shock (thick solid line) for model 5q02. The thick dashed line is the initial matching surface. 

Open with DEXTER  
In the text 
Figure 8: Logarithm of the density for the final numerical solutions of models 1q01, 3q05, 4q01, 6q05 (left to right). The two leftmost cases have a weaker stellar component compared to the two rightmost. 

Open with DEXTER  
In the text 
Figure 9: Logarithmic poloidal velocity and streamlines (dashed lines) for the unmixed ADO (top left) and ASO (top right) models separately. In the lower panel, models 7B05 ( ; left) and 8B5 ( ;right) are shown. The maximum values of the poloidal velocity of both of these twocomponent cases are km s^{1}, despite the misleading colorbar, which was accordingly chosen to match that of the top panel. 

Open with DEXTER  
In the text 
Figure 10: Logarithmic pressure contours for model 10d1 on the left and 11d4 on the right. 

Open with DEXTER  
In the text 
Figure 11: Magnetic fieldlines (red) and logarithmic density contours at t = 250 for models 1S1 (left) and 3S100 (right). 

Open with DEXTER  
In the text 
Figure 12: Logarithmic density contours (thin blue lines) and the magnetic field (red lines) for model 3S10^{2} at t = 50. The red crosses and the thick blue lines denote the critical poloidal velocity surfaces of the ADO solution and those of the final equilibrium, respectively. This plot is equivalent to Fig. 3. 

Open with DEXTER  
In the text 
Figure 13: Logarithm of the density along with magnetic fieldlines (red lines) at t = 50 for models 6X1 and 8X10^{2}. 

Open with DEXTER  
In the text 
Figure 14: The quantity , which is roughly related with emissivity, is plotted for the low frequency variability models examined in large scales. In the left, models 4aS10^{3} (top) and 4bS10^{3} (bottom) are shown, whereas in the right model 5S10^{4} is displayed. Note that for model 4aS10^{3}, for model 4bS10^{3} and for model 5S10^{4}. However, the colorbars in the online version use a lower maximum value in order to enhance the displayed features. The length code unit corresponds to 1 AU. 

Open with DEXTER  
In the text 
Figure 15: Density, pressure and sonic Mach number plotted along the jet axis for model 4aS10^{3}. 

Open with DEXTER  
In the text 
Copyright ESO 2009