MAGNETIC NESTED-WIND SCENARIOS FOR BIPOLAR OUTFLOWS: PREPLANETARY AND YSO NEBULAR SHAPING^*

We execute a series of axisymmetric, radiatively cooled simulations of two coaxial, or "nested," MHD and HD winds flowing simultaneously into a rectangular domain from a source located at the lower left boundary. The long (z) axis is chosen to correspond to a physical size equal to $2\times 10^4 \rm \ AU,$ while the short (r) axis is $8\times 10^3\rm \ AU$. The direction of flow is predominantly along z. A diagram of the boundary condition assumed in the region where the winds are launched is shown in Figure 1. Two radii, r_d and r_s, are defined to delineate the regions where the inner "stellar" wind and the outer "disk" wind operate. In the range r < r_s of the launch region, the parameters specifying the stellar wind determine the properties of the flow there, while for r_s < r < r_d the flow properties are determined by the parameters specifying the disk wind. The parameters chosen to specify the winds include the disk-wind and stellar-wind mass-loss rates: $\dot{M}_{\rm d}$, and $\dot{M}_{\rm s}$; their velocities: v_d and v_s; their opening angles: θ_d and θ_s; and two dimensionless parameters: β_m and σ to be discussed below. As indicated in Figure 1, the wind opening angles are upper limits on the degree to which the flow directions depart from being parallel with the z-axis. For a cell in the launch region, a distance r from the origin, the velocity vector for the flow emerging from that element makes an angle (r/r_w)θ_w with respect to the z-axis, where θ_w is one of θ_s or θ_d, and r_w is one of r_s or r_d depending on whether r < r_s or r>r_s, respectively. The opening angles are both nonzero with the disk wind opening angle shallower than that of the stellar wind. The intent is to choose opening angles here that accord with the notion that the disk wind is launched magnetocentrifugally and "flung" out along poloidal field lines, while the stellar wind mass loss is more nearly isotropic. The angles chosen for the simulations presented here represent our best effort to model this scenario while respecting the technical limitations imposed upon us by the code.

Zoom In Zoom Out Reset image size

Throughout the launch region, at each time step, a toroidal magnetic field is embedded in the winds. To characterize the strength and dynamical significance of this field, we introduce the independent parameters σ and β_m, where σ is the ratio of wind magnetic energy density to wind kinetic energy or,

$$\begin{equation} \sigma =\frac{B^2/8\pi }{\rho _{\rm w}v_{\rm w}^2}, \end{equation} \tag{ 1 }$$

and β_m is a particular value of the ratio of thermal pressure to magnetic pressure, chosen to be characteristic of the wind, i.e.,

$$\begin{equation} \beta _{\rm m}=8\pi P_{\rm w}\big/B_{\rm m}^2, \end{equation} \tag{ 2 }$$

where B_m is a maximum value for the magnetic field strength. By specifying the values of σ and β and constraining the value of wind mass density, ρ_w, by requiring a relation of the form

$$\begin{equation} \dot{M} = \Omega \rho _{\rm w}v_{\rm w}r_{\rm w}^2 \end{equation} \tag{ 3 }$$

(where Ω is the solid angle of the wind) to hold among the mass-loss rates, velocities, and radii of the stellar and disk winds⁵, we infer a characteristic value for the thermal pressure of the wind P_w by way of the relation

$$\begin{equation} P_{\rm w}=\rho _{w}v_{\rm w}^2\sigma \beta, \end{equation} \tag{ 4 }$$

and fix the value of B_m, from the definition of β_m. That is,

$$\begin{equation} B_{\rm m}=\sqrt{8\pi P_{\rm w}/\beta _{\rm m}}. \end{equation} \tag{ 5 }$$

Finally, using the values of P_w and B_m thus obtained we model the magnetic and pressure profiles B(r) and P(r) for the winds after the form first introduced by Lind et al. (1989). For the magnetic field profile, we write

$$\begin{equation} B\left(r\right)=\Biggl \lbrace \begin{array}{@{}ll} B_{\rm m}\frac{r}{r_{\rm m}}, & 0\le r<r_{\rm m}, \\ \ms B_{\rm m}\frac{r_{\rm m}}{r}, & r_{\rm m}\le r<r_{\rm d}, \\ \ms 0, & r\ge r_{\rm d}, \end{array} \end{equation} \tag{ 6 }$$

