# Massive Compact Disks around FU Orionis-type Young Eruptive Stars Revealed by ALMA

Á. Kóspál, F. Cruz-Sáenz de Miera, J. White, P. Ábrahám, L. Chen, T. Csengeri, R. Dong, M. Dunham, O. Fehér, J. Green, et al.

► To cite this version:

Á. Kóspál, F. Cruz-Sáenz de Miera, J. White, P. Ábrahám, L. Chen, et al. Massive Compact Disks around FU Orionis-type Young Eruptive Stars Revealed by ALMA. The Astrophysical Journal Supplement, 2021, 256 (2), pp.30. 10.3847/1538-4365/ac0f09. insu-03439571

HAL Id: insu-03439571

https://insu.hal.science/insu-03439571v1

Submitted on 2 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Massive Compact Disks around FU Orionis-type Young Eruptive Stars Revealed by ALMA

Á. Kóspál<sup>1,2,3 Ⓓ</sup>, F. Cruz-Sáenz de Miera<sup>1 Ⓓ</sup>, J. A. White<sup>1,4,18 Ⓓ</sup>, P. Ábrahám<sup>1,3 Ⓓ</sup>, L. Chen<sup>1 Ⓓ</sup>, T. Csengeri<sup>5 Ⓓ</sup>, R. Dong<sup>6 Ⓓ</sup>, M. M. Dunham<sup>7 Ⓓ</sup>, O. Fehér<sup>1,8 Ⓓ</sup>, J. D. Green<sup>9,10 Ⓓ</sup>, J. Hashimoto<sup>11 Ⓓ</sup>, Th. Henning<sup>2 Ⓓ</sup>, M. Hogerheijde<sup>12,13 Ⓓ</sup>, T. Kudo<sup>14 Ⓓ</sup>, H. B. Liu<sup>15 Ⓓ</sup>, M. Takami<sup>15 Ⓓ</sup>, and E. I. Vorobyov<sup>16,17 Ⓓ</sup>

<sup>1</sup> Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary; kospal@konkoly.hu

<sup>2</sup> Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

<sup>3</sup> ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

<sup>4</sup> National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

<sup>5</sup> Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>6</sup> Department of Physics & Astronomy, University of Victoria, Victoria, BC, V8P 1A1, Canada

<sup>7</sup> Department of Physics, State University of New York at Fredonia, Fredonia, NY 14063, USA

<sup>8</sup> IRAM, 300 Rue de la piscine, F-38406 Saint-Martin-d'Hères, France

<sup>9</sup> The University of Texas at Austin, Department of Astronomy, 2515 Speedway, Stop C1400, Austin, TX 78712, USA

<sup>10</sup> Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 20138, USA

<sup>11</sup> Astrobiology Center of NINS, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan

<sup>12</sup> Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands

<sup>13</sup> Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

<sup>14</sup> Subaru Telescope, National Astronomical Observatory of Japan, National Institutes of Natural Sciences (NINS), 650 North Aohoku Place, Hilo, HI 96720, USA
<sup>15</sup> 11F of AS/NTU Astronomy-Mathematics Building, No.1, Section 4, Roosevelt Road, Taipei 10617, Taiwan, R.O.C.

<sup>16</sup> Department of Astrophysics, The University of Vienna, Vienna, A-1180, Austria

<sup>17</sup> Research Institute of Physics, Southern Federal University, Stachki 194, Rostov-on-Don, 344090, Russia

Received 2021 April 10; revised 2021 June 25; accepted 2021 June 26; published 2021 September 23

## Abstract

The FU Orionis-type objects (FUors) are low-mass pre-main-sequence stars undergoing a temporary but significant increase of mass accretion rate from the circumstellar disk onto the protostar. It is not yet clear what triggers the accretion bursts and whether the disks of FUors are in any way different from the disks of nonbursting young stellar objects. Motivated by this, we conducted a 1.3 mm continuum survey of 10 FUors and FUor-like objects with ALMA, using both the 7 m array and the 12 m array in two different configurations to recover emission at the widest possible range of spatial scales. We detected all targeted sources and several nearby objects as well. To constrain the disk structure, we fit the data with models of increasing complexity from 2D Gaussian to radiative transfer, enabling comparison with other samples modeled in a similar way. The radiative transfer modeling gives disk masses that are significantly larger than what is obtained from the measured millimeter fluxes assuming optically thin emission, suggesting that the FUor disks are optically thick at this wavelength. In comparison with samples of regular class II and class I objects, the disks of FUors are typically a factor of 2.9–4.4 more massive and a factor of 1.5–4.7 smaller in size. A significant fraction of them (65%–70%) may be gravitationally unstable.

Unified Astronomy Thesaurus concepts: Young stellar objects (1834); FU Orionis stars (553); Circumstellar disks (235); Submillimeter astronomy (1647); Interferometry (808)

Supporting material: data behind figure

## 1. Introduction

The FU Orionis-type objects (FUors) are low-mass pre-main-sequence stars undergoing large outbursts in visible light, attributed to temporarily enhanced accretion (Hartmann & Kenyon 1996). Episodes of high accretion are believed to occur at all early stages of star formation after the prestellar core phase but become observable at optical or infrared wavelengths only as the circumstellar envelope thins (Audard et al. 2014, and references therein). During outburst, the accretion rate from the circumstellar disk to the star may reach up to $10^{-4} M_{\odot} \text{yr}^{-1}$, several orders of magnitude higher than in quiescence or in normal T Tauri stars. Repetitive FUor outbursts, each lasting ~100 yr and dumping as much as $10^{-2} M_{\odot}$ onto the star, significantly contribute to the buildup of the final stellar mass

and offer a possible solution for the long-standing “protostellar luminosity problem” (e.g., Dunham et al. 2014a). The physical mechanism of FUor outbursts is not yet known, but several processes have been proposed (Audard et al. 2014).

The FUor outbursts may eventually clear up the protostar’s environment at the end of the embedded (class I) phase. By this time, envelopes are depleted to a level comparable to the mass accreted in a single outburst, and UV/X-ray radiation and outflows accompanying the eruptions may disperse them altogether. Thus, individual FUor eruptions might modify the structure of the circumstellar environment considerably, both the envelope and the disk itself (Vorobyov et al. 2020). Following a sequence of outbursts (the “FUor phase”), the envelope disappears, and the young star enters the disk-only class II phase. The FUor eruptions may also change the mineralogy and chemistry of the disk (Molyarova et al. 2018). Therefore, FUor outbursts mark a short but important

<footnote n="18"><sup>18</sup> Jansky Fellow of the National Radio Astronomy Observatory.</footnote>
transitional period in the evolution of the circumstellar material (Takami et al. 2018, 2019). Numerical simulations support the notion that most FUors occur in the class I or early class II phases (Vorobyov & Basu 2015). Exploring this period via characterizing FUors, in particular their disks where the outburst happens, is the main topic of the present paper.

Millimeter continuum emission has long been used to study the dusty environment of FUors (Henning et al. 1998; Sandell & Weintraub 2001). Interferometric observations of the millimeter dust continuum are the standard way to estimate basic disk properties like the mass and size or characteristic radius. Using the Atacama Large Millimeter/submillimeter Array (ALMA), Cieza et al. (2018) observed three FUors, four EXors (objects showing smaller, shorter, repetitive accretion outbursts), and one intermediate-type object. They found that FUor disks are significantly more massive than EXor disks, and they are more compact for a given disk mass than what is typical for normal T Tauri disks. Liu et al. (2018) reached a similar conclusion in a Submillimeter Array (SMA) study of 29 FUors, EXors, and FUor-like objects. Overall, the observed targets showed a systematically higher millimeter luminosity distribution than those of normal T Tauri stars in nearby star-forming regions, and within the sample, all EXors had smaller 1.33 mm fluxes (when scaled to the same distance) than any of the FUors. This suggests that relatively massive circumstellar disks may provide optimal conditions for accretion outbursts.

Not all FUors are isolated objects; some of them are parts of small clusters of young stars (like AR 6A/6B, HBC 722, or V960 Mon; Aspin & Reipurth 2003; Kóspál et al. 2015, 2016), or they may even be close binaries (like FU Ori or L1551 IRS 5; Rodríguez et al. 1998; Wang et al. 2004). Multiplicity affects disk formation and evolution (Tobin et al. 2018). In the SMA sample of eruptive young stars, Liu et al. (2018) found that binaries or triple systems are systematically fainter than the rest of the sources.

In this paper, we present new ALMA 1.33 mm continuum observations of 10 FUors and FUor-like objects, including some multiple systems. In Section 2, we explain the target selection and describe the observations and data reduction. In Section 3, we show our ALMA images and derive basic disk properties by Gaussian fitting in the image space and analytical disk model fitting in the visibility space. We note that due to the finite sensitivity of our observations, and given that ALMA filters out the emission on large spatial scales, our data mostly trace structures on the disks' spatial scale, even though some of our targets are embedded in envelopes. In Section 4, we detail our radiative transfer (RT) disk modeling. We explain our new results on individual objects in the context of the literature in Section 5, discuss the whole sample and draw our conclusions in Section 6, and summarize our main findings in Section 7.

## 2. Observations and Data Reduction

### 2.1. Target Selection

For our project, we selected those FUors from the eruptive star list of Audard et al. (2014) that were not included in any previous ALMA project but still observable from the ALMA site. These were AR 6A/6B, Bran 76, HBC 687, Haro 5a IRS, OO Ser, V346 Nor, V900 Mon, and L1551 IRS 5. We added two extra FUors that were discovered after the publication of Audard et al. (2014), V899 Mon and V960 Mon. Basic data (distances and coordinates) for our targets are shown in

Table 1. Recently, Connelley & Reipurth (2018) completed a near-infrared spectroscopic survey of FUors and classified them based on various criteria as bona fide FUors (whose spectra display the necessary near-IR characteristics and an eruption was observed), FUor-like objects (whose spectra display the necessary near-IR characteristics but no eruption was observed), and peculiar objects (which have or had some spectral similarities to FUors but are different enough that a FUor classification is not warranted). Based on their classification, from our sample, V900 Mon and V960 Mon are bona fide FUors; L1551 IRS 5, Haro 5a IRS, Bran 76, and HBC 687 are FUor-like objects; and V346 Nor and AR 6A/6B are peculiar objects. V899 Mon and OO Ser were not studied in Connelley & Reipurth (2018).

### 2.2. ALMA Observations and Data Processing

We observed our targets with ALMA in Band 6 between 2016 October 9 and 2018 April 5. A detailed log of observations can be found in Table 1. For each object, we took observations with the 7 m array and two different configurations of the 12 m array (TM2 being the more compact C40–3 configuration and TM1 being the more extended C40–6 configuration). We used a standard frequency setup that covered the $J=2-1$ line of $^{12}\text{CO}$, $^{13}\text{CO}$, and $\text{C}^{18}\text{O}$ in 0.58 GHz wide spectral windows and defined two 1.875 GHz wide windows centered at 234.0 and 216.5 GHz to cover the emission from the thermal dust continuum. Here we focus on the continuum results, while the molecular line data will be analyzed in subsequent papers. For the purpose of atmospheric, bandpass, flux, phase, and pointing calibration, the ALMA observatory mostly selected different bright quasars. In a few cases, the atmospheric and flux calibrator was Neptune, Titan, or Uranus.

