# An APEX Study of Molecular Outflows in FUor-type Stars

Fernando Cruz-Sáenz de Miera, Ágnes Kóspál, Péter Ábrahám, Timea

Csengeri, Orsolya Fehér, Rolf Güsten, Thomas Henning

► To cite this version:

Fernando Cruz-Sáenz de Miera, Ágnes Kóspál, Péter Ábrahám, Timea Csengeri, Orsolya Fehér, et al. An APEX Study of Molecular Outflows in FUor-type Stars. The Astrophysical Journal, 2023, 945 (1), pp.80. 10.3847/1538-4357/acb332. hal-04298978

HAL Id: hal-04298978

https://hal.science/hal-04298978v1

Submitted on 23 Nov 2023

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# An APEX Study of Molecular Outflows in FUor-type Stars

Fernando Cruz-Sáenz de Miera<fnref n="1" />,<fnref n="2" /> ID, Ágnes Kóspál<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> ID, Péter Ábrahám<fnref n="1" />,<fnref n="2" />,<fnref n="4" /> ID, Timea Csengeri<fnref n="5" /> ID, Orsolya Fehér<fnref n="6" />,<fnref n="7" /> ID,
Rolf Güsten<fnref n="8" />, and Thomas Henning<fnref n="3" /> ID

Received 2022 August 8; revised 2023 January 7; accepted 2023 January 7; published 2023 March 8

## Abstract

The FU Orionis-type objects (FUors) are low-mass pre-main-sequence objects that go through a short-lived phase (~100 yr) of increased mass accretion rate (from $10^{-8}$ to $10^{-4} M_{\odot}$ yr$^{-1}$). These eruptive young stars are in the early stages of stellar evolution and thus still deeply embedded in a massive envelope that feeds material to the circumstellar disk that is then accreted onto the star. Some FUors drive molecular outflows, i.e., low-velocity wide-angle magnetohydrodynamical winds, that inject energy and momentum back to the surrounding envelopes and help clear the material surrounding the young star. Here we present a $^{12}$CO (3–2), $^{13}$CO (3–2), and $^{12}$CO (4–3) survey of 20 FUor-type eruptive young stars observed with APEX. We use our $^{13}$CO (3–2) observations to measure the masses of the envelopes surrounding each FUor and find an agreement with the FUor evolutionary trend found from the 10 $\mu$m silicate feature. We find outflows in 11 FUors, calculate their masses and other kinematic properties, and compare these with those of outflows found around quiescent young stellar objects gathered from the literature. This comparison indicates that outflows in FUors are more massive than outflows in quiescent sources, and that FUor outflows have a higher-ratio outflow mass with respect to the envelope than the quiescent sample, indicating that the eruptive young stars have lower star-forming efficiencies. Finally, we find that the outflow forces in FUors are similar to those of quiescent young stellar objects, indicating that their accretion histories are similar or that the FUor outflows have lower velocities.

*Unified Astronomy Thesaurus concepts: Low mass stars (2050); Stellar winds (1636); Eruptive variable stars (476); FU Orionis stars (553); Circumstellar envelopes (237)*

*Supporting material: figure sets*

## 1. Introduction

Jets and molecular outflows are ubiquitous phenomena in the process of star formation. The former are highly collimated gas streams at high velocities (≥100 km s$^{-1}$), and the latter have wider opening angles and velocities between 1 and 50 km s$^{-1}$ in the case of low-mass stars. Jets are detected with optical, near-infrared, radio molecular lines, and radio continuum, while the slower outflows are typically detected with molecular line tracers (Frank et al. 2014; Bally 2016).

Both types of mass ejection events are driven by accretion; thus, the physical properties of the outflows depend on the accretion history of the star. Indeed, evidence has shown that Class 0 objects (i.e., younger protostars with higher mass accretion rates) have elevated outflow mass-loss rates and higher outflow forces compared to more evolved Class I or II objects (Mottram et al. 2017). The mass accretion rates from protostellar disks to protostars are expected to undergo episodic variations. Detailed analysis of jet knots (e.g., Ellerbroek et al. 2014; Lee et al. 2017; Garufi et al. 2019) and molecular outflow shells (Plunkett et al. 2015; Zhang et al. 2019; Nony et al. 2020; Vazzano et al. 2021)

shows how the study of outflows can shed light on the accretion history of the protostars that drive them.

The FU Orionis-type objects (FUors) are examples of the episodic nature of accretion (Hartmann & Kenyon 1996; Audard et al. 2014; Fischer et al. 2022). These eruptive young stars are low-mass protostars characterized by a sudden increase in their mass accretion rate, going from typical values of ~$10^{-8}$ up to ~$10^{-4} M_{\odot}$ yr$^{-1}$. These events are typically detected as a 3–5 mag brightening at optical and near-infrared wavelengths and are expected to last up to a century, meaning that they increase the final stellar mass by a significant amount. The FUor-type events generally occur in Class I objects. Accretion outbursts have been detected in earlier stages, e.g., the Class 0 HOPS 383 (Safron et al. 2015), and later stages, e.g., the Class II Gaia20eae (Cruz-Sáenz de Miera et al. 2022); however, these are not considered FUors. This differentiation is because to classify an object as a FUor, the near-infrared spectrum of the protostar must also present the spectral signatures found in the prototypical FUors (Connelley & Reipurth 2018). In the case where a protostar shows these signatures and the photometric outburst was not detected, the source is considered “FUor-like.” And if a Class I protostar shows an outburst and no or a minimal number of spectral signatures, it is considered “peculiar.”

Outflows play an important role in the star formation process as they remove angular momentum from the accretion disk,

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), Konkoly-Thege Miklós út 15–17, 1121 Budapest, Hungary; cruzsaenz.fernando@csfk.org</footnote>

<footnote n="2">² CSFK, MTA Centre of Excellence, Konkoly Thege Miklós út 15–17, 1121, Budapest, Hungary</footnote>

<footnote n="3">³ Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="4">⁴ ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="5">⁵ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>

<footnote n="6">⁶ IRAM, 300 Rue de la piscine, F-38406 Saint-Martin-d'Hères, France</footnote>

<footnote n="7">⁷ School of Physics and Astronomy, Cardiff University, Queen's Buildings, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote n="8">⁸ Max Planck Institute for Radioastronomy, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>
inject mass and energy into their surroundings, and clear material from the envelope (Arce & Sargent 2006). The circumstellar envelopes are the remains of the parent molecular cloud core that surround the protostar, and their properties (i.e., mass and extension) are deeply connected with how evolved a young star is, with younger objects having more massive and larger envelopes than their evolved counterparts (Andre & Montmerle 1994). Therefore, if the elevated accretion rates during the outbursts can inject more momentum into the envelopes via outflows, then these episodic events must play an important role in the evolution of their protostellar system. Indeed, it is expected that after an eruption, the inner circumstellar disk becomes depleted and will be replenished by the surrounding envelope (Vorobyov & Basu 2006) until the system can erupt again (Bell & Lin 1994; Takami et al. 2018). Eventually, the repetitive outbursts will clear out the envelope, and the young system will move to its next evolutionary phase, from Class I to Class II (Green et al. 2006; Quanz et al. 2007; Green et al. 2013).

Previous observations of CO rotational transitions have shown the presence of outflows in some known FUors: V1057 Cyg (Rodríguez et al. 1990), V1735 Cyg (Evans et al. 1994), L1551 IRS 5 (Wu et al. 2009, and references therein), V883 Ori (Ruíz-Rodríguez et al. 2017a), Reipurth 50 (Ruíz-Rodríguez et al. 2017b), FU Ori (Hales et al. 2015), V1647 Ori (Principe et al. 2018), V2775 Ori (Zurlo et al. 2017), V346 Nor (Kóspál et al. 2017b), and V900 Mon (Takami et al. 2019). In other cases, optical and near-infrared spectroscopy has shown an indication of high-velocity jets: Z CMa (Poetzel et al. 1989), V899 Mon (Ninan et al. 2015), iPTF 15afq (Hillenbrand 2019), and V346 Nor (Kóspál et al. 2020b).

Each of the aforementioned studies focused on a single FUor-type object or a few of them, preventing a statistical analysis of their properties. In this paper, we present a systematic study of the envelopes surrounding FUors, and we search for outflows among our full sample. We then compare our results with outflows found in young stellar objects (YSOs) that are currently quiescent and for which it is unknown whether they experienced an outburst or not. As outflows found at thousands of astronomical units are an indication of the accretion history of a protostar, this comparison allows us to examine how comparable the histories of the FUors are with those of the quiescent sample. The triggering mechanism behind the FUor-type outbursts is still not understood; however, it is possible that an examination of the differences between the two samples might hint that FUors are protostars with intrinsic differences that caused the outburst. Alternatively, it could show that the samples are similar; thus, we cannot rule out that quiescent sources experienced FUor-type outbursts in the past.

The structure of the paper is as follows. The observed sample is briefly introduced in Section 2, while in Section 3, we describe the observations and the data reduction. In Section 4, we present the distribution of the gas in the environment surrounding the FUors and the properties of the integrated line profiles. The main goal of this paper is to study the properties of the circumstellar gas; this analysis is found in Section 5, including the characterization of the molecular outflows, where detected. In Section 6, we compare the outflows found in the FUor sample with nonoutbursting sources and draw our conclusions. Finally, in Section 7, we summarize our work and present our main findings.

## 2. Sample

Our sample is composed of 20 eruptive young stars, including most of the known FUors accessible from the APEX site (Audard et al. 2014; Connelley & Reipurth 2018). We note that not all of the targets in our list are considered FUors, as some objects are cataloged as FUor-like objects. This subclassification is used when the photometric outburst was not detected but the near-infrared spectrum shows features similar to those of the prototypical FUors (e.g., BBW 76; Connelley & Reipurth 2018). The target list also includes eruptive young stars with peculiar accretion histories (e.g., V1647 Ori) and a massive star with a powerful accretion outburst (V723 Car). The non-FUor objects were included because of the sudden increases of their mass accretion rate: thus, the properties of their outflows could be affected. The first part of the sample, composed of eight targets, was analyzed by Kóspál et al. (2017a), where they found outflows in three objects: HBC 494, Haro 5a IRS, and V346 Nor. Here we will analyze the full sample, including a reprocessing of the target list presented in Kóspál et al. (2017a). The sample presented in this paper includes ~50% of the currently known FUors and FUor-like objects (Connelley & Reipurth 2018). The full target list is presented in Table 1.

## 3. Observations and Data Reduction

We carried out two programs with the FLASH⁺ receiver (Klein et al. 2014) at the APEX telescope (Güsten et al. 2006) to measure the ¹²CO (3–2), ¹³CO (3–2), and ¹²CO (4–3) lines toward our targets. Program 094.F-9508 was observed between 2014 August 23 and 28 and program 098.F-9505 between 2016 August 25 and 2016 September 10. Both programs used the same technical setup and reduction process. The lower-frequency channel was tuned to 344.2 GHz in USB to cover ¹³CO (3–2) at 330.588 GHz and ¹²CO (3–2) at 345.796 GHz. The higher-frequency channel was tuned to the ¹²CO (4–3) line at 461.041 GHz in USB. We used the XFFTS backends providing a nominal 38 kHz spectral resolution for the *J* = 3–2 lines and 76 kHz for the *J* = 4–3 line; these resulted in spectral resolutions of ~34.5, ~32.9, and ~49.4 m s⁻¹ for the ¹³CO (3–2), CO (3–2), and CO (4–3) lines, respectively. For each target, 120″ × 120″ on-the-fly maps were obtained at 6″ s⁻¹ using a relative reference off position 1000″ away in R.A.

We removed a first-order baseline from the spectra and calibrated the data using a main beam efficiency of 0.73 and 0.60 at 352 and 464 GHz, respectively, and the values were converted to janskys using 41 and 48 Jy K⁻¹ at 352 and 464 GHz, respectively. We calculated the noise levels of each CO line by first selecting the first and last 100 channels of each cube (individually confirmed to be free of line emission), calculating the noise levels for each FUor using these channels, and then calculating the median noise level of all FUors to obtain representative values. The rms noise levels at the native spectral resolution mentioned earlier are 3.6 Jy for ¹³CO (3–2), 3.7 Jy for ¹²CO (3–2), and 9.4 Jy for ¹²CO (4–3). The telescope’s half-power beamwidth is 19″/2, and it is 15″/3 at the corresponding frequencies. As mentioned earlier, the first half of the survey has already been published by Kóspál et al. (2017a), and here we use their calibrated data for our analyses.
<table><caption>**Table 1**
Objects Observed</caption><thead><tr><th>Name</th><th>Coordinates</th><th>Class<sup>a</sup></th><th>FUor Classification</th><th>v<sub>LSR</sub><sup>b</sup><br>(km s<sup>-1</sup>)</th><th>Distance<sup>c</sup><br>(pc)</th><th>L<sub>bol</sub><sup>d</sup><br>(L<sub>☉</sub>)</th></tr></thead><tbody><tr><td>L1551 IRS 5</td><td>04:31:34.07 +18:08:04.9</td><td>I</td><td>FUor-like</td><td>6.46</td><td>147</td><td>25</td></tr><tr><td>V582 Aur</td><td>05:25:51.97 +34:52:30.0</td><td>F</td><td>Bona fide FUor</td><td>-10.85</td><td>1320</td><td>146</td></tr><tr><td>Haro 5a IRS</td><td>05:35:26.75 -05:03:55.1</td><td>I</td><td>FUor-like</td><td>10.90</td><td>391</td><td>50</td></tr><tr><td>V883 Ori</td><td>05:38:18.09 -07:02:25.9</td><td>F</td><td>Bona fide FUor</td><td>4.10</td><td>388</td><td>400</td></tr><tr><td>Reipurth 50<sup>e</sup></td><td>05:40:27.45 -07:27:30.0</td><td>I</td><td>Peculiar</td><td>3.76</td><td>460</td><td>300</td></tr><tr><td>FU Ori</td><td>05:45:22.36 +09:04:12.2</td><td>I/II</td><td>Bona fide FUor</td><td>11.96</td><td>402</td><td>420</td></tr><tr><td>V1647 Ori</td><td>05:46:13.13 -00:06:04.8</td><td>I/II</td><td>Peculiar</td><td>10.06</td><td>388</td><td>39</td></tr><tr><td>V2775 Ori</td><td>05:42:48.48 -08:16:34.7</td><td>I</td><td>Bona fide FUor</td><td>3.08</td><td>428</td><td>25</td></tr><tr><td>V899 Mon</td><td>06:09:19.24 -06:41:55.8</td><td>F/II</td><td>Peculiar</td><td>9.57</td><td>785</td><td>419</td></tr><tr><td>AR 6A<sup>f</sup></td><td>06:40:59.30 +09:35:52.3</td><td>I</td><td>Peculiar</td><td>5.02</td><td>890</td><td>450</td></tr><tr><td>V900 Mon</td><td>06:57:22.22 -08:23:17.6</td><td>I</td><td>Bona fide FUor</td><td>13.77</td><td>1130</td><td>106</td></tr><tr><td>V960 Mon</td><td>06:59:31.58 -04:05:27.7</td><td>II</td><td>Bona fide FUor</td><td>23.81</td><td>2068</td><td></td></tr><tr><td>Z CMa</td><td>07:03:43.15 -11:33:06.2</td><td>I</td><td>FUor-like</td><td>13.91</td><td>1150</td><td>500</td></tr><tr><td>iPTF 15afq<sup>g</sup></td><td>07:09:21.39 -10:29:34.4</td><td>I/F</td><td>Peculiar</td><td>14.04</td><td>920</td><td>7.2</td></tr><tr><td>BBW 76<sup>e</sup></td><td>07:50:35.59 -33:06:23.9</td><td>I</td><td>FUor-like</td><td>17.64</td><td>1040</td><td>287</td></tr><tr><td>V723 Car</td><td>10:43:23.44 -59:33:55.3</td><td>I</td><td>Peculiar</td><td>-19.58</td><td>2500</td><td>4000</td></tr><tr><td>GM Cha</td><td>11:09:28.55 -76:33:28.1</td><td>I/II</td><td>Peculiar</td><td>4.86</td><td>160</td><td>1.5</td></tr><tr><td>V346 Nor</td><td>16:32:32.19 -44:55:30.7</td><td>0/I</td><td>Peculiar</td><td>-3.08</td><td>700</td><td>135</td></tr><tr><td>OO Ser</td><td>18:29:49.13 +01:16:20.6</td><td>I</td><td>Peculiar</td><td>8.36</td><td>311</td><td>31</td></tr><tr><td>HBC 687<sup>h</sup></td><td>19:29:00.87 +09:38:42.9</td><td>II</td><td>FUor-like</td><td>16.98</td><td>400</td><td>10</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Here “Class” refers to the Classification based on the shape of its SED (Lada 1987; Andre et al. 1993; Greene et al. 1994).

