Atomic Shocks in the Outflow of L1551 IRS 5 Identified with SOFIA-upGREAT Observations of [O i]

To combine two [O i] spectra taken with two setups of local oscillator frequencies ("m50" and "p50;" see Section 2), we took the weighted average of the two spectra using on-source time as the weight. Figure 2 shows the [O i] spectra observed toward the "red," "center," and "blue" positions. For the baseline fitting and subtraction, we selected velocity ranges of −150 km s⁻¹ < v_lsr < −110 km s⁻¹ and 20 km s⁻¹ < v_lsr <130 km s⁻¹. The baseline noise is 0.10 K, 0.15 K, and 0.09 K for the "red," "center, " and "blue" positions, respectively. The "center" spectrum shows a clear detection of the [O i] 63 μm line, with a peak at around 0–5 km s⁻¹ and a broad blueshifted tail. We detect no emission at redshifted velocities at the "center" position. We also do not detect any significant [O i] emission in the "red" or "blue" positions in the outflows. To better characterize the broad blueshifted feature at the center position, we further smoothed the [O i] spectra to 5 km s⁻¹ resolution. At the "center" position, the smoothed [O i] spectrum shows two tentative local peaks at −53 and −31 km s⁻¹ in addition to the prominent peak around the source velocity. The intensities of these two peaks are roughly equal to the half of the maximum intensity. If these two peaks belong to the same broad feature, the FWHM of the [O i] line is then ∼100 km s⁻¹. However, if these two peaks are separate features, the primary line profile centered at the source velocity would have an FWHM of ∼40 km s⁻¹ along with two lines at −53 and −31 km s⁻¹ with a width of ∼10 km s⁻¹ each.

Zoom In Zoom Out Reset image size

Figure 3 shows the [C ii] spectra toward the "red," "center," and "blue" positions, where the baseline noise is 0.12 K, 0.09 K, and 0.10 K, respectively. We used the central 80% of the spectral window during baseline fitting to avoid large baseline variation toward the edge of the spectral window. We detect a narrow feature near the source velocity in the "center" and "red" spectra (see Figure 4 for a zoomed-in view). These narrow features have a low signal-to-noise ratio (S/N). By fitting a Gaussian profile, we constrain the narrow feature at 4.8 ± 0.3 km s⁻¹ ("center" position) to a width of 2.3 ±0.6 km s⁻¹, a peak temperature of 0.39 ± 0.09 K, and an S/N of 4.4, while, at the "red" position, the narrow feature is centered at 6.9 ± 0.3 km s⁻¹ and is fitted with a width of 2.8 ± 0.6 km s⁻¹, a peak temperature of 0.40 ± 0.07 K, and an S/N of 3.4. We also note that the feature in the "red" spectrum has an asymmetric line profile with the blue side of the line having a steeper profile than the red side.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using Herschel/PACS, Green et al. (2016) measured line strengths of (4.2 ± 0.1) × 10⁻¹² erg s⁻¹ cm⁻² for the [O i] line from a spatial pixel (hereafter spaxel) centered on L1551 IRS 5, which has a size of 94 × 94. For these SOFIA/upGREAT observations, we derived an integrated line flux at 63 μm of (2.2 ± 0.2) × 10⁻¹² erg s⁻¹ cm⁻² using a symmetric Gaussian beam (${\rm{\Omega }}={\theta }_{\mathrm{HPBW}}^{2}\pi /4\mathrm{ln}2;$ θ_HPBW = 63). The greater line flux observed with Herschel suggests that the emission extends beyond the SOFIA beam. If we naively assume a uniform brightness distribution and scale the SOFIA line flux to the square PACS spaxel (94 × 94), the line flux becomes (4.4 ± 0.4) × 10⁻¹² erg s⁻¹ cm⁻². The consistency indicates that redshifted emission is also missing in the Herschel observations (see further discussion in Section 4.1).

Using the Kuiper Airborne Observatory (KAO), Cohen et al. (1988) measured an [O i] line flux of (5.0 ± 1.2) × 10⁻¹² erg s⁻¹ cm⁻² from a 44'' aperture, consistent with the measurements using Herschel and SOFIA. Green et al. (2016) measured a total [O i] 63 μm line flux of (7.1 ± 0.1) × 10⁻¹² erg s⁻¹ cm⁻² over the entire 47'' × 47'' field of view, consistent with the KAO line flux given the larger aperture.

Green et al. (2016) mapped the [C ii] line flux in the 47'' × 47'' region around the protostar using Herschel. From the Herschel intensity map of the [C ii] emission (Figure 1, right), we can estimate the [C ii] line flux within the SOFIA beam of 141, finding 1.1 × 10⁻¹³ and 1.2 × 10⁻¹³ erg s⁻¹ cm⁻² in the "center" and "red" positions, respectively. The narrow feature in our SOFIA observations has a line flux of 3.5 × 10⁻¹⁴ and 3.1 ×10⁻¹⁴ erg s⁻¹ cm⁻² in the "center" and "red" positions, respectively, which only represents 32% and 26% of the total line flux detected with Herschel. This discrepancy hints at a broad [C ii] emission undetected in our SOFIA observations or that SOFIA chopped out extended emission more than Herschel did (i.e., observations of SOFIA have more significant [C ii] emission at the chopped position than that of Herschel). We further discuss the line width of the undetected [C ii] emission in Section 4.3.

The [O i] emission at the "central" position shows only blueshifted emission and no detectable redshifted emission. However, we expected to detect the emission from both lobes of outflows, resulting in a more symmetric line profile centered at the source velocity. Given this puzzling asymmetric line profile, we consider a few scenarios that may explain our observations. In the case of massive protostars, colder [O i] gas in the surrounding envelope can absorb the emission from warm [O i] originating in the outflow, resulting in absorption features around the envelope velocity, which is similar to the source velocity (Karska et al. 2014a; Leurini et al. 2015). However, in this scenario, the emission should show a rather narrow absorption feature centered at the source velocity, which appears in our observations, and the high-velocity redshifted [O i] emission should have been detected (like the detected blueshifted emission). Thus, this scenario is unlikely. Asymmetric outflows where the blueshifted outflow is much stronger than the redshifted outflow would lead to an asymmetric [O i] line profile, such as the outflow of TMC 1A (Bjerkeli et al. 2016). However, the [O i] emission observed by Herschel appears quite symmetric in both outflow lobes (Figure 1, left), so this scenario is also unlikely.

