An APEX Study of Molecular Outflows in FUor-type Stars

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 ⊙ 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 12CO (3–2), 13CO (3–2), and 12CO (4–3) survey of 20 FUor-type eruptive young stars observed with APEX. We use our 13CO (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 μ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.


INTRODUCTION
Jets and molecular outflows are a ubiquitous phenomenon 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 km s −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 Corresponding author: Fernando Cruz-Sáenz de Miera cruzsaenz.fernando-at-csfk.orgcontinuum, 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 Class II objects (Mottram et al. 2017).The mass accretion rates from protostellar disks to protostars are expected to undergo episodic variations.De-tailed 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) show how the study of outflows can shed light on the accretion history of the protostars that drive them.
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 M yr −1 up to ∼ 10 −4 M yr −1 .These events are typically detected as a 3 − 5 magnitude brightening at optical and nearinfrared wavelengths, and are expected to last up to a century, meaning that these events increase the final stellar mass by a significant amount.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 in 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 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 none, or a minimal amount, of the spectral signatures, then it is considered as Peculiar.
Outflows play an important role in the star formation process as they remove angular momentum from the accretion disk, 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).
Each of the aforementioned studies focused on a single FUor-type object or on 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 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 is the histories of the FUors with those of the quiescent sample.The triggering mechanisms 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 protostar with intrinsic differences that caused the outburst.Alternatively, it could show the samples are similar and, 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 non-outbursting sources and draw our conclusions.Finally, in Section 7 we summarize our work and present our main findings.

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 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 their 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 their sudden increases of their mass accretion rate, and 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 re-processing 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.

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 12 CO (3-2), 13 CO (3-2), and 12 CO (4-3) lines towards our targets.Program 094.F-9508 was observed between 2014 August 23-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 the 13 CO (3-2) at 330.588 GHz, and the 12 CO (3-2) at 345.796 GHz, respectively.The higher frequency channel was tuned to the 12 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 m s −1 , ∼32.9 m s −1 and ∼49.4 m s −1 for the 13 CO (3-2), CO (3-2) and CO (4-3) lines, respectively.. For each target, 120 ×120 on-the-fly (OTF) maps were obtained at 6 s −1 , using a relative reference off position 1000 away in right ascension.
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 GHz and 464 GHz, respectively, and the values were converted to Jansky using 41 Jy K −1 and 48 Jy K −1 at 352 GHz and 464 GHz, respectively.We calculated the noise levels of each CO line by first selecting the first and the last 100 channels of each cube (individually confirmed to be free of line emission), calculated the noise levels for each FUor using these channels, and then we calculated 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 13 CO (3-2), 3.7 Jy for 12 CO (3-2), and 9.4 Jy for 12 CO (4-3).The telescope's half-power beam-width is 19.2, and 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.