<sup>b</sup> See Section 4.2.

<sup>c</sup> See Section 5.1.

<sup>d</sup> Obtained from the literature (e.g., Audard et al. 2014; Connelley & Reipurth 2018, and references therein).

<sup>e</sup> Reipurth 50 and BBW 76 are labeled as HBC 494 and Bran 76 in Kóspál et al. (2017a).

<sup>f</sup> Also known as V912 Mon.

<sup>g</sup> Also known as Gaia19fct.

<sup>h</sup> Also known as Parsamian 21.

# 4. Results

## 4.1. Distribution of Gas

We constructed velocity-integrated emission maps (moment zero) for $^{12}$CO (3–2) using all channels in our data cubes. The resulting maps are presented in Figure 1. Some young eruptive stars are still deeply embedded; therefore, it is possible that the observed CO emission originates from the remaining material in their surrounding envelope. In order to verify that our CO detections come from the FUors, we compared our moment zero maps with the dust continuum emission. We searched for 250 µm continuum maps taken with Herschel/SPIRE maps in the Herschel Science Archive$^9$ and found data at an angular resolution of 17".6 for 18 sources. For the two remaining sources (V900 Mon and Z CMa), we searched the Canadian Astronomy Data Centre$^{10}$ for archival 850 µm observations taken with the James Clerk Maxwell Telescope$^{11}$ (JCMT) with an angular resolution of 14".5. In the cases of L1551 IRS 5, Haro 5a IRS, V883 Ori, Reipurth 50, V899 Mon, V960 Mon, Z CMa, V346 Nor, GM Cha, and HBC 687, the peaks of both the dust and gas emission are located at the position of the protostar. For five of our targets (V582 Aur, AR 6A, iPTF 15afq, V723 Car, and OO Ser), the brighter peaks of

both gas and dust emission are offset from the position of the protostar. The continuum peaks to these five sources are 27" to the southeast, 12" to the east, 28" to the northwest, 22" to the west, and 6" to the southwest, respectively. In the case of OO Ser, there is continuum emission toward the position of the FUor; however, the brightest peak is the one previously mentioned. We find that for FU Ori, V1647 Ori, V2775 Ori, and V900 Mon, the dust emission is located at the position of the protostar; however, the peak of the gas emission is offset. Finally, in the case of BBW 76, the dust emission peaks at the position of the source; however, the CO map shows that the emission is weak and extended.

## 4.2. Systemic Velocity

As shown below, to estimate the kinematic properties of the outflows, we need a reliable estimate of the systemic velocity for each target so that we can measure the velocity of the outflow relative to the protostar. The systemic velocities of our targets were estimated by fitting a Gaussian function to the line profile of $^{13}$CO, extracted using a circular aperture with a radius of 10,000 au, and using the center of the best-fitting Gaussian as the systemic velocity. Due to the proximity of L1551 IRS 5 and the field of view of our observations, we had to use a smaller aperture of 8000 au for this FUor. The line profiles and best-fit Gaussians are shown in Figure 2.

A number of sources have complicated line profiles that could not be fitted by a single Gaussian. Some of these show asymmetric dips around the peak of the emission (e.g., V1647 Ori, GM Cha, and OO Ser), which can be due to the

<footnote>$^9$ http://archives.esac.esa.int/hsa whsa/</footnote>

<footnote>$^{10}$ https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/</footnote>

<footnote>$^{11}$ The JCMT has historically been operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the National Research Council of Canada, and the Netherlands Organisation for Scientific Research.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Integrated intensity (moment zero) maps of our targets for the ¹²CO (3–2) line observed with APEX (orange contours). The moment zero maps were generated by integrating the full spectral cube in order to produce an unbiased map. The purple contours are the 250 µm continuum emission from Herschel for most of our targets; the two exceptions are V900 Mon and Z CMa, where we show the contours of the 850 µm continuum emission from the JCMT (see Section 4.1). The CO and dust contours are plotted with levels at 0.3, 0.4, ..., 0.9 of the peak intensity and are meant to be representative. The stars indicate the nominal positions of the protostars. The bars at the bottom of each panel represent 10,000 au.</figcaption></figure>

self-absorption of the envelope or the rotation of the gas. The former scenario is more likely based on the inspection of the channel maps of ¹³CO; thus, for these sources, we discarded the velocity range of the dip and fitted the Gaussian function using the remaining velocities. The line profiles of other objects show asymmetric shapes out to the wings of the line profiles (e.g., V582 Aur, Reipurth 50, and V346 Nor), an indication of multiple components (e.g., envelope, outflows, Keplerian disk, or unrelated gas in the same line of sight) showing emission at ¹³CO. For these objects, we fitted a combination of two or three Gaussian functions to the line profile and used the best-fit mean of the Gaussian with the highest amplitude to estimate the systemic velocity.

To verify our estimated systemic velocities, we searched for the velocity around which the line profile is most symmetrical. As expected, we found that the sources with asymmetrical line profiles have the largest differences, but still less than 1 km s⁻¹. For the symmetrical sources, the differences are less than 0.3 km s⁻¹. In addition, we examined the channel maps of each target to confirm our systemic velocity estimates.

As a final step, we compared our estimates with those from the literature. The differences between our estimates and those

obtained from previous observations are 0.04 km s⁻¹ for L1551 IRS 5 (Wu et al. 2009), 0.02 km s⁻¹ for V2775 Ori (Zurlo et al. 2017), 0.04 km s⁻¹ for GM Cha (Mottram et al. 2017), 0.20 km s⁻¹ for V883 Ori (Ruíz-Rodríguez et al. 2017a), 0.06 km s⁻¹ for V1647 Ori (Principe et al. 2018), 0.35 km s⁻¹ for V582 Aur (Ábrahám et al. 2018), 0.27 km s⁻¹ for V900 Mon (Takami et al. 2019), and 0.56 km s⁻¹ for FU Ori North (Pérez et al. 2020, who resolved the binary system with an angular resolution of 0.05 using the Atacama Large Millimeter/ submillimeter Array, ALMA). In the cases of Haro 5a IRS, AR 6A, BBW 76, OO Ser, and HBC 687, the estimates by Kóspál et al. (2017a) are in agreement within 0.30 km s⁻¹. For V899 Mon, V960 Mon, Z CMa, iPTF 15afq, and V723 Car, these are the first estimates of their systemic velocity. Two of our measurements deviate from those determined by interferometric observations of C¹⁸O: V346 Nor and Reipurth 50. In the case of the former, Kóspál et al. (2017b) found that the line profile peaks at −3.55 km s⁻¹, indicating a difference of 0.47 km s⁻¹ from our estimate, and in the case of the latter FUor, Ruíz-Rodríguez et al. (2017b) determined a systemic velocity of 4.6 km s⁻¹, a value 0.77 km s⁻¹ different from ours. It is likely that the larger
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Line profiles of $^{13}$CO (black lines) extracted using a circular aperture with a radius of 10,000 au to determine the systemic velocities. The blue lines indicate the best-fit Gaussian when using the full velocity range, and the green lines indicate where the fit was done without including velocities close to the peak. The gray horizontal line at 0 Jy indicates the range of velocities excluded from this second fit. The red lines show the best fit when using two or three Gaussians. The vertical dashed line indicates the systemic velocity of each FUor. In the cases of V883 Ori, V2775 Ori, and V723 Car, i.e., FUors with known emission from other sources in the same line of sight, we did not use additional Gaussians to fit the additional components because they can be easily separated from the single Gaussian fit. See Section 4.2 for details.</figcaption></figure>

differences are due to the interferometric observations resolving out emission from the extended envelopes. The final values for the systemic velocities are presented in Table 1.

### 4.3. Line Profiles

In order to examine the outflows using their line profiles, we must select apertures that cover the gas emission. We began by exploring the channel maps of the two $^{12}$CO transitions and checking which channels and regions show emission above the $3\sigma$ contour level. The channels with wide extended emission that showed little variation from channel to channel were considered as envelope emission. Then we inspected the blueshifted and redshifted channel maps for emission similar to what is found in outflows, i.e., emission whose redshifted channels are on the opposite side from the blueshifted channels with respect to the expected position of the star and that is generally more extended in the channels with velocities closer to the systemic velocity and more compact toward higher velocities. Finally, we created a polygon whose shape would cover this emission in both transitions. For the targets where

the CO emission does not follow the morphology described above, the spectra were extracted using a 10,000 au aperture centered on the nominal position of the protostar. The only exception is L1551 IRS 5, where we used a circular aperture with a radius of 8000 au, due to the proximity of this source (see below) and the size of our CO map. For each target, we used the same aperture in the three CO maps. The aperture used for each target can be seen in the channel maps in Appendix A, and the spectral line profiles integrated over these apertures for all three observed CO lines are presented in Figure 3.

The line profiles and channel maps show contamination caused by faint extended emission in four of the FUors: V582 Aur (~9 km s⁻¹), V883 Ori (~5.5 km s⁻¹), V2775 Ori (~5.9 km s⁻¹), and V723 Car (~24 km s⁻¹). The peaks in the V883 Ori profiles were reported by White et al. (2019), and the blueshifted broad feature in V582 Aur was discussed by Ábrahám et al. (2018).

Objects HBC 687, BBW 76, FU Ori, and V883 Ori show the narrowest lines in our sample. The line profiles of V1647 Ori, V900 Mon, and Z CMa are slightly wider and do not show
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The CO line profiles of our targets observed with APEX. The vertical dotted line is the systemic velocity. The vertical dashed lines are the range of velocities of the CO (3–2) outflows. The line profiles have been smoothed for presentation purposes.</figcaption></figure>

obvious indications of wings caused by high-velocity outflows. The remainder of the sources exhibit much wider profiles with a clear indication of line wings and possible outflows, mainly in the $^{12}$CO (4–3) and $^{12}$CO (3–2) transitions. The $^{12}$CO (4–3) line is the strongest line for most sources, except for AR 6A and V960 Mon, where both transitions are equally strong. Indeed, for most FUors, the ratio between line profiles, ($J=4-3$)/($J=3-2$), is <1.5 at the systemic velocity of each object. The two exceptions are FU Ori, where the $J=3-2$ transition almost reaches zero due to strong self-absorption, and V582 Aur, where there is a ratio of ~4.5. For the latter, this ratio suggests different excitation conditions, which can be explained by the intense radiation from two early B-type stars within 30 pc of V582 Aur that are exciting the region surrounding the FUor (Kun et al. 2017).

Some of our line profiles are different from those presented by Kóspál et al. (2017a). These discrepancies are because of differences in the distance to the FUors and the shape of the apertures. An example of the former is BBW 76, for which they used a distance 660 pc larger than ours (see below); thus, their aperture covered fewer pixels, causing a difference in the integrated flux of a factor of ~3. A similar scenario applies to AR 6A. Concerning the different shapes of the apertures,

Kóspál et al. (2017a) used a 10,000 au circular aperture for all targets, while we tailored the shape of our apertures. Haro 5a IRS, Reipurth 50, and V346 Nor are examples of this, where our apertures produced higher integrated fluxes by a factor of ~3.

Unsurprisingly, we find that the $^{13}$CO (3–2) transition produces the faintest line in all targets. Its line profiles are single-peaked for most FUors, with the maximum at velocities close to the systemic velocity (see below). Objects GM Cha, iPTF 15afq, and OO Ser are double-peaked with slightly less emission at the systemic velocity, a possible indication that $^{13}$CO (3–2) is optically thick at the line center.

# 5. Analysis

## 5.1. Distances

The estimation of gas masses is dependent on the distance to the target. For FU Ori, V899 Mon, AR 6A, V900 Mon, V960 Mon, and BBW 76, we used photogeometric distances from the Gaia Early Data Release 3 (Bailer-Jones et al. 2021). In the case of the more embedded objects (i.e., undetected by Gaia), L1551 IRS 5, V883 Ori, V1647 Ori, and V2775 Ori, we used the distances estimated from the distance to their molecular clouds
<table><caption>**Table 2**
Envelope Masses of the FUors Based on ¹³CO (3–2)</caption><thead><tr><th>Name</th><th>v<sub>min</sub><br>(km s<sup>-1</sup>)</th><th>v<sub>max</sub><br>(km s<sup>-1</sup>)</th><th>Int. Flux<br>(Jy km s<sup>-1</sup>)</th><th>M<sub>env</sub><br>(M<sub>☉</sub>)</th><th>10 µm Feature</th><th>Si Reference</th></tr></thead><tbody><tr><td>L1551 IRS 5</td><td>3.11</td><td>9.03</td><td>4920.70</td><td>0.193</td><td>Absorption</td><td>Quanz et al. (2007)</td></tr><tr><td>V582 Aur</td><td>-11.73</td><td>-8.34</td><td>43.35</td><td>0.137</td><td>Emission</td><td>Kóspál et al. (2020a)</td></tr><tr><td>Haro 5a IRS</td><td>8.92</td><td>13.21</td><td>3714.50</td><td>1.032</td><td>Absorption</td><td>Postel et al. (2019)</td></tr><tr><td>V883 Ori</td><td>3.00</td><td>4.80</td><td>506.04</td><td>0.138</td><td>Absorption</td><td>Quanz et al. (2007)</td></tr><tr><td>Reipurth 50</td><td>1.28</td><td>7.61</td><td>1567.74</td><td>0.603</td><td>Absorption</td><td>Quanz et al. (2007)</td></tr><tr><td>FU Ori</td><td>11.20</td><td>12.73</td><td>323.91</td><td>0.095</td><td>Emission</td><td>Quanz et al. (2007)</td></tr><tr><td>V1647 Ori</td><td>7.95</td><td>11.79</td><td>1903.05</td><td>0.521</td><td>Absorption</td><td>Quanz et al. (2007)</td></tr><tr><td>V2775 Ori</td><td>0.72</td><td>4.53</td><td>1363.13</td><td>0.454</td><td>Absorption</td><td>Kim et al. (2016)</td></tr><tr><td>V899 Mon</td><td>7.92</td><td>11.27</td><td>196.46</td><td>0.220</td><td>Emission</td><td>Kóspál et al. (2020a)</td></tr><tr><td>AR 6a</td><td>0.24</td><td>8.40</td><td>899.39</td><td>1.295</td><td>Unknown</td><td></td></tr><tr><td>V900 Mon</td><td>12.26</td><td>14.50</td><td>58.73</td><td>0.136</td><td>Emission</td><td>Kóspál et al. (2020a)</td></tr><tr><td>V960 Mon</td><td>21.82</td><td>25.70</td><td>80.29</td><td>0.624</td><td>Emission</td><td>Kóspál et al. (2020a)</td></tr><tr><td>Z CMa</td><td>11.96</td><td>15.73</td><td>499.94</td><td>1.201</td><td>Absorption</td><td>Quanz et al. (2007)</td></tr><tr><td>iPTF 15afq</td><td>12.27</td><td>14.63</td><td>73.41</td><td>0.113</td><td>Unknown</td><td></td></tr><tr><td>BBW 76</td><td>17.36</td><td>17.98</td><td>8.39</td><td>0.016</td><td>Emission</td><td>Quanz et al. (2007)</td></tr><tr><td>V723 Car</td><td>-21.70</td><td>-16.13</td><td>294.20</td><td>3.341</td><td>Absorption</td><td>Kóspál et al. (2020a)</td></tr><tr><td>GM Cha</td><td>3.31</td><td>6.39</td><td>4099.17</td><td>0.191</td><td>Absorption</td><td>Manoj et al. (2011)</td></tr><tr><td>V346 Nor</td><td>-6.60</td><td>-0.86</td><td>443.84</td><td>0.395</td><td>Absorption</td><td>Quanz et al. (2007)</td></tr><tr><td>OO Ser</td><td>5.22</td><td>12.38</td><td>3071.60</td><td>0.540</td><td>Absorption</td><td>Quanz et al. (2007)</td></tr><tr><td>HBC 687</td><td>16.71</td><td>17.40</td><td>34.74</td><td>0.010</td><td>Emission</td><td>Quanz et al. (2007)</td></tr></tbody></table>

