# The Circumstellar Disk and Asymmetric Outflow of the EX Lup Outburst System

A. S. Hales, S. Pérez, M. Saito, C. Pinte, L. B. G. Knee, I. de Gregorio-Monsalvo, B. Dent, C. López, A. Plunkett, P. Cortés, et al.

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

## To cite this version:

A. S. Hales, S. Pérez, M. Saito, C. Pinte, L. B. G. Knee, et al., The Circumstellar Disk and Asymmetric Outflow of the EX Lup Outburst System. The Astrophysical Journal, 2018, 859, 10.3847/1538-4357/aac018. insu-03693555

HAL Id: insu-03693555

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

Submitted on 19 Aug 2025

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

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

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

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

# The Circumstellar Disk and Asymmetric Outflow of the EX Lup Outburst System

A. S. Hales<sup>1,2</sup>, S. Pérez<sup>3</sup>, M. Saito<sup>4,5</sup>, C. Pinte<sup>6</sup>, L. B. G. Knee<sup>7</sup>, I. de Gregorio-Monsalvo<sup>1,8</sup>, B. Dent<sup>1,8</sup>, C. López<sup>1</sup>,
A. Plunkett<sup>8</sup>, P. Cortés<sup>1,2</sup>, S. Corder<sup>1,2</sup>, and L. Cieza<sup>9,10</sup>

Received 2018 February 10; revised 2018 April 19; accepted 2018 April 23; published 2018 May 30

## Abstract

We present Atacama Large Millimeter/submillimeter Array (ALMA) observations at 0.′′3 resolution of EX Lup, the prototype of the EXor class of outbursting pre-main-sequence stars. The circumstellar disk of EX Lup is resolved for the first time in 1.3 mm continuum emission and in the *J* = 2–1 spectral line of three isotopologues of CO. At the spatial resolution and sensitivity achieved, the compact dust continuum disk shows no indications of clumps, fragments, or asymmetries above the 5σ level. Radiative transfer modeling constrains the characteristic radius of the dust disk to 23 au and the total dust mass to $1.0 \times 10^{-4} M_\odot$ (33 $M_\oplus$), similar to other EXor sources. The $^{13}$CO and C$^{18}$O line emissions trace the disk rotation and are used to constrain the disk geometry, kinematics, and a total gas disk mass of $5.1 \times 10^{-4} M_\odot$. The $^{12}$CO emission extends out to a radius of 200 au and is asymmetric, with one side deviating from Keplerian rotation. We detect blueshifted, $^{12}$CO arc-like emission located 0.′′8 to the northwest and spatially disconnected from the disk emission. We interpret this extended structure as the brightened walls of a cavity excavated by an outflow, which are more commonly seen in FUor sources. Such outflows have also been seen in the borderline FU/EXor object V1647 Ori, but not toward EXor objects. Our detection provides evidence that the outflow phenomenon persists into the EXor phase, suggesting that FUor and EXor objects are a continuous population in which outflow activity declines with age, with transitional objects such as EX Lup and V1647 Ori.

**Key words:** circumstellar matter – protoplanetary disks – stars: formation – stars: individual (EX Lup) – stars: pre-main sequence

## 1. Introduction

Most young stars are now thought to undergo phases of episodic accretion during which their accretion rates can increase by several orders of magnitude (Hartmann & Kenyon 1996). Although short-lived compared to the total life of the accretion disks, these bursts are believed to play a key role in the building up of the final stellar mass (Hartmann 2008; Hartmann et al. 2016). These episodic accretion events may alleviate the well-known luminosity problem, in which low-mass stars appear underluminous compared to the predictions of steady-state accretion models (Kenyon et al. 1990; Evans et al. 2009). The enhancement in accretion luminosity can affect the dust and gas properties, and the variable snowlines can influence the disk chemistry, dust grain properties, recondensation of ices, and possibly planet formation (Ábrahám et al. 2009; Banzatti et al. 2015; Cieza et al. 2016; Hubbard 2017; Rab et al. 2017).

Outburst sources have been historically divided into two classes, FUors (named after their prototype FU Ori) and EXors (named after their prototype EX Lup). The class distinction is observational: FUors have large (>5 mag), long-lived (years to decades) bursts (Herbig 1966), whereas EXors have moderate (2–4 mag), shorter (days to months) episodes of high accretion. Recent observations have questioned the need for two separate classifications, leading to the suggestion that both types of

outburst may have a common origin (Audard et al. 2014, and references therein).

Different possible triggering mechanisms have been proposed for the episodic outburst: (1) thermal instability (Bell & Lin 1994), (2) coupling of the gravitational instability (GI) and magnetorotational (MRI) instability (Armitage et al. 2001; Zhu et al. 2009), (3) disk fragmentation (Vorobyov & Basu 2015), and (4) tidal interactions between the disk and either planets (Lodato & Clarke 2004) or stellar companions (Bonnell & Bastien 1992). D'Angelo & Spruit (2010) showed that disk material near corotation radius can become trapped, which could also produce variable accretion rates. Recently, Armitage (2016) proposed that EXor outbursts could be produced by the interaction between stellar dynamo cycles and the inner regions of the disk.

Some proposed mechanisms for FUor and EXor outbursts (e.g., GI and disk fragmentation) predict the presence of distinct morphological features that might be detectable in the dust continuum of massive disks, yet to date none of these signatures have been found in high-resolution interferometric images of FUor/EXor disks (Hales et al. 2015; Cieza et al. 2016; Kóspál et al. 2017b; Ruíz-Rodríguez et al. 2017; Zurlo et al. 2017; Principe et al. 2018). Even in V883 Ori, the most massive FUor disk known and thus the best candidate for a GI to occur, ALMA observations at 0.′′03 (12 au) resolution were unable to detect the predicted signatures of disk instability

<footnote n="1"><sup>1</sup> Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura 7630355, Santiago, Chile; antonio.hales@alma.cl</footnote>

<footnote n="2"><sup>2</sup> National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA</footnote>

<footnote n="3"><sup>3</sup> Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago 8330015, Chile</footnote>

<footnote n="4"><sup>4</sup> National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="5"><sup>5</sup> Department of Astronomy, The Graduate University for Advanced Studies (SOKENDAI), Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="6"><sup>6</sup> Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="7"><sup>7</sup> NRC Herzberg Astronomy and Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="8"><sup>8</sup> European Southern Observatory, Avenida Alonso de Córdova 3107, Vitacura 7630355, Santiago, Chile</footnote>

<footnote n="9"><sup>9</sup> Nucleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile</footnote>

<footnote n="10"><sup>10</sup> Millenium Nucleus Protoplanetary Discs in ALMA Early Science, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile</footnote>
or fragmentation (Cieza et al. 2016). Observations of FUor/EXor sources show that the continuum disks are compact, with characteristic radii $R_c < 20–40$ au (Cieza et al. 2017), and most likely optically thick in their inner regions—making the total dust mass difficult to estimate (Cieza et al. 2016; Liu et al. 2017).

Complementary to dust continuum observations are spectroscopic studies of the gas content of the disk, which provides information on both the structure and kinematics of the system. For example, deviations from Keplerian motion may be indicative of disk instabilities or enhanced accretion (Ilee et al. 2011; Evans et al. 2015). Several physical mechanisms can produce non-Keplerian radial velocity profiles, such as (1) radial inflows, (2) a physical warp in the inner disk, (3) circumbinary disk interactions, and (4) a perturbation from a stellar or planetary companion (Rosenfeld et al. 2014; Casassus et al. 2015; Takakuwa et al. 2017). All these processes can be linked to enhanced stellar accretion; therefore, detecting such signatures in spectral line images can provide important clues to the physical processes that drive episodic accretion in young stars. Furthermore, observation of mass inflow onto the star may help constrain estimates of the quantity of mass being transferred during episodic accretion events.

Spectral line observations show that all studied FUor sources have active circumstellar environments characterized by strong outflow activity (Kóspál et al. 2017b; Ruíz-Rodríguez et al. 2017; Zurlo et al. 2017). On the other hand, EXor sources do not show detectable outflows, with the possible exception of V1647 Ori—a system with an ambiguous FUor/EXor classification (Principe et al. 2018). There is also observational evidence that most FUor sources are surrounded by large envelopes that are still transferring material onto the disk (e.g., Fehér et al. 2017; Kóspál et al. 2017a, and references therein). It has been proposed that the imbalance between mass transferred from the envelope to the disk and the accretion of disk mass into the star may trigger the disk instability responsible for the outburst (Bell & Lin 1994; Goto et al. 2011; Kóspál et al. 2017b). Studying the large-scale structure (outer disk, envelope, outflows) and the substructures (disk inner regions, clumps) of eruptive sources is therefore crucial for understanding the nature of the FUor/EXor outbursts.

EX Lup is a young M0 star located in the Lupus 3 cloud ($d = 155$ pc; Gras-Velázquez & Ray 2005; Comerón 2008; Lombardi et al. 2008) and is the prototype of the EXor class of young eruptive stars (Herbig 1950, 1989, 2008). EX Lup has undergone several outbursts, with the two strongest observed in 1955 (Herbig 1977) and 2008 (Jones 2008). These outbursts are usually ascribed to episodes of mass infall up to three orders of magnitude higher than typical quiescent accretion rates (Audard et al. 2014). The inner disk of EX Lup has been studied thoroughly using spectroscopic techniques, which has shown strong variations in the physical conditions of the gaseous and silicate components in the inner $<0.4$ au of the disk during and post-outburst (Goto et al. 2011; Kóspál et al. 2011; Juhász et al. 2012; Sicilia-Aguilar et al. 2012, 2015). Studies of ionized metal lines have been used to probe accretion processes in EX Lup. These line profiles appear to arise in a hot, asymmetric accreting region located in the inner disk regions 0.1–0.2 au from the star (Sicilia-Aguilar et al. 2012), providing evidence for a strong inner disk wind that increases during the periods of high accretion. Observations of the vibrationally excited mid-infrared CO lines also indicate a strong inner disk asymmetry, attributed to the