The circumbinary disk can preferentially block the redshifted emission if most of the [O i] emission would originate only from the inner ∼300 au region (Cruz-Sáenz de Miera et al. 2019; Takakuwa et al. 2020). The observed [O i] line profile has a low-velocity (≲20 km s⁻¹) and a high-velocity (∼20–80 km s⁻¹) component (Figure 2; see also Section 4.2). In this scenario, both components would have a scale of ≲300 au. Disk wind may occur at the disk scale (r ∼ 150 au) and could produce the low-velocity component (Ercolano & Owen 2016; Gressel et al. 2020) but not the high-velocity component, which may be due to jets. While the jets remain spatially unresolved with our SOFIA observations, the comparison of the [O i] line profile observed by Herschel-PACS and SOFIA finds predominantly blueshifted emission between the inner 62 and 94 regions (Figure 5), suggesting that the jet extends beyond the disk scale. The missing redshifted emission at high velocity indicates that the obscuration is, in fact, beyond the disk scale; thus, the circumbinary disk is unlikely to be the main obscuration mechanism.

Zoom In Zoom Out Reset image size

Dust in the protostellar envelope can produce differential extinction in the blue- and redshifted outflows if the outflows are inclined with respect to the plane of the sky. The dust opacity from the envelope would impact the [O i] 63 μm line more than the [O i] 145 μm line because the dust opacity is typically proportional to λ^−β , where β is often assumed as ∼1.8, making the ratio of these two lines an indicator of differential extinction. From the Herschel observations, the ratio of the two [O i] lines increases toward the blueshifted outflow and decreases toward the redshifted outflow, indicative of significant dust extinction toward the redshifted outflow (Figure 6). Nisini et al. (2015) performed a similar analysis on other low-mass protostellar outflows but found no obvious variation of the [O i] line ratio. The disappearance of the redshifted outflow in the Two Micron All Sky Survey (2MASS) image of L1551 IRS 5 (Skrutskie et al. 2006) also corroborates the scenario of differential dust extinction (Figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantitatively test the effect of differential extinction, we use the envelope model of L1551 IRS 5 presented in Kristensen et al. (2012), which has a density profile

$$\begin{eqnarray}&&{n}_{\mathrm{gas}}={n}_{\mathrm{in}}{\left(\displaystyle \frac{r}{{r}_{\mathrm{in}}}\right)}^{-1.8},\mathrm{where}\ {r}_{\mathrm{in}}\lt r\lt {r}_{\mathrm{out}}.\end{eqnarray} \tag{ 1 }$$

r_in and r_out are the inner and outer radius of the envelope, and n_in is the gas number density at r_in. By fitting the spectral energy distribution and the radial brightness profiles at 450 and 850 μm, Kristensen et al. (2012) constrained these parameters to n_in ∼ 6.9 × 10⁸ cm⁻³, r_in ∼ 28.9 au, and r_out ∼ 14,000 au. To model the optical depth toward the blue- and redshifted outflows, we assume triangular outflow cavities and adopt geometric parameters from the literature.

For example, we adopt an outflow half-opening angle (θ_cav) of 22° (Wu et al. 2009) and an inclination angle (θ_incl) of 30° with respect to the plane of the sky (Chou et al. 2014). Then, to derive the column density (N_gas) at any impact parameter (b), we can simply integrate the density profile (Equation (1)) from either the inner edge of the envelope, $\sqrt{{r}_{\mathrm{in}}^{2}-{b}^{2}}$, or the outflow cavity wall, $b\,\tan ({\theta }_{\mathrm{cav}}+{\theta }_{\mathrm{incl}})$, to the outer edge of the envelope, $\sqrt{{r}_{\mathrm{out}}^{2}-{b}^{2}}$. The θ_cav is positive for the blueshifted outflow cavity and negative for the redshifted outflow cavity. With an estimate of the column density in hand, we can estimate the dust optical depth from

$$\begin{eqnarray}&&{\tau }_{\mathrm{dust}}={\kappa }_{63\,\mu {\rm{m}}}{N}_{\mathrm{gas}}\mu {m}_{{\rm{H}}}0.01,\end{eqnarray} \tag{ 2 }$$

where κ₆₃ μm = 2.1 × 10² cm⁻² g⁻¹ is the dust opacity at 63 μm taken from Table 1, Column 5 in Ossenkopf & Henning (1994), μ = 2.8 is the mean molecular weight (Kauffmann et al.2008), m_H is the mass of a hydrogen atom, and 0.01 is the dust-to-gas mass ratio.

Because the outflow is inclined, the integrated dust optical depth toward the redshifted outflow is greater than that toward the blueshifted outflow (Figure 8). Assuming both outflows have the same intrinsic [O i] line flux, the observed flux ratio between two lobes can be approximated by the ratio of Σe^−τ Δb, where τ is the dust optical depth in each outflow lobe as presented in Figure 8. The results are given in Figure 9, which shows that the dust extinction from the envelope can significantly reduce the flux ratio to ∼0.55 ± 0.15 with a SOFIA beam size of 63. The observed flux ratio upper limit agrees with the modeled ratio within the uncertainty, but is lower than the modeled ratio. If the [O i] emission originates from shocked gas, the [O i] line intensity would be higher toward the protostar, resulting in a lower flux ratio in the model that assumes a uniform brightness. With a resolution of 94, the Herschel/PACS observations (Figure 1) cannot resolve the inner ∼6'' where the redshifted outflow is mostly extincted (τ > 1 in Figure 8). Beyond the inner 94 region, the dust extinction has negligible impact on the [O i] emission, resulting in a rather symmetric distribution of [O i] emission in Figure 1. Thus, extinction due to the dusty envelope surrounding L1551 IRS 5 is the most likely reason for the absence of redshifted emission in the resolved [O i] line, although the circumbinary disk may block some low-velocity redshifted [O i] emission. The asymmetric line profiles of CO J = 16 → 15 and H₂O 2₀₂ → 1₁₁, which are discussed in Section 4.2, also corroborate with the extinction due to the envelope. As the wavelength increases, from the CO line to the H₂O line, the asymmetry becomes less obvious (Figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We further test the extinction from the envelope model using observations of [Fe ii] jets in L1551 IRS 5. Hayashi & Pyo (2009) detected a faint [Fe ii] jet in the redshifted outflow, appearing only at separations >10'' from the protostar. They also measured a visual extinction, A_V , of 20–30 mag for this faint jet. We estimate the A_V using the dust opacity at 0.55 μm (2.1 × 10⁴ cm² g⁻¹) and taking ${A}_{V}=2.5{\mathrm{log}}_{10}({e}^{{\tau }_{\lambda }})=1.086{\tau }_{V}$. Figure 10 shows the A_V along the red- and blueshifted outflows. At an offset >10'', the A_V toward the redshifted outflows is consistent with the A_V measured from the [Fe ii] jet, supporting the scenario of envelope extinction for the [O i] line.

Zoom In Zoom Out Reset image size

The broad line width and significant velocity offset of the [O i] line profile suggests an outflow origin. From observations of far-infrared emission in outflows, a composite scenario emerges from high-J CO and H₂O lines (e.g., Kristensen et al. 2012, 2017b; Mottram et al. 2017; Karska et al. 2018; Yang et al. 2018). In this scenario, the outflow-envelope interaction leads to shock-excited emission of high-J CO, H₂O, and [O i] lines. So-called "cavity" shocks at the envelope-outflow interface can produce broad (typical FWHM of >15 km s⁻¹) emission centered at the source velocity, while "spot" shocks occurring within the outflow (typical FWHM of =5–15 km s⁻¹) result in emission that is offset from the source velocity, and often blueshifted. Magneto-hydrodynamical disk winds can also produce similar features as the cavity shocks (Yvart et al. 2016). Some sources show a component similar to that of spot shocks but with significant velocity offset, and are often called "extremely high-velocity gas" or "bullets." Both components are thought to originate from spot shocks, likely tracing dissociative J-shocks (Mottram et al. 2014; van Dishoeck et al. 2021). Colder molecular gas entrained by outflows along the edges of the cavity mostly dominates the emission of low-J CO lines, such as CO J = 3 → 2, showing a narrow (FWHM < 5 km s⁻¹) line profile centered on the source velocity (Yıldız et al. 2013, Figure 11).

Zoom In Zoom Out Reset image size

The [O i] 63 μm line profile in the "center" position seems to have two components as well, including a component centered at the source velocity with its redshifted part blocked by the dusty envelope (see discussion in Section 4.1) and a broad component at high blueshifted velocity. Assuming a two-component model, we fitted two Gaussian profiles with one component fixed at the source velocity, 6.2 km s⁻¹, and constrained another broad component at −30.0 km s⁻¹ (Figure 12). Table 2 lists the fitted parameters of the two Gaussian profiles. This composite model successfully reproduces the observed line profile, suggesting two distinct components of [O i] gas; however, the fitted parameters of the extremely broad component have ∼35%–60% uncertainties.

Zoom In Zoom Out Reset image size

Kristensen et al. (2017b) decomposed the Herschel/HIFI CO J = 16 → 15 emission in L1551 IRS 5 into two components tracing different physical components. The component that is centered around the source velocity has a narrow line width of 3.6 km s⁻¹, and traces the entrained gas. Another component at 3.6 km s⁻¹ has a width of 14.9 km s⁻¹ and traces the cavity shocks (Figure 12). In comparison, the [O i] emission has no narrow component in emission but does show a narrow absorption near the source velocity (at ∼5.4 km s⁻¹). This narrow absorption may come from the colder [O i] gas in the quiescent or the infalling envelope. The broad component in CO J = 16 → 15 has a similar line profile to the [O i] component centered at the source velocity. The centroid of the CO broad component is 3.5 km s⁻¹, which is consistent with the [O i] component centered at the source velocity given the uncertainty in the [O i] spectrum, suggesting that both lines trace the component.

Four H₂O lines were detected by Herschel (Kristensen et al. 2012; Mottram et al. 2014). The 2₀₂ → 1₁₁ and 2₁₁ → 2₀₂ lines show a narrow component (FWHM = 4.3 ± 0.8 km s⁻¹) centered at the source velocity with a blueshifted line wing (the H₂O 2₀₂ → 1₁₁ emission is shown in Figure 11). However, the other two H₂O lines, 1₁₁ → 0₀₀ and 1₁₀ → 1₀₁, show a broad (FWHM = 26.1 ± 6.3 km s⁻¹) redshifted component at 19.2 ± 3.0 km s⁻¹ but no narrow component. Figure 12 also shows the decoupled components derived from H₂O lines (Mottram et al. 2014) and the blueshifted component from a double-Gaussian fit to the H₂O 2₀₂ → 1₁₁ line (Figure 12). At blueshifted velocities, the H₂O lines in L1551 IRS 5 support a similar scenario as the CO J = 16 → 15 line, where the emission comes from entrained gas and cavity shocks traced by the blueshifted line wing. At redshifted velocities, the H₂O lines show a velocity component at 19.2 ± 3.0 km s⁻¹ with a width of 26.1 ± 6.3 km s⁻¹, suggesting an origin in cavity shocks (Mottram et al. 2014). The reason why the four H₂O lines exhibit different kinematics remains unclear. A slightly blueshifted narrow absorption also appears in the 1₁₁ → 0₀₀ H₂O line, indicative of absorption by a colder envelope (Mottram et al. 2014).

Among the spectra of high-J CO, H₂O, and [O i], only the [O i] emission has the extremely broad component with a width of ∼90 km s⁻¹. Combining the insights learned from the [O i] 63 μm and CO J = 16 → 15 lines, a revised picture of the L1551 IRS 5 outflow emerges: the narrow (FWHM = 3.6 km s⁻¹) component traces entrained material that appears in emission for CO J = 16 → 15 and H₂O (Mottram et al. 2014). The quiescent envelope appears as narrow absorption seen in both the H₂O 1₁₁ → 0₀₀ line and the [O i] 63 μm line. The broad component (FWHM ∼ 20 km s⁻¹) comes from the cavity shocks where the outflow drives shocks into the cavity walls. The extremely broad component (FWHM ∼ 90 km s⁻¹) that is only seen in [O i] traces J-type spot shocks in the outflow cavity or the jet. The line profiles of all three species, [O i], high-J CO, and H₂O, show asymmetries that agree with the extinction due to the envelope (see Section 4.1).

Green et al. (2016) detected [O i] emission at 63 and 145 μm using Herschel/PACS, further constraining the excitation condition of [O i] in L1551 IRS 5. The ratio of [O i] 63 μm/145 μm is 32.0 ± 4.4 for the central spaxel and 25.4 ± 2.5 for the 1D spectrum extracted with a 248 aperture, which best matches the the Spectral and Photometric Imaging Receiver spectrum (Yang et al. 2018, see also Figure 6 for the map of ratios). Nisini et al. (2015) modeled the ratio of the [O i] 63 to 145 μm lines, which increases with the density of collision partners, suggesting a H₂ density of ∼10 ^5.5 cm⁻³ for the [O i] outflow.

The study of the R1 position in the outflow of NGC 1333 IRAS 4A is a notable comparison to L1551 IRS 5 (Kristensen et al. 2017a). In R1, the [O i] 63 μm emission observed by SOFIA has a similar line profile as that of H₂O and high-J CO lines, which are tracers of shocks, suggesting that these emission lines mostly come from the outflows. In contrast, the H₂O line in L1551 IRS 5 has a different line profile compared to the [O i] spectrum, suggesting that the [O i] emission uniquely trace an atomic outflow and/or jet.

Both shocked gas in outflows and photodissociation regions (PDRs) can contribute to the [O i] emission (Hollenbach 1985; Tielens & Hollenbach 1985). Although our SOFIA observations did not detect significant [O i] emission at "red" and "blue" outflow positions, Green et al. (2016) detected extended [O i] emission with Herschel/PACS on ∼20'' scales (Figure 1). Interpolating from the Herschel [O i] emission map, we estimate [O i] velocity-unresolved fluxes of 7.7 × 10⁻¹³ and 7.6 × 10⁻¹³ erg s⁻¹ cm⁻² in the "red" and "blue" positions, respectively. If the [O i] flux observed by Herschel comes from a single line, our SOFIA observations have a better sensitivity to detect a narrow (FWHM ≲ 5 km s⁻¹) line from PDRs compared to a broad line due to shocks. The non-detection of narrow [O i] emission in all three pointings allows us to constrain the maximum contribution to the [O i] emission from PDRs by comparing the non-detection to the interpolated [O i] flux from Herschel/PACS. [C ii] emission is a common tracer of PDRs (Hollenbach & Tielens 1999; Kaufman et al. 1999). If we take the mean FWHM of the narrow [C ii] lines, 2.5 km s⁻¹, as the characteristic line width of gas in PDRs, the upper limit of the [O i] line flux is then 2.4 ×10⁻¹⁴ erg s⁻¹ cm⁻² at both outflow positions and 3.7 × 10⁻¹⁴ erg s⁻¹ cm⁻² at the central position. Comparing to the Herschel [O i] line fluxes estimated within the "central" (from Section 3.3), "red," and "blue" positions, the upper limit of the PDR contribution to the [O i] flux is 3.1% in "red" and "blue" positions, while the contribution to the [O i] flux is 2.6% in the "central" position, suggesting that the shocked gas (i.e., broad emission) dominates the [O i] flux. Assuming that a single Gaussian line profile can describe all of the non-PDR [O i] emission, [O i] emission at "red" and "blue" positions should have FWHM greater than ∼80 km s⁻¹ had it been detected, which is similar to the line width of the extremely broad [O i] component at the "central" position.

The narrow line width of the observed [C ii] emission indicates a PDR origin. Extracting the line flux from the Herschel intensity map (Figure 1, right) using the SOFIA beam (141 at the [C ii] line), the detected narrow feature makes up 30% and 26% of the total Herschel [C ii] flux at the "center" and "red" positions, respectively. If the [C ii] emission from the large-scale cloud is removed by beam-chopping (see Section 2), the detected [C ii] line flux then traces the emission from local PDRs. The PDRs contribute ∼30% of the total [C ii] flux but only represent ≲3% of the [O i] flux, suggesting that [O i] emission is indeed a good tracer of shocked gas, and the contribution from PDRs is negligible. Furthermore, the spatial extent of the [O i] emission coincides with the blueshifted [Fe ii] emission that unambiguously traces dissociative shocks in the jets (Pyo et al. 2009). If the discrepant [C ii] emission between observations of SOFIA and Herschel is instead due to excessive baseline subtraction during data reduction, we can assume a Gaussian profile for the undetected [C ii] emission and estimate the lower limit of its line width given the spectral noise. We estimate that the [C ii] flux missing in the SOFIA observations is 7.7 × 10⁻¹⁴ and 8.9 ×10⁻¹⁴ erg s⁻¹ cm⁻², resulting in a lower limit of 23.3 and 20.2 km s⁻¹ for the line width at the "center" and "red" positions, respectively. Comparing to the lower limit of the [O i] line width, the undetected [C ii] emission could be narrower than the undetected [O i] emission.

With the PDR contribution constrained toward the protostar, we can further use the line ratio between [C ii] and [O i] to infer the UV radiation incident on the PDR. Karska et al. (2018) calculated the ratio between the [C ii] 158 μm and the [O i] 63 μm lines as a function of G₀ from 10¹ to 10^6.5 and densities from 10³ to 10⁷ cm⁻³. For the PDR contribution in the "center" position, we measure a [C ii] line flux of 3.3 × 10⁻¹⁴ erg s⁻¹ cm⁻² in a 141 beam, and an upper limit of the [O i] line flux of 3.7 × 10⁻¹⁴ erg s⁻¹ cm⁻² in a 63 beam. Thus, if we scale the [O i] upper limit to the [C ii] beam size, we have a lower limit of 0.17 on the ratio of [C ii] and [O i]. The corresponding [C ii] integrated line intensity is 9.0 × 10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹ and, combining with the [C ii] intensity, we estimate that the incident UV radiation, G₀, is at most ∼10 for a density of ∼10⁵ cm⁻³, which is similar to the density estimated from the ratio of [O i] 63 μm and 145 μm lines (see Section 4.2; Nisini et al.2015).

A strength of ∼10 for the UV radiation is consistent with prior estimates for low-mass protostars, which assumed that 90% of unresolved [O i] emission originates from shocks and 10% from PDRs (Karska et al. 2018). The velocity-resolved line profiles from SOFIA/upGREAT verify such assumptions by identifying the narrow and broad components in [O i] and [C ii] lines.

Besides being a good tracer of energetic outflows and jets, atomic oxygen is a major carrier of volatile oxygen in protostellar shocks. The Herschel Key Program, Water In Star-forming regions with Herschel (WISH), demonstrated that CO, H₂O, and atomic O represent most of the volatile oxygen budget in protostellar outflows and shocks, but up to 50% of the total oxygen abundance remains unaccounted for (van Dishoeck et al. 2021, and references therein). Velocity-resolved spectra of [O i], CO J = 16 → 15, and H₂O give us a unique opportunity to identify the emission from shocks as well as to quantify the column density of each oxygen carrier in shocks. The two components of the [O i] emission (Table 2) represent two types of shocks. The broad component (FWHM ∼ 15–25 km s⁻¹) seen in [O i], CO J = 16 → 15, and H₂O trace cavity shocks, while the extremely broad component only seen in [O i] probes the shocked gas in jets. The absence of an extremely broad component in the high-J CO and H₂O lines indicates fully dissociative shocks. Tracing the ionized jet, the [Fe ii] emission from the ionized jet exhibits an extremely broad (FWHM = 150–180 km s⁻¹) component centered at ∼−120 km s⁻¹ and a narrower (FWHM = 40 km s⁻¹) component at −300 to −400 km s⁻¹ (Pyo et al. 2009). While our observations only show the far-infrared [O i] line up to ∼−100 km s⁻¹, the optical [O i] line at 6364 Å has a broad (FWHM ∼ 200 km s⁻¹) line at ∼−100 to −300 km s⁻¹ similar to the characteristics of the [Fe ii] emission (White & Hillenbrand 2004), suggesting that both [O i] and [Fe ii] trace the jet. The brightness of this extremely high-velocity component likely falls below the sensitivity of our observations.

We can further estimate the oxygen budget in the shocks. To estimate the [O i] column density for the broad component centered at the source velocity (Figure 13, orange line), we used the non-LTE radiative transfer code, Radex (van der Tak et al. 2007) and the collision rates from the Leiden Atomic and Molecular Database (LAMDA; Schöier et al. 2005). While the ratio of [O i] 63 and 145 μm line fluxes indicates an H₂ density of 10^5.5 cm⁻³ (Nisini et al. 2015), both the kinetic temperature (T_kin) and the column density of [O i] remains unknown. Using the line fitting results listed in Table 2, we modeled an [O i] column density of (1–5) × 10¹⁶ cm⁻² with a kinetic temperature ranging from 500 to 8000 K (Figure 13). The modeling used the [O i] atomic data from LAMDA (Schöier et al. 2005) and assumed a uniform geometry for calculating the escape probability to reproduce the observed [O i] line flux within 1%. The collision rates from LAMDA only cover the first three states of [O i]. To test the necessity of including higher states that produce the 6300 and 6364 Å lines, we further include the ¹D₂ and ¹S₀ states using Table F.3 and Equation (2.27) of Draine (2010), which employs the data from Pequignot (1990). Assuming a typical electron density of 2 × 10⁴ cm⁻³ in protostellar jets (Hartigan et al. 2019), the calculations with the updated collision rates show that these excited states have negligible impact on the atomic oxygen column density derived from the ³P ₁ → ³ P₂ transition. The uncertainty of the [O i] column density is ∼40% due to the uncertainty in the measured flux. The modeled [O i] column density only varies by a factor of ∼2 in 8000 K > T_kin > 750 K. Because both atomic and molecular lines are present in the cavity shocks, the temperature has to be moderate to prevent dissociation. Therefore, we assume that the cavity shocks in L1551 IRS 5 have the same temperature as the shocks in NGC 1333 IRAS 4A R1, 750 K (Kristensen et al. 2017a), which then gives a [O i] column density of (2.0 ± 0.7) × 10¹⁶ cm⁻². Assuming a temperature of 750 K and an H₂ density of 10^5.5 cm⁻³, we can also model the broad component of the CO and H₂O emission using Radex. For the CO J = 16 → 15 line, the broad (FWHM = 14.9 km s⁻¹) component has a peak intensity of 0.36 K (Figure 12), resulting in a column density of (2.7 ± 0.2) × 10¹⁵ cm⁻². For H₂O, we modeled the red- and blueshifted components separately. The modeling suggests a column density of (5.1 ± 2.6) × 10¹² cm⁻² for the redshifted component in the 1₁₁ → 0₀₀ line with an FWHM of 26.1 km s⁻¹ and a peak intensity of 0.02 K; the modeling meanwhile gives a column density of (5.0 ± 2.3) ×10¹² cm⁻² for the blueshifted component in the 2₀₂ → 1₁₁ line with an FWHM of 8.4 km s⁻¹ and a peak intensity of 0.03 K.

Zoom In Zoom Out Reset image size

To derive the abundance of dominant oxygen carriers, we need to estimate the mass of each species in the same location. Because of the different beam sizes for [O i], CO, and H₂O observations, we examine two scenarios to derive the range of possible abundances. First, we assume that the shocks are spatially unresolved for all three lines. In this case, we can derive the mass of each species by multiplying the column density with the respective beam size. The SOFIA beam for [O i] is 63 and the Herschel beams for CO J = 16 → 15 and H₂O lines are 115 and 22'', respectively. The relative abundance of each species can be written as

$$\begin{eqnarray}&&\displaystyle \frac{X({\rm{Y}})}{X({{\rm{O}}}_{\mathrm{total}})}=\displaystyle \frac{N({\rm{Y}}){{\rm{\Omega }}}_{{\rm{Y}}}}{N({\rm{O}}){{\rm{\Omega }}}_{{\rm{O}}}+N(\mathrm{CO}){{\rm{\Omega }}}_{\mathrm{CO}}+N({{\rm{H}}}_{2}{\rm{O}}){{\rm{\Omega }}}_{{{\rm{H}}}_{2}{\rm{O}}}},\end{eqnarray} \tag{ 3 }$$

where Ω_Y is the beam area for the observations of each species. In the second scenario, we assume the shocks are uniformly spread out over the beam. A larger beam would probe more shocked gas. Therefore, we derive the abundance of each species using the smallest beam size (i.e., 63 for [O i]) instead of their respective beams. In both cases, we assume that all three species trace the same shocks. Furthermore, we only consider the gaseous oxygen carriers and the refractory oxygen remains locked on dust grains because of the SiO non-detection (Codella et al. 1999). Any refractory oxygen released from dust grains will make the relative abundances in Table 3 upper limits. For cavity shocks and/or disk winds, atomic oxygen is the primary carrier of oxygen, representing 69% and 88% of the gaseous oxygen abundance in the extended and unresolved scenarios, respectively, with CO as the secondary oxygen carrier (Table 3).

To estimate the absolute abundances of these major oxygen carriers, we need to make assumptions about the total gaseous oxygen abundance. There is a long-standing problem of missing oxygen, where the elemental oxygen abundance in the interstellar medium, 5.8 × 10⁻⁴ (Przybilla et al. 2008), is consistently higher than all of the observable oxygen carriers in star-forming regions (e.g., van Dishoeck et al. 2021). This so-called undetected depleted oxygen (UDO) abundance with respect to H is ∼1–2 × 10⁻⁴ in star-forming regions (van Dishoeck et al. 2021). If we assume a UDO abundance of 2 × 10⁻⁴, the remaining volatile oxygen abundance is 3.8 × 10⁻⁴, similar to the value of (3.4 ± 0.5) × 10⁻⁴ measured from using translucent sight lines toward stars (Sofia et al. 2004). Furthermore, refractory dust also carries oxygen. Following the assumption made in van Dishoeck et al. (2021), we assume a refractory abundance of 1.4 × 10⁻⁴ derived from the Mg, Si, and Fe abundances in Przybilla et al. (2008). Henceforth, we assume a total gaseous oxygen abundance of 2.0 × 10⁻⁴. From Equation (3), we then derive an atomic oxygen abundance relative to H for the two scenarios described above ("extended" and "unresolved" in Table 3). The dominant presence of [O i] is expected in L1551 IRS 5; Karska et al. (2018) found the [O i] luminosity dominates the line luminosity of major coolants, such as H₂O, CO, and OH. In L1551 IRS 5, the [O i] luminosity takes up 55% of the total line luminosity, while the warm and hot CO luminosity takes up 38%. The X(O) is shockingly high compared to the oxygen abundance in other observations of protostellar shocks. In the R1 shock knot of NGC 1333 IRAS 4A, for instance, Kristensen et al. (2017a) analyzed the velocity-resolved [O i] emission and estimated an X(O) of 5–6 × 10⁻⁵, taking up only ∼15% of the total oxygen volatile abundance.

The ratio of the [O i] to H₂O emission can constrain the velocity where H₂O becomes dissociated. For example, in irradiated C-shock models, the O abundance decreases with shock velocity, while the H₂O abundance increases, resulting in a decreasing [O i]/H₂O ratio (Melnick & Kaufman 2015; Godard et al. 2019), the opposite to this SOFIA observation shows. A high UV radiation would promote photodissociation, efficiently destroying H₂O; however, Lehmann et al. (2020) found that only J-shocks can produce sufficient UV radiation to generate significant photodissociation. Smooth, steady-state disk winds can reproduce the Herschel H₂O observations that trace outflows and contain little dust to block the UV radiation due to the accreting protostar (Yvart et al. 2016). Thus, the cavity shock component can be a C-type shock irradiated by an exceptional amount of UV radiation from sources other than the shock itself, a J-type shock that shows increasing O abundance with velocities (Lehmann et al. 2020), or disk winds. In Figure 14, we show the [O i]/H₂O intensity ratio as a function of velocity. The ratio has a minimum at ∼9 km s⁻¹, a value of ∼5 at the source velocity, 6.2 km s⁻¹, and then increases to ∼60 at −10 km s⁻¹.

Zoom In Zoom Out Reset image size

Only the [O i] emission shows an extremely broad component, suggesting that the shocks are fully dissociative. The modeled column density of this component is (2.7 ± 1.6) ×10¹⁶ cm⁻², assuming a kinetic temperature of 750 K. The models of far-UV-irradiated C-shocks break down at velocity ∼20 km s⁻¹, assuming a pre-shock density of 10⁵ cm⁻³, becoming dissociative at larger velocities (Melnick & Kaufman 2015). The lack of H₂O emission beyond v ≲ −15 km s⁻¹ is consistent with its photodissociation in a J-type shock. Hollenbach (1985) presented a simple J-shock model, suggesting that [O i] emission dominates the cooling at T ≲ 5000 K. If we take T = 5000 K as an upper limit of the [O i] temperature for the extremely broad component, we can estimate a lower limit of [O i] column density of 2 × 10¹⁶ cm⁻² from the non-LTE radiative transfer modeling (Figure 13).

A long-standing issue with molecular outflows is whether they can be driven by stellar winds or jets. The basic idea is that the wind emerges from the star and/or disk at high speeds and sweeps out a cavity in the infalling core. The swept-up gas is observed in low rotational transitions of CO. The highly radiative nature of the shocks in the molecular gas implies that only the momentum of the wind is available to sweep up material (so-called "momentum-conserving" situations). In this picture, the momentum in the CO flow must be matched by the momentum of the wind,

$$\begin{eqnarray}&&{P}_{{\rm{w}}}={P}_{\mathrm{CO}}={M}_{{\rm{w}}}{v}_{{\rm{w}}}={M}_{\mathrm{CO}}{v}_{\mathrm{CO}}.\end{eqnarray} \tag{ 4 }$$

If this is correct, an intrinsic mass-loss rate in the wind can be derived from

$$\begin{eqnarray}&&{\dot{M}}_{\star }={P}_{\mathrm{CO}}/({t}_{\mathrm{CO}}{v}_{{\rm{w}}})={F}_{\mathrm{CO}}/{v}_{{\rm{w}}},\end{eqnarray} \tag{ 5 }$$

where F_CO is the average force that must have been applied to the swept-up gas, v_CO is the velocity of the swept-up gas, and v_w is the velocity of the stellar wind. The intrinsic mass-loss rate is averaged over the age of the flow, t_CO. This basic idea was developed by Levreault (1988) and applied to CO maps of young stellar objects, including L1551 IRS 5. Maps of CO are readily obtained, while observations of the wind itself have proved challenging. Since the latter inevitably average over much shorter times, comparison of the two provides a potential means to measure changes, whether secular or episodic, in the intrinsic mass-loss rate (${\dot{M}}_{\star }$), which in turn, is likely to be closely related to the mass accretion rate onto the star, a quantity of considerable interest. The binary protostars of L1551 IRS 5 drive two jets, as can be seen in [Fe ii] (Fridlund & Liseau 1998; Pyo et al. 2009), but remain unresolved with [O i] and other tracers. Thus, we are measuring the total intrinsic mass-loss rate instead of rates for individual jets.

There have been many maps of different transitions of CO toward L1551 IRS 5, with a wide range of inferred momenta. For the most stringent test of the constancy of the mass-loss rate, we prefer the largest, most complete maps of the lower J transitions of CO, which turn out to be the oldest. The outflow was the first bipolar flow discovered (Snell et al. 1980), but better maps were made in the J = 1 → 0 transition by Snell & Schloerb (1985) and the J = 2 → 1 transition by Levreault (1988) both covering about 30'.

Determination of the force from the observations can take various paths, with issues of how to handle the fact that the outflow is inclined to the plane of sky (θ_incl). Some definitions of the inclination angle use the angle between the outflow axis and the plane of the sky, while others use the line of sight to the observer as the reference. We adopt the former. The L1551 IRS 5 outflow is mostly in the plane of the sky, and we adopt θ_incl = 30°, based on a determination of the disk angle (90-θ_incl) of 60° ± 5° (Chou et al. 2014).

Following Dunham et al. (2014), we make no correction to the observed P_CO (see also Downes & Cabrit 2007) but correct t_CO because both length and velocity are affected. For our convention on inclination angle,

$$\begin{eqnarray}&&{t}_{\mathrm{CO}}={t}_{\mathrm{obs}}\times \tan ({\theta }_{\mathrm{incl}}).\end{eqnarray} \tag{ 6 }$$

This approach is similar to Method M7 of van der Marel et al. (2013). Therefore, forces determined from observations, with no corrections, should be multiplied by $1/\tan ({\theta }_{\mathrm{incl}})=1.73$. We added forces for the blue and red lobes. The result for the J = 1 → 0 map of Snell & Schloerb (1985), after correction of age by inclination angle, is F_CO = 1.2 × 10⁻⁴ M_⊙ yr⁻¹ km s⁻¹ averaged over an age of 4.6 × 10⁴ yr. The result from the J = 2 → 1 map from Levreault (1988) is F_CO = 2.1 × 10⁻⁴ M_⊙ yr⁻¹ km s⁻¹ averaged over an age of 2.9 × 10⁴ yr corrected by inclination angle.

The wind velocity has been estimated from spectra of [Fe ii] λ = 1.644 μm emission. Pyo et al. (2005) find velocities as high as 300 km s⁻¹, with blueshifted components at 100 km s⁻¹ to 279 ± 10 km s⁻¹. Giovanardi et al. (1992) also observed velocities up to 150 km s⁻¹ in H i. Correcting for inclination by $1/\sin ({\theta }_{\mathrm{incl}})=2$ implies v_w = 200 to 560 km s⁻¹ for ionized gas and v_w = 260 km s⁻¹ for the neutral gas. Such high velocities strongly indicate an origin for the wind near the star, rather than a disk wind. The SOFIA [O i] observations only detect [O i] up to ∼75 km s⁻¹ in blueshift, which may be due to insufficient sensitivity. We used v_w = 300 km s⁻¹ as a compromise. The values for ${\dot{M}}_{\star }$ are then 4.0 × 10⁻⁷ M_⊙ yr⁻¹ for J = 1 → 0 and 7.0 × 10⁻⁷ M_⊙ yr⁻¹ for J = 2 → 1.

Yıldız et al. (2015) also mapped the outflow of L1551 IRS 5 using the CO J = 3 → 2 line to derive the force and dynamical age of the blue- and redshifted lobes; however, the map size is only ∼5' × 5'. We use their Equation (4) to retrieve the momentum without any correction, 0.11 and 0.40 M_⊙ km s⁻¹, for the blue- and redshifted outflows, respectively. They also measured a dynamical age of 8.3 × 10³ and 6.7 × 10³ yr uncorrected for inclination in the blue- and redshifted outflows, respectively. Using Equation (5) to correct the averaged age for inclination angle, and v_w = 300 km s⁻¹, we derive an intrinsic mass-loss rate of 4.0 × 10⁻⁷ M_⊙ yr⁻¹, consistent with the estimates from low-J CO emission.

In the momentum-conserving scenario, atomic winds are thought to drive the molecular outflows because the ionized wind/jet has too little momentum (Mundt et al. 1987). Giovanardi et al. (1992) observed the high-velocity wings of H i line profiles using the Arecibo Observatory. Assuming optically thin emission and mass conservation of H i, they modeled the line profiles with a decelerating wind model to get ${\dot{M}}_{\star }=2.4$ × 10⁻⁶ M_⊙ yr⁻¹. However, their model assumes that all H i gas is due to the wind, while photodissociation in the envelope can produce H i. Moreover, the H i spectra suffer from intervening structured clouds at the blueshifted velocity. Thus, our SOFIA [O i] observations would better measure ${\dot{M}}_{\star }$ in atomic winds and so, in the following discussion, we focus on the extremely broad component of [O i].

[O i] is a good tracer of atomic outflows in the model proposed by Hollenbach (1985). The model has a fast (v_w) wind from the star–disk system driving a shock into the infalling molecular gas, pushing it away from the forming star at speeds of a few to tens of kilometers per second (v_m). The slowing of the wind in the molecular shock then drives a reverse shock into the wind at v_wsh = v_w − v_m ∼ v_w km s⁻¹. In this model, the molecular shock is formed in the ambient material, producing the emission from CO, H₂, H₂O, and OH. The rapid deceleration of the wind as it encounters the molecular material produces the wind shocks, resulting in a J-type shocked region, with T_K = 10³–10⁵ K. Such a high temperature would dissociate most molecules, making the fine-structure [O i] 63 μm line the dominant coolant in the shocked wind. The shocks are strongly radiative, and only the momentum of the wind is assumed to be available to accelerate the molecular outflow. Thus, in this model, [O i] 63 μm emission can conveniently be related to the intrinsic mass-loss rate by (Hollenbach 1985)

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{[{\rm{O}}\,{\rm\small{I}}]\,63\,\mu {\rm{m}}}}{{L}_{\odot }}=\displaystyle \frac{0.1{\dot{M}}_{\star }}{{10}^{-5}{M}_{\odot }\,{\mathrm{yr}}^{-1}}.\end{eqnarray} \tag{ 7 }$$

Using the flux of the extremely broad [O i] component that traces high-velocity gas (Table 2), we can derive ${\dot{M}}_{\star }=3.4\,\pm 1.8\,\times \,{10}^{-7}$ M_⊙ yr⁻¹ using the observed [O i] luminosity, 3.4 ± 1.8 × 10⁻³ L_⊙, assuming that the missing redshifted [O i] emission has the same luminosity as the observed blueshifted part. If we adopt the total [O i] 63 μm luminosity from Herschel observations, 4.8 ± 0.1 × 10⁻³ L_⊙ (Green et al. 2016), ${\dot{M}}_{\star }$ becomes 4.8 ± 0.1 × 10⁻⁷ M_⊙ yr⁻¹. Sperling et al. (2021) also estimated ${\dot{M}}_{\star }$ using the [O i] emission observed with SOFIA/FIFI-LS. They employed two methods, one assuming collisionally excited [O i] emission from purely atomic oxygen gas (Sperling et al. 2020) and the other one using the relation in Hollenbach (1985), deriving ${\dot{M}}_{\star }$ of 5.8–11.8 × 10⁻⁷ and 4.9–5.4 × 10⁻⁷ M_⊙ yr⁻¹, respectively.

Equation (7) provides a zeroth-order estimate of the intrinsic mass-loss rate and there are important caveats. This method assumes the entire wind shock radiates before any gas gets shock-heated or entrained, which only occurs when measuring the jet close to its source (Hartigan et al. 2019). Second, if the observations probe several unresolved shocks instead of a single shock, Equation (7) may overestimate the intrinsic mass-loss rate (Nisini et al. 2015). Finally, if parts of the wind propagate without intervention, this method would underestimate the intrinsic mass-loss rate (Cohen et al. 1988). With the velocity-resolved SOFIA [O i] observations, we can directly estimate the mass outflow rate traced by [O i] to partially circumvent these caveats, especially the first one. If the extremely broad component of [O i] traces the atomic wind that carries most momentum in the wind, we can simply estimate ${\dot{M}}_{\star }$ by taking the ratio of the mass traced by [O i] and the dynamical time. With a column density of 2.7 ±1.6 × 10¹⁶ cm⁻², we can estimate the mass outflow rate using

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\star } & = & \displaystyle \frac{\mu {m}_{{\rm{H}}}{N}_{[{\rm{O}}\,{\rm\small{I}}]}{{\rm{\Omega }}}_{{\rm{b}}}/{X}_{{\rm{O}}}}{{t}_{\mathrm{dyn}}};\\ {t}_{\mathrm{dyn}} & = & \displaystyle \frac{{R}_{\mathrm{beam}}}{{v}_{[{\rm{O}}\,{\rm\small{I}}]}}\,\tan ({\theta }_{\mathrm{incl}}),\end{array}\end{eqnarray} \tag{ 8 }$$

