RADIATION-HYDRODYNAMIC SIMULATIONS OF PROTOSTELLAR OUTFLOWS: SYNTHETIC OBSERVATIONS AND DATA COMPARISONS

Identifying the gas associated with outflowing material remains a key observational challenge. Generally, outflow gas is higher velocity, hotter, and less dense than the accreting material. It is difficult to use temperature as a means to distinguish outflow gas. In nearby well-resolved regions, it is possible to image outflow cavities in scattered light (Seale & Looney 2008). However, observers typically use low-density tracers such as ¹²CO to map the outflow gas, since these data supply both direct velocity information and indirect estimates of the gas mass. Within molecular line datacubes, line-of-sight velocity remains the primary means of identification. In AS06, outflows are defined by selectively integrating over the ¹²CO(1–0) emission in a range of high-velocity channels determined by visual inspection.

To identify the outflows in the simulations, we set a minimum outflow gas velocity of 2 km s⁻¹. We then generate a column density map of cells above this cutoff, which is analogous to an observed intensity map integrated over selected velocity channels.

Figure 2 shows the angle distributions of integrated emission for schematic elliptical outflows. In this two-dimensional space, we treat each outflow lobe independently and characterize outflow morphology using five parameters: opening angle, length, width, shape, and inclination. To measure the opening angle, we calculate the angle with respect to the z-axis for each pixel that is above a given velocity minimum. In the analysis, we weight all such pixels equally. The opening angle is then defined to be the full width at quarter-maximum of the angle distribution. This definition is empirically selected to correspond to opening angles derived from intensity maps by eye (and protractor). We derive the outflow inclination with respect to the vertical axis by fitting the angle distribution with either a Gaussian or higher order polynomial. Note that an outflow may also be inclined along the line of sight toward (or away from) the observer. However, our method only fits for the inclination in the plane of the sky. For an elliptical outflow, the inclination is the location of the distribution maximum, while for a more conically shaped outflow, it is the local minimum. Note that since the simulation inclination is obtained from the projected distribution of cells it may not identically correspond to the direction in which the outflow is launched. In addition, interaction between the outflow and the envelope may affect the inferred outflow inclination.

Zoom In Zoom Out Reset image size

For a given distribution of connected high-velocity cells, the physical length of the outflow is the longest extent measured from the center outward along the outflow major axis. The width is the maximum extent along the axis perpendicular to the outflow length. The Gaussian parameter, Gauss, is the χ² value of a Gaussian fit to the angle distribution. As illustrated by Figure 2, the distribution of pixels for elliptically shaped outflows is well described by a Gaussian. Our opening angle definition is formulated to be consistent with an elliptical outflow such that arctan(Width/Length) ≃ θ. As we show in Section 4.2, the pixels close to the protostar are likely to be dominated by the beam characteristics, and thus are not suitable for an opening angle determination.

We first characterize an idealized, symmetric outflow geometry. Figure 2 shows an ellipse of length 0.1 pc, which has been superimposed upon the AMR grid cell hierarchy of a simulation output. The effect of truncating the outflow length is apparent in the angle distribution evolution shown on the right. The inferred angle increases with increasing truncation as the wide-angle cells at the base of the ellipse dominate the distribution.

Figure 3 shows the opening angle and major axis inclination for a series of truncated elliptical outflows. The error bars indicate the $\sqrt{N}$ error, where N is the number of beams contained in the ellipse assuming that the outflow is placed 250 pc away and observed with a beam of 4''. Figure 3(a) illustrates that the error bars do not necessarily contain the inferred opening angle for the outflow with full information (i.e., no truncation) even though the uncertainty increases as the outflow length shrinks. This indicates that a poorly resolved outflow may be measured to have a fundamentally different opening angle than the same outflow with complete information. The outflow major axis remains fixed during the truncation. Figure 3(b) demonstrates that the derived direction of the outflow axis is fairly insensitive to the truncation.

Zoom In Zoom Out Reset image size