presence of a disk hot spot, as well as the presence of a CO disk wind (Goto et al. 2011).

However, very little is known about the outer disk of EX Lup since it has never been spatially resolved at any wavelength. There is evidence for a $0.025 M_\odot$ and 150 au disk from spectral energy distribution (SED) modeling (Sipos et al. 2009). The SED of EX Lup resembles that of a Class II object, suggesting that the circumstellar environment consists of a disk with possibly a thin envelope. In that sense, despite its outburst activity and young age (0.5 Myr; Frasca et al. 2017), EX Lup appears to be relatively evolved among eruptive young stellar objects. Recently, Kóspál et al. (2016b) detected molecular emission from EX Lup in the $^{12}\text{CO}$ (3–2), $^{12}\text{CO}$ (4–3), and $^{13}\text{CO}$ (3–2) lines. They detected single-peaked line profiles that they modeled as a Keplerian disk. The total disk gas mass inferred from the $^{13}\text{CO}$ data ($2.3 \times 10^{-4} M_\odot$) is significantly lower than the one derived from the continuum data (assuming a canonical gas-to-dust ratio of 100), indicative of CO depletion. Despite these efforts, the basic disk parameters of the dust and gas disk in EX Lup remain poorly constrained in the absence of spatially resolved images.

In this work we present ALMA observations that resolve the disk around EX Lup in both dust and gas. Section 2 describes the details of the ALMA observations and data reduction. In Section 3 we present the results and an analysis of the data. A comparison of these results with radiative transfer modeling is discussed in Section 4, with a full discussion of our results in Section 5. Our conclusions are summarized in Section 6.

## 2. Observations

ALMA observations of EX Lup using the facility Band 6 ($\sim 230$ GHz) receivers were acquired on 2016 July 25 (during ALMA Cycle 3) utilizing interferometric observations on 40 antennas in the main 12 m antenna array, providing baselines ranging from 15.1 m to 1.12 km. For Band 6 this antenna configuration yields an angular resolution of $\sim 0.0''3$ ($\sim 45$ au) and a Maximum Recoverable Scale (MRS; see Schieven 2017) of $\sim 3''2$ ($\sim 500$ au). During the observations, the precipitable water vapor column in the atmosphere was stable at 0.8–0.9 mm with clear sky conditions, resulting in a median system temperature of 70 K.

Four spectral windows were positioned to target the $^{12}\text{CO}$ (2–1), $^{13}\text{CO}$ (2–1), $\text{C}^{18}\text{O}$ (2–1), and the vibrationally excited $^{13}\text{CO}$ (2–1, $v=1$) transitions of carbon monoxide (rest frequencies of 230.538, 220.399, 219.560, and 218.437 GHz, respectively). The ALMA correlator was configured in Frequency Division Mode (FDM) to provide spectral resolutions of 61.035 kHz for $^{12}\text{CO}$ (2–1) and 122.070 kHz for the other transitions, corresponding to velocity resolutions of 0.092 and 0.192 km s$^{-1}$, respectively. A fifth spectral window configured in Time Division Mode (TDM) was positioned in a region devoid of line emission for detecting the continuum dust emission (centered at 232.477 GHz/1.29 mm). The effective bandwidth for each spectral window was 117.187 MHz for $^{12}\text{CO}$ (2–1), $^{13}\text{CO}$ (2–1), and $\text{C}^{18}\text{O}$ (2–1), 234.375 MHz for $^{13}\text{CO}$ (2–1, $v=1$), and 1.875 GHz for the continuum.

The primary flux calibrator was the quasar J1427–4206, which was also used as the bandpass calibrator. The flux density of the amplitude calibrator was obtained as a result of interpolation between Band 3 and Band 7 monitoring observations, both acquired 2 days before the target observation (the monitor observations used Callisto as the main flux calibrator).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left panel: ALMA 1.29 mm continuum image of the EX Lup disk. Middle panel: CLEAN image of the best-fit model to the EX Lup continuum data (as described in Section 4.1) after filtering through the observed visibilities using CASA task FT (version 5.0). Right panel: residual image obtained by subtracting the model from the data. The color maps in the left and middle panels are normalized to an asinh scale. The residual image on the right has maximum and minimum values equal to 0.24 and -0.18 mJy beam⁻¹, respectively. Contours show residual levels at -3, -2, -1, 0, 1, 2, and 3 times the rms of 38 µJy beam⁻¹. The residuals suggest the presence of structure not accounted for by our axisymmetric model, at the 5σ level.</figcaption></figure>

The Band 3 to Band 7 spectral index used for the interpolation was -0.54 (uncertainty ±0.05). The quasar J1610-3958 was observed for phase calibration. Observations of the phase calibrator, located 1.4 from the target, were alternated with the science target every 6 minutes to calibrate the time-dependent variation of the complex gains. The total time spent on-source was 28.3 minutes. A secondary phase calibrator, J1604-4228, was observed regularly as a check source. The purpose of the check source was to assess the quality of the phase transfer (i.e., to estimate the coherence of the calibrated data during the Quality Assurance data processing).

All data were calibrated using the ALMA Science Pipeline (version r36660 of Pipeline-Cycle3-R4-B release) in *Common Astronomy Software Applications*<fnref n="11" />; McMullin et al. 2007) by the East Asia ALMA ARC Staff. The pipeline uses CASA tasks wherever possible to perform the data reduction and calibration in a standard fashion, which included offline Water Vapor Radiometer (WVR) calibration, system temperature correction, and bandpass, phase, and amplitude calibrations. Online flagging and nominal flagging, such as for shadowed antennas and band edges, were applied for calibration. In addition, some scans of the phase calibrator for six antennas (DA55, DA62, DV16, DV19, DV21, DV25) and the adjacent target scans were flagged owing to large phase scatter and/or discontinuous phase of the calibrator.

Imaging of the continuum and molecular emission lines was performed using the CLEAN task in CASA. The 1.875 GHz wide, line-free, spectral window centered at 232.477 GHz was imaged using Briggs weighting with a robust parameter of 0.5, resulting in a synthesized beam size of 0.32 × 0.26 at position angle PA = 61.8 in the final continuum image. A single iteration of phase-only self-calibration was performed to improve the coherence in the final image. The signal-to-noise ratio (S/N) improved by a factor of 1.8 after applying self-calibration. A second continuum image was produced at 218.474 GHz (1.37 mm) using the 234.375 MHz wide ¹³CO (2-1, ν = 1) spectral window using CLEAN parameters similar to those for the 1.29 mm image.

Continuum subtraction in the visibility domain was performed prior to imaging of each molecular line. Self-calibration tables

derived from the continuum were applied to the FDM spectral windows before imaging the CO lines. During CLEANing, the spectral axis was binned to 0.2 km s⁻¹ per channel for ¹²CO (2-1), to 0.3 km s⁻¹ for ¹³CO (2-1), and to 0.8 km s⁻¹ for C¹⁸O (2-1) and ¹³CO (2-1, ν = 1), to increase the S/N. CLEAN was run with natural weighting to maximize sensitivity. An rms noise level of 4.5 mJy beam⁻¹ per 0.2 km s⁻¹ channel was reached in the ¹²CO (2-1) line, with a 0.39 × 0.29 synthesized beam (PA = 62.7). The rms noise per channel in the ¹³CO (2-1) line was found to range between 3 and 4 mJy beam⁻¹, with a mean value of 3.6 mJy beam⁻¹. For the C¹⁸O (2-1) and ¹³CO (2-1, ν = 1) cubes, the rms ranged between 2 and 3 mJy beam⁻¹, with mean values of 2.2 and 2.0 mJy beam⁻¹ respectively. A second ¹²CO (2-1) line image was created using Briggs weighting with a robust parameter of 0.5, which provided a resolution of 0.32 × 0.26 (PA = 61.8).

# 3. Results and Analysis

## 3.1. Dust Continuum Emission

The disk around EX Lup is detected and resolved in the dust continuum. Figure 1 shows the 1.29 mm compact source, which is detected at an S/N of 226 (peak flux 8.8 mJy beam⁻¹ with noise rms 0.038 mJy beam⁻¹). The peak of the continuum emission coincides with the near-infrared position of EX Lup (as listed in the Two Micron All Sky Survey catalog, after correcting for the proper motion; Cutri et al. 2003). The parameters of the compact continuum emission can be determined from elliptical Gaussian fitting in the image domain using the CASA task IMFIT. The parameters of the elliptical Gaussian component determined using IMFIT are detailed in Table 1. A 0.45 × 0.37 (FWHM) disk at PA = 63.1 ± 1.2 is resolved, with a total 1.29 mm flux of 17.37 ± 0.15 mJy (uncertainties are derived from the IMFIT task). An estimate of the deconvolved disk size is 0.31 × 0.26 at PA = 64.6 ± 3.3. There is no evidence of residual emission above 5σ after subtracting the best-fit elliptical Gaussian from the continuum image.

Assuming that the deprojected EX Lup disk is circular, the major/minor-axis ratio implies a disk inclination angle of 32.4 ± 0.9. This inclination is in reasonable agreement with the value of 40°-50° derived by Goto et al. (2011) and with the