where μ is the mean molecular weight of 2.8 (Kauffmann et al. 2008), Ω_b (1.13 × 63²) is the area of the Gaussian beam with FWHM = 63, X_O is the atomic oxygen abundance, t_dyn is the inclination-corrected dynamical time of the [O i] gas, R_beam is the radius of the beam, and v_{[O I]} is the velocity of [O i] gas. We estimate the dynamical age of the wind as the time required for the gas to cross the beam radius with correction of inclination angle. The t_dyn is 35 yr. The atomic oxygen abundance is a key parameter to relate the [O i] mass flux to the intrinsic mass outflow rate. Because neither the H₂O nor high-J CO emission shows an extremely broad component, this component is likely to be fully dissociated. Thus, if we assume a total oxygen volatile abundance (X_O), 3.4 × 10⁻⁴, as discussed in Section 4.4, and an equal rate in each outflow lobe, the intrinsic mass outflow rate becomes 1.3 ± 0.8 × 10⁻⁶ M_⊙ yr⁻¹. Considering the uncertainty, the intrinsic mass outflow rate from [O i] agrees with not only the estimate from H i but also our estimate from the CO outflows after inclination correction, unambiguously confirming the momentum-driven outflows. In this scenario, both atomic and ionized winds belong to the mass ejected near the protostar, but the atomic wind carries most of the momentum, driving the molecular outflows. Given its spectral profile in the optical and near-infrared, some studies consider L1551 IRS 5 as an FU Orionis-like object (e.g., Mundt et al. 1985; Connelley & Greene 2010), which implies it should have a strongly varying time-dependent mass-loss rate. But the similar ${\dot{M}}_{\star }$ found over 3–5 × 10⁴ yr (low-J CO) and 35 yr ([O i]) timescales suggests no evidence of varying ${\dot{M}}_{\star }$ within a factor of three.