© ESO 2018

1 Introduction

Protoplanetary disks are considered to be the sites of planet formation (e.g., Williams & Cieza 2011). Statistical studies on properties and evolution of protoplanetary disks are essential to shed light on the diverse properties of exoplanetary systems (e.g., Winn & Fabrycky 2015). With the unprecedented sensitivity of the Atacama Large Millimeter/submillimeter Array (ALMA), several tens of protoplanetary disks can be imaged in continuum at sub-arcsecond resolutions and at high signal-to-noise ratios (S/Ns) in a few hours, providing large samples for statistical studies. Such ALMA surveys have been conducted toward a few star-forming regions, including Upper Scorpius (Carpenter et al. 2014; Barenfeld et al. 2016), Lupus (Ansdell et al. 2016), Chamaeleon I (Pascucci et al. 2016), and σ Orionis (Ansdell et al. 2017). With these large samples, the relations between the disk dust mass and the stellar mass in the individual star-forming regions have been revealed, and there is a trend that the relation becomes steeper with age. This can be explained with grain growth, drift, and fragmentation (Pascucci et al. 2016; Ansdell et al. 2017). Nevertheless, the disk dust mass estimated from (sub-)millimeter continuum emission can be underestimated by a factor of a few due to uncertainties in the optical depth of the continuum emission and the maximum grain size (Dunham et al. 2014; Tsukamoto et al. 2017). The high-resolution and high-sensitivity continuum data also allow detailed analyses of physical structures of disks, and disks in older star-formingregions tend to be less massive and larger than those in younger star-forming regions (Tazzari et al. 2017). In addition, in synergy with the ALMA and VLT/X-shooter surveys, the disk dust mass and the mass accretion rate estimated from the X-shooter spectra in the Lupus and Chamaeleon I regions are found to be correlated (Manara et al. 2016; Mulders et al. 2017). The overall properties of dusty disk populations appear to be consistent with some general features predicted by viscous accretion disk models under the condition that the viscous timescale is of the order of ~1 Myr (Lodato et al. 2017; Mulders et al. 2017; Rosotti et al. 2017). Detailed investigations of the gaseous disk properties are essential to verify these initial findings.

To investigate gas components in disks, CO isotopolog lines are also observed simultaneously in these ALMA surveys, but less than one third of the disks are detected in the CO isotopolog lines (Ansdell et al. 2016, 2017; Long et al. 2017). The total fluxes of the CO isotopolog lines are measured and compared with the grid of physical–chemical models of protoplanetary disks to estimate thedisk gas mass (Williams & Best 2014; Miotello et al. 2016). A correlation between the disk gas mass and the stellar mass is also observed, but the correlation is less significant compared to that between the disk dust mass and the stellar mass due to the smaller number of the detections in the CO isotopolog lines and large uncertainties in deriving the total gas mass from line fluxes (Ansdell et al. 2016; Miotello et al. 2017; Long et al. 2017). On the other hand, no correlation was found between thedisk gas mass and the mass accretion rate, in contrast to the expectation from viscously evolving disks (Manara et al. 2016). In addition, the gas-to-dust mass ratios derived with the disk gas mass estimated from the CO isotopolog lines are typically between 1 and 10 in these surveys (Ansdell et al. 2016; Miotello et al. 2017; Long et al. 2017). These results hint at the uncertainty in estimating the disk gas mass with CO isotopologs, possibly due to the carbon depletion (e.g., Miotello et al. 2017). Results of chemical models of protoplanetary disks also suggest that the disk gas mass can be underestimated by a factor of a few to two orders of magnitude with CO isotopolog lines when the conventional interstellar medium (ISM) CO abundance is adopted (Miotello et al. 2014, 2016; Yu et al. 2016, 2017; Molyarova et al. 2017). Nevertheless, a larger sample of disks detected in molecular lines is needed to study the relations between disk gas and dust masses and mass accretion rate.

With an on-source integration time of few minutes in the current ALMA surveys on protoplanetary disks, the gas distributions and velocity structures of the disks traced by the molecular lines can still not be imaged at high S/Ns, except for a few very bright disks. In most of the disks, the CO isotopolog lines are detected after integrating the emission over spatial and velocity ranges (Ansdell et al. 2016, 2017; Long et al. 2017). The central stellar mass and disk geometry of the majority of the disks in the surveys cannot be constrained from the gas kinematics because of the limited S/Ns in the velocity channel maps. Therefore, currently it is not possible to compare stellar masses dynamically determined from gas kinematics with spectroscopically determined masses in stellar evolutionary models of a large sample. Such a comparison is, nevertheless, important to understand stellar evolution of young stellar objects (e.g., Rizzuto et al. 2016; Simon et al. 2017). Applying further techniques, such as Keplerian masking, is required to better reveal distributions of molecular lines and to constrain gas kinematics (e.g., Matrà et al. 2015; Marino et al. 2016; Salinas et al. 2017).

In order to extract more information from molecular-line data of these ALMA surveys and to study disk properties traced by molecular lines, we apply in this work the velocity-aligned stacking method described in Yen et al. (2016) on the ¹³CO (3–2) and C¹⁸O (3–2) data of the Lupus survey (Ansdell et al. 2016). With this method, we enhance the S/Ns of the data and obtain better detections of the disks in the ¹³CO and C¹⁸O lines. We can then constrain disk orientation and stellar mass and re-estimate the disk gas mass with new measurements of the line flux. This paper is organized as follows: Sect. 2 introduces the observations and the data. Section 3 describes and demonstrates our method to measure disk parameters in the ¹³CO (3–2) line. Section 4 presents our measurements of disk orientation, disk gas mass, and stellar mass. Section 5 discusses the relations between dynamical and spectroscopic stellar masses and between the disk gas and dust masses and the mass accretion rate found in our results.

2 Observations

The data analyzed here were retrieved from the ALMA archive (project code: 2013.1.00220.S). In this ALMA project, 88 young stellar objects (YSOs) with Class II or flat infrared spectra in the Lupus star-forming region were observed with an on-source integration time of approximately 1-min per source. 0.9 mm continuum, ¹³CO (3–2; 330.587965 GHz), and C¹⁸O (3–2; 329.330552 GHz) were observed simultaneously. The details of the observations have been described in Ansdell et al. (2016). In this work, we re-analyzed the ¹³CO and C¹⁸O (3–2) data. The raw visibility data were calibrated using the standard reduction script for the cycle-2 data, which uses tasks in Common Astronomy Software Applications (CASA; McMullin et al. 2007) of version 4.2.2. The image cubes of the ¹³CO and C¹⁸O (3–2) lines were generated with the briggs weighting with a robust parameter of +0.5 and cleaned with the CASA task “clean” at a spectral resolution of 122 kHz, corresponding to a velocity resolution of 0.11 km s⁻¹. The typical angular resolution achieved is ~0.′′34, and the typical noise levels per channel are 70 and 67 mJy beam⁻¹ for the ¹³CO (3–2) and C¹⁸O (3–2) emission, respectively.

3 Method

We applied the velocity-aligned stacking method described in Yen et al. (2016) on the ALMA ¹³CO and C¹⁸O (3–2) data of the Lupus sample. This method aligns spectra at different positions in a protoplanetary disk by shifting them by the projected Keplerian velocities at their positions and then stacks them. With this alignment, the signals at the different positions are coherently added, and the total flux is accumulated in a narrower velocity range. As a result, the S/N of an aligned stacked spectrum is enhanced. We adopted this method to enhance the S/Ns of the ALMA ¹³CO and C¹⁸O data and to estimate the stellar mass of the targets and the orientations of their associated disks.

With a known stellar position and the assumption that its disk traced by the ¹³CO and C¹⁸O (3–2) lines is geometrically thin, the velocity pattern of Keplerian rotation of this disk can be described with four parameters: stellar mass (M_⋆), position angle of the major axis (PA), inclination angle (i), and systemic velocity (V_sys). A detailed formulation is given in Yen et al. (2016). In this work, for each source, we generated a series of aligned stacked spectra with different combinations of these four parameters. Each aligned stacked spectrum was generated by integrating emission over the entire disk area (i.e., azimuthally over 2π and radially from the center to the disk outer radius). Subsequently, we searched for the parameter set that maximized the S/Ns of the aligned stacked spectra, and adopted this parameter set as our measurements of M_⋆, PA, i, and V_sys of that source. In other words, we measured M_⋆, PA, i, and V_sys of each source by maximizing the auto-correlation between the data and the various generated Keplerian rotational patterns.