Note. Here $v_{min}$ and $v_{max}$ indicate the velocity range used to integrate and calculate the line fluxes. The emission or absorption of the 10 µm silicate feature is indicated when known, as is the reference used for each target.

and the positions of the FUors within them (Connelley & Reipurth 2018, and references therein). For V582 Aur, we used the distance estimated by Kun et al. (2017) under the assumption that the FUor is related to the Aur OB1 association. We followed Tapia et al. (2015) and used the mean distance to the Great Carina Nebula (NGC 3372) for V723 Car. The distance to iPTF 15afq was estimated by Park et al. (2022) after comparing different distance estimates based on kinematics, Gaia parallax, and the distance to the CMa OB1 association, to which this object belongs. For Reipurth 50, Z CMa, GM Cha, V346 Nor, and OO Ser, we used distances compiled from the literature (Audard et al. 2014, and references therein). We compared our distances to those used by Kóspál et al. (2017a) and found that four FUors have different distance estimates: Haro 5a IRS, AR 6A, V900 Mon, and BBW 76. The differences between our and their estimates are -79, 90, 30, and -660 pc, respectively. If we had used the same apertures and velocity integration ranges as Kóspál et al. (2017a), these differences in distance would translate to a difference in mass of factors of 0.69, 1.24, 1.06, and 0.37, respectively.

## 5.2. Envelope Masses

We used the ¹³CO (3–2) emission to calculate the masses of the envelopes surrounding the FUors. To calculate the integrated fluxes, we used the line profiles defined in Section 4.3 and integrated the channels that had emission above $3\sigma$. We assumed local thermodynamical equilibrium and used an excitation temperature of 20 K, a ¹³CO/¹²CO abundance ratio of 69 (Wilson 1999), and a ¹²CO/H₂ abundance ratio of $10^{-4}$ (Bolatto et al. 2013). The velocity range used to calculate the line fluxes, the resulting line fluxes, and the envelope mass estimates are presented in Table 2. To test the impact of our choice on gas temperature, we did the calculations using 10 or 50 K and found that our estimated envelope masses would change by a factor of ~0.94 or ~2.55,

respectively. As we show later, some FUors have optically thick emission at velocities close to the systemic velocity; therefore, the estimated masses are lower limits for these sources.

## 5.3. Outflow Detection

Here we explain the process we followed to determine whether a FUor had an outflow detection. We inspected all of the sources in our sample, including the ones that Kóspál et al. (2017a) considered as not having an outflow.

As mentioned earlier, high-velocity wings in the line profiles of ¹²CO are a common indicator of outflows, and these are present in some of our FUors. However, this feature by itself is not enough. Thus, we examined the ¹²CO channel maps for each target (found in Appendix A) to verify the existence of the outflows via a visual inspection of the different distribution of the gas between the channels close to the systemic velocity and outward to higher velocities. The envelope emission dominates the velocities closest to the systemic, which are approximately the same velocities covered by the ¹³CO emission (Table 2), so we focused on velocities beyond these. If the emission is not detected at velocities higher than the ones overrun by the envelope, then we consider the FUor as not having an outflow. As outflows originate from the protostars, it is expected that at lower velocities (with respect to the systemic), the outflow is extended, and its position is closer to the FUor. So, in the case where there is emission in the maps beyond the envelope velocities, we checked the separation of this gas with respect the position of the FUor. Should the emission be close to the FUor, we consider it to be an outflow, and in the case where they are separated, we do not. When compared to the outflow detections of Kóspál et al. (2017a), our methodology resulted in almost the same detections and nondetections, with the only difference begin V900 Mon. In this case, the line profiles do
<table><caption>**Table 3**
Position Angles, Velocities, Extensions, and Dynamical Times of Outflows</caption><thead><tr><th rowspan="2">Target</th><th rowspan="2">Inc. (deg)</th><th rowspan="2">PA (deg)</th><th rowspan="2">Side</th><th colspan="5">CO (3–2)</th><th colspan="5">CO (4–3)</th></tr><tr><th>v<sub>in</sub> (km s<sup>-1</sup>)</th><th>v<sub>out</sub> (km s<sup>-1</sup>)</th><th>|v<sub>max</sub>| (km s<sup>-1</sup>)</th><th>R<sub>lobe</sub> (10<sup>3</sup> au)</th><th>τ<sub>d</sub> (10<sup>3</sup> yr)</th><th>v<sub>in</sub> (km s<sup>-1</sup>)</th><th>v<sub>out</sub> (km s<sup>-1</sup>)</th><th>|v<sub>max</sub>| (km s<sup>-1</sup>)</th><th>R<sub>lobe</sub> (10<sup>3</sup> au)</th><th>τ<sub>d</sub> (10<sup>3</sup> yr)</th></tr></thead><tbody><tr><td>L1551 IRS 5</td><td>70</td><td>65</td><td>Blue*</td><td>5.03</td><td>-0.27</td><td>6.73</td><td>11.7</td><td>8.2</td><td>5.64</td><td>0.93</td><td>5.54</td><td>11.7</td><td>10.0</td></tr><tr><td rowspan="2">Haro 5a IRS</td><td rowspan="2">50</td><td rowspan="2">70</td><td>Red*</td><td>7.57</td><td>12.80</td><td>6.33</td><td>11.7</td><td>8.7</td><td>7.67</td><td>12.24</td><td>5.77</td><td>11.7</td><td>9.6</td></tr><tr><td>Blue*</td><td>9.20</td><td>5.50</td><td>5.35</td><td>27.2</td><td>24.1</td><td>9.53</td><td>5.30</td><td>5.55</td><td>25.3</td><td>21.6</td></tr><tr><td rowspan="2">Reipurth 50</td><td rowspan="2">70</td><td rowspan="2">150</td><td>Red*</td><td>13.40</td><td>16.63</td><td>5.78</td><td>26.6</td><td>21.8</td><td>12.80</td><td>15.50</td><td>4.65</td><td>22.9</td><td>23.3</td></tr><tr><td>Blue*</td><td>1.03</td><td>-0.98</td><td>4.81</td><td>18.8</td><td>18.5</td><td>2.60</td><td>1.01</td><td>2.82</td><td>9.2</td><td>15.5</td></tr><tr><td rowspan="2">V2775 Ori</td><td rowspan="2">10</td><td rowspan="2">...</td><td>Red*</td><td>7.91</td><td>8.97</td><td>5.14</td><td>32.0</td><td>29.5</td><td>5.18</td><td>8.01</td><td>4.18</td><td>23.6</td><td>26.7</td></tr><tr><td>Blue</td><td>1.36</td><td>-5.58</td><td>8.66</td><td>6.3</td><td>3.4</td><td>1.61</td><td>-1.62</td><td>4.70</td><td>10.7</td><td>10.8</td></tr><tr><td rowspan="2">V899 Mon</td><td rowspan="2">50</td><td rowspan="2">60</td><td>Red*</td><td>8.60</td><td>6.92</td><td>2.62</td><td>46.1</td><td>83.3</td><td>8.93</td><td>7.99</td><td>1.55</td><td>24.4</td><td>74.6</td></tr><tr><td>Blue*</td><td>10.85</td><td>12.47</td><td>2.93</td><td>23.0</td><td>37.2</td><td>10.71</td><td>11.66</td><td>2.12</td><td>18.4</td><td>41.2</td></tr><tr><td rowspan="2">V900 Mon<sup>†</sup></td><td rowspan="2">30</td><td rowspan="2">80</td><td>Red</td><td>12.02</td><td>11.22</td><td>2.53</td><td>21.5</td><td>40.4</td><td>13.44</td><td>12.49</td><td>1.26</td><td>21.5</td><td>80.9</td></tr><tr><td>Blue</td><td>14.46</td><td>14.73</td><td>0.98</td><td>25.9</td><td>125.7</td><td>13.44</td><td>14.23</td><td>0.48</td><td>25.9</td><td>256.2</td></tr><tr><td rowspan="2">V960 Mon</td><td rowspan="2">10</td><td rowspan="2">...</td><td>Red</td><td>21.42</td><td>17.32</td><td>6.48</td><td>14.5</td><td>10.6</td><td>22.12</td><td>20.93</td><td>2.87</td><td>14.5</td><td>23.9</td></tr><tr><td>Blue</td><td>27.14</td><td>28.80</td><td>5.00</td><td>14.5</td><td>13.7</td><td>24.61</td><td>27.04</td><td>3.24</td><td>14.5</td><td>21.2</td></tr><tr><td rowspan="2">Z CMa</td><td rowspan="2">30</td><td rowspan="2">45</td><td>Red</td><td>12.13</td><td>10.71</td><td>3.15</td><td>29.9</td><td>45.0</td><td>12.44</td><td>11.25</td><td>2.61</td><td>21.6</td><td>39.2</td></tr><tr><td>Blue</td><td>16.00</td><td>19.34</td><td>5.48</td><td>47.0</td><td>40.6</td><td>15.56</td><td>19.18</td><td>5.32</td><td>40.8</td><td>36.3</td></tr><tr><td rowspan="2">iPTF 15afq<sup>†</sup></td><td rowspan="2">50</td><td rowspan="2">135</td><td>Red*</td><td>11.16</td><td>9.84</td><td>4.20</td><td>28.1</td><td>31.7</td><td>13.00</td><td>11.36</td><td>2.89</td><td>28.1</td><td>49.7</td></tr><tr><td>Blue*</td><td>17.97</td><td>20.65</td><td>6.61</td><td>60.4</td><td>43.3</td><td>13.00</td><td>16.22</td><td>1.97</td><td>60.4</td><td>131.3</td></tr><tr><td rowspan="2">GM Cha</td><td rowspan="2">70</td><td rowspan="2">100</td><td>Red</td><td>3.48</td><td>2.82</td><td>2.02</td><td>2.8</td><td>6.5</td><td>3.65</td><td>3.00</td><td>1.84</td><td>1.3</td><td>3.3</td></tr><tr><td>Blue</td><td>6.59</td><td>10.00</td><td>5.16</td><td>6.8</td><td>6.2</td><td>6.08</td><td>9.80</td><td>4.96</td><td>6.4</td><td>6.1</td></tr><tr><td rowspan="2">V346 Nor</td><td rowspan="2">30</td><td rowspan="2">45</td><td>Red</td><td>-5.41</td><td>-11.16</td><td>8.18</td><td>22.7</td><td>13.1</td><td>-4.58</td><td>-9.59</td><td>6.61</td><td>18.0</td><td>12.9</td></tr><tr><td>Blue*</td><td>-0.25</td><td>5.94</td><td>8.92</td><td>39.6</td><td>21.1</td><td>-1.26</td><td>4.38</td><td>7.36</td><td>30.0</td><td>19.3</td></tr></tbody></table>

Note. The inclination angles here are the values used for the inclination correction in Section 6.5; see text for details. The position angles were estimated by hand using the CO (3–2) integrated emission maps. The asterisk indicates the lobes that extend beyond the field of view of our observations; thus, their values of R<sub>lobe</sub> and τ<sub>d</sub> are lower limits. The † labels the two FUors with tentative outflow detections.

not show high-velocity wings and thus resulted in a nondetec-
tion for them.

After we detected an outflow, we determined the velocity ranges at which it was present in the blueshifted and redshifted sides. First, we found the velocity channel on which the envelope is not dominant and considered it as the “inner” velocity, v<sub>in</sub>, of the outflow. Next, we located the velocity channel where a 3σ detection was not found and considered it as the “outer” velocity, v<sub>out</sub>. We then calculated the maximum velocity of each lobe of the outflow with respect to the systemic velocity as v<sub>max</sub> = v<sub>out</sub> − v<sub>sys</sub>. The list of FUors with outflows and these three velocities is presented in Table 3. In Figure 3, we marked with vertical dotted lines the velocity ranges where outflows are detected. For two FUors, V900 Mon and iPTF 15afq, the <sup>12</sup>CO (4–3) emission at velocities close to the systemic does not appear to be dominated by the envelope; thus, we used the systemic as v<sub>in</sub> for both lobes.

Finally, we used these velocities to produce blueshifted and redshifted integrated emission maps of the *J* = 3–2 and *J* = 4–3 transitions of <sup>12</sup>CO, and the contour maps are presented in Figures 4 and 5. We used these maps to estimate the position angle of the outflow, reported in Table 3, and the extension of each lobe (see below).

## 5.4. Outflow Properties

One of the goals of this work is to compare the outflows emanating from eruptive young stars to those from quiescent YSOs. We carried out our calculations for the blueshifted and redshifted parts of the spectra separately, and here we describe how we carried out these calculations. The results for the

*J* = 3–2 and *J* = 4–3 transitions are presented in Tables 4 and 5, respectively.

### 5.4.1. Mass, Momenta, and Energy

The outflow masses were calculated assuming the wind emission is in local thermodynamical equilibrium with an excitation temperature of 75 K (van Kempen et al. 2009; Yıldız et al. 2015) and assuming a CO abundance of 10<sup>-4</sup> with respect to H<sub>2</sub> (Bolatto et al. 2013). We calculated the mass (M<sub>v</sub>) for each velocity channel (v) for all pixels above 3σ. Afterward, we calculated the momentum and kinematic energy for each channel with P<sub>v</sub> = M<sub>v</sub> × v and E<sub>v</sub> = 0.5M<sub>v</sub> × v<sup>2</sup>, respectively. Finally, we integrated the three properties over the same velocity range to obtain the total values (M<sub>of</sub>, P<sub>of</sub>, and E<sub>of</sub>) for the blueshifted and redshifted lobes of the outflow.

### 5.4.2. Force and Luminosity

The outflow force and luminosity are calculated as F<sub>of</sub> = P<sub>of</sub>/τ<sub>d</sub> and L<sub>of</sub> = E<sub>of</sub>/τ<sub>d</sub>, respectively, where τ<sub>d</sub> is the dynamical time of the outflow. The dynamical time is defined as τ<sub>d</sub> = R<sub>lobe</sub>/v<sub>max</sub>, where R<sub>lobe</sub> is the projected extension of the outflow lobe, and v<sub>max</sub> is the maximum velocity of the outflow.