and for the pressure profile, we write

$$\begin{equation} P\left(r\right)=\left \lbrace \begin{array}{@{}ll} \displaystyle \left[\alpha +\frac{2}{\beta _{\rm m}}\left(1-\frac{r^2}{r_{\rm m}^2}\right) \right]P_{\rm w}, & 0\le r<r_{\rm m}, \cr \alpha P_{\rm w}, & r_{\rm m}\le r<r_{\rm d}, \cr P_{\rm w}, & r \ge r_{\rm d}. \end{array}\right. \end{equation} \tag{ 7 }$$

The quantity α is a constant related to β_m according to

$$\begin{equation} \beta _{\rm m}^{-1}=\left(1-\alpha \right)\left(r_{\rm d}/r_{\rm m}\right)^2, \end{equation} \tag{ 8 }$$

and r_m is the value of r at which B = B_m, and is chosen in all of the simulations presented here to be equal to the "stellar" radius r_s. One may verify that with the profiles defined as in Equations (6) and (7), the pressure and magnetic field satisfy the condition of magnetostatic equilibrium, i.e.,

$$\begin{equation} \frac{dP}{dr}=-\frac{B}{4\pi r}\frac{d(rB)}{dr}, \end{equation} \tag{ 9 }$$

everywhere in the interior of the winds except at r = r_m, where the value of P_w (and therefore B_m) changes. The winds are launched into a homogeneous, unmagnetized ambient medium whose total number density n_a and temperature T_a are chosen independently. As a consequence, the winds are not pressure-matched with respect to their environments in the simulations. We allow this freedom so that we may explore how shaping is affected by changes in the circumstellar environment (see Section 3.3).

Lastly, we note that since all of the parameters σ, β_m, $\dot{M}_{\rm s}$, $\dot{M}_{\rm d}$, v_d, v_s, r_d, and r_s are set independently, the simulations are not controlled with respect to the total power in the winds. However, because the field strength in the winds has been kept relatively low, the magnetic contribution to the outflow power, $\mathcal {P}_{\rm B}$, is small ($\mathcal {P}_{\rm B}\lesssim \sigma \mathcal {P}_{\rm tot}$) in all of the simulations presented, facilitating the comparisons below and arguing against the possibility that the differences seen between the MHD and HD runs are the result of a dominant energy effect.

We have executed nine, radiatively cooled, axially symmetric simulations using "AstroBEAR." AstroBEAR is an AMR Hydro/MHD code based on the conservative form of the MHD equations and designed for use with high-resolution shock capturing methods for sets of nonlinear hyperbolic equations. We use AstroBEAR to solve either the ideal MHD equations for the magnetized winds or the equations of inviscid hydrodynamics (i.e., the Euler equations) for the unmagnetized winds. In both cases, we assume that the field variables do not depend on the azimuthal angle so that our solutions are axially symmetric. This allows us to reduce the problem to a two-dimensional calculation. Because the solutions that result represent "slices" through the axis of symmetry of the full three-dimensional axisymmetric solution, they are said to be 2.5 dimensional (2.5D). For a description of AstroBEAR and the equation set it solves, see Cunningham et al. (2009).

Short descriptions of each simulation and their parameterizations are given in Tables 1 and 2. Table 1 lists the values of those parameters that are common to all cases, and Table 2 lists each run and its corresponding parameterization. All but one of these simulations was carried out on a grid of base resolution 640 × 256 and with two levels of refinement for an effective resolution 2560 × 1024. One simulation (run A2) was carried out on a fixed grid of resolution 2560 × 1024. In every case the discretization used corresponds to a resolution of 64 cells per disk-wind radius or 16 cells per stellar-wind radius. The standard case, in which two MHD winds interact as suggested by Figure 1, is presented as run A1. For purposes of comparison, we also present an HD case, cases for which the outermost of the two winds is absent, and cases for which the stellar wind is "lighter" than the disk wind (i.e., the stellar mass-loss rate is an order of magnitude smaller). Two more cases are also presented for nested HD and MHD winds in order to compare the morphologies obtained in a warm (${\sim} 10^4 \rm \ K$) circumstellar environment with those obtained in a cool (${\sim} 10^2 \rm \ K$) environment.