Distribution of gas
We constructed velocity integrated emission maps (Moment 0) 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 0 maps with the dust continuum emission.We searched for 250 µm continuum maps taken with Herschel /SPIRE maps in the Herschel Science Archive 1 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 Centre2 for archival 850 µm observations taken with the James Clerk Maxwell Telescope3 (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 emission and the 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, a Here we refer as Class to the classification based on the shape of its SED (Lada 1987;Andre et al. 1993;Greene et al. 1994).b See Section 4.2.
c See Section 5.1.
f Also known as V912 Mon.
g Also known as Gaia19fct.
h Also known as Parsamian 21.
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.

Systemic velocity
As it is 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 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 Moment 0 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 contours of the 850 µm continuum emission from the JCMT (see Section 4.1).The CO and the dust contours are plotted with levels at 0.3, 0.4, . . ., 0.9 of the peak intensity, and are meant to be representative.The star symbols indicate the nominal positions of the protostars.The bars at the bottom of each panel represent 10 000 au.
The line profiles and the 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 self-absorption of the envelope or due to the rotation of the gas.The former scenario is more likely based on the inspection of the channel maps of 13 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, out-flows, Keplerian disk, or unrelated gas in the same lineof-sight) showing emission at 13 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 −1 .For the symmetrical sources, the differences are less than 0.3 km s −1 .In addition, we examined the channel maps of each target to confirm our systemic velocity estimates.Figure 2. Line profiles of 13 CO (black lines) extracted using a circular aperture with 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 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 separate from the single Gaussian fit.See Section 4.2 for details.
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 −1 for L1551 IRS 5 (Wu et al. 2009), 0.02 km s −1 for V2775 Ori (Zurlo et al. 2017), 0.04 km s −1 for GM Cha (Mottram et al. 2017), 0.20 km s −1 for V883 Ori (Ruíz-Rodríguez et al. 2017a), 0.06 km s −1 for V1647 Ori (Principe et al. 2018), 0.35 km s −1 for V582 Aur ( Ábrahám et al. 2018), 0.27 km s −1 for V900 Mon (Takami et al. 2019), 0.56 km s −1 for FU Ori North (Pérez et al. 2020, who resolved the binary system with an angular resolution of 0.05" using ALMA).In the cases of Haro 5a IRS, AR 6A, BBW 76, OO Ser, and HBC 687 the esti-mates by (Kóspál et al. 2017a) are in agreement within 0.30 km s −1 .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 18 O: V346 Nor and Reipurth 50.In the case of the former, Kóspál et al. (2017b) found the line profile peaks at −3.55 km s −1 , indicating a difference of 0.47 km s −1 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 −1 , a value 0.77 km s −1 different from ours.It is likely that the larger 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.

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 which regions show emission above the 3σ contour level.The channels with wide extended emission that showed little variations from channel to channel were considered as envelope emission.Then we inspected the blue-and red-shifted channel maps for emission similar to what is found in outflows, i.e. emission whose red-shifted channels is in the opposite side from the blue-shifted channels with respect to the expected position of the star, and emission that is generally more extended in the channels with velocities closer to the systemic velocity and more compact towards 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 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 their 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.
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 obvious indications of wings caused by high velocity outflows.The remainder of the sources exhibit much wider profiles with 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 0 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 in 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).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.

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 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 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 from 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 difference between our and their estimates are −79 pc, 90 pc, 30 pc 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 masses of factors of 0.69, 1.24, 1.06 and 0.37, respectively.

Envelope masses
We used the 13 CO (3-2) emission to calculate the masses of the envelopes surrounding the FUors.To calculate the integrated fluxes, we used the line pro-files defined in Section 4.3, and integrated the channels that had emission above 3σ.We assumed local thermodynamical equilibrium, used an excitation temperature of 20 K, a 13 CO/ 12 CO abundance ratio of 69 (Wilson 1999) and a 12 CO/H 2 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 K or 50 K, and found 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.

Outflow detection
Here we explain the process we followed to determine whether a FUor had an outflow detection.We inspected all 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 12 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 12 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 outwards to higher velocities.The envelope emission dominates the velocities closest to the systemic, which are approximately the same velocities covered by the 13 CO emission (Table 2), so we focused on velocities beyond these.If there emission is not detected at velocities higher the ones overran 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 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 them to be an outflow, and in the case 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 non-detections with the only difference begin V900 Mon.In this case, the line profiles do not show high velocity wings and thus resulted in a non-detection 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 in , of the outflow.Next, we located the velocity channel where a 3σ detection was not found and considered it as the "outer" velocity, v out .We then calculated the maximum velocity of each lobe of the outflow with respect to the systemic velocity as v max = v out − v sys .The list of FUors with outflows and these three velocities are 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 12 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 in 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 12 CO, and their contour maps are presented in Figure 4 and Figure 5.We used these maps to estimate the position angle of the outflow, reported in Table 3, and to estimate the extension of each lobe (see below).