For the projected extension of the outflows, we used the integrated emission maps (Figures 4 and 5) to measure the separation between the position of the star and the maximum length at which the outflow is above 3σ for each transition separately. However, the outflows around some of our targets extend beyond the field of view of our observations, so our extension measurements are only a lower limit. The lobes for which this is the case are indicated with an asterisk in Table 3. We estimated the maximum outflow velocity for each lobe
<figure><img src="image_7.png" /><figcaption>Figure 4. The red and blue contours show redshifted and blueshifted CO (J = 3–2) emission integrated in the velocity ranges indicated in Figure 3. The stars mark the stellar position as given in Table 1. The hatched circle in the bottom right panel is the APEX beam size, and the arrows indicate the orientation of the outflow. In the cases of V2775 Ori and V960 Mon, the outflow appears to be expanding in the direction of the line of sight. For L1551 IRS 5, the blueshifted contours are 3, 13, 24, 35, 45, 56, 67, and 78σ for σ = 0.38 K km s⁻¹, while the redshifted contours are 3, 18, 33, 48, 63, 78, 93, and 109σ for σ = 0.33 K km s⁻¹. For Haro 5a IRS, the blueshifted contours are 3, 7, 12, 16, 21, 25, 30, and 35σ for σ = 0.61 K km s⁻¹, while the redshifted contours are 3, 6, 10, 14, 1, 22, 26, and 30σ for σ = 0.48 K km s⁻¹. For Reipurth 50, the blueshifted contours are 3, 5, 7, 9, 11, 13, 15, and 17σ for σ = 0.48 K km s⁻¹, while the redshifted contours are 3, 6, 10, 14, 17, 21, 25, and 29σ for σ = 0.43 K km s⁻¹. For V2775 Ori, the blueshifted contours are 3, 12, 22, 31, 41, 50, 60, and 70σ for σ = 0.73 K km s⁻¹, while the redshifted contours are 3, 7, 12, 17, 21, 26, 31, and 36σ for σ = 0.68 K km s⁻¹. For V899 Monσ, the blueshifted contours are 3, 5, 7, 9, 11, 13, 15, and 17σ for σ = 0.52 K km s⁻¹, while the redshifted contours are 3, 4, 6, 8, 9, 11, 13, and 15σ for σ = 0.43 K km s⁻¹. For V900 Mon, the blueshifted contours are 3, 4, 5, 6, 7, 8, 9, and 10σ for σ = 0.62 K km s⁻¹, while the redshifted contours are 3, 4, 5, 6, 7, and 8σ for σ = 0.47 K km s⁻¹. For V960 Mon, the blueshifted contours are 3, 6, 9, 13, 16, 20, 23, and 27σ for σ = 0.44 K km s⁻¹, while the redshifted contours are 3, 7, 11, 15, 19, 23, 27, and 32σ for σ = 0.25 K km s⁻¹. For Z CMa, the blueshifted contours are 3, 4, 6, 8, 9, 11, 13, and 15σ for σ = 0.52 K km s⁻¹, while the redshifted contours are 3, 8, 13, 19, 24, 30, 35, and 41σ for σ = 0.39 K km s⁻¹. For iPTF 15afq, the blueshifted contours are 3, 5, 7, 9, 11, 13, 15, and 17σ for σ = 0.44 K km s⁻¹, while the redshifted contours are 3, 6, 10, 13, 17, 20, 24, and 28σ for σ = 0.32 K km s⁻¹. For GM Cha, the blueshifted contours are 3, 5, 7, 9, 11, 13, 15, and 18σ for σ = 0.43 K km s⁻¹, while the redshifted contours are 3, 9, 16, 23, 29, 36, 43, and 50σ for σ = 0.39 K km s⁻¹. For V346 Nor, the blueshifted contours are 3, 12, 22, 32, 41, 51, 61, and 71σ for σ = 0.52 K km s⁻¹, while the redshifted contours are 3, 11, 20, 29, 38, 47, 56, and 65σ for σ = 0.56 K km s⁻¹.</figcaption></figure>

independently by calculating the difference between the systemic velocity and the minimum/maximum velocity where there is blueshifted/redshifted emission. In both R<sub>lobe</sub> and v<sub>max</sub>, the sensitivity and spatial resolution of the observations directly affect their measured values.

### 5.4.3. Caveats

Due to the nature of our observations and methodology, it is important to understand the limitations of our estimated values.

First, large uncertainties exist in the estimations of the outflow masses due to the presence of envelopes. This surrounding material dominates emission at velocities close to the systemic velocity; thus, we have calculated the outflow properties using only channels where the outflow is the predominant flux contributor. Therefore, we are knowingly underestimating the outflow masses by not integrating the emission at low velocities to prevent this contamination. However, it is likely that some of the envelope emission is still included in our calculations. Indeed, in a couple of cases (V883 Ori and V1647 Ori), we are not able to separate their known outflows due to their emission being at velocities comparable to those of the surrounding cloud.

Second, the sensitivity of the observations puts strong constraints on the maximum velocities where outflows are detected. For example, in the case of L1551 IRS 5, Yıldız et al. (2015) reported a sum of maximum velocities of 21.5 km s⁻¹, while we obtained ~13.1 km s⁻¹ (see Table 3). Their JCMT observations had typical noise values of ~0.1 K, while our APEX observations have a noise value of ~0.38 K.

Finally, the extension of the outflow can reach beyond the field of view of our observations. An extreme example is L1551 IRS 5, whose outflow extends out to ~20' (e.g., Stojimirović et al. 2006), and our field of view is less than 1' (see Figures 4 and 5).

Therefore, both the outflow masses and their dynamical ages should be considered as lower limits, while the outflow forces and luminosities must be considered as highly uncertain.

### 5.5. Optical Depth Correction

We calculated the outflow parameters assuming the ¹²CO lines are optically thin; however, this isotopologue is typically optically thick. One way to correct for this optical depth issue is by using the ¹³CO emission, under the assumption that that isotopologue is optically thin, to correct the fluxes of ¹²CO. In this section, we present our methodology to correct the emission of the J = 3–2 transition of ¹²CO.

This correction was done following the procedure presented by Dunham et al. (2014). We assumed that both CO isotopologues are in local thermodynamical equilibrium at the same excitation temperature and with identical beam filling factors. Under these conditions, the brightness temperature ratio between the two isotopologues is given by

$$ \frac{T_{\text{mb,12}}}{T_{\text{mb,13}}} = \frac{1 - e^{-\tau_{12}}}{1 - e^{-\tau_{13}}}, \qquad (1) $$

where T<sub>mb,12</sub> and T<sub>mb,13</sub> are the brightness temperatures of ¹²CO and ¹³CO, respectively, and τ<sub>12</sub> and τ<sub>13</sub> are their respective
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Similar to Figure 4 but for $^{12}$CO (4–3). In the case of V900 Mon, the emission at this transition is not significant enough to be seen in this map, and the emission detected is in the outer parts of the field of view. In the iPTF 15afq map, the emission of the outflow is weak compared to the surrounding gas. For L1551 IRS 5, the blueshifted contours are 4, 7, 11, 15, 19, 23, 27, and 31σ for σ = 0.96 K km s⁻¹, while the redshifted contours are 4, 8, 12, 17, 21, 26, 30, and 35σ for σ = 0.84 K km s⁻¹. For Haro 5a IRS, the blueshifted contours are 4, 5, 7, 8, 10, 11, 13, and 15σ for σ = 2.09 K km s⁻¹, while the redshifted contours are 4, 5, 6, 8, 9, 11, 12, and 14σ for σ = 1.65 K km s⁻¹. For Reipurth 50, the blueshifted contours are 5, 6, 7, 8, and 9σ for σ = 1.22 K km s⁻¹, while the redshifted contours are 5, 7, 10, 13, 16, 19, 22, and 25σ for σ = 1.13 K km s⁻¹. For V2775 Ori, the blueshifted contours are 5, 8, 11, 15, 18, 22, 25, and 29σ for σ = 1.72 K km s⁻¹, while the redshifted contours are 5, 6, 7, 8, 9, 10, 11, and 12σ for σ = 1.62 K km s⁻¹. For V899 Mon, the blueshifted contours are 3, 4, 5, 6, 7, 8, and 9σ for σ = 1.50 K km s⁻¹, while the redshifted contours are 3, 4, 5, 6, 7, and 8σ for σ = 1.22 K km s⁻¹. For V900 Mon, the blueshifted contours are 4, 5, 6, 7, and 8σ for σ = 2.16 K km s⁻¹, while the redshifted contours are 4, 5, and 6σ for σ = 1.64 K km s⁻¹. For V960 Mon, the blueshifted contours are 4, 5, 6, and 7σ for σ = 1.38 K km s⁻¹, while the redshifted contours are 4, 6, 9, 12, 14, 17, 20, and 23σ for σ = 0.82 K km s⁻¹. For ZCMa, the blueshifted contours are 4, 7, 10, 13, 16, 19, 22, and 25σ for σ = 0.94 K km s⁻¹. For iPTF 15afq, the blueshifted contours are 4, 5, 6, 8, 9, 11, 12, and 14σ for σ = 1.15 K km s⁻¹, while the redshifted contours are 4, 5, 7, 9, 10, 12, 14, and 16σ for σ = 0.88 K km s⁻¹. For GM Cha, the blueshifted contours are 3, 4, 5, 6, and 7σ for σ = 1.21 K km s⁻¹, while the redshifted contours are 4, 6, 9, 12, 14, 17, 20, and 23σ for σ = 1.09 K km s⁻¹. For V346 Nor, the blueshifted contours are 4, 8, 13, 18, 23, 28, 33, and 38σ for σ = 1.66 K km s⁻¹, while the redshifted contours are 4, 8, 12, 16, 20, 24, 28, and 33σ for σ = 1.77 K km s⁻¹.</figcaption></figure>

<table><caption>**Table 4**
Outflow Properties from $^{12}$CO (3–2) Observations, Assuming They Are Optically Thin</caption><thead><tr><th>Target</th><th>Side</th><th>M<sub>of</sub><br>(M<sub>☉</sub>)</th><th>P<sub>of</sub><br>(M<sub>☉</sub> km s<sup>-1</sup>)</th><th>E<sub>of</sub><br>(erg)</th><th>F<sub>of</sub><br>(M<sub>☉</sub> yr<sup>-1</sup> km s<sup>-1</sup>)</th><th>L<sub>of</sub><br>(L<sub>☉</sub>)</th></tr></thead><tbody><tr><td rowspan="2">L1551 IRS 5</td><td>Blue</td><td>6.3 × 10<sup>-3</sup></td><td>1.9 × 10<sup>-2</sup></td><td>6.8 × 10<sup>41</sup></td><td>2.4 × 10<sup>-6</sup></td><td>6.9 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>6.8 × 10<sup>-3</sup></td><td>1.7 × 10<sup>-2</sup></td><td>4.9 × 10<sup>41</sup></td><td>1.9 × 10<sup>-6</sup></td><td>4.7 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">Haro 5a IRS</td><td>Blue</td><td>1.6 × 10<sup>-2</sup></td><td>3.6 × 10<sup>-2</sup></td><td>8.7 × 10<sup>41</sup></td><td>1.5 × 10<sup>-6</sup></td><td>3.0 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>8.6 × 10<sup>-3</sup></td><td>2.7 × 10<sup>-2</sup></td><td>8.6 × 10<sup>41</sup></td><td>1.2 × 10<sup>-6</sup></td><td>3.2 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">Reipurth 50</td><td>Blue</td><td>3.5 × 10<sup>-3</sup></td><td>1.1 × 10<sup>-2</sup></td><td>3.6 × 10<sup>41</sup></td><td>6.0 × 10<sup>-7</sup></td><td>1.6 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>7.2 × 10<sup>-3</sup></td><td>3.3 × 10<sup>-2</sup></td><td>1.5 × 10<sup>42</sup></td><td>1.1 × 10<sup>-6</sup></td><td>4.4 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">V2775 Ori</td><td>Blue</td><td>4.6 × 10<sup>-2</sup></td><td>1.7 × 10<sup>-1</sup></td><td>7.4 × 10<sup>42</sup></td><td>4.9 × 10<sup>-5</sup></td><td>1.8 × 10<sup>-2</sup></td></tr><tr><td>Red</td><td>3.5 × 10<sup>-2</sup></td><td>9.5 × 10<sup>-2</sup></td><td>2.8 × 10<sup>42</sup></td><td>2.2 × 10<sup>-5</sup></td><td>5.3 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">V899 Mon</td><td>Blue</td><td>1.3 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-2</sup></td><td>3.2 × 10<sup>41</sup></td><td>2.4 × 10<sup>-7</sup></td><td>3.2 × 10<sup>-5</sup></td></tr><tr><td>Red</td><td>2.2 × 10<sup>-2</sup></td><td>3.8 × 10<sup>-2</sup></td><td>6.8 × 10<sup>41</sup></td><td>1.0 × 10<sup>-6</sup></td><td>1.5 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">V900 Mon<sup>†</sup></td><td>Blue</td><td>6.0 × 10<sup>-3</sup></td><td>1.2 × 10<sup>-2</sup></td><td>2.5 × 10<sup>41</sup></td><td>3.1 × 10<sup>-7</sup></td><td>5.3 × 10<sup>-5</sup></td></tr><tr><td>Red</td><td>1.3 × 10<sup>-2</sup></td><td>1.0 × 10<sup>-2</sup></td><td>8.0 × 10<sup>40</sup></td><td>7.9 × 10<sup>-8</sup></td><td>5.2 × 10<sup>-6</sup></td></tr><tr><td rowspan="2">V960 Mon</td><td>Blue</td><td>1.1 × 10<sup>-1</sup></td><td>4.2 × 10<sup>-1</sup></td><td>1.7 × 10<sup>43</sup></td><td>3.9 × 10<sup>-5</sup></td><td>1.4 × 10<sup>-2</sup></td></tr><tr><td>Red</td><td>5.9 × 10<sup>-2</sup></td><td>2.3 × 10<sup>-1</sup></td><td>8.7 × 10<sup>42</sup></td><td>1.6 × 10<sup>-5</sup></td><td>5.2 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">Z CMa</td><td>Blue</td><td>1.9 × 10<sup>-2</sup></td><td>4.2 × 10<sup>-2</sup></td><td>9.3 × 10<sup>41</sup></td><td>9.4 × 10<sup>-7</sup></td><td>1.7 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>4.7 × 10<sup>-2</sup></td><td>1.5 × 10<sup>-1</sup></td><td>4.8 × 10<sup>42</sup></td><td>3.5 × 10<sup>-6</sup></td><td>9.7 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">iPTF 15afq<sup>†</sup></td><td>Blue</td><td>1.8 × 10<sup>-2</sup></td><td>5.8 × 10<sup>-2</sup></td><td>1.9 × 10<sup>42</sup></td><td>1.8 × 10<sup>-6</sup></td><td>5.0 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>3.8 × 10<sup>-2</sup></td><td>1.8 × 10<sup>-1</sup></td><td>8.9 × 10<sup>42</sup></td><td>4.2 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">GM Cha</td><td>Blue</td><td>7.4 × 10<sup>-4</sup></td><td>1.2 × 10<sup>-3</sup></td><td>1.8 × 10<sup>40</sup></td><td>1.8 × 10<sup>-7</sup></td><td>2.3 × 10<sup>-5</sup></td></tr><tr><td>Red</td><td>1.2 × 10<sup>-3</sup></td><td>3.3 × 10<sup>-3</sup></td><td>9.8 × 10<sup>40</sup></td><td>5.3 × 10<sup>-7</sup></td><td>1.3 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">V346 Nor</td><td>Blue</td><td>1.7 × 10<sup>-2</sup></td><td>6.4 × 10<sup>-2</sup></td><td>2.7 × 10<sup>42</sup></td><td>4.8 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-3</sup></td></tr><tr><td>Red</td><td>4.6 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-1</sup></td><td>9.8 × 10<sup>42</sup></td><td>9.8 × 10<sup>-6</sup></td><td>3.9 × 10<sup>-3</sup></td></tr></tbody></table>