For comparison, we re-analyze the AS06 data to obtain the opening angles quantitatively. Figure 4 shows the measured angles for the upper and lower lobes and the mean of the two. The slope of the fit we find to the data, 0.15, is quite close to the value of 0.16 ± 0.4 found by AS06. The intercept of 1.22 is also within error of the previously published value of 1.1 ± 0.2. This is encouraging since it suggests that our method is a good quantitative alternative to fitting the outflow angles by eye. One caveat to the fitting method is that for small numbers of pixels the algorithm becomes sensitive to the bin size. We indicate this uncertainty using vertical error bars on the plot, which are twice the bin size used in the fit.

Zoom In Zoom Out Reset image size

As shown in the bottom panel of Figure 4, the angle distribution tends to evolve from an elliptical shape (Gauss (χ²) ⩽20) to a more conical or irregular shape (Gauss (χ²) >20). This trend may be partially an artifact of how the outflow is observationally sampled at late times rather than an indication of how the geometry actually changes. This is supported by the apparent decrease of outflow length with time, which is most likely because the older gas has blended with the ambient gas.

Figure 5 shows a volume rendering of the simulation velocities for both runs. The outflow axis in R1 is almost directly aligned with the z-axis so that v_z traces the fastest moving gas, which ranges up to 70 km s⁻¹. R2 forms two protostars with fairly distinct outflows. The secondary protostar arises from fragmentation in the core rather than disk fragmentation, so that the outflow axes of the two are somewhat misaligned. However, the two R2 protostars are sufficiently close that separating the outflow lobes in projection is difficult. Consequently, in the following analysis we will use R1 to explore outflow evolution for individual sources, while R2 will illustrate a case in which observations are confused by the presence of a second, unidentified outflow source.

Zoom In Zoom Out Reset image size

All three outflows shown in Figure 5 exhibit a significant amount of asymmetry between the lobes. The sub-grid model launches the outflow gas from the grid symmetrically about the protostars. (On the scale of the figure the outflow launching region is contained in the central pixel.) This means that any asymmetry that arises is directly the result of the interaction between the outflow and the turbulent envelope. For example, the R1 protostar forms in a more filamentary structure and the net angular momentum direction is almost directly aligned with the filament axes. As a result, the positive z outflow lobe is diverted from the launching axis and confined by the dense filament gas. The outflow lobe along the negative z-axis breaks out of the filament and extends further. In R2 asymmetry in the core also leads to mismatched outflow lobe sizes.

Figure 6 shows the R1 outflow evolution as a function of time. The outflow gas is identified using a velocity cutoff of 2 km s⁻¹, where the full three-dimensional information is used to determine the outflow and envelope masses. These cells are projected along the x-direction so that the outflow opening angle, width, and length are calculated when the outflow axis is nearly parallel to the z-axis. A fit to the opening angles demonstrates that the lobes widen with time similarly to observations. The earliest time only contains a small number of cells in the outflow so that it is artificially broadened, but a strong trend is apparent for 3.7 ⩽ log (t) ⩽ 4.7. (See Section 5 for discussion of the origin of this broadening.) However, the different outflow lobes in Figure 6 differ from one another and display significant shorter-time variation. For example, the upper outflow lobe growth stalls at late times. This suggests that even if outflows generally evolve with time, variation in individual outflow behavior can be large.

Zoom In Zoom Out Reset image size

Considering the large uncertainties in observed outflow ages, poorly constrained line-of-sight inclinations, and instrument limitations, the agreement between the simulations and the observations shown in Figure 4 is striking. This supports the finding of AS06 that outflow angles evolve with time. However, the protostellar masses in these calculations are ${\lesssim} 0.1 {\,M}_{\mathord \odot }$ for the time of comparison, so they are still very young. Outflow sizes are ≲ 0.1 pc, although this is similar to the sizes of those in AS06, which are ∼0.005–0.1 pc. Figure 6 illustrates the length of only the connected cells with velocities ⩾2 km s⁻¹. These cells are also very low density and have a total mass of ${\sim} 0.001 {\,M}_{\mathord \odot }$.