<footnote n="11">¹¹ http://casa.nrao.edu/</footnote>
<table><caption>**Table 1**
Continuum and Line Imaging Parameters</caption><thead><tr><th>Source Properties</th><th>1.29 mm</th><th>1.37 mm</th><th><sup>12</sup>CO</th><th><sup>13</sup>CO</th><th>C<sup>18</sup>O</th></tr></thead><tbody><tr><td>R.A.<sup>a</sup></td><td>16:03:05.479</td><td>0.000</td><td>-0.015</td><td>-0.013</td><td>-0.014</td></tr><tr><td>Decl.<sup>a</sup></td><td>-40:18:25.789</td><td>0.000</td><td>+0.010</td><td>+0.006</td><td>+0.019</td></tr><tr><td>Major-axis FWHM (arcsec)<sup>b</sup></td><td>0.45 ± 0.003</td><td>0.45 ± 0.007</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Minor-axis FWHM (arcsec)<sup>b</sup></td><td>0.37 ± 0.002</td><td>0.38 ± 0.005</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Position angle (deg)<sup>c</sup></td><td>63.1 ± 1.2</td><td>62.1 ± 3.0</td><td>80 ± 2</td><td>76 ± 3</td><td>71 ± 4</td></tr><tr><td>Peak intensity<sup>d</sup></td><td>8.8</td><td>8.2</td><td>605.9</td><td>94.7</td><td>34.4</td></tr><tr><td>Integrated flux<sup>e</sup></td><td>17.37 ± 0.15</td><td>15.16 ± 0.32</td><td>5011 ± 510</td><td>341 ± 37</td><td>77 ± 10</td></tr><tr><td>Line centroid<sup>f</sup> (km s<sup>-1</sup>)</td><td>...</td><td>...</td><td>3.7 ± 0.05</td><td>3.9 ± 0.2</td><td>5.0 ± 0.6</td></tr><tr><td colspan="6"><strong>Beam Properties and Image rms</strong></td></tr><tr><td>Major axis (arcsec)</td><td>0.32</td><td>0.34</td><td>0.39</td><td>0.41</td><td>0.41</td></tr><tr><td>Minor axis (arcsec)</td><td>0.26</td><td>0.28</td><td>0.29</td><td>0.31</td><td>0.32</td></tr><tr><td>Position angle (deg)</td><td>61.8</td><td>61.43</td><td>62.7</td><td>63.9</td><td>63.4</td></tr><tr><td>rms (mJy beam<sup>-1</sup>)</td><td>0.038</td><td>0.1</td><td>4.5</td><td>3.6</td><td>2.2</td></tr><tr><td>rms moment 0 (mJy beam<sup>-1</sup> km s<sup>-1</sup>)</td><td>...</td><td>...</td><td>8.9</td><td>7.3</td><td>5.3</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Centroid coordinates for 1.37 mm, <sup>12</sup>CO, <sup>13</sup>CO, and C<sup>18</sup>O relative to the 1.29 mm centroid, in units of arcseconds.

<sup>b</sup> Convolved with the beam.

<sup>c</sup> PA for continuum images derived from IMFIT. PA for line images derived from spectro-astrometry.

<sup>d</sup> Units are mJy beam<sup>-1</sup> for continuum and mJy beam<sup>-1</sup> km s<sup>-1</sup> for moment 0 images.

<sup>e</sup> Units are mJy for continuum and mJy km s<sup>-1</sup> for line images. Integrated fluxes for spectral line images are obtained by integrating over apertures of 1"/7, 0"/75, and 0"/55 for <sup>12</sup>CO, <sup>13</sup>CO, and C<sup>18</sup>O, respectively.

<sup>f</sup> Derived from Gaussian fits to the line profiles.

40° inclination determined by Kóspál et al. (2014). It should be noted that CASA documentation indicates that current limitations in IMFIT generally result in it underestimating the uncertainty in output axis size and PA fit values.

The total flux measured in the 234.375 MHz wide continuum window at 1.37 mm is 15.16 ± 0.32 mJy (the rms noise was 0.1 mJy beam<sup>-1</sup>). This results in a spectral index ($F_i \propto \nu^{\alpha}$) between 1.30 and 1.37 mm of 2.19 ± 0.47. Since the Band 3 and Band 7 observations of the flux calibrator were obtained simultaneously, and only 2 days apart from the observations of EX Lup, the effect of uncertainties in the spectral slope of the calibrator on the in-band spectral index of EX Lup is minimal. The derived spectral index is in agreement with previous estimates of the 1–3 mm slope of 2.6 ± 0.4 made with ATCA (Ubach et al. 2012). Assuming that all the flux detected at 0.87 mm with APEX/LABOCA comes from the disk (41 ± 10 mJy; Juhász et al. 2012), the derived 0.87/1.3 mm spectral index is 2.18, almost identical to the 1.3 mm spectral index we measure in-band.

## 3.2. Spectral Line Emission

Rotational line emission ($J = 2–1$) from all three main isotopologues of CO was detected. The vibrationally excited <sup>13</sup>CO (2–1, $\nu = 1$) line was not detected and provides a $1\sigma$ upper limit of 4.7 mJy km s<sup>-1</sup> (assuming a line width of 9 km s<sup>-1</sup>). Figure 2 shows the integrated emission (moment 0) maps for each detected transition with continuum contours overlaid. The <sup>12</sup>CO integrated line emission is 5011 ± 510 mJy km s<sup>-1</sup> (integrated over a circular aperture of 1"/7 centered at the stellar position). The <sup>13</sup>CO and C<sup>18</sup>O total integrated emissions are 341 ± 37 mJy km s<sup>-1</sup> and 77 ± 10 mJy km s<sup>-1</sup> (integrated over circular apertures of 0"/75 and 0"/55, respectively). The quoted error includes both the statistical noise and absolute flux uncertainty, added in quadrature.

The <sup>12</sup>CO-emitting region is notably more extended than the dust continuum, <sup>13</sup>CO-emitting, and C<sup>18</sup>O-emitting regions. The <sup>12</sup>CO can be traced above $3\sigma$ out to ~1"/3 (~200 au) from the stellar position, while <sup>13</sup>CO and C<sup>18</sup>O extend out to ~0"/5 and ~0"/3, respectively. We note that the full extension of the observed <sup>12</sup>CO emission (~3"/4 at $1\sigma$) is close to the MRS allowed by the observations.

Figures 3 and 4 show the intensity-weighted velocity fields (moment 1) and the integrated spectrum for the three CO isotopologues. The moment 1 images are indicative of rotation around the continuum disk, although only <sup>13</sup>CO (2–1) and C<sup>18</sup>O (2–1) seem to trace mostly Keplerian rotation at a position angle of approximately 70°–80°. We performed astrometry of the channel maps of <sup>13</sup>CO and C<sup>18</sup>O (presented in the Appendix) in order to determine the position angle of the disk by measuring the centroid position of the redshifted and blueshifted emission of each line independently. Only channels where the emission has an S/N > 5 were included in the fit. This resulted in seven data points for <sup>13</sup>CO and three data points for C<sup>18</sup>O. All data points align well with the peak of the continuum (Figure 4, right panel). By performing a $\chi^2$ fit, we determined the PA to be 76° ± 3° and 71° ± 3° for <sup>13</sup>CO and C<sup>18</sup>O, respectively. Fitting both the <sup>13</sup>CO and C<sup>18</sup>O data simultaneously yielded a PA = 74° ± 3°.

Like the <sup>13</sup>CO and C<sup>18</sup>O first moments, the <sup>12</sup>CO first moment (Figure 3) has a velocity gradient along the disk major axis typical of disk rotation, although with deviations from a pure Keplerian pattern (and at a scale larger than for the <sup>13</sup>CO and C<sup>18</sup>O emission). The <sup>12</sup>CO first moment also shows a velocity gradient perpendicular to the disk’s major axis, consistent with the presence of an outflow. The <sup>12</sup>CO channel maps (Figure 5) reveal an asymmetric, rotating structure, near a systemic velocity somewhere between +3.6 and +4.6 km s<sup>-1</sup>. There is a “twist” at systemic velocities seen in the <sup>12</sup>CO
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left panel: ALMA $^{12}$CO (2–1) integrated intensity (moment 0) image of EX Lup. Contours start at $3\sigma$ with steps of $10\sigma$ ($\sigma = 9$ mJy beam$^{-1}$ km s$^{-1}$). Middle panel: integrated intensity of $^{13}$CO (2–1). Contours start at $3\sigma$ with steps of $5\sigma$ ($\sigma = 7$ mJy beam$^{-1}$ km s$^{-1}$). Right panel: C$^{18}$O (2–1). Contours start at $3\sigma$ with steps of $1\sigma$ ($\sigma = 5$ mJy beam$^{-1}$ km s$^{-1}$). The $^{12}$CO (2–1) data have been displayed using a sqrt stretch in order to highlight faint extended emission. The color bar units are Jy beam$^{-1}$ km s$^{-1}$. The star symbol in the center of each panel represents the stellar position, which coincides with the peak of the continuum emission. The velocity range used to make each map is annotated in the top right. The dashed white circle shows the region used to compute the integrated line emissions listed in Table 1.</figcaption></figure>

first-moment image (Figure 3). This is similar to the infalling-
rotating motion from a tenuous envelope into a disk observed in
Class 0/I sources (e.g., Yen et al. 2014; Oya et al. 2016) and
recently discovered in Class II sources (Yen et al. 2017). The red
side of the disk ($v > +4.4$ km s$^{-1}$) exhibits the “butterfly”
pattern typical of a rotating Keplerian disk seen at intermediate
inclination. This is in agreement with previous works that have
estimated the disk inclination angle as between about 20° and 50°
(Sipos et al. 2009; Goto et al. 2011; Kóspál et al. 2014; Sicilia-
Aguilar et al. 2015) and is consistent with our estimate of 32°
from the continuum image. Optical absorption-line measurements
(Kóspál et al. 2014) suggest an optical systemic local standard of
rest (LSR) velocity of $\sim +4.4$ km s$^{-1}$, in reasonable agreement
with centroid velocities of the CO line profiles, particularly for
the more optically thin and spatially compact $^{13}$CO and C$^{18}$O
lines (Table 1). The line centroids determined from Gaussian
fitting (Table 1) are virtually identical to the line positions
measured by Kóspál et al. (2016b) using single-dish APEX
spectra (3.7, 3.9, and 5.0 km s$^{-1}$ for ALMA versus 3.6, 3.9, and
$5.1$ km s$^{-1}$ for APEX).