We reduced the delivered data with the ALMA pipeline using the Common Astronomy Software Applications (CASA) package version 5.1.1 (McMullin et al. 2007). The observational executions with a high signal-to-noise ratio (S/N) were self-calibrated by using 30 s intervals for the creation of the phase-calibration tables and the total integration time for the amplitude calibration. We checked the continuum spectral windows for line emission and selected the uncontaminated channels to extract the continuum using the uvcontsub task in CASA. After verifying that the common baselines between the three different observing configurations provided consistent visibilities, we concatenated all observations and produced a single measurement set for each target, which recovers emission in a wide range of spatial scales up to 5"/9, the maximum recoverable scale for the C40–3 configuration in Band 6. We extracted the continuum from the combined data sets and created uniformly weighted images using the CLEAN algorithm implemented in the tclean task for CASA (version 5.4.0). The reference frequency for all maps was fixed at 224.55 GHz. We ran tclean noninteractively and set the flux density threshold to three times the noise level of the dirty image. We applied primary beam correction to the final clean images. The primary beam sizes at this frequency are 28" for the 12 m antennas and 48" for the 7 m antennas. The rms of our non-primary-beam-corrected cleaned maps and the sizes and position angles (PAs) of our resulting beams are listed in Table 2.
<table><caption>**Table 1**
Log of Our ALMA Continuum Observations</caption><thead><tr><th>Name</th><th>Distance<sup>a</sup><br>(pc)</th><th>R.A.<sub>12000</sub><br>(h:m:s)</th><th>Decl.<sub>12000</sub><br>(deg:arcmin:arcsec)</th><th>Date</th><th>Config.<sup>b</sup></th><th>Ang. Res.<br>(arcsec)</th><th>T<sub>exp</sub><br>(minutes)</th><th>PWV<br>(mm)</th><th>Baselines<br>(m)</th></tr></thead><tbody><tr><td rowspan="3">L1551 IRS 5</td><td rowspan="3">140 ± 10</td><td rowspan="3">04:31:34.08</td><td rowspan="3">+18:08:04.9</td><td>2016 Nov 14</td><td>7 m</td><td>5.2</td><td>4.6</td><td>...</td><td>8–44</td></tr><tr><td>2018 Apr 5</td><td>TM2</td><td>0.68</td><td>3.6</td><td>0.8</td><td>15–500</td></tr><tr><td>2017 Jul 24</td><td>TM1</td><td>0.12</td><td>11.3</td><td>0.2</td><td>30–3637</td></tr><tr><td rowspan="5">Haro 5a IRS</td><td rowspan="5">391 ± 24</td><td rowspan="5">05:35:26.71</td><td rowspan="5">−05:03:54.9</td><td>2016 Oct 10</td><td>7 m</td><td>5.1</td><td>4.1</td><td>...</td><td>8–48</td></tr><tr><td>2017 Apr 6</td><td>7 m</td><td>5.1</td><td>4.1</td><td>...</td><td>8–48</td></tr><tr><td>2017 Aug 23</td><td>TM2</td><td>0.45</td><td>3.1</td><td>2.1</td><td>15–783</td></tr><tr><td>2017 Jul 7</td><td>TM1</td><td>0.19</td><td>9.3</td><td>0.8</td><td>16–2647</td></tr><tr><td>2017 Jul 26</td><td>TM1</td><td>0.11</td><td>9.3</td><td>0.4</td><td>16–3696</td></tr><tr><td rowspan="5">V899 Mon</td><td rowspan="5">769<sup>+36</sup><sub>−33</sub></td><td rowspan="5">06:09:19.28</td><td rowspan="5">−06:41:55.4</td><td>2016 Oct 9</td><td>7 m</td><td>5.1</td><td>4.1</td><td>...</td><td>8–48</td></tr><tr><td>2017 Apr 29</td><td>TM2</td><td>0.67</td><td>3.1</td><td>2.0</td><td>15–459</td></tr><tr><td>2018 Mar 22</td><td>TM2</td><td>0.46</td><td>3.1</td><td>2.2</td><td>15–783</td></tr><tr><td>2017 Jul 5</td><td>TM1</td><td>0.19</td><td>9.3</td><td>0.6</td><td>21–2647</td></tr><tr><td>2017 Jul 27</td><td>TM1</td><td>0.11</td><td>9.3</td><td>0.5</td><td>16–3696</td></tr><tr><td rowspan="3">AR 6A/6B</td><td rowspan="3">719 ± 16</td><td rowspan="3">06:40:59.31</td><td rowspan="3">+09:35:52.0</td><td>2016 Oct 9</td><td>7 m</td><td>5.1</td><td>4.1</td><td>...</td><td>8–48</td></tr><tr><td>2017 Apr 20</td><td>TM2</td><td>0.70</td><td>3.1</td><td>2.3</td><td>16–459</td></tr><tr><td>2017 Jul 22</td><td>TM1</td><td>0.13</td><td>10.3</td><td>1.0</td><td>30–3636</td></tr><tr><td rowspan="3">V900 Mon</td><td rowspan="3">1100 ± 100</td><td rowspan="3">06:57:22.22</td><td rowspan="3">−08:23:17.6</td><td>2016 Oct 11</td><td>7 m</td><td>4.7</td><td>8.1</td><td>...</td><td>8–47</td></tr><tr><td>2017 Apr 20</td><td>TM2</td><td>0.70</td><td>6.1</td><td>2.2</td><td>16–459</td></tr><tr><td>2017 Jul 27</td><td>TM1</td><td>0.11</td><td>18.6</td><td>0.5</td><td>16–3696</td></tr><tr><td rowspan="3">V960 Mon</td><td rowspan="3">1574<sup>+173</sup><sub>−143</sub></td><td rowspan="3">06:59:31.58</td><td rowspan="3">−04:05:28.1</td><td>2016 Oct 11</td><td>7 m</td><td>4.7</td><td>8.1</td><td>...</td><td>8–47</td></tr><tr><td>2017 Apr 20</td><td>TM2</td><td>0.70</td><td>6.1</td><td>2.2</td><td>16–459</td></tr><tr><td>2017 Jul 27</td><td>TM1</td><td>0.11</td><td>18.6</td><td>0.5</td><td>16–3696</td></tr><tr><td rowspan="3">Bran 76</td><td rowspan="3">1060<sup>+29</sup><sub>−28</sub></td><td rowspan="3">07:50:35.60</td><td rowspan="3">−33:06:23.9</td><td>2016 Oct 21</td><td>7 m</td><td>4.7</td><td>8.3</td><td>...</td><td>8–47</td></tr><tr><td>2017 Apr 20</td><td>TM2</td><td>0.70</td><td>6.2</td><td>2.3</td><td>16–459</td></tr><tr><td>2017 Jul 30</td><td>TM1</td><td>0.10</td><td>20.7</td><td>0.8</td><td>16–3637</td></tr><tr><td rowspan="3">V346 Nor</td><td rowspan="3">700<sup>+600</sup><sub>−350</sub></td><td rowspan="3">16:32:32.19</td><td rowspan="3">−44:55:30.7</td><td>2017 May 15</td><td>7 m</td><td>5.2</td><td>4.1</td><td>...</td><td>8–47</td></tr><tr><td>2018 Mar 16</td><td>TM2</td><td>0.46</td><td>3.1</td><td>2.0</td><td>15–783</td></tr><tr><td>2017 Aug 25</td><td>TM1</td><td>0.10</td><td>9.8</td><td>0.4</td><td>21–3696</td></tr><tr><td rowspan="4">OO Ser</td><td rowspan="4">438 ± 11</td><td rowspan="4">18:29:49.13</td><td rowspan="4">+01:16:20.7</td><td>2017 Apr 30</td><td>7 m</td><td>5.2</td><td>4.1</td><td>...</td><td>8–48</td></tr><tr><td>2017 May 14</td><td>7 m</td><td>5.2</td><td>4.1</td><td>...</td><td>8–47</td></tr><tr><td>2018 Mar 27</td><td>TM2</td><td>0.46</td><td>3.1</td><td>1.9</td><td>15–783</td></tr><tr><td>2017 Jul 22</td><td>TM1</td><td>0.12</td><td>9.8</td><td>0.5</td><td>16–3696</td></tr><tr><td rowspan="3">HBC 687</td><td rowspan="3">400</td><td rowspan="3">19:29:00.86</td><td rowspan="3">+09:38:42.9</td><td>2016 Nov 18</td><td>7 m</td><td>5.2</td><td>9.3</td><td>...</td><td>8–44</td></tr><tr><td>2017 Apr 30</td><td>TM2</td><td>0.67</td><td>7.3</td><td>1.6</td><td>13–459</td></tr><tr><td>2017 Jul 22</td><td>TM1</td><td>0.12</td><td>22.8</td><td>0.5</td><td>16–3696</td></tr></tbody></table>

**Notes.**

**3. Results**

Figure 1 shows our 1.33 mm images corresponding to the thermal dust continuum emission. All of our targets were detected with a high peak S/N of at least 23. Some of our maps show multiple sources. In the map of the AR 6A/6B system, we detect three sources that can be identified as 6A, 6B, and 6C (see Aspin & Reipurth 2003). In addition, we also detected a bright source to the southwest of AR 6B, which we have labeled AR 6D and is outside of the area plotted in Figure 1. In the map of Haro 5a IRS, there is a faint source near our target that had not been previously reported in the literature. In the following, we will call this source Haro 5a IRS B. Around V960 Mon, there are three nearby sources, two of which can be

identified as V960 Mon N and V960 Mon SE (see Kóspál et al. 2015). In the following, we will call the third one V960 Mon E. We detected V370 Ser in our OO Ser map at a separation of ~11"5 to the north of the FUor (not visible in Figure 1, as it falls outside of the plotted region). We detected the continuum source to the northwest of V346 Nor that had previously been reported by Kóspál et al. (2017a). At a separation of 15"2, this source also falls outside of the region plotted in Figure 1.

The brightness distribution of our targets typically has a centrally peaked, bright, rather compact central part, in some cases surrounded by extended faint emission. To derive basic parameters such as coordinates, sizes, inclinations, PAs, and fluxes, we used CASA's IMFIT tool to find the best-fit 2D Gaussian models for each FUor disk. The radii were estimated

<footnote><sup>a</sup> References for the distances are as follows. L1551 IRS 5: Kenyon et al. (1994), Haro 5a IRS: Großschedl et al. (2018), V899 Mon: Bailer-Jones et al. (2018), AR 6A/6B: Maíz Apellániz (2019), V900 Mon: Reipurth et al. (2012), V960 Mon: Bailer-Jones et al. (2018), Bran 76: Bailer-Jones et al. (2018), V346 Nor: Graham & Frogel (1985), OO Ser: Herczeg et al. (2019), HBC 687: Kóspál et al. (2008).</footnote>