The R1 opening angle extrapolated at t = 0 is lower than that found by AS06, although the positive z fit is still within the observational error. Smearing due to the observational beam likely accounts for some of this discrepancy. Figure 7 shows the measured opening angles for R1 as a function of convolution with different beam sizes. For a 4'' beam, which is comparable to the resolution of AS06, some angles may be broadened by up to ∼80%. Broadening becomes more significant for larger beam sizes. The earliest times in the simulation, for which the outflow is already broadened by the small number of cells, appear least affected by the beam smearing.

Zoom In Zoom Out Reset image size

While comparison between the simulation and observed ¹²CO data in Section 4.2 is suggestive, it is important to compare the data directly with synthetic observations in ¹²CO. In this section we use the radiative transfer code MOLLIE (Keto & Rybicki 2010) to model the emission from the first seven ¹²CO and ¹³CO transitions. We assume a standard abundance of 5.6 × 10^{−5 12}CO relative to H₂ and an abundance of 7.3 × 10^{−5 13}CO relative to H₂, which are comparable to CO estimated abundances in low-mass envelopes and outflows (Carolan et al. 2008, 2009). However, there are a wide range of inferred CO abundances among different cores and star-forming regions, values of which are often individually uncertain by a factor of two or more. These uncertainties complicate the estimation of gas mass beyond the already challenging problem of identifying the outflow gas and accounting for projection.

We perform the radiative transfer calculation on a cube of side 0.15 pc centered on the protostar. We regrid the AMR data to produce a grid of 128³ with 130 AU cell resolution nested within a grid of 128³ and 260 AU cell resolution.

Figure 8 shows integrated ¹²CO(1–0) emission maps from observations of R1 at different inclinations with respect to the line of sight. This includes all gas with velocities in the range ±20 km s⁻¹, which encompasses the minimum and maximum outflow velocity at 130 AU resolution.⁴ The outflow structure is clearly visible in the integrated maps. Such structure may not be apparent in observed integrated maps due to emission from ambient cloud gas between the source and the observer. The right panels show the map an observer would make when integrating over gas in higher velocity channels (|v| ⩾ 2 km s⁻¹).

Zoom In Zoom Out Reset image size

The lower outflow lobe shows substructure that could be interpreted as the result of accretion bursts. However, the accretion history, and hence the outflow mass ejection rate, is fairly smooth (see Figure 12). Instead, the discontinuous morphology results from the outflow interacting with the turbulent and asymmetric core gas.

Figure 9 shows the opening angle versus inclination measured from the ¹²CO emission maps with no beam smearing. The opening angle appears slightly broader as the outflow axis becomes more parallel to the line of sight, but there is not a strong dependence on the inclination. At low inclinations, the highest velocities are mainly perpendicular to the line of sight and the angle is measured using only a small sample of cells. When the outflow motion is mostly along the line of sight, the lobes become shorter and rounder. Observationally, outflow orientations with completely overlapping lobes are not well suited to measuring opening angles compared to those with intermediate inclinations. Although the outflow inclination is difficult to infer from projection, the axes of the outflows investigated by AS06 are likely inclined between 30° and 60° with respect to the line of sight.

Zoom In Zoom Out Reset image size

Figure 10 shows the opening angle versus time using the CO intensity maps to calculate the angle. Like Figure 6, the angle increases on average as a function of time. Linear fits to the average opening angles in Figure 10 give θ = 0.80 + 0.18log (t) and θ = 0.53 + 0.25log (t) for inclinations of 30° and 45°, respectively. However, like Figure 6, the two lobes show different trends with time. For the projected data, the lower outflow lobe is nearly flat as a function of time, while the upper lobe widens significantly. This is mainly due to the concurrent widening and elongation of the lower lobe. Despite the variation, this suggests that trends inferred from molecular ppv information are reflective of the actual ppp information. In this case, the similarity is particularly strong because CO is an efficient tracer of the low-density gas within the outflow cavity.