To further highlight the differences between the blue and red
sides of the $^{12}$CO (2–1) emission, Figure 6 shows the integrated
intensity maps at systemic-, intermediate-, and high-velocity
ranges. For comparison, moment 0 maps for $^{13}$CO (2–1) and
C$^{18}$O (2–1) in identical velocity ranges are presented in the
Appendix. The blue side of the disk is strikingly different from
the red side. Instead of mirroring the redshifted side, the
blueshifted emission is dominated by a bright, compact region
at systemic and intermediate velocities. The blue side shows an
additional extended spiral-shaped component shown in
Figure 6 up to 1.2 in extent and also visible in the +2.6 and
+2.8 km s$^{-1}$ channel maps. At higher blueshifted velocities
($v < +1.6$ km s$^{-1}$) there is emission near the central star as a
counterpart to the high-velocity emission seen on the red side.
Another feature of the line maps is an extended component at
systemic and intermediate redshifted velocities, extending up to
1.5 from the star.

The compact high-velocity emission on the blue side is
fainter than its redshifted counterpart, and the redshifted
emission extends to higher relative velocities than does the
blueshifted counterpart (assuming a systemic velocity of
$\sim +4.1$ km s$^{-1}$). The redshifted high-velocity gas, which is

located very close to the central star, can be traced down to
+9.4 km s⁻¹, whereas the blueshifted compact emission is
detected only out to +0.2 km s⁻¹. The overall rotation follows
a slightly larger PA than derived for $^{13}$CO and C$^{18}$O. We
performed an astrometric analysis for $^{12}$CO (2–1) similar to that
for the two isotopologues, but selecting high-velocity channels
between +0.6 and +1.6 km s⁻¹ for the blue side and between
+6.8 and +8.6 km s⁻¹ for the red side. This resulted in a PA of
80° ± 2° (Figure 4).

Between +1.2 and +1.6 km s⁻¹, we note extended, arc-like
emission located northwest of the star. It is located 0.′′8 (125 au)
away from the disk plane, in the direction of the rotational axis of
the system, where molecular outflows are typically seen. We
interpret this arc-like feature as the interaction of a molecular
outflow with remnant ambient material. Such arc-like shapes
have been discovered recently in other outbursting sources
(Kóspál et al. 2017b; Ruíz-Rodríguez et al. 2017; Zurlo et al.
2017; Principe et al. 2018), although mostly around FUors and
one transitional FUor/EXor source (Principe et al. 2018).

Given the strong asymmetry between the redshifted and
blueshifted sides, determining the exact systemic velocity is
difficult. The typical “hourglass” shape from an inclined
Keplerian disk seen at the systemic velocity is normally used to
identify the central velocity using channel maps. In the EX Lup
channel maps, channels between +3.6 and +4.4 km s⁻¹ could
plausibly be attributed to the systemic velocity, with
+3.6 km s⁻¹ perhaps slightly favored. However, not only is
the hourglass shape slightly twisted, but the center of the
hourglass is not located at the position of the continuum as is
typical of circumstellar disks (see Section 5.4). The double-
lobed pattern only becomes centered with the continuum
centroid at velocities of +4.0 to +4.2 km s⁻¹, which is closer
to the optically measured system velocity. The $^{12}$CO spectrum
(Figure 4) does not exhibit a double-peaked profile as expected
from a purely rotating, inclined, symmetric disk. Instead, it
shows a single-peaked profile centered at +3.7 km s⁻¹. The
$^{12}$CO (2–1) line shape is similar to the $^{12}$CO (3–2) line profile
from APEX single-dish observations (Kóspál et al. 2016b),
which also peaks at +3.6 km s⁻¹.

Position–velocity (PV) diagrams for all three isotopologues
were constructed by extracting the averaged spectra along the
disk major axis (assuming a PA of 78°). The features of the
<figure><img src="image_6.png" /><figcaption>Figure 3. Left to right: moment 1 maps for $^{12}$CO (2–1), $^{13}$CO (2–1), and C$^{18}$O (2–1), respectively. Only pixel values above 3σ have been included. Contour levels are 0.9, 1.8, 2.8, 3.7, 4.7, 5.6, 6.5, and 7.5 km s$^{-1}$.</figcaption></figure>

$^{12}$CO (2–1) maps are evident in the PV diagram (Figure 7): the redshifted emission is slightly fainter than the blue and shows compact high-velocity emission near the star. The blueshifted emission from part of the high-velocity arc-like feature is evident at $r > 100–200$ au. The emission on the redshifted side seems consistent with a Keplerian disk of inclination close to 40° (i.e., the value adopted by Kóspál et al. [2016b] for modeling the single-dish CO spectra). In Figure 7, we extract the centroid velocities as determined from a Gaussian fit to the $^{12}$CO (2–1) line profiles on both sides of the disk and plot them against position. The velocities from the redshifted side can be fit with an $r^{-0.53}$ power law, consistent with Keplerian rotation. However, the blueshifted velocities deviate significantly from Keplerian rotation and increase with radius at $r > 55$ au. In Figure 8 we use the redshifted data points from the inner $r < 100$ au from Figure 7 to fit the inclination angle, assuming Keplerian rotation around a star of mass 0.5 $M_{\odot}$. The inclination angle derived from $\chi^2$ fitting is 38° ± 4°, thus consistent with the value derived by Kóspál et al. (2016b).

The $^{13}$CO and C$^{18}$O data do not show any clear asymmetries in the PV diagrams (Figure 9) and are consistent with a Keplerian disk with a systemic velocity between +4.0 and +4.4 km s$^{-1}$. There is no evidence of outflows/inflows or additional emission from the arc feature in $^{13}$CO or C$^{18}$O, suggesting that the arc CO gas has a relatively low column density compared to the rest of the disk.

In the following section we use continuum and line radiative transfer modeling to characterize the physical properties of the EX Lup disk.

# 4. Disk Modeling

## 4.1. Radiative Transfer Modeling of Dust Continuum

The EX Lup continuum emission presented here has no resolved substructures and can be approximated as a continuous axisymmetric disk. In order to derive disk structural parameters and the dust mass, the continuum was modeled using the radiative transfer code RADMC-3D (Dullemond et al. 2012). We adopt the standard tapered-edge model to describe the surface density profile (Andrews et al. 2009; Williams & Cieza 2011) parameterized by the total dust mass $M_{d}$, a power-law slope $\gamma$, a characteristic radius $R_{c}$, and a surface density $\Sigma$ that is tapered with an exponential decay

beyond $R_{c}$:

$$ \Sigma = \Sigma_c \left(\frac{R}{R_c}\right)^{-\gamma} \exp\left[-\left(\frac{R}{R_c}\right)^{2-\gamma}\right]. \quad (1) $$

The vertical scale height of the disk as a function of radius is also a power law $H(r) = H_{100}(r/100 \text{ au})^{\psi}$, with a flaring index $\psi = 1.09$ as determined by Sipos et al. (2009). A stellar spectrum with $T_{\text{star}} = 3859$ K and a radius of 0.8 $R_{\odot}$ were adopted (Frasca et al. 2017). The model includes a standard power-law distribution of dust grain sizes $a$, given by $n(a) \propto a^{-3.5}$ and extending from $a_{\text{min}} = 0.1 \mu\text{m}$ to $a_{\text{max}} = 3.0 \text{ mm}$. The optical properties of the dust correspond to a mix of amorphous carbon grains (Li & Greenberg 1997) and astro-silicate grains (Draine & Lee 1984), which were combined using Bruggeman’s rules. We used a standard mix of 30% amorphous carbons and 70% astro-silicates. The opacities of the mix were computed using the “Mie Theory” code written by Bohren & Huffman (1983). The mass opacity at 1.3 mm is thus $\kappa_{\text{abs}} = 2.2 \text{ cm}^2 \text{ g}^{-1}$.

The inner radius was set to approximately the expected dust sublimation radius of ~0.05 au. This is smaller than the value used by Sipos et al. (2009), who found that an inner radius smaller than 0.2 au caused an excess emission in the near- and mid-infrared part of the SED. We caution that this radius depends on disk parameters for an optically thick disk, but our current data do not constrain this radius. Figure 1 shows that our model fits the resolved image of the disk at 1.3 mm, and it can be reconciled with the SED at shorter wavelengths if the small and large grains follow different spatial distributions (as also noted by Sicilia-Aguilar et al. 2015). Large grains accumulate in the inner regions (as predicted by radial drift models; e.g., Pinte & Laibe 2014), while the inner <0.2 au region is depleted of smaller grains.

### 4.1.1. Fitting Procedure