In our analysis, the stellar position of each source is adopted to be the center of its 0.9 mm continuum emission observed in the same observations or the pointing center of the observations if the continuum is not detected. The coordinates were obtained from Ansdell et al. (2016). The distance is adopted to be 200 pc for the sources in the Lupus III cloud and to be 150 pc for the sources in the Lupus I, II, and IV clouds (Comerón 2008). The parameter ranges that we searched to maximize the S/Ns of the aligned stacked spectra are 0.1 M_⊙≤ M_⋆ ≤ 3 M_⊙ in steps of 0.1 M_⊙ or an increase by 10% when M_⋆ ≥ 1.5 M_⊙, 0° ≤PA ≤ 355° in steps of 5°, 5° ≤ i ≤ 85° in steps of 5°, and 1.7 km s⁻¹ ≤ V_sys ≤ 6.3 km s⁻¹ in steps of 0.1 km s⁻¹. The position angle of the major axis is defined as the direction from the blue- to red-shifted emission. From the stacked “moment 0” map of the ¹³CO emission of all the sources detected in the ¹³CO (3–2) line made by Ansdell et al. (2016), we found that the ¹³CO emission is primarily within a radius of 1′′. Therefore, in our parameter search, all the aligned stacked spectra were generated by integrating the area within a de-projected radius of 1′′. Our integrated area varies with i, as the projected area of a circular disk on the plane of the sky changes.

As demonstrated in Yen et al. (2016), an aligned stacked spectrum is centered at V_sys and is approximately symmetric with respect to V_sys, when data are properly aligned. Therefore, for each aligned stacked spectrum, we fitted a Gaussian line profile. When the fitted line center is close to V_sys within the 1σ Gaussian line width, we computed the S/N of its integrated intensity within the 1σ Gaussian line width. The integration was weighted by the fitted Gaussian line shape, ∫ I_ν × G_νdν ∫ G_νdν, where I_ν is the observed spectrum and G_ν is the fitted Gaussian profile with the peak scaled to be one, such that channels closer to V_sys have larger weights. For a Gaussian-like spectrum, its weighted integrated intensity is proportional to its unweighted integrated intensity divided by its line width. Therefore, for aligned stacked spectra with the same total flux, the one that is more symmetric and has a narrower line width has a higher S/N than the others. In addition, as shown in Yen et al. (2016), the noise of aligned stacked spectra changes with different alignments even if the integrated area is the same. This is caused by de-correlation of pixels within one synthesized beam due to the alignment. Therefore, when we compared the S/N ratios of stacked spectra aligned with different parameters, their noise levels are computed from the original non-aligned data cube with the same area and channel ranges adopted to generate these aligned stacked spectra. We note that the area adopted to generate aligned stacked spectra changes with i, and the integrated channel range to compute S/Ns depends on the line profiles of spectra. Better aligned spectra tend to have narrower line widths and higher integrated fluxes, leading to higher S/Ns. In this case, our comparison of the S/Ns is not biased by the effects of the de-correlation.

After the S/Ns of the weighted integrated intensity of all the stacked spectra aligned with different parameters were computed, we calculated the S/N-weighted means and dispersions (i.e., moment 1 and 2) of PA, i, and V_sys, and the number distribution of S/Ns. Three examples of the distributions of the S/Ns in the parameter space are presented in Appendix A. The aligned stacked spectra with the S/Ns of their weighted integrated intensity below 3σ are excluded, and we trimmed the outliers in the number distribution at the high-S/N end because they are found to be isolated in the parameter space and are false signals caused by shifting high-noise channels to V_sys. We then searched for the maximum S/N near the means and within the dispersions of the parameters. The parameter set resulting in the highest S/N is adopted as our measurements.

To estimate the uncertainties of our measurements for each source, we searched for its stacked spectra that are aligned with different parameters but have the line profiles consistent with the best-aligned stacked spectrum within the uncertainty. For an aligned stacked spectrum, its line profile was quantified with its peak intensity, center, and width measured from Gaussian fitting. With these three quantities, the line profile of an aligned stacked spectrum generated with different parameters was compared with the best-aligned stacked spectrum. When their peak intensities, centers, and widths are all consistent within the uncertainties, and the difference in their S/Ns of the integrated fluxes is less than $\sqrt{2}$, we claim that this spectrum is consistent with the best-aligned stacked spectrum. We adopted the ranges of the parameters resulting in the spectra consistent with the best-aligned stacked spectrum as the uncertainties of our measurements of that source.

We performed this method on the ¹³CO (3–2) data of 88 sources observed with ALMA. A detection of the ¹³CO (3–2) emission is claimed if (1) the peak intensity of the best-aligned stacked spectrum is above 4σ and (2) there are at least three channels above 3σ within the 2σ Gaussian line width. Then we applied the best parameters of the alignment on the C¹⁸O (3–2) data because the S/Ns of the C¹⁸O (3–2) data are not sufficiently high to constrain the parameters even with this method. We have tested the method on 30 synthetic ALMA images of a blank field, which were generated with the CASA simulator and have the same resolution and noise level as the observations. Our method indeed did not find detections in these synthetic images. Therefore, our criteria of detection are sufficiently stringent to rule out the possibility that our method coherently adds noise and mimics a signal by coincidence.

3.1 Demonstration with observational data

In this ALMA dataset, we found only a few sources where the ¹³CO emission is sufficiently bright to directly image its distribution and velocity structures. In Fig. 1, we demonstrate our method with the ¹³CO data of one of those bright sources, J16083070−3828268. The ¹³CO spectra at different positions in the disk around J16083070−3828268 are centeredat different velocities (Figs. 1b–d) because of different projected Keplerian velocities at these positions. Consequently, when these spectra are directly stacked together, the emission is not coherently added (Fig. 1e), resulting in a stacked spectrum with double peaks and a wide line width of ~10 km s⁻¹ (Fig. 1g). In contrast, our method aligns these spectra with different centroid velocities first before stacking them (Fig. 1f). Thus, the emission originated from different positions in the disk is coherently added, and the stacked spectrum has a single peak and a narrower line width of ~6 km s⁻¹. The peak intensity also increases by almost a factor of two, and the S/N of the stacked spectrum is enhanced.

For the bright sources in this Lupus sample, we also analyzed their data with conventional methods to test the robustness of our velocity-aligned stacking method. Below, we present three cases, RY Lup, J16083070−3828268, and Sz 83. Figures 2a–c presents the total integrated intensity (moment 0) and intensity-weighted mean velocity (moment 1) maps of the ¹³CO emission in RY Lup, J16083070−3828268, and Sz 83. The ¹³CO components in RY Lup and J16083070−3828268 are elongated along the northwest–southeast direction, and exhibit clear velocity gradients along the elongation. The ¹³CO component in Sz 83 also exhibits a clear velocity gradient but does not show any obvious elongation. We fitted a two-dimensional (2D) Gaussian function to the intensity distributions of the ¹³CO emission. The deconvolved full-width-half-maximum (FWHM) size is measured to be 0.′′ 85 ± 0.11 × 0.′′ 34 ± 0.04 with a PA of the major axis of 287° ± 6° in RY Lup. Those in Sz 83 are measured to be 0.′′ 7 ± 0.11 × 0.′′66 ± 0.04 and 116° ± 80°, respectively. The PA of the major axis in Sz 83 is poorly constrained because the aspect ratio is almost unity. The ¹³CO emission in J16083070−3828268 shows two peaks, so its intensity distribution cannot be fitted with a 2D Gaussian function. The PA of the major axis in J16083070−3828268 is measured from the axis passing through the two peaks to be 106°±7°. The size of the intensity distribution is 1.′′27 ± 0.′′ 03 × 0.′′ 52 ± 0.′′ 02 from the FWHM widths of the intensity profiles extracted along and perpendicular to the elongation.

We extracted the position–velocity (PV) diagrams along the major axes and passing through the continuum peak positions in these sources (Figs. 2d–f). The PV diagram of Sz 83 was extracted along the PA of the velocity gradient, 140°. The data cube of RY Lup was binned every four channels, and those of J16083070−3828268 and Sz 83 were binned every two channels to increase the S/Ns per channel. We measured the peak positions at different velocity channels in the PV diagrams following the method described in Yen et al. (2013). The distances from the peak positions to the stellar positions were adopted as rotational radii (R_rot), and the relative velocities with respect to V_sys at those velocity channels were adopted as rotational velocities (V_rot). We fitted Keplerian rotational profiles, $V_{\text{rot}}\propto R_{\text{rot}}{}^{-0.5}$, to these data points obtained from the PV diagrams, and there were two free parameters in the fitting, V_rot at a representative radius and V_sys. The measured Keplerian velocities are 2.2 ± 0.1 km s⁻¹ at a radius of 150 AU in RY Lup, 2.3 ± 0.1 km s⁻¹ at a radius of 200 AU in J16083070−3828268, and 0.56 ± 0.04 km s⁻¹ at a radius of 150 AU in Sz 83. On the assumption that the disks around these sources observed in the ¹³CO emission are circular and geometrically thin, the inclination angles, i, are estimated to be 63° ± 8° in RY Lup, 66° ± 6° in J16083070−3828268, and 19° ± 12° in Sz 83 from the aspect ratios of the major and minor axes. Then, with the measured Keplerian velocities and inclination angles, we derived M_⋆ of RY Lup, J16083070−3828268, and Sz 83. The estimated M_⋆, PA, i, V_sys in these three sources with the analysis described above are listed in Table 1. We note that Sz 83 is close to face on, so its estimated M_⋆, which is proportional to 1∕ sini², is very uncertain.