Zoom In Zoom Out Reset image size

Molecular hydrogen has no dipole moment and is observationally invisible unless the gas is strongly shocked (Yu et al. 1999). Consequently, CO, which is the next most abundant molecular species, serves as the primary tracer of molecular gas. In this section we consider two methods for estimating molecular gas mass from CO emission.

For an assumed abundance of CO, a conversion factor can be used to relate CO intensity to total gas mass. This quantity, known as the X-factor, is defined as

$$\begin{equation} X_{\rm CO} = \frac{N({\rm H_2})}{W(^{12}{\rm CO})}, \end{equation} \tag{ 2 }$$

where W(¹²CO) is the ¹²CO(1–0) integrated intensity and N(H₂) is the column density of molecular hydrogen. The abundance of CO depends upon the local gas density, temperature, and UV field, so that in reality, the X-factor varies widely over an individual cloud (Pineda et al. 2008; Glover & Mac Low 2011) However, this complicated chemistry is usually reduced to an average X-factor. As discussed in Section 4.3, we adopt a typical mean CO abundance from observations. Here, we use the simulated outflows to assess how accurately the standard X-factor and other more complicated methods recover outflow mass.

We adopt X_CO = 1.8 × 10²⁰ cm⁻² K⁻¹ km s⁻¹, which is the mean value for molecular clouds in the solar neighborhood (Dame et al. 2001). Pineda et al. (2008) found that this factor overestimated the gas mass by a factor of 45% compared to the estimated extinction mass. However, the typical mean abundances inferred by Pineda et al. (2008) are a factor of ∼10–30 higher than the value we assume for our molecular line transfer calculation, where the difference is due to lower expected abundances for dense core and outflow gas versus lower density cloud gas (e.g., Carolan et al. 2008).

Pineda et al. (2008) found that the X-factor is most reliable when the gas is both optically thin and A_v > 4. While the emission from the outflow gas alone is marginally optically thin (with increasing optical depth as the inclination approaches 90° and lower velocities), the integrated ¹²CO(1–0) emission through the core envelope is optically thick with τ > 10. This means that the emission is saturated and X_CO, which assumes a linear relationship between emission and gas mass, will underestimate the outflow gas mass.

Bally et al. (1999) and Yu et al. (1999) present a more nuanced approach for obtaining the outflow gas mass from CO emission. Their methods are similar to that of Arce & Goodman (2001, henceforth AG01) which we follow here. All three techniques improve upon the estimation above by correcting for optical depth effects and including more of the low-velocity outflow gas mass. We outline the AG01 procedure below and post-process the simulations accordingly (see AG01 and references therein).

While the ¹²CO(1–0) emission may be quite optically thick, it is possible to utilize the more optically thin ¹³CO(1–0) emission to calculate the true optical depth and correct the ¹²CO(1–0) mass estimate. In the first step of the procedure, we spatially average over all the ¹²CO(1–0) and ¹³CO(1–0) spectra and derive the ratio of the two lines, R_12/13, as a function of velocity. The resultant parabolic line ratio shape indicates that the opacity is not constant with velocity as many studies assume.

We then perform a polynomial fit of R_12/13(v), which is constrained to reach a minimum at the cloud velocity. Following AG01, who limit the velocity range to exclude a secondary coincident cloud, we limit the fitted velocity range to within 1.5 km s⁻¹ of the minimum. Using this fit, we extrapolate to the high-velocity wings, where observationally ¹³CO is too weak to be detected. (Even in the noiseless synthetic observations, the ¹³CO emission is negligible at velocities above a few km s⁻¹.)

We next derive an effective "¹³CO main beam temperature," T¹³_mb(x, y, v), to estimate the ¹³CO opacity. Observationally, if the ¹³CO(1–0) emission is greater than twice the rms noise, then T¹³_mb(x, y, v) can be used directly. Otherwise, the ¹³CO main beam temperature must be inferred from ¹²CO: T¹³_mb(x, y, v) = T_mb¹²(x, y, v)/R_12/13(v_i). We adopt an rms noise value of 0.06 K per 0.2 km s⁻¹ channel, similar to that used by AG01.

