A Two-Temperature Model of Magnetized Protostellar Outflows Liang-Yao Wang (王亮堯)1,2 , Hsien Shang (尚賢)2 , Ruben Krasnopolsky2 , Tzu-Yang Chiang (江 子揚)2 arXiv:1512.02609v1 [astro-ph.SR] 8 Dec 2015 October 14, 2018 ABSTRACT We explore kinematics and morphologies of molecular outflows driven by young protostars using magnetohydrodynamic simulations in the context of the unified wind model of Shang et al. The model explains the observed high-velocity jet and lowvelocity shell features. In this work we investigate how these characteristics are affected by the underlying temperature and magnetic field strength. We study the problem of a warm wind running into a cold ambient toroid by using a tracer field that keeps track of the wind material. While an isothermal equation of state is adopted, the effective temperature is determined locally based on the wind mass fraction. In the unified wind model, the density of the wind is cylindrically stratified and highly concentrated toward the outflow axis. Our simulations show that for a sufficiently magnetized wind, the jet identity can be well maintained even at high temperatures. However, for a high temperature wind with low magnetization, the thermal pressure of the wind gas can drive material away from the axis, making the jet less collimated as it propagates. We also study the role of the poloidal magnetic field of the toroid. It is shown that the wind-ambient interface becomes more resistant to corrugation when the poloidal field is present, and the poloidal field that bunches up within the toroid prevents the swept-up material from being compressed into a thin layer. This suggests that the ambient poloidal field may play a role in producing a smoother and thicker swept-up shell structure in the molecular outflow. Subject headings: ISM: jets and outflows – stars: winds, outflows 1. INTRODUCTION Bipolar molecular outflows from young stellar objects have been found to accompany the star formation process since the early, embedded phase. Understanding their properties may ultimately 1 Graduate Institute of Astronomy and Astrophysics, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan; lywang@asiaa.sinica.edu.tw 2 Academia Sinica, Institute of Astrophysics (ASIAA), and Theoretical Institute for Advanced Research in Astrophysics (TIARA), Academia Sinica, 11F of Astronomy-Mathematics Building, AS/NTU. No.1, Sec. 4, Roosevelt Rd, Taipei 10617, Taiwan –2– help shed light on the central compact protostar-disk systems that is believed to drive the energetic phenomenon. Historically, wind-blown and jet-driven models have been raised to account for the different observed morphologies of a wider opening shell or a collimated jet-like structure (Lee et al. 2000). Later, it is found in a young group of (Class 0) outflows that both features are simultaneous present (e.g., IRAS 04166+2706, Santiago-Garcı́a et al. 2009). It is therefore important for models aiming to provide a comprehensive picture of outflow phenomenon to be capable of explaining such dual characteristics. One model that provides a natural reconciliation between the jet and wind nature is the unified wind model proposed by Shang et al. (2006) (hereafter paper I). It consists of a wide angle wind which has a stratified density structure resembling that predicted by the X-wind theory (Shu et al. 2000). The wind is assumed to first run into an ambient toroid structure expected to form from a dense core that is partially supported by an ordered magnetic field which it inherited from the parent cloud. The dual features can be naturally reproduced in this scenario where the axially concentrated part of the wind forms a collimated jet while the more tenuous diverging part of the wind impacts upon the ambient and forms a low-velocity shell. The model has been successfully demonstrated using magnetohydrodynamic (MHD) simulations in paper I. Outflow properties can vary from source to source on top of the common characteristics of bipolarity and jet/shell features. The diversity in cavity shape, opening angle, kinematic pattern, knot structure, etc. can in principal reflect the underlying physical conditions as well as the evolution history. An important key to extract information from the rich phenomenon is to build the connection between physical factors and their consequences on a larger observable scale. Numerical experiments of outflow models are particularly valuable in this respect. In paper I it has been demonstrated that the unified wind model can provide a natural explanation of the jet and shell features, but the roles of different physical parameters remain to be explored. The magnetic field is core to the wind model and its role is naturally of great interests. On the other hand, molecular line studies have suggested a possible warm gas content in the molecular jet (Nisini et al. 2007), which also motivates an extension of the cold simulations in paper I to a higher temperature regime. Therefore in this work we examine the consequences of varying the magnetic field strength and introducing a different gas temperature in context of the unified wind model of paper I. Meanwhile, the poloidal field of the toroid may also play a role in the wind-ambient interaction and we also explore its effects. We briefly introduce the theoretical background and the numerical setup in section 2 and 3, and present our results in section 4. We discuss the implication in section 5 and conclude in section 6. –3– 2. 2.1. MODEL BACKGROUND Asymptotic Structure of Magnetocentrifugal Winds Magnetocentrifugal wind is probably the most promising mechanism to date to explain the launching of bipolar outflows from the young stellar objects. As in paper I, we adopt the same limiting case of the asymptotic behavior predicted by the X-wind theory. The density, toroidal magnetic field, and velocity of the wind take the following forms: D0 B0 , Bφ (̟) = , vw = V0 , (1) 2 ̟ ̟ where ̟ = r sin θ is the cylindrical radius from the axis, and the wind is in radial direction. The density constant D0 is related to and is determined by the adopted (two-sided) wind mass-loss rate ṀW through ṀW = 2πD0 V0 ln |(1 + cos θ1 )/(1 − cos θ1 )|. The limiting angle θ1 is taken to be half the first cell size of the θ-grid. The magnetic field strength constant B0 is related to the √ density constant through B0 = 4πD0 V0 /MA (in Gaussian convention), where MA is the Alfvén Mach number of the wind. Thus for a given mass-loss rate (a given density constant), the magnetic field strength is inversely proportional to MA . We adopt a constant velocity of 50 km s−1 for the wind, and with the chosen velocity value we further fix the wind density constant D0 for a wind mass-loss rate of 3×10−6 M⊙ yr−1 . The initial toroidal magnetic field strengths is ∼ 0.75/MA gauss at ̟ = 1 AU from the outflow axis given the adopted values. ρ(̟) = 2.2. Singular Isothermal Toroids For the young outflows from the embedded system, the wind will first interact with the protostellar envelopes, which are part of the star forming core yet to be accreted onto the central protostar-disk. When dynamically important magnetic fields are present during core formation, anisotropic mass distribution is expected to form due to the easier settling of mass along the field lines than across. This problem was studied by Li & Shu (1996) and an equilibrium self-similar model was proposed to describe the core density ρ(r, θ) and magnetic flux function Φ(r, θ). The radial and angular dependence of the functions are assumed to take a separable form of a2 4πa2 r R(θ), Φ = φ(θ), (2) 2πGr2 G1/2 where R(θ) and φ(θ) are dimensionless angular functions. The poloidal field is obtained through 1 Φ B(r, θ) = 2π ∇ × ( r sin θ eφ ). The possible solutions are characterized by a single parameter H0 representing the fractional overdensity supported by the magnetic field above that by the thermal pressure. An equivalent parameter n (n = 4H0 ) is sometimes used in place of H0 , and a larger n corresponds to a more flattened core structure. We adopt the density distribution of n = 2 toroid as the initial condition. Simulations are carried out both with and without including the poloidal magnetic field threading the matter. Such ambient field was ignored in paper I for simplicity, and we will examine its role here. ρ(r, θ) = –4– 3. SIMULATION SETUP The problem of a warm magnetized wind running into an ambient toroid is studied using magnetohydrodynamic (MHD) simulations. We evolve the 2D axisymmetric problem in spherical coordinate (r, θ, φ) with the outflow axis lying at θ = 0◦ . The usual set of ideal MHD equations are solved using the ZEUS-TW code (Krasnopolsky et al. 2010). A scaled grid of 1600 points is used to cover the computation domain from r = 1.5 AU to 2 × 105 AU. The finest grid spacing is 0.08 AU at the innermost cell. We only consider one side of the bipolar outflow and the θ grid is in uniform step of 0.125◦ across θ = 0◦ to 90◦ . The initial condition is set up with a toroid with zero velocity. The density distribution is always initialized with the toroid solution while the poloidal field threading it is ignored in the first part of our exploration. We later include this poloidal field component to gain insight to its role. The boundary conditions on the outflow axis and at the equator are set by symmetry. Outflow conditions are set for the outer r boundary. A wind is imposed at the inner r boundary by setting the density and toroidal magnetic field values with the wind solutions in every time step. The inner r boundary conditions of the V × B are set by continuity to avoid anchoring the poloidal filed lines. There are several differences in terms of numerical setup between this work and paper I. For example, we have moved from the ZEUS-2D code (Stone & Norman 1992) to ZEUS-TW in order to take advantage of the much faster production speed of parallel computing. This allows us to finish high resolution tests within affordable time. We have also adopted a spherical coordinate system instead of a cylindrical one to accommodate the radial wind velocity more easily. The wind zone, a box that was used to adapt the spherical wind into a cylindrical grid, is abandoned, and the winds are imposed through boundary conditions (the ghost zones). We checked that the new numerical implementation produces nearly identical output for the same problem setup. Hence the results in this work should be consistent with that of paper I. For our study of a warm gas propagating into a cold medium, two characteristic temperatures are specified. For the ambient toroid we adopt a value of aambient = 2 × 104 cm s−1 (∼ 10 K), which is also the global isothermal sound speed used in paper I. We adapt the code by introducing a second sound speed awind to evolve the wind while retaining an isothermal equation of state. This is done by allowing an effective temperature to be determined locally based on the presence of wind or ambient material at each position. An evolving tracer field (ρwind ) is used to keep track of the wind material and facilitate this prescription. The gas is cold if it comes from the ambient toroid and warm if injected as a wind. With these considerations, we come to adopt the formula of q a = f a2wind + (1 − f ) a2ambient for the effective sound speed. The parameter f is the wind-originated mass fraction defined as f ≡ ρwind /ρ. Such formula allows the pure toroid and wind material to evolve under their own isothermal sound speed, while a simple linear interpolation between the two temperatures is used for region of mixing. Note the simulations essentially reduce to a globally isothermal condition when awind is set equal to aambient . As mentioned, to enable calculation of this two-temperature –5– isothermal scheme, we keep track of the wind material by using a tracer field. We initialize the tracer to zero at t = 0 in the simulations when no wind material is yet injected. The same wind boundary condition used to set the density at the inner radius is also applied to the tracer field in every time step. Since it simply follows the matter transport passively, it keeps a copy of all wind-originated material. 4. NUMERICAL RESULTS In the following subsections we present our simulation results. First we demonstrate the general features of the warm magnetized wind simulation with one example case. We then explore the consequences of varying the wind temperature and toroidal magnetic field strength of the wind by examining a group of runs with different awind and MA . Probing the relative importance of thermal and magnetic pressures in each cases using for instance the plasma beta parameter (βplasma ≡ 8πa2 ρ/B 2 ) may help provide insights. Note in the first part of our explorations (sections 4.1 and 4.2) the poloidal magnetic field threading the ambient medium is not included. In section 4.3 we initialize the poloidal magnetic field with the proper toroid solution of Li & Shu (1996). 4.1. Basic Features As a demonstration of the general results, we present snapshots of a typical simulation at t = 100 and 1000 yr in the upper and lower row panels of Figure 1, respectively. For this particular case the wind sound speed is awind = 1.2 × 105 cm s−1 (∼ 400 K) and the Alfvén Mach number is MA = 30. From left to right columns are maps of the density, magnitude of poloidal velocity (in km s−1 ), the tracer field ρwind , and the wind fraction f . The density and wind density maps in this work are plotted as number density of molecules assuming a mean molecular weight of 3.88 × 10−24 g. The basic characteristic of jet and shell features can be clearly identified in the density map. The collimated structure concentrated on the z-axis is a result of the intrinsic wind density profile (equation 1). Although the initial velocity is in radial direction, one does not see a significant density drop (near the axis) expected for a spherically expanding wind. This is due to magnetic collimation of the toroidal field in the wind which helps maintain this cylindrically stratified density structure. Moving away from the central axis, the flow density decreases with cylindrical radius until reaching a higher density structure at the wind-ambient interface, which is the swept-up shell. If we further compare the upper and lower panels, one can see the clear resemblance between the swept-up shell at different epochs of t = 100 and 1000 yr. As pointed out in paper I, this self-similar behavior is not surprising since both the asymptotic wind and the toroid solutions are scale-free. The wind propagates mostly near its initial speed of 50 km s−1 until dropping abruptly at the windtoroid interface, as shown by the velocity map. This confirms that the wind-ambient interaction –6– still lies within the supersonic/superalfvénic regime in our new simulation as in the cold flow of paper I. This example case basically demonstrates that the overall dynamic of the unified wind model is qualitatively unchanged when we increase the wind temperature from ∼10 K (in paper I) to ∼400 K. In the third column of Figure 1 are snapshots of the tracer field. This field is a critical component in our simulations which facilitates the calculation of effective temperatures and sound speeds in our two-temperature prescription. Moreover, it keeps a record of the distribution of wind-originated material which can not be seen in the density maps. The figures show that the primary wind is confined within a limited volume. The values of the tracer is essentially zero or negligible outside this region, but appear to be very similar or identical to the density maps (in the first column) within the region. Note while the central jet component is clearly seen in the tracer field maps, the shell feature is not. This implies different compositions and origins of the two characteristic features. While the jet structure is essentially part of the primary wind, the shell consists mostly of the swept-up ambient material which does not show up in the tracer field maps. The boundary where the density of wind material drops to near zero is close to the inner edge of the swept-up shell. Another way to visualize the wind distribution is to examine the relative mass composition at each position. We show the results in the fourth column of Figure 1. The wind fraction of 10% and 90% are highlighted with white and gray contours. It is clear from the resulting maps that the majority part of the simulation are either wind-dominated (within the gray contour, f > 90%) or ambient-dominated (outside the white contour, f < 10%). The transition from wind to toroid usually happens rather quickly with only little intermediate region. Hence we may understand the two-temperature simulations as follow: the major part of the wind evolves isothermally under the higher temperature specified for the wind, while the toroid retains its original temperature of ∼ 10 K. Materials in between the two contours are subjected to notable degree of mixing, and are the region that will take on the interpolated temperature value in the simulation. While the division between wind and ambient appears obvious, there are structures with notable degree of mixing located within the wind-dominated region. For example, at a height of z ≃ 350, 450, and 700 AU in the upper panel of wind fraction map, there are blobs of gas with f ≃ 0.8. The velocity map at the corresponding positions shows that these patches have slower velocities compared to the wind surrounding them. Since they do not propagate as fast as the major part of the wind, they should lag behind the main portion of the flow. This is indeed seen in the lower panel of the wind fraction map which is a snapshot at a later time frame. The development of these mixing structures shows that the ambient material can be entrained at the strongly sheared wind-ambient interface where Kelvin-Helmholtz instability is active. Occasionally, larger blobs of ambient materials can also get incorporated into the wind and subsequently get ejected along the flow. Given the much higher density of the inner toroid than the wind, the entrained ambient mass can dominate the local mass fraction for some time before being diluted by the wind. We discuss this in more detail in section 5.2. –7– 4.2. The Joint Effect of Wind Temperature And Bφ A group of simulations has been carried out to investigate effects of the wind temperature and the wind magnetic field on the outflow morphology and dynamics. We present the resulting density snapshots at t = 1000 yr in Figure 2. Within this three-by-three mosaic, the wind temperature increases from ∼10 K in column 1 to ∼400 K and ∼2000 K in column 2 and 3, corresponding to awind = 2 × 104 , 1.2 × 105 , and 2.67 × 105 cm s−1 , respectively. Note that the sound speed of the toroid is 2 × 104 cm s−1 , so the three cases in column 1 are in fact globally isothermal. The strength of the toroidal magnetic field in the wind, on the other hand, decreases form MA = 6 on the top row (a) to MA = 30 and 90 in the second and the third row (b and c). The collection of awind —MA combinations are meant to explore a wide range of physical conditions. With a high field strength and a low temperature, the upper-left case a1 represents the magnetic-dominated situation with negligible gas pressure. The lower-right panel c3 shows the contrary case where thermal energy dominates. For a temperature of ∼ 400 K and MA = 30 in case b2 (βplasma ≃ 1), both magnetic field and thermal pressure play a role. Case b2 is also the demonstrative case discussed in the previous section. We directly list the value of awind and MA above each panel for easy reference. How raising the temperature can affect the wind is best illustrated with the three cases of relatively weak magnetic fields, namely, those in the third row of Figure 2. When the wind temperature increases from ∼10 K (c1) to ∼2000 K (c3), the change in the jet component is evident. The initial density distribution of the wind is cylindrically stratified and takes the form of ρ ∝ ̟−2 (see equation 1), which is strongly peaked toward the outflows axis. In case c1 this jet-like structure is clearly seen and is well maintained along the flow. At higher temperatures, however, we begin to see how the strongly peaked density profile gradually flattens with distance. In c2 the jet component becomes rather flat after reaching z ∼ 4000 AU, and in c3 it flattens so rapidly that beyond z ∼ 2000 AU the jet identity is barely recognizable. The rapid drop of the axial density peak in c2 and c3 is a result of thermal pressure gradient force which pushes the gas from high density region toward the lower ones. The strength of the force scales linearly with temperature and increases over two order of magnitude from 10 K to 2000 K, resulting in the strongly modified jet density profile. Another notable change is that the opening angle of the shell component becomes larger as temperature increases. As the figure shows, the lateral size of the swept-up shell is about twice as large in c3 than in c1. On the other hand, the effects of varying toroidal magnetic field strength in the wind are most evident when we compare the simulations at a same high wind temperature of T ∼ 2000 K, namely, the three cases in the third column of Figure 2. The weakly magnetized case c3 presents a jet component that flattens rapidly with distance due to the strong pressure gradient force. In case b3 where the wind is more strongly magnetized, we tentatively see that the jet is maintained to a larger distance, although the map is complicated by some substructures. In the most strongly magnetized case a3 where the toroidal field is more than one order of magnitude stronger than in c3, the initial jet density structure turns out to be well-maintained despite the high wind temperature of ∼ 2000 K. Since the three simulations start from the same initial wind density profile and evolve –8– under the same temperature, the pressure gradient force responsible for smoothing the jet in c3 must also be active in a3. It is only because the stronger toroidal magnetic field in a3 has countered the effect of thermal pressure that the stratified density profile of the jet is maintained. This same collimation effect has been shown in paper I for the cold flow, and for a higher wind temperature the effect is even more dramatic. We have shown that the temperature and magnetic field each plays a role in determining the fate of jet-like density structure. For example, in contrast to the clear difference between c1 and c3 due to the increased temperature, the counterpart pairs of a1 and a3 appear unaffected by the temperature because of their stronger toroidal field. The joint effect of the two factors can be better understood with an overview of all nine cases in Figure 2. While the temperature varies column by column and the toroidal field strength row by row, a general trend is found in direction of upper-left to lower-right. As we move from a magnetic-dominated regime to a thermal-dominated one, the peak density of the jet-like structure drops more rapidly with distance. To compare the different degrees of flattening, in the upper panel of Figure 3 we present on a log-log scale the density profiles of all nine cases cut perpendicular to the jet at a height of z = 200 AU. The three profiles of the most strongly magnetized cases (MA = 6, solid lines) are very similar despite the different temperatures. Case a1 is almost indistinguishable from a2, and a3 only shows a small discrepancy from the other two. They all preserve the two order of magnitude drop in density from ̟ = 1 to 10 AU obtained from the initial wind solution where ρ ∝ ̟−2 . For the less magnetized cases of MA = 30 (dashed lines) and MA = 90 (dotted lines), the flattened density profiles are clearly seen for the 400 K and 2000 K cases, but not for 10 K. In fact, all three cases at 10 K (a1, b1, and c1) as well as the strongly magnetized cases (a2 and a3) all show similar density profiles in this figure. We also note that a less magnetized wind with a modest temperature of ∼ 400 K (the intermediate case b2) shows only a mild density drop. Although the density profile is modified by the thermal pressure, the change is slow enough that its jet identity is still clearly seen up to a distance of at least 10000 AU from the source (see Figure 1). The jet identity is strongly modified only in case of a poorly magnetized wind at high temperature. The jet structure is clearly related to both the temperature and the toroidal magnetic field strength. How the relative importance of the two factors is related to the density drop is examined in the lower panel of Figure 3. We use the wind density at ̟ ≃1 AU in each profile as a proxy to the degree of jet flattening and arrange the results using their βplasma values. Different symbols are used for different magnetic field strengths, while different wind temperatures are represented by different colors. Beside the nine cases that have been discussed, we also include the results of non-magnetized wind at the three temperatures explored (the cross symbol on the right boundary of the plot). The three extra points can be thought of as asymptotic results one would obtain when the magnetic field strength is decreased indefinitely for a fixed wind temperature. The plot shows that for magnetic-dominated cases the axially concentrated profiles are well maintained and the density values are about the same (a1, b1, c1, a2, and a3). In the intermediate regime the drop in jet density becomes apparent, and the value gradually approaches what one would obtain in –9– the hydrodynamic limit as thermal pressure becomes more and more dominant. In absence of the toroidal field, the jet follows a simply trend where it flattens more significantly for a higher wind temperature, as shown by the values of the three cross symbols. For different temperatures, the toroidal magnetic field strengths sufficiently high to collimate the flow are different. But once the field is strong enough the jet density can be kept at a similar high value. Although the scattering of the data points shows that jet structure is not simply determined by the relative strength of the thermal and magnetic energy in the wind, with this ratio we are able to show the opposing roles of the two. The opposing roles of the magnetic and thermal pressure on the stratified density structure is closely related to the properties of the wind solution. The density distribution of the wind takes the form of ρ(̟) = D0 /̟2 , which naturally produces a pressure gradient force pointing in the positive ̟ direction. Since its strength is proportional to both the gas temperature and the density gradient, it is strongest near the axis where the gradient is high and is more significant at higher temperatures. On the other hand, with the initial configuration of Bφ (̟) ∝ ̟−1 , the toroidal field exerts no Lorentz force. But when the thermal pressure pushes the centrally concentrated gas outward, it also modifies the toroidal field distribution and as a consequence a magnetic force can arise. In a strongly magnetized case a slight modification in the field configuration can result in a strong magnetic force, making it more capable of countering the flattening of density peak caused by the thermal pressure. 4.3. The Role of Ambient Poloidal Fields The ordered poloidal magnetic field assumed to support and shape the toroid in the model of Li & Shu (1996) has been ignored in the simulations so far. As in paper I, no poloidal field is set up and only the wind-originated toroidal field is evolved for simplicity. It has, however, been anticipated that the presence of a poloidal field in the ambient could help resist lateral expansion of the wind and therefore aid in the shell collimation. In this section we present simulation results conducted with the poloidal magnetic fields initialized to the appropriate toroid solution. Both toroidal and poloidal fields are evolved. In Figure 4 we present the density snapshots at t = 1000 yr of a series of simulations with ambient poloidal fields. The adopted combinations of wind temperatures and toroidal field strengths are the same as those in Figure 2 so that for each case in Figure 4 there is a counterpart in Figure 2. A comparison between the two sets of simulations thus straightforwardly shows the effects of the poloidal field. The most evident change is clearly seen in the swept-up shell: a much smoother and thicker structure is formed in place of the original corrugated interface. Such effect is quite general since it occurs in all nine cases from magnetically to thermally dominated and does not depend on the wind properties. On the other hand, we find that the ambient field has little or negligible effects on the evolution of axial density concentration of the jets. This is expected since the central jet component propagates almost freely within the outflow cavity and knows little about the ambient. – 10 – The shell, on the contrary, is itself a result of wind-ambient interaction and is directly affected by the ambient poloidal field. We examine the outflow structure in more detail again using case b2 as an example. In Figure 5 we compare maps of density and poloidal velocity produced either with ambient field ignored (upper row ) or included (lower row ). The poloidal field lines are shown with thin white contours in the lower panels. To illustrate the region within which the wind is confined, we overlay the contours of 1% wind fraction with thick solid lines (dark in the left panels and red in the right). Beyond this line the wind have negligible mass fraction, and we take this as a proxy of the wind boundaries. By comparing the upper and lower density maps we find only tiny differences in the enclosed regions. At a given outflow height z, the cylindrical radius ̟ of the boundary is marginally larger when the ambient poloidal field is absent. Although this appears compatible with the expectation that an ambient field can contribute to the lateral confinement of the wind, the effect is not significant for the provided poloidal field strength. The modified shell structure in the presence of an ambient poloidal field can also be examined more clearly in Figure 5. In the upper density map we can easily recognize the thin, high-density layer lying just outside the wind boundary (dark solid line) as the swept-up shell. It consists mainly of piled up ambient material that has received momentum from the wind. In the lower panels, however, a smooth parabolic structure is formed instead of a corrugated thin boundary. This shell-like density structure also consists mainly of ambient material, but it is now thicker and appears to be threaded by the poloidal magnetic field (white thin lines). Therefore, by contrasting the two different sets of simulations, we find that the poloidal field can help prevent corrugation and results in a smoother interface. This agrees with the expectation that the presence of a magnetic field along a sheared flow can suppress the growth of instabilities and results in a more organized flow behavior (see e.g., Chandrasekhar 1961; Frank et al. 1996, and references therein). In addition, when the ambient field lines bunch up in the toroids as the wind pushes them aside, the enhanced magnetic pressure tends to resist compression and prevents the density layer from being all squeezed together. The field thus acts like a cushion within the toroid. Finally, we examine the magnitude of poloidal velocity in the right column of Figure 5. The red dashed lines are contours of |vpoloidal | = aambient , and the flow directions are shown with yellow arrows for positions with velocity magnitude higher than this value. The same contours are overlaid on the density maps with dark dashed lines as well. It turns out that this contour well delineate the boundary between the shell structure and the rest of the toroid. The reason could be that material not yet affected by the wind would remain nearly static in our simulations, while that already affected can have a finite speed. In the lower right panel, one finds that the velocity arrows in the shell point roughly in a direction perpendicular to the poloidal fields. This is consistent with the idea that the gas in the expanding shell is accelerated by the magnetic force. In summary, although the poloidal field threading the toroid has only limited effect in terms of confining the wind, it plays a significant role in wind-ambient interaction by modifying the resulting shell structure. This plot also shows more clearly that the strong shear happens near the wind boundary (thick solid lines) – 11 – at the inner edge of the compressed shell. 5. 5.1. DISCUSSION Implications on The Jet And Shell Accumulating observations have revealed details in many molecular outflows. While bearing the same basic characteristics of jet and shell, the variety of cavity shapes, collimation degree, velocity patterns, etc. suggests some underlying differences among different sources. The rich morphologies and kinematics could be potential probes to the underlying physical conditions or evolution histories that are not well understood. Nevertheless, we need to understand the possible influences of different physical factors in order to learn from these clues. In view of its success in reproducing the basic jet and shell characteristics, we adopt the framework of the unified wind model to explore possible signatures of temperatures and magnetic fields. Several findings from the series of numerical experiments presented in the previous section may help illuminate the direction. The central portion of the cylindrically stratified wind in the model is identified with the extremely high velocity jet component (Bachiller 1996) observed in young molecular outflows . Our simulations suggest that this density profiles can be modified by thermal pressure when the temperature is high and the toroidal magnetic field is not too strong. In the most extreme case such as c3 in Figure 2, the jet feature flattens very quickly. It is difficult to identify a collimated jet beyond a few thousand AU from the origin. Therefore, the presence of a clear collimated jet would suggest either the flow is cold enough that pressure is negligible, or the toroidal field is strong enough to suppress the effect of gas pressure if the wind is hot. It may be worth noting that since a direct detection of the magnetic field in the wind is difficult, a highly collimated jet-like structure found to persist at a high temperature could be a clue of the presence of dynamically important toroidal field. More quantitative prediction would rely on comparisons between model and observations, which is however expected to be complicated by the knotty nature of the real sources. The simulations also pose some interesting questions on the shell structure. In the shell model of Shu et al. (1991), molecular outflows are the ambient material driven into expansion by the momentum of a radial wind. The resulting morphology of the shell is determined by a bipolarity function which considers the directional input of the wind and the anisotropic distribution of the ambient. Paper I has numerically demonstrated the formation of such swept-up shells, and the results appear to be in good agreement with the low velocity component observed in young outflows such as HH 211 (Gueth & Guilloteau 1999), L1448C (Hirano et al. 2010), IRAS 04166+2706 (Santiago-Garcı́a et al. 2009; Wang et al. 2014) etc. With much improved spatial resolution over paper I, the new simulations can resolve instabilities on a smaller scale. While the overall size and extent of the outflow are not affected by the resolutions, we find corrugated shells in many cases as a result of the finer grids as shown in Figure 2. Since the real world wind and toroid density – 12 – distributions are unlikely as smooth as the analytic description adopted in the model, structures of the swept-up shell can in principal be more irregular than what is shown in the current simulations. This, however, is based on the results obtained before including the ambient poloidal field in the simulations. In Figure 4, it is clearly shown that the presence of such field component can provide support to the boundary and maintain a much smoother interface between the wind and the ambient. It then appears that a smooth boundary of the low velocity shell could actually suggest the presence of the ambient poloidal field, since otherwise a more corrugated structure may be expected. To tell whether these implications may be probed by observations would require further predictions of observational features from the model. For example, synthetic position-velocity diagrams could be used to interpret the kinematics in a more observational fashion. As an example, in Figure 6 we show the CO J = 3–2 position-velocity diagram cut along the outflow axis of case b2 (without poloidal field). A gas temperature of 400 K and a condition of local thermal equilibrium is assumed while calculating the emission. At an inclination angle of 30◦ from the plane of the sky, the axially concentrated jet-like structure shows up with a projected radial velocity of ∼25 km s−1 in the plot. A wider velocity distribution is seen near the base (the origin), which is a characteristic of a wide angle wind. On the other hand, the swept-up material that constitutes the shell has a much lower velocity and its emission is largely found near the system velocity which is taken to be zero. How the different physical conditions may introduce observable signatures would need to be examined more carefully and will be addressed in an upcoming paper. 5.2. Instability And Mixing Between Wind And Toroid One question about molecular outflows is whether they are a part of the primary wind or are made of entrained ambient material. With the tracer field, we can examine the wind fraction at each position to address this problem. We have seen in Figure 1 that the major portion of the flow (including the central jet-like structure) consists of wind-originated material. On the other hand, the shell structure formed at the wind-ambient interface has been shown to consist mainly of swept-up ambient material. Therefore in the context of the unified wind model, the simple answer to the question would be that a molecular outflow is a primary-wind-originated jet surrounded by a swept-up shell of ambient material. For a more detailed picture we need to look into the mixing region. As described in section 4.1, mixing happens within two kinds of regions: a thin transition layer between the wind and the ambient, and blobs and feather-like substructures. The mixing substructures are mostly a result of dense toroid material entrained and carried by the wind, and we illustrate their development in Figure 7 using case b2 (without ambient poloidal field) as an example. We plot the wind fraction at the inner few hundred AU for three different epochs t = 500, 550, and 600 yr. Stochastic, corrugated structures are found to evolve at the dense base of the toroid as well as on the interface. Blobs of ambient gas, big or small, can be entrained by the wind and subsequently carried by the – 13 – flow. Given the strong velocity shear at the wind-ambient interface, the corrugated boundary is likely related to the Kelvin-Helmholtz instability. We have also mentioned in section 4.3 that the presence of poloidal field can help resist the corrugation of the swept-up shell. In Figure 8 we present wind fraction snapshots similar to those of Figure 7 for the simulations including poloidal field. The field lines are plotted with white contours. A comparison to Figure 7 again shows how corrugation is suppressed at the interface and a much smoother boundary is found. Note, however, that even in the presence of a poloidal field, there are still feather-like mixing structures growing inside the wind region. One reason could be that the poloidal field which helps stabilize the shear instabilities is only present in the toroid and not in the wind. Another could be the much higher density of the toroids as compared with the wind. Any wind material running into the toroids will not go too far before losing its momentum to the dense swept-up shell. The wind boundary therefore appears largely well defined despite the complexity of the substructures. Various types of models have been proposed to explain bipolar molecular outflows, and among which wind-driven and jet entrainment are two popular mechanisms (Cabrit et al. 1997). In the jet-driven scenario, one key question is how a collimated jet can drive a less collimated molecular outflow. Since a turbulent mixing layer between the jet and the ambient is likely thin (e.g., Cantó & Raga 1991), entrainment through a jet bowshock (e.g., Raga & Cabrit 1993) appears to be more promising. Here we would like to note that the setup of our wind simulation is sufficiently different from the jet bowshock scenario so that bowshock entrainment is irrelevant here. In the jet bowshock entrainment model, the molecular outflow is identified with the environmental gas entrained into the wake of the bowshock, whereas in the wind scenario ambient material is swept up by the primary wind. There is no prominent bowshock in our wind simulations, which is partly due to the lack of ambient material in the polar region of the toroid solution. With the tracer field, we have also seen that the mixed material within the wind region mainly results from entrainment and/or ejection at the wind boundary, especially near the inner dense toroid. We do not find significant ambient material being incorporated or entrained near the head (the polar region) of the outflow. In our wind simulations, the wide opening angle of the shell is a natural consequence of the wide-angle wind and the toroidal distribution of the ambient density, and has little to do with the entrainment process going on at the sheared interface. With improved resolution, in this work we have resolved corrugated swept-up shells as well as small blobs of mixed material in the wind region that were not seen in paper I. The swept-up shells are much smoother in the simulations including poloidal field, which is related to the stabilization of shear instability by the magnetic field. To check how an even higher resolution can affect the results, we have run simulations for case b2 at a resolution of 3200×990. In the upper panel of Figure 9 we present the density map at t = 1000 yr for the run without ambient poloidal field. A close comparison to its ordinary resolution counterpart (see the upper left panel of Figure 5) shows that the details of the mixing structures are quite different in the two runs. More complex and fragmented mixing structures are generated in the higher resolution simulation, which indicates that the process producing these structures is not completely resolved at the current resolution – 14 – (see, e.g., Frank et al. 1996). Despite the different fine structures, the overall outflow shape and the jet structure in the two simulations are in good agreement. In the lower panel of Figure 9 we present the same maps for the run with ambient poloidal field. The smoother and thicker shell seen here is almost identical to its ordinary resolution counterpart (see the lower left panel of Figure 5), and similar feather-like structures are also present near the wind boundary. There is also no sign of intense stochastic behavior like that seen in the upper panel of Figure 9. This suggests that the stabilization of the wind-ambient interface by the poloidal field is quite robust, and this conclusion holds well at least for the resolution we have used. Finally, it should be noted that the wind in the current simulations fans out in all radial directions even at the equator where an accreting inner envelope or disk is instead expected. Whether this setup could have anything to do with the growth of stochastic mixing structures near the inner toroid will need to be resolved in future simulations with more realistic setup of boundary conditions. 6. CONCLUSION In this work, we aim to understand the relation between outflow characteristics and their underlying physical factors using MHD numerical simulations in context of the unified wind model (Shang et al. 2006). In the wind model, the jet and shell characteristics common to young molecular outflows are results of the intrinsic stratification of wind density profile and the swept-up materials during the wind-ambient interaction. Their detailed morphologies and kinematics may contain information about the underlying physical conditions of the wind. We study the consequences of varying temperature and toroidal magnetic field in the wind using the ZEUS-TW code. A series of runs are carried out using a tracer field technique to keep track of the wind material and to enable the two-temperature prescription adopted. While keeping the ambient toroids unchanged, different combinations of wind temperatures (sound speed awind = 2 × 104 , 1.2 × 105 , 2.67 × 105 cm s−1 ) and wind Alfvén Mach numbers (MA = 6, 30, 90) have been investigated. We also examine the role of poloidal magnetic field of the toroid by comparing runs with and without including the ambient poloidal fields. Our findings are as follows: 1. For a less magnetized wind, the cylindrically stratified density profile that forms the central jet component can be modified by the thermal pressure at a high wind temperature. In that case, as the wind temperature increases, the wind gas pressure can drive material away from the dense axis so that the density profile will gradually flatten as the wind propagates. In extreme cases, the jet-like structure is barely identifiable beyond several thousand AU (e.g., case c3 in Figure 2). 2. The toroidal magnetic field, on the other hand, tends to maintain the initial density structure of the wind and appears to oppose the change made by the thermal pressure. For a strongly magnetized wind (such as MA = 6), it dominates over thermal pressure even at a temperature as high as ∼ 2000 K, and the jet structure appears unaltered as it propagates away from the – 15 – source (case a3 in Figure 2). 3. Therefore in a magnetized wind it is the magnetic field and thermal pressure that jointly dictate the evolution of the wind, especially the jet-like density profile. For a sufficiently magnetized wind, the jet identity is well maintained by the magnetic collimation effect. If the field is not strong enough to oppose the effects of thermal pressure, the density profile will evolve and become less axially concentrated as the flow propagates. Such effect is more significant for a higher wind temperature. Arranging the results in terms of the relative importance of thermal and the magnetic pressure may help understand the opposing roles of the two factors (Figure 3). For the intermediate case b2, the change in the density profile is mild and the jet-like structure is still clearly seen up to a distance of at least 10000 AU. 4. The poloidal magnetic field that threads the toroid is found to play a role in the windambient interaction. When threaded by the field lines, the toroids become more resistant to corrugation and a smoother wind-ambient interface is formed. Also, as the wind pushes the toroids aside, field lines bunch up and the magnetic pressure is enhanced, making the region more resistant to compression. The swept-up shell appears thicker as compared to the thin layer obtained in the absence of poloidal field. The cushioning effect of poloidal field thus leads to a smoother and thicker swept-up shell (Figure 5). In short, our explorations suggest that the density contrast of the jet is closely related to the relative strengths of the magnetic field and the thermal pressure in the wind. We also find that an ambient poloidal field could be the key factor to form a smooth and thicker swept-up shell. The variety of results also demonstrate the potential capability of the model in accommodating the rich outflow phenomenon observed. While the basic principals are clearly demonstrated, whether they can help constrain physical conditions of molecular outflows in the real world will depend on more detailed predictions of observable diagnostics. There are several possible directions to pursue in the future work. For example, it would also be interesting to check whether non-ideal MHD effects will result in different behaviors, especially for the wind-ambient interactions. This work was supported by funds from the Theoretical Institute for Advanced Research in Astrophysics (TIARA) through the Academia Sinica and from the Ministry of Science and Technology of Taiwan by MOST 102-2119-M-001-008-MY3. REFERENCES Bachiller, R. 1996, ARA&A, 34, 111 Cabrit, S., Raga, A., & Gueth, F. 1997, in IAU Symp. 182, Herbig-Haro Flows and the Birth of Stars, ed. B. Reipurth & C. Bertout (Dordrecht: Kluwer), 163 – 16 – Fig. 1.— Maps of (left to right) density (cm−3 ), poloidal velocity magnitude (km s−1 ), tracer field (wind density ρwind ) (cm−3 ), and wind fraction for the reference case of MA = 30 and awind = 1.2 × 105 cm s−1 (T∼ 400 K). The upper and lower rows are snapshots at t = 100 and t = 1000 yr, respectively. The white and gray lines overlaid on the wind fraction maps are 10% and 90% contours. – 17 – Fig. 2.— A three-by-three mosaic of density snapshots at t = 1000 yr for nine different cases. The wind temperature increases from ∼ 10 K in column 1 to ∼ 400 K and ∼ 2000 K in column 2 and 3, corresponding to awind = 2 × 104 , 1.2 × 105 , and 2.67 × 105 cm s−1 , respectively. The strength of the toroidal magnetic field in the wind decreases form the strongest MA = 6 on the top row (a) to MA = 30 and 90 in the second and the third row (b and c). The results demonstrate the effect of varying magnetic field and wind temperature in the wind. The jet structure is well-maintained in the magnetic-dominated regime, and flattens with distance in magnetic-dominated regime. – 18 – Fig. 3.— Top: density profile cut perpendicular to the outflow axis at z = 200 AU for the nine cases presented in Figure 2. The modification of the cylindrically stratified density profile by thermal pressure is clearly seen. Note that the lines for cases a1, b1, c1, a2, and a3 are very close or even overlap with each other in this plot. Bottom: wind density sampled at ̟ ≃ 1 AU plotted against the βplasma value for each case. Data points with the same wind temperature are connected using dashed lines. We also include results of non-magnetized wind simulations at the three temperatures explored (the three cross symbols on the right boundary of the plot), which can be thought of as βplasma = ∞. In magnetic-dominated regime, the initial density is well maintained by the magnetic collimation. In thermal-dominated regime, the density approaches what one would obtain in a non-magnetized wind. – 19 – Fig. 4.— A three-by-three mosaic of density snapshots at t = 1000 yr for nine simulations with ambient poloidal magnetic fields. The same combinations of awind and MA as in Figure 2 are used. The swept-up shell is clearly affected by the presence of an ambient poloidal field. – 20 – Fig. 5.— Density and velocity snapshots at t = 1000 yr comparing simulations ignoring (upper ) and including (lower ) the ambient poloidal fields. In this example the Alfvén Mach number MA is 30 and the wind sound speed awind is 1.2 × 105 cm s−1 (∼ 400 K). The thick solid lines are contours of 1% wind fraction, and the poloidal magnetic fields are plotted with thin white contours in the lower panels. The thick dashed lines are velocity contours of |vpoloidal | = 2 × 104 cm s−1 , which is the ambient sound speed aambient . The flow directions are indicated by arrows (only showing those with velocity greater than aambient ). The modification of shell structures by the presence of poloidal magnetic field is clearly seen. – 21 – Fig. 6.— A synthetic position-velocity diagram cut along the outflow axis of case b2 (without poloidal field). The calculation of CO J = 3–2 line emission assumes a gas temperature of 400 K and local thermal equilibrium. With an inclination angle of 30◦ from the plane of the sky, the dense axially concentrated jet material mainly has a projected velocity of ∼25 km s−1 . On the other hand, the swept-up material at low velocity forms another feature near the system velocity which is taken to be zero. – 22 – Fig. 7.— Snapshots of wind fraction map of case b2 (without ambient poloidal magnetic field) at t = 500, 550, and 600 yr to demonstrate the process of mixing between the wind and the ambient. Stochastic structures are found to develop in the inner (z < 100 AU) dense portion of the toroid during interaction with the wind. The blobs are subsequently ejected along the flow. Given the strong shear between wind and ambient, Kelvin-Helmholtz instability could aid in producing the corrugated interface. – 23 – Fig. 8.— Snapshots of wind fraction map of case b2 (with ambient poloidal magnetic field included) at t = 500, 550, and 600 yr to demonstrate the process of mixing between the wind and the ambient. The poloidal field lines are plotted with white lines. Stochastic structures are still seen to develop in the inner dense portion of the toroid, but appear to be less significant. With the presence of poloidal magnetic field in the toroid, the interface is more stable and less corrugated compared with Figure 7. – 24 – Fig. 9.— Density maps of high resolution simulations for case b2 at t = 1000 yr. The resolution is 3200×990, which is higher than the ordinary 1600×720 used throughout this work. The upper and lower panels show the run without and with ambient poloidal field, respectively. The poloidal magnetic field is plotted with thin white lines in the lower panel. In the absence of poloidal field, details of the stochastic structures appear quite different from its ordinary resolution counterpart (upper left panel of Figure 5). The run with ambient poloidal magnetic field, on the other hand, appears quite similar to its ordinary resolution counterpart (lower left panel of Figure 5). – 25 – Cantó, J., & Raga, A. C. 1991, ApJ, 372, 646 Chandrasekhar, S. 1961, International Series of Monographs on Physics, Oxford: Clarendon, 1961, Frank, A., Jones, T. W., Ryu, D., & Gaalaas, J. B. 1996, ApJ, 460, 777 Gueth, F., & Guilloteau, S. 1999, A&A, 343, 571 Hirano, N., Ho, P. P. T., Liu, S.-Y., et al. 2010, ApJ, 717, 58 Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2010, ApJ, 716, 1541 Lee, C.-F., Mundy, L. G., Reipurth, B., Ostriker, E. C., & Stone, J. M. 2000, ApJ, 542, 925 Li, Z.-Y., & Shu, F. H. 1996, ApJ, 472, 211 Nisini, B., Codella, C., Giannini, T., et al. 2007, A&A, 462, 163 Raga, A., & Cabrit, S. 1993, A&A, 278, 267 Santiago-Garcı́a, J., Tafalla, M., Johnstone, D., & Bachiller, R. 2009, A&A, 495, 169 Shang, H., Allen, A., Li, Z.-Y., et al. 2006, ApJ, 649, 845 Shu, F. H., Najita, J. R., Shang, H., & Li, Z.-Y. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Tucson, AZ: Univ. Arizona Press), 789 Shu, F. H., Ruden, S. P., Lada, C. J., & Lizano, S. 1991, ApJ, 370, L31 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791 Wang, L.-Y., Shang, H., Su, Y.-N., et al. 2014, ApJ, 780, 49 This preprint was prepared with the AAS LATEX macros v5.2.