The model parameters that we aim to determine are $M_{d}$, $\gamma$, $R_{c}$, and $H_{100}$. Before exploring the parameter space, the centroid shift, the inclination angle $i$, and the PA of the continuum emission were determined using the CASA task IMFIT. The disk structural parameters were then constrained using a Bayesian approach. The posterior distribution for each parameter was recovered using the Goodman & Weare (2010) affine-invariant MCMC ensemble sampler as implemented by
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left panel: integrated spectrum for $^{12}$CO (2–1) (top), $^{13}$CO (2–1) (middle), and C$^{18}$O (2–1) (bottom). Spectra were computed by integrating the line images over apertures of 1′′7, 0′′75, and 0′′55 for $^{12}$CO, $^{13}$CO, and C$^{18}$O, respectively. The red lines show the spectra for the gas disk model presented in Section 4.2. The vertical dotted line represents the system’s velocity, while the horizontal blue, orange, and green lines represent the velocity ranges used to construct the moment 0 maps at different velocity intervals from Figure 6. Middle panel: azimuthally averaged surface brightness profiles for $^{12}$CO, $^{13}$CO, and C$^{18}$O moment 0 images. The dashed and dotted lines represent the 3σ and 1σ levels of each map. $^{13}$CO and C$^{18}$O have been multiplied by a factor of 5 for display purposes. Right panel: spectro-astrometry measurements of the $^{13}$CO (2–1) and C$^{18}$O (2–1) channel maps (open and filled black circles, respectively). Each data point corresponds to the centroid of the line emission relative to the center of the continuum emission. Only channels with emission above an S/N of 5 have been included. The black line corresponds to a χ² fit to the combined $^{13}$CO (2–1) and C$^{18}$O (2–1) data. The resulting position angle is annotated in the image, together with the derived uncertainty from the fit. The error bars correspond to the astrometric accuracy of the data. Red points show the spectro-astrometry measurements of the $^{12}$CO (2–1) high-velocity channel maps relative to the center of the continuum emission. Velocity channels between +0.6 and +1.6 km s⁻¹ were selected for the blue side, while channels between +6.8 and +8.6 km s⁻¹ were selected for the red side. The red dashed line shows the best fit to the $^{12}$CO (2–1) data.</figcaption></figure>

Foreman-Mackey et al. (2013). We used the latter’s publicly available PYTHON module EMCEE to sample the parameter space and maximize the likelihood function. The likelihood function is proportional to $exp(-χ^2/2)$, where χ² is the sum over the squared difference of the model and measured visibilities divided by their variance. The model and the data (the self-calibrated Measurement Set including only the continuum) were compared in the uv plane. The model visibilities were obtained by taking the fast Fourier transform of model images and interpolating to the same uv points as the observations (S. Marino et al. 2018, in preparation).

Our priors for the free parameters came from assuming uniform distributions, specifically

$$ 
\begin{array}{l}
M_d \in [0.001, 0.1] \times 10^{-2} M_{\odot} \\
\gamma \in [-2.0, 2.0] \\
R_c \in [5.0, 100.0] \text{ au} \\
H_{100} \in [0.1, 10.0] \text{ au.}
\end{array}
$$

The best-fitting parameters and their uncertainties were obtained after running 1000 iterations (~10 times the autocorrelation time) with 500 walkers. The posterior distributions of $M_{\text{dust}}$, $\gamma$, $R_c$, and $H_{100}$ are presented in Figure 10.

## 4.2. Comparison with Axisymmetric CO Gas Disk Model

In this section we compare the EX Lup line observations to an axisymmetric CO gas disk model, with the purpose of highlighting the nonsymmetric $^{12}$CO (2–1) line emission, as well as reproducing the $^{13}$CO (2–1) and C$^{18}$O (2–1) emission. For this purpose, we use the 3D Monte Carlo radiative transfer code MCFOST (Pinte et al. 2006, 2009) to produce model images of $^{12}$CO (2–1), $^{13}$CO (2–1), and C$^{18}$O (2–1) for a disk model that yields reasonable agreement with the observed $^{13}$CO and C$^{18}$O emission. We assumed that the low- J CO level populations are in LTE at a temperature of $T_{\text{gas}} = T_{\text{dust}}$ in order to produce line emission data cubes. The abundance of $^{12}$CO relative to H₂ was set to $10^{-4}$ and assumed constant throughout the disk. Freeze-out of CO molecules onto dust grains was assumed to occur at temperatures below 20 K with a depletion factor of 1000.

We used a genetic algorithm (e.g., Mathews et al. 2013; Pinte et al. 2018) to fit the $^{13}$CO and C$^{18}$O moment 0 maps simultaneously. The CO disk model does not attempt to fit the nonaxisymmetric $^{12}$CO (2–1) line emission, which is beyond the scope of this paper. Each generation was composed of 100 models, and the genetic algorithm was run for 50 generations. The following parameters were allowed to vary during the fitting process: $M_{\text{gas}}$, $\gamma$, $R_{\text{in}}$, and $H_{100}$, as well as the disk’s flaring exponent. Initially $R_c$ and the disk’s outer radius $R_{\text{out}}$ were also part of the variable parameters, but we found that better solutions were obtained if we restricted $R_c$ to 75 au and the outer radius to 150 au. To determine these values, a range of $R_c$ between 50 and 90 au and a range of $R_{\text{out}}$ between 100 and 200 au were explored by eye, in steps of 5 and 10 au, respectively. Furthermore, because the algorithm does not fully sample the parameter space, the solution found by the genetic algorithm may not be unique, and therefore one must be careful when interpreting the model parameters presented in Table 2.

The CO disk model is different from the dust disk model presented in Section 4.1: both the dust and CO disk models assume the same total dust mass from Section 4.1, but the radial distribution of dust in the CO disk model is assumed to be more radially extended, as described by the characteristic radius parameter $R_c$ for the gas model. Differences in the radial distribution millimeter-sized dust grains and CO gas have been reported in protoplanetary disk systems (de Gregorio-Monsalvo et al. 2013; Pinte et al. 2016) and can be explained by faster grain growth in the central regions and/or inward radial migration of larger dust grains (Barrière-Fouchet et al. 2005; Pinte & Laibe 2014).

The model images were used to create synthetic visibilities using CASA task SIMOBSERVE, using the same integration time, spectral setup, and antenna configuration used for the observations, as well as injecting the appropriate amount of thermal noise. The resulting model visibilities were then imaged using the same CLEAN parameters used for the real data. Figures 11 and 12 show the resulting Moment 0 images for the data, model, and residuals for $^{13}$CO (2–1) and C$^{18}$O (2–1), respectively. The model reproduces well the observed
<figure><img src="image_8.png" /><figcaption>**Figure 5.** $^{12}$CO channel maps toward EX Lup. The velocity of the channels is shown in the LSR frame, centered at the rest frequency of $^{12}$CO (2–1). The data have been binned to a velocity resolution of 0.2 km s$^{-1}$. Contours of the $^{12}$CO emission are overlaid. Contour levels are 0.018, 0.036, 0.054, 0.072, 0.09, 0.11, 0.13, and 0.15 Jy beam$^{-1}$. The star symbol in the center of each panel represents the stellar position, coinciding with the peak of the continuum emission.</figcaption></figure>

$^{13}$CO (2–1) and C$^{18}$O (2–1) integrated intensity images. We find that a vertical optical depth of 1 is reached at a radius of 60 and 15 au for $^{13}$CO (2–1) and C$^{18}$O (2–1), respectively. $^{12}$CO (2–1) is very optically thick throughout the disk and has an optical depth of 1.5 even at the outer radius of 150 au.

# 5. Discussion

## 5.1. Dust Disk

The ALMA observations of EX Lup resolve the dust disk for the first time. The total dust mass derived from the radiative transfer model fitting process is $1.0 \times 10^{-4} M_{\odot}$ (33.3 $M_{\odot}$). Our

model predicts, however, that the disk becomes optically thick at $r < 10$ au, so the total dust mass could be higher by a factor of 2–5. Deriving the dust mass using the standard assumption of optically thin emission, a dust temperature of 20 K, and the opacities of Beckwith et al. (1990) underestimates the dust mass by a factor of at least 3 compared to our model. The model predicts a 0.87 mm flux of ~50 mJy, consistent with the flux measured with LABOCA (Juhász et al. 2012). The 1.3 mm spectral index measured in-band and also that between 0.87 and 1.3 mm are both close to the value of 2.3 measured toward FU Ori itself (Liu et al. 2017). The low spectral index of FU Ori can be explained by the presence of a compact optically
<figure><img src="image_9.png" /><figcaption>**Figure 6.** $^{12}$CO (2–1) moment 0 maps at different velocity ranges. The high-velocity channels (left panel) have been integrated between −1.6 and +1.7 km s$^{-1}$ (blue) and between +6.5 and +9.8 km s$^{-1}$ (red). The gray scale shows 1.29 mm continuum emission from the compact dust disk (in units of Jy beam$^{-1}$). Intermediate-velocity channels (middle panel) have been integrated between +1.8 and +3.5 km s$^{-1}$ (blue) and between +4.7 and +6.4 km s$^{-1}$ (red). Right panel: integrated emission at systemic velocities (+3.6 to +4.6 km s$^{-1}$). Contour levels are the same as for Figure 5.</figcaption></figure>

thick inner region surrounded by a more extended optically thin disk (Liu et al. 2017). A similar spectral behavior has been reported in the bright FUor system V883 Ori (Cieza et al. 2016), which shows an axisymmetric optically thick inner disk within 42 au surrounded by an optically thin halo extending out to 125 au. The modeling procedure applied to EX Lup is similar to the one we performed for other FUor/EXor sources (Cieza et al. 2017), allowing direct comparison. EX Lup has an intermediate dust mass compared to other FUor and EXor objects (Liu et al. 2016; Cieza et al. 2017). The ~23 au characteristic radius we measure is consistent with the compact radii measured toward other FUor/EXor sources. The γ value of 0.34 is lower than measured toward most FUor objects observed by Cieza et al. (2017) and closer to the values seen around the FUor/EXor object V1647 Ori and toward T Tauri stars (Andrews et al. 2010). We find that the disk around EX Lup follows the trend that EXor stars are surrounded by compact dust disks that become optically thick in their inner regions (Cieza et al. 2016, 2017; Liu et al. 2017), which could be the reservoir feeding the accretion.

We find no evidence of strong asymmetries in the dust continuum, which have been suggested to occur as a direct consequence of gravitational instability (e.g., Zhu et al. 2012). Comparison with an axisymmetric model suggests the presence of structure not accounted for at the 5σ level. Assuming a gas-to-dust ratio of 100, the total disk mass is $0.011 M_\odot$, close to the $0.025 M_\odot$ derived with SED fitting by Sipos et al. (2009). Given that gravitational instability and disk fragmentation require $M_D/M_* > 0.1$ (e.g., Armitage et al. 2001; Vorobyov & Basu 2015), the total disk mass derived from the continuum is consistent with the nondetection of asymmetries in the dust disk triggered by gravitational instability.