**Note.** The † labels the two FUors with tentative outflow detections.
<table><caption>**Table 5**
Outflow Properties Estimated from the ¹²CO (4–3) Observations</caption><thead><tr><th>Target</th><th>Side</th><th><i>M</i><sub>of</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>P</i><sub>of</sub><br>(<i>M</i><sub>☉</sub> km s<sup>-1</sup>)</th><th><i>E</i><sub>of</sub><br>(erg)</th><th><i>F</i><sub>of</sub><br>(<i>M</i><sub>☉</sub> yr<sup>-1</sup> km s<sup>-1</sup>)</th><th><i>L</i><sub>of</sub><br>(<i>L</i><sub>☉</sub>)</th></tr></thead><tbody><tr><td rowspan="2">L1551 IRS 5</td><td>Blue</td><td>3.8 × 10<sup>-3</sup></td><td>8.4 × 10<sup>-3</sup></td><td>2.3 × 10<sup>41</sup></td><td>8.4 × 10<sup>-7</sup></td><td>1.9 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>2.9 × 10<sup>-3</sup></td><td>6.9 × 10<sup>-3</sup></td><td>1.9 × 10<sup>41</sup></td><td>7.2 × 10<sup>-7</sup></td><td>1.6 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">Haro 5a IRS</td><td>Blue</td><td>1.8 × 10<sup>-2</sup></td><td>3.3 × 10<sup>-2</sup></td><td>6.7 × 10<sup>41</sup></td><td>1.6 × 10<sup>-6</sup></td><td>2.6 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>1.3 × 10<sup>-2</sup></td><td>2.9 × 10<sup>-2</sup></td><td>6.9 × 10<sup>41</sup></td><td>1.2 × 10<sup>-6</sup></td><td>2.4 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">Reipurth 50</td><td>Blue</td><td>4.3 × 10<sup>-3</sup></td><td>7.3 × 10<sup>-3</sup></td><td>1.3 × 10<sup>41</sup></td><td>4.6 × 10<sup>-7</sup></td><td>6.8 × 10<sup>-5</sup></td></tr><tr><td>Red</td><td>1.9 × 10<sup>-2</sup></td><td>4.5 × 10<sup>-2</sup></td><td>1.2 × 10<sup>42</sup></td><td>1.7 × 10<sup>-6</sup></td><td>3.7 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">V2775 Ori</td><td>Blue</td><td>9.9 × 10<sup>-3</sup></td><td>2.4 × 10<sup>-2</sup></td><td>6.4 × 10<sup>41</sup></td><td>2.2 × 10<sup>-6</sup></td><td>4.9 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>1.0 × 10<sup>-2</sup></td><td>2.3 × 10<sup>-2</sup></td><td>5.7 × 10<sup>41</sup></td><td>2.7 × 10<sup>-6</sup></td><td>5.5 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">V899 Mon</td><td>Blue</td><td>3.2 × 10<sup>-3</sup></td><td>3.3 × 10<sup>-3</sup></td><td>3.7 × 10<sup>40</sup></td><td>4.5 × 10<sup>-8</sup></td><td>4.2 × 10<sup>-6</sup></td></tr><tr><td>Red</td><td>2.9 × 10<sup>-3</sup></td><td>4.5 × 10<sup>-3</sup></td><td>7.3 × 10<sup>40</sup></td><td>1.1 × 10<sup>-7</sup></td><td>1.4 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">V900 Mon<sup>†</sup></td><td>Blue</td><td>1.8 × 10<sup>-2</sup></td><td>1.3 × 10<sup>-2</sup></td><td>1.1 × 10<sup>41</sup></td><td>1.7 × 10<sup>-7</sup></td><td>1.2 × 10<sup>-5</sup></td></tr><tr><td>Red</td><td>1.3 × 10<sup>-2</sup></td><td>2.8 × 10<sup>-3</sup></td><td>8.1 × 10<sup>39</sup></td><td>1.1 × 10<sup>-8</sup></td><td>2.5 × 10<sup>-7</sup></td></tr><tr><td rowspan="2">V960 Mon</td><td>Blue</td><td>4.0 × 10<sup>-2</sup></td><td>8.2 × 10<sup>-2</sup></td><td>1.7 × 10<sup>42</sup></td><td>3.4 × 10<sup>-6</sup></td><td>6.0 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>1.7 × 10<sup>-1</sup></td><td>2.6 × 10<sup>-1</sup></td><td>4.5 × 10<sup>42</sup></td><td>1.2 × 10<sup>-5</sup></td><td>1.7 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">Z CMa</td><td>Blue</td><td>1.6 × 10<sup>-2</sup></td><td>2.8 × 10<sup>-2</sup></td><td>4.9 × 10<sup>41</sup></td><td>7.2 × 10<sup>-7</sup></td><td>1.1 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>4.3 × 10<sup>-2</sup></td><td>1.1 × 10<sup>-1</sup></td><td>3.2 × 10<sup>42</sup></td><td>3.0 × 10<sup>-6</sup></td><td>7.2 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">iPTF 15afq<sup>†</sup></td><td>Blue</td><td>2.2 × 10<sup>-2</sup></td><td>3.6 × 10<sup>-2</sup></td><td>6.3 × 10<sup>41</sup></td><td>7.2 × 10<sup>-7</sup></td><td>1.1 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>2.6 × 10<sup>-2</sup></td><td>2.2 × 10<sup>-2</sup></td><td>2.7 × 10<sup>41</sup></td><td>1.7 × 10<sup>-7</sup></td><td>1.7 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">GM Cha</td><td>Blue</td><td>3.8 × 10<sup>-4</sup></td><td>5.3 × 10<sup>-4</sup></td><td>7.6 × 10<sup>39</sup></td><td>1.6 × 10<sup>-7</sup></td><td>1.9 × 10<sup>-5</sup></td></tr><tr><td>Red</td><td>9.1 × 10<sup>-4</sup></td><td>2.1 × 10<sup>-3</sup></td><td>5.6 × 10<sup>40</sup></td><td>3.4 × 10<sup>-7</sup></td><td>7.5 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">V346 Nor</td><td>Blue</td><td>1.2 × 10<sup>-2</sup></td><td>3.7 × 10<sup>-2</sup></td><td>1.3 × 10<sup>42</sup></td><td>2.8 × 10<sup>-6</sup></td><td>8.4 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>2.5 × 10<sup>-2</sup></td><td>8.5 × 10<sup>-2</sup></td><td>3.3 × 10<sup>42</sup></td><td>4.5 × 10<sup>-6</sup></td><td>1.4 × 10<sup>-3</sup></td></tr></tbody></table>

Note. The † labels the two FUors with tentative outflow detections.

opacities. Assuming that ¹³CO is optically thin, Equation (1) can be rewritten as

$$ \frac{T_{\text{mb,12}}}{T_{\text{mb,13}}} = \frac{[^{12}\text{CO}]}{[^{13}\text{CO}]} \frac{1 - e^{-\tau_{12}}}{\tau_{12}}, \qquad (2) $$

where $[^{12}\text{CO}]/[^{13}\text{CO}]$ is the abundance ratio, for which we use a value of 69 (Wilson 1999).

We began the estimation of the correction factor $(1 - \exp(-\tau_{12}))/\tau_{12}$ by calculating $T_{\text{mb,12}}/T_{\text{mb,13}}$ for each channel where both isotopologues were detected above $6\sigma$. In some low-velocity channels for a few FUors, the ¹³CO appears to be optically thick; therefore, we dropped these points from the fitting.

We then fitted a parabola,

$$ \frac{T_{\text{mb,12}}}{T_{\text{mb,13}}} = A + B(v - v_{\text{sys}}) + C(v - v_{\text{sys}})^2, \qquad (3) $$

which will allow us to correct for the velocity channels where the ¹³CO emission was not detected. We fixed $B=0$ to keep the parabola symmetric with respect to the systemic velocity and prevent overcorrecting one side of the outflow. Finally, the correction factor selected for each channel was the lower value between the fitted parabola and the expected abundance ratio of 69. The plots and values of the fitted parabolas for each target are presented in Appendix B. We note that in the case of iPTF 15afq, due to the complex emission of ¹³CO, we only used blueshifted points to fit the parabola (see Appendix B).

## 6. Discussion

Here we present our discussion about the envelope masses and their relationship with the FUor evolutionary scheme presented by Quanz et al. (2007). Then we discuss the FUors for which we detected outflows, and we comment on the

sources for which an outflow was not detected. Finally, we make a statistical comparison between the properties of the outflows in our FUor sample and those from other works in the literature focused on quiescent protostars.

### 6.1. Envelope Masses

Quanz et al. (2007) targeted 14 FUor-type objects and obtained mid-infrared spectra. They found that the silicate feature at 10 µm could be present in either absorption or emission and suggested that when the feature is in absorption, it is an indication of a higher content of mass in the envelope surrounding the FUor and thus of the object being younger. In Figure 6, we compare our estimations of envelope masses to the emission/absorption of the silicate feature based on the references for each object listed in Table 2.

We expanded on the work presented by Kóspál et al. (2017b), who analyzed the first half of the FUor sample, and we found that the FUors with the least massive envelopes show the silicate feature in emission, while those with more massive envelopes show it in absorption. We found two exceptions to this trend: V899 Mon and V960 Mon. For the latter, there could be two explanations. As mentioned below, there are three YSOs inside the beam of our observations; thus, we could be significantly overestimating the amount of material in the line of sight to this FUor. Alternatively, if we consider that the outflow is indeed driven by the FUor then the moment zero maps indicate that the direction of the outflow is aligned with the line of sight; therefore, it could be that the outflow has already cleared the line of sight to the FUor, allowing the detection of the silicate feature in emission while maintaining a high envelope mass. It is harder to explain the case of V899 Mon, as our observations indicate that the direction of the outflow is perpendicular to the line of sight. Under the assumption that the outflows are perpendicular to the
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison between envelope masses and the emission/absorption of the silicate feature.</figcaption></figure>

inclination of the disks, we tried to verify the geometry of the systems using ALMA continuum observations (Kóspál et al. 2021). However, both disks were barely resolved; thus, the inclinations of the uncertainties are large enough to allow the scenarios of almost edge-on and almost face-on geometries. Observations with higher angular resolution and sensitivity are needed to determine the geometry of these systems and understand this discrepancy between envelope mass and the silicate feature.

The transition between absorption and emission appears to occur between 0.1 and $0.2 M_{\odot}$. Indeed, V900 Mon, V582 Aur, and V883 Ori have comparable envelope masses, with only the latter FUor having the silicate feature in absorption. Here it is not clear if the geometry of the system could explain this difference. The inclination of V582 Aur is unknown because continuum observations have not resolved the disk (Ábrahám et al. 2018), and the latter two FUors have comparable inclination angles (Cieza et al. 2018; Kóspál et al. 2021).

We do not have 10 µm data for two sources: AR 6A and iPTF 15afq. Based on its high envelope mass, we could expect the silicate feature around AR 6A to be in absorption. However, the peak of its CO emission is off-center (Figure 1), so the direct line of sight to our target could have less material and show the feature in emission. In the case of iPTF 15afq, the peak of the CO is also slightly off-center; however, its mass envelope falls in the intermediate range of masses, so we expect this to depend on the geometry of the system.

This suggests that HBC 687 is the most evolved FUor in our sample. The case for the least evolved FUor is less clear, as V723 Car is a massive young star and thus this evolutionary trend might not apply to it, and ZCMa is a binary with one of its stars being an intermediate-mass star (Koresko et al. 1991). Therefore, we consider Haro 5a IRS to be the youngest FUor in our sample, as it is the one with the most massive envelope with the silicate feature in absorption.

## 6.2. FUors with Outflows

L1551 IRS 5—This Class I protostar was among the first detections of bipolar outflows from YSOs (Snell et al. 1980). Later observations recovered the blueshifted and redshifted

lobes of the bipolar outflow in CO (2–1) (Moriarty-Schieven et al. 2006; Wu et al. 2009) and CO (3–2) (Yıldız et al. 2015). Based on our maps, the molecular outflows have the same geometry as seen in those previous works (e.g., Wu et al. 2009), and the position angle of the outflow (~45°) is almost perpendicular to the position angle of the circumstellar disks in the system (~160°; Lim et al. 2016; Cruz-Sáenz de Miera et al. 2019). Comparing our estimated outflow properties to those calculated by Yıldız et al. (2015), we find that our mass estimate is in agreement with their result, while our force and luminosity are a factor of ~6 lower than theirs, even when taking into account the inclination correction factor the authors applied. However, this difference can be due to the higher sensitivity (their $v_{max}$ is higher for both lobes) and the larger field of view of their observations.

Haro 5a IRS—This Class I protostar is located in the Orion star-forming region and was identified as a FUor-like object by Reipurth et al. (2012). Previous CO observations of the source revealed its outflow (Takahashi et al. 2006, 2008) with the same geometry as what we detected, including the slight overlap between the redshifted and blueshifted emission. Kóspál et al. (2017a) presented the $J=4–3$ and $J=3–2$ $^{12}$CO and $J=3–2$ $^{13}$CO observations of this FUor and found narrow outflow in an almost east-west direction. Our analysis, based on the same observations as them, recovered the same morphology. Their estimates for outflow masses are higher than ours by less than a factor of 2, which can be explained by the difference in distances, excitation temperatures, and the velocity ranges used in the calculation of the outflow properties. Tobin et al. (2020) and Kóspál et al. (2021) presented continuum observations at millimeter wavelengths with data from ALMA and the Very Large Array, and both reported that this FUor is also a proto-binary star.