Table 1. Parameters Common to All Simulations

Effective Resolution	16 cells/rs
σ...	0.1
β...	1.0
na...	$5\times 10^3\rm \ cm^{-3}$
θs...	30°
θd...	15°
rs...	$125\rm \ AU$
rd...	$500\rm \ AU$

Download table as:  ASCIITypeset image

Table 2. Simulation Parameters

Run	Description	$\dot{M}_{\rm s}$	$\dot{M}_{\rm d}$	vs	vd	Ta	Run Time
		(M☉ yr−1)	(M☉ yr−1)	(km $\rm \ s^{-1}$)	(km $\rm \ s^{-1}$)	(K)	($10^3\rm \ yr$)
A1	MHD, both winds on	10−7	10−7	150	50	104	1.3
B1	HD, both winds on	10−7	10−7	150	50	104	3.1
A2	MHD, disk wind off	10−7	0	150	0	104	0.9
B2	HD, disk wind off	10−7	0	150	0	104	3.0
A3	Cool MHD, both winds on	10−7	10−7	150	50	102	1.3
B3	Cool HD, both winds on	10−7	10−7	150	50	102	3.9
A4	MHD, light stellar wind	10−8	10−7	150	50	104	2.6
A5	MHD, light stellar wind, eq. vel.'s	10−8	10−7	100	100	104	1.7
A6	MHD, light stellar wind, eq. vel.'s	10−8	10−7	100	100	102	1.6

Download table as:  ASCIITypeset image

The domain of computation is scaled so that one computational unit of length corresponds to $500 \rm \ AU$. The physical size of the domain thus corresponds to a length along z of $2\times 10^4\rm \ AU$ and a length along r of $8\times 10^3\rm \ AU$. The simulations are each run for the length of time required for the flow to reach the right end of the domain. These times are given in physical units in Table 2.

We note that the use of a constant density ambient medium was a response to our desire to make the simulations as general as possible, allowing these models to address issues relevant to YSO jets, PNs, and other disk/central source systems. Not including a 1/r² density fall off appropriate to the AGB wind will affect the results in terms of timescales (as the winds/jets will see lower momentum densities in the ambient medium as they propagate outward). If the AGB wind is spherically symmetric on nebular scales, we do not, however, expect dramatic changes in morphology. This point can be explored in further studies that are beyond the scope of the current work. We note that some aspects of this problem have been covered by Akashi (2007) and Akashi & Soker (2008) for the nonmagnetic case.

We also note that the velocity scales are chosen to be appropriate for either YSOs (a wind from the inner edge of a disk) or the preplanetary nebular phase. For the preplanetary case, we take this speed to be indicative of the bridge between the AGB wind ($10 \rm \ km\ s^{-1}$) and the circumstellar wind ($1000 \rm \ km\ s^{-1}$). In addition, we note that we are interested in the long term of evolution of the morphology and so we have had to make certain choices based on computational expediency and our desire to provide simulations that are generally relevant to nested-wind systems. Thus, the size scale of the disk boundary condition is larger than should be expected in some systems. Future studies will be needed to connect behavior at the smallest scales where the disk winds are launched and disk and stellar winds interact (Garcia-Arrendondo & Frank 2006) and the largest scales where full nebular morphology has been established.

While our initial ambient temperature is appropriate for mature PNs (see Akashi & Soker 2008), it is too high for PPNs and YSOs except irradiated YSO jets where T = 10⁴ K is a reasonable choice for the ambient medium. As we shall demonstrate, however, in the case of MHD winds, which are our principal concern, the choice of the ambient temperature is not a significant factor in determining the morphology. We have included a discussion early in the paper on the distinction between YSO and PN temperatures.

We use relatively wide jets (θ>10°) in our simulations. The presence of such wide jets or wide-angle winds has been conjectured for sometime in YSO systems (Shu et al. 1994) and can be seen as a diagnostic for launch mechanisms. In PN systems, Soker (2004) has presented analytical arguments for the existence of such wide winds/jets. Such wide outflow systems are likely to be important for creating wider lobes in the observed nebulae in both classes of bipolar nebulae (YSOs and PNs).