## 5.2. CO Gas Disk

We detect widespread $^{12}$CO (2–1) emission, of which some fraction can be attributed to a Keplerian gas disk. Williams & Best (2014) developed a grid of protoplanetary disk models that can be used to derive disk masses based on the $^{13}$CO (2–1)/$C^{18}$O (2–1) integrated line ratios (the models take into account photodissociation, CO freeze-out, and basic CO chemistry). Comparing the EX Lup $^{13}$CO/$C^{18}$O line ratios to those from the models of Williams & Best (2014), we obtain a total disk mass of $5.4 \times 10^{-4} M_\odot$. This result is remarkably

similar to the $5.1 \times 10^{-4} M_\odot$ we derived by modeling of the $^{13}$CO and $C^{18}$O moment 0 images in Section 4.2. This is an order of magnitude lower than the $(1-3) \times 10^{-3} M_\odot$ total disk mass derived by Sicilia-Aguilar et al. (2015) from revised SED fitting, but close to the value of $2.3 \times 10^{-4} M_\odot$ derived by Kóspál et al. (2016b) based on the single-dish detection of the $^{13}$CO (3–2) line. We note that since in our model the $^{13}$CO and $C^{18}$O become optically thick at radii of 60 and 15 au, respectively, the resulting gas mass estimate may only represent a lower limit to the total gas mass. Combined with the dust modeling results, our gas mass estimates imply a gas-to-dust ratio of 4.5, which is similar to the value reported by Kóspál et al. (2016b). Gas-to-dust ratios lower than the typical interstellar medium value of 100 using Williams & Best (2014) seem to be common around disks in Lupus and Chameleon (Ansdell et al. 2016; S. Villenave et al. 2018, in preparation). Depletion of CO and other molecular species by up to an order of magnitude during accretion outbursts in EX Lup has been detected through mid-infrared spectroscopy (Banzatti et al. 2015) and could possibly be related to the low gas-to-dust ratios we derived based on the $^{13}$CO and $C^{18}$O lines.

## 5.3. Widespread $^{12}$CO (2–1) Gas Emission