Reipurth 50—This Class I protostar is also referred to as HBC 494. Ruíz-Rodríguez et al. (2017b) presented high angular resolution observations with ALMA in which they traced the emission from the outflow in $^{12}$CO ($J=2–1$) and the envelope emission with the same transition of $^{13}$CO and $^{18}$O. The extension of the $J=2–1$ outflow obtained with ALMA is smaller than the size of our beam. This can be explained by the maximum recoverable scale of their ALMA configuration (11")., which is comparable to the extended emission seen in their channel maps (e.g., the 6 km s$^{-1}$ channel in their Figure 3); thus, it is likely that they have resolved out most of the extended emission of the outflow. Indeed, most of the emission they recovered with ALMA originates from the dense cavity walls of the bipolar outflow. Nevertheless, the position angle obtained from the high-resolution interferometric observations (~145°) is comparable to our estimation of the position angle (150°). The outflow mass from the ALMA observations, calculated assuming an excitation temperature of 50 K, is a factor of 60 higher than the one we determined using the $J=3–2$ transition. If we adjust for the higher temperature used in our calculations (see Section 6.6), the mass estimated from ALMA measurements is still a factor of 50 higher. This difference hints that most of the mass of the outflow of Reipurth 50 is located in the narrow cavity walls, which are severely diluted by our single-dish beam.

V2775 Ori—The first detection of a molecular outflow on this object was done in the $J=2–1$ transitions of $^{12}$CO, $^{13}$CO, and $^{18}$O with ALMA (Zurlo et al. 2017). The authors found that the system is almost face-on with an inclination angle of
~14°. Our observations recovered a similar orientation of the outflow (see Appendix A). In addition, we found significant extended emission at both the systemic velocity and redshifted velocities (+3 km s⁻¹; see also Figure 3). Zurlo et al. (2017) reported different velocity ranges for ¹²CO and C¹⁸O (see their Table 2). The velocities of the ¹²CO match those of the redshifted excess emission (peaking at ~6 km s⁻¹; see Figure 3), and the velocities of C¹⁸O match those of the systemic emission we report in Table 1. The redshifted cloud emission appears at velocities where the outflow is still detected. Therefore, for all of the analyses of the outflow in this FUor, we removed the cloud's contamination. Similar to the case of Reipurth 50, the difference in beam sizes and sensitivities complicates the comparison between our estimated physical properties and those from Zurlo et al. (2017). However, we found that the masses from the *J* = 3–2 transition are higher by a factor of ~8 than those estimated from the ALMA observations, which is even more surprising due to the lower excitation temperature used by Zurlo et al. (2017). We suggest that, contrary to the case of Reipurth 50, the extended emission, which is likely resolved out by their interferometric observations, contains more of the mass of the outflow of V2775 Ori than the narrow cavity walls.

**V899 Mon**—This source with a flat or Class II spectral energy distribution (SED) was originally reported as a FUor by Wils et al. (2009). Follow-up observations indicated that the source was dimming, which was interpreted as a decrease in accretion rate by Ninan et al. (2015). The authors also recovered a PCygni profile for several forbidden lines, indicating the presence of outflows. Our observations are the first to recover an indication of a bipolar molecular outflow in CO, which follows a northeast–southwest direction. However, the outflow position angle disagrees with the position angle of the disk (Kóspál et al. 2021) and of the jets detected at optical wavelengths (Park et al. 2021); thus, the analysis of outflows with higher angular resolution is needed to resolve this discrepancy. Based on its channel maps (Appendix A), there is also significant extended emission at low velocities due to the envelope.

**V960 Mon**—Based on its pre-outburst SED, this is a Class II object (Kóspál et al. 2015), and our observations are the first to study the gas surrounding the system. The ¹²CO line profiles show high-velocity wings (Figure 3), which we interpreted as an indication of a bipolar molecular outflow. The integrated emission maps (Figures 4 and 5) and channel maps (Appendix A) show the two outflow lobes overlapping, an indication of the outflow having a direction along the line of sight. High angular resolution ALMA observations barely resolve the FUor disk and indicate a disk inclination between 16° and 60°, depending on the method used (Kóspál et al. 2021). Therefore, for the rest of the analysis, we assumed the lower inclination angle for the outflow. Kóspál et al. (2015) detected two sources close to this FUor (one to the north and one to the southeast), and Kóspál et al. (2021) found a third one to the east. As these sources are located within our beams, our observations contain emission from these neighboring sources, and it is possible that a source other than the FUor drives the outflow. Therefore, our results for the outflow around this FUor must be taken with caution, as an analysis of higher angular resolution observations is needed.

**ZCMa**—This source is a binary composed of the FUor and a Herbig Be star with a separation of 0.1". Levreault (1988) did

not detect an outflow in the *J* = 1–0 transition of ¹³CO and the *J* = 2–1 and *J* = 1–0 transitions of ¹²CO. Evans et al. (1994) and Liljeström & Olofsson (1997) detected the bipolar outflow emanating from this FUor in the *J* = 3–2 and *J* = 1–0 transitions of CO, respectively. Our observations recovered emission from the outflow with a northwest–southeast orientation (similar to that found in previous works), and we find that the outflow is compact and has low velocities. Our estimations of the outflow properties are in general agreement with those of Evans et al. (1994). It is unknown which of the two binary components drives the outflow, however; since both sources drive jets (Whelan et al. 2010), it is possible that both sources drive outflows.

**GM Cha**—The outflows around this Class I/II object had been previously reported using a single-dish antenna (e.g., Mottram et al. 2017) and ALMA (Hales et al. 2020). We recovered the east–west outflow orientation found by these authors. Comparing our results with those of Mottram et al. (2017), we find that the redshifted lobe is more extended than the blueshifted side. Our estimation of the mass for the blueshifted lobe is higher than their estimations, which is explained by us integrating lower-velocity fluxes compared to them. The redshifted mass and other outflow properties are comparable to those by Mottram et al. (2017).

**V346 Nor**—Kóspál et al. (2017a) presented single-dish observations of the *J* = 3–2 and *J* = 4–3 transitions, and Kóspál et al. (2017b) presented ALMA Cycle 2 observations of the *J* = 2–1 transition. Our analysis uses the same observations as Kóspál et al. (2017a), and the properties of the outflow have the same values within 10%. The small differences are due to slight differences in the methodology, such as different apertures and systemic velocities. The orientation of the outflow is the same as that obtained at high angular resolution (Kóspál et al. 2017b). Based on the ¹²CO/¹³CO ratio used in the optical depth correction, it appears that even the rarer isotopologue is optically thick at velocities close to the systemic.

### 6.3. FUors with Tentative Detections

Below, we present the two FUors for which we can only make a tentative detection of their outflows; thus, we consider that these two sources require follow-up observations.

**V900 Mon**—One of the most recently discovered FUors, it is a Class I source bordering on Class II (Reipurth et al. 2012). Kóspál et al. (2017b) used the same data and carried out a similar analysis as us, and they did not find outflow emission. However, Takami et al. (2019) presented high angular resolution ALMA observations of the *J* = 2–1 transition of ¹²CO, ¹³CO, and C¹⁸O, where they identified a bipolar outflow where the redshifted and blueshifted lobes are in the east and west directions, respectively. Following their results, we searched for the velocity ranges that could be integrated in the *J* = 3–2 transition for which we could find emission that follows that detected in the *J* = 2–1 observations. We found bipolar emission only in the *J* = 3–2 transition that follows a similar east–west alignment (see Figure 4) using the velocity range indicated in Table 3; thus, we considered this source to drive an outflow and estimated its properties. The *J* = 4–3 transition does not show significant emission (see Figure 5), which prompted us to consider this as only a tentative detection.

**iPTF 15afq**—This Class I object is one of the latest discovered FUors. It showed an ~2.5 mag brightening in
<table><caption>**Table 6**
Outflow Properties from ¹²CO (3–2) Observations after Optical Depth Correction</caption><thead><tr><th>Target</th><th>Side</th><th><i>M</i><sub>of</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>P</i><sub>of</sub><br>(<i>M</i><sub>☉</sub> km s<sup>-1</sup>)</th><th><i>E</i><sub>of</sub><br>(erg)</th><th><i>F</i><sub>of</sub><br>(<i>M</i><sub>☉</sub> yr<sup>-1</sup> km s<sup>-1</sup>)</th><th><i>L</i><sub>of</sub><br>(<i>L</i><sub>☉</sub>)</th></tr></thead><tbody><tr><td rowspan="2">L1551 IRS 5</td><td>Blue</td><td>2.2 × 10<sup>-2</sup></td><td>5.2 × 10<sup>-2</sup></td><td>1.4 × 10<sup>42</sup></td><td>6.4 × 10<sup>-6</sup></td><td>1.5 × 10<sup>-3</sup></td></tr><tr><td>Red</td><td>3.5 × 10<sup>-2</sup></td><td>6.6 × 10<sup>-2</sup></td><td>1.4 × 10<sup>42</sup></td><td>7.6 × 10<sup>-6</sup></td><td>1.4 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">Haro 5a IRS</td><td>Blue</td><td>9.5 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-1</sup></td><td>4.2 × 10<sup>42</sup></td><td>8.2 × 10<sup>-6</sup></td><td>1.5 × 10<sup>-3</sup></td></tr><tr><td>Red</td><td>2.8 × 10<sup>-2</sup></td><td>8.1 × 10<sup>-2</sup></td><td>2.4 × 10<sup>42</sup></td><td>3.7 × 10<sup>-6</sup></td><td>9.1 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">Reipurth 50</td><td>Blue</td><td>2.9 × 10<sup>-2</sup></td><td>9.2 × 10<sup>-2</sup></td><td>2.9 × 10<sup>42</sup></td><td>4.9 × 10<sup>-6</sup></td><td>1.3 × 10<sup>-3</sup></td></tr><tr><td>Red</td><td>3.4 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-1</sup></td><td>7.2 × 10<sup>42</sup></td><td>5.4 × 10<sup>-6</sup></td><td>2.0 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">V2775 Ori</td><td>Blue</td><td>9.1 × 10<sup>-2</sup></td><td>2.7 × 10<sup>-1</sup></td><td>9.6 × 10<sup>42</sup></td><td>7.8 × 10<sup>-5</sup></td><td>2.3 × 10<sup>-2</sup></td></tr><tr><td>Red</td><td>9.3 × 10<sup>-2</sup></td><td>2.1 × 10<sup>-1</sup></td><td>5.3 × 10<sup>42</sup></td><td>4.8 × 10<sup>-5</sup></td><td>9.9 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">V899 Mon</td><td>Blue</td><td>6.8 × 10<sup>-2</sup></td><td>8.9 × 10<sup>-2</sup></td><td>1.3 × 10<sup>42</sup></td><td>1.1 × 10<sup>-6</sup></td><td>1.3 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>8.5 × 10<sup>-2</sup></td><td>1.4 × 10<sup>-1</sup></td><td>2.2 × 10<sup>42</sup></td><td>3.6 × 10<sup>-6</sup></td><td>4.9 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">V900 Mon<sup>†</sup></td><td>Blue</td><td>2.5 × 10<sup>-2</sup></td><td>5.1 × 10<sup>-2</sup></td><td>1.0 × 10<sup>42</sup></td><td>1.3 × 10<sup>-6</sup></td><td>2.1 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>2.4 × 10<sup>-1</sup></td><td>1.8 × 10<sup>-1</sup></td><td>1.4 × 10<sup>42</sup></td><td>1.4 × 10<sup>-6</sup></td><td>9.2 × 10<sup>-5</sup></td></tr><tr><td rowspan="2">V960 Mon</td><td>Blue</td><td>1.9 × 10<sup>-1</sup></td><td>6.5 × 10<sup>-1</sup></td><td>2.4 × 10<sup>43</sup></td><td>6.2 × 10<sup>-5</sup></td><td>1.9 × 10<sup>-2</sup></td></tr><tr><td>Red</td><td>8.2 × 10<sup>-2</sup></td><td>3.1 × 10<sup>-1</sup></td><td>1.2 × 10<sup>43</sup></td><td>2.3 × 10<sup>-5</sup></td><td>7.1 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">Z CMa</td><td>Blue</td><td>8.1 × 10<sup>-2</sup></td><td>1.7 × 10<sup>-1</sup></td><td>3.5 × 10<sup>42</sup></td><td>3.8 × 10<sup>-6</sup></td><td>6.6 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>1.2 × 10<sup>-1</sup></td><td>3.2 × 10<sup>-1</sup></td><td>9.2 × 10<sup>42</sup></td><td>7.8 × 10<sup>-6</sup></td><td>1.9 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">iPTF 15afq<sup>†</sup></td><td>Blue</td><td>2.2 × 10<sup>-2</sup></td><td>7.1 × 10<sup>-2</sup></td><td>2.3 × 10<sup>42</sup></td><td>2.2 × 10<sup>-6</sup></td><td>6.0 × 10<sup>-4</sup></td></tr><tr><td>Red</td><td>3.8 × 10<sup>-2</sup></td><td>1.8 × 10<sup>-1</sup></td><td>8.9 × 10<sup>42</sup></td><td>4.2 × 10<sup>-6</sup></td><td>1.7 × 10<sup>-3</sup></td></tr><tr><td rowspan="2">GM Cha</td><td>Blue</td><td>1.3 × 10<sup>-3</sup></td><td>2.0 × 10<sup>-3</sup></td><td>3.1 × 10<sup>40</sup></td><td>3.1 × 10<sup>-7</sup></td><td>3.9 × 10<sup>-5</sup></td></tr><tr><td>Red</td><td>1.2 × 10<sup>-3</sup></td><td>3.4 × 10<sup>-3</sup></td><td>1.0 × 10<sup>41</sup></td><td>5.4 × 10<sup>-7</sup></td><td>1.3 × 10<sup>-4</sup></td></tr><tr><td rowspan="2">V346 Nor</td><td>Blue</td><td>1.9 × 10<sup>-2</sup></td><td>6.9 × 10<sup>-2</sup></td><td>2.8 × 10<sup>42</sup></td><td>5.2 × 10<sup>-6</sup></td><td>1.8 × 10<sup>-3</sup></td></tr><tr><td>Red</td><td>4.6 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-1</sup></td><td>9.9 × 10<sup>42</sup></td><td>9.8 × 10<sup>-6</sup></td><td>3.9 × 10<sup>-3</sup></td></tr></tbody></table>

Note. The † labels the two FUors with tentative outflow detections.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison of outflow masses determined from the *J* = 3–2 and *J* = 4–3 transitions. The dashed line indicates a ratio of 1.</figcaption></figure>

2015, which lasted a few months (Miller et al. 2015), and follow-up brightenings in 2018 and 2019 (Hillenbrand 2019). The 2019 outburst lasted until early 2021 and was followed by another outburst that was ongoing as of this writing.<fnref n="12" /> Hillenbrand (2019) presented high-resolution (*R* = 37,000) spectra taken during outburst and found that Hα and the Ca II triplet showed a P Cygni profile, an indicator of high-velocity winds. Our observations are the first submillimeter-wavelength

observations of this object. Its CO line profiles (Figure 3) show high-velocity line wings, particularly on the redshifted side. Based on its *J* = 3–2 integrated emission maps (Figure 4) and channel maps (Appendix A), there appears to be an outflow whose blueshifted and redshifted lobes are in the southeast and northwest directions, respectively. The blueshifted component is broader and with lower velocities than its redshifted counterpart. However, we consider this FUor as only a tentative detection because the emission is heavily dominated by the envelope; thus, it is hard to confirm the morphology seen in the *J* = 3–2 transition as an outflow.

## 6.4. FUors without Outflow Detections

ALMA observations of the *J* = 2–1 transition showed outflow emission for two FUors: V883 Ori (Ruíz-Rodríguez et al. 2017a) and V1647 Ori (Principe et al. 2018). There could be multiple causes behind our lack of detection: the combination of the higher sensitivity and angular resolution in the ALMA observations, the possible low temperatures in the system, and the low velocities of the ALMA outflows. Indeed, in the case of V883 Ori, White et al. (2019) found that the emission of ¹³CO *J* = 3–2 was a combination of the outflow at low velocities and a spherical-like envelope. In addition, when considering interferometric observations, it is possible that they have resolved out the contribution from the envelope, which our single-dish observations did not; therefore, the envelope emission dominates in the low-velocity channels of our observations. A similar case was found for FU Ori, for which previous observations reported that it did not drive an outflow (Levreault 1988), but ALMA observations of *J* = 2–1 hint toward an outflow, thus making its detection uncertain (Pérez et al. 2020). We detect emission in the northeast–southwest direction, which is perpendicular to the position angles of the resolved disks (Pérez et al. 2020); however, the angular

<footnote n="12">¹² http://gsaweb.ast.cam.ac.uk/alerts/alert/Gaia19fct/</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the outflows from quiescent sources in the literature with the FUor outflows. Shown are the outflow forces ($J = 3–2$) plotted against envelope masses (panel (a)), outflow masses ($J = 3–2$; panel (b)), ratio between outflow mass and envelope mass (panel (c)), and bolometric luminosities (panel (d)). In case of the FUors, all bolometric luminosities are during outburst. The two FUors with tentative outflow detections are marked with open diamonds. The $L_{bol}$ of iPTF 15afq is unknown and thus not shown in panel (d).</figcaption></figure>

resolution of our observations prevents us from determining if it is a bipolar outflow, so we do not consider this as a detection.

For the remaining FUors (AR 6A, BBW 76, V582 Aur, V723 Car, and OO Ser), we did not detect outflows, and we did not find previous publications that reported outflows.

### 6.5. Outflow Parameters

We detect clear outflow emission in ~55% of the FUors in our sample (10 out of the 18), which is lower than the 92% found in Class 0 and I objects (e.g., Mottram et al. 2017). Even including the two cases where ALMA detected outflows when we did not (V883 Ori and V1647 Ori) and the possible outflow in FU Ori, we would only find outflows in ~73% of the FUors of our sample. However, this is not surprising, considering that some of the FUors in our sample are classified as flat-spectrum or Class II objects, and the outflows in these evolved stages might be harder to detect due to the lower densities of the enveloping material (Arce & Sargent 2006). Indeed, only two FUors of Class II had evidence of outflows: V960 Mon and GM Cha.

After the optical depth correction, the outflow masses increased by a median factor of 3, with values ranging between 1 (V346 Nor) and 14 (V900 Mon). The values of the kinematic properties (e.g., momentum, energy) after the optical depth correction are also a factor of a few higher than without the correction. Even after applying this correction, we can still expect the outflow properties presented in Table 6 to be underestimated, as explained in Section 5.4.3.

In Figure 7, we show a comparison of the outflow masses determined for the FUors that had outflow detection in both transitions of $^{12}$CO. For this comparison, we estimated the outflow mass uncertainties by multiplying the number of pixels used when calculating the mass by the rms of each data cube and then converted these fluxes to masses using the same assumptions as the outflows. We find that within these uncertainties, most outflows have comparable masses in both transitions. Reipurth 50 is the only FUor in which the mass estimate is higher in the $J = 4–3$ transition than in $J = 3–2$ by a factor of ~2. The outflows of V2775 Ori and V899 Mon are more massive in the lower transition by factors of 4 and 5, respectively; thus, this could be an indication of different excitation properties causing the lower transition to be stronger. However, these comparisons are limited by differences in the observations (i.e., angular resolution and sensitivities) and

<table><caption>**Table 7**
Ratio between Outflow Masses and Envelope Masses using the $J = 3–2$ Transition of the Two Observed CO Isotopologues and the Core-to-star Formation Efficiency, $\epsilon$</caption><thead><tr><td>Target</td><td><i>M</i><sub>outflow</sub>/<i>M</i><sub>envelope</sub></td><td>ε</td></tr></thead><tbody><tr><td>L1551 IRS 5</td><td>0.068</td><td>-2.91</td></tr><tr><td>Haro 5a IRS</td><td>0.024</td><td>0.50</td></tr><tr><td>Reipurth 50</td><td>0.018</td><td>0.69</td></tr><tr><td>V2775 Ori</td><td>0.178</td><td>-19.23</td></tr><tr><td>V899 Mon</td><td>0.159</td><td>0.05</td></tr><tr><td>V900 Mon<sup>†</sup></td><td>0.140</td><td>0.44</td></tr><tr><td>V960 Mon</td><td>0.271</td><td>-8.89</td></tr><tr><td>Z CMa</td><td>0.055</td><td>0.39</td></tr><tr><td>iPTF 15afq<sup>†</sup></td><td>0.496</td><td>-4.73</td></tr><tr><td>GM Cha</td><td>0.010</td><td>0.23</td></tr><tr><td>V346 Nor</td><td>0.159</td><td>-2.77</td></tr></tbody></table>

Note. The † labels the two FUors with tentative outflow detections.

images (i.e., pixel size and field of view) and by using the same excitation for the two transitions in all of the FUors. A large-scale program to target multiple CO transitions under comparable conditions would alleviate these limitations and provide more insight on the masses of the outflows.

### 6.6. Comparison with Quiescent YSOs

We put our outflow properties into context by comparing them with the values of similar studies based on quiescent sources. This comparison is not straightforward due to the differences in the observational properties (i.e., angular resolution and sensitivity) and methodology (i.e., choosing velocities for integration, optical depth correction, and inclination correction), which have significant effects on the resulting values of the outflow properties. It is expected that FUor outbursts last for up to 100 yr, and the dynamical ages of the outflows are on the order of thousands of years (see Table 3); thus, the outflows we have detected around FUors are not related to the current outbursts. This means that we are comparing the histories of the two samples, which could provide hints toward the nature behind the outbursts.

We compared our sample with the values of the outflow properties published in the following studies: Dunham et al. (2014), Yıldız et al. (2015), and Mottram et al. (2017). The
<table><caption>**Table 8**

p-values of the Three Statistical Tests Done for the Envelope Masses, Outflow Masses and Forces, and the Ratio between Outflow and Envelope Masses</caption><thead><tr><th rowspan="2">Test</th><th colspan="2">Envelope Masses</th><th colspan="2">Outflow Masses</th><th colspan="2">O/E Mass Ratio</th><th colspan="2">Outflow Forces</th></tr><tr><th>Class 0</th><th>Class I</th><th>Class 0</th><th>Class I</th><th>Class 0</th><th>Class I</th><th>Class 0</th><th>Class I</th></tr></thead><tbody><tr><td>K-S</td><td>0.011</td><td>0.605</td><td>0.220</td><td>1 × 10<sup>-3</sup></td><td>2 × 10<sup>-3</sup></td><td>2 × 10<sup>-4</sup></td><td>0.070</td><td>0.247</td></tr><tr><td>MWU</td><td>2 × 10<sup>-3</sup></td><td>0.667</td><td>0.129</td><td>0.002</td><td>5 × 10<sup>-4</sup></td><td>2 × 10<sup>-4</sup></td><td>0.272</td><td>0.726</td></tr><tr><td>kAD</td><td>&lt;0.001</td><td>&gt;0.250</td><td>0.200</td><td>0.002</td><td>&lt;0.001</td><td>&lt;0.001</td><td>0.097</td><td>0.208</td></tr></tbody></table>

Note. The p-values from anderson_ksamp are capped between 0.001 and 0.250; thus, these values are upper and lower limits, respectively.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Cumulative histograms of the outflow properties and outflow-to-envelope mass ratio for the Class 0 and I objects from the literature (Dunham et al. 2014; Yıldız et al. 2015; Mottram et al. 2017) and the FUors from this work. The bin widths for each histogram were selected using the Freedman–Diaconis rule.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Ratio between outflow mass and envelope mass plotted against the distance to each target. The open diamonds are the two FUors with tentative outflow detections.</figcaption></figure>

three studies cover a combination of Class 0 and Class I objects, and our calculations followed similar methods to theirs. We included quiescent Class 0 objects even when the vast majority of FUors were Class I objects (see Table 1) because we wanted to compare how the FUor outbursts compare to the different stages of the star formation process. We compared the outflow masses and forces, as those were the

only two properties presented by all three studies from the literature.

Dunham et al. (2014) presented outflow properties with and without optical depth correction, while Yıldız et al. (2015) and Mottram et al. (2017) did not calculate this correction; thus, we used the optically thin values for this comparison.

Dunham et al. (2014) assumed $T_{ex} = 50$ K for their calculations, while Yıldız et al. (2015) and Mottram et al. (2017) used the same temperature as we did in our analysis, $T_{ex} = 75$ K. To test the effect of using the lower temperature, we calculated the outflow properties of the FUors using $T_{ex} = 50$ K and found that the values were a factor of ~1.19 higher when using the higher temperature. Thus, we multiplied the outflow properties of Dunham et al. (2014) by this factor to minimize the differences in methodology. The outflow forces estimated by Yıldız et al. (2015) and Mottram et al. (2017) were corrected due to the inclination of the systems based on Downes & Cabrit (2007). These correction values were estimated for Class 0 objects and are not recommended to correct Class I objects; however, for the sake of a comparison between our sample and the ones from the literature, we applied this correction factor to the FUor outflows even if they were at later stages (Class I or II). We estimated the inclination of the outflows by assuming that the inclination of the outflow was perpendicular to the inclination of the disk. For most sources, we used the inclination of the FUor disks obtained from high angular resolution observations with ALMA (Cieza et al. 2018; Hales et al. 2020; Kóspál et al. 2021), while for ZCMa, we used the estimate by Antoniucci et al. (2016) from an analysis with data from optical interferometry, and in the case of iPTF 15afq, we assumed an inclination of 45°. However, the works studying the quiescent sample used a
coarse correction table, e.g., Table A.5 in Mottram et al. (2017). Thus, we rounded our inclinations to the closest values in the inclination correction table, and the angles assumed are listed in Table 3. We combined the quiescent samples from the literature into one, divided it by Class, and compared the two subsamples with the FUors.

In Figure 8, we plot the outflow forces calculated from the $J=3-2$ transition against the envelope masses calculated from the $^{13}$CO emission (panel (a)), the outflow masses also calculated from the $^{12}$CO $J=3-2$ transition (panel (b)), the ratio between the outflow mass and envelope mass (panel (c)), and the bolometric luminosities obtained from the literature (panel (d); Table 1).

Panels (a) and (b) of Figure 8 show that FUor outflows follow the same trends as the quiescent sources from the literature; i.e., higher envelope and outflow masses indicate higher outflow forces. In the outflow-to-envelope mass ratio subplot, panel (c) of Figure 8, the FUor sample is offset from the quiescent samples. This ratio has been used to discuss the core-to-star formation efficiency in the quiescent sample (Mottram et al. 2017); thus, it hints that FUors are less efficient at driving mass from the envelope onto the star, and this relationship will be discussed further below. The values for this ratio are presented in Table 7.

The correlation between outflow force and bolometric luminosity, panel (d) of Figure 8, has been well studied for quiescent sources (Cabrit & Bertout 1992; Bontemps et al. 1996; Yıldız et al. 2015; Mottram et al. 2017), and it would appear that FUors do not follow this correlation. However, the FUor bolometric luminosities were estimated from photometry taken while in outburst, and none have sufficient pre-outburst photometric data to estimate their pre-outburst luminosities. Even if we do not have sufficient information about the individual FUors, when in quiescence, the protostars are expected to be low-mass and low-luminosity objects (except for V723 Car), and our measured outflow parameters are consistent with this. Thus, our results suggest that FUors, when in quiescence, produce molecular outflows with forces comparable to those from outflows in quiescent stars.

In order to get a better estimate of how similar FUor outflows are to their quiescent counterparts, we present cumulative histograms comparing different properties (Figure 9), and we carried out three complementary statistical tests to examine whether the samples of the quiescent young stars were drawn from the same sample as the FUors. The first was a two-sided Kolmogorov–Smirnov (K-S) test that compares the shapes of the distributions; the second was a Mann–Whitney U-test (MWU), which is more sensitive to the mean of the two samples rather than the shape of both distributions; and the third was a k-sample Anderson–Darling test (kAD), which is more sensitive to the tails of the distributions. The three tests were done using the SciPy functions `kstest`, `mannwhitneyu`, and `anderson_k samp`, respectively. The results of the statistical tests are presented in Table 8.

With a significance level of 5%, we found that the distribution of FUor envelope masses is similar to that of Class I and different from Class 0s (Figure 9(a)), the outflow masses of FUors are different from those of Class I objects and comparable with those of Class 0s (Figure 9(b)), the outflow-to-envelope mass ratio is different in the FUors when compared to either sample (Figure 9(c)), and the outflow forces of FUors are

comparable with the two quiescent samples (Figure 9(d)). Our tests did not include the two tentative detections. We ran the tests including these two FUors and found that our statistical results would not change.

The envelope mass result is not surprising because, based on their SEDs, the FUors are also Class I. Following the same logic, the result for the outflow masses is surprising, as the outflow masses are close to those of the Class 0s. This could be interpreted as an indication that FUors are in the very early part of their Class I stage, and their outflows have not had sufficient time to dissipate. However, the masses of the outflows are a combination of the material that has passed through the accretion disk and is now being driven away from the star by the outflow and the material in the envelope that has been entrained by the outflow. The FUor outflows have higher outflow masses but similar envelope masses compared to the Class I, pointing toward FUor outflows having a higher percentage of material that was ejected from the accretion disk, i.e., material that was not accreted onto the star. This can be seen in the ratio between the outflow and envelope mass in Figures 8 and 9.

The separation between the quiescent sample and the FUors might be biased due to the distance of the targets. The YSOs of the quiescent samples are all within 500 pc from the Sun, while half of the FUor sample is beyond this distance. Therefore, if we consider that in most cases, the extension of the outflows is larger than that of the envelopes, our analysis might be biased toward the FUors, because it is likely that we are measuring the full extension of their outflows in comparison to the quiescent sample. In Figure 10, we plot the outflow envelope mass ratio versus the distance for each target. The distribution of points indicates that indeed there might be a positive correlation between the mass ratio and the distance. However, the maps in Yıldız et al. (2015) and Mottram et al. (2017) indicate that 40% of the outflows extend beyond the areas of the sky covered by their respective observations. Therefore, the mass ratios for those sources are lower limits, thus raising the question of whether this relationship is real or not. An in-depth observational program covering outflows at a wide range of distances should shed some light on this matter.

Mottram et al. (2017) used the outflow envelope mass ratio to examine the core-to-star efficiency in a group of quiescent young stars. They assumed that during the whole duration of the Class 0 and I phases, a star has an outflow rate with small enough variations that it can be approximated by a constant value, and calculated the core-to-star formation efficiency, $\epsilon$, as follows:

$$ \epsilon = 1 - \frac{M_{\text{of}}}{M_{\text{env}}} \frac{\tau_{0+1}}{\tau_{\text{d}}}, \quad (4) $$

where $\tau_{0+1}$ is the total duration of the Class 0 and I phases, i.e., 0.5 Myr. The authors used their typical values of $M_{\text{of}}/M_{\text{env}} = 10^{-2}$ and $\tau_{\text{d}} = 10^4$ yr and found a core-to-star formation efficiency of 0.5, which is in agreement with the literature. We used our derived masses and dynamical times to calculate $\epsilon$ for the FUors. The values can be found in Table 7.

As can be seen from our results, five of the FUors have negative $\epsilon$ values. These can be explained by the underestimation of the dynamical age because the outflow either appears to be face-on (V2775 Ori and V960 Mon) or extends beyond the field of view of our observations (L1551 IRS 5, iPTF 15afq, and V346 Nor). Three of the six FUors with
positive values extend beyond our field of view (Haro 5a IRS, Reipurth 50, and V899 Mon); thus, their $\epsilon$ values are uncertain because it is unknown how much of the mass and the extension of the outflow is beyond our field of view. The three remaining FUors (V900 Mon, Z CMa, and GM Cha) with positive $\epsilon$ and the full outflow inside the field of view have lower efficiencies than the quiescent sample. These values would suggest that a significant amount of material that was fed from the envelope onto the disk was not accreted onto the star but instead was driven outward by the outflow. However, two of these have strong caveats. For V900 Mon, its outflow was only tentatively detected, and follow-up observations might reveal a different geometry than ours, which would lead to a different value of $\epsilon$. Object GM Cha is a Class II object, i.e., in a more evolved stage than most FUors, and the equation used to calculate the efficiencies was created for Class I objects whose accretion rate is orders of magnitude higher than that for Class IIs (Fiorellino et al. 2023) and have higher envelope masses, which means that our estimated value is not accurate.

As a final point, we address the similar distributions of the outflow forces between the two quiescent samples and the FUors. The outflow force is a property that is commonly associated with the accretion history of a YSO (Bontemps et al. 1996). This is because, as can be seen in panel (a) of Figure 8, there is a positive correlation between the outflow force and the envelope mass, and the more evolved young stars have lower envelope masses and mass accretion rates. Here we present some scenarios that can explain the lack of separation between the FUors and quiescent sample.

First, if we assume that the similar distributions are because the accretion histories of the two samples are the same, then this can be interpreted as either the “quiescent” sample had outbursts that were undetected, or the current outbursts in the FUors are the first ones in their accretion histories. If these are indeed the first outbursts in each of the FUors, then their effects would be undetected because the angular resolution of our observations is insufficient to resolve the inner parts of the outflows where the effects of a <100 yr old outburst could be detected. Few FUor outbursts (~20; Connelley & Reipurth 2018) have been detected, so the incidence rate of these events is unknown; as such, it is difficult to separate these two scenarios.

Second, if we assume that the accretion histories are different between quiescent and outbursting samples, then the likeness between distributions should be because of a physical property of the outflows. The outflow force is calculated using the outflow momentum and the dynamical age of the outflow, and both of these properties depend on the distribution of velocities in the outflow. If the FUor outflows have masses comparable to Class 0 outflows and similar velocity ranges as the outflows from the literature, we would expect the FUor outflow forces to be close to those of the Class 0 objects from the other samples. Therefore, the divergence between Class 0 and FUor outflow forces indicates that the latter have lower velocities. This had already been mentioned by Principe et al. (2018) when comparing the V1647 Ori outflow to other outflows observed with ALMA.

However, as mentioned earlier, the outflow forces presented here are highly uncertain because they are calculated as the ratio of two lower limits, the outflow momentum and the outflow dynamical age, and because of the understudied effects of the inclination of the outflow with respect to the plane of the

sky. As such, there is a need for a thorough program to study these molecular outflows.

## 7. Summary and Conclusions

We presented APEX observations of 20 FUors or FUor-like objects from which we estimated the envelope mass and searched for outflows. Using a combination of line profiles and inspection of channel maps, we detected outflows in 45% of our sample. These include the possible first detections of molecular outflows in V899 Mon and V960 Mon, although these should be observed with higher angular resolution to corroborate them. We also found two tentative detections in V900 Mon and iPTF 15afq that require follow-up observations to confirm. In the case of V883 Ori, V1647 Ori, and possibly FU Ori, we did not detect the outflows that have been observed by ALMA.

Based on our $^{13}$CO measurements, envelopes with masses higher/lower than 0.1–0.2 $M_{\odot}$ show the silicate feature at 10 $\mu$m in absorption/emission. If the envelope mass is close to this threshold level, the geometry of the system determines whether the spectral feature is in emission or absorption. The most significant outlier of this trend is V960 Mon, which shows the 10 $\mu$m feature in emission despite having an envelope of ~0.6 $M_{\odot}$.

The masses of outflows estimated from the $^{12}$CO 3–2 and 4–3 transitions are in agreement, except for two FUors: V900 Mon and V960 Mon. We suggest that these two sources could be colder than the rest of the sample; thus, the higher transition is dimmer.

As V960 Mon is an outlier in both trends, we proposed another possible explanation. This FUor has three companion YSOs in its proximity, with separations smaller than the sizes of our beams. Therefore, we suggest that these additional sources move this object away from the trend seen in the rest of the sample.

The kinematic outflow properties (momenta, energies, forces, and luminosities) are higher when estimated from the $J=3-2$ transition than from $J=4-3$. We attribute this to the higher sensitivity of the lower transition, which causes a difference in the range of velocities in which we detected outflow emission.

After applying an optical depth correction to the $J=3-2$ transition using the $^{13}$CO emission, we found that the mass of the outflows increased by a median factor of 3 and up to an order of magnitude. The minimum improvement, seen in a few cases, showed that the outflow mass increased by only a few percent.

We compared the outflows found in our FUor sample with three works from the literature and found that outflows emanating from FUors are more massive than those from quiescent Class I sources but with masses comparable to outflows in Class 0 sources. We found that FUors have a higher outflow envelope mass ratio than the quiescent sample, although this result could be biased by the distance. We calculated the core-to-star efficiencies of the FUors, and although our results are severely constrained by the geometry of the outflows, it could indicate that a significant portion of the material that was deposited into the accretion disk from the envelope is not accreted onto the star but is instead driven back to the envelope by the outflow. Finally, we found that outflow forces from the FUor sample are comparable to the two
quiescent sources, which can be interpreted as similar accretion histories or very low velocities in the FUor outflows.

This study focuses on the outflow histories of the FUors observable from the APEX site. The dynamical ages of the detected outflows indicate that they are much older than any of the ongoing outbursts, which are less than 100 yr old. Indeed, any outflow emission directly related to the current outburst would be detected at high velocities close to the protostar and would have small spatial scales that would be diluted by the beam of our single-dish observations. Our comparison between the outflow properties of FUors and of other quiescent objects should be taken with caution due to the varying quality of the individual observations and the methodology used by each research group. A complete survey of all known FUors in both hemispheres with similar observational setups and sensitivities and a control sample of multiple quiescent YSOs at different evolutionary stages would greatly improve our analysis.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 716155 (SACCRED). T.Cs. has received financial support from the French State in the framework of the IdEx Université de Bordeaux Investments for the Future Program. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX) under program IDs 094.F-9508 and 098.F-9505. APEX is a collaboration between the Max-

Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory.