We also measured the M_⋆, PA, i, V_sys of these three sources with our velocity-aligned stacking method. The aligned stacked spectra with maximum S/Ns are shown in Figs. 2g–i. The measurements from the velocity-aligned stacking method are also listed in Table 1 for comparison. The distributions of the highest S/N achieved with a given pair of parameters are presented in Appendix A to demonstrate the covariance between these parameters. We found that there is a clear correlation in the S/N distribution between M_⋆ and i. Therefore, the uncertainty in i can propagate to M_⋆, resulting in poorly constrained M_⋆, as in the case of Sz 83. On the other hand, no correlation is found between other pairs of parameters, so there is no significant covariance between these parameters. This demonstration of our velocity-aligned stacking method with the observational data shows that the results from the conventional method and our new method are consistent within the uncertainties, and our new method indeed provides robust measurements of M_⋆, PA, i, and V_sys with high-S/N data.

Fig. 1 Demonstration of our method with the 13CO (3–2) data ofJ16083070−3828268. Panel a: moment 0 map of the 13CO (3–2) emission in J16083070−3828268 obtained with the ALMA observations. Panels b–d: 13CO Spectra at the same de-projected radius of 0.′′5 but at different position angles in the disk, shown as red, black, and blue dots in panel a. Panels e and f: demonstration of stacking the three spectra with and without aligning them with the Keplerian velocities at their positions in the disk. Panels g and h final stackedspectra integrated over the entire disk with and without the alignment.

3.2 Demonstration with models

To test our method on faint disks, we generated a series of synthetic images of simple model disks. The model disks were made with different intensity distributions, stellar masses, inclination angles, and total integrated fluxes. The details of the model generation are described in Appendix B, and the model parameters are listed in Table C.1. We then simulated ALMA observations with a one-minute integration time on our model image cubes. Our synthetic images have the same angular resolution and noise level as the observations. We applied our velocity-aligned stacking method on the synthetic data to measure M_⋆, PA, i, and V_sys of these model disks. The measurements from our method are listed and compared with model inputs in Table C.1. In this test, we detected 17 out of 24 model disks. The detailed results and discussion of this demonstration with the model disks are presented in Appendix C. Figure 3 presents two examples of the synthetic moment 0 and 1 maps and stacked spectra with and without alignment of our model disks. There is no clear detection in either the moment 0 maps or the stacked spectra without alignment. For comparison, we also generated synthetic spectra without simulating observational noise (Figs. 3e and f). As shown in the synthetic spectra without noise, the intensity of the model disks is below the 1σ noise level of the observations, and the intensity becomes higher than the noise level after applying the alignment (Figs. 3e and f). Therefore, after applying our method, the detections are obtained in the synthetic aligned stacked spectra. This test demonstrates that even if the intensity distribution and velocity structure of a disk are not detected in its moment 0 and 1 maps and spectrum, our method can still find the signals and measure its M_⋆, PA, i, and V_sys. The results of our test show that we are able to detect disks with a total integrated flux as low as 400 mJy kms⁻¹. The detection limit changes with M_⋆, i, and intensity profiles, and it can be even lower than 400 mJy km s⁻¹.

The measurements from the synthetic images are typically consistent with the model inputs within the 1–2σ uncertainties, with exceptions where the difference is more than 3σ, for example, model 22 and 23 (Table C.1). The inclination is more difficult to constrain because it requires more spatial information along the direction of the minor axis. The detailed comparison between the measurements from the synthetic images and the model inputs is described in Appendix C. This demonstration with the synthetic data shows that our method can provide robust measurements for disks with an integrated line flux of ≳400 mJy km s⁻¹. However, when the total flux is as low as 400 mJy km s⁻¹, M_⋆ can be overestimated because of the limitation in constraining a disk orientation (e.g., model 22 in Table C.1). As discussed in Sect. 5.1, these possibly overestimated M_⋆ can be identified from their flux densities. Nevertheless, that does not affect the robustness of the measured flux because our method conserves the flux (Yen et al. 2016). In this demonstration, on average our method detects 80% of the total fluxes of the model disks (Table C.1). There are only two cases where only 50% of the total fluxes are detected, and only one case where the measured flux is 70% higher than the input flux. Therefore, the uncertainty in the measured flux due to the limited sensitivity with this method is typically 20%, and can be a factor of two in certain circumstances.

Fig. 2 ALMA observational results of RY Lup (top row), J16083070−3828268 (middle row), and Sz 83 (bottom row). Panels a–c: moment 0 map (contour) overlaid on the moment 1 map (color; in units of km s−1) of the 13CO emission. Green lines denote the axes where the PV diagrams are extracted. Contour levels start from 3σ in steps of 3σ, where 1σ is 77 mJy km s−1 in panel a, 74 mJy km s−1 in panel b, and 44 mJy km s−1 in panel c. Panels d–f: PV diagrams of the 13CO emission along the major axes and passing through the stellar positions. Blue and red data points denote the measured peak positions in the high-velocity channels, and green curves present the fitted Keplerian rotational profiles to the data points. Contour levels start from 2σ in steps of 2σ, where 1σ is 35 mJy km s−1 in panel d and 50 mJy in panels e and f. Panels g–i: stacked 13CO spectra with (red) and without (black) alignment integrated over the disk area. Zero velocity refers to the measured Vsys. Red and black horizontal dotted lines denote ±1σ levels with and without alignment, respectively.

Fig. 3 Synthetic images and spectra of disk model 3 (top row) and 6 (bottom row). The model parameters are listed in Table C.1. Panels a and b: moment 0 maps (contour) overlaid on the moment 1 maps (color; in units of km s−1) of the model disks. The moment 0 and 1 maps were computed by integrating the velocity range of VLSR = 3.7 ± 2.3 km s−1, the same as in Ansdell et al. (2016). Contour levels start from 2σ in steps of 2σ, where 1σ is 56 mJy km s−1. Panels c and d: synthetic stacked spectra with the alignment using the measured parameters (red) overlaid on the ones without alignment (black). Both are integrated over the disk area. Zero velocity refers to the measured Vsys. Red and black horizontal dotted lines denote ±1σ levels with and without alignment, respectively. Panels e and f: as in panelsc and d but for synthetic spectra without noise. Black histograms present the stacked spectrum without alignment, and red histograms present the stacked spectrum with the alignment using the input parameters.

4 Measured stellar mass and disk properties

With the velocity-aligned stacking method described in Sect. 3, we detect 41 out of 88 disks in the ¹³CO emission. Eleven of them were not detected in the ¹³CO emission in Ansdell et al. (2016). There are six disks detected in the ¹³CO emission in Ansdell et al. (2016) which could not be identified in our analysis. Our method could not find parameters to align their data and obtain a detection at more than a 4σ level in aligned stacked spectra. We note that our criteria to claim a detection are more stringent than in Ansdell et al. (2016), where they require that the integrated intensity is larger than 3σ level. All the measured M_⋆, PA, i, and V_sys of the 41 disks are listed in Table 2, and their velocity-aligned stacked spectra are shown in Appendix D. As demonstrated in Sect. 3.2, M_⋆ can be overestimated when the integrated flux is lower than 400 mJy km s⁻¹, and those sources with possibly overestimated M_⋆ are identified in Sect. 5.1 and are labeled in Table 2.

Among our detected disks in ¹³CO, there are 23 disks which are resolved with their PA and i also measured from the continuum emission (Ansdell et al. 2016; Tazzari et al. 2017). Figures 4a and b present the comparison of PA and i measured from our method and from the continuum emission. There is a clear correlation in PA from the two different methods, except for two outliers where the difference in PA is more than 45°. i from the two methods are also correlated but have a larger scatter, as the uncertainties in i are larger. The mean difference of the disk inclination, i, measured from our method and from the continuum emission is 15°, and i for 18 out of 23 disks are consistent within 25°. This comparison shows that our method detects independently and consistently the disk orientations, and therefore, can also reveal a disk orientation in the absence of a clear direct detection in continuum or lines.