In all cases, optically thin, atomic line cooling based upon a cooling curve is assumed. A temperature "floor," T_f = 9.0 × 10³, is set so that only material heated to temperatures T ⩾ T_f is subject to radiative energy loss. No attempt has been made to follow the ionization dynamics of the flow. Note that adiabatic cooling continues to operate below the floor for radiative cooling. Because our focus in this paper is restricted to the morphological features that arise from the interaction of nested winds and from their interaction with their common environment, this approximation is not expected to materially affect our conclusions.

We present the results of our simulations in Figures 2–7. In Figure 2, we present late-time density maps for runs A1, A2, B1, and B2.⁶ The intent here is to contrast the case of simultaneous nested winds with the case of the stellar wind in the absence of the disk wind, for both the MHD (runs A1 and A2) and HD (runs B1 and B2) cases. Collimation is evident in all four simulations. For the HD cases, this is a consequence of the ram pressure of the ambient medium. For the MHD flows, the hoop stresses associated with the toroidal field will also contribute.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The presence of the slower disk wind does not appear to have as significant an effect on the resulting structure of the flow in the case of the HD winds. In both cases, the stellar wind is refocused toward the axis via the shock at the wind/wind or wind/medium interface over comparable length scales with the focusing length slightly smaller for the case of the single-wind simulation. This small reduction is likely due to the fact that mixing of the disk wind with the post-shock stellar wind will impart additional z-directed momentum flux to the flow. Otherwise, the resulting shock structures appear very similar in the HD cases.

The differences between the nested-wind and single-wind cases for the MHD cases, on the other hand, are quite striking. We note first that in the absence of the disk wind, a large rarefied region (a cocoon) opens up along the axis of the wind, while a relatively short MHD nose cone evolves at the working surface of the wind shock with several vortices appearing along the bow shock. When both winds operate, however, there is considerable mixing between the two flows giving rise to the filimentary structures appearing in the region near the axis in place of rarefaction. In addition, the bow shock has developed a "shoulder," below which an extended and somewhat flattened nose cone protrudes with what appears to be a refocusing event similar to those seen in the HD cases occurring in its interior. The shoulder begins to develop early in the simulation soon after the formation of the nose cone and is built up smoothly over the course of the simulation as separate shedding events occurring near the head of the nose cone combine, cool, and expand.

Next, we compare maps of the total energy density with the magnetic energy density and with maps of plasma-β for both the nested-wind and single-wind cases. Results for the nested-wind case are shown in Figure 3, while those of the single-wind case are presented in Figure 4. The top two panels of Figure 3 compare total energy density and magnetic energy density for the nested-wind case. We first note that the map of the total energy density is quite similar in appearance to the map of the density. This is to be expected since most of the energy in the system is either thermal or kinetic and thus closely traces the density. The map of the magnetic energy density has a different appearance. Here we note that in the region bounded on the right by the bow shock "shoulder," the magnetic energy density is distributed in two roughly uniform layers—distinguished by their typical values—throughout the area enclosed and falls off steeply only near the wind/medium interface. The energy densities characteristic of the inner layer are evidently of order $E_{\mathrm B}\gtrsim 10^{-6}{\rm \ erg\rm \ cm^{-3}}$ while those of the outer layer are characterized (very roughly) by energy densities in the range $0{\rm \ erg\rm \ cm^{-3}}\lesssim E_{\mathrm B}\lesssim 10^{-10}{\rm \ erg\rm \ cm^{-3}}$. It is interesting to note that in the extended conical region (an MHD nose cone) to the right of the shoulder, the typical magnetic energy density is of the order of the highest values found in the outer of the two layers to the left of the shoulder. This result is confirmed in the maps of β given in the bottom two panels of Figure 3 which indicate values of β ≲ 1.0 throughout much of the inner layer, while in both the outer layer to the left of the shoulder and the conical region beyond, we find β ≳ 10. Note that in the maps of β, the color black is serving double duty and indicates very low values of β in the interior of the flow and "infinite" values exterior to it. For this reason, we have presented two maps for each simulation: one in which the full range of values of β are mapped and another "overexposed" map in which values of β are restricted to the range 0 ⩽ β ⩽ 1. In this way we are able to show that the very dark interior regions in the panel with unrestricted range are regions of low-β plasma. Note that the magnetic field in the cocoon/bubble is not, in fact, very strong and does not play an important role in determining the morphology there. This can be seen by examining the plots of β. In the body of the jet downstream of the interaction region we see that β can be of order 1, and in those regions the field is contributing more significantly to the dynamics.