Facilities: APEX, Herschel, JCMT.

Software: NumPy (Harris et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020).

## Appendix A: Channel Maps

Here we present the channel maps for the three observed transitions of L1551 IRS 5 (Figure A1), and the complete figure set with the rest of the targets in the FUor sample is available in the online journal. The minimum and maximum velocities in the channel maps are those when the gas emission starts or finishes being significant. The purple contours are used for all of the ¹³CO channel maps and for the ¹²CO maps when outflows were not detected. The blue and red contours show the blueshifted and redshifted emission of outflows, and the green contours show the envelope emission. All channel maps show the aperture used for the calculation of the envelope mass in the case of ¹³CO and the outflow properties in the case of both ¹²CO transitions. The velocities shown in the plots were chosen so that the maximum and minimum velocities are shown within 27 frames, which can cause some irregular velocity steps in the plot. However, these differences are of one channel, and we do not expect to see significant changes in the distribution of CO between two continuous channels.
<figure><img src="image_14.png" /><figcaption>**Figure A1.** The L1551 IRS 5 channel maps: ¹³CO (3–2) with contours at 3σ, 5σ, 7σ, 9σ, 11σ, 13σ, 15σ, 17σ, 19σ, and 21σ with σ = 0.32 K; CO (3–2) with contours at 3σ, 6σ, 9σ, 12σ, 15σ, 18σ, 21σ, 24σ, 27σ, 30σ, 33σ, 36σ, and 39σ with σ = 0.37 K; and CO (4–3) with contours at 3σ, 5σ, 7σ, 9σ, 11σ, 13σ, 15σ, 17σ, 19σ, and 21σ with σ = 0.65 K. The black circles indicate the 8000 au aperture used to extract the line profiles used to analyze the outflow.</figcaption></figure>