<footnote><sup>b</sup> TM1 means the C40–6 configuration of the 12 m array, while TM2 means the C40–3 configuration.</footnote>
<table><caption>**Table 2**
CLEAN Beam and Noise</caption><thead><tr><th>Target</th><th>Beam Size<br>(arcsec)</th><th>Beam PA<br>(deg)</th><th>rms<br>(mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td>L1551 IRS 5</td><td>0.14 × 0.10</td><td>85.7</td><td>0.374</td></tr><tr><td>Haro 5a IRS</td><td>0.18 × 0.18</td><td>-85.5</td><td>0.138</td></tr><tr><td>V899 Mon</td><td>0.15 × 0.15</td><td>84.8</td><td>0.023</td></tr><tr><td>AR 6A/6B</td><td>0.17 × 0.17</td><td>-87.9</td><td>0.046</td></tr><tr><td>V900 Mon</td><td>0.13 × 0.13</td><td>85.6</td><td>0.030</td></tr><tr><td>V960 Mon</td><td>0.14 × 0.14</td><td>85.5</td><td>0.028</td></tr><tr><td>Bran 76</td><td>0.13 × 0.13</td><td>-77.5</td><td>0.027</td></tr><tr><td>V346 Nor</td><td>0.13 × 0.13</td><td>82.1</td><td>0.031</td></tr><tr><td>OO Ser</td><td>0.15 × 0.15</td><td>85.7</td><td>0.064</td></tr><tr><td>HBC 687</td><td>0.14 × 0.14</td><td>87.8</td><td>0.042</td></tr></tbody></table>

from the deconvolved best-fit 2D Gaussian; we took the FWHM of the major axis and divided it by $2\sqrt{2\ln 2}$ to obtain $\sigma$, which we use as a proxy for the radius of the disk based on their primary beam-corrected images. We also estimated their masses by assuming optically thin emission, a dust opacity of $\kappa_{\nu}$, and a dust temperature of $T_{\text{dust}}$ and using the following equation:

$$M_{\text{dust}} = \frac{F_{\nu} d^2}{\kappa_{\nu} B_{\nu} (T_{\text{dust}})} \quad (1)$$

The dust opacity was estimated following Beckwith et al. (1990), where $\kappa_{\nu} = 0.1(\nu/10^{12}\text{Hz})^{\beta}$, assuming an opacity power-law index of $\beta = 1$, which resulted in $\kappa_{\nu} = 2.254 \text{ cm}^2 \text{ g}^{-1}$ at 1.33 mm.

The dust temperature was calculated following Equation (2) of Tobin et al. (2020):

$$T_{\text{bol}}^{\text{dust}} = 43 \text{ K} \left( \frac{L_{\text{bol}}}{1 L_{\odot}} \right)^{0.25} \quad (2)$$

In the case of the known FUors, we used the bolometric luminosities measured during their outburst phase (Audard et al. 2014, and references therein). For the nonoutbursting sources in our fields of view, we used the median luminosity estimated from 230 protostars (Dunham et al. 2013), $L_{\text{bol}} = 1.3 L_{\odot}$, which resulted in ~47 K, and we rounded to $T_{\text{bol}} = 50$ K. We note that other works have used lower temperatures, e.g., 20 K (Ansdell et al. 2016, 2017), which would result in higher mass estimates.

We then converted the obtained dust masses to total masses assuming the usual gas-to-dust mass ratio of 100. We emphasize that the millimeter emission is likely (partially) optically thick, so these values should be considered as lower limits for the true mass (see Dunham et al. 2014b; Liu 2019; Zhu et al. 2019). The obtained parameters are shown in Table 3. With the 2D Gaussian fitting procedure, V899 Mon, AR 6A, and AR 6C appear to be unresolved; therefore, we could not determine radius, inclination, and PA for them. However, we provide upper limits for their radii.

We measured the in-band spectral index $\alpha$ of our targets by constructing maps from the two continuum spectral windows separately, at $\nu_1 = 216.9$ and $\nu_2 = 232.2$ GHz, and measuring the targets' fluxes ($F$) at these two frequencies. Then we

calculated $\alpha$ using

$$\alpha = \frac{\log F_{\nu_1} - \log F_{\nu_2}}{\log \nu_1 - \log \nu_2} \quad (3)$$

and its uncertainties ($\sigma$) with

$$\sigma_{\alpha}^{2} = \left( \frac{1}{\ln \nu_{1} - \ln \nu_{2}} \right)^{2} \left( \frac{\sigma_{\nu_{1}}^{2}}{F_{\nu_{1}}^{2}} + \frac{\sigma_{\nu_{2}}^{2}}{F_{\nu_{2}}^{2}} \right) \quad (4)$$

If the dust emission comes from large grains that emit as blackbodies or if the emission is optically thick, $\alpha$ is expected to be $\alpha = 2$ (Pavlyuchenkov et al. 2019). If the emission is optically thin and the grain size distribution is typical for the interstellar matter, then the millimeter spectral slope is expected to be steeper, with $\alpha \approx 3.7$ (Draine 2006). The measured indices are presented in Table 3. We obtained spectral indices between 0.27 and 4.95; however, these values are highly uncertain due to the small difference between frequencies $\nu_1$ and $\nu_2$. We did not make further analysis of these values.

### 3.1. Fitting in the Visibility Space

Fitting disk properties in the image space is sensitive to the deconvolution procedure used. Therefore, to obtain results independently of the deconvolution, we fitted the disk emission in the visibility space using a procedure similar to that described in Manara et al. (2019). Based on the continuum images, we assumed that the disks were axisymmetrical and modeled each of them using a radial profile

$$I(r) = I_0 r^{-\gamma_1} \exp\left(-\frac{r}{R_{\text{vis}}}\right)^{\gamma_2} \quad (5)$$

where $I_0$ is the peak emission of the disk in janskys per steradian, $R_{\text{vis}}$ is the characteristic radius of the disk in arcseconds, $\gamma_1$ is the exponent of the power law and describes the flatness of the inner parts of the brightness radial profile, and $\gamma_2$ is the exponent of the exponential cutoff that dictates the cutoff of the outer parts of the radial profile. We used the Galario<fnref n="19" /> code (Tazzari et al. 2018) to apply the appropriate inclination and position angle ($i$ and PA, respectively) and the offsets with respect to the phase center ($\Delta R.A.$ and $\Delta \text{decl.}$) and transform the radial profile model into visibilities ($V_m$). Therefore, each disk has been fitted with eight free parameters: $I_0$, $R_{\text{vis}}$, $\gamma_1$, $\gamma_2$, $i$, PA, $\Delta R.A.$, and $\Delta \text{decl.}$ We explored the parameter space with a Goodman & Weare affine invariant Markov Chain Monte Carlo (MCMC) using the Python ensemble sampler emcee<fnref n="20" /> (Foreman-Mackey et al. 2013). Our likelihood function was also calculated with Galario following

$$\chi^2 = \sum_{j=1}^{N} \left[ \begin{aligned} & (\text{Re}(V_{o,j}) - \text{Re}(V_{m,j}))^2 \\ & + (\text{Im}(V_{o,j}) - \text{Im}(V_{m,j}))^2 \end{aligned} \right] w_j \quad (6)$$

where $N$ is the number of visibility points and $w_j$ is the weight corresponding to each point.

For the cases where our field of view includes more than one source (i.e., Haro 5a IRS, AR 6, V960 Mon, V346 Nor, and

<footnote n="19">¹⁹ https://mtazzari.github.io/galario/index.html</footnote>

<footnote n="20">²⁰ https://emcee.readthedocs.io/en/stable/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The 1.33 mm cleaned non-primary-beam-corrected ALMA images of our sample. Each panel is 4'' × 4''. The synthesized beam is indicated in the lower left corner of each panel. The target name and contour levels are indicated in the upper left corner of the panels. The beam sizes and rms noise ($\sigma$) values are listed in Table 2.</figcaption></figure>

(The data used to create this figure are available.)
<table><caption>**Table 3**
Basic Parameters of the Detected Sources from Gaussian Fits</caption><thead><tr><th>Name</th><th>Radius<br>(au)</th><th><i>i</i><br>(deg)</th><th>PA<br>(deg)</th><th>Int. Flux<br>(mJy)</th><th>Mass (50 K)<br>(M<sub>☉</sub>)</th><th><i>L</i><sub>bol</sub><br>(L<sub>☉</sub>)</th><th><i>T</i><sub>bol</sub><sup>dust</sup><br>(K)</th><th>Mass (<i>T</i><sub>bol</sub><sup>dust</sup>)<br>(M<sub>☉</sub>)</th><th>α</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">FUors</td></tr><tr><td>L1551 IRS 5 N</td><td>10.6 ± 0.7</td><td>44 ± 7</td><td>158 ± 12</td><td>248.7 ± 9.4</td><td>1.48 × 10<sup>-2</sup></td><td>30 (a)</td><td>101</td><td>6.96 × 10<sup>-3</sup></td><td>2.08 ± 0.93</td></tr><tr><td>Haro 5a IRS</td><td>39 ± 2</td><td>55 ± 1</td><td>169 ± 2</td><td>166.4 ± 1.7</td><td>7.73 × 10<sup>-2</sup></td><td>50 (b)</td><td>114</td><td>3.18 × 10<sup>-2</sup></td><td>1.41 ± 0.36</td></tr><tr><td>V899 Mon</td><td>&lt;37</td><td>...</td><td>...</td><td>...</td><td>1.6 ± 0.1</td><td>150 (c)</td><td>151</td><td>8.87 × 10<sup>-4</sup></td><td>3.78 ± 1.15</td></tr><tr><td>AR 6A</td><td>&lt;38</td><td>...</td><td>...</td><td>...</td><td>1.5 ± 0.2</td><td>450 (d)</td><td>198</td><td>5.48 × 10<sup>-4</sup></td><td>4.95 ± 2.89</td></tr><tr><td>AR 6B</td><td>31 ± 4</td><td>54 ± 27</td><td>12 ± 27</td><td>4.9 ± 0.4</td><td>7.70 × 10<sup>-3</sup></td><td>450 (d)</td><td>198</td><td>1.79 × 10<sup>-3</sup></td><td>2.23 ± 2.01</td></tr><tr><td>V900 Mon</td><td>20 ± 2</td><td>28 ± 20</td><td>169 ± 73</td><td>8.4 ± 0.1</td><td>3.09 × 10<sup>-2</sup></td><td>106 (e)</td><td>138</td><td>1.04 × 10<sup>-2</sup></td><td>2.62 ± 0.51</td></tr><tr><td>V960 Mon</td><td>37 ± 28</td><td>61 ± 28</td><td>71 ± 20</td><td>0.7 ± 0.1</td><td>5.27 × 10<sup>-3</sup></td><td>...</td><td>...</td><td>...</td><td>3.51 ± 4.11</td></tr><tr><td>Bran 76</td><td>24 ± 3</td><td>45 ± 11</td><td>79 ± 18</td><td>4.8 ± 0.1</td><td>1.64 × 10<sup>-2</sup></td><td>287 (f)</td><td>177</td><td>4.28 × 10<sup>-3</sup></td><td>2.13 ± 0.75</td></tr><tr><td>V346 Nor</td><td>42 ± 28</td><td>38 ± 10</td><td>105 ± 23</td><td>11.7 ± 0.7</td><td>1.74 × 10<sup>-2</sup></td><td>160 (g)</td><td>153</td><td>5.29 × 10<sup>-3</sup></td><td>3.44 ± 1.92</td></tr><tr><td>OO Ser</td><td>26 ± 2</td><td>50 ± 6</td><td>60 ± 10</td><td>12.8 ± 0.5</td><td>7.46 × 10<sup>-3</sup></td><td>36 (h)</td><td>102</td><td>3.48 × 10<sup>-3</sup></td><td>2.79 ± 1.38</td></tr><tr><td>HBC 687</td><td>20 ± 2</td><td>64 ± 1</td><td>77 ± 1</td><td>30.5 ± 0.4</td><td>1.48 × 10<sup>-2</sup></td><td>10 (i)</td><td>77</td><td>9.33 × 10<sup>-3</sup></td><td>2.40 ± 0.31</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">Neighboring Sources</td></tr><tr><td>L1551 IRS 5 S</td><td>9.5 ± 1.1</td><td>25 ± 23</td><td>160 ± 77</td><td>153.2 ± 9.3</td><td>9.12 × 10<sup>-3</sup></td><td>10 (a)</td><td>77</td><td>5.74 × 10<sup>-3</sup></td><td>3.52 ± 1.6</td></tr><tr><td>Haro 5a IRS B</td><td>34 ± 6</td><td>56 ± 7</td><td>97 ± 10</td><td>13.4 ± 1.3</td><td>6.22 × 10<sup>-3</sup></td><td>...</td><td>...</td><td>...</td><td>0.27 ± 2.26</td></tr><tr><td>AR 6C</td><td>&lt;38</td><td>...</td><td>...</td><td>...</td><td>1.10 × 10<sup>-3</sup></td><td>...</td><td>...</td><td>...</td><td>4.77 ± 3.88</td></tr><tr><td>AR 6D</td><td>50 ± 7</td><td>56 ± 14</td><td>43 ± 14</td><td>4.1 ± 0.3</td><td>6.44 × 10<sup>-3</sup></td><td>...</td><td>...</td><td>...</td><td>2.50 ± 1.60</td></tr><tr><td>V960 Mon N</td><td>63 ± 16</td><td>34 ± 41</td><td>33 ± 57</td><td>1.2 ± 0.1</td><td>9.03 × 10<sup>-3</sup></td><td>...</td><td>...</td><td>...</td><td>0.95 ± 1.95</td></tr><tr><td>V960 Mon E</td><td>126 ± 29</td><td>65 ± 8</td><td>96 ± 11</td><td>1.5 ± 0.2</td><td>1.13 × 10<sup>-2</sup></td><td>...</td><td>...</td><td>...</td><td>0.28 ± 3.03</td></tr><tr><td>V960 Mon SE</td><td>227 ± 25</td><td>73 ± 2</td><td>170 ± 2</td><td>7.6 ± 0.3</td><td>5.72 × 10<sup>-2</sup></td><td>...</td><td>...</td><td>...</td><td>1.74 ± 1.00</td></tr><tr><td>V346 Nor B</td><td>84 ± 57</td><td>67 ± 4</td><td>145 ± 4</td><td>17.8 ± 1.1</td><td>2.65 × 10<sup>-2</sup></td><td>...</td><td>...</td><td>...</td><td>2.69 ± 1.36</td></tr><tr><td>V370 Ser</td><td>57 ± 6</td><td>55 ± 8</td><td>15 ± 9</td><td>5.9 ± 0.6</td><td>3.44 × 10<sup>-3</sup></td><td>...</td><td>...</td><td>...</td><td>1.23 ± 1.96</td></tr></tbody></table>

**Notes.** References for $L_{bol}$: (a) Liseau et al. (2005), (b) Reipurth & Aspin (1997), (c) Ninan et al. (2015), (d) Aspin & Reipurth (2003), (e) Reipurth et al. (2012), (f) Green et al. (2006), (g) Kóspál et al. (2017b), (h) Kóspál et al. (2007), (i) Kóspál et al. (2008). Disk masses were obtained from the integrated fluxes using Equation (1) and should be considered as lower limits. Disk radii are deconvolved sizes.

OO Ser), we used the same model for them, calculated the visibilities individually, and added them together before calculating $\chi^2$. The initial values of the exploration of the parameter space were based on the results of our 2D Gaussian fitting (Table 3). We utilized flat priors for all variables. Depending on the number of fitted disks, we used a different number of walkers and steps. In cases where only one disk is in the field of view, we used 150 walkers and at least 10,000 steps. In cases where more than one disk was detected, we used 250 walkers and between 25,000 and 35,000 steps. For all cases, we discarded a few thousand steps as the “burn-in” phase and built the posterior distributions with at least 5000 steps. The best-fit values were taken from the maximum-likelihood estimate of the posterior distribution for each parameter and are presented in Table 4. The uncertainties for each parameter were derived from the 16th and 84th percentiles of their respective posterior distributions.

We extracted the sizes of the disks from the best-fit values following Manara et al. (2019). The disk radius ($R_{disk}$) contains 95% of the continuum flux, and the effective radius ($R_{eff}$) contains 68% of the continuum flux. The two radii are defined as follows:

$$ \int_{0}^{R_{\text{disk}}} 2\pi I(r) r dr / F_{\text{tot}} = 0.95, \qquad (7) $$

$$ \int_{0}^{R_{\text{eff}}} 2\pi I(r) r dr / F_{\text{tot}} = 0.68, \qquad (8) $$

where $I(r)$ was calculated using the best-fit parameters and $F_{tot}$ is the total flux of the continuum emission, defined from

$$ F_{\text{tot}} = I_0 \cos i \int_{0}^{\infty} 2\pi r r^{-\gamma_1} \exp \left( -\frac{r}{R_{\text{vis}}} \right)^{-\gamma_2} dr. \qquad (9) $$

Finally, we corrected $F_{tot}$ based on the inclination of the disk. The calculated radii and $F_{tot}$ values are also presented in Table 4.

In the case of barely resolved disks, degeneracies would cause the MCMC to not converge into a solution. Hence, for the three sources unresolved in the Gaussian fitting (i.e., V899 Mon, AR 6A, and AR 6C), we fixed the exponents to $\gamma_1 = 0$ and $\gamma_2 = 50$ to describe those disks with a flat surface profile whose emission decreases steeply after a certain radius. We examined the effects of this decision by fitting one of the well-resolved sources (Bran 76) with the exponents fixed to these values and compared our results to those obtained with the run where all parameters were free. We found that $i$, PA, and $F_{tot}$ are the same within their uncertainties, and $R_{eff}$ and $R_{disk}$ are smaller for the run with the fixed exponents. Therefore, we expect that in the case of the unresolved sources, the true disk radius lies between the $R_{eff}$ from the visibility fitting and the upper limit estimated from the Gaussian fitting.

In the case of Haro 5a IRS, the fitting would not converge into a solution. To overcome this, we fixed both angles ($i$ and PA) for the two disks to the values obtained in the Gaussian fitting and obtained a solution.

The cleaned images for the data, model, and residuals, as well as the observed and model visibilities as a function of uv distance, are plotted in Figures 5–14 in Appendix A.
<table><caption>**Table 4**
Results of Visibility Fitting</caption><thead><tr><th>Name</th><th>log<sub>10</sub>I<sub>0</sub><br>(Jy sr<sup>-1</sup>)</th><th>R<sub>vis</sub><br>(arcsec)</th><th>i<br>(deg)</th><th>PA<br>(deg)</th><th>γ<sub>1</sub></th><th>γ<sub>2</sub></th><th>F<sub>tot</sub><br>(mJy)</th><th>R<sub>eff</sub><br>(au)</th><th>R<sub>disk</sub><br>(au)</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">FUors</td></tr><tr><td>Haro 5a IRS</td><td>17.17<sup>+0.10</sup><sub>-0.11</sub></td><td>3.05<sup>+0.27</sup><sub>-0.38</sub></td><td>55</td><td>169</td><td>-0.86<sup>+0.02</sup><sub>-0.02</sub></td><td>54.86<sup>+4.92</sup><sub>-6.80</sub></td><td>160.65<sup>+0.17</sup><sub>-0.15</sub></td><td>33<sup>+2</sup><sub>-2</sub></td><td>38<sup>+2</sup><sub>-2</sub></td></tr><tr><td>V899 Mon</td><td>10.86<sup>+0.16</sup><sub>-0.19</sub></td><td>0.90<sup>+0.25</sup><sub>-0.21</sub></td><td>58<sup>+23</sup><sub>-38</sub></td><td>44<sup>+38</sup><sub>-21</sub></td><td>0</td><td>50</td><td>1.59<sup>+0.10</sup><sub>-0.50</sub></td><td>16<sup>+5</sup><sub>-4</sub></td><td>19<sup>+5</sup><sub>-5</sub></td></tr><tr><td>AR 6A</td><td>11.06<sup>+0.33</sup><sub>-0.27</sub></td><td>0.69<sup>+0.30</sup><sub>-0.26</sub></td><td>69<sup>+16</sup><sub>-27</sub></td><td>67<sup>+47</sup><sub>-37</sub></td><td>0</td><td>50</td><td>1.27<sup>+0.25</sup><sub>-0.86</sub></td><td>12<sup>+5</sup><sub>-4</sub></td><td>14<sup>+6</sup><sub>-5</sub></td></tr><tr><td>AR 6B</td><td>14.13<sup>+3.56</sup><sub>-3.24</sub></td><td>4.36<sup>+2.31</sup><sub>-1.34</sub></td><td>72<sup>+4</sup><sub>-6</sub></td><td>7<sup>+4</sup><sub>-4</sub></td><td>-0.52<sup>+0.47</sup><sub>-0.53</sub></td><td>106.93<sup>+48.88</sup><sub>-41.26</sub></td><td>4.87<sup>+0.14</sup><sub>-0.15</sub></td><td>43<sup>+4</sup><sub>-3</sub></td><td>49<sup>+5</sup><sub>-4</sub></td></tr><tr><td>V900 Mon</td><td>9.12<sup>+6.42</sup><sub>-5.53</sub></td><td>0.64<sup>+0.29</sup><sub>-0.24</sub></td><td>34<sup>+15</sup><sub>-16</sub></td><td>181<sup>+27</sup><sub>-21</sub></td><td>0.34<sup>+0.71</sup><sub>-0.93</sub></td><td>45.33<sup>+39.52</sup><sub>-19.90</sub></td><td>7.24<sup>+0.52</sup><sub>-0.16</sub></td><td>15<sup>+2</sup><sub>-2</sub></td><td>18<sup>+5</sup><sub>-3</sub></td></tr><tr><td>V960 Mon</td><td>14.95<sup>+9.85</sup><sub>-8.50</sub></td><td>0.83<sup>+0.79</sup><sub>-0.36</sub></td><td>58<sup>+15</sup><sub>-23</sub></td><td>88<sup>+41</sup><sub>-25</sub></td><td>-0.59<sup>+1.18</sup><sub>-1.52</sub></td><td>63.19<sup>+18.35</sup><sub>-19.93</sub></td><td>0.74<sup>+0.14</sup><sub>-0.13</sub></td><td>30<sup>+22</sup><sub>-15</sub></td><td>35<sup>+25</sup><sub>-16</sub></td></tr><tr><td>Bran 76</td><td>-0.04<sup>+3.20</sup><sub>-1.69</sub></td><td>1.75<sup>+1.94</sup><sub>-1.01</sub></td><td>52<sup>+9</sup><sub>-11</sub></td><td>52<sup>+12</sup><sub>-10</sub></td><td>1.53<sup>+0.23</sup><sub>-0.43</sub></td><td>29.11<sup>+37.64</sup><sub>-20.16</sub></td><td>5.84<sup>+2.13</sup><sub>-0.70</sub></td><td>23<sup>+5</sup><sub>-7</sub></td><td>46<sup>+18</sup><sub>-14</sub></td></tr><tr><td>V346 Nor</td><td>-0.39<sup>+0.10</sup><sub>-0.09</sub></td><td>4.24<sup>+12.57</sup><sub>-3.17</sub></td><td>25<sup>+7</sup><sub>-7</sub></td><td>2<sup>+5</sup><sub>-2</sub></td><td>1.57<sup>+0.01</sup><sub>-0.01</sub></td><td>0.30<sup>+0.76</sup><sub>-0.24</sub></td><td>46.90<sup>+56.17</sup><sub>-17.77</sub></td><td>707<sup>+752</sup><sub>-567</sub></td><td>882<sup>+690</sup><sub>-652</sub></td></tr><tr><td>OO Ser</td><td>6.56<sup>+1.04</sup><sub>-1.24</sub></td><td>3.12<sup>+1.90</sup><sub>-1.07</sub></td><td>31<sup>+7</sup><sub>-7</sub></td><td>32<sup>+7</sup><sub>-8</sub></td><td>0.58<sup>+0.18</sup><sub>-0.16</sub></td><td>43.30<sup>+29.88</sup><sub>-15.98</sub></td><td>15.57<sup>+0.46</sup><sub>-0.32</sub></td><td>30<sup>+3</sup><sub>-3</sub></td><td>37<sup>+5</sup><sub>-4</sub></td></tr><tr><td>HBC 687</td><td>24.01<sup>+0.44</sup><sub>-0.41</sub></td><td>0.35<sup>+0.03</sup><sub>-0.03</sub></td><td>56<sup>+0</sup><sub>-0</sub></td><td>74<sup>+0</sup><sub>-0</sub></td><td>-1.77<sup>+0.06</sup><sub>-0.06</sub></td><td>20.58<sup>+1.62</sup><sub>-1.59</sub></td><td>31.19<sup>+0.03</sup><sub>-0.03</sub></td><td>13<sup>+1</sup><sub>-1</sub></td><td>14<sup>+1</sup><sub>-1</sub></td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">Neighboring Sources</td></tr><tr><td>Haro 5a IRS B</td><td>16.99<sup>+0.17</sup><sub>-0.27</sub></td><td>2.17<sup>+0.52</sup><sub>-0.45</sub></td><td>56</td><td>97</td><td>-1.04<sup>+0.04</sup><sub>-0.03</sub></td><td>32.29<sup>+7.76</sup><sub>-6.59</sub></td><td>12.32<sup>+0.16</sup><sub>-0.16</sub></td><td>42<sup>+3</sup><sub>-3</sub></td><td>47<sup>+3</sup><sub>-3</sub></td></tr><tr><td>AR 6C</td><td>10.33<sup>+0.63</sup><sub>-0.43</sub></td><td>0.94<sup>+0.56</sup><sub>-0.50</sub></td><td>39<sup>+29</sup><sub>-26</sub></td><td>117<sup>+44</sup><sub>-60</sub></td><td>0</td><td>50</td><td>0.72<sup>+0.11</sup><sub>-0.17</sub></td><td>16<sup>+9</sup><sub>-8</sub></td><td>19<sup>+11</sup><sub>-10</sub></td></tr><tr><td>AR 6D</td><td>40.10<sup>+1.13</sup><sub>-2.51</sub></td><td>5.97<sup>+3.37</sup><sub>-2.48</sub></td><td>84<sup>+3</sup><sub>-2</sub></td><td>44<sup>+1</sup><sub>-2</sub></td><td>-4.77<sup>+0.40</sup><sub>-0.17</sub></td><td>126.97<sup>+74.33</sup><sub>-56.13</sub></td><td>3.94<sup>+0.22</sup><sub>-0.18</sub></td><td>104<sup>+5</sup><sub>-5</sub></td><td>109<sup>+5</sup><sub>-5</sub></td></tr><tr><td>V960 Mon N</td><td>11.11<sup>+14.32</sup><sub>-9.21</sub></td><td>4.20<sup>+2.75</sup><sub>-1.41</sub></td><td>41<sup>+9</sup><sub>-10</sub></td><td>22<sup>+13</sup><sub>-13</sub></td><td>-0.79<sup>+1.63</sup><sub>-2.11</sub></td><td>38.13<sup>+24.41</sup><sub>-18.66</sub></td><td>1.19<sup>+0.35</sup><sub>-0.19</sub></td><td>237<sup>+178</sup><sub>-70</sub></td><td>273<sup>+255</sup><sub>-82</sub></td></tr><tr><td>V960 Mon E</td><td>27.32<sup>+5.83</sup><sub>-6.67</sub></td><td>4.69<sup>+2.15</sup><sub>-1.97</sub></td><td>60<sup>+10</sup><sub>-13</sub></td><td>70<sup>+10</sup><sub>-11</sub></td><td>-3.48<sup>+1.21</sup><sub>-1.20</sub></td><td>72.76<sup>+33.67</sup><sub>-34.27</sub></td><td>2.04<sup>+0.39</sup><sub>-2.04</sub></td><td>239<sup>+124</sup><sub>-53</sub></td><td>255<sup>+134</sup><sub>-58</sub></td></tr><tr><td>V960 Mon SE</td><td>37.50<sup>+2.55</sup><sub>-9.37</sub></td><td>6.03<sup>+0.78</sup><sub>-0.46</sub></td><td>80<sup>+1</sup><sub>-1</sub></td><td>152<sup>+1</sup><sub>-1</sub></td><td>-4.53<sup>+1.51</sup><sub>-0.41</sub></td><td>84.03<sup>+9.63</sup><sub>-16.34</sub></td><td>6.62<sup>+0.28</sup><sub>-0.24</sub></td><td>342<sup>+47</sup><sub>-38</sub></td><td>361<sup>+52</sup><sub>-40</sub></td></tr><tr><td>V346 Nor B</td><td>38.95<sup>+0.15</sup><sub>-0.42</sub></td><td>13.63<sup>+12.81</sup><sub>-2.24</sub></td><td>85<sup>+0</sup><sub>-0</sub></td><td>150<sup>+10</sup><sub>-10</sub></td><td>-4.97<sup>+0.07</sup><sub>-0.02</sub></td><td>115.13<sup>+99.19</sup><sub>-18.18</sub></td><td>10.39<sup>+0.34</sup><sub>-10.39</sub></td><td>290<sup>+286</sup><sub>-212</sub></td><td>307<sup>+303</sup><sub>-222</sub></td></tr><tr><td>V370 Ser</td><td>1.75<sup>+4.86</sup><sub>-9.90</sub></td><td>4.13<sup>+2.62</sup><sub>-2.37</sub></td><td>81<sup>+6</sup><sub>-9</sub></td><td>17<sup>+17</sup><sub>-9</sub></td><td>0.70<sup>+0.65</sup><sub>-3.60</sub></td><td>33.69<sup>+16.25</sup><sub>-22.23</sub></td><td>1.11<sup>+0.71</sup><sub>-1.11</sub></td><td>52<sup>+136</sup><sub>-27</sub></td><td>69<sup>+186</sup><sub>-38</sub></td></tr></tbody></table>

Note. The indicated parameters are the most probable values, while the uncertainties for each parameter were derived from the 16th and 84th percentiles of their respective posterior distributions.

## 4. RT Modeling

To constrain the disk parameters (especially the mass) of each system more accurately than with a simple Gaussian or analytical disk model fitting, we adopted an RT model-fitting approach similar to Cieza et al. (2018). We use the RT code RADMC-3D 0.41 (Dullemond et al. 2012) with the Python interface radmc3dPy<fnref n="21" /> to set the code input parameters for a given model. As the luminosity of the central heating source in young outbursting systems may not be well represented by the stellar photosphere temperature, we use the same approach as Cieza et al. (2018); we set the stellar blackbody temperature to 10,000 K and allow the stellar radius, $R_*$, and therefore the luminosity, to be a free parameter in the model-fitting procedure.<fnref n="22" /> This will allow us to directly compare our results with those that Cieza et al. (2018) obtained for V883 Ori, HBC 494, V1647 Ori, and V2775 Ori.

The density profile of a trial protoplanetary disk model is given by

$$ \rho = \frac{\Sigma(r, \phi)}{H_p \sqrt{(2\pi)}} \exp\left(-\frac{z^2}{2H_p^2}\right), \qquad (10) $$

where $\Sigma$ is the surface density profile, and $H_p$ is the pressure scale height. The surface density profile includes a power-law inner disk and an exponential taper at large disk radii

(Andrews et al. 2009),

$$ \Sigma(r) = \Sigma_c \left(\frac{r}{R_c}\right)^\gamma \exp\left\{-\left(\frac{r}{R_c}\right)^{2+\gamma}\right\}, \qquad (11) $$

where $R_c$ is the characteristic radius of the disk, $\gamma$ is the power-law exponent of the radial surface density profile, and $\Sigma_c$ is a surface density normalization factor with

$$ \Sigma_c = (2 - \gamma) \frac{M_{\text{disk}}}{2\pi R_c^2}. \qquad (12) $$

The pressure scale height is defined as

$$ H_p = h_c \left(\frac{r}{R_c}\right)^\psi, \qquad (13) $$

where $h_c$ is the ratio of the pressure scale height over the radius at $R_c$ and $\psi$ is the degree of flaring for the disk. We used OpacityTool<fnref n="23" /> (Toon & Ackerman 1981; Woitke et al. 2016) to generate the dust absorption and scattering parameters. We assumed isotropic scattering. We include 0.1–3000 µm silicates and a size distribution of –3.5 as input to calculate the dust opacities with a volume mixture of 60% amorphous silicates (e.g., Dorschner et al. 1995) and 15% amorphous carbon (e.g., Zubko et al. 1996) and a 25% porosity. We further assume a canonical dust-to-gas mass ratio of 0.01.

To converge on the best-fit model, we use a Metropolis-Hastings MCMC model-fitting approach. The free parameters

<footnote n="21">21 https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/manual_rmcpy/index.html</footnote>

<footnote n="22">22 For AR 6C, L1551 IRS 5 N, V899 Mon, V900 Mon, and V960 Mon (all components except S), we found that using a stellar temperature of 5000 K allowed for a better convergence on a best-fit model.</footnote>

<footnote n="23">23 The OpacityTool Software was obtained from https://dianaproject.wp.st-andrews.ac.uk/data-results-downloads/fortran-package/.</footnote>
<table><caption>**Table 5**
Results of RADMC-3D Fitting</caption><thead><tr><th>Name</th><th>R<sub>*</sub><br>(R<sub>☉</sub>)</th><th>M<sub>disk</sub><br>(M<sub>☉</sub>)</th><th>ψ</th><th>R<sub>c</sub><br>(au)</th><th>γ</th><th>h<sub>c</sub><br>(R<sub>c</sub>)</th><th>i<br>(deg)</th><th>PA<br>(deg)</th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-weight:bold;">FUors</td></tr><tr><td>L1551 IRS 5 N</td><td>3.57<br>[3.38, 5.05]</td><td>0.24<br>[0.16, 0.28]</td><td>1.99<br>[1.61, 1.99]</td><td>16<br>[15, 18]</td><td>-0.00038<br>[-0.39, -0.00018]</td><td>0.97<br>[0.82, 0.98]</td><td>79<br>[69, 83]</td><td>62<br>[56, 108]</td></tr><tr><td>Haro 5a IRS</td><td>7.35<br>[4.36, 7.37]</td><td>0.38<br>[0.36, 0.61]</td><td>0.52<br>[0.25, 1.34]</td><td>69<br>[57, 79]</td><td>-0.31<br>[-1.08, -0.17]</td><td>0.49<br>[0.44, 0.78]</td><td>68<br>[54, 76]</td><td>155<br>[150, 170]</td></tr><tr><td>V899 Mon</td><td>6.39<br>[2.42, 6.57]</td><td>0.030<br>[0.018, 0.056]</td><td>0.21<br>[0.21, 1.52]</td><td>21<br>[14, 40]</td><td>-0.32<br>[-1.56, -0.28]</td><td>0.20<br>[0.20, 0.77]</td><td>41<br>[14, 53]</td><td>75<br>[46, 115]</td></tr><tr><td>AR 6A</td><td>2.61<br>[2.02, 3.04]</td><td>0.080<br>[0.038, 0.10]</td><td>0.63<br>[0.39, 1.27]</td><td>29<br>[27, 51]</td><td>-0.63<br>[-1.45, -0.29]</td><td>0.13<br>[0.11, 0.71]</td><td>75<br>[22, 77]</td><td>83<br>[55, 90]</td></tr><tr><td>AR 6B</td><td>3.01<br>[1.83, 7.67]</td><td>0.011<br>[0.0052, 0.024]</td><td>1.37<br>[0.45, 1.73]</td><td>20<br>[14, 48]</td><td>-0.42<br>[-1.34, -0.23]</td><td>0.73<br>[0.15, 0.81]</td><td>84<br>[14, 79]</td><td>71<br>[26, 162]</td></tr><tr><td>V900 Mon</td><td>3.66<br>[3.43, 6.76]</td><td>0.30<br>[0.24, 0.47]</td><td>0.53<br>[0.37, 1.55]</td><td>22<br>[18, 40]</td><td>-0.11<br>[-1.52, -0.04]</td><td>0.82<br>[0.44, 0.85]</td><td>14<br>[10, 46]</td><td>89<br>[66, 121]</td></tr><tr><td>V960 Mon</td><td>1.69<br>[1.68, 6.45]</td><td>0.17<br>[0.066, 0.27]</td><td>0.63<br>[0.31, 1.60]</td><td>46<br>[32, 89]</td><td>-0.64<br>[-1.39, -0.32]</td><td>0.48<br>[0.21, 0.80]</td><td>15<br>[12, 69]</td><td>96<br>[34, 139]</td></tr><tr><td>Bran 76</td><td>1.09<br>[1.05, 2.14]</td><td>0.16<br>[0.089, 0.23]</td><td>1.58<br>[0.42, 1.65]</td><td>27<br>[19, 41]</td><td>-0.077<br>[-0.70, -0.068]</td><td>0.29<br>[0.27, 0.73]</td><td>20<br>[17, 68]</td><td>102<br>[46, 119]</td></tr><tr><td>V346 Nor</td><td>1.95<br>[1.59, 2.21]</td><td>0.18<br>[0.15, 0.23]</td><td>0.30<br>[0.098, 0.61]</td><td>49<br>[48, 104]</td><td>-0.15<br>[-0.48, -0.082]</td><td>0.41<br>[0.33, 0.84]</td><td>20<br>[10, 30]</td><td>105<br>[90, 107]</td></tr><tr><td>OO Ser</td><td>2.80<br>[2.42, 5.10]</td><td>0.025<br>[0.022, 0.048]</td><td>0.21<br>[0.20, 1.59]</td><td>25<br>[22, 46]</td><td>-0.1<br>[-1.40, -0.19]</td><td>0.78<br>[0.51, 0.90]</td><td>74<br>[16, 74]</td><td>62<br>[-34, 177]</td></tr><tr><td>HBC 687</td><td>2.61<br>[2.76, 6.52]</td><td>0.053<br>[0.040, 0.063]</td><td>1.36<br>[0.53, 1.69]</td><td>17<br>[13, 24]</td><td>-0.21<br>[-1.48, -0.20]</td><td>0.67<br>[0.32, 0.80]</td><td>37<br>[21, 70]</td><td>66<br>[40, 130]</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">Neighboring Sources</td></tr><tr><td>L1551 IRS 5 S</td><td>3.58<br>[2.11, 6.74]</td><td>0.028<br>[0.012, 0.039]</td><td>0.32<br>[0.30, 1.62]</td><td>11<br>[9, 20]</td><td>-1.26<br>[-1.64, -0.33]</td><td>0.16<br>[0.15, 0.79]</td><td>44<br>[13, 50]</td><td>121<br>[64, 122]</td></tr><tr><td>V960 Mon N</td><td>5.24<br>[1.84, 6.65]</td><td>0.13<br>[0.081, 0.30]</td><td>1.37<br>[0.51, 1.69]</td><td>77<br>[50, 156]</td><td>-0.42<br>[-1.52, -0.32]</td><td>0.63<br>[0.19, 0.80]</td><td>60<br>[18, 74]</td><td>12<br>[2, 140]</td></tr><tr><td>V960 Mon S</td><td>7.20<br>[4.25, 12.73]</td><td>0.29<br>[0.20, 0.52]</td><td>1.52<br>[0.85, 1.67]</td><td>292<br>[132, 318]</td><td>-1.65<br>[-1.78, -1.18]</td><td>0.29<br>[0.23, 0.62]</td><td>54<br>[44, 68]</td><td>168<br>[165, 183]</td></tr><tr><td>V960 Mon E</td><td>2.09<br>[1.32, 6.63]</td><td>0.19<br>[0.11, 0.28]</td><td>1.26<br>[0.38, 1.69]</td><td>74<br>[58, 76]</td><td>-0.63<br>[-1.49, -0.32]</td><td>0.70<br>[0.124, 0.87]</td><td>19<br>[12, 66]</td><td>91<br>[73, 110]</td></tr><tr><td>V346 Nor B</td><td>2.85<br>[2.76, 6.52]</td><td>0.55<br>[0.33, 0.77]</td><td>0.62<br>[0.26, 1.16]</td><td>145<br>[136, 163]</td><td>-1.05<br>[-1.32, -0.45]</td><td>0.48<br>[0.34, 0.79]</td><td>66<br>[57, 74]</td><td>142<br>[134, 154]</td></tr></tbody></table>

Note. The indicated parameters are the most probable values, while the ranges in the square brackets show the 68% credible region for each parameter.

considered in the modeling are stellar radius, $R_*$; total disk mass, $M_{\text{disk}}$; degree of disk flaring, $\psi$; characteristic radius, $R_c$; power-law exponent of the surface density profile, $\gamma$; scale height ratio, $h_c$; disk inclination, $i$; and disk PA. Given the detailed calculations in RADMC-3D, a given trial model can become relatively computationally expensive, depending on the trial disk parameters. Therefore, we perform the MCMC modeling in the image plane, as we can more quickly converge to the most probable values of the free parameters (see, e.g., Booth et al. 2016; White et al. 2016). After a trial model is computed in RADMC-3D, a 1.3 mm continuum image is then produced and projected to a trial inclination and PA. The trial model is then attenuated by the primary beam and convolved with the synthetic beam for a given observational setup. To assess the likelihood of a given model, a $\chi^2$ is calculated as

$$ \chi^2 = \frac{(\text{Data} - \text{Model})^2}{\sigma^2}, \qquad (14) $$

where $\sigma$ is the observed $\sigma_{\text{rms}}$ for a given observation multiplied by the synthetic beam size in pixels (see Booth et al. 2016). A trial model is accepted if a random number drawn from a

uniform distribution [0, 1] is less than $\alpha$, where

$$ \alpha = \min(e^{\frac{1}{2}(\chi_i^2 - \chi_{i+1}^2)}, 1). \qquad (15) $$

For each system, we ran 10 chains with 1000 links each (minus 100 for burn-in). The most probable values and 68% credible regions are summarized in Table 5. If the posterior distributions were Gaussian, these would correspond to the $1\sigma$ uncertainties. The data minus model residual plots for each system are in Figures 15–25 in Appendix B. We fit all FUors in our sample and most of the nearby sources detected in our images, with the exception of Haro 5a IRS B, AR 6C, AR 6D, and V370 Ser due to their low S/N. In Figures 26–28 in Appendix C, we show some examples for the posterior distributions.

## 5. Results on Individual Objects

In the following, we discuss each of our targets individually. We first summarize the earlier results in relation to their circumstellar properties, then describe our new ALMA results in the context of the literature. To compare the results of the different types of modeling presented in this paper, we refer to Appendix D. The geometrical parameters from the different
<table><caption>**Table 6**
Mass Fraction of the RT-generated Disk that Is Optically Thick ($\tau_{1.3mm} > 1.0$) and Midplane Temperature at the Characteristic Radius</caption><thead><tr><th>Name</th><th><i>M</i><sub>disk</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>τ>1.0</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>τ>1.0</sub>/<i>M</i><sub>disk</sub></th><th><i>T</i><sub>RC</sub><br>(K)</th></tr></thead><tbody><tr><td>L1551 IRS 5 N</td><td>0.24</td><td>0.12</td><td>0.48</td><td>130</td></tr><tr><td>Haro 5a IRS</td><td>0.38</td><td>0.10</td><td>0.27</td><td>125</td></tr><tr><td>V899 Mon</td><td>0.030</td><td>0.0093</td><td>0.31</td><td>51</td></tr><tr><td>AR 6A</td><td>0.080</td><td>0.041</td><td>0.51</td><td>60</td></tr><tr><td>AR 6B<sup>a</sup></td><td>0.011</td><td>...</td><td>&lt;0.01</td><td>210</td></tr><tr><td>V900 Mon</td><td>0.30</td><td>0.19</td><td>0.64</td><td>110</td></tr><tr><td>V960 Mon</td><td>0.17</td><td>0.0093</td><td>0.055</td><td>33</td></tr><tr><td>Bran 76</td><td>0.16</td><td>0.065</td><td>0.41</td><td>82</td></tr><tr><td>V346 Nor</td><td>0.18</td><td>0.050</td><td>0.28</td><td>68</td></tr><tr><td>OO Ser<sup>a</sup></td><td>0.025</td><td>...</td><td>&lt;0.01</td><td>220</td></tr><tr><td>HBC 687</td><td>0.053</td><td>0.014</td><td>0.26</td><td>240</td></tr></tbody></table>

**Note.**

<sup>a</sup> These RT-generated disks were completely optically thin, likely due to the combination of large scale heights and lower relative mass; therefore, we write the optically thick mass fraction as <0.01.

models are generally in good agreement within the uncertainties (Figure 29). However, there are significant differences between the masses obtained from the RT modeling (Table 5) and the estimates that simply convert the flux to mass (Table 3).

To understand and quantify the differences, we calculated the mass that is within the optically thick part ($\tau > 1.0$) of the disk from the best-fitting RT models for each of our targets, along with the optically thick mass fraction and the midplane temperature at the characteristic radius (Table 6). The RT masses for the FUor sample are a factor of 1.4–34 higher than the optically thin masses, i.e., masses obtained from the fluxes assuming 50 K, and a factor of 5.7–146 higher when using the temperatures from the bolometric luminosities. While there is no one-on-one correspondence between the optically thick mass fraction and the ratio of the RT masses to optically thin masses, they do generally correlate; the disks with the highest/lowest optically thick mass fraction have the highest/lowest ratio between the RT and optically thin masses. These indicate that (at least part of) the disks are optically thick at 1.33 mm, and the RT masses are higher because the ray tracing takes into account optical depth effects. Therefore, the RT masses are probably closer to the real disk masses, and in the following, we will use them in the discussion.

### 5.1. L1551 IRS 5

L1551 IRS 5 is a class 0/I FUor-like embedded object in the Taurus star-forming region. The system is a 365 mas (51 au) protobinary surrounded by circumstellar disks, circumbinary material (Cruz-Sáenz de Miera et al. 2019, and references therein), and a massive envelope (Chou et al. 2014). It drives a bipolar outflow (Wu et al. 2009), and two jets are detected at radio (Rodríguez et al. 1998) and optical/near-infrared wavelengths (Fridlund et al. 2005; Pyo et al. 2005, and references therein).

Our observations recovered the two protostellar disks and resolved the circumbinary emission. The different components of the system were already analyzed by Cruz-Sáenz de Miera et al. (2019), and here we will use their reported disk properties and compare them with the rest of our sample. For the Gaussian fitting, Cruz-Sáenz de Miera et al. (2019) separated the

northern disk into two components to avoid contamination from circumbinary material, and their fitted parameters for both disks match those from 7 mm Very Large Array observations with higher angular resolution (Lim et al. 2016). For the RT fitting, we used a single component to model each protostellar disk, and, in the case of the northern disk, this causes differences in the inclination and PA between the Gaussian fitting and the RT. For the total disk masses in Table 3, we made different assumptions on the disk temperature than Cruz-Sáenz de Miera et al. (2019). They used the brightness temperature as a proxy for the disk temperature and obtained 160 and 100 K for the northern and southern disks, respectively, whereas here we estimated the disk temperature from the bolometric luminosity and obtained 101 and 77 K, respectively.

Taking into account both our and the literature data, we can conclude that the two circumstellar disks have the same inclinations and PAs, and our best estimates for them are $i = 44^{\circ}$ and $PA = 158^{\circ}$. These disks are the smallest in our sample. This is not surprising if we consider that this is the closest separation binary in our sample and dynamical effects might have truncated the circumstellar disks. From our RT modeling, we obtained a disk mass of $0.24 M_{\odot}$ for the disk of the northern component. Interestingly, the southern component, which is not outbursting, has a disk that is almost an order of magnitude less massive.

### 5.2. Haro 5a IRS

Haro 5a IRS is an embedded FUor-like object in the Orion Nebula Cluster, located between two sides of a bipolar nebula, Haro 5a to the west and Haro 6a to the east (Reipurth & Aspin 1997). Haro 5a and 6a are separated by a high extinction ridge that hides the young stellar object (YSO) at optical wavelengths, but Haro 5a IRS is well detected at longer wavelengths, e.g., in the mid- and far-infrared (IRAS 05329 $-0505$, MIR 10, Nielbock et al. 2003; HOPS 84, Fischer et al. 2010) and the (sub)millimeter (OMC3 MMS 7, Chini et al. 1997; OMC3 SMM 11, Takahashi et al. 2013).

Haro 5a IRS is surrounded by a large amount of circumstellar material. ALMA 3 mm continuum observations by Kainulainen et al. (2017) suggest a total mass of $1.9 M_{\odot}$ (using 20 K, $\kappa_{\nu} = 0.1(\nu/1000 \text{ GHz})^{\beta} \text{ cm}^2 \text{ g}^{-1}$, $\beta = 1.5$, and a gas-to-dust ratio of 100). Using optically thin $^{13}$CO $J=3-2$ line observations, Kóspál et al. (2017b) obtained $1.2 M_{\odot}$. Considering the uncertainties in dust opacities and optical depth, these numbers are fairly consistent with each other.

In accordance with the bipolar nebula nature of the object, it is not surprising that our ALMA data show that the disk of Haro 5a IRS is inclined ($55^{\circ}-68^{\circ}$). Based on our RT model, we obtained a total disk mass of $0.38 M_{\odot}$, while its characteristic radius is 69 au. Both of these parameters are the largest in our FUor sample.

We detected a nearby source at a distance of $0.770$ to the southwest ($PA = 209^{\circ}$) of Haro 5a IRS, which has not been reported earlier in the literature. We call this source Haro 5a IRS B. Our ALMA line data indicate CO emission from this source at approximately the same velocities as Haro 5a IRS itself (Á. Kóspál et al. 2021, in preparation), suggesting that they form a physical pair. Their projected separation corresponds to 273 au. Haro 5a IRS B is fainter and not as well resolved as Haro 5a IRS; therefore, it must correspond to a less massive disk.
### 5.3. V899 Mon

Ninan et al. (2015) classified V899 Mon as a young eruptive star with photometric and spectroscopic properties that fall between those of classical FUors and EXors. They compiled the spectral energy distribution (SED) of the object and modeled it using the SED fitter tool of Robitaille et al. (2007), considering data points only up to 70 µm. They found a disk mass between $10^{-5}$ and $10^{-3} M_\odot$ and an envelope mass between $10^{-7}$ and $10^{-5} M_\odot$, depending on whether they used quiescent or outburst photometry. Because the SED of V899 Mon shows a curious dip at 70 µm and increases again until 160 µm, Ninan et al. (2015) interpreted this as a separate clump, possibly an envelope component with a large inner cavity. They obtained 20–23 $M_\odot$ for the mass of this component. Following Kóspál et al. (2017b), we used optically thin $^{13}$CO line observations from APEX to convert the observed line fluxes to total gas mass assuming local thermodynamic equilibrium and 20 K temperature, integrated over a 10,000 au radius area. We obtained 0.24 $M_\odot$, which is between the values published for the disk and the far-infrared envelope component by Ninan et al. (2015).

Our ALMA observations suggest that V899 Mon is among the less massive FUor disks in our sample; only AR 6B and OO Ser have smaller disk masses. It has a moderate characteristic radius. Although the Gauss fitter in casaviewer found the object to be a point source, our visibility space fitting and RT modeling could measure its size. As opposed to Ninan et al. (2016), who derived an inclination of 10°, i.e., an almost face-on geometry from optical spectroscopy, our ALMA data show an intermediate inclination (41°–58°). The small mass we measure with ALMA is in contrast with the massive envelope suggested by the far-infrared and submillimeter SED shape. Although our ALMA data recover all emission up to 5"/9 (or about 4500 au), we do not detect the suspected envelope component. This may be a sensitivity issue if the envelope has a very low surface brightness, or this component may be too extended and is filtered out even with the ALMA 7 m array. Indeed, Ninan et al. (2015) wrote that the source's FWHM is 50" in the Herschel/SPIRE 500 µm image. Since the FWHM of SPIRE at this wavelength is only 35"/2, this hints at an envelope that extends at least about 35" or several times 10,000 au.

### 5.4. AR 6A/6B

Although no eruption was observed for them, both AR 6A and AR 6B were classified as FUors based on their spectroscopic similarities to FU Ori and PP 13S* by Aspin & Reipurth (2003). However, Connelley & Reipurth (2018) noticed that their spectrum of AR 6A, obtained more than 10 yr after that of Aspin & Reipurth (2003), shows fewer of the spectral characteristics required for a FUor-like classification, suggesting significant spectral variability. Connelley & Reipurth (2018) did not find convincing evidence to classify AR 6B as a FUor-like object either. They classify both AR 6A and AR 6B as peculiar objects.

The objects are located in the NGC 2264 cluster, which is part of the Mon OB1 East giant molecular cloud. Forbrich et al. (2010) presented a Spitzer/IRS spectrum for AR 6A and found ice absorption features and strong polycyclic aromatic hydrocarbon emission in the 5–38 µm range. They constructed the SED of the object and tried to fit it using the SED fitter tool of

Robitaille et al. (2007), but they found that it was not well fitted by the models and did not give model parameters. On the other hand, Gramajo et al. (2014) could fit both AR 6A and AR 6B with Robitaille et al. (2007) models and found the following disk and envelope masses: for AR 6A, $M_{disk} = 0.34 M_\odot$, $M_{env} = 0.20 M_\odot$; for AR 6B, $M_{disk} = 0.37 M_\odot$, $M_{env} = 0.01 M_\odot$. Using APEX CO observations, Kóspál et al. (2017b) found that there is no distinct peak in the CO emission at the stellar position, concluded that the CO emission is most likely not associated with these sources, and could only give an upper limit of 1.3 $M_\odot$ for the total gas mass using the optically thin $^{13}$CO(3–2) line.

We clearly detect both AR 6A and AR 6B in our ALMA images. At 1.33 mm, AR 6B is the brightest target; then comes AR 6A, and even AR 6C is detected at an ~7σ level. The disks of AR 6A and AR 6B are not particularly large or massive, compared to the rest of our FUor sample. In fact, AR 6B has the least massive disk within our FUor sample. We measure separations of 2.7"/1 (A–B), 0.8"/9 (A–C), and 2.2"/18 (B–C). For comparison, Aspin & Reipurth (2003) measured a separation of 2.2"/8 between A and B and 0.2"/85 between A and C in Subaru near-infrared *JHK* images. We detected an additional source 4.4"/43 to the southwest (PA = 244°) of AR 6B at an ~40σ level, which we have labeled AR 6D.

### 5.5. V900 Mon

V900 Mon went into eruption some time between 1953 and 2009 (Reipurth et al. 2012). Based on this, as well as various spectroscopic characteristics, it is classified as a bona fide FUor by Connelley & Reipurth (2018). According to Reipurth et al. (2012), V900 Mon is a class I–type YSO embedded in a cool envelope. Its SED, if corrected for $A_V = 13$ mag extinction, is remarkably similar to that of FU Ori. Kóspál et al. (2017b) measured the object with APEX in different CO lines but found no distinct CO peak toward the stellar position, suggesting a rather low-mass envelope. Therefore, an upper limit of 0.1 $M_\odot$ for the total mass of the circumstellar material was given. Based on modeling the SED with Robitaille et al. (2007) models, Gramajo et al. (2014) found a disk mass of 0.2 $M_\odot$ and an envelope mass of 0.072 $M_\odot$.

Based on ALMA data, Takami et al. (2019) analyzed the morphology and kinematics of the circumstellar matter around V900 Mon at a $10^4$ au spatial scale. The $^{12}$CO emission indicated an extended bipolar outflow in the east–west direction. The $^{13}$CO emission traced compressed gas in the cavity wall. The C$^{18}$O emission revealed a rotating envelope.

We clearly detected V900 Mon in our ALMA continuum image. No other millimeter sources are present in the field of view of the primary beam. With a mass of 0.3 $M_\odot$ and characteristic radius of 22 au, it has a small and massive disk. The disk is seen under a modest inclination between 14° and 34°, depending on the type of model fitting. The Gaussian and visibility space fitting suggest that the orientation of the major axis of the disk is close to north–south, which means that the molecular outflow seen by Takami et al. (2019) is perpendicular to the disk.

### 5.6. V960 Mon

V960 Mon went into outburst in 2014 and shows all of the necessary spectral properties of a bona fide FUor (Kóspál et al. 2015 and references therein; Connelley & Reipurth 2018). Based on periodic oscillations in the outburst light curve,
Hackstein et al. (2015) suggested that it is a close binary. Kóspál et al. (2015) studied the SED of the progenitor and found that it indicates the presence of a moderately flared disk and an outer envelope in the system, whose total mass is between 0.1 and $0.7 M_\odot$, if we scale it to the new Gaia distance of 1574 pc from the previously assumed 450 pc.

We detect several millimeter sources in our ALMA image. Object V960 Mon itself is clearly detected, just like V960 Mon N (5.53 away at PA = 15°) and V960 Mon SE (8.45 away at PA = 152°), two nearby sources identified previously in infrared images (Kóspál et al. 2015). None of the T Tauri candidates listed in Kóspál et al. (2015) are detected in the millimeter. One more millimeter source to the east of V960 Mon, V960 Mon E (5.36 away at PA = 97°), has no known infrared counterpart. Caratti o Garatti et al. (2015) reported a companion at a separation of 227 mas and PA of 131° detected using VLT/SINFONI at near-infrared wavelengths, as well as a bump to the southwest of V960 Mon. In our ALMA image, we do not detect the bump. In the residuals of the visibility fitting for this target (Figure 10), there are residuals over $5\sigma$ at the same PA but slightly closer (155 mas) than the companion discovered by Caratti o Garatti et al. (2015). Considering that our ALMA observations were taken within 2 yr of the SINFONI data, it is unlikely that the millimeter emission in our residual map can be associated with this companion.

Our modeling suggests that within our FUor sample, V960 Mon has a rather large disk with an average mass. The disk's PA is close to east-west, and its inclination is between 15° and 61°, depending on the type of model fitting. APEX observations of the CO emission from this source suggest an outflow oriented along the line of sight (F. Cruz-Sáenz de Miera et al. 2021, in preparation), which favors a more face-on geometry for the disk.

## 5.7. Bran 76

Bran 76 was classified as a FUor by Eisloffel et al. (1990) based on high-resolution optical spectra. Connelley & Reipurth (2018) confirmed that its infrared characteristics are also FUor-like, although no outburst was ever observed. Reipurth et al. (2002) noted that it has been slowly fading since 1983, although photographic plates from 1900 and 1927 suggest that the fading rate must have been lower earlier. Its SED is consistent with a steady accretion disk model shortward of 10 µm. At longer wavelengths, there is excess emission due to an envelope, whose luminosity is intermediate between that of FU Ori and V1057 Cyg. Earlier APEX observations indicate a moderate-mass envelope in the system ($0.02 M_\odot$; Kóspál et al. 2017b) and no molecular outflow.

We clearly detected Bran 76 in the 1.33 mm continuum with ALMA. Our modeling suggests that it has an average disk size and mass within the observed FUor sample. The disk has an intermediate inclination. We do not see any other source in the primary beam, but the residuals in Figure 11 hint at some extended emission at a several hundred au scale, possibly from the envelope.

## 5.8. V346 Nor

The fifth object that was discovered to undergo a FUor outburst (Reipurth 1985) was V346 Nor. It brightened some time between 1976 and 1980, reached peak brightness in 1992,

showed a deep minimum in 2010–11, and is currently brightening again (Kraus et al. 2016; Kóspál et al. 2017c, 2020b). In Kóspál et al. (2017a), we analyzed earlier 1" resolution ALMA 1.3 mm continuum, as well as $^{12}$CO, $^{13}$CO, and $^{18}$O J = 2–1 line observations of V346 Nor and found that it is surrounded by a 10,000 au envelope with an outflow cavity. The central ~700 au part of the circumstellar matter forms a flattened pseudodisk where material is infalling with conserved angular momentum. Within about 350 au, the velocity profile is more consistent with a disk in Keplerian rotation. The infall rate from the envelope onto the disk is $6 \times 10^{-6} M_\odot \text{ yr}^{-1}$, which is a factor of a few higher than the quiescent accretion rate from the disk onto the star. This hints at a mismatch between the infall and accretion rates as the cause of the eruption.

From our RT modeling, we obtained a characteristic radius of 49 au, which suggests that the detected emission is dominated by the inner part of the Keplerian disk. With its $M_{\text{disk}} = 0.18 M_\odot$, the disk of V346 Nor is among the most massive FUor disks. This disk mass is on the same order of magnitude as the protostellar mass.

An earlier, lower spatial resolution ALMA continuum image (Figure 1 in Kóspál et al. 2017a) revealed faint extended emission around the star that seems to follow the outflow cavity walls. We detect extended emission in our new ALMA images, although with low S/N, especially to the northwest of the source. The extension of this faint emission causes our single-component visibility space fitting to flatten out (i.e., low value of $\gamma_2$), causing $R_{\text{eff}}$ and $R_{\text{disk}}$ to reach hundreds of au. To estimate the radius of the compact disk that can be compared with the modeling done in the image space, we run another fit using a two-component model. This model gives $R_{\text{eff}} = 21$ and $R_{\text{disk}} = 25$ au, somewhat smaller values than either the Gaussian or the RT fits. The extended component reaches out to ~300 au, which matches the ~0.5" extension seen in the continuum maps (Figures 1 and 12). We also detected V346 Nor B, a deeply embedded protostar Kóspál et al. (2017a) discovered in the earlier ALMA data set at a separation of 15.2".

## 5.9. OO Ser

OO Ser is a deeply embedded pre-main-sequence star that went into outburst in 1995 and faded back to quiescence on an ~8 yr timescale (Hodapp et al. 1996, 2012; Kóspál et al. 2007). Kóspál et al. (2007) compiled its broadband SED for different epochs and found data until 200 µm, although the longer-wavelength points suffer from source confusion due to nearby sources like V370 Ser, V371 Ser, EC 38, and SMM 9. Deep ice absorption features in the infrared spectrum suggest the presence of a thick cold envelope (Kóspál et al. 2007). The optical depth of the 9.7 µm silicate feature implies $A_V = 42$ mag. OO Ser is clearly visible in a JCMT/SCUBA 850 µm image by Yoo et al. (2017), but they did not give its flux or determine its mass. Kóspál et al. (2017b) detected CO line emission with APEX and determined a total circumstellar mass of $0.6 M_\odot$ using the optically thin $^{13}$CO J = 3–2 line flux. Although the CO line profiles showed high-velocity wings, the spatial resolution of APEX was not sufficient to disentangle which sources drive molecular outflows in this region. On the other hand, Hodapp (1999) found $H_2 S(1)$ emission to the east of OO Ser, suggesting jet activity.
We detect two millimeter sources in our image, OO Ser and V370 Ser (at a distance of 11"/78 to the north of OO Ser at a PA of 7°/6). OO Ser has a relatively low-mass, small disk within our FUor sample. It has the second-smallest disk mass after AR 6B. The disk has an intermediate inclination and a PA between 32° and 62°, depending on the type of model fitting, and OO Ser has a bipolar reflection nebula (Hodapp et al. 1996) extending in the east-west direction. If we interpret this as the outflow direction, and if the disk is perpendicular to that, then the scattered-light morphology supports the lower PA for the disk.

## 5.10. HBC 687

Although it displayed no outburst, HBC 687 (also known as Parsamian 21) was first classified as a FUor-like object based on spectroscopic properties by Staude & Neckel (1992), later confirmed by Connelley & Reipurth (2018). Kóspál et al. (2008) obtained near-infrared polarimetric images of the scattered-light nebula in the system and found a circumstellar envelope with a polar cavity and an edge-on disk with a PA of 78° ± 4°. The disk seems to be geometrically flat and extends approximately 48–360 au from the star. The SED of HBC 687 is consistent with this picture (see also Gramajo et al. 2014). APEX observations by Kóspál et al. (2017b) revealed a relatively tenuous envelope around the system, with a total mass of 0.01 $M_\odot$ measured within 10,000 au of the central star using the optically thin $^{13}$CO(3–2) line.

With ALMA, we detected a single source at the stellar position. Our modeling suggests a modest disk mass and rather small disk size within our FUor sample. While the millimeter-emitting disk seems to be more inclined than perfectly edge-on, its PA is consistent with what is seen in scattered light (Kóspál et al. 2008).

# 6. Discussion

## 6.1. Diversity of the Circumstellar Matter of FUors

The FUor envelopes are predicted to exhibit large diversity, depending on how far they are in the dispersal process during the FUor phase. Indeed, Quanz et al. (2007) and Kóspál et al. (2020a) found both FUors showing 10 µm silicate absorption (younger, more embedded objects) and FUors with silicate emission (objects with envelopes that are already partly opened up, enabling a direct view of the surface of the accretion disk). A similar evolutionary sequence was outlined by Green et al. (2006) based on the amount of far-infrared excess, which traces the dust mass in the envelope. Concerning the gas component, Kóspál et al. (2017b) found a similar division between younger and more evolved FUors based on single-dish millimeter CO data.

Envelopes play a significant role in the outburst physics by both replenishing the disk material after each outburst and sustaining gravitational instability, the key ingredient of several outburst models (Vorobyov & Basu 2005, 2010, 2015; Zhu et al. 2009; Kuffmeier et al. 2018). ALMA already made an important contribution in mapping the envelope structure of a few selected FUors, resulting in the discovery of the widest outflow cavity in a class I object known to date (HBC 494; Ruíz-Rodríguez et al. 2017) and the first clear detection of the transition from an infalling envelope to a rotating disk in a FUor (V346 Nor; Kóspál et al. 2017a). However, in such a diverse class of objects, general conclusions cannot be drawn

from data on a limited number of objects, and a complete, deep, and homogeneous survey is indispensable.

The FUor disks are also expected to display a large diversity while evolving from a class I massive accretion disk toward the passive protoplanetary phase. While it is challenging to observationally separate disks from remnant circumstellar envelopes, earlier studies in the literature based on (sub) millimeter interferometry indeed revealed very different disk masses in FUors (in the range of 0.007–0.3 $M_\odot$; see L1551 IRS 5 in Takakuwa et al. 2004, PP 13S* in Pérez et al. 2010, V883 Ori in Cieza et al. 2016, HBC 494 in Ruíz-Rodríguez et al. 2017, FU Ori in Hales et al. 2015 and Pérez et al. 2020, or V2775 Ori in Zurlo et al. 2017). As a new example, our group is studying HBC 722, a FUor where the disk remained undetected with SMA and IRAM (Dunham & Vorobyov 2012; Kóspál et al. 2016) but could be detected with ALMA, obtaining a total disk mass of 0.01 $M_\odot$ (Xi et al. 2020). This large variety of disk properties might imply different outburst mechanisms for objects in different evolutionary stages or with different disk structures.

## 6.2. Differences between FUor and Class I/II Disks

Episodic accretion during pre-main-sequence stellar evolution may be common (e.g., Audard et al. 2014). This is supported by the observational fact that accretion-related outbursts can be observed for YSOs of various evolutionary states from very embedded to disk-only. Also, accretion outbursts are observed for YSOs of various protostellar masses from $\lesssim 0.1$ (V346 Nor in Kóspál et al. 2017a or OO Ser in Hodapp et al. 1996) to $20 M_\odot$ (S2551RIRS 3 in Caratti & Garatti et al. 2017). Therefore, it is a highly relevant question whether the disk properties of eruptive YSOs are in any way different from the disks of normal YSOs not showing signs of present or past outburst. Our current study provides an opportunity to inspect this further by comparing the disk parameters we obtained from our RT modeling of FUors to samples of regular class I or II disks modeled in a similar way.

To compare disk properties between FUors and regular YSOs, we take the class II comparison sample from Andrews et al. (2010), who observed 16 T Tauri disks with the SMA in the Ophiuchus star-forming region. We chose this sample because they were fitted with the same disk profile and RT models as our sample of FUors (Section 4).

For a class I comparison sample, we took the results of Sheehan & Eisner (2014, 2017a, 2017b), who observed objects mostly in Taurus and one source in Ophiuchus. Their disk modeling employed a pure power-law surface density profile without an exponential taper. We numerically compared their surface density profiles with ours (Equation 11) and found that, using the outer disk radius $R_{\text{disk}}$ from Sheehan & Eisner (2017a) as our characteristic radius $R_c$, we would get very similar surface density profiles and identical total disk masses. Therefore, for lack of any work for class I disks using exactly the same model as here, the parameters from Sheehan & Eisner (2014, 2017a, 2017b) are the most comparable to ours.

To extend the FUor sample, we searched the literature for other FUors modeled in a way similar to ours. We found that PP 13S* from Pérez et al. (2010); V883 Ori, HBC 494, V1647 Ori, and V2775 Ori from Cieza et al. (2018); and FU Ori from Pérez et al. (2020) can be added to our sample, although for PP 13S*, we first had to calculate $R_c$, because Pérez et al. (2010) normalized the surface brightness profile in a
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Total disk mass (in $M_c$) as a function of characteristic radii. The orange points are the FUors presented here. The green points are FUors from Pérez et al. (2010, 2020) and Cieza et al. (2018). The black points are T Tauri disks from Andrews et al. (2010). The blue points are class I disks from Sheehan & Eisner (2014, 2017a, 2017b). Linear fits are $M_{disk} \propto R_c^{1.40}$ for class II and $M_{disk} \propto R_c^{0.99}$ for class I.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Histograms and cumulative distributions of disk masses and radii from RT modeling. FUors include the sample presented here and those from Pérez et al. (2010, 2020) and Cieza et al. (2018). Class II and I disks are the same as in Figure 2.</figcaption></figure>

<table><caption>**Table 7**

P-values of Various Two-sample Statistical Tests</caption><thead><tr><th>Statistical Test</th><th>FUors/Class II<br>Disk Masses</th><th>FUors/Class I<br>Disk Masses</th><th>Class II/Class I<br>Disk Masses</th><th>FUors/Class II<br>Disk Radii</th><th>FUors/Class I<br>Disk Radii</th><th>Class II/Class I<br>Disk Radii</th></tr></thead><tbody><tr><td>Kolmogorov–Smirnov test</td><td>0.0254</td><td>0.0214</td><td>0.9997</td><td>0.1967</td><td>0.0348</td><td>0.3738</td></tr><tr><td>Wilcoxon rank sum test</td><td>0.0069</td><td>0.0090</td><td>0.9410</td><td>0.0660</td><td>0.0016</td><td>0.0641</td></tr><tr><td>Welch two-sample t-test</td><td>0.0136</td><td>0.0123</td><td>0.8594</td><td>0.0252</td><td>0.0291</td><td>0.0872</td></tr></tbody></table>

different way (their $R_f = 13$ au corresponds to $R_c = 30$ au). In the following, we work with these 17 FUors, 16 class II objects, and 11 class I objects.

In Figure 2, we show the total disk mass as function of the characteristic radius. The FUors presented here, along with FUors from the literature, show a clear separation both from typical T Tauri disks and from typical class I disks in that they are more massive and compact. For normal class I/II disks, there is a significant correlation between disk mass and radius (Pearson, Spearman, and Kendall correlation statistics give coefficients in the range of 0.42–0.82, with very low probabilities of 0.0005–0.05 of the null hypothesis of no correlation). To guide the eye, we fitted a line using an ordinary least-squares bisector regression method for these comparison samples. Class II disks follow a mass–radius relation of $M_{disk} \propto R_c^{1.40}$, and class I disks are $M_{disk} \propto R_c^{0.99}$. All FUors in Figure 2 are located above these average trends for the comparison samples. In Figure 3, we plotted histograms and cumulative distributions for the masses and radii of the three samples. These graphs also suggest that FUor disks are generally more massive and smaller in size than the disks of either class II or class I objects.

To quantify whether the mass or radius of the different samples follows the same distribution (null hypothesis), we run various

two-sample tests using the R software (R Core Team 2021). Table 7 shows the resulting probabilities of the null hypothesis for the two-sample Kolmogorov–Smirnov test, the Wilcoxon rank sum exact test, and the Welch two-sample t-test. The results suggest that the mass distribution of FUor disks is significantly different from that of both class II and class I disks. The mass distributions of class II and class I disks, on the other hand, seem to be similar. This is supported by their cumulative distributions in the lower left panel of Figure 3 as well. Concerning the disk radii, all three samples are different, with the FUor disks being the most compact and the Class I disk being the largest (see also the lower right panel of Figure 3).

Our results suggest that FUor disks differ from class I or II disks in the distribution of their masses and radii. For the merged FUor sample used in Figures 2 and 3, the median FUor disk radius is 27 au, while the median FUor disk mass is 0.08 $M_⊙$. The respective values for the comparison samples are as follows: the median disk radius for class I disks is 127 au, that for class II disks is 40 au; the median disk mass for class I disks is 0.02 $M_⊙$, and that for class II disks is 0.03 $M_⊙$. Therefore, while individual FUor disks may be large and not very massive, the whole sample, on average, is more compact (by a factor of 1.5–4.7) and massive (by a factor of 2.9–4.4) than the disks around class I or II objects.
<table><caption>**Table 8**
Gravitational Instability as Calculated by Four Different Metrics</caption><thead><tr><th>Name</th><th>M<sub>disk</sub><br>(M<sub>☉</sub>)</th><th>M<sub>*</sub><br>(M<sub>☉</sub>)</th><th>M<sub>disk</sub>/M<sub>*</sub> > 0.1</th><th>M<sub>disk</sub>/M<sub>*</sub> > 0.25</th><th>M<sub>disk</sub>/M<sub>*</sub> > H/r</th><th>Q(r) < 1.4</th><th>Q<sub>min</sub></th><th>r(Q<sub>min</sub>)<br>(au)</th></tr></thead><tbody><tr><td>L1551 IRS 5 N</td><td>0.24</td><td>0.25</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>0.71</td><td>19</td></tr><tr><td>Haro 5a IRS</td><td>0.38</td><td>0.6</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>0.75</td><td>95</td></tr><tr><td>V899 Mon</td><td>0.030</td><td>2.0</td><td>N</td><td>N</td><td>N</td><td>N</td><td>7.4</td><td>17</td></tr><tr><td>AR 6A</td><td>0.080</td><td>0.6</td><td>Y</td><td>N</td><td>Y</td><td>N</td><td>2.5</td><td>24</td></tr><tr><td>AR 6B</td><td>0.011</td><td>0.6</td><td>N</td><td>N</td><td>N</td><td>N</td><td>27</td><td>22</td></tr><tr><td>V900 Mon</td><td>0.30</td><td>0.6</td><td>Y</td><td>Y</td><td>N</td><td>Y</td><td>0.74</td><td>29</td></tr><tr><td>V960 Mon</td><td>0.17</td><td>0.75</td><td>Y</td><td>N</td><td>N</td><td>Y</td><td>0.90</td><td>46</td></tr><tr><td>Bran 76</td><td>0.16</td><td>0.6</td><td>Y</td><td>Y</td><td>N</td><td>Y</td><td>1.3</td><td>27</td></tr><tr><td>V346 Nor</td><td>0.18</td><td>0.1</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>0.52</td><td>55</td></tr><tr><td>OO Ser</td><td>0.025</td><td>&lt;0.1</td><td>Y</td><td>Y</td><td>N</td><td>N</td><td>&lt;5.4</td><td>29</td></tr><tr><td>HBC 687</td><td>0.053</td><td>0.6</td><td>N</td><td>N</td><td>N</td><td>N</td><td>5.6</td><td>17</td></tr><tr><td>Total number of gravitationally unstable disks</td><td></td><td></td><td>8</td><td>6</td><td>3</td><td>6</td><td></td><td></td></tr></tbody></table>

### 6.3. Gravitational Instability in FUor Disks

The more massive and moderately more compact nature of the FUor disks could yield insight into the source of the episodic accretion that drives the outbursts. While the exact mechanisms that cause enhanced accretion of disk material onto the protostar are still unknown, it may be triggered if the disk becomes gravitationally unstable, possibly in combination with magnetorotational instability (Armitage et al. 2001; Lodato & Clarke 2004; Vorobyov & Basu 2006, 2010; Boley & Durisen 2008; Zhu et al. 2009; Martin & Lubow 2011; Bae et al. 2014; Kadam et al. 2019). This scenario most likely occurs in relatively massive disks when $M_{disk}/M_* \ge 0.1$ (Kratter & Lodato 2016). The actual criterion is closer to $M_{disk}/M_* \ge H/r$, i.e., if the $M_{disk}/M_*$ ratio is greater than or equal to the disk aspect ratio (the disk scale height over the characteristic radius). More accurately, gravitational instability occurs if the value of Toomre's $Q$ parameter is $Q = c_s \Omega / \pi G \Sigma \lesssim 1.4$, where $c_s$ is the sound speed, $\Omega$ is the epicyclic frequency for a Keplerian disk, and $\Sigma$ is the column density at a given radius in the disk. Gravitationally unstable disks are expected to develop substructures like spiral arms (e.g., Durisen et al. 2007; Kratter & Lodato 2016). Numerical simulations by Hall et al. (2019) indicate that gravitational instability-induced spiral structures are only detectable when $M_{disk}/M_* \ge 0.25$.

To check these criteria for our sample, we collected stellar mass estimates for FUors from the literature. We used $<0.1 M_☉$ for OO Ser (Hodapp et al. 1996), $0.3 M_☉$ for L1551 IRS 5 (Osorio et al. 2003), $0.75 M_☉$ for V960 Mon (Kóspál et al. 2015), $2 M_☉$ for V899 Mon (Ninan et al. 2015), $0.1 M_☉$ for V346 Nor (Kóspál et al. 2017a), and $0.6 M_☉$ for FU Ori (Pérez et al. 2020). For objects without a published stellar mass, we assumed the same value as for FU Ori, $0.6 M_☉$. For all calculations, the RT disk masses were used. To calculate $c_s$ and $\Omega$, we used the RT temperature and density profiles. Table 8 shows whether a disk satisfies any of the four criteria mentioned above and the number of disks that are gravitationally unstable based on these metrics.

We plotted the stellar masses as a function of the disk mass in Figure 4. Here we complemented our sample with objects from the literature. The $M_{disk}/M_* \ge 0.1$ and $0.25$ criteria (for likely instability and the appearance of substructures, respectively) are indicated by dashed and dotted curves in this figure. Interestingly, while some FUors clearly satisfy the stricter criterion, we do not detect such structures in our targets.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Stellar mass as a function of disk mass (both in $M_☉$). Orange dots indicate FUors from this work, and green dots are objects from the literature. Light blue shading indicates the area where the disk is relatively massive compared to the host star ($M_{disk}/M_* \ge 0.1$) and thus most likely to be gravitationally unstable. The dotted curve marks the disk-to-stellar mass ratio above which gravitational instability-induced substructures are readily detectable. V899 Mon is actually outside of the plotted region at $M_* = 2.0 M_☉$.</figcaption></figure>

reason for this may be that the spatial resolution of our ALMA observations is insufficient to resolve the expected substructures in these relatively distant objects. Alternatively, the ongoing burst might have sent out an expansion wave through the disk and smeared out the preexisting spiral structure, as recently proposed by Vorobyov et al. (2020).

Figure 4 shows that 11 out of a sample of 17 FUors (65%) are in the gravitational instability region. Even if we only consider objects with published stellar masses, there are seven FUors out of 10 (70%) that may have gravitationally unstable disks, suggesting that gravitational instability may be an important factor in driving the outbursts. On the other hand, V899 Mon, FU Ori, and possibly AR 6A and HBC 687 seem to harbor gravitationally stable disks, and triggering their outbursts may need some other mechanisms.
## 6.4. Implications for the Ubiquity of the FUor Phenomenon

We discussed previously that FUor disks are, on average, more compact and massive than the disks around class I or II objects. This conclusion has important implications and raises interesting questions. One possible interpretation of these results is that there is a phase of the circumstellar disk evolution when the disk structure temporarily deviates from the typical one, forms a new structure in which it can produce outbursts, then relaxes back to a normal disk structure.

Alternatively, FUor eruptions may happen in normal disks, but the outburst itself (the increased accretion heating, enhanced wind activity, etc.) transforms the disk structure from normal into the special one we observe around eruptive stars. In these cases, it is a question what physical mechanism can make the disks smaller and their masses larger (especially in later evolutionary phases when replenishment from the envelope is no longer significant).

Finally, it may be possible that some disks are born with a smaller radius and larger mass than usual, and later they follow a particular evolutionary path that leads to episodic eruptions. Numerical simulations suggest that this may be the consequence of delayed infall from a slowly rotating envelope (Vorobyov et al. 2015). If FUor disks are born special and follow a special evolutionary path, we would need to accept that not all YSOs go through FUor eruptions.

While it is often assumed that all YSOs undergo FUor eruptions, our results raise the possibility that they instead form a special class of YSOs defined by their unusual disk structure and enhanced disk accretion. If so, then members of this special class are relatively rare among young stars. Contreras Peña et al. (2019) analyzed a large sample of 15,000 YSOs and identified only nine FUors, suggesting that at any given time, only 0.06% of young stars are in FUor outburst. A somewhat higher rate was reported for the population of YSOs in the Orion molecular clouds by Fischer et al. (2019), who found two outbursting sources out of 319 protostars (0.6%). The true incidence rate of FUors can be much higher if we take into account their duty cycle, i.e., the fraction of the life of a FUor-type object spent in outburst.

One way to determine the true fraction of eruptive stars among YSOs is to monitor them and discover new outbursts. All-sky monitoring surveys like Gaia and the Wide-field Infrared Survey Explorer gave a new impetus to this field of study (e.g., Hillenbrand et al. 2018, 2019, 2021; Szegedi-Elek et al. 2020). Another way would be to explore the accretion properties and structure of disks falling above the class I/II relations in Figure 2.

Less than 1% is a rather low rate for actual outbursting FUors, but these objects may still be important players in the star formation process. Considering that their accretion rate is ~1000 times higher than that of normal YSOs ($10^{-4}$ instead of $10^{-7}$; Hartmann & Kenyon 1996), the matter accreted in the young eruptive stars may be comparable to or even higher than the total mass accreted in the whole known YSO population.

## 7. Summary and Conclusions

We presented an ALMA millimeter continuum study of FUors, a class of low-mass pre-main-sequence stars undergoing accretion-related outbursts. Our sample of 10 objects contains a large part of all FUors visible from the ALMA site. We mapped the thermal dust continuum emission at 1.33 mm using the 7 m

array and two different configurations of the 12 m array. This way, we recover emission at the widest possible range of spatial scales, from the typical beam size of 0.′′15 to the largest recoverable scale of 5.′′9. We detected all of the targeted FUors and several neighboring sources as well. In most cases, the detected emission is bright, centrally peaked, and rather compact, and we attribute them to the circumstellar disks. To obtain the basic geometrical parameters and integrated millimeter fluxes of the sources, we first fitted 2D Gaussians to the images. As a next step, we fitted an analytical disk model of a power-law brightness distribution with an exponential taper in the visibility space using an MCMC method. Finally, we computed RT disk models with RADMC-3D and fitted them in the image space again with an MCMC method.

The geometrical parameters from the different models are generally in good agreement within the uncertainties. However, the RT modeling gives disk masses significantly larger than what is obtained from the measured millimeter fluxes assuming optically thin emission, suggesting that the FUor disks are optically thick at 1.33 mm. From the RT modeling, we obtained disk masses between 0.011 and $0.38 M_{\odot}$ and disk characteristic radii between 16 and 69 au. We complemented our sample with FUors from the literature and compared their masses and radii with comparison samples of normal, noneruptive class I and II objects modeled in a similar way with RT codes. The comparison revealed that FUor disks are typically a factor of 2.9–4.4 more massive and a factor of 1.5–4.7 smaller in size than the disks of class I/II objects. Using the criterion of $M_{\text{disk}}/M_{*} \ge 0.1$, about 65%–70% of the FUor disks may be gravitationally unstable.

It remains unclear whether all normal protostellar disks could develop properties that we observe for the FUor sample discussed here, which would imply that episodic accretion during the pre-main-sequence stellar evolution may be common. Our main finding that FUor disks are generally massive and compact may imply that they form a special class or evolutionary path.

There are several possible future directions that may provide better insights into the physics of young eruptive stars than what is possible at the moment. The RT modeling of larger samples of YSO disks will enable more robust comparison statistics than we could do here. Deeper and higher spatial resolution millimeter interferometric observations may better trace the disk–envelope connection and might reveal sub-structures in the disks related to gravitational instability. Observations at multiple bands, including longer wavelengths with the planned Band 1 and 2 receivers, will help constrain dust properties and optical depth effects. Targeting the appropriate molecular tracers will allow the study of possible outflow activity and interesting chemical processes related to the enhanced luminosity during the outburst. Thanks to several ongoing and planned monitoring surveys, new eruptive star candidates will be discovered whose follow-up may hold the key to understanding episodic accretion and its effects.

We thank the anonymous reviewer for useful comments that helped us to improve the paper. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 716155 (SACCRED). This paper makes use of the following ALMA data: ADS/JAO. ALMA#2016.1.00209.S. ALMA is a partnership of the ESO
(representing its member states), NSF (USA), and NINS (Japan), together with the NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by the ESO, AUI/NRAO, and NAOJ. On behalf of the SACCRED project, we are thankful for the usage of MTA Cloud (https://cloud.mta.hu/), which helped us significantly in achieving the results published in this paper. M.T. is supported by the Ministry of Science and Technology (MoST) of Taiwan (grant Nos. 106-2119-M-001-026-MY3, 109-2112-M-001-019). H.B.L. and M.T. are supported by MoST of Taiwan 108-2923-M-001-006-MY3 for the Taiwanese-Russian collaboration project. E.V. acknowledges support from the

Russian Fund for Fundamental Research, Russian–Taiwanese project 19-52-52011.

*Facility: ALMA.*

*Software:* CASA (McMullin et al. 2007), Galario (Tazzari et al. 2018), emcee (Foreman-Mackey et al. 2013), RADMC-3D (Dullemond et al. 2012), radmc3dPy.

## Appendix A
Visibility Space Fitting Results

Figures 5–14 show the observed clean maps, best-fit models convolved with the clean beam, residual maps, and observed and model visibilities as a function of uv distance.

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

**Figure 5.** Maps and visibility fitting of Haro 5a IRS. The first panel is the clean map with uniform weighting, the second panel is the best-fit model convolved with a 2D Gaussian representing the clean beam, and the third panel is the residual map. For the three panels in the image space, the contours are at −10σ, −8σ, −7σ, −5σ, −3σ, 3σ, 9σ, 31σ, 103σ, and 334σ with σ = 138 µJy. In the fourth panel, the black points show the observed visibility components plotted against the uv distance, and the blue line shows the best-fit model. The uv distances of the observed and modeled visibilities have been phase shifted to the best-fit ΔR.A. and Δdecl. and deprojected by *i* and PA.

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

**Figure 6.** Same as Figure 5 but for V899 Mon. The contours are at −5σ, −4σ, −3σ, 3σ, 6σ, 14σ, 30σ, and 64σ with σ = 23 µJy.

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

**Figure 7.** Same as Figure 5 but for AR 6A. The contours are at −5σ, −4σ, −3σ, 3σ, 5σ, 9σ, 17σ, and 30σ with σ = 46 µJy.
<figure><img src="image_10.png" /><figcaption>Figure 8. Same as Figure 5 but for AR 6B. The contours are at −7σ, −6σ, −5σ, −4σ, −3σ, 3σ, 7σ, 14σ, 32σ, and 70σ with σ = 46 µJy.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9. Same as Figure 5 but for V900 Mon. The contours are at −6σ, −5σ, −4σ, −3σ, 3σ, 9σ, 26σ, 78σ, and 232σ with σ = 30 µJy.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10. Same as Figure 5 but for V960 Mon. The contours are at −4σ, −3σ, 3σ, 5σ, 8σ, 14σ, and 23σ with σ = 28 µJy.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Same as Figure 5 but for Bran 76. The contours are at −5σ, −4σ, −3σ, 3σ, 8σ, 21σ, 55σ, and 146σ with σ = 27 µJy.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12. Same as Figure 5 but for V346 Nor. The contours are at −7σ, −6σ, −5σ, −4σ, −3σ, 3σ, 8σ, 23σ, 63σ, and 174σ with σ = 31 µJy.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13. Same as Figure 5 but for OO Ser. The contours are at −3σ, 3σ, 7σ, 17σ, 42σ, and 101σ with σ = 64 µJy.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14. Same as Figure 5 but for HBC 687. The contours are at −19σ, −15σ, −11σ, −7σ, −3σ, 3σ, 11σ, 38σ, 135σ, and 478σ with σ = 42 µJy.</figcaption></figure>
Appendix B
RADMC Modeling Results

The RT modeling used an image plane fitting approach. The
observed dirty maps, best-fit models convolved with the dirty

beam, and data minus model residuals for each system are
shown in Figures 15–25.

<figure><img src="image_17.png" /><figcaption>Figure 15. The RADMC-3D model fitting of L1551 IRS 5 N and S. The left panel is the dirty map, the middle panel is the best-fit model convolved with the dirty beam, and the right panel is the residual map.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16. Same as Figure 15 but for Haro 5A IRS.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17. Same as Figure 15 but for V899 Mon.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Same as Figure 15 but for AR 6A.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Same as Figure 15 but for AR 6B.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Same as Figure 15 but for V900 Mon.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Same as Figure 15 but for V960 Mon.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Same as Figure 15 but for Bran 76.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 23.** Same as Figure 15 but for V346 Nor.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 24.** Same as Figure 15 but for OO Ser.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 25.** Same as Figure 15 but for HBC 687.</figcaption></figure>
Appendix C

Posterior Distribution Functions for the RT Modeling

Figures 26–28 show examples of the posterior distributions of the fitted parameters in the RT modeling.

<figure><img src="image_28.png" /><figcaption>Figure 26. Posterior distributions for L1551 IRS 5 N.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure 27. Posterior distributions for Bran 76.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 28. Posterior distributions for V346 Nor.</figcaption></figure>
Appendix D
Comparison of Geometrical Disk Parameters

Figure 29 shows a comparison of the geometrical disk parameters (radius, inclination, and PA) we obtained from the

different types of modeling. Here 2DG stands for 2D Gaussian fitting in the image space, RT stands for RT modeling in the image space with RADMC-3D (Section 4), and Vis stands for disk model fitting in the visibility space (Section 3.1).

<figure><img src="image_31.png" /><figcaption>(a) Haro 5a IRS</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(c) AR 6A</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(e) V900 Mon</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(g) Bran 76</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(i) OO Ser</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(b) V899 Mon</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(d) AR 6B</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(f) V960 Mon</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(h) V346 Nor</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(j) HBC 687</figcaption></figure>

**Figure 29.** Comparison of the best-fit geometrical parameters (radius, inclination, and PA) of the disks with three different methods. The blue circles show the 2D Gaussian fitting in the image space, the red squares are the results from the RT modeling in the image space with RADMC-3D, and the diamonds show the disk model fitting in the visibility space. For the latter, purple indicates $R_{\text{eff}}$, while green indicates $R_{\text{disk}}$. Uncertainties correspond to $1\sigma$.
ORCID iDs

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

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

J. A. White https://orcid.org/0000-0001-8445-0444

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

L. Chen https://orcid.org/0000-0003-2835-1729

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

R. Dong https://orcid.org/0000-0001-9290-7846

M. M. Dunham https://orcid.org/0000-0003-0749-9505

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

J. D. Green https://orcid.org/0000-0003-1665-5709

J. Hashimoto https://orcid.org/0000-0002-3053-3575

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

M. Hogerheijde https://orcid.org/0000-0001-5217-537X

T. Kudo https://orcid.org/0000-0002-9294-1793

H. B. Liu https://orcid.org/0000-0003-2300-2626

M. Takami https://orcid.org/0000-0001-9248-7546

E. I. Vorobyov https://orcid.org/0000-0002-6045-0359

References

Andrews, S. M., Wilner, D., Hughes, A., Qi, C., & Dullemond, C. 2009, ApJ, 700, 1502

Andrews, S. M., Wilner, D., Hughes, A., Qi, C., & Dullemond, C. 2010, ApJ, 723, 1241

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241

Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705

Aspin, C., & Reipurth, B. 2003, AJ, 126, 2936

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

Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Boley, A. C., & Durisen, R. H. 2008, ApJ, 685, 1193

Booth, M., Jordán, A., Casassus, S., et al. 2016, MNRAS, 460, L10

Caratti o Garatti, A., Garcia Lopez, R., Ray, T. P., et al. 2015, ApJL, 806, L4

Caratti o Garatti, A., Steclum, B., Garcia Lopez, R., et al. 2017, NatPh, 13, 276

Chini, R., Reipurth, B., Ward-Thompson, D., et al. 1997, ApJL, 474, L135

Chou, T.-L., Takakuwa, S., Yen, H.-W., Ohashi, N., & Ho, P. T. P. 2014, ApJ, 796, 70

Cieza, L. A., Casassus, S., Tobin, J., et al. 2016, Natur, 535, 258

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

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

Contreras Peña, C., Naylor, T., & Morrell, S. 2019, MNRAS, 486, 4590

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

Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503

Draine, B. T. 2006, ApJ, 636, 1114

Dullemond, C., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, http://ascl.net/1202.015

Dunham, M. M., Arce, H. G., Allen, L. E., et al. 2013, AJ, 145, 94

Dunham, M. M., Stutz, A. M., Allen, L. E., et al. 2014a, Protostars and Planets VI (Tucson, AZ: Univ. Arizona Press), 195

Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52

Dunham, M. M., Vorobyov, E. I., & Arce, H. G. 2014b, MNRAS, 444, 887

Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 607

Eisloffel, J., Hessman, F. V., & Mundt, R. 1990, A&A, 232, 70

Fischer, W. J., Megeath, S. T., Ali, B., et al. 2010, A&A, 518, L122

Fischer, W. J., Safron, E., & Megeath, S. T. 2019, ApJ, 872, 183

Forbrich, J., Tappe, A., Robitaille, T., et al. 2010, ApJ, 716, 1453

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fridlund, C. V. M., Liseau, R., Djupvik, A. A., et al. 2005, A&A, 436, 983

Graham, J. A., & Frogel, J. A. 1985, ApJ, 289, 331

Gramajo, L. V., Rodón, J. A., & Gómez, M. 2014, AJ, 147, 140

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

Großschedl, J. E., Alves, J., Meingast, S., & Hasenberger, B. 2018, arXiv:1812.08024

Hackstein, M., Haas, M., Kóspál, Á., et al. 2015, A&A, 582, L12

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

Hall, C., Dong, R., Rice, K., et al. 2019, ApJ, 871, 228

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

Henning, T., Burkert, A., Launhardt, R., Leinert, C., & Steclum, B. 1998, A&A, 336, 565

Herczeg, G. J., Kuhn, M. A., Zhou, X., et al. 2019, ApJ, 878, L11

Hillenbrand, L. A., Contreras Peña, C., Morrell, S., et al. 2018, ApJ, 869, 146

Hillenbrand, L. A., De, K., Hankins, M., et al. 2021, AJ, 161, 220

Hillenbrand, L. A., Reipurth, B., Connelley, M., Cutri, R. M., & Isaacson, H. 2019, AJ, 158, 240

Hodapp, K. W. 1999, AJ, 118, 1338

Hodapp, K. W., Chini, R., Watermann, R., & Lemke, R. 2012, ApJ, 744, 56

Hodapp, K.-W., Hora, J. L., Rayner, J. T., Pickles, A. J., & Ladd, E. F. 1996, ApJ, 468, 861

Kadam, K., Vorobyov, E., Regály, Z., Kóspál, Á., & Ábrahám, P. 2019, ApJ, 882, 96

Kainulainen, J., Stutz, A. M., Stanke, T., et al. 2017, A&A, 600, A141

Kenyon, S. J., Dobrzycka, D., & Hartmann, L. 1994, AJ, 108, 1872

Kóspál, Á., Ábrahám, P., Acosta-Pulido, J. A., et al. 2016, A&A, 596, A52

Kóspál, Á., Ábrahám, P., Apai, D., et al. 2008, MNRAS, 383, 1015

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

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

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

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

Kóspál, Á., Ábrahám, P., Prusti, T., et al. 2007, A&A, 470, 211

Kóspál, Á., Ábrahám, P., Westhues, C., & Haas, M. 2017c, A&A, 597, L10

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

Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

Kraus, S., Caratti o Garatti, A., Garcia-Lopez, R., et al. 2016, MNRAS, 462, L61

Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642

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

Liseau, R., Fridlund, C. V. M., & Larsson, B. 2005, ApJ, 619, 959

Liu, H. B. 2019, ApJL, 877, L22

Liu, H. B., Dunham, M. M., Pascucci, I., et al. 2018, A&A, 612, A54

Lodato, G., & Clarke, C. J. 2004, MNRAS, 353, 841

Maíz Apellániz, J. 2019, A&A, 630, A119

Manara, C. F., Tazzari, M., Long, F., et al. 2019, A&A, 628, A95

Martin, R. G., & Lubow, S. H. 2011, ApJL, 740, L6

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Molyarova, T., Akimkin, V., Semenov, D., et al. 2018, ApJ, 866, 46

Nielbock, M., Chini, R., & Müller, S. A. H. 2003, A&A, 408, 245

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

Ninan, J. P., Ojha, D. K., & Philip, N. S. 2016, ApJ, 825, 65

Osorio, M., D'Alessio, P., Muzerolle, J., Calvet, N., & Hartmann, L. 2003, ApJ, 586, L148

Pavlyuchenkov, Y., Akimkin, V., Wiebe, D., & Vorobyov, E. 2019, MNRAS, 486, 3907

Pérez, L. M., Lamb, J. W., Woody, D. P., et al. 2010, ApJ, 724, 493

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

Pyo, T.-S., Hayashi, M., Kobayashi, N., et al. 2005, ApJ, 618, 817

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

R Core Team 2021, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, https://www.R-project.org/

Reipurth, B. 1985, A&A, 143, 435

Reipurth, B., & Aspin, C. 1997, AJ, 114, 2700

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

Reipurth, B., Hartmann, L., Kenyon, S. J., Smette, A., & Bouchet, P. 2002, AJ, 124, 2194

Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328
Rodríguez, L. F., D’Alessio, P., Wilner, D. J., et al. 1998, Natur, 395, 355

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

Sandell, G., & Weintraub, D. A. 2001, ApJS, 134, 115

Sheehan, P. D., & Eisner, J. A. 2014, ApJ, 791, 19

Sheehan, P. D., & Eisner, J. A. 2017a, ApJ, 851, 45

Sheehan, P. D., & Eisner, J. A. 2017b, ApJL, 840, L12

Staude, H. J., & Neckel, T. 1992, ApJ, 400, 556

Szegedi-Elek, E., Ábrahám, P., Wyrzykowski, Ł., et al. 2020, ApJ, 899, 130

Takahashi, S., Ho, P. T. P., Teixeira, P. S., Zapata, L. A., & Su, Y.-N. 2013, ApJ, 763, 57

Takakuwa, S., Ohashi, N., Ho, P. T. P., et al. 2004, ApJL, 616, L15

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

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

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43

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

Toon, O. B., & Ackerman, T. 1981, ApOpt, 20, 3657

Vorobyov, E. I., & Basu, S. 2005, ApJL, 633, L137

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

Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896

Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115

Vorobyov, E. I., Elbakyan, V. G., Takami, M., & Liu, H. B. 2020, A&A, 643, A13

Vorobyov, E. I., Lin, D. N. C., & Guedel, M. 2015, A&A, 573, A5

Wang, H., Apai, D., Henning, T., & Pascucci, I. 2004, ApJL, 601, L83

White, J. A., Boley, A., Dent, W., Ford, E., & Corder, S. 2016, MNRAS, 466, 4201

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

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

Xi, Y., Dunham, M. M., Arce, H. G., et al. 2020, RNAAS, 4, 155

Yoo, H., Lee, J.-E., Mairs, S., et al. 2017, ApJ, 849, 69

Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009, ApJ, 701, 620

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18

Zubko, V., Kreiowski, J., & Wegner, W. 1996, MNRAS, 283, 577

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