The noticeable decrease in magnetic energy density in the extended conical region is particularly interesting in light of the observation made earlier of what appears to be the occurrence of a HD refocusing event in the interior of the extended conical region. These results suggest that while the inertia and overpressuring of the winds initially lead to an expansion into the ambient medium of the winds, hoop stresses—coupled with the gradual accumulation of shed material—eventually create a bottleneck for the more highly magnetized material in the interior, leading to the ejection of relatively less magnetized material, which then exhibits refocusing similar to the HD wind within this extended region beyond the shoulder.

The single-wind case shown in Figure 4 tells a very different story. While the map of total energy density shows some layering indicative of larger energy densities in regions with more material, the magnetic energy is distributed uniformly throughout the region of the flow including in the nose cone at the tip. This is again confirmed in the bottom two panels of Figure 4 showing the maps of β for this case.

Most of the simulations run for this paper assumed a warm ($T_a=10^4 \rm \ K)$ ambient medium. Since the type of flows studied here are hypothesized to occur in the nonionized mediums of post-AGB stars, we present in Figure 5 the density maps for a pair of nested-wind simulations for the cases of MHD and HD flow, respectively, with ambient temperature set at $T=100\rm \ K$. The results shown are from runs A3 and B3 in Table 1, respectively. We note that in the case of the MHD nested wind, the same bow shock shoulder, filamentary structures, and extended and flattened nose cone (with some refocusing again evident) are all present in run A3 as well. A comparison with the top panel of Figure 2 also indicates that in spite of the now overpressured jet impinging upon an ambient environment that has had its pressure lowered by two orders of magnitude, the flow suffers very little additional expansion into the medium, suggesting that the hoop stresses play the predominant role in maintaining collimation for the case of MHD nested winds. The HD nested wind, on the other hand, suffers a significant amount of expansion when the pressure of the ambient environment is reduced. The rarified region along the axis quadruples its radius and is now bordered by a relatively dense "beam" of material consisting of thoroughly mixed post-shock stellar and disk material, with a characteristic width much greater than the original dimensions of the flow. This layer is in turn bordered by an even denser layer of ambient material plowed up by the flow which, in spite of the shock that forms at its outer edge, fails to warm sufficiently to give rise to significant radiative cooling. The resulting appearance is that of an adiabatic wind "piggybacked" upon a radiative wind. This impression is further suggested by the presence of a nose cone near the axis and a Mach disk above and to its left.

It is expected for both the post-AGB and YSO systems that both the stellar and disk winds are expected to have—to the order of magnitude—comparable velocities and mass-loss rates. Still, it is important to investigate how differences in these quantities within the expected range of variation are likely to affect the appearance and physics of the flows. With this in mind we have produced the three simulations for which we present late-time density maps in Figure 6. In the top panel of this figure we show the results of run A4, in which the velocities of the flows—$150 \rm \ km\ s^{-1}$ for the stellar wind and $50 \rm \ km\ s^{-1}$ for the disk wind—are unchanged from our standard run A1, but the mass-loss rate for the stellar wind has been reduced by one order of magnitude relative to the mass-loss rate for the disk wind. Because the temperature of the stellar wind is determined by its velocity and the two parameters σ and β_m, all of which are unchanged from run A1, the effect is to lower the mean density and therefore the pressure of the stellar wind. This allows material from the slower disk wind to expand into, mix with, and disrupt the flow of the stellar wind material, creating a turbulent, filamentary field in the region near the axis which, in run A1, is relatively evacuated.