If ¹³CO(1–0) is indeed optically thin, then its opacity may be estimated as follows:

$$\begin{equation} \tau _{13}(x,y,v) = - {\rm ln} \left(1- \frac{{T_{\rm mb}^{13}(x,y,v)}}{T_0/[{\rm exp}(T_0/T_{\rm ex})-1]-0.87} \right), \end{equation} \tag{ 3 }$$

where T₀ = hν/k = 5.29 and T_ex is the excitation temperature under the assumption that ¹²CO(1–0) is optically thick:

$$\begin{equation} T_{\rm ex} = \frac{5.53}{\rm {ln}[1+5.53/{(\it T}_{\rm peak}+0.82)]}. \end{equation} \tag{ 4 }$$

Here, T_peak is the peak temperature for the ¹²CO emission. We find T_ex ∼ 10–12 K, in good agreement with the known simulation ambient gas temperature of 10 K.

Finally, we derive the column density of ¹³CO for each pixel:

$$\begin{equation} N_{13}(x,y,v) = 2.42 \times 10^{14}(T_{\rm ex}\, +\, 0.88) \frac{\tau _{13}(x,y,v) dv}{1-{\rm exp}(-T_0/T_{\rm ex})}, \end{equation} \tag{ 5 }$$

where dv is the channel width in km s⁻¹.⁵ The total mass per pixel of area A is given by

$$\begin{equation} M(x,y,v) = m_{\rm H_2}N_{\rm H_2}(x,y,v)A, \end{equation} \tag{ 7 }$$

where $m_{\rm H_2} = 2.72 m_{\rm H}$ is the mean molecular weight and the conversion between the ¹³CO column density and the molecular hydrogen column density is assumed to be $N_{\rm H_2} = 1.4 \times 10^6 N_{13}$. Note that the exact value depends upon the local X-factor. When modeling the HH300 outflow in the Taurus star-forming region, AG01 adopt $N_{\rm H_2} = 7 \times 10^5 N_{13}$, a generic value for dense gas in Taurus that was obtained by Frerking et al. (1982).

An additional step is necessary to distinguish between cloud mass and outflow mass. Summing over the total area, we obtain M(v) as shown in Figure 11. AG01 found that the peak, which includes the bulk of the core mass, was well fitted by a Gaussian. We find that a Gaussian provides a satisfactory fit in the case of R1, where there is only a single outflow, while the R2 data, in which there is a second "hidden" outflow, a Gaussian does not well describe all of the high-density gas (see Figure 11). We then subtract this fit from M(v) and integrate over all velocities to get the outflow mass, i.e., the shaded area shown in Figure 11.

Zoom In Zoom Out Reset image size

Figure 12 shows the simulated outflow mass compared with that derived using the X-factor and the AG01 method. For R2, Figure 12 displays the total mass ejected by both stars. Although the axes of the two R2 outflows are not aligned, projection, beam resolution, close proximity, and the low velocities of the outflows would prohibit individual identification of the two outflows observationally.

Zoom In Zoom Out Reset image size

In the figure, the "total ejected mass" is the total mass launched by the protostar according to the outflow model. This is an essentially numerical quantity representing the total mass deposited in the zones around the protostar as described in Section 2. This quantity is distinct from what an observer would identify as the "outflow mass," which is defined on the basis of velocity and thus also includes entrained gas. At any given time, the launched mass includes gas that goes on to comprise a low-velocity, less collimated component of the outflow and gas that mixes with the ambient gas becoming indistinguishable even with full information.⁶ Since the launched material is placed on the grid near the protostar by design, this total does not include entrained envelope gas. Once the wind is deposited in the grid it is impossible to differentiate launched gas from entrained, accelerated gas. In principle, due to the addition of entrained gas, the total mass of high-velocity gas, i.e., the outflow mass, could exceed the mass numerically ejected by the protostars.