Outflow properties
One of the goals of this work is to compare the outflows emanating from eruptive young stars to those from quiescent young stellar objects.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 Table 4 and Table 5, respectively.

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 −4 with respect to H 2 (Bolatto et al. 2013).We calculated the mass (M v ) for each velocity channel (v) for all pixels above 3σ.Afterwards, we calculated the momentum and kinematic energy for each channel with P v = M v ×v and E v = 0.5M v ×v 2 , respectively.Finally, we integrated the three properties over the same velocity range to obtain the total values (M of , P of , E of ) for the blueshifted and redshifted lobes of the outflow.

Force and luminosity
The outflow force and luminosity are calculated as F of = P of /τ d and L of = E of /τ d , respectively, where τ d is the dynamical time of the outflow.The dynamical time is defined as τ d = R lobe /v max , where R lobe is the projected extension of the outflow lobe and v max is the maximum velocity of the outflow.
For the projected extension of the outflows, we used the integrated emission maps (Figure 4 and Figure 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 independently by calculating the difference between the systemic velocity and the minimum/maximum velocity where there is blueshifted/redshifted emission.In both R lobe and v max , the sensitivity and spatial resolution of the observations directly affect their measured values.

Caveats
Due to the nature of our observations and methodology, it is important to understand the limitations of our estimated values.
First, the large uncertainties 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 into 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 put 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 −1 , while we obtained ∼13.1 km s −1 (see Table 3).Their JCMT observations had typical noise values 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 that 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 Figure 4 and Figure 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.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 lobe and τ d are lower limits.The † labels the two FUors with tentative outflow detections.

Optical depth correction
We calculated the outflow parameters assuming the 12 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 13 CO emission, under the assumption that that isotopologue is optically thin, and use it to correct the fluxes of 12 CO.In this section we present our methodology to correct the emission of the J=3-2 transition of the 12 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 where T mb,12 and T mb,13 are the brightness temperatures of 12 CO and 13 CO, respectively, and τ 12 and τ 13 are their respective opacities.Assuming that 13 CO is optically thin, Equation 1 can be re-written as where [ 12 CO]/[ 13 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 (−τ 12 ))/τ 12 by calculating T mb,12 /T mb,13 for each channel where both isotopologues were detected above 6σ.In some low-velocity channels for a few FUors, the 13 CO appears to be optically thick, therefore, we dropped these points from the fitting.
We then fitted a parabola which will allow us to correct for the velocity channels where the 13 CO emission was not detected.We fixed B = 0 to keep the parabola symmetric with respect to the systemic velocity and prevent over-correcting 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 the 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 13 CO, we only used blueshifted points to fit the parabola (see Appendix B).

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 higher content of mass in the envelope surrounding the FUor and, thus, an indication 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 young stellar objects inside the beam of our observations, and 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 based on the Moment 0 maps, we found that the direction of the outflow is aligned with the line of sight, and, 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 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 and thus the inclination of 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.
Haro 5a IRS -This is a Class I protostar is located in the Orion star forming region and it was identified as a FUor-like object by Reipurth et al. (2012).Previous CO observations of source releaved its outflow (Takahashi et al. 2006(Takahashi et al. , 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 Table 6.Outflow properties from 12 CO (3-2) observations after optical depth correction.
masses are higher than ours by less than a factor of 2, which can be explained by the difference in distances and excitation temperatures, and by the difference in 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 VLA, and both reported that this FUor is also a proto-binary star.
Reipurth 50 -A Class I protostar 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 C 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 its is likely they have resolved out most of the extended emis-sion 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 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 C 18 O with ALMA (Zurlo et al. 2017).
The authors found 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 at redshifted velocities (+3 km s −1 , see also Figure 3).Zurlo et al. (2017) reported different velocity ranges for 12 CO and C 18 O (see their Table 2).The velocities of the 12 CO match those of the redshifted excess emission (peaking at ∼6 km s −1 , see Figure 3), and the velocities of C 18 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 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 complicate 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 on 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 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 P Cygni 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 12 CO line profiles show high-velocity wings (Figure 3), which we interpreted as an indication of a bipolar molecular outflow.The integrated emission maps (Figure 4 and Figure 5) and the 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 re-solved 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.
Z CMa -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 13 CO, and in the J=2-1 and J=1-0 transitions of 12 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 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 the redshifted lobe is more extended than the blueshifted side.Our estimation of 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 the other outflow properties are comparable to those by Mottram et al. (2017).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 12 CO/ 13 CO ratio used in the optical depth correction, it appears even the rarer isotopologue is optically thick at velocities close to the systemic.

FUors with tentative detections
Below we present the two FUors for which we can only make a tentative detection of their outflows, and thus we consider that these two sources require followup 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 as us and carried out a similar analysis as us, and did not find outflow emission.However, Takami et al. (2019) presented high-angular resolution ALMA observations of the J=2-1 transition of 12 CO, 13 CO and C 18 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 following the one 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, and 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 a ∼2.5 mag brightening in 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 which is ongoing as of this writing4 .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 sub-millimeter wavelength observations of this object.Its CO line profiles (Figure 3) show high velocity line wings, in particular on the redshifted side.Based on its J=3-2 integrated emission maps (Figure 4) and channels maps (Appendix A), there appears to be an outflow whose blueshifted and redshifted lobes are on 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, and thus, it is hard to confirm that the morphology seen in the J=3-2 transition as an outflow.

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 13 CO J=3-2 was a combination of the outflow at low velocities and of 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 at dominates in the lowvelocity channels of our observations.A similar case was found for FU Ori, for which previous observations reported it did not drive an outflow (Levreault 1988) but ALMA observations of J=2-1 hint towards 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 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.

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 Class 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% 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), in the images (i.e.pixel size and field of view), and by using the same excitation for the two transitions in all 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.

Comparison with quiescent young stellar objects
We put into context our outflow properties 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 a hundred years and the dynamical ages of the outflows are in the other 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 towards the nature behind the outbursts.
We compared our sample with the values of the outflow properties published in the following studies: Dun-   2014), Yıldız et al. (2015) and Mottram et al. (2017).The 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 are Class I objects (see Table 1) because we want 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), Mottram et al. (2017) used the same temperature 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 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 are at later stages (Class I or II).We estimated the inclination of the outflows by assuming the inclination of the outflow is 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 Z CMa we used the estimate by Antoniucci et al. (2016) estimated 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 correcion 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 sub-samples with the FUors.
In Figure 8 we plotted the outflow forces calculated from the J=3-2 transition against the envelope masses calculated from the 13 CO emission (panel a), against the outflow masses also calculated from the 12 CO J=3-2 transition (panel b), against the ratio between the outflow mass and envelope mass (panel c), and against the bolometric luminosities obtained from the literature (panel d; Table 1).
Panels a and b of Figure 8 show that FUors outflows follow the same trends as the quiescent sources from the literature, i.e. higher envelope masses and higher outflow masses indicate higher outflow forces.In the outflow mass 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 Figure 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 with 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 ksamp, 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 Is and different from the Class 0s (Figure 9, panel a), the outflow masses of FUors are different to those of Figure 9. Cumulative histograms of outflow properties, and the outflow mass to envelope mass ratio, for the Class 0 and Class 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.Class I objects and comparable with those of Class 0s (Figure 9, panel b), the outflow mass to envelope mass ratio is different in the FUors when compared to either sample (Figure 9, panel c), and the outflow forces of FUors are comparable with the two quiescent samples (Figure 9, panel 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 Is.Following the same logic, the result for the outflow masses is surprising as the outflow masses are close to the Class 0s.This could be interpreted as an indication that FUors are in the very early stages of their Class I stage and their outflows have not had sufficient time to dissipate.However, the masses of 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.FUor outflows have higher outflow masses but similar envelope masses compared to the Class Is, pointing towards 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 mass and the envelope mass in Figure 8 and Figure 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 withing 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 towards the FUors as it is likely that for we are measuring the full extension of their outflows in comparison to the quiescent sample.In Figure 10 we plotted 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 their outflows extend beyond the areas of the sky covered by their respective observations.Therefore, the mass ratios for those sources are lower limits, and thus it raises the question 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 calculate the core-to-star formation efficiency, , as follows: where τ 0+I is the total duration of the Class 0 and I phases, i.e. 0.5 Myr.The authors used their typical values of M of /M env = 10 −2 and τ d = 10 4 years, 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 for the FUors.The values can be found in Figure 7.
As can be seen from our results, five of the FUors have negative values.These can be explained by the underestimation of the dynamical age because either the outflow appears to be face-on (V2775 Ori and V960 Mon), or it 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) and thus their 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 and with 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 outwards 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 .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 for Class IIs (Fiorellino et al. 2022) 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 to the accretion history of a young stellar object (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 lower 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 year old outburst could be detected.Few FUor outbursts (∼20 Connelley & Reipurth 2018) have been detected so the incidence rate of these events is unknown, and as such it is difficult to separate between 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 outflow forces 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 those of others 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.

SUMMARY & 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.2M show the silicate feature at 10 µ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 µm feature in emission despite having an envelope of ∼0.6 M .
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 and, thus, the higher transition is dimmer.
V960 Mon is an outlier in both trends, thus 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 those 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.

APPENDIX
A. CHANNEL MAPS Here we present the channel maps for the three observed transitions of L1551 IRS 5, 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 channels maps are those when the gas emission starts or finished being significant.The purple contours are used for all the 13 CO channel maps and for the 12 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 13 CO, and the outflow properties in the case of both 12 CO transition.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.In Table B1 we present the parameters of the best-fitted parabola used to determine the optical depth correction for the sources with outflows, and in Figure B1 we present the parabolic fit and the line profiles used in the fitting., 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. CO (4-3) with contours at 3,5,7,9,11,13,15,17,19

Figure 1 .
Figure1.Integrated intensity (Moment 0) maps of our targets for the 12 CO (3-2) line observed with APEX (orange contours).The Moment 0 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 contours of the 850 µm continuum emission from the JCMT (see Section 4.1).The CO and the dust contours are plotted with levels at 0.3, 0.4, . . ., 0.9 of the peak intensity, and are meant to be representative.The star symbols indicate the nominal positions of the protostars.The bars at the bottom of each panel represent 10 000 au.

Figure 3 .
Figure 3. 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.

V346
Nor -Kóspál et al. (2017a)  presented single dish observations of the J=3-2 and J=4-3 transitions, andKóspál et al. (2017b) presented ALMA Cycle 2 observations of the J=2-1 transition.Our analysis uses the same observations asKóspál et al. (2017a) and the properties of outflow have the same values within 10%.

Figure 7 .
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.

Figure 8 .
Figure 8.Comparison of the outflows from quiescent sources in the literature with the FUor outflows.Outflow forces (J=3-2) plotted against envelope masses (panel a), outflow masses (J=3-2, panel b), the 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 empty diamond symbols.The L bol of iPTF 15afq is unknown and thus it is not shown in panel d.

Figure 10 .
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.

Figure B1 .
Figure B1.Optical depth correction for L1551 IRS 5.In the left panel we show the line profiles, where the green, purple and black colors represent the 13 CO, the observed 12 CO and the corrected 12 CO, respectively, and the vertical dashed lines indicate the rage of velocities used in the parabola fit.The right panel shows the ratio of main beam temperatures (TMB), where the light blue crosses are all the values of the ratio for each velocity channel, in pink dots with errorbars are the points used in the parabolic fitting, and the black line is the resulting best-fitted parabola.The complete figure set (10 images) is available in the online journal.

Table 1 .
Objects observed.FUor classification vLSR b Distance c L bol d

Table 2 .
Envelope masses of the FUors based on 13 CO (3-2).vmin and vmax 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, and the reference used for each target.

Table 3 .
Position angles, velocities, extensions and dynamical times of outflows

Table 7 .
Ratio between outflow masses and envelope masses using the J=3-2 transition of the two observe CO isotopologues, and the core-t-star formation efficiency, .

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.

Table B1 .
Parameters of best fitted parabolas used for optical depth correction.