(The complete figure set (19 images) is available.)
## Appendix B
Optical Depth Correction

In Table B1, we present the parameters of the best-fitted parabola used to determine the optical depth correction for the

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

**Figure B1.** Optical depth correction for L1551 IRS 5. In the left panel, we show the line profiles, where green, purple, and black represent the ¹³CO, observed ¹²CO, and corrected ¹²CO, respectively, and the vertical dashed lines indicate the range of velocities used in the parabola fit. The right panel shows the ratio of main beam temperatures ($T_{MB}$), where the light blue crosses are all of the values of the ratio for each velocity channel, the pink dots with error bars are the points used in the parabolic fitting, and the black line is the resulting best-fit parabola.

(The complete figure set (10 images) is available.)

<table><caption>**Table B1**
Parameters of the Best-fit Parabolas Used for Optical Depth Correction</caption><thead><tr><th>Target</th><th>A</th><th>C</th></tr></thead><tbody><tr><td>L1551 IRS 5</td><td>2.388</td><td>2.683</td></tr><tr><td>Haro 5a IRS</td><td>1.327</td><td>2.194</td></tr><tr><td>Reipurth 50</td><td>2.548</td><td>0.643</td></tr><tr><td>V2775 Ori</td><td>1.825</td><td>4.320</td></tr><tr><td>V899 Mon</td><td>1.138</td><td>6.189</td></tr><tr><td>V900 Mon</td><td>1.531</td><td>3.946</td></tr><tr><td>V960 Mon</td><td>1.409</td><td>3.357</td></tr><tr><td>Z CMa</td><td>0.924</td><td>3.350</td></tr><tr><td>iPTF 15afq</td><td>-6.751</td><td>6.224</td></tr><tr><td>GM Cha</td><td>-4.909</td><td>19.505</td></tr><tr><td>V346 Nor</td><td>-2.419</td><td>8.246</td></tr></tbody></table>

sources with outflows, and in Figure B1, we present the parabolic fit and line profiles used in the fitting.

## ORCID iDs

Fernando Cruz-Sáenz de Miera https://orcid.org/0000-0002-4283-2185

Ágnes Kóspál https://orcid.org/0000-0001-7157-6275

Péter Ábrahám https://orcid.org/0000-0001-6015-646X

Timea Csengeri https://orcid.org/0000-0002-6018-1371

Orsolya Fehér https://orcid.org/0000-0003-3453-4775

Thomas Henning https://orcid.org/0000-0002-1493-300X

## References

* Ábrahám, P., Kóspál, Á., Kun, M., et al. 2018, ApJ, 853, 28

* Andre, P., & Montmerle, T. 1994, ApJ, 420, 837

* Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122

* Antoniucci, S., Podio, L., Nisini, B., et al. 2016, A&A, 593, L13

* Arce, H. G., & Sargent, A. I. 2006, ApJ, 646, 1070

* Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

* Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

* Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 387

* Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

* Bally, J. 2016, ARA&A, 54, 491

* Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987

* Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

* Bontemps, S., Andre, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858

* Cabrit, S., & Bertout, C. 1992, A&A, 261, 274

* Cieza, L. A., Ruíz-Rodríguez, D., Perez, S., et al. 2018, MNRAS, 474, 4347

* Connelley, M. S., & Reipurth, B. 2018, ApJ, 861, 145

* Cruz-Sáenz de Miera, F., Kóspál, Á., Ábrahám, P., et al. 2022, ApJ, 927, 125

* Cruz-Sáenz de Miera, F., Kóspál, Á., Ábrahám, P., Liu, H. B., & Takami, M. 2019, ApJL, 882, L4

* Downes, T. P., & Cabrit, S. 2007, A&A, 471, 873

* Dunham, M. M., Arce, H. G., Mardones, D., et al. 2014, ApJ, 783, 29

* Ellerbroek, L. E., Podio, L., Dougados, C., et al. 2014, A&A, 563, A87

* Evans, N. J., I., Balkum, S., Levreault, R. M., Hartmann, L., & Kenyon, S. 1994, ApJ, 424, 793

* Fiorellino, E., Tychoniec, L., Cruz-Sáenz de Miera, F., et al. 2023, ApJ, 944, 135

* Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2022, arXiv:2203.11257

* Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 451

* Garufi, A., Podio, L., Bacciotti, F., et al. 2019, A&A, 628, A68

* Green, J. D., Evans, N. J., I, Kóspál, Á., et al. 2013, ApJ, 772, 117

* Green, J. D., Hartmann, L., Calvet, N., et al. 2006, ApJ, 648, 1099

* Greene, T. P., Wilking, B. A., Andre, P., Young, E. T., & Lada, C. J. 1994, ApJ, 434, 614

* Güsten, R., Booth, R. S., Cesarsky, C., et al. 2006, Proc. SPIE, 6267, 626714

* Hales, A. S., Corder, S. A., Dent, W. R. D., et al. 2015, ApJ, 812, 134

* Hales, A. S., Pérez, S., Gonzalez-Ruilova, C., et al. 2020, ApJ, 900, 7

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

* Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207

* Hillenbrand, L. A. 2019, ATel, 13321, 1

* Hunter, J. D. 2007, CSE, 9, 90

* Kim, K. H., Watson, D. M., Manoj, P., et al. 2016, ApJS, 226, 8

* Klein, T., Ciechanowicz, M., Leinz, C., et al. 2014, ITTST, 4, 588
Koresko, C. D., Beckwith, S. V. W., Ghez, A. M., Matthews, K., & Neugebauer, G. 1991, AJ, 102, 2073

Kóspál, Á., Ábrahám, P., Carmona, A., et al. 2020a, ApJL, 895, L48

Kóspál, Á., Ábrahám, P., Csengeri, T., et al. 2017a, ApJ, 836, 226

Kóspál, Á., Ábrahám, P., Csengeri, T., et al. 2017b, ApJ, 843, 45

Kóspál, Á., Ábrahám, P., Moór, A., et al. 2015, ApJL, 801, L5

Kóspál, Á., Cruz-Sáenz de Miera, F., White, J. A., et al. 2021, ApJS, 256, 30

Kóspál, Á., Szabó, Z. M., Ábrahám, P., et al. 2020b, ApJ, 889, 148

Kun, M., Szegedi-Elek, E., & Reipurth, B. 2017, MNRAS, 468, 2325

Lada, C. J. 1987, in IAU Symp. 115, Star Forming Regions, ed. M. Peimbert & J. Jugaku (Dordrecht: Reidel), 1

Lee, C.-F., Ho, P. P., Li, Z.-Y., et al. 2017, NatAs, 1, 0152

Levreault, R. M. 1988, ApJ, 330, 897

Liljeström, T., & Olofsson, G. 1997, ApJ, 478, 381

Lim, J., Yeung, P. K. H., Hanawa, T., et al. 2016, ApJ, 826, 153

Manoj, P., Kim, K. H., Furlan, E., et al. 2011, ApJS, 193, 11

Miller, A. A., Hillenbrand, L. A., Bilgi, P., et al. 2015, ATel, 7428, 1

Moriarty-Schieven, G. H., Johnstone, D., Bally, J., & Jenness, T. 2006, ApJ, 645, 357

Mottram, J. C., van Dishoeck, E. F., Kristensen, L. E., et al. 2017, A&A, 600, A99

Ninan, J. P., Ojha, D. K., Baug, T., et al. 2015, ApJ, 815, 4

Nony, T., Motte, F., Louvet, F., et al. 2020, A&A, 636, A38

Park, S., Kóspál, Á., Ábrahám, P., et al. 2022, ApJ, 941, 165

Park, S., Kóspál, Á., Cruz-Sáenz de Miera, F., et al. 2021, ApJ, 923, 171

Pérez, S., Hales, A., Liu, H. B., et al. 2020, ApJ, 889, 59

Plunkett, A. L., Arce, H. G., Mardones, D., et al. 2015, Natur, 527, 70

Poetzel, R., Mundt, R., & Ray, T. P. 1989, A&A, 224, L13

Postel, A., Audard, M., Vorobyov, E., et al. 2019, A&A, 631, A30

Principe, D. A., Cieza, L., Hales, A., et al. 2018, MNRAS, 473, 879

Quanz, S. P., Henning, T., Bouwman, J., et al. 2007, ApJ, 668, 359

Reipurth, B., Aspin, C., & Herbig, G. H. 2012, ApJL, 748, L5

Rodriguez, L. F., Hartmann, L. W., & Chavira, E. 1990, PASP, 102, 1413

Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2017a, MNRAS, 468, 3266

Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2017b, MNRAS, 466, 3519

Safron, E. J., Fischer, W. J., Megeath, S. T., et al. 2015, ApJL, 800, L5

Snell, R. L., Loren, R. B., & Plambeck, R. L. 1980, ApJL, 239, L17

Stojimirović, I., Narayanan, G., Snell, R. L., & Bally, J. 2006, ApJ, 649, 280

Takahashi, S., Saito, M., Ohashi, N., et al. 2008, ApJ, 688, 344

Takahashi, S., Saito, M., Takakuwa, S., & Kawabe, R. 2006, ApJ, 651, 933

Takami, M., Chen, T.-S., Liu, H. B., et al. 2019, ApJ, 884, 146

Takami, M., Fu, G., Liu, H. B., et al. 2018, ApJ, 864, 20

Tapia, M., Roth, M., & Persi, P. 2015, MNRAS, 446, 4088

Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130

van Kempen, T. A., van Dishoeck, E. F., Hogerheijde, M. R., & Gusten, R. 2009, A&A, 508, 259

Vazzano, M. M., Fernández-López, M., Plunkett, A., et al. 2021, A&A, 648, A41

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956

Whelan, E. T., Dougados, C., Perrin, M. D., et al. 2010, ApJL, 720, L119

White, J. A., Kóspál, Á., Rab, C., et al. 2019, ApJ, 877, 21

Wils, P., Greaves, J., Catelan, M., et al. 2009, ATel, 2307, 1

Wilson, T. L. 1999, RPPh, 62, 143

Wu, P.-F., Takakuwa, S., & Lim, J. 2009, ApJ, 698, 184

Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2015, A&A, 576, A109

Zhang, Y., Arce, H. G., Mardones, D., et al. 2019, ApJ, 883, 1

Zurlo, A., Cieza, L. A., Williams, J. P., et al. 2017, MNRAS, 465, 834