For both runs, the mass derived from the raw simulation data with |v| ⩾ 2 km s⁻¹ and outflow orientation of 45° is significantly less than the total ejected mass by a factor of 5–10. This reinforces the point that the majority of outflow gas in an inclined outflow is simply not observable even with a perfect method for converting between CO emission and total gas mass. However, the integrated ejected outflow mass does not include gas that is entrained and accelerated by the outflowing material, so in principle the simulated and observed gas could be larger than the total ejected mass.

The gas mass estimated from the CO emission using the X-factor appears to track the mass of the high-velocity gas fairly well, although it is generally a factor of 5–10 lower. This discrepancy is due to the high mean optical depth of the core (τ > 10 for ¹²CO(1–0) and τ ∼ a few for ¹³CO(1–0)). The AG01 method, which subtracts off the core mass, appears to perform somewhat better for most simulation outputs. However, the accuracy of the mass estimation strongly depends upon how well M(v) can be fitted by a simple function. For example, at late times the R1 M(v) is slightly asymmetric, leading to a negative derived outflow mass. Note that Figure 11 shows the log of the mass, whereas the fit is dominated by a few points around v ≃ 0 km s⁻¹. Differences between the gas mass and fit at these low speeds far exceed the total mass contributed by high-velocity gas in the distribution wings. To account for this, we only subtract the curves beyond the full width at half-maximum of the Gaussian fit and omit velocities where the difference is negative. However, since the mass in the core is so much higher than the mass in the high-velocity channels, the mass estimates strongly reflect the symmetry and Gaussianity of the mass distribution. The success of AG01 in estimating the mass of outflow HH300 relies upon the goodness of fit of a Gaussian and the relatively larger outflow mass relative to the core mass.

If the outflow mass is instead defined as a sum over M(v) where |v| ⩾ 2 km s⁻¹ then the mass estimate follows a similar trend to the mass derived using the X-factor. However, this suggests that the more complicated AG01 procedure still underestimates the total outflow mass in the optically thick case. The factor of ∼5 discrepancy between the two methods is due to different assumptions about CO abundance.

Observationally, ¹²CO(1–0) and ¹²CO(2–1) may be used in combination to derive the excitation temperature, T_ex. Temperature variation may be used to discriminate between different outflow models (Arce & Goodman 2002). For example, shocked outflow gas should be higher temperature than ambient gas, and the gas temperature should rise with the maximum outflow velocity (Lee et al. 2001). For optically thin gas, the excitation temperature may be related to the ratio of the line intensities by

$$\begin{equation} R_{21/10} =4 e^{-11/T_{\rm ex}}, \end{equation} \tag{ 8 }$$

where R_21/10 is the ratio of the ¹²CO(2–1) to ¹²CO(1–0) lines (Arce & Goodman 2002). Observers typically assume that the high-velocity gas that constitutes the outflows is optically thin. We find that in many cases line ratios in the channels with |v| ⩾ 2 km s⁻¹ are above 1.0, the line ratio limit in the optically thick case, indicating that the gas is at most marginally optically thick in the outflow.

Figure 13 shows that the estimated excitation temperatures are generally higher than the ambient gas temperature of 10 K. The excitation temperatures also increase with time, which is expected if the shock strength and, hence, post-shock temperature increase with time. This effect also occurs in the simulations, although a rising gas temperature also results from increased protostellar heating (see OKMK09). The two causes are difficult to distinguish from a single average temperature value. However, T_ex does not appear to be well correlated with the actual gas temperatures in the simulation. The disagreement at later times is likely due to the increasing optical depth of the high-velocity gas. Hatchell et al. (1999) show that gas with τ = 1 and inferred T_ex = 30 K under the assumption that the gas is optically thin will actually have T_ex = 40 K. The optically thin approximation increasingly underestimates temperatures for higher optical depths and gas temperatures.

Zoom In Zoom Out Reset image size