With the measured M_⋆, PA, i, and V_sys, for each source, we then generated aligned stacked spectra of the ¹³CO emission in a series of radial bins in steps of half of the synthesized beam. We fitted Gaussian line profiles to these spectra and identified the outer radius where the detection was below 2σ. We adopted that radius as the detected disk radius of that source. Subsequently, we generated the aligned stacked spectrum of the ¹³CO emission integrated over the area within the detected disk radius, and we measured the total integrated flux within twice of the FWHM line widths in the aligned stacked spectra. We applied the same parameters to align the C¹⁸O data and generated aligned stacked spectra. Except for a few bright disks, the C¹⁸O emission is still not clearly identified in the aligned stacked spectra. Nevertheless, we assume that the ¹³CO and C¹⁸O lines trace the same regions, and we adopted the same area and velocity range as the ¹³CO emission to measure the total integrated flux of the C¹⁸O emission. Although the distribution of C¹⁸O is affected by isotope-selective effects more than ¹³CO (Miotello et al. 2016), our selected spatial and velocity ranges for the integration from the ¹³CO emission are expected to fully cover the distribution of the C¹⁸O emission. We define that the C¹⁸O emission is detected when the integrated flux is above the 4σ threshold. There are 18disks detected in the C¹⁸O line, and 9 of them are new detections compared to Ansdell et al. (2016). Their velocity-aligned stacked spectra are shown in Fig. D.2. The detected disk radii and the measured integrated ¹³CO and C¹⁸O fluxes are listed in Table 2, and the 4σ upper limit of the integrated C¹⁸O flux is listed for the C¹⁸O non-detected disks.

Figure 5 compares our measured integrated fluxes of the ¹³CO and C¹⁸O emission with the fluxed reported in Ansdell et al. (2016). Our measured fluxes are tightly correlated with the measurements in Ansdell et al. (2016). Specifically, for bright disks with integrated fluxes larger than 1000 mJy km s⁻¹, the difference between the fluxes from our method and from Ansdell et al. (2016) ranges from a few percent to 30%, and the mean difference is 10%. On the other hand, for faint disks, the difference in the measured fluxes ranges from 10% to as large as a factor of five, and several more disks with integrated fluxes of 100–500 mJy km s⁻¹ are detected with our method. We note that the area and velocity ranges for integration to measure the integrated fluxes are different between our method and Ansdell et al. (2016). In Ansdell et al. (2016), for disks that are not detected in the velocity channel maps, the integrated velocity ranges were all adopted to have V_LSR of 1.4–6 km s⁻¹, and the integrated area was determined by a curve-of-growth method on moment 0 maps. In contrast, our measurements show that V_sys of the detected disks range from V_LSR of 1.8 to 5.4 km s⁻¹ (Table 2). With a typical line width of a few kilometres per second for Keplerian disks (e.g., Fig. 2), part of the emission can be located beyond V_LSR of 1.4–6 km s⁻¹ and is therefore not included in the calculation by Ansdell et al. (2016). In addition, our method is expected to exclude more velocity channels with only noise in the integration because the integrated pixels and channels are selected based on a derived Keplerian rotation.

Figure 6 presents the comparison between the measured ¹³CO and C¹⁸O line fluxes and the continuum fluxes. There is a clear correlation between the continuum and line fluxes. Disks that are bright in the continuum tend to have higher line fluxes. To estimate the disk gas mass from the ¹³CO and C¹⁸O line fluxes, we followed the method in Ansdell et al. (2016) and Miotello et al. (2017), and compared the measured ¹³CO and C¹⁸O line luminosities with the grid of physical–chemical models of protoplanetary disks.

Figures 7 and 8 present our measurements and the expected ¹³CO and C¹⁸O line luminosities of protoplanetary disks with different masses from models by Williams & Best (2014)and Miotello et al. (2016). In Williams & Best (2014), two sets of models with different C¹⁸O abundances were made, one with the typical ISM abundance (Fig. 7, left) and one with a three times lower abundance to approximate selective photodissociation of C¹⁸O (Fig. 7, right). For each set of models, there are a few measurements that do not match any model line luminosities. Nevertheless, when the two sets of models are combined, all the measured line luminosities can be reproduced with the models by Williams & Best (2014). On the other hand, Miotello et al. (2016) modeled detailed disk chemistry and solved temperature structures with radiative transfer, but the temperature ranges covered by the grid of their models are not as wide as those in Williams & Best (2014). As discussed in Miotello et al. (2016), the temperature effect and spatial variation of the C¹⁸O abundance are the primary difference between their work and Williams & Best (2014). As there are fewer models and the parameter ranges are narrower in Miotello et al. (2016), there are more measurements that cannot be reproduced with these models.

We estimated disk masses with both models. As discussed in Sect. 3, there can be a typical systematic uncertainty of 20% in the measured integrated flux. Therefore, in our estimate of disk mass, we included this 20% uncertainty when the uncertainty of the measured integrated flux due to the noise is less than 20%. For disks detected in both ¹³CO and C¹⁸O lines, we searched for the model line luminosities consistent with our measurements within the uncertainties. The estimated disk mass is then adopted to be the mean mass of those model disks weighted by the difference between the model and observed line luminosities. The uncertainty of the estimated disk mass is adopted to be the maximum and minimum disk mass of those models. Six out of the 18 disks detected in both the ¹³CO and C¹⁸O lines can be reproduced with the models by Miotello et al. (2016). For the remaining 12 disks, we did not estimate their disk masses with the models by Miotello et al. (2016). For disks only detected in the ¹³CO line but not in the C¹⁸O line, we placed the upper limit of disk mass with the models by Williams & Best (2014). We searched for the disk models with the ¹³CO line luminosities consistent with our measurements within the uncertainties and the C¹⁸O line luminosities below our measured upper limit. We then adopted the maximum disk mass among these models to be the upper limit. In addition, we also estimated their disk masses with the fitting function of disk mass versus ¹³CO line luminosity from the grid of the disk models in Miotello et al. (2017). There are different fitting functions for face-on (i = 10°) and edge-on (i = 80°) disks. We measured i of all the detected disks. Therefore, the fitting function was selected for each disk based on whether i is smaller or larger than 70°, the same as Miotello et al. (2017). For a given line luminosity, the difference in the disk gas masses estimated with the two fitting functions for face-on and edge-on disks is ~50%. In our sample, we found two disks whose gas masses were estimated with the fitting function, and that could suffer from this ambiguity in the disk inclination, Sz 72 and Sz 129. Therefore, their uncertainties in the disk gas mass are larger.

Figure 9 compares the disk masses estimated with the models by Williams & Best (2014) and Miotello et al. (2016). For disks that are detected in both lines and whose line luminosities can be reproduced with both models, the disk gas masses estimated from the two different models are consistent within the uncertainties, and their difference ranges from a factor of one to five. For disks only detected in the ¹³CO line, with the models by Williams & Best (2014), we can only place a high upper limit of 10⁻³ to 3 × 10⁻² M_⊙. In contrast, with the fitting function by Miotello et al. (2017), the disk masses of those disks are estimated to be on the order of 10⁻⁵ to 10⁻⁴ M_⊙. All the disk masses estimated with the models by Miotello et al. (2016) and Williams & Best (2014) are listed in Table 3.

Fig. 4 Comparison of position angles (upper panel) and inclination angles (lower panel) estimated from the continuum emission (horizontal axis) in the literature and from the 13CO emission (vertical axis) with our method. In the upper panel, green and red dashed lines denote an angle difference of ±45° and ±90°, respectively. In the lower panel, green and red dashed lines denote an angle difference of ±30° and ±60°, respectively.

Fig. 5 Comparison of integrated 13CO (3–2) and C18O (3–2) fluxes from our method (vertical axis) and in the literature (horizontal axis; Ansdell et al. 2016).

Fig. 6 Comparison of the integrated 13CO (blue) and C18O (red) fluxes with the continuum fluxes of our sample disks. Triangles show the 4σ upper limit for the C18O non-detected disks.

Fig. 7 13CO (horizontal axis) and C18O (vertical axis) line luminosities from the physical–chemical models of protoplanetary disks (color dots) by Williams & Best (2014) overlaid on the measured line luminosities from the velocity-aligned stacking method (data points with error bars). The C18O non-detected disks are presented with open circles. Red, orange, yellow, green, blue, purple, and gray dots denote the models with disk masses of 10−1, 3 × 10−2, 10−2, 3 × 10−3, 10−3, 3 × 10−4, and 10−4 M⊙. The left panel is for the typical ISM abundance of C18O, while the right panel is for the C18O abundance that is three times lower the typical value.