We detect $^{12}$CO (2–1) emission up to the full extent of the MRS (~3"/2), suggesting that there could be larger-scale emission filtered out by the interferometer. Unfortunately, there are no $^{12}$CO (2–1) single-dish observations available that could have been used to estimate the amount of resolved-out emission.

High accretion rates in EXor systems must be accompanied by high mass-loss rates (Hartmann 2008; Hartmann et al. 2016). Prominent collimated outflows, as seen around FUor objects (Fehér et al. 2017; Kóspál et al. 2017b; Ruíz-Rodríguez et al. 2017; Zurlo et al. 2017; Principe et al. 2018), are not evident in the EX Lup data. However, we identify an arc-like feature in $^{12}$CO blueshifted emission, 100–200 au northwest of the star (see Figures 5 and 6). This emission is detected between +0.6 and +1.8 km s$^{-1}$, i.e., displaced to the blue of the systemic velocity by at least ~2.0 km s$^{-1}$. This gas and its corresponding velocity gradient are both in a direction perpendicular to the disk’s major axis, and we conclude that this blue emission corresponds to a molecular outflow, possibly interacting with ambient material. We expect this CO emission
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left panel: ¹²CO position–velocity diagram along the disk major axis (assumed PA = 78°). The horizontal dashed white line represents the possible systemic velocity of +4.1 km s⁻¹, while the vertical dotted black line represents the stellar position. The red and blue markers show measured centroid positions in different velocity channels for the red- and blueshifted sides of the disk. The green curves show the Keplerian velocity for gas at the tangential point at an inclination angle of 38°. The mass of the central object was assumed to be 0.5 M☉ (Frasca et al. 2017). For comparison, the Keplerian velocities for a 0.2 and 0.8 M☉ central star are also shown (inner and outer dashed white lines, respectively). Right panel: rotation profile of the ¹²CO redshifted (red) and blueshifted (blue) emission. The dashed line represents the results from least-squares fitting to the inner r < 100 au redshifted data points, which results in an r⁻⁰.⁵³ slope. The blueshifted emission contrarily has a positive velocity gradient at r > 50 au.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Best fit to the inclination of angle to the inner r < 100 au redshifted data points from Figure 7. The stellar mass was assumed to have a value of 0.5 M☉.</figcaption></figure>

more likely to pertain to swept-up molecular material rather than material launched directly from the source during the recent burst, because the latter scenario would require velocities of 100 km s⁻¹ in order for the material to have traveled 150 au in 8 yr.

Accretion-related winds have been detected through optical spectroscopy in EX Lup during the periods of enhanced accretion and are believed to originate in the inner r < 0.4 au region of the disk (Goto et al. 2011; Kóspál et al. 2011; Sicilia-Aguilar et al. 2012). This probable detection of an outflow provides evidence that the bipolar outflow phenomenon persists into the EXor phase and that accretion/outflow is still clearing molecular gas from the areas above and below the disk plane. EXor/FUor sources are therefore likely to be valuable for studying episodic infall/ accretion/outflows near the end stages of the outflow

phenomenon. These larger-scale outflow phenomena are still relevant to the evolution of the disk, as they regulate the accretion process via removal of angular momentum.

We do not detect a clear redshifted counterpart to the blueshifted outflow. There is redshifted emission at high velocities close to the star in excess of our Keplerian disk model, which is not seen on the blue side. This could be material either outflowing or infalling toward the star, although we believe that it is most likely the redshifted counterpart of the molecular outflow.

This interpretation places the disk’s far (near) side to the northwest (southeast), with the exposed surface rotating clockwise as seen by the observer. This is consistent with the observed disk kinematics, best evidenced in the C¹⁸O first moment map (Figure 3). The proposed geometry for the disk and molecular outflow is shown in Figure 13. This geometry is consistent with the inner disk configuration proposed by Sicilia-Aguilar et al. (2012).

Molecular outflows and complex circumstellar structures have been observed in FUor sources and in the borderline FUor/EXor object V1647 Ori, but not in EXors: this has led to the suggestion that FUor objects may be precursors to EXors, with EXors corresponding to a slightly more evolved stage (Hartmann 2008; Cieza et al. 2017). EX Lup itself seems to have a moderately complex circumstellar environment, with evidence of recent outflowing material still impacting the tenuous remnants of the ambient cloud and extended material near systemic velocities possibly still spiraling into the system. Compared with the rich, active environments of FUor sources and the more quiescent environments of the few surveyed EXor objects, EX Lup itself may be a system in a transitional phase between classical FUors and EXors. Ironically, EX Lup and FU Ori seem more similar to each other in terms of disk size and dust mass than FU Ori is to other FUor sources or EX Lup is to other EXors (Hales et al. 2015; Liu et al. 2016; Cieza et al. 2017; Liu et al. 2017).

Fainter extended material at low velocity is detected on the redshifted side of the disk (Figure 7). This extended material has
<figure><img src="image_12.png" /><figcaption>Figure 9. Position–velocity diagram for the observed $^{13}$CO (left panel) and C$^{18}$O (right panel) along the disk major axis (assuming a PA of 78°). The horizontal dashed white line represents the possible systemic velocity of +4.1 km s$^{-1}$, while the vertical dotted black line represents the stellar position. The green curves show the Keplerian velocity for gas at the tangential point at an inclination angle of 38° and a stellar mass of 0.5 $M_{\odot}$. The Keplerian velocities for a 0.2 and 0.8 $M_{\odot}$ central star are also shown (inner and outer dashed black lines, respectively).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10. Posterior distributions of each dust disk parameter, including their marginalized distributions for EX Lup. The vertical dashed lines represent the 16th, 50th, and 84th percentiles. Contours correspond to the 68%, 95%, and 99.7% confidence regions.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11. Left to right: $^{13}$CO (2–1) moment 0 images for data, model, and residual, respectively. Contour levels are the same as in Figure 2.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 12. Left to right: C$^{18}$O (2–1) moment 0 images for data, model, and residual, respectively. Contour levels are the same as in Figure 2.</figcaption></figure>

<table><caption>Table 2
EX Lup CO Disk Model</caption><thead><tr><th>Physical Parameter</th><th></th><th>Reference</th></tr></thead><tbody><tr><td colspan="3">Stellar properties</td></tr><tr><td>Primary star mass, M<sub>*</sub> (M<sub>☉</sub>)</td><td>0.5</td><td>Frasca et al. (2017)</td></tr><tr><td>Effective temperature, T<sub>*</sub> (K)</td><td>3859</td><td>Frasca et al. (2017)</td></tr><tr><td>Stellar radius, R<sub>*</sub> (R<sub>☉</sub>)</td><td>0.8</td><td>Frasca et al. (2017)</td></tr><tr><td>Distance, d (pc)</td><td>155</td><td>Lombardi et al. (2008)</td></tr><tr><td colspan="3">Disk structure<sup>a</sup></td></tr><tr><td>Disk dust mass, M<sub>d</sub> (M<sub>☉</sub>)</td><td>1.1 × 10<sup>-4</sup></td><td>1.1<sup>+0.3</sup><sub>-0.3</sub> × 10<sup>-4</sup></td></tr><tr><td>Disk total mass, M<sub>D</sub> (M<sub>☉</sub>)</td><td>5.1 × 10<sup>-4</sup></td><td>100 × M<sub>d</sub></td></tr><tr><td>Inner rim radius, R<sub>in</sub> (au)</td><td>0.27</td><td>0.05</td></tr><tr><td>Characteristic radius, R<sub>c</sub> (au)</td><td>75</td><td>22.9<sup>+1.2</sup><sub>-0.9</sub></td></tr><tr><td>Outer radius, R<sub>out</sub> (au)</td><td>150</td><td>...</td></tr><tr><td>Characteristic height at 100 au, H<sub>c</sub> (au)</td><td>7.6</td><td>2.2<sup>+0.3</sup><sub>-0.2</sub></td></tr><tr><td>Surface density exponent, γ</td><td>0.32</td><td>0.31<sup>+0.09</sup><sub>-0.10</sub></td></tr><tr><td>Flaring exponent, ψ</td><td>1.0</td><td>1.1</td></tr><tr><td>Inclination angle, i (deg)</td><td>38</td><td>32</td></tr><tr><td>Position angle, PA (deg)</td><td>78</td><td>65</td></tr><tr><td>Systemic velocity, V<sub>LSR</sub> (km s<sup>-1</sup>)</td><td>+4.1</td><td></td></tr></tbody></table>

**Note.**

a Uncertainties are only given for parameters fit in the MCMC (see Section 4).

a radial velocity profile $\propto r^{-2.9}$, much steeper than the nearly $r^{-0.5}$ Keplerian dependence found in the inner $r < 100$ au, and also steeper than the $r^{-1}$ gradient characteristic of flattened, infalling envelopes (e.g., Kóspál et al. 2017b, and references therein). It must be noted that the redshifted emission at larger radii is fainter

(at the limit of the MRS), and thus the fit to the centroid positions is less accurate. This diffuse extended emission is likely associated with a more tenuous rotating envelope, which has to be low density to be consistent with the low-reddening and disk-like SED of EX Lup.

The first moment map also shows a twist around the systemic velocity, similar to those found in sources where material is still being transferred from the envelope into the disk (Yen et al. 2017). This twist corresponds to a spiral-shaped structure that seems to cross the disk at velocities between +2.6 and +3.4 km s$^{-1}$. The two sides of the spiral cross the disk at the center of the blue-side asymmetry. However, interpretation of structure near the systemic velocity may be confused by extended cloud emission larger than the MRS, causing problems during the image reconstruction.

## 5.4. Nature of the Asymmetry and Connection with the EXor Outburst

The $^{12}$CO (2–1) channel maps reveal clear deviations from Keplerian kinematics when compared to an axisymmetric disk model (Figure 14). The “butterfly” pattern that is clearly seen on the redshifted side of the disk should be symmetric on the blueshifted side with respect to the systemic velocity. If the systemic velocity is assumed to be around +4.1 km s$^{-1}$, then for an axisymmetric disk the channel maps at +3.0 and +5.0 km s$^{-1}$ should be mirror images of each other. Instead, they are very different, with the +3.0 km s$^{-1}$ map being more compact than its +5.0 km s$^{-1}$ red-side counterpart. The velocity channels at +3.2, +3.4, and +3.6 km s$^{-1}$ show an elongated, hourglass shape that seems to cross the disk in the north–south direction. The hourglass shape does not cross through the disk’s center (i.e., the
<figure><img src="image_16.png" /><figcaption>Figure 13. Schematic diagram of the EX Lup system.</figcaption></figure>

stellar position and centroid of the continuum), but at ~0.08 to the southwest (~1/3 of the beam). Figure 15 shows a blowup of the channel maps near the systemic velocity of the $^{12}$CO (2–1) image cube created using Briggs weighting. The slightly higher resolution compared to natural weighting shows that the hourglass shape may consist of at least three unresolved clumps connected by filamentary emission.

One possible scenario to explain the morphology seen on the blue side is that we are observing an asymmetry related to the outflow, the offset corresponding to the launching radius of the disk wind. The spiral seen near the systemic velocity seems to “land” at this position, consistent with the asymmetry coinciding with the landing site of infalling material from the remnant envelope and possibly to the launching site of a disk outflow. This is similar to what has been reported by Alves et al. (2017) for the Class I object BHB07-11. Although more massive ($1.7 M_\odot$) and younger than EX Lup, ALMA observations of BHB07-11 detect bipolar outflows launched symmetrically from the disk at a radius of 80 au from the star, where larger-scale spirals and the disk intersect. In the case of EX Lup, the center of the two-lobed structure seen at $+3.6 \text{ km s}^{-1}$ is located at 12–15 au, requiring velocities of ~10 km s$^{-1}$ to escape the system.

Alternatively, twisted iso-velocities similar to the morphology seen at $+3.6 \text{ km s}^{-1}$ are associated with enhanced accretion onto the central star (Rosenfeld et al. 2014; Casassus et al. 2015). In this scenario the twisted iso-velocities should be centered at the stellar position; however, in our data the center of the twisted hourglass is shifted by about 1/3 of a beam from the star. Such a shift in the accretion center could plausibly be the result of an unseen companion. Many efforts have been made to identify a possible binary companion around EX Lup (Ghez et al. 1997; Bailey 1998; Kóspál et al. 2014). Kóspál et al. (2014) detect a radial velocity signal that can be described by a companion with an $m \sin i = 14.7 \pm 0.7 M_{\text{Jup}}$ ($0.018 M_\odot$ assuming an inclination of 50° located 0.06 au from the primary. Sicilia-Aguilar et al. (2015) however, show that the same signal can also be reproduced by a rotating, line-dependent veiling from an accretion column, without the need of a stellar companion. In addition, a Jupiter-mass companion at 0.06 au may be too close-in to explain the shift in the $^{12}$CO emission centroid.

A similar twisted feature has been proposed as a signature of ongoing planet formation by Perez et al. (2015). As an accreting giant planet develops its own circumplanetary disk, it distorts the global velocity field owing to its own distinct kinematics and accretion streams. This circumplanetary disk also forces some of the gas to follow horseshoe orbits, distorting the Keplerian pattern along its full orbit, producing a large-scale skewed pattern in the iso-velocities, similar to the one seen in EX Lup's systemic velocity channels (S. Perez et al. 2018, in preparation).

Asymmetries in the inner disk of EX Lup have been reported previously (Goto et al. 2011; Sicilia-Aguilar et al. 2012), associated with a hot spot in the inner 0.4 au of the disk. Our observations indicate that the gas disk is also asymmetric at outer radii (>100 au), suggesting a possible link between the outer disk structure and the accretion from the inner regions into the star. In this picture, the asymmetry observed in the blue regions corresponds to the region of the disk where low-density material from the envelope falls into the disk and piles up until it can be efficiently transferred inward. The infall dragging of the magnetic field lines can trigger the launching of disk winds (Zhao et al. 2016), which may enhance the asymmetry between redshifted and blueshifted sides. A mismatch between infall from the envelope and accretion onto the star has been suggested as the possible cause for the piling up of material in the disk, leading to disk instabilities and enhanced accretion (Bell & Lin 1994). Observational evidence of this has recently been reported in the FUor star V346 Nor by Kóspál et al. (2017b) using ALMA observations incorporating Morita Array short-baseline and single-dish observations. This combination of data from different array families allows characterization of the circumstellar environment on different spatial scales and to constrain better the envelope infall rates. Our ALMA 12 m array CO observations of EX Lup suggest the presence of extended emission, but observations targeting larger-scale emission are needed to quantify the balance between envelope mass infall and how material from the asymmetric outer disk is transferred into the inner disk and ultimately onto the star.

The clumpy substructure detected in CO inside the region of the asymmetry could represent the clumps and spiral filaments predicted by models of fragmentation in gravitationally unstable disks (Vorobyov & Basu 2015). The models predict that these features should be detectable in the gas and that the interaction between disk and clump fragments could be responsible for the nonperiodic accretion outbursts. The current total disk mass estimates suggest, however, that the disk is not massive enough for gravitational instabilities to occur. Moreover, the dust disk appears less asymmetric than might be expected. Detailed hydrodynamic modeling and higher-resolution observations of the EX Lup disk are required to quantify whether the observed structure could be explained by fragmentation driven by gravitational instability.

## 6. Conclusions

We present the first millimeter-wave resolved observations of the prototype EXor outbursting system EX Lup using ALMA. We detect and resolve emission from EX Lup's disk in both continuum and the CO isotopologues $^{12}$CO (2–1), $^{13}$CO (2–1), and $C^{18}O$ (2–1). Our main findings are as follows:
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Channel maps of the data, axisymmetric model, and residuals for the ¹²CO emission. Contour levels and intensity scale are identical to those for Figure 5, but including negative contours to highlight negative residuals. The star symbol in the center of each panel represents the stellar position. The beam is shown in the top left corner of the leftmost panels.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Blowup of the central velocity channels, highlighting the clumpy, non-Keplerian ¹²CO emission on the blueshifted side. The star symbol in the center of each panel represents the stellar position, coinciding with the continuum peak.</figcaption></figure>

1. The continuum disk is compact, with a characteristic radius of 23 au and a total dust mass of $1.0 \times 10^{-4} M_{\odot}$ (33 $M_{\odot}$). Our modeling shows that the 1.3 mm continuum is optically thin beyond ~10 au, but the data show no asymmetries or substructures above the 5$\sigma$ level.

2. The ¹²CO disk is highly asymmetric, with the blue side deviating significantly from Keplerian rotation. On the other hand, ¹³CO and C¹⁸O trace disk rotation and are used to estimate a disk gas mass of $5.1 \times 10^{-4}$. The red side of the disk follows Keplerian rotation. Assuming a stellar mass of 0.5 $M_{\odot}$, we derive an inclination angle of 38° ± 4°.

3. The overall ¹²CO emission shows various components additional to the disk, such as extended emission at systemic velocities and a blueshifted molecular outflow at 100–200 au. The detection of a molecular outflow around EX Lup provides evidence that episodic accretion occurs throughout the early stages of star formation and the amount of activity declines gradually with age of the YSOs; therefore, the distinction between FUors and EXors is one of degree, not of kind.

The bright asymmetry on the blueshifted side shows a clumpy, twisted structure that is likely associated with the launching site of the molecular outflow. However, ruling out other possible explanations such as disk instabilities, disk-companion interactions, or optical depth effects requires further modeling. Follow-up observations targeting larger-scale emission are crucial to understanding the mass transport balance from the remnant envelope onto the disk.

We thank the anonymous referee for their insightful review. We thank E. Fomalont and J. Williams for useful discussions and suggestions during the analysis of the data. This paper makes use of the following ALMA data: ADS/JAO.ALMA.2015.1.00200.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research made use of ASTROPY, a community-developed core PYTHON package for astronomy (Astropy Collaboration et al. 2013). This research used the facilities of the Canadian Astronomy
<figure><img src="image_19.png" /><figcaption>**Figure 16.** $^{13}$CO (2–1) channel maps toward EX Lup. The velocity of the channels is shown in the LSR frame, centered at the rest frequency of $^{13}$CO (2–1). The data have been binned to a velocity resolution of 0.3 km s$^{-1}$. All maps share the same linear color scale. Contour levels are 0.012, 0.016, 0.024, 0.032, and 0.04 Jy beam$^{-1}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** C$^{18}$O (2–1) channel maps toward EX Lup. The velocity of the channels is shown in the LSR frame, centered at the rest frequency of C$^{18}$O (2–1). The data have been binned to a velocity resolution of 0.8 km s$^{-1}$. All maps share the same linear color scale. Contour levels are 0.055, 0.010, and 0.015 Jy beam$^{-1}$.</figcaption></figure>

Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. S.P. and L.C. acknowledge support from the Millennium Science Initiative (Chilean Ministry of Economy) through grant Nucleus RC130007. S.P. acknowledges CONICYT-Gemini grant 32130007. L.C. acknowledges support from CONICYT-FONDECYT grant 1171246. This paper used the Brelka cluster, financed by Fondequip project EQM140101. C.P. acknowledges funding from the Australian Research Council (ARC) under the Future Fellowship no. FT170100040. This work was undertaken with the SOLA (Soul of Lupus with ALMA) collaboration centered at the Joint ALMA Observatory.

*Facility:* Atacama Large Millimeter/Submillimeter Array.

*Software:* Common Astronomy Software Applications (McMullin et al. 2007), RADMC-3D (Dullemond et al. 2012), MCFOST (Pinte et al. 2006, 2009), Astropy (Astropy Collaboration et al. 2013).

## Appendix

### Channel Maps

Figures 16 and 17 show the channel maps for $^{13}$CO (2–1) and C$^{18}$O (2–1), respectively. Figure 18 shows the corresponding moment 0 maps in different velocity ranges.
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Moment 0 maps in different velocity ranges. $^{13}$CO (2–1) and C$^{18}$O (2–1) are shown in the top and bottom rows, respectively. The high-velocity channels (left panels) have been integrated between −1.6 and +1.7 km s$^{-1}$ (blue) and between +6.5 and +9.8 km s$^{-1}$ (red). Intermediate-velocity channels (middle panels) have been integrated between +1.8 and +3.5 km s$^{-1}$ (blue) and between +4.7 and +6.4 km s$^{-1}$ (red). Right panels: integrated emission at systemic velocities (+3.6 to +4.6 km s$^{-1}$). Contour levels for $^{13}$CO and C$^{18}$O are the same as in Figures 16 and 17, respectively.</figcaption></figure>

## ORCID iDs

A. S. Hales https://orcid.org/0000-0001-5073-2849

M. Saito https://orcid.org/0000-0003-0769-8627

C. Pinte https://orcid.org/0000-0001-5907-5179

L. B. G. Knee https://orcid.org/0000-0002-9342-9003

A. Plunkett https://orcid.org/0000-0002-9912-5705

L. Cieza https://orcid.org/0000-0002-2828-1153

## References

Ábrahám, P., Juhász, A., Dullemond, C. P., et al. 2009, *Natur*, **459**, 224

Alves, F. O., Girart, J. M., Caselli, P., et al. 2017, *ApJ*, **603**, L3

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

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

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

Armitage, P. J. 2016, *ApJL*, **833**, L15

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

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

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

Bailey, J. 1998, *MNRAS*, **301**, 161

Banzatti, A., Pontoppidan, K. M., Bruderer, S., Muzerolle, J., & Meyer, M. R. 2015, *ApJL*, **798**, L16

Barrière-Fouchet, L., Gonzalez, J.-F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, *A&A*, **443**, 185

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

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

Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (New York: Wiley)

Bonnell, I., & Bastien, P. 1992, *A&A*, **401**, L31

Casassus, S., Marino, S., Pérez, S., et al. 2015, *ApJ*, **811**, 92

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

Cieza, L. A., Ruiz-Rodríguez, D., Perez, S., et al. 2017, arXiv:1711.08693

Comerón, F. 2008, in Handbook of Star Forming Regions, Vol. II, ed. B. Reipurth (San Francisco, CA: ASP), 295

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, *yCat*, 2246, 0

D'Angelo, C. R., & Spruit, H. C. 2010, *MNRAS*, **406**, 1208

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, *A&A*, **557**, A133

Draine, B. T., & Lee, H. M. 1984, *ApJ*, **285**, 89

Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D, Astrophysics Source Code Library, ascl:1202.015

Evans, M. G., Ilee, J. D., Boley, A. C., et al. 2015, *MNRAS*, **453**, 1147

Evans, N. J., II, Dunham, M. M., Jørgensen, J. K., et al. 2009, *ApJS*, **181**, 321

Fehér, O., Kóspál, Á., Ábrahám, P., Hogerheijde, M. R., & Brinch, C. 2017, *A&A*, **607**, A39

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

Frasca, A., Biazzo, K., Alcalá, J. M., et al. 2017, *A&A*, **602**, A33

Ghez, A. M., McCarthy, D. W., Patience, J. L., & Beck, T. L. 1997, *ApJ*, **481**, 378

Goodman, J., & Weare, J. 2010, *Comm. Appl. Math. Comp. Sci.*, 5, 65

Goto, M., Regály, Z., Dullemond, C. P., et al. 2011, *ApJ*, **728**, 5

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

Hartmann, L. 2008, Accretion Processes in Star Formation (Cambridge: Cambridge Univ. Press)

Hartmann, L., Herczeg, G., & Calvet, N. 2016, *ARA&A*, **54**, 135

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

Herbig, G. H. 1950, *PASP*, **62**, 211

Herbig, G. H. 1966, *VA*, **8**, 109

Herbig, G. H. 1977, *ApJ*, **217**, 693
Herbig, G. H. 1989, in Proc. ESO Workshop on Low Mass Star Formation and Pre-main Sequence Events, Vol. 33, ed. B. Reipurth (Garching: ESO), 233

Herbig, G. H. 2008, AJ, 135, 637

Hubbard, A. 2017, MNRAS, 465, 1910

Ilee, J. D., Boley, A. C., Caselli, P., et al. 2011, MNRAS, 417, 2950

Jones, A. F. A. L. 2008, CBET, 1217, 1

Juhász, A., Dullemond, C. P., van Boekel, R., et al. 2012, ApJ, 744, 118

Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869

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

Kóspál, Á., Ábrahám, P., Csengeri, T., et al. 2016b, ApJL, 821, L4

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

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

Kóspál, Á., Ábrahám, P., Goto, M., et al. 2011, ApJ, 736, 72

Kóspál, Á., Mohler-Fischer, M., Sicilia-Aguilar, A., et al. 2014, A&A, 561, A61

Li, A., & Greenberg, J. M. 1997, A&A, 323, 566

Liu, H. B., Galván-Madrid, R., Vorobyov, E. I., et al. 2016, ApJL, 816, L29

Liu, H. B., Vorobyov, E. I., Dong, R., et al. 2017, A&A, 602, A19

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

Lombardi, M., Lada, C. J., & Alves, J. 2008, A&A, 480, 785

Mathews, G. S., Pinte, C., Duchêne, G., Williams, J. P., & Ménard, F. 2013, A&A, 558, A66

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, adass XVI, 376, 127

Oya, Y., Sakai, N., López-Sepulcre, A., et al. 2016, ApJ, 824, 88

Perez, S., Casassus, S., Ménard, F., et al. 2015, ApJ, 798, 85

Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJL, 811, L5

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., & Laibe, G. 2014, A&A, 565, A129

Pinte, C., Ménard, F., Duchêne, G., et al. 2018, A&A, 609, A47

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

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

Rab, C., Elbakyan, V., Vorobyov, E., et al. 2017, A&A, 604, A15

Gras-Velázquez, À., & Ray, T. P. 2005, A&A, 443, 541

Rosenfeld, K. A., Chiang, E., & Andrews, S. M. 2014, ApJ, 782, 62

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

Schieven, G. (ed.) 2017, Observing with ALMA A Primer, ALMA Doc. 5.1, ver. 2

Sicilia-Aguilar, A., Fang, M., Roccatagliata, V., et al. 2015, A&A, 580, A82

Sicilia-Aguilar, A., Kóspál, Á., Setiawan, J., et al. 2012, A&A, 544, A93

Sipos, N., Ábrahám, P., Acosta-Pulido, J., et al. 2009, A&A, 507, 881

Takakuwa, S., Saigo, K., Matsumoto, T., et al. 2017, ApJ, 837, 86

Ubach, C., Maddison, S. T., Wright, C. M., et al. 2012, MNRAS, 425, 3137

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

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Yen, H.-W., Takakuwa, S., Chu, Y.-H., et al. 2017, A&A, 608, A134

Yen, H.-W., Takakuwa, S., Ohashi, N., et al. 2014, ApJ, 793, 1

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050

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

Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110

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