In the second panel of Figure 6 we show results from run A5, in which we have now set the velocities to equal values ($100 \rm \ km\ s^{-1}$). We see that because of the reduced pressure of the stellar wind, the disk wind is still able to expand into the region of the stellar wind and mix with the material there, but because the velocities are equal, this expansion is less disruptive to the flow there. We see a relatively more uniform density field in this region. From comparing these panels with the first panel of Figure 2, in which we also noted the presence of filamentary density structures, we may conclude that differences in both mass-loss rates and flow velocities contribute to determining the degree to which these filamentary density structures arise, and to the characteristic length scales associated with them, in the region of the flow near the axis. In the bottom panel of Figure 6, we present the result of run A6, which is identical to run A5 with the exception that the temperature of the ambient medium has been reduced from $10^4 \rm \ K$ to $100\rm \ K$. It is evident in this panel that the effect of this reduction is to allow for greater expansion of the entire overpressured flow into the surrounding medium with otherwise little qualitative alteration in the appearance of the flow.

In an attempt to approximate roughly how our interacting winds might appear on the sky, we present in Figure 7 synthetic maps of emission for both our nested-wind and single-wind cases as obtained in the MHD and HD runs. The intensity shown, which does not distinguish among cooling lines, was determined according to

$$\begin{equation} I_{i,j,k} = \Sigma _k n_{i,j,k}^2\Lambda (T_{i,j,k}), \end{equation} \tag{ 10 }$$

where i, j, and k refer to the x, y, and z directions, respectively, in the final data cube created by rotating n(r, z) and T(r, z) about the axis of symmetry, and Λ is the cooling function. A projection angle of 20° is assumed in all panels. The left column shows the results for the MHD runs with the nested winds shown in the top panel and the single wind shown below. The column on the right shows the corresponding results for the HD case. The maps were generated by revolving the 2.5D data about the axis of symmetry, "tilting" the symmetry axis by a projection angle of 20°, and subsequently summing along lines of sight and projecting the results on to the image plane. The results are for the total radiative emission. No attempt was made to distinguish among lines of emission. The images show a marked distinction between the MHD and HD cases. While the HD simulations give rise to smooth conical segments of emission that broaden with distance, the MHD results instead give rise to emission divided by primarily on-axis features and ring-like structures centered on the symmetry axis. We also note an interesting difference that appears in the on-axis emission of the nested winds as contrasted with the stellar-wind case. The emission features on both axes are both quite narrow; but while the on-axis emission for the stellar wind is quite smooth, the corresponding MHD emission shows small but distinct knots of emission distributed along the left half of the axis and vanishing approximately where we earlier noted a bottleneck in magnetic energy arose. While it is possible that the smooth emission seen in both panels could be, in part, an artifact of the numerical method employed for handling the cylindrical symmetry of the problem, the appearance of the knots in the nested-wind case and their absence in the stellar-wind case suggests that these features are "real." Given that well-aligned knots of emission are routinely observed in optical HH jets, these results, while preliminary, provide cause for asking if these knots of emission are a consequence of interacting nested winds in the environments of YSOs.

The feature that most distinguishes the MHD nested wind from the other scenarios considered here is the development of the bow shock shoulder mentioned above. This feature was plainly evident in both the warm and cool nested-wind MHD simulations presented here, and would serve as a marker for determining whether an observed object has been formed from simultaneously operating magnetized disk and stellar winds. Before we can robustly employ such a feature as a means of identifying objects as candidates for formation by nested winds, it will be necessary to study the formation of the shoulder in a greater detail. For example, it will be necessary to determine how the formation and appearance of the shoulder varies with such parameters as the disk opening angle, ratios of disk and stellar flow velocities, and ratios of mass-loss rates. (In particular, it is evident from Figure 6 that the shoulder is much more subtle in appearance when the mass-loss rate of the star is significantly lower than that of the disk.) But to illustrate what we have in mind, we present as a case in point the HST image of the planetary nebula Hen 2–320 in Figure 8. Though other interpretations are possible, the right lobe of this object exhibits features that are qualitatively similar to our MHD nested-wind simulations. We see—indicated by the arrows labeled "shoulders"—an expanded region of the flow nearer to the nebular core and—indicated by the arrows labeled "nose cones"—a narrower and conical extended region to its right. We also note a brightening of emission at the location where the expanded region gives way to the narrow conical region. This brightening is reminiscent of the large ring of (synthetic) emission seen in the upper-left panel of Figure 7, which, we note, coincides with the location of the bow shock shoulder that arises in this simulation.

Zoom In Zoom Out Reset image size