Fig. 8 13CO (horizontal axis) and C18O (verticalaxis) line luminosities from the physical–chemical models of protoplanetary disks (color dots) by Miotello et al. (2016) overlaid on the measured line luminosities from the velocity-aligned stacking method (data points with error bars). The C18O non-detected disks are presented with open circles. Red, orange, yellow, green, blue, purple, and gray dots denote the models with disk masses of 10−1, 10−2, 5 × 10−3, 10−3, 5 × 10−4, and 10−5, 10−4 M⊙.

5 Discussion

5.1 Spectroscopic and dynamical stellar mass

We estimated the stellar masses of the 41 YSOs from their gas kinematics traced by the ¹³CO emission with the velocity-aligned stacking method. The stellar masses of 37 of them have also been estimated by comparing their stellar effective temperatures and luminosities with the evolutionary models by Siess et al. (2000) as discussed by Alcalá et al. (2014, 2017). Figure 10a compares the stellar masses estimated with the dynamical information from this work and the spectroscopic information in the literature. There is a group of sources showing a correlation between dynamically and spectroscopically determined stellar masses. In 21 out of 37 (~60%) sources, the two masses are consistent within the 2σ uncertainties, or their difference is less than 50%. We note that there are several sources that have dynamical stellar masses significantly larger than their spectroscopic stellar masses by a factor of several (upper left corner in Fig. 10a). We compare the estimated stellar masses as a function of integrated ¹³CO flux (Fig. 10b). We find that high estimated stellar masses tend to be associated with lower integrated ¹³CO fluxes. S/Ns per unit area and per unit velocity (i.e., flux density) are proportional to integrated flux over disk size and line width. The projected disk area on the plane of the sky is proportional to cos i, and the line width of disk rotation is proportional to $\sin i\sqrt{M_{\star }}$. Therefore, we compute flux densities of the disks as $F{}^{13}\text{CO}/(\cos i\sin i\sqrt{M_{\star }}$) and compare with our estimated stellar masses in Fig. 10d. A sub-group with high estimated stellar mass and low flux density is clearly seen. Therefore, their stellar masses can be overestimated due to their low fluxes, as discussed in Sect. 3.2. In Fig. 10c, we exclude the sub-group with flux density less than 250 mJy arcsec⁻². As a result, the correlation between our dynamically estimated stellar masses and the spectroscopically determined stellar masses becomes clear.

The main uncertainty in our estimated stellar mass is from the constraint on the inclination angle. As discussed in Sect. 3, it requires more spatial information to better constrain an inclination angle of a disk. The inclination angles of a subsample of the disks in this Lupus survey have also been measured from the continuum emission (Fig. 4), which has much higher S/Ns to constrain disk orientations as compared to the ¹³CO data. Therefore, we have also estimated the stellar masses of this subsample with our method but adopted PA and i from the continuum results (Ansdell et al. 2016; Tazzari et al. 2017). There are therefore only two free parameters, M_⋆ and V_sys, in this analysis. The re-estimated stellar masses are typically consistent with the original values within the 1 to 2σ uncertainties because our method can trace the disk orientation reasonably well (Fig. 4), and they tend to become closer to the spectroscopic stellar masses. The estimated stellar masses with the fixed PA and i are presented and compared with their spectroscopic stellar masses in Table 4 and Fig. 11. In this subsample, the dynamical stellar masses are in very good agreement with the spectroscopic stellar masses. Thedynamical and spectroscopic stellar masses in 10 out of 16 sources are consistent within the 1σ uncertainties, and they are consistent within 2σ in all the sources. The mean difference between the dynamical and spectroscopic stellar masses in this subsample is 0.15 M_⊙.

Our results are similar to a recent report by Simon et al. (2017). They analyzed and modeled Keplerian rotation of disks around 25 pre-main sequence stars observed with ALMA at angular resolutions of 0.′′ 2 to 0.′′ 8. With the high-S/N data and detailed modeling, the achieved 1σ uncertainties in dynamical stellar mass in their studies are less than 0.1 M_⊙, and the mean uncertainty is 0.03 M_⊙. They found that in 7 out of 25 (30%) pre-main sequence stars, their dynamical stellar masses range from 1.0 to 2.3 M_⊙, and their luminosities (log L = −1.57–0.37 L_⊙) and effectivetemperatures (logT_eff = 3.51–3.64 K) are inconsistent with the expectation for stars with such masses from the stellar evolutionary models. They suggest that these sources are unresolved binary or multiple systems. For the remaining stars with masses ranging from 0.1 to 1.1 M_⊙, their estimated dynamical masses are consistent with their spectroscopic masses within the uncertainties. Although the accuracy of our estimated stellar mass is not sufficiently high to calibrate stellar evolution models due to the low S/Ns of the data, the observed trend between the dynamical and spectroscopic stellar masses of 30 sources in Fig. 10c demonstrates the robustness of the velocity-aligned stacking method to estimate stellar masses, and our estimated dynamical masses of the subsample of 16 sources, where the inclination angles are better constrained, are in good agreement with the stellar evolutionary models of the mass range from 0.1 to 2 M_⊙.

Fig. 9 Disk gas mass estimated from the 13CO and C18O line luminosities with the disk models by Williams & Best (2014) and Miotello et al. (2016), shown in the vertical and horizontal axes, respectively.

Fig. 10 Panel a: stellar masses estimated with the spectroscopic method (horizontal axis; Alcalá et al. 2014, 2017) and with the dynamical information traced by the 13CO emission from this work (vertical axis). Filled and open symbols present the disks with the integrated 13CO fluxes above and below 400 mJy km s−1, respectively.Black, green, and red dashed lines denote the regions where the differences in mass are 0%, 25%, and 50%, respectively. Panel b: stellar masses estimated with the dynamical information from this work (vertical axis) compared with measured integrated 13CO fluxes (horizontal axis). Panel c: as in panel a but with a subsample of lowest 13CO flux density discarded. Panel d: as in panel b but the horizontal axis presents 13CO flux density. A vertical dashed line in (panel d) denotes our selection criterion for the subsample shown in panel c.

Fig. 11 As inFig. 10a but fixed PA and i were adopted from the continuum results (Ansdell et al. 2016; Tazzari et al. 2017) when we estimated the stellar masses with the 13CO emission.

5.2 Gas-to-dust ratio and disk mass

Our analysis enabled more measurements of disk gas masses from ¹³CO and C¹⁸O data, especially for disks with an integrated ¹³CO flux lower than 500 mJy km s⁻¹, as compared to Ansdell et al. (2016). Figure 12 presents our disk gas mass estimated with the ¹³CO and C¹⁸O emission and gas-to-dust mass ratios compared with the disk dust mass. In Fig. 12, the gas mass of the disks detected in the ¹³CO line but not in the C¹⁸O line is estimated with the fitting function in Miotello et al. (2017). For the disks detected in both the ¹³CO and C¹⁸O lines, their disk gas masses are estimated with the grid of the disk models by Miotello et al. (2016) when their line luminosities can be reproduced with these disk models. Otherwise their disk gas masses are estimated with the grid of the disk models by Williams & Best (2014). As shown in Fig. 9, the estimated gas masses with the two models are consistent within the uncertainties, when the observed line luminosities can be reproduced by both models. The disk dust mass shown in Fig. 12 is estimated with the continuum flux in Ansdell et al. (2016) and the fitting function in Miotello et al. (2017). We note that there are a few disks with estimated gas masses of a few × 10⁻⁵ M_⊙ but very low dust masses on the order of 10⁻⁷ M_⊙, and all their derived gas-to-dust mass ratios are a factor of several hundred. On the other hand, the distribution of the gas-to-dust ratios of the disks with their dust masses higher than 10⁻⁶ M_⊙ is relatively flat and has a median of 12. When this median gas-to-dust mass ratio is applied to disks with very low dust masses, assuming that the gas-to-dust mass ratio does not strongly depend on disk dust mass (as shown below), we expect their disk gas masses traced by the ¹³CO emission to be on the order of 10⁻⁶ M_⊙. That is ten times lower than the minimum disk mass of the grid of the disk models by Miotello et al. (2016), which is 10⁻⁵ M_⊙. In addition, Fig. 4 in Miotello et al. (2017) shows that the gas mass of a disk with a ¹³CO line luminosity lower than 7 × 10⁴ Jy km s⁻¹ pc² can be potentially lower than 10⁻⁵ M_⊙, but that low mass regime was not explored. We therefore consider that the gas mass estimates of those disks with dust masses lower than 10⁻⁶ M_⊙ may not be reliable, and we exclude those disks in the following discussion.

Figure 12a shows that high disk gas masses estimated with the ¹³CO and C¹⁸O emission tend to be associated with high disk dust masses estimated with the continuum. We performed a linear fit to the logarithms of the disk gas and dust masses by using the routine of a Bayesian approach with errors of both axes written by Kelly (2007), linmix_err. The error bars of the disk masses are estimated from the mass ranges of the grid of the disk models showing consistent line and continuum luminosities with the observations within the uncertainties (Figs. 1 and 4 in Miotello et al. 2017), so they are not conventional 1σ uncertainties with normal distributions. This is different from the assumption in linmix_err, which considers that errors are normally distributed. For simplicity and to be conservative, we adopted the logarithms of the ratios of the upper limits to the measurements as 1σ uncertainties in the linear fit, and therefore the uncertainty in our linear fit is unlikely to be an underestimate. The fitting result shows that the logarithms of the disk gas and dust masses are possibly correlated with a slope of 1.1 ± 0.5 and a correlation coefficient of 0.8 ± 0.2 in the logarithmic scale. As a consequence, the logarithms of the gas-to-dust mass ratio and the dust mass do not have any significant correlation but with the possible flat value of 12, despite a large scatter. We did not include the disks which are not detected in the ¹³CO emission in the linear fit because it is not straightforward to estimate the upper limits of their line fluxes and their resulting disk gas masses without knowing their line widths, systemic velocities, and disk sizes. Nevertheless, as demonstrated in Sect. 3, we are able to detect disks with an integrated flux of a few hundred mJy km s⁻¹ at a distance of 150 pc, depending on the actual disk parameters. Therefore, the upper limits of their ¹³CO line luminosities are most likely on the order of 10⁵ Jy km s⁻¹ pc², corresponding to upper limits of their disk gas masses of a few ×10⁻⁵ M_⊙ with the fitting function in Miotello et al. (2017). On the other hand, the disk dust masses of the disks undetected in ¹³CO are mostly less than 10⁻⁵ M_⊙. Only one undetected disk has a dust mass of more than 10⁻⁴ M_⊙, and three have a mass of more than 10⁻⁵ M_⊙. Therefore, when we include these disks and adopt an upper limit of the disk gas mass of 10⁻⁴ M_⊙ in the linear fit¹, our results remain unchanged. We have also verified that our results do not depend on whether the disk dust masses are adopted from Ansdell et al. (2016) or from Miotello et al. (2017), and that the results are not affected by including an additional 10% uncertainty due to the absolute flux calibration in the analysis.

However, we note that the scatter in Fig. 12 is large with a 1σ intrinsic scatter of 0.2 ± 0.2 estimated from linmix_err. The estimated gas masses and the derived gas-to-dust mass ratios span over three to four orders of magnitude at a given bin of the disk dust masses. We subtracted the derived correlation from the distribution of the disk gas and dust masses and measured the standard deviation of the residual gas masses to be 0.9 dex. Such a large scatter and wide range of gas-to-dust mass ratios likely suggests that additional important physics and chemistry are not yet captured in the current disk mass estimates, leading to the inaccuracy in the estimated disk gas masses. As discussed in Miotello et al. (2017), the estimates of the disk gas masses with the models by Williams & Best (2014) and Miotello et al. (2016) are based on the assumption that there is no depletion of volatile carbon and oxygen, but observations have shown hints of carbon depletion in protoplanetary disks (Schwarz et al. 2016; McClure et al. 2016), which can lead to inaccurate estimates of the disk gas masses.

In addition, we note that although the line fluxes are correlated with the continuum fluxes (Fig. 6), the presence of the possible correlation between the disk gas and dust masses is subject to the choice of disk model adopted to estimate the gas masses from the line fluxes. The possible correlation between the disk gas and dust masses is only present when we include the data points of the disk gas masses estimated with the fitting function in Miotello et al. (2017) for those disks detected in the ¹³CO line but not in the C¹⁸O line and having gas masses around a few ×10⁻⁵ M_⊙. For these disks that are undetected in the C¹⁸O line, when we estimate their disk gas masses with the disk models by Williams & Best (2014), we can only obtain upper limits of 10⁻³ M_⊙ or above (Fig. 9). When we performed the linear fit by adopting the upper limits of the gas masses estimated with the disk models by Williams & Best (2014) instead of adopting the gas masses estimated with the fitting function in Miotello et al. (2017) or by excluding this subsample, we found that the correlation between disk gas and dust masses becomes insignificant. Therefore, measuring gas masses of disks detected in ¹³CO but not in C¹⁸O is essential to constrain the possible correlation between the disk gas and dust masses. Deeper observations to detect their C¹⁸O emission are needed to compare with the disk models and examine the current gas mass estimate.

Correlations between mass accretion rates and disk dust masses have been observed in this Lupus sample (Manara et al. 2016) and in the Chamaeleon I sample (Mulders et al. 2017). The presence of a correlation between mass accretion rates and disk total masses is expected in the theory of viscously evolving disks (e.g., Hartmann et al. 1998; Jones et al. 2012; Rosotti et al. 2017; Lodato et al. 2017). However, no significant correlation between the mass accretion rate and the disk gas mass estimated from the ¹³CO and C¹⁸O lines was found in these studies. Here, with more measurements of the disk gas mass, we revisit the relation between the mass accretion rate and the disk gas mass (Fig. 13). The mass accretion rates are obtained from Alcalá et al. (2014, 2017). In Manara et al. (2016), several disks only had upper limits of gas masses ranging from 5 × 10⁻⁴ to 10⁻³ M_⊙. Here, their gas masses are re-estimatedto be on the order of 10⁻⁵ M_⊙ with the fitting function in Miotello et al. (2017). However, after adding these new measurements to the analysis, we still find no significant correlation between the mass accretion rate and the disk gas mass, regardless of whether we include the data points with only the upper limits of the disk gas mass or the mass accretion rate, similarly to the results in Manara et al. (2016) and Miotello et al. (2017), even though high disk gas masses tend to be associated with high disk dust masses (Fig. 12a). We found that even if the disk gas masses are correlated with the disk dust masses, still no correlation between the mass accretion rate and the disk gas mass is expected because of the large uncertainties in deriving the total disk gas masses from the ¹³CO and C¹⁸O emission (see Miotello et al. 2017 for a detailed discussion). We have generated simulated samples with random dust masses in the range 10⁻³–10⁻⁶ M_⊙, and we computed their disk gas masses and mass accretion rates with the observed correlations and randomly added the observed scatters (this work and Mulders et al. 2017). We subsequently compared the simulated disk gas masses and mass accretion rates, and indeed did not find any significant correlation. With the current data, disks with gas masses less than 10⁻⁵ M_⊙ cannot be probed. To investigate the correlation between the disk gas mass and the mass accretion rate with the current uncertainties in deriving the disk gas mass, measurements of disk gas masses at this low-mass end (<10⁻⁵ M_⊙) to expand the range of disk masses in the sample are needed.

Fig. 12 Panel a: disk gas mass as a function of disk dust mass. Panel b: gas-to-dust mass ratio as a function of disk dust mass. The disk dust masses are estimated with the fitting function in Miotello et al. (2017). A dashed line presents the best linear fit. Only the disks with their dust masses larger than 10−6 M⊙ were included in the linear fit.

Fig. 13 Mass accretion rate as a function of disk gas mass. The mass accretion rates are obtained from Alcalá et al. (2014, 2017).

6 Summary

In order to extract more information from molecular-line data of the ALMA surveys on protoplanetary disks and to study disk properties in molecular lines, we re-analyzed the ALMA ¹³CO (3–2) and C¹⁸O (3–2) data of 88 Lupus YSOs with the velocity-aligned stacking method. With this method, we search for detection and estimate the dynamical stellar masses of the YSOs and the position angle of the major axis, inclination angle, and systemic velocity of their associated disks with the ¹³CO data by maximizing the auto-correlation between the data and the various generated Keplerian rotational patterns. After applying this method, the S/Ns of these data are enhanced, and we obtain ¹³CO detections in 41 disks and C¹⁸O detections in 18 disks. Eleven of the 41 disks detected in ¹³CO were not detected in the previous study, and there are six disks detected in the ¹³CO emission in the previous study but not in our analysis (Ansdell et al. 2016). Our main results are summarized below.

We have tested the velocity-aligned stacking method to estimate stellar mass and disk orientation with observational and synthetic data. We measured the disk orientations of RY Lup, J16083070−3828268, and Sz 83 from their moment 0 maps of the ¹³CO emission and their stellar masses and systemic velocities from the velocity structures in the PV diagrams of the ¹³CO emission. The estimates from our velocity-aligned stacking method are consistent with these measurements from the direct images within the uncertainties. Our test with synthetic data shows that with these ALMA observations, this method can detect disks and obtain robust measurements when the integrated line flux is ≳400 mJy km s⁻¹. When the integrated line flux is as lowas 400 mJy km s⁻¹, the stellar mass can be overestimated because of the limitation in constraining disk orientation. Nevertheless, the measurement of the integrated line flux remains robust. The uncertainty in the measured flux due to the limited sensitivity with this method is typically 20% and exceptionally it can be a factor of two.

The position angles of the disks estimated from our method are clearly correlated with those from the continuum emission, except for two outliers where the difference is more than 45°. The inclination angles of the disks estimated from our method and from the continuum emission are also correlated but have a larger scatter due to the larger uncertainties in the inclination angles. Therefore, our method can detect independently and consistently the disk orientations and also reveal a disk orientation in the absence of a clear direct detection in continuum or lines. Our measured fluxes are consistent with the measurements in Ansdell et al. (2016) within a few percent to 30% for bright disks with integrated fluxes larger than 1000 mJy km s⁻¹. The difference in the measured fluxes ranges from 10% to a factor of five for fainter disks. Compared to Ansdell et al. (2016), our method additionally detects several more disks with integrated fluxes of 100–500 mJy km s⁻¹.

The stellar masses of 37 of our ¹³CO detected disks have been spectroscopically determined with stellar evolutionary models. After excluding a subsample of disks with a low ¹³CO flux density, where their dynamical stellar masses can be overestimated with our method, a correlation between the dynamical and spectroscopic stellar masses of 30 sources is seen, demonstrating the robustness of the velocity-aligned stacking method to estimate stellar masses. In a subsample of 16 sources, where the inclination angles are better constrained with the continuum data, we have also estimated their stellar masses with our method but adopted PA and i from the continuum results. In these sources, our estimated dynamical masses are consistent with their spectroscopic masses within the uncertainties with a mean difference between the two masses of 0.15 M_⊙. Therefore, the dynamical masses in this subsample are in good agreement with the stellar evolutionary models of the mass range from 0.1 to 2 M_⊙.

We estimate the total gas masses of our detected disks by comparing the measured ¹³CO and C¹⁸O line luminosities with the grid of the disk models by Williams & Best (2014) and Miotello et al. (2016). We find that high disk gas masses estimated with the ¹³CO and C¹⁸O emission tend to be associated with high disk dust masses estimated with the continuum, hinting at a linear relation, and therefore the derived gas-to-dust ratios show no dependence on the disk dust mass. However, we note that the distribution of the disk gas anddust masses has a large scatter, and the derived gas-to-dust mass ratios span over three to four orders of magnitude at a given bin of the disk dust mass, likely suggesting that additional important physics and chemistry are not yet captured in the current disk mass estimates, leading to the inaccuracy in the estimated disk gas masses. The 1σ dispersion in the distribution of the disk gas and dust masses is estimated to be 0.9 dex. We find that with such a large scatter, even if the disk gas masses are correlated with the disk dust masses, it is still expected that there is no correlation between the disk gasmass and the mass accretion rate seen in this sample. With the current uncertainties in deriving the disk gasmass, deeper observations, which can detect molecular lines in disks with a gas mass lower than 10⁻⁵ M_⊙ and expand the range of disk masses in the sample, are needed to investigate the correlation between the disk gas mass and the mass accretion rate.

Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00220.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatoryis operated by ESO, AUI/NRAO and NAOJ. We thank all the ALMA staff supporting this work. HWY, CFM, AM acknowledge ESO fellowship. PMK acknowledges support from the Ministry of Science and Technology through grant MOST 104-2119-M-001-019-MY3 and from an Academia Sinica Career Development Award. This work was partly supported by by the Italian Ministero dell’Istruzione, Università e Ricerca through the grant Progetti Premiali 2012 iALMA (CUP C52I13000140001), and by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation)–Ref no. FOR 2634/1 TE 1024/1-1. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize, and by the European Union A-ERC grant 291141 CHEMPLAN.

Appendix A Distributions of S/Ns in parameter space

Figures A.1–A.3 present the distributions of the highest S/Ns of the weighted integrated intensity of the velocity aligned stacked spectra achieved with a given pair of the parameters by varying the other two parameters. The results of the analysis of the three Lupus sources, RY Lup, J16083070−3828268, and Sz 83 are presented. Red data points with error bars denote our measurements and the uncertainties. We note that in our analysis, when estimating the uncertainties of our measurements, the information of the line profiles is also taken into account (Sect. 3). The plots here only show the S/Ns of the integrated intensity. Therefore, the estimated uncertainties can be smaller than the parameter ranges delineated by the contour of the maximum S/N − 1 of the integrated intensity.

Fig. A.1 Highest S/Ns of the weighted integrated intensity of the aligned stacked spectra achieved with a given pair of the parameters by varying the other two parameters for RY Lup. Contours from higher (dark) to lower (light) levels delineate the parameter space with achieved S/Ns of the maximum S/N − 1, − 2, and − 3. Red data points with error bars show our measurements and uncertainties reported in Table 1. The uncertainties are smaller than the parameter ranges delineated by the maximum S/N − 1 contour because of the additional constraints from the line profiles (see Sect. 3).

Fig. A.2 As in Fig. A.1 but for J16083070−3828268.

Fig. A.3 As in Fig. A.1 but for Sz 83.

In the case of RY Lup, there are isolated regions showing high S/Ns; for example, when the parameter M_⋆ is lower than 0.5 M_⊙. The presence of these regions is caused by shifting high-noise channels to V_sys with these particular sets of parameters when aligning the spectra, and they are false signals. Therefore, such isolated regions showing high S/Ns are trimmed in our analysis.

The distributions of the S/Ns also show a clear correlation between the parameters i and M_⋆, where higher inclination angles tend to be associated with lower stellar masses. Therefore, the uncertainty in i can propagate to M_⋆. This is clearly seen in the case of Sz 83. The disk around Sz 83 is close to face on, and therefore its inclination angle is poorly constrained (Sect. 3.1). This leads to a larger uncertainty in the estimated stellar mass. On the other hand, there is no correlation found between other pairs in the parameter set. Therefore, they do not have any significant covariance, even for the pair of i and PA. This is different from estimating i and PA from the intensity distributions of the continuum or molecular lines, where i and PA can have a significant covariance. In particular, when the inclination is close to face on, the disk orientation cannot be constrained from the intensity distributions. In our analysis, the disk orientation PA is constrained from the direction of the velocity gradient due to disk rotation, and therefore, it does not have a significant covariance with the inclination i.

Appendix B Simple disk models

The model disks were made with two intensity distributions, uniform and Gaussian distributions, two stellar masses, 0.2 and 0.5 M_⊙, two inclination angles, 30° and 60°, and three different total integrated fluxes, 400, 600, and 800 mJy km s⁻¹. The position angles of their disk major axes were allfixed at 135°. For those uniform disks, the disk radii were adopted to be 1′′, and for those Gaussian disks, the FWHM widths were set to be 1′′. V_sys was adopted to be V_LSR of 3.7 km s⁻¹, the mean V_sys of Lupus YSOs (Ansdell et al. 2016). The distance was assumed to be 150 pc, the same as that to the Lupus I, II, and IV regions. We first used the Mirad tasks, imgen and velmodel, to generate intensity distributions and velocity patterns of model disks, respectively. From those, model image cubes were produced using Mirad task, velimage, assuming a constant Gaussian 1σ line width of 0.2 km s⁻¹. We then used CASA to simulate ALMA observations with the array configuration of C34–5 and a one-minute integration time on our model image cubes, resulting in an angular resolution of 0.′′ 35 × 0.′′ 29 and a noise level of 72 Jy Beam⁻¹ per channel, the same as the observations.

Appendix C Detailed results of demonstration with models

We applied our method on the synthetic images of the model disks to measure their M_⋆, PA, i, and V_sys. All the results are listed in Table C.1. For the detected model disks with the integrated fluxes of 600 and 800 mJy km s⁻¹, the measured M_⋆, PA, and V_sys are all consistent with the model inputs within the 1σ uncertainties, except for model 24 and 11, where the measured M_⋆ and PA are within the 2σ uncertainties, respectively. The inclination is more difficult to constrain because it requires more spatial information along the direction of the minor axis. In 9 out of 14 detected model disks with the integrated fluxes larger than 600 mJy km s⁻¹, the measured and input i are consistent within the 1σ uncertainties, and all of them are within the 3σ uncertainties, except for model 23, where the measured i significantlydeviates from the input. Three of the model disks with the integrated fluxes of 400 mJy km s⁻¹ were detected with our method (model 4, 7, and 22). Even though their integrated fluxes are low, all the measured parameters in model 4 are consistent with the inputs within the 1σ uncertainties, and in model 7, only the measured PA significantly deviates from the inputs, while the remaining parameters are all consistent with the inputs within the 1σ uncertainties. The only exception among these detected model disks is model 22. For this model, all the measured parameters are different from the inputs. In Fig. C.1, we compare the stacked spectra

Fig. C.1 Aligned stacked spectra of model 22 in Table C.1. Black and red histograms present the stacked spectra aligned with the input and measured parameters, respectively.

of model 22 aligned with the measured and input parameters. The one aligned with the measured parameters has a higher S/N, which could be caused by its particular distribution of noise by coincidence. We have verified this by performing the same analysis on model 22 with different seeds of random noise, and there are cases with the measured parameters consistent with the inputs within the 1σ uncertainties. Although in Fig. C.1 the parameters adopted to align the data are different, their measured fluxes are consistent within the uncertainties because our method conserves the flux (Yen et al. 2016).

Appendix D Aligned stacked spectra of the full sample

Figures D.1 to D.2 present the velocity-aligned stacked spectra of the ¹³CO and C¹⁸O lines of all the detected disks in comparison with the original non-aligned stacked spectra.

Fig. D.1 Stacked 13CO spectra with (red) and without (black) alignment integrated over the disk area. Zero velocity refers to the measured Vsys. Red and black horizontal dotted lines denote ±1σ levels with and without alignment, respectively.

Fig. D.2 As in Fig. D.1 but for the C18O emission.

References

All Tables

All Figures

Fig. 1 Demonstration of our method with the 13CO (3–2) data ofJ16083070−3828268. Panel a: moment 0 map of the 13CO (3–2) emission in J16083070−3828268 obtained with the ALMA observations. Panels b–d: 13CO Spectra at the same de-projected radius of 0.′′5 but at different position angles in the disk, shown as red, black, and blue dots in panel a. Panels e and f: demonstration of stacking the three spectra with and without aligning them with the Keplerian velocities at their positions in the disk. Panels g and h final stackedspectra integrated over the entire disk with and without the alignment.

In the text

Fig. 2 ALMA observational results of RY Lup (top row), J16083070−3828268 (middle row), and Sz 83 (bottom row). Panels a–c: moment 0 map (contour) overlaid on the moment 1 map (color; in units of km s−1) of the 13CO emission. Green lines denote the axes where the PV diagrams are extracted. Contour levels start from 3σ in steps of 3σ, where 1σ is 77 mJy km s−1 in panel a, 74 mJy km s−1 in panel b, and 44 mJy km s−1 in panel c. Panels d–f: PV diagrams of the 13CO emission along the major axes and passing through the stellar positions. Blue and red data points denote the measured peak positions in the high-velocity channels, and green curves present the fitted Keplerian rotational profiles to the data points. Contour levels start from 2σ in steps of 2σ, where 1σ is 35 mJy km s−1 in panel d and 50 mJy in panels e and f. Panels g–i: stacked 13CO spectra with (red) and without (black) alignment integrated over the disk area. Zero velocity refers to the measured Vsys. Red and black horizontal dotted lines denote ±1σ levels with and without alignment, respectively.

In the text

Fig. 3 Synthetic images and spectra of disk model 3 (top row) and 6 (bottom row). The model parameters are listed in Table C.1. Panels a and b: moment 0 maps (contour) overlaid on the moment 1 maps (color; in units of km s−1) of the model disks. The moment 0 and 1 maps were computed by integrating the velocity range of VLSR = 3.7 ± 2.3 km s−1, the same as in Ansdell et al. (2016). Contour levels start from 2σ in steps of 2σ, where 1σ is 56 mJy km s−1. Panels c and d: synthetic stacked spectra with the alignment using the measured parameters (red) overlaid on the ones without alignment (black). Both are integrated over the disk area. Zero velocity refers to the measured Vsys. Red and black horizontal dotted lines denote ±1σ levels with and without alignment, respectively. Panels e and f: as in panelsc and d but for synthetic spectra without noise. Black histograms present the stacked spectrum without alignment, and red histograms present the stacked spectrum with the alignment using the input parameters.

In the text

Fig. 4 Comparison of position angles (upper panel) and inclination angles (lower panel) estimated from the continuum emission (horizontal axis) in the literature and from the 13CO emission (vertical axis) with our method. In the upper panel, green and red dashed lines denote an angle difference of ±45° and ±90°, respectively. In the lower panel, green and red dashed lines denote an angle difference of ±30° and ±60°, respectively.

In the text

	Fig. 5 Comparison of integrated 13CO (3–2) and C18O (3–2) fluxes from our method (vertical axis) and in the literature (horizontal axis; Ansdell et al. 2016).
In the text

	Fig. 6 Comparison of the integrated 13CO (blue) and C18O (red) fluxes with the continuum fluxes of our sample disks. Triangles show the 4σ upper limit for the C18O non-detected disks.
In the text

Fig. 7 13CO (horizontal axis) and C18O (vertical axis) line luminosities from the physical–chemical models of protoplanetary disks (color dots) by Williams & Best (2014) overlaid on the measured line luminosities from the velocity-aligned stacking method (data points with error bars). The C18O non-detected disks are presented with open circles. Red, orange, yellow, green, blue, purple, and gray dots denote the models with disk masses of 10−1, 3 × 10−2, 10−2, 3 × 10−3, 10−3, 3 × 10−4, and 10−4 M⊙. The left panel is for the typical ISM abundance of C18O, while the right panel is for the C18O abundance that is three times lower the typical value.

In the text

Fig. 8 13CO (horizontal axis) and C18O (verticalaxis) line luminosities from the physical–chemical models of protoplanetary disks (color dots) by Miotello et al. (2016) overlaid on the measured line luminosities from the velocity-aligned stacking method (data points with error bars). The C18O non-detected disks are presented with open circles. Red, orange, yellow, green, blue, purple, and gray dots denote the models with disk masses of 10−1, 10−2, 5 × 10−3, 10−3, 5 × 10−4, and 10−5, 10−4 M⊙.

In the text

	Fig. 9 Disk gas mass estimated from the 13CO and C18O line luminosities with the disk models by Williams & Best (2014) and Miotello et al. (2016), shown in the vertical and horizontal axes, respectively.
In the text

Fig. 10 Panel a: stellar masses estimated with the spectroscopic method (horizontal axis; Alcalá et al. 2014, 2017) and with the dynamical information traced by the 13CO emission from this work (vertical axis). Filled and open symbols present the disks with the integrated 13CO fluxes above and below 400 mJy km s−1, respectively.Black, green, and red dashed lines denote the regions where the differences in mass are 0%, 25%, and 50%, respectively. Panel b: stellar masses estimated with the dynamical information from this work (vertical axis) compared with measured integrated 13CO fluxes (horizontal axis). Panel c: as in panel a but with a subsample of lowest 13CO flux density discarded. Panel d: as in panel b but the horizontal axis presents 13CO flux density. A vertical dashed line in (panel d) denotes our selection criterion for the subsample shown in panel c.

In the text

	Fig. 11 As inFig. 10a but fixed PA and i were adopted from the continuum results (Ansdell et al. 2016; Tazzari et al. 2017) when we estimated the stellar masses with the 13CO emission.
In the text

	Fig. 12 Panel a: disk gas mass as a function of disk dust mass. Panel b: gas-to-dust mass ratio as a function of disk dust mass. The disk dust masses are estimated with the fitting function in Miotello et al. (2017). A dashed line presents the best linear fit. Only the disks with their dust masses larger than 10−6 M⊙ were included in the linear fit.
In the text

	Fig. 13 Mass accretion rate as a function of disk gas mass. The mass accretion rates are obtained from Alcalá et al. (2014, 2017).
In the text

Fig. A.1 Highest S/Ns of the weighted integrated intensity of the aligned stacked spectra achieved with a given pair of the parameters by varying the other two parameters for RY Lup. Contours from higher (dark) to lower (light) levels delineate the parameter space with achieved S/Ns of the maximum S/N − 1, − 2, and − 3. Red data points with error bars show our measurements and uncertainties reported in Table 1. The uncertainties are smaller than the parameter ranges delineated by the maximum S/N − 1 contour because of the additional constraints from the line profiles (see Sect. 3).

In the text

	Fig. A.2 As in Fig. A.1 but for J16083070−3828268.
In the text

	Fig. A.3 As in Fig. A.1 but for Sz 83.
In the text

	Fig. C.1 Aligned stacked spectra of model 22 in Table C.1. Black and red histograms present the stacked spectra aligned with the input and measured parameters, respectively.
In the text

	Fig. D.1 Stacked 13CO spectra with (red) and without (black) alignment integrated over the disk area. Zero velocity refers to the measured Vsys. Red and black horizontal dotted lines denote ±1σ levels with and without alignment, respectively.
In the text

	Fig. D.2 As in Fig. D.1 but for the C18O emission.
In the text