CHRONOLOGY OF EPISODIC ACCRETION IN PROTOSTARS - AN ALMA SURVEY OF THE CO AND H2O
SNOWLINES
Tien-Hao Hsieh1, Nadia M. Murillo2, Arnaud Belloche3, Naomi Hirano4, Catherine Walsh5, Ewine F. van Dishoeck2,6, Jes K., Jørgensen7, Shih-Ping Lai1,8
1Institute of Astronomy and Astrophysics, Academia Sinica, P.O. Box 23-141, Taipei 106, Taiwan 2Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA, Leiden, the Netherlands
3Max-Planck-Institut f¨ur Radioastronomie, Auf dem H¨ugel 69, 53121 Bonn, Germany 4Institute of Astronomy and Astrophysics, Academia Sinica, P.O. Box 23-141, Taipei 106, Taiwan
5School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK
6Max-Planck-Institut f¨ur extraterrestrische Physik, Giessenbachstraße 1, 85748, Garching bei M¨unchen, Germany 7Niels Bohr Institute, University of Copenhagen, Øster Voldgade 5–7, DK 1350 Copenhagen K., Denmark and
8Institute of Astronomy, National Tsing Hua University (NTHU), Hsinchu 30013, Taiwan Draft version September 9, 2019
ABSTRACT
Episodic accretion has been used to explain the wide range of protostellar luminosities, but its origin and influence on the star forming process are not yet fully understood. We present an ALMA survey of N2H+(1 − 0) and HCO+ (3 − 2) toward 39 Class 0 and Class I sources in the Perseus molecular cloud. N2H
+
and HCO+are destroyed via gas-phase reactions with CO and H2O, respectively, thus tracing
the CO and H2O snowline locations. A snowline location at a much larger radius than that expected from the current luminosity suggests that an accretion burst has occurred in the past which has shifted the snowline outward. We identified 18/18 Class 0 and 9/10 Class I post-burst sources from N2H
+
, and 7/17 Class 0 and 1/8 Class I post-burst sources from HCO+. The accretion luminosities during the past bursts are found to be ∼ 10 − 100 L . This result can be interpreted as either evolution of
burst frequency or disk evolution. In the former case, assuming that refreeze-out timescales are 1000 yr for H2O and 10,000 yr for CO, we found that the intervals between bursts increases from 2400 yr in the Class 0 to 8000 yr in the Class I stage. This decrease in the burst frequency may reflect that fragmentation is more likely to occur at an earlier evolutionary stage when the young stellar object is more prone to instability.
Subject headings: Star formation, Interstellar medium, Protostars, Astrochemistry
1. INTRODUCTION
Episodic accretion plays an important role in star for-mation (Audard et al. 2014). In the episodic accretion scenario, a protostellar system stays in a quiescent ac-cretion phase most of the time, and acac-cretion bursts oc-casionally occur to deliver material onto the central pro-tostar. Because the accretion luminosity dominates the stellar luminosity at the early embedded phase (Hart-mann & Kenyon 1996), this behavior leads to a low pro-tostellar luminosity for the majority of the time. Such low luminosities have been revealed by recent surveys in star-forming regions with large statistical samples (Evans et al. 2009; Enoch et al. 2009; Kryukova et al. 2012; Hsieh & Lai 2013; Dunham et al. 2014), yet it is still un-known how the accretion luminosity evolves and affects the luminosity distribution (Offner & McKee 2011).
The origin of episodic accretion is still unknown due to the difficulty in directly observing the accretion process. The most plausible explanation is disk instability which can originate from several mechanisms such as thermal instability (Lin et al. 1985; Bell & Lin 1994), gravita-tional instability (Vorobyov & Basu 2005, 2010; Boley & Durisen 2008; Machida et al. 2011), magnetorotational instability (Armitage et al. 2001; Zhu et al. 2009; Zhu et al. 2010a,b). Stellar (or planet) encounters have also been proposed to explain accretion bursts (Clarke & Syer 1996; Lodato & Clarke 2004; Forgan & Rice 2010). Furthermore, Padoan et al. (2014) suggest that
the mass accretion rate is controlled by the mass infall from the large-scale turbulent cloud. These possibilities make episodic accretion a key mechanism in star forma-tion because it is associated with the timeline of disk fragmentation, planet formation, and multiplicity. Ob-servationally, Liu et al. (2016) and Takami et al. (2018) found large-scale arms and arc structures in the disks to-ward four out of five FU Ori-type objects (Herbig 1966, 1977, see below) with Lbol∼ 100−600 L . This supports
the hypothesis that bursts are triggered by fragmentation due to gravitational instabilities. However, the remain-ing source, V1515 Cyg, hosts a smooth and symmetric disk. From an evolutionary point of view, Vorobyov & Basu (2015) find that the outburst should preferentially occur during the Class I stage after the disk has accreted sufficient material to fragment. However, Hsieh et al. (2018) found that accretion bursts with a few to a few tens of L have occurred in Very Low Luminosity
Ob-jects (VeLLOs), which are extremely young or very-low mass protostars and thus unlikely hosts of massive disks. Episodic accretion alters the star formation process by regulating the radiative feedback (Offner et al. 2009). The change in the thermal structure of the disk can directly affect the chemical composition of gas and ice (Cieza et al. 2016; Wiebe et al. 2019). For exam-ple, Taquet et al. (2016) found that complex organic molecules could be formed via gas-phase reactions in the hot region (T &100 K) during the outburst phase. In a continuous accretion process, the radiative feedback
could suppress fragmentation by keeping the disk and/or cloud core warm (Offner et al. 2009; Yıldız et al. 2012, 2015; Krumholz et al. 2014). On the contrary, episodic accretion can moderate this effect, and during the qui-escent phase, the disk has sufficient time to cool down and fragment (Stamatellos et al. 2012). Such a process can be associated with the formation of binary/multiple systems and the formation of substellar objects, affecting the multiplicity and initial mass function (Kratter et al. 2010; Stamatellos et al. 2011; Mercer & Stamatellos 2016; Riaz et al. 2018). Therefore, it is crucial to reveal the time intervals and the magnitude of outbursts in or-der to study how this radiative feedback affects the star formation process.
Variations in protostellar luminosity have been found in the past decades (i.e., FU Orionis and EX Orionis events; Herbig 1966, 1977), which are considered to arise directly from episodic accretion. Given the time intervals between bursts (∼ 5×103−5×104yr, Scholz et al. 2013),
there are only a handful of cases in which luminosity vari-ability has been reported to date (V1647 Ori: ´Abrah´am et al. 2004; Andrews et al. 2004; Acosta-Pulido et al. 2007; Fedele et al. 2007; Aspin et al. 2009, OO Serpentis: K´osp´al et al. 2007, [CTF93]216-2: Caratti o Garatti et al. 2011, VSX J205126.1: Covey et al. 2011; K´osp´al et al. 2011, HOPS 383: Safron et al. 2015). An outburst has also been detected toward the high-mass star forming region, S255IR-SMA1 (Caratti o Garatti et al. 2016; Liu, S. et al. 2018). However, among these sources, HOPS 383 is the only source at an early stage (near the end of the Class I phase) due to the difficulty of infrared/optical observations to probe the embedded phase (Safron et al. 2015). At longer wavelengths, Liu, H. et al. (2018) found variations in millimeter flux of 30 − 60% toward 2 out of 29 sources using SMA. The James Clerk Maxwell Telescope (JCMT) transient sur-vey monitored 237 sources at submillimeter wavelengths for 18 months (Herczeg et al. 2017; Johnstone et al. 2018), and they identified only one burst, from the Class I source, EC53.
Chemical tracers are sensitive to the thermal history and can be used to probe past luminosity outbursts over a much longer timescale than direct observations of the luminosity (Lee 2007; Kim et al. 2011, 2012; Visser & Bergin 2012; Visser et al. 2015). Jørgensen et al. (2013) found a ring-like structure of H13CO+
surround-ing the extended CH3OH emission toward IRAS 15398-3359. They found that this anti-correlation highlights the H2O snowline location because CH3OH has a sub-limation temperature similar to that of H2O and H2O can destroy H13CO+ via gas-phase reactions (Visser et al. 2015). Given the current luminosity, they suggested that IRAS 15398-3359 has experienced a past luminos-ity outburst, sublimating H2O over a larger region. This result has later been confirmed with HDO by Bjerkeli et al. (2016) that directly revealed the radial extent of H2O emission. Such an anti-correlation between H2O and H13CO+ has also been found by van ’t Hoff et al. (2018b) in NGC 1333-IRAS2A, supporting that H13CO+ is a good tracer of the water snowline. An alternative method is to look directly at the CO snowline, which is at larger distances and can be more easily resolved us-ing extended C18O emission. Jørgensen et al. (2015)
and Frimann et al. (2017) studied 16 and 24 embed-ded protostars, respectively, and found that 20 − 50 % of the Class 0/I sources have experienced recent accretion bursts. Assuming a CO refreeze-out time of ∼10,000 yr, they estimated that the time interval between accretion bursts is 2−5×104yr. Later, Hsieh et al. (2018) derived
the CO snowline radius using the spatial anti-correlation of CO and N2H+ in VeLLOs, and found that 5 out of 7
sources are post-burst sources. On the other hand, us-ing CO and N2H+ to trace the snowline, Anderl et al.
(2016) found no evidence for past luminosity outbursts in four Class 0 protostars.
Here we present a survey of HCO+ and N2H +
, chem-ical tracers of past accretion bursts, in 39 protostars in Perseus that include 22 Class 0 sources and 17 Class I sources. In Section 2 we describe the sample and the ob-servations. The observational results are shown in Sec-tion 3, and the detailed analysis and the modeling are given in Section 4. Finally, we discuss the implications of our findings on episodic accretion, and summarize these results in Sections 5 and 6, respectively.
2. OBSERVATIONS
2.1. Sample
We selected 39 protostars from Dunham et al. (2015) in Perseus (d = 293 pc, Ortiz-Le´on et al. 2018), a star-forming region containing sufficient Class 0 and Class I sources for a statistical survey. These targets are lo-cated in the Western Perseus region (Hsieh & Lai 2013) mostly near the NGC1333 and B1 regions. Our sam-ple includes 22 Class 0 and 17 Class I protostars with Tbol = 25 − 490 K (Table 1). The bolometric
luminosi-ties were taken from Dunham et al. (2015) and were scaled with the new measured distance (250 pc → 293 pc), which yields Lbol = 0.1 − 45 L . These 39 sources
are selected because they are detected at 850 or 1120 µm continuum emission by single dish observations (COM-PLETE survey, Ridge et al. 2006); the presence of the continuum emission suggests that the envelope has not yet dissipated, which might be a marker of stronger line emission. All targets are included in the sample of the “Mass Assembly of Stellar Systems and their Evolution with the SMA (MASSES)” survey (Lee et al. 2015, 2016; Stephens et al. 2018), and we use the naming convention Per-emb-XX denoted by Enoch et al. (2009).
2.2. N2H+ (1 − 0) observation
We observed the N2H+ (1 − 0) line emission toward 36
out of the 39 targets from 2018 March to 2018 April using ALMA (Cycle 5 project, 2017.1.01693.S, PI: T. Hsieh). The N2H+(1−0) data for the remaining three targets are
taken from an earlier ALMA project (2015.1.01576.S), and the results of which were reported in Hsieh et al. (2018). With the C43-4 configuration, the resulting beam size was ∼ 2.004 × 1.005 using natural weighting. The largest scale covered is ∼ 1200. The channel width was 30 kHz (∼0.1 km s−1 at a frequency of 93 GHz). The on-source time toward each on-source was ∼ 7.5 min, resulting in an rms noise level of ∼ 13 mJy beam−1 at a spec-tral resolution of 0.1 km s−1. The gain calibrator was
TABLE 1 Targets
Source Other name R.A. Dec F1.2 mm Tbol Lbol P.A. Reference
(J2000) (J2000) mJy T L degree Per-emb-2 IRAS 03292 + 3029 03h32m17.92s +30d49m47.85s 702.0±34.8 25 1.8 127 1,2,3 Per-emb-3 03h29m00.58s +31d12m00.17s 52.6±1.1 30 0.9 277 1,4 Per-emb-4 DCE065 03h28m39.11s +31d06m01.66s 0.9±0.2 28 0.3 ∗50 1,5 Per-emb-5 IRAS 03282 + 3035 03h31m20.94s +30d45m30.24s 279.4±4.5 32 1.6 125 1,2,6 Per-emb-6 DCE092 03h33m14.41s +31d07m10.69s 10.7±0.4 34 0.9 53 1,7 Per-emb-7 DCE081 03h30m32.70s +30d26m26.47s 8.1±1.0 34 0.2 165 1,5
Per-emb-9 IRAS 03267 + 3128, Perseus5 03h29m51.83s +31d39m05.85s 11.1±1.3 39 0.7 63 1
Per-emb-10 03h33m16.43s +31d06m52.01s 21.9±0.6 26 1.4 230 1 Per-emb-14 NGC 1333 IRAS4C 03h29m13.55s +31d13m58.10s 97.6±1.7 35 1.2 95 1,8 Per-emb-15 RNO15-FIR 03h29m04.06s +31d14m46.21s 6.3±0.9 17 0.9 145 1,4 Per-emb-19 DCE078 03h29m23.50s +31d33m29.12s 16.8±0.4 60 0.5 335 1,7 Per-emb-20 L1455-IRS4 03h27m43.28s +30d12m28.78s 9.9±1.5 54 2.3 295 1 Per-emb-22 L1448-IRS2 03h25m22.41s +30d45m13.21s 51.4±5.0 52 2.7 318 1,2,8 03h25m22.36s +30d45m13.12s - - - -Per-emb-24 03h28m45.30s +31d05m41.66s 3.9±0.3 62 0.6 281 1,8 Per-emb-25 03h26m37.51s +30d15m27.79s 120.5±1.5 64 1.2 290 1,9,10 Per-emb-27 NGC 1333 IRAS2A 03h28m55.57s +31d14m36.98s 247.6±12.4 54 30.2 204 1,2,4 03h28m55.57s +31d14m36.42s - - - -Per-emb-29 B1-c 03h33m17.88s +31d09m31.78s 133.1±5.5 41 4.8 110 1 Per-emb-30 03h33m27.31s +31d07m10.13s 47.9±0.8 62 1.8 109 1,11 Per-emb-31 DCE064 03h28m32.55s +31d11m05.04s 2.1±0.4 52 0.4 345 1,7 Per-emb-34 03h30m15.17s +30d23m49.19s 11.4±0.5 93 1.9 45 1,10
Per-emb-35 NGC 1333 IRAS1, Per-emb-35A 03h28m37.09s +31d13m30.76s 27.8±1.1 100 13.0 290 1,2,6
Per-emb-35B 03h28m37.22s +31d13m31.73s - - - -Per-emb-36 NGC 1333 IRAS2B 03h28m57.38s +31d14m15.74s 154.6±3.0 100 7.3 204 1,2,4 Per-emb-38 DCE090 03h32m29.20s +31d02m40.75s 26.0±0.7 120 0.7 250 1,7 Per-emb-39 03h33m13.82s +31d20m05.11s 2.0±0.4 59 0.1 - 1 Per-emb-40 B1-a 03h33m16.67s +31d07m54.87s 16.9±0.5 100 2.2 280 1,2 Per-emb-41 B1-b 03h33m20.34s +31d07m21.32s 11.2±0.4 47 0.8 210 1,6 B1-bS 03h33m21.36s +31d07m26.37s - - -
-Per-emb-44 SVS 13A, Per-emb-44-B 03h29m03.75s +31d16m03.77s 313.6±22.1 170 45.3 130 1,6
Per-emb-44-A 03h29m03.77s +31d16m03.78s - - - -SVS 13A2 03h29m03.39s +31d16m01.58s - - - -SVS 13B 03h29m03.08s +31d15m51.70s - - - -Per-emb-45 03h33m09.58s +31d05m30.94s 1.4±0.2 210 0.1 - 1 Per-emb-46 03h28m00.42s +30d08m00.97s 4.3±0.4 230 0.3 315 1 Per-emb-48 L1455-FIR2 03h27m38.28s +30d13m58.52s 4.0±0.4 260 1.1 295 1 Per-emb-49 Per-emb-49-A 03h29m12.96s +31d18m14.25s 21.9±1.9 240 1.4 207 1,6 Per-emb-49-B 03h29m12.98s +31d18m14.34s - - - -Per-emb-51 03h28m34.51s +31d07m05.25s 85.4±4.3 150 0.2 110 1 Per-emb-52 03h28m39.70s +31d17m31.84s 4.3±0.4 250 0.2 25 1 Per-emb-54 NGC 1333 IRAS6 03h29m01.55s +31d20m20.48s 3.1±0.4 230 11.3 310 1 Per-emb-58 03h28m58.43s +31d22m17.42s 4.6±0.2 240 1.3 5 1 Per-emb-59 03h28m35.06s +30d20m09.44s 1.6±0.1 49 0.5 - 1 Per-emb-63 03h28m43.27s +31d17m32.90s 24.6±0.5 490 2.2 ∗20 1,9 03h28m43.36s +31d17m32.69s - - - -03h28m43.57s +31d17m36.31s - - - -Per-emb-64 03h33m12.85s +31d21m24.00s 39.7±0.6 480 4.0 ∗70 1 Per-emb-65 03h28m56.32s +31d22m27.75s 35.6±0.7 440 0.2 140 1
Note. — The source coordinates are obtained by a Gaussian fitting for the 1.2 mm images and the fluxes are listed in the next column (Figure 1). The bolometric temperature (Tbol) and luminosity (Lbol) are taken from Dunham et al. (2015), in which Lbol is scaled to the new measured distance of Perseus (250 pc → 293 pc). We use Tbol = 70 K (Evans et al. 2009) as a boundary to classify Class 0 and Class I sources. The position angle (P.A.) of the outflow axis is taken from the corresponding references.
∗ The presumed P.A. of the source comes from disk or envelope structures rather than outflows.
References: (1) Stephens et al. (2018),(2) Tobin et al. (2016),(3) Schnee et al. (2012),(4) Plunkett et al. (2013),(5) Hsieh et al. (2018),(6) Lee et al. (2016),(7) Hsieh et al. (2017),(8) Tobin et al. (2015),(9) Segura-Cox et al. (2018),(10) Dunham et al. (2014),(11) Davis et al. (2008)
2.3. HCO+ (3 − 2) and 1.2 mm continuum observations
The HCO+ (3 − 2) data were taken simultaneously
with the continuum and CH3OH (20,2− 1−1,1) data
to-ward the 39 targets in 2018 September during the same project (2017.1.01693.S). The integration time is ∼ 12 min for each source. The array configuration was C43-5, resulting in a beam size of 0.0045 × 0.0030 using natural weighting and 0.0033 × 0.0022 using uniform weighting; for each source, we choose the weighting based on the S/N of the image: uniform weighting is used for better detected sources. The largest scale covered is ∼ 2.006. The chan-nel widths for both HCO+ and CH3OH were 30.5 kHz
(∼ 0.03 km s−1) and were averaged to 0.1 km s−1 when imaging. The rms noise level at a spectral resolution of
0.1 km s−1 is ∼ 8 − 13 mJy beam−1 depending on the weighting. The window for continuum emission was cen-tered at 268 GHz with a bandwidth of ∼ 1.85 GHz. For all executions, the flux and band pass calibrators were J0237+2848, and the phase calibrator was J0336+3218.
3. RESULTS
3.1. Continuum emission at 1.2 mm
Per-emb-2 1" Per-emb-3 1" Per-emb-4 1" Per-emb-5 1" Per-emb-6 1" No vandam Per-emb-7 1" Per-emb-9 1" Per-emb-10 1" Per-emb-14 1" Per-emb-15 1" Per-emb-19 1" Per-emb-20 1" Per-emb-22 1" Per-emb-24 1" Per-emb-25 1" Per-emb-27 1" Per-emb-29 1" Per-emb-30 1" Per-emb-31 1" Per-emb-34 1" Per-emb-35 1" Per-emb-36 1" Per-emb-38 1" Per-emb-39 1" Per-emb-40 1" Per-emb-41 1" Per-emb-44 1" Per-emb-45 1" Per-emb-46 1" Per-emb-48 1" Per-emb-49 1" Per-emb-51 1" Per-emb-52 1" Per-emb-54 1" Per-emb-58 1" Per-emb-59 1" Per-emb-63 1" Per-emb-64 1" Per-emb-65 1"
1.2 mm Continuum
Per-emb-2 Per-emb-3 Per-emb-4 Per-emb-5 Per-emb-6 Per-emb-7
Per-emb-9 Per-emb-10 Per-emb-14 Per-emb-15 Per-emb-19 Per-emb-20
Per-emb-22 Per-emb-24 Per-emb-25 Per-emb-27 Per-emb-29 Per-emb-30
Per-emb-31 Per-emb-34 Per-emb-35 Per-emb-36 Per-emb-38 Per-emb-39
Per-emb-40 Per-emb-41 Per-emb-44 Per-emb-45 Per-emb-46 Per-emb-48
Per-emb-49 Per-emb-51 Per-emb-52 Per-emb-54 Per-emb-58 Per-emb-59
Per-emb-63 Per-emb-64 Per-emb-65
HCO
+3
−2
and N
2
H
+1
−0
sources at 1.2 mm as the system centers. The compan-ions of Per-emb-40 and 48 found by the VLA Nascent Disk and Multiplicity Survey (VANDAM) at 8 mm (To-bin et al. 2016) are not detected at 1.2 mm probably due to insufficient sensitivity.
3.2. N2H+ maps
Figure 2 presents the N2H+ (1 − 0) integrated
inten-sity maps and the velocity ranges over which they were integrated are listed in Table 2. In order to maximize the signal-to-noise ratio, we integrated the emission from all seven hyperfine components in N2H+ (1 − 0). All
targets are detected except for Per-emb-38, 45, 58, 59, and 64. For Per-emb-41, the N2H
+
emission is likely associated with B1b-S located in the north east. The N2H+emission in the Per-emb-39, 49, and 65 maps likely traces background cloud structures rather than the en-velopes of the sources. These non-detected and ambigu-ous sources are thus excluded in the following analysis. As a result, 19 Class 0 and 11 Class I sources are re-maining. In most of the sources, the N2H
+
emission peaks are offset from the continuum source, and are anti-correlated with HCO+ which is also shown in Figure 2. This suggests that the N2H+ emission is suppressed in the warmer regions where it is destroyed through reac-tions with CO sublimating off dust grains (Mauersberger & Henkel 1991; Jørgensen et al. 2004; van ’t Hoff et al. 2017). Several targets show strong negative contours in the map, likely caused by the spatial filtering of large-scale emission. This is not expected to affect our analy-sis in which we measure the peak position of the N2H+ emission and compare it with model predictions.
3.3. HCO+ maps
The blow-ups of the HCO+ integrated intensity maps are shown in Figure A1. We integrated the line emission excluding the optically thick region near the systemic ve-locity (Table 2, see Appendix A). This effect of optical depth is discussed in section 4.4.1. The velocity ranges of integration are listed in Table 2. We removed the fol-lowing sources in the analysis because the emission does not seem to reflect the H2O snowline in the envelope: (1) Per-emb-4, 41, 49, and 59 show no detection near the source, (2) Per-emb-52 shows weak emission and ambigu-ous structures, and (3) For Per-emb-2, 36, 51, 64 and 65, the HCO+ emission peaks near the outflows axes, which are likely associated with the outflows rather than the envelopes. We discuss if and how exclusion of these targets might introduce a bias. Most of these targets (except for Per-emb-4, see section 4.1 and Hsieh et al. 2018) in categories (1) and (2) have no N2H+ detec-tions, and in category (3), Per-emb-36, 64, and 65 also show no (or ambiguous) N2H+detections. This suggests that their envelopes are almost dissipated; for Per-emb-36, 64, and 65 in category (3), outflows might dominate the emission with low envelope densities. Since the dissi-pation is likely associated with the evolution rather than one accretion outburst, we speculate that the removal of these sources does not add a selection bias. However, this narrows down our sampling evolutionary phase in the more evolved Class I stage with Tbol & 250 K.
Per-emb-2 (Class 0) and Per-emb-51 (Class I) from category
(3) can introduce a bias but the small number only af-fects the statistical results by 5% for the Class 0 stage and by 11% for the Class I stage given the number of remaining sources. We note that in order to minimize the outflow contamination, our analysis focusses on the emission roughly along the axis perpendicular to the out-flow. As a result, 18 Class 0 and 11 Class I sources are remaining for the following analysis.
3.4. CH3OH (20,2− 1−1,1) maps
CH3OH (20,2 − 1−1,1) at 254.015377 GHz is detected
toward the source center in six targets (Figure 3). The emission most likely traces the region where CH3OH
sublimates due to central heating. All these detected sources show an anti-correlation between CH3OH and HCO+ except for Per-emb-20 and 22 which have very weak CH3OH emission. Because CH3OH shares a similar
sublimation temperature with H2O (∼ 100 K, Collings
et al. 2004), these anti-correlations could be used to con-firm the radii of the H2O snowlines (van ’t Hoff et al.
2018a). For Per-emb-20 and 22, CH3OH and HCO
+
emission have a similar peak position at the center (see also Appendix C), which likely comes from unresolved structures at the current resolution. However, the non-detections of CH3OH do not indicate a temperature less than ∼ 100 K due to the unknown CH3OH abundance
and probably the low Einstein coefficient of 1.9×10−5s−1 for the transition. The origin and the presence or absence of CH3OH emission will be discussed in a separate paper
(Murillo et al., in prep.).
4. ANALYSIS
In order to identify the post-burst sources, we model the line emission at different central luminosities. We compare the N2H+ and HCO+ peak positions derived from the integrated intensity maps with those from the models. Here we describe how we measured the peak positions from the observed images (Section 4.1) and how we construct the models (Section 4.2). Then, we discuss the identification of sources that have likely experienced a past burst, i.e., post-burst sources, in section 4.3, and the caveats in section 4.4.
4.1. The peak radii in the integrated intensity maps The difficulty in measuring the radius of the emission peak is that the observed peak position is not always located at the equatorial plane and such a plane is not necessarily perpendicular to the outflow. These misalign-ments can be reproduced by simulations (Offner et al. 2016) and are widely seen in observations (Lee et al. 2016; Hsieh et al. 2016; Stephens et al. 2018). To derive the radius of the peak emission, we employ a biconical mask to filter out the outflow contaminated region; depend-ing on sources, the region with a position angle smaller than 35 − 85◦ from the outflow axis is excluded. We identify the local maxima from the resulting maps as the peak position of the emission on one or two sides. The mask in use and the selection of peak introduce artificial effects; taking the N2H
+
TABLE 2 N2H
+
and HCO+integrated intensities
N2H+ HCO+
source beam vel. range rms beam vel. range rms
00 km s−1 mJy beam−1km s−1 00 km s−1 mJy beam−1km s−1
Per-emb-2 2.5 × 1.5 5.9 − 7.9 25.4 0.33 × 0.22 3.9 − 6.1, 7.9 − 9.3 10.0 Per-emb-3 2.5 × 1.5 6.2 − 8.0 24.1 0.33 × 0.22 5.1 − 6.4, 7.8 − 9.5 9.6 Per-emb-4 3.2 × 2.0 6.7 − 7.3 9.8 0.46 × 0.30 5.9 − 8.0 5.4 Per-emb-5 2.5 × 1.5 6.5 − 7.3 21.3 0.34 × 0.22 5.0 − 6.4, 8.2 − 9.5 9.1 Per-emb-6 2.5 × 1.5 5.6 − 6.2 16.3 0.43 × 0.29 5.1 − 5.4, 7.8 − 8.2 3.2 Per-emb-7 3.1 × 2.0 5.9 − 6.4 14.3 0.33 × 0.22 3.1 − 5.5, 6.6 − 8.4 10.8 Per-emb-9 2.5 × 1.5 7.7 − 8.3 17.8 0.34 × 0.22 5.0 − 7.6, 8.6 − 10.2 11.0 Per-emb-10 2.5 × 1.5 6.2 − 6.9 19.3 0.34 × 0.22 4.6 − 5.6, 7.4 − 8.7 8.4 Per-emb-14 2.5 × 1.5 6.5 − 8.7 25.3 0.34 × 0.22 5.0 − 6.7, 8.7 − 10.3 9.6 Per-emb-15 2.5 × 1.5 6.3 − 7.0 19.3 0.38 × 0.25 2.6 − 5.9, 8.7 − 9.0 7.8 Per-emb-19 2.5 × 1.5 7.2 − 7.7 18.5 0.38 × 0.25 6.0 − 6.9, 8.0 − 8.8 5.7 Per-emb-20 2.4 × 1.5 4.4 − 5.2 21.0 0.33 × 0.22 2.5 − 4.1, 6.2 − 8.3 10.6 Per-emb-22 2.4 × 1.5 3.5 − 4.7 23.0 0.33 × 0.22 1.2 − 3.3, 4.9 − 7.9 12.7 Per-emb-24 2.5 × 1.5 7.0 − 7.5 17.6 0.38 × 0.25 5.6 − 6.6, 8.1 − 9.1 6.1 Per-emb-25 2.3 × 1.5 4.7 − 5.6 20.8 0.33 × 0.22 2.7 − 4.6, 6.1 − 8.1 11.1 Per-emb-27 2.5 × 1.5 6.4 − 8.3 25.5 0.34 × 0.22 4.7 − 5.9, 8.8 − 10.8 9.6 Per-emb-29 2.5 × 1.5 5.4 − 7.2 24.4 0.34 × 0.22 3.9 − 5.7, 7.3 − 8.9 9.9 Per-emb-30 2.5 × 1.5 6.6 − 7.4 19.9 0.34 × 0.22 −0.1 − 6.1, 7.5 − 15.0 22.2 Per-emb-31 3.2 × 2.0 6.7 − 7.6 19.0 0.38 × 0.25 5.5 − 6.4, 7.9 − 9.4 7.0 Per-emb-34 2.4 × 1.5 5.6 − 6.4 21.3 0.33 × 0.22 1.6 − 5.1, 6.7 − 9.6 14.5 Per-emb-35 2.5 × 1.5 6.8 − 7.4 19.2 0.34 × 0.22 4.9 − 6.6, 8.0 − 9.6 10.5 Per-emb-36 2.5 × 1.5 7.0 − 7.6 18.3 0.34 × 0.22 1.1 − 5.8, 9.3 − 10.9 13.6 Per-emb-38 2.5 × 1.5 6.3 − 6.7 16.1 0.43 × 0.29 5.6 − 8.3 6.4 Per-emb-39 2.5 × 1.5 6.6 − 7.3 18.4 0.46 × 0.30 6.3 − 6.7 2.4 Per-emb-40 2.5 × 1.5 5.7 − 7.2 22.6 0.34 × 0.22 1.5 − 5.0, 7.8 − 11.5 14.8 Per-emb-41 2.5 × 1.5 6.4 − 6.9 17.4 0.46 × 0.30 5.8 − 7.6 5.2 Per-emb-44 2.5 × 1.5 7.7 − 9.1 23.3 0.34 × 0.22 6.6 − 8.3, 9.3 − 10.8 9.9 Per-emb-45 2.5 × 1.5 6.7 − 7.3 17.8 0.43 × 0.29 5.5 − 7.7 5.7 Per-emb-46 2.4 × 1.5 4.7 − 5.3 18.1 0.33 × 0.22 3.3 − 4.4, 5.3 − 5.8 7.0 Per-emb-48 2.4 × 1.5 3.9 − 5.2 22.1 0.43 × 0.28 1.2 − 3.5, 5.0 − 5.7 6.8 Per-emb-49 2.5 × 1.5 8.1 − 8.7 17.7 0.46 × 0.30 7.8 − 9.5 4.9 Per-emb-51 2.4 × 1.5 6.5 − 7.2 20.0 0.34 × 0.22 5.1 − 6.2, 7.1 − 7.9 7.7 Per-emb-52 2.5 × 1.5 7.6 − 8.2 18.1 0.46 × 0.30 6.9 − 7.7 3.4 Per-emb-54 2.5 × 1.5 7.7 − 8.4 19.5 0.38 × 0.25 3.0 − 7.6, 9.0 − 14.0 14.0 Per-emb-58 2.5 × 1.5 7.6 − 8.4 19.8 0.34 × 0.22 4.4 − 6.4, 9.1 − 10.5 10.4 Per-emb-59 2.4 × 1.5 5.2 − 5.8 18.7 0.45 × 0.30 5.6 − 7.6 5.7 Per-emb-63 2.5 × 1.5 7.5 − 8.3 19.7 0.38 × 0.25 6.9 − 8.6 5.6 Per-emb-64 2.5 × 1.5 4.1 − 4.6 17.5 0.34 × 0.22 −3.1 − 6.1, 7.7 − 15.0 26.2 Per-emb-65 2.5 × 1.5 8.4 − 9.2 21.0 0.46 × 0.30 6.5 − 9.1 6.2
identified, the uncertainty is taken as the half-beam size. If two peaks are found, the peak radius is defined as the average of their distances to the primary source and the uncertainty is the taken as the difference of that but with a minimum value as the half-beam size (Figures 2, A1, and Table 3).
We check using the images in Figure 2 if the N2H+ peaks are located beyond the HCO+ emitting region; HCO+ is indeed expected to form from CO and traces the region where CO returns to the gas phase. We find that most of the HCO+emitting areas lie reasonably in-side the N2H
+
peak radii, at least along the major axis of the source. For several sources, the elongated HCO+ emission is likely tracing the outflow and thus extends beyond the N2H+ peak radius (e.g. Per-emb-5, 7, 9, 22, 27, 29, 48, 51, and 54). This result suggests that the measured N2H+peaks reflect the CO snowline radii. Per-emb-4 and 52 show N2H+ depletion without detections of HCO+ toward the center. Per-emb-52 has extended C18O emission toward the center (Hsieh, T., in prep.). On the other hand, the CO isotopologues, 13CO, C18O, and C17O (1 − 0) and (2 − 1) are not (or marginally) detected in Per-emb-4 (Hsieh et al. 2018, DCE065 in the paper). Thus, for Per-emb-4, we cannot exclude the
possibility that N2H+ is absent due to freeze out of N2 in the central dense region (Belloche & Andr´e 2004).
We check if the HCO+ peaks are located beyond the CH3OH emitting regions for the six sources with CH3OH detections (Figure 3). The spatial extent of CH3OH emission is broadly within the radius of the measured HCO+ peak. This suggests that HCO+ reflects well the H2O sublimation region and that the measured radii are reasonable. Unfortunately, most of the sources have no CH3OH detection; such non-detections do not rule out the hypothesis that HCO+ probes the location of the H2O snowline but prevent us to confirm it.
Figure 4 shows the intensity profiles of N2H +
Per-emb-20 200 au Per-emb-22 200 au Per-emb-27 200 au Per-emb-29 200 au Per-emb-35 200 au Per-emb-44 200 au
CH
3OH 2
0, 21
1, 1Fig. 3.— CH3OH integrated intensity map (blue contours) over-laid on that of HCO+ (green scale). The contour levels are 3σ, 5σ, 10σ, 30σ, and 70σ. The purple circles showing the radii of the measured HCO+ peak are the same as described in Figure A1.
4.2. MHW19 models
To help the interpretation of the observational data, we construct a model framework for studying the rela-tion between the luminosity and the emission peak radii of N2H+ and HCO+ (Murillo, N. M, in prep.,
here-after MHW19). MHW19 have built a grid of density
structures by varying the disk and envelope geometry. For each density structure, RADMC3D1 (Dullemond et
al. 2012) is employed to calculate the temperature pro-file with a given central luminosity from 0.01 to 200 L . These physical conditions are later used for
deriv-ing the molecular abundance profiles. The model starts from initially icy H2O and CO, and a chemical network based on the UMIST database for Astrochemistry, ver-sion RATE12 (McElroy et al. 2013) is used to calcu-late the static chemical abundance profiles. The abun-dances of CO, N2 and H2O relative to H2 are set to 2 × 10−4, 1 × 10−4, and 2 × 10−4, respectively. A con-stant cosmic-ray ionization rate of appropriate for the interstellar medium, 1.2 × 10−17 s−1, is adopted. It is noteworthy that Padovani et al. (2016) found that the
1 http://www.ita.uni-heidelberg.de/∼dullemond/software/radmc-3d/
shocked gas from protostellar jets enhances the ioniza-tion rate by a few orders of magnitude, which can affect the abundances of HCO+and N2H+(Gaches et al. 2019). However, this effect is less significant in the midplane of a flattened envelope which is the region in which we are interested.
Due to central heating, CO and H2O are sublimated from the dust grains and destroy N2H+ and HCO+ via the gas-phase reactions
CO + N2H+−−→ HCO++ N2 and
H2O + HCO+−−→ H3O++ CO,
respectively (see Section 1). Thus, N2H+and HCO+are expected to be depleted in the central regions where the temperature is higher than the sublimation temperatures of CO (∼20 K) and H2O (∼100 K), respectively. Finally, MHW19 make line-emission images using the line radia-tive transfer code available as part of RADMC3D us-ing molecular transition data from the Leiden database (Sch¨oier et al. 2005) (N2H
+
from Green 1975 and HCO+ from Flower 1999). These images perform as a reference for comparison with observations with different burst lu-minosities Lburst.
In this work, we adopt the model with a rotationally flattened envelope given by Ulrich (1976), i.e., with-out a disk component, from MHW19. We include an outflow cavity with the edge, z(r) following the func-tion z(r) ∝ r1.5with an opening angle of 50◦at z = 50 au
(Whitney et al. 2003; Robitaille et al. 2006). For a quan-titative comparison between the model and observation, we measure the molecular peak radii from the models with different luminosities at different inclination angles (see Appendix B). Here we show an example comparing the observed N2H
+
and HCO+integrated intensity maps of Per-emb-5 to that of the model (Figure 5); the modeled image is generated using the CASA task “simobserve” for a source with a luminosity of 30 L at an inclination
an-gle of 45◦. Figure 6 compares the modeled and observed molecular peak radii as a function of luminosity. The modeled peak radius is affected by the inclination due to the emission contributed from the inner or outer envelope (see Appendix B). Besides, because a disk could shield the outer region from the central radiative heating, the modeled luminosity should be considered as a lower limit; if a disk exists, a higher luminosity is needed to shift the peak position outward to match the observation (see the disk model in Figure 6).
4.3. Identification of post-burst sources
10 0
10
0.0
0.5
1.0
I/I
peakPer-emb-44
L
bol= 45.3 L
burst
Per-emb-44
L
bol= 45.3 L
burst
1
0
1
offset (arcsec)
0.0
0.5
1.0
I/I
peak10 0
10
Per-emb-29
L
bol= 4.8 L
< 1000 yr
Per-emb-29
L
bol= 4.8 L
< 1000 yr
1
0
1
offset (arcsec)
10 0
10
Per-emb-20
L
bol= 2.3 L
1000 10, 000 yr
Per-emb-20
L
bol= 2.3 L
1000 10, 000 yr
1
0
1
offset (arcsec)
10 0
10
Per-emb-51
L
bol= 0.2 L
> 10,000 yr
Per-emb-51
L
bol= 0.2 L
> 10,000 yr
N2H+ HCO+1
0
1
offset (arcsec)
CH3OH HCO+Fig. 4.— Intensity profiles along cuts across the source center and the local minima from the integrated intensity maps toward four sources. The top panel shows the N2H+and smoothed HCO+profiles at large scales, and the bottom panel shows the HCO+and CH3OH profiles at small scales. The profiles are normalized to their peaks. The bolometric luminosity of the target is labeled in the upper left corner. The purple area indicates the radius of the measured peak for N2H+ (top) and HCO+(bottom) in Table 3. The vertical and the horizontal dashed lines represent the source position and the intensity zero level, respectively.
the peak luminosity in the past, i.e., the burst luminos-ity (Lburst, Table 3), and identify the post-burst sources
with Lburst > Lbol. However, this estimated luminosity
is degenerate with the inclination angle for the case of HCO+ at small scales (Figure 6). This degeneracy be-comes severe near the pole-on case (see Appendix B). Fortunately, most of our targets show clear bipolar out-flows (Stephens et al. 2018), suggesting that they are not pole-on sources. Statistically, the nearly pole-on proba-bility (θinc< 25◦) is less than 10%. Thus, we derive the
burst luminosity by comparing the measured peak radius to the model at an inclination angle of 45◦ and use the angle from 25◦− 75◦ as the uncertainty (Table 3).
As a result, we found that with N2H+, almost all Class 0 and Class I sources are identified as post-burst sources. With HCO+, 10/17 Class 0 sources and 2/10 Class I sources are identified as post-burst sources. The sources with a peak radius less than the half-beam size are not classified as post-burst sources; they are classi-fied as sources without a past burst or sources where CO or H2O have refrozen onto the dust grains after the last burst.
It is noteworthy that Per-emb-4 has no detection of HCO+ nor CO isotopologues (Hsieh et al. 2018) in the N2H+ depletion region. Thus, we are not able to exclude the possibility that the N2H+ depletion comes from freeze out of its parent molecule, N2.
4.4. Caveat in identification of the post-burst sources 4.4.1. Effect of optical depth for probing the H2O snowline
The optical depth from the continuum or line emis-sion could affect the measured peak positions. If the line emission and dust continuum emission are optically thin, the integrated intensity maps should properly re-flect the snowline locations. N2H+ is expected to be
less affected by this issue, because it is usually optically thin throughout the outer envelope due to its relatively low abundance. However, for the H2O snowline traced by HCO+ in the inner dense region, the effects of opti-cal depth need to be addressed. Here we discuss how it influences the measured snowline radii in two cases, opti-cally thick HCO+ line emission and optically thick dust continuum emission:
(a) If the HCO+ emission is optically thick, it prevents us to probe the inner dense region. An HCO+ hole would be seen in the HCO+integrated intensity map due to line self-absorption and/or continuum sub-traction. In order to reduce this effect, we integrated the spectra avoiding the optically thick regions near the systemic velocity (Table 2) (see Appendix A). Excluding the low-velocity channels should not af-fect the measured snowline location because it is ex-pected to locate at the inner region where the ve-locity is high. However, it is unclear what fraction of the continuum emission is absorbed by the fore-ground HCO+ gas at the selected velocity ranges for integration, resulting in an “over subtraction” in the continuum subtraction process. This issue leads to a mis-identification of HCO+ depletion, e.g., as for the case of Per-emb-49 in Figure A1. The current data cannot completely rule out this possibility except for sources with CH3OH detections (see below).
TABLE 3
Identification of the past burst from N2H +
and HCO+
name Lbol RN2H+ peak Lburst,CO M˙acc RHCO+ peak Lburst,H2O M˙acc Last burst
L 00 L 10−6M yr−1 00 L 10−6M yr−1 yr Per-emb-2 1.8 4.8 ± 0.9 18.1+1.9−0.1 6.9+0.7−0.1 - - - < 10 000 Per-emb-3 0.9 1.7 ± 0.9 3.7+0.1−0.3 1.4−0.1+0.1 0.16 ± 0.13 1.6+1.1−1.0 0.6+0.4−0.4 < 1000 Per-emb-4 0.3 4.9 ± 1.3 20.0+0.1−1.9 7.7+0.1−0.7 - - - ∗< 10 000 Per-emb-5 1.6 5.2 ± 1.0 22.1+0.1−2.1 8.5−0.8+0.1 0.35 ± 0.14 49.2+17.2−22.2 18.9+6.6−8.5 < 1000 Per-emb-6 0.9 5.0 ± 1.0 20.0+0.1−0.1 7.7−0.1+0.1 0.19 ± 0.18 6.7+4.3−5.7 2.6+1.6−2.2 < 1000 Per-emb-7 0.2 2.0 ± 1.3 4.5+0.1−0.1 1.7+0.1−0.1 0.04 ± 0.13 <0.9 <0.3 1000 − 10 000 Per-emb-9 0.7 1.1 ± 1.0 2.2+0.1−0.2 0.8+0.1−0.1 0.06 ± 0.14 <0.9 <0.3 1000 − 10 000 Per-emb-10 1.4 3.8 ± 1.0 12.1+1.3−0.1 4.6+0.5−0.1 0.10 ± 0.14 <0.9 <0.3 1000 − 10 000 Per-emb-14 1.2 2.2 ± 1.0 4.9+0.1−0.1 1.9−0.1+0.1 0.31 ± 0.14 40.3+14.1−22.2 15.4+5.4−8.5 < 1000 Per-emb-15 0.9 5.7 ± 0.9 24.4+0.1−0.1 9.4+0.1−0.1 0.03 ± 0.15 <1.6 <0.6 1000 − 10 000 Per-emb-19 0.5 5.9 ± 1.0 27.0+0.1−0.1 10.4+0.1−0.1 0.14 ± 0.16 <1.8 <0.7 1000 − 10 000 Per-emb-20 2.3 4.4 ± 0.9 16.4+0.1−1.6 6.3+0.1−0.6 0.09 ± 0.14 <0.9 <0.3 1000 − 10 000 Per-emb-22 2.7 7.1 ± 0.9 36.4+0.1−0.1 14.0−0.1+0.1 0.20 ± 0.14 8.1+6.7−6.9 3.1+2.6−2.6 < 1000 Per-emb-24 0.6 7.1 ± 0.9 36.4+0.1−0.1 14.0−0.1+0.1 0.26 ± 0.15 24.4+12.0−19.5 9.4+4.6−7.5 < 1000 Per-emb-25 1.2 3.1 ± 0.9 9.0+0.1−0.9 3.5−0.3+0.1 0.29 ± 0.13 36.4+12.8−24.3 14.0+4.9−9.3 < 1000 Per-emb-27 30.2 5.5 ± 1.0 24.4+0.1−2.3 9.4−0.9+0.1 0.48 ± 0.14 99.1+0.1−25.7 38.0+0.1−9.9 0 Per-emb-29 4.8 3.1 ± 1.0 9.0+0.1−0.9 3.5−0.3+0.1 0.33 ± 0.14 44.5+15.6−22.4 17.1+6.0−8.6 < 1000 Per-emb-30 1.8 2.0 ± 1.0 4.5+0.1−0.4 1.7+0.1−0.2 0.13 ± 0.14 <0.9 <0.3 1000 − 10 000 Per-emb-31 0.4 4.6 ± 1.3 18.1+0.1 −1.7 6.9 +0.1 −0.7 0.34 ± 0.15 49.2 +17.2 −24.8 18.9 +6.6 −9.5 < 1000 Per-emb-34 1.9 2.0 ± 0.9 4.0+0.4−0.1 1.5−0.1+0.2 0.18 ± 0.14 5.4+2.7−4.5 2.1+1.0−1.7 < 1000 Per-emb-35 13.0 7.3 ± 0.9 40.3+0.1 −0.1 15.4 +0.1 −0.1 0.16 ± 0.14 1.6 +1.1 −1.0 0.6 +0.4 −0.4 1000 − 10 000 Per-emb-36 7.3 12.4 ± 0.9 109.5+0.1−0.1 42.0+0.1−0.1 - - - < 10 000 Per-emb-38 0.7 - - - 0.05 ± 0.18 <6.7 <2.6 > 1000 Per-emb-39 0.1 - - - 0.23 ± 0.19 13.4+11.0 −11.6 5.1+4.2−4.4 < 1000 Per-emb-40 2.2 3.6 ± 1.0 11.0+1.2−0.1 4.2+0.5−0.1 0.08 ± 0.14 <0.9 <0.3 1000 − 10 000 Per-emb-41 0.8 - - - -Per-emb-44 45.3 6.6 ± 1.0 33.0+0.1 −0.1 12.7+0.1−0.1 0.41 ± 0.14 73.4+7.7−24.2 28.1+3.0−9.3 0 Per-emb-45 0.1 - - - 0.07 ± 0.18 <6.7 <2.6 > 1000 Per-emb-46 0.3 2.4 ± 0.9 5.4+0.1−0.1 2.1+0.1−0.1 0.12 ± 0.13 <0.9 <0.3 1000 − 10 000 Per-emb-48 1.1 1.3 ± 0.9 2.7+0.3−0.3 1.0−0.1+0.1 0.39 ± 0.17 66.4+14.7−26.1 25.5+5.6−10.0 < 1000 Per-emb-49 1.4 - - - -Per-emb-51 0.2 0.3 ± 1.0 <2.2 <0.8 - - - > 10 000 Per-emb-52 0.2 4.2 ± 0.9 14.8+0.1−1.4 5.7+0.1−0.5 - - - < 10 000 Per-emb-54 11.3 9.2 ± 0.9 60.1+6.3−0.1 23.0+2.4−0.1 0.11 ± 0.15 <1.6 <0.6 1000 − 10 000 Per-emb-58 1.3 - - - 0.05 ± 0.14 <0.9 <0.3 > 1000 Per-emb-59 0.5 - - - -Per-emb-63 2.2 5.1 ± 0.9 20.0+2.1−0.1 7.7+0.8−0.1 0.15 ± 0.15 <1.6 <0.6 1000 − 10 000 Per-emb-64 4.0 - - - -Per-emb-65 0.2 - - -
-Note. — Col. (1): Source name. Col. (2): Bolometric luminosity from Dunham et al. (2014) scaled from 230 pc to 293 pc. Col. (3): Radius of the measured N2H+ peak. Col. (4): Luminosity corresponding to the radius of the measured N2 H+ peak from models, i.e., the accretion luminosity during the past burst. Col. (5): Mass accretion rate estimated from Col. (4). Col. (6) − (8): Same as Col. (3) − (5) but with the numbers measured from HCO+ . Col. (9): Time after the last burst. ∗ Per-emb-4 has no HCO+ and CO detection toward the center. Thus, it is unclear if the N2H+ depletion comes from destruction by CO or freeze-out of the parent molecule, N2.
to effective destructive collisions and higher dust den-sities (Banzatti et al. 2015; Cieza et al. 2016). If this is the case, the HCO+ depletion toward the center may still reflect the H2O snowline locations assum-ing that the optically thick dust region and an associ-ated HCO+depleted region are due to a dust opacity change within the water snowline.
Since CH3OH has a sublimation temperature similar
to that of H2O (∼ 100 K), CH3OH line emission can be used as a proxy for the location of the H2O snowline. For those sources with CH3OH detections, the measured HCO+ peak radii broadly agree with the CH3OH emis-sion extents (Figure 3), which can resolve the issues of optical depths mentioned above. This indicates that the estimates of the H2O snowline locations from HCO+ are reasonable at least for these six sources. Thus, we
spec-ulate that HCO+ is a good tracer of the H2O snowline, but future observations with a resolution sufficient to re-solve the continuum source or more warm-gas tracers are required to completely rule out the optical-depth issue.
4.4.2. Dependence of the physical and chemical models
The binding energy used in the chemical model deter-mines the sublimation temperature, the decisive param-eter for the snowline locations (Collings et al. 2003). The binding energy of pure CO is found between ∼ 850−1000 K (Bisschop et al. 2006) but can be increased to ∼ 1200 − 1700 K depending on the substrate (Fayolle et al. 2016), resulting in a CO sublimation temperature of ∼ 17 − 33 K at a gas density of ∼ 107 cm−3. To make a
Observation
N
2H
+(1 0)
Per-emb-5
2000 au
Model
2000 au
HCO
+(3 2)
Per-emb-5
100 au
100 au
10
110
210
310
4R (au)
10
1510
1310
1110
910
710
510
3[X
]/[
H
2]
CO gas
N
2H
+H2O gas
HCO
+Fig. 5.— An example showing the comparison of the observa-tions and models for Per-emb-5. The top row shows the observed N2H+ map (left) and HCO+map (right) in the same contours as in Figures 2 and A1. These images are at different sizes with scale bars shown in the bottom right corner. The images of a model with 30 L are shown in the second row with arbitrary intensity scales and contours to emphasize the emission peak. The purple circles indicate the radii of the N2H+and HCO+peaks in the correspond-ing map from observations. The bottom plot shows the molecular abundance profiles in the equatorial plane for the model.
water ice) and 4820 K for H2O (Sandford & Allamondola 1993; Fraser et al. 2001). Figure 6 shows the modeled curve with the CO binding energy of 1150 K (Collings et al. 2004, Tsub∼ 21 K) and 1307 K (Tsub∼ 25 K).
How-ever, the binding energy is degenerate with the density structures (see below), and the latter are not necessarily the same in different sources.
Although a massive unstable disk is presumed to trig-ger the accretion burst, we perform our analysis using models without a disk. A protostellar disk can shield the envelope or itself from the central radiation, chang-ing the temperature structure mainly along the equa-torial plane (Murillo et al. 2015). Figure 6 (bottom) shows the modeled curve from a disk model in compar-ison with the no-disk model that is used for the HCO+
10
110
010
110
210
310
4R
NH2 +pe ak(a
u)
Class 0
Class I
VeLLO (Hsieh+)
Per27
Per44
T
sub21
K
T
sub25
K
10
110
010
1L
bol(L )
10
110
2R
HCO +pe ak(a
u)
Class 0
Class I
Per27
Per44
Disk
No disk
Fig. 6.— Measured N2H+(top) and HCO+peak radii (bottom) as a function of the bolometric luminosities of the sources. The peak positions measured in this work are shown by open circles (Class 0) and filled triangles (Class I), and that from Hsieh et al. (2018) are in red. Per-emb-27 and 44 have been classified as sources undergoing an accretion burst and are plotted in black. We note that the peak position in Hsieh et al. (2018) is measured using abundance profiles, which could result in a slightly higher value. The modeled peak radius as a function of luminosity is shown with the blue lines for an inclination angle θinc= 45◦. The dark and light blue area represent a range of θinc= 25 − 65◦and θinc= 15 − 75◦, respectively, for which the upper boundary has a smaller inclination angle (0◦ for pole-on). The dashed black line shows a model at an inclination angle of 45◦with a CO sublimation temperature of 21 K (top) and with the presence of a disk (bottom). The grey area corresponds to a half-beam size, the limitation to resolve the depletion in the observations.
peak radii. A higher central luminosity is required to heat the envelope and shift the H2O snowline outward for the disk model compared with the no-disk model. If a massive disk is included in the model, the num-ber of post-burst sources and the burst luminosities are expected to increase. However, a disk model is much more complicated since the disk density, geometry, and grain size distribution all affect the temperature struc-ture. To keep the model simple, we chose to use the no-disk model. The no-no-disk model could be considered as a conservative approach for identifying the occurrence of a past burst. Specifically, the modeled snowline radii are an upper limit, and correspondingly the estimated burst luminosities, Lburst,H2O and Lburst,CO, are lower limits.
the temperature structures and snowline locations.
5. DISCUSSION
5.1. Sources in the burst phase
Our sample includes two sources, Per-emb-27 (NGC 1333 IRAS2A) and Per-emb-44 (SVS 13A), which are likely undergoing an accretion burst given the current huge Lbol= 30 L and Lbol= 45 L , respectively.
Several complex organic molecules are detected toward Per-emb-27 in the inner 40−100 au region where the dust temperature is > 100 K (Maret et al. 2014; Maury et al. 2014; Taquet et al. 2015). Codella et al. (2014) found a knotty jet driven by Per-emb-27 with a dynamical time of < 30 − 90 yr, which is considered to be a signature of episodic accretion (Vorobyov et al. 2018). Per-emb-44 contains also a central hot region with a detection of glycolaldehyde (De Simone et al. 2017; Bianchi et al. 2019). Multiple components are found in continuum ob-servations by Tobin et al. (2016, 2018), including a close binary with a separation of 0.003 (70 au). Furthermore, Lef`evre et al. (2017) speculate the existence of a compan-ion with a separatcompan-ion of 20−30 au in order to explain the knotty jets with a period of ∼ 300 yr, a putative compan-ion that triggers the burst episodes in the close perihelcompan-ion approach with an eccentric orbit. These results suggest that Per-emb-27 and 44 are in the accretion-burst phase. These burst-phase sources can be used as calibrators for the model. If we assume that the current bolometric luminosity determines the observed CO and H2O peak radii, the modeled curve in Figure 6 should go through the data points of Per-emb-27 and 44. As a result, the adopted model without a disk component and a CO sub-limation temperature of ∼ 25 K looks reasonable (see Section 4.4.2). However, the assumption is not necessar-ily true because the luminosity variation during a burst can be very large (Elbakyan et al. 2016; Vorobyov et al. 2018). Furthermore, the density structures should be quite different for each source. Thus, this calibration provides only a rough confirmation.
5.2. Chronology of episodic accretion in protostars In this section we discuss the history of the episodic accretion process in a statistical way from Class 0 to Class I. We derive frequencies of accretion bursts in Sec-tion 5.2.1 and mass accreSec-tion rates during burst phases in Section 5.2.2. Then, based on these results, we discuss the mass accumulation history of protostars in episodic accretion in Section 5.2.3
5.2.1. Evolution of burst frequency
Estimation of the outburst frequency has been done by monitoring a sample of protostars in the more evolved Class II stage (Scholz et al. 2013; Contreras Pe˜na et al. 2019). These studies require a large survey with a long baseline in time (Hillenbrand, & Findeisen 2015). The chemical probes extend the baseline to 1000-10,000 yr by considering the refreeze-out time scales.
The refreeze-out time scales of CO and H2O are
differ-ent because their snowlines are located at differdiffer-ent radii with different densities. Thus, these two chemical trac-ers provide complementary information to constrain the time since a past burst. The refreeze-out time scales can
10
010
110
2L
burst ,C O/L
bolPer27
Per44
Burst rate: 27/28 (96±3%)
Class 0: 18/18 (100±0%)
Class I: 9/10 (90±9%)
Class 0
Class I
10
110
2T
bol(K)
10
110
010
110
2L
burst ,H2 O/L
bolPer27
Per44
Burst rate: 12/27 (44±24%)
Class 0: 10/17 (58±24%)
Class I: 2/10 (20±16%)
Fig. 7.— Lburst/Lbol as a function of the bolometric tempera-ture Tbolfrom N2H+ (top, Lburst,CO/Lbol) and HCO+ (bottom,
Lburst,H2O/Lbol). The markers shown in this plot is the same as
that in Figure 6 for Class 0 and Class I sources. The vertical dashed lines indicate the bolometric temperature of 70 K, the boundary between Class 0 and Class I. The horizontal dashed lines represent
Lburst> Lbolfor classification of post-burst sources.
be expressed as a function of gas density (nH2) and dust
temperature (Tdust), i.e.,
τfr = 1 × 104 yr r 10 K Tdust 106 cm−3 nH2 (1)
from Visser & Bergin (2012) and Visser et al. (2015). Here we assume τfr is 10,000 yr for CO and 1000 yr for
H2O (Visser et al. 2015). As a result, we categorize
these targets into: (1) post-burst source from HCO+. The burst has occurred in the past 1000 yr. (2) post-burst source from N2H+but not from HCO+. The burst has occurred during the past 1000 to 10,000 yr. (3) no signature of a burst neither in N2H+ nor HCO+. No burst has occurred during the past 10,000 yr (Table 3). Figure 4 shows the intensity profiles of the standard cases for these three categories together with that during an outburst. Table 3 lists the time since the last burst for each source.
Figure 7 shows Lburst/Lbol as a function of the
evo-lutionary indicator, bolometric temperature (Tbol), from
both CO and H2O. Sources with Lburst> Lbol are
experienced a past burst within the refreeze-out time. Excluding the two sources in the burst phase (Per-emb-27 and 44, section 5.1), there are 28 and (Per-emb-27 sources for the following statistical analyses with N2H+and HCO+, respectively. Lburst is defined as an upper limit for those
sources whose measured peak radius is less than the beam size, and the upper limit is obtained using the half-beam size as the peak radius. As a result, we cannot identify the chemical signature of a past burst for Per-emb-51 on the basis of its N2H
+
map, and for seven sources on the basis of their HCO+ maps. For Per-emb-51 with N2H+, the upper limit of Lburst,CO, 2.2 L ,
sug-gests that it unlikely experienced a past burst, or at least a strong burst. For those seven sources from Lburst,H2O,
depending on the weighting of the map, the upper limits of Lburst,H2Oare 0.9 L for three sources, 1.8 L for two
sources, and 6.7 L for two sources. These upper
lim-its are generally smaller than the burst luminosities of the post-burst sources, a median of 31.1 L and a
stan-dard deviation of 20.8 L (Figure 8). Therefore, we do
not consider these sources to be post-burst sources. As a result, we find 100 ± 0% of Class 0 objects are burst sources and 90 ± 9% of Class I objects are post-burst sources from N2H+ alone, where the uncertainty is derived using binomial statistics. This result implies that the burst interval is <10,000 yr for Class 0 objects and ∼11,000 yr (10,0000.9 yr) for Class I objects. From HCO+ alone, we identify 58 ± 24% of Class 0 objects and 20 ± 16% of Class I objects are post-burst sources. This gives us burst intervals of ∼1700 yr for Class 0 ob-jects and ∼5,000 yr for Class I obob-jects. The combination of these results suggests that the burst frequency is de-creasing from the Class 0 to the Class I stage.
The criterion used to define post-burst sources, Lburst > Lbol, may not well portray sources during
an accretion burst especially for the estimate of burst frequencies; intuitively, large bursts occur rarely com-pared with small bursts. For example, several post-burst sources (Per-emb-3, 9, 30, etc.) have their Lburst,CO of
∼ 2.2 − 4.5 L , only ∼ 3 − 4 times larger than their
Lbol. These sources may have experienced a small
ac-cretion outburst ( ˙Macc ∼ (1 − 2) × 10−6 M yr−1, see
section 5.2.2) that is expected to occur more frequently. In addition, the observed bolometric luminosity can be affected by the viewing angle of a disk-outflow system. A small Lbolmay result from a nearly edge-on configuration
(Offner et al. 2012). This can lead to a misidentification of a post-burst source if Lbol is very small and if Lbust
is only slightly larger. Thus, we decide to focus on those post-burst sources robustly identified. If we consider only large outbursts with Lburst > Lbol and Lburst > 10 L
(Enoch et al. 2009), 56 ± 24% of Class 0 objects and 50 ± 25% of Class I objects are post-burst sources from N2H
+
, implying an interval of ∼18,000 yr and ∼20,000 yr for Class 0 and Class I, respectively (Figure 8). From HCO+, there are 41 ± 24% of post-burst sources in Class 0 with an interval of ∼2,400 yr and 12 ± 10% of post-burst sources in Class I with an interval of ∼8,000 yr for Class I (note that two sources with Lbol > 10 L
are excluded given the new criterion). The inconsisten-cies between the burst intervals traced with N2H+ and HCO+is difficult to explain but this result still suggests
a decrease of burst frequency from the Class 0 stage to the Class I stage. However, from an evolutionary point of view, if a disk is significantly denser/larger at the Class I stage than that at the Class 0 stage, it might shrink the emission peak inward (Figure 6). The disk evolution may thus lead us to underestimate the number of post-burst sources at the Class I stage, and our conclusion that the burst frequency decreases from the Class 0 to the Class I stage might in turn not be robust.
Based on N2H+ observations, Hsieh et al. (2018) found that episodic accretion can start at a very early evolutionary stage. Here we find that the burst fre-quency is higher in the Class 0 stage than in the Class I stage. The accretion outbursts are believed to be associ-ated with a massive and large disk with the gravitational and/or magnetorotational instability (Vorobyov & Basu 2015; Zhu et al. 2010b). Therefore, the burst frequency may reflect the disk formation and/or evolution. The onset of disk formation is still unclear, but disks have been found in some Class 0 sources (Tobin et al. 2012; Murillo & Lai 2013; Ohashi et al. 2014; Lee et al. 2017, 2018; Aso et al. 2017; Hsieh et al. 2019; Maury et al. 2019). Besides, numerical simulations suggest that an initially unstable cloud core can promote disk formation and tend to have a higher mass accretion rate (Machida et al. 2016). Vorobyov & Basu (2013) further support that the episodic accretion process is highly dependent on the core initial conditions; with a higher ratio of rota-tional to gravitarota-tional energy, the strength of the burst is increased.
If accretion bursts are triggered by infalling fragments in an gravitationally unstable disk (Vorobyov & Basu 2005), the decreasing burst-frequency implies that, at an earlier stage, either disk fragmentation occurs more fre-quently or that the fragments tend to fall more often onto the central source. It is also noteworthy that Reg´aly & Vorobyov (2017) find that the gravitational instability in Vorobyov & Basu (2005) is overestimated with a fixed central source because the disk angular momentum can translate into the orbital motion of the central source. The gravitational instability can be controlled by infall onto the disk and its thermodynamics (Kratter & Lodato 2016). A high mass infall rate is crucial for sustaining the gravitational instability in disks (Vorobyov & Basu 2005; Kratter et al. 2008), which might explain the high burst frequency in the Class 0 stage. In addition, frag-mentation is suggested to occur in cold regions which re-quire a sufficient cooling time (Vorobyov & Basu 2010; Kratter et al. 2010; Kratter & Murray-Clay 2011). Ob-servations of multiple systems (Murillo et al. 2016; Tobin et al. 2018) in the cold disk/envelope support such frag-mentations at an early stage.
5.2.2. Mass accretion rate
The derived burst luminosities give us indications on the mass accretion rate during the past burst phase. If we assume that Lburst = Lacc, the accretion luminosity,
we derive the mass accretion rate with Lacc=
GMstarM˙acc
R (2)
10
010
110
2L
burst ,C O(L
)
Per27
Per44
Burst rate: 14/25 (56±24%)
Class 0: 10/17 (58±24%)
Class I: 4/8 (50±25%)
Class 0
Class I
10
110
2T
bol(K)
10
010
110
2L
burst ,H2 O(L
)
Per27
Per44
Burst rate: 8/25 (32±21%)
Class 0: 7/17 (41±24%)
Class I: 1/8 (12±10%)
10
610
510
4M
acc(M
yr
1)
10
610
510
4M
acc(M
yr
1)
Fig. 8.— Luminosity during the past burst obtained from N2H+ (top) and HCO+ (bottom) as a function of bolometric tempera-ture. The luminosity is derived by modeling the line-emission peak offset, and sources with an offset less than the half-beam are la-beled with the half-beam as an upper limit. The horizontal dashed lines indicate the burst luminosity Lburst> 10 L . The right axis shows the mass accretion rate derived assuming that the luminosity is dominated by accretion luminosity.
and Mstar is the mass of the central source (assumed
to be 0.25 M , Evans et al. 2009, half of the average
stellar mass). Figure 8 shows the inferred burst-phase mass accretion rate as a function of the evolutionary indicator Tbol. This figure aims to reveal the
evolu-tion of episodic accreevolu-tion while the ratios of the post-burst sources (> 10 L ) indicate the burst frequency,
and ˙Macc represents their burst strength. We find mass
accretion rates during the burst phase (Lburst > Lbol
and Lburst > 10 L ) between 4.2 × 10−6 M yr−1 and
4.2 × 10−5 M yr−1 from N2H +
and between 5.1 × 10−6M yr−1and 2.5×10−5M yr−1from HCO+
(Ta-ble 3). The median is ∼ 7.6 × 10−6M yr−1from N2H +
and ∼ 1.6 × 10−5M
yr−1from HCO+with a standard
deviation of ∼ (5.8 − 8.8) × 10−6 M yr−1. This
system-atic shift might come from the adopted model parameters such as the CO binding energy. From an evolutionary point of view, the median does not change from Class 0 to Class I from both N2H+ and HCO+. However, the estimation has a strong bias in selection as the analyses only identify those sources with strong bursts. Besides, the assumption of Mstar = 0.25 M , can be unrealistic
as the stellar mass must increase from the Class 0 to the Class I stage. Thus, Mstarshould be an increasing value
as a function of time, and given the similar accretion lu-minosity during outbursts from Class 0 to Class I, ˙Macc
is subsequently expected to decrease with time. 5.2.3. Mass accumulation of protostars
We have derived the mass accretion rates and the inter-vals between accretion episodes, which allows us to probe the growing process of the central stars. Considering a lifetime of 0.15 − 0.24 Myr (Dunham et al. 2015) and an interval of 2400 yr for the Class 0 stage, accretion bursts would occur 63 − 100 times during this stage. Similarly in Class I, with a lifetime of 0.31 − 0.48 Myr and an in-terval of 8,000 yr, accretion bursts would occur 39 − 60 times; note that our sample includes only one Late Class I protostar (Per-emb-63) with Tbol> 300 K (Evans et al.
2009), and the burst frequency likely keeps decreasing in the Class I stage; thus, the total burst number might be overestimated in Class I. Besides, if the lifetimes of Class 0 and Class I are 30% shorter as suggested by Carney et al. (2016), the total number of bursts would be revised downward by 30%.
Statistically, Enoch et al. (2009) found that ∼ 5% of the embedded protostars have a high luminosity with Lbol> 10 L , or ˙Macc> 10−5 M yr−1. If these sources
are during the burst phase, each burst would last for 2400 × 0.05 = 120 yr at Class 0 stage, which is consistent with the duration of 100 − 200 yr predicted from simu-lations of Vorobyov & Basu (2005). Given the median burst-phase accretion rate (7.6 − 16.2) × 10−6M yr−1,
each burst would thus deliver ∼ 9.1 − 19.4 × 10−4 M
onto the central star. Thus, the protostar would ac-cumulate ∼ 0.06 − 0.19 M at the Class 0 stage and
∼ 0.04 − 0.12 M at the Class I stage. Assuming a mean
stellar mass of 0.5 M (Evans et al. 2009), this result
implies that only ∼ 20 − 60% (0.1 − 0.3 M ) of mass is
accumulated during burst phases. A simple explanation for the discrepancy is that the final stellar mass is in fact smaller; for example, the peak of the initial mass func-tion is ∼ 0.3 M (Muench et al. 2002; Alves et al. 2007).
Alternatively, we propose three possibilities to comple-ment the remaining mass accumulated: (1) an under-estimation of the burst-phase mass accretion rate. The accretion rate is estimated considering a model without a disk component such that the derived accretion lumi-nosity is likely a lower limit. (2) non-negligible mass accumulation during the quiescent phase. As our esti-mated accumulated mass of 0.1 − 0.3 M in Class 0/I
stage arises solely from the accretion burst, it requires . 80% (0.4 M ) of the mass accreted in quiescent phase
with ˙Macc. (5.6 − 8.7) × 10−7 M yr−1 to build a star
with an average mass of 0.5 M . (3) the existence of
super bursts (i.e. FU ori-type). Offner & McKee (2011) estimate that ∼ 25% of mass is accreted during the FU ori events with M˙acc ∼ 10−5− 10−4 M yr−1. This is
consistent with the highest mass accretion rate estimated from the N2H
+
observations as ∼ 4.2 × 10−5 M yr−1.
If such super bursts last for longer, they could deliver significant material onto the central sources.
6. SUMMARY
CO and H2O and destroys N2H +
and HCO+, respec-tively. We compare such derived luminosity to the bolo-metric luminosity (the current luminosity), thus identi-fying the sources that experienced a past accretion burst with Lburst > Lbol, i.e., the post-burst sources. Our
re-sults are summarized as follows:
1. N2H+ and HCO+ peak positions can be used to trace the CO and H2O sublimation regions, respec-tively, and in turn to estimate the luminosity dur-ing the past burst. While N2H+ at large scale is less affected by the system geometry, HCO+in the inner regions is sensitive to the inclination angle but is crucial to trace the past burst over a shorter timescale.
2. We find that 7/17 Class 0 and 1/8 Class I are post-burst sources from HCO+. This decrease of the fraction of post-burst sources may result from the evolution of burst frequency, but we cannot exclude the possibility that the snowline radius is shrunk due to the increase of disk density/size from the Class 0 to the Class I stage. If the disk evolution is not the main factor, then we can draw the following conclusions about the mass accumulation history from the Class 0 to the Class I stage.
3. We derive the intervals between accretion episodes of ∼ 2, 400 yr for Class 0 sources and ∼ 8, 000 yr for Class I sources, suggestive of a decrease in burst frequency during the embedded phase. If the ac-cretion outburst is triggered by disk fragmentation due to gravitational instability, our result suggests that the fragmentation occurs more frequently at an earlier evolutionary stage. Alternatively, the fragment has a higher probability to fall onto the central star at such a stage.
4. We estimate the mass accretion rates at the burst-phase to be (7.6 − 16.2) × 10−6 M yr−1. From an
evolutionary point of view, the burst magnitude is likely unchanged from Class 0 to Class I.
5. Based on the estimate of mass accretion rate and interval between episodes, we derive an accumu-lated mass of 0.06 − 0.19 M at the Class 0 stage
and ∼ 0.04 − 0.12 M at the Class I stage, in
total 0.1 − 0.3 M during burst phases. This
value is smaller than the typical stellar mass of 0.3 − 0.5 M . More material needs to be accreted
to build the star during the quiescent phase or per-haps the star can accumulate mass via a few super accretion bursts.
We are thankful for the referee for many insightful comments for the discussion that helped to improve this paper significantly. The authors thank Merel van ’t Hoff, Jeong-Eun Lee, and Marc Audard for providing valuable discussions. This paper makes use of the fol-lowing ALMA data: ADS/JAO.ALMA#2017.1.01693.S. ALMA is a partnership of ESO (representing its mem-ber states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Re-public of Chile. The Joint ALMA Observatory is oper-ated by ESO, AUI/NRAO and NAOJ. T.H.H. and N. H. acknowledges the support by Ministry of Science and Technology of Taiwan (MoST) 107-2119-M-001-041 and 108-2112-M-001-048. N.H. acknowledges a grant from MoST 108-2112-M-001-017. C.W. acknowledges finan-cial support from the University of Leeds and from the Science Facilities and Technology Council (STFC), under grant number ST/R000549/1. J.K.J. acknowledges sup-port from the European Research Council (ERC) under the European Union’s Horizon 2020 research and inno-vation programme (grant agreement No 646908). S.P.L. acknowledges the support from the Ministry of Science and Technology (MOST) of Taiwan with grant MOST 106-2119-M-007-021-MY3
Software: CASA (McMullin et al. 2007), RADMC-3D (Dullemond et al. 2012), Astropy (Astropy Collaboration et al. 2013), APLpy (Robitaille, & Bressert 2012)
REFERENCES Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17
´
Abrah´am, P., K´osp´al, ´A., Csizmadia, S., et al. 2004, A&A, 419, L39
Acosta-Pulido, J. A., Kun, M., ´Abrah´am, P., et al. 2007, AJ, 133, 2020
Anderl, S., Maret, S., Cabrit, S., et al. 2016, A&A, 591, A3 Andrews, S. M., Rothberg, B., & Simon, T. 2004, ApJ, 610, L45 Aspin, C., Reipurth, B., Beck, T. L., et al. 2009, ApJ, 692, L67 Aso, Y., Ohashi, N., Aikawa, Y., et al. 2017 ApJ, 850, 2 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al.
2013, Astronomy and Astrophysics, 558, A33
Armitage, P. J., Livio, M.,& Pringle J. E. 2001, MNRAS, 324, 705 Audard, M., ´Abrah´am, P., Dunham, M. M., et al. 2014, arXiv:
1401.3368
Banzatti, A., Pinilla, P., Ricci, L., et al. 2015, ApJ, 815, L15 Bell, K. R.,& Lin, D. N. C. 1994, ApJ, 427, 987
Belloche, A., & Andr´e, P. 2004, A&A, 418, L35
Bianchi, E., Codella, C., Ceccarelli, C., et al. 2019, MNRAS, 483, 1850
Bisschop, S. E., Fraser, H. J., ¨Oberg, K. I., et al. 2006, A&A, 449, 1297
Bjerkeli, P., Jrgensen, J. K., Bergin, E. A., et al. 2016, A&A, 595, 39
Boley, A. C.,& Durisen, R. H. 2008, ApJ, 685, 1193
Caratti o Garatti, A., Garcia Lopez, R., Scholz, A., et al. 2011, A&A, 526, L1
Caratti o Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2016, NatPh, 13, 276
Carney, M. T., Yıldız, U. A., Mottram, J. C., et al. 2016, A&A, 586, A44
Chiang, H.-F., Looney, L. W., & Tobin, J. J. 2012, ApJ, 756, 168 Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Nature, 535, 258 Clarke, C. J., & Syer, D., 1996, MNRAS, 278, L23
Codella, C., Maury, A. J., Gueth, F., et al. 2014, A&A, 563, L3 Collings, M. P., Dever, J. W., Fraser, H. J., et al. 2003, ApJ, 583,
1058
Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133
Contreras Pe˜na, C., Naylor, T., & Morrell, S. 2019, MNRAS, 486, 4590
Covey, K. R., Hillenbrand, L. A., Miller, A. A., et al. 2011, AJ, 141, 40
Davis, C. J., Scholz, P., Lucas, P., Smith, M. D., & Adamson, A. 2008, MNRAS, 387, 954
De Simone, M., Codella, C., Testi, L., et al. 2017, A&A, 599, A121
Dunham, M. M., Evans, N. J., Bourke, T. L., et al. 2010, ApJ, 721, 995
Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, ascl:1202.015
