© ESO, 2017

1. Introduction

Planets form in the circumstellar disks orbiting young pre-main-sequence stars. In the last decade the number of known exoplanetary systems has increased exponentially, uncovering a large diversity in their architectures as well as in the physical properties – mass, size and average density – of the single planets (Lissauer et al. 2014; Winn & Fabrycky 2015). The number and location of planets that can form around a star strongly depend not only on the total mass of the parent protoplanetary disk but also on how it is spatially distributed (Mordasini et al. 2012; Alibert et al. 2011).

In the core accretion theory (Safronov 1972), planet formation starts with the growth of the sub-micron interstellar dust grains that initially populate the disk to mm/cm-sized pebbles via pair-wise collisions and subsequently to km-sized planetesimals via gravitational interaction (Testi et al. 2014; Birnstiel et al. 2016, as reviews). The population of planetesimals is then assembled into a rocky core that rapidly accrete gaseous material. The efficiency of the formation of planetesimals is affected by the local conditions within the disk and is tightly linked to the available mass in solids (Chiang & Youdin 2010, and references therein). The characterization of how solids are spatially distributed in protoplanetary disks thus provides an excellent probe of the initial conditions of the planet formation process.

Optical and near- to mid-infrared observations have been used to study the dust emission in protoplanetary disks (Bouwman et al. 2001; van Boekel et al. 2003; Juhász et al. 2010; Miotello et al. 2012) but they effectively trace only the emission of the micron-sized grains in the upper layers of the disk structure. In these regions the dust emission is largely optically thick and therefore provides us with a measurement of the disk surface temperature rather than of the total dust mass. Measurement of the disk mass in solids can be obtained from millimeter and sub-millimeter continuum observations, which are sensitive to the thermal emission of mm-sized grains located in the disk midplane. At these wavelengths the continuum emission is optically thin and it can thus be used – given an assumption on the dust opacity and temperature – to infer the dust mass (Beckwith et al. 1990; Beckwith & Sargent 1991).

In the last two decades, the development of sub-mm/mm interferometers provided us with the first spatially resolved images of protoplanetary disks which have typical angular sizes of 1–2′′ at the distance of the nearby star forming regions (150–200 pc). Using sub-mm/mm resolved observations and a proper modeling of temperature and opacity it is possible to infer the radial profile of the dust surface density (Williams & Cieza 2011, and references therein). The radial profile generally used is a power law Σ(R) = Σ₀(R/R₀)^{− γ}, which can be either truncated at a radius R_out or – following the arguments of accretion theory (Lynden-Bell & Pringle 1974; Hartmann et al. 1998) – exponentially tapered with an e-folding radius R_c. The temperature profile is usually parametrized as a power-law or in some cases derived self-consistently with a simplified physical model for irradiated disks (Isella et al. 2009). Simultaneous fit of photometric infrared fluxes is also used to provide additional constraint on the vertical structure of the disk (Andrews et al. 2009).

The first sub-arcsecond constraints on the dust surface density profiles obtained for disks in the Taurus-Auriga (Isella et al. 2009; Guilloteau et al. 2011) and Ophiuchus (Andrews et al. 2009, 2010) clouds relied on interferometric observations at 870 μm and 1.3 mm and indicated that disks are described well from an exponentially tapered profile and generally appear to have flat interiors (γ ~ 1) and sharp outer edges. These studies provided the first constraints on the dust distribution on scales of 30–40 au, but are limited in terms of sample size (10–15 objects) and result difficult to compare due to the different modeling techniques. The determination of a profile able to describe the overall dust surface density in protoplanetary disks still needs an homogeneous analysis of larger samples of disks and is now actively investigated.

In the last years, the advent of the Atacama Large Millimeter and sub-millimeter Array (ALMA) marked the beginning of a new era for protoplanetary disk studies. On the one hand, thanks to its unprecedented angular resolution, it is revealing structures in disks at ~au scales that provide excellent testbeds for our understanding of the physical processes involved in planet formation (ALMA Partnership et al. 2015; Pérez et al. 2016; Andrews et al. 2016; Isella et al. 2016). On the other hand, its increased collecting area allows us to target large samples of disks at high resolution and sensitivity in a modest amount of time. The recent surveys at 890 μm in the Lupus (Ansdell et al. 2016), Upper Scorpius (Barenfeld et al. 2016) and Chamaeleon I (Pascucci et al. 2016) star forming regions targeted hundreds of disks with comparable sensitivity and resolution (0.3–0.5′′, 15–25 au in radius). These observational data sets constitute an exceptionally large and homogeneous sample that can be used to infer fundamental physical properties of protoplanetary disks.

In this work we focus on the determination of the physical properties of the disks in the Lupus star forming region, analyzing the ALMA observations presented in Ansdell et al. (2016), which is a near complete (96% completeness) survey of Class II disks in the Lupus I–IV clouds. By modeling the continuum emission with a physical disk model for passively irradiated disks (Chiang & Youdin 2010; Dullemond et al. 2001) we determine the surface density profile of more than 20 disks and we put constraints on their midplane temperature. We perform the fit of the interferometric visibilities directly in the uv-plane and we notice that the surface brightness profiles of most disks can be described remarkably well with a smooth exponentially tapered surface density profile, with the exception of two disks (IM Lup and Sz 98) that exhibit significant ring-like excesses (van Terwisga et al., in prep.).

By comparing the inferred properties of disks in star forming regions of different mean age, it is possible to trace their time evolution. In particular, we are interested at the time evolution of the distribution of solids, which in first instance can be characterized in simple terms of radial extent and total mass. In this work we use the inferred surface density cut-off radius R_c and the 890 μm integrated flux as proxies for these two quantities. By comparing the mass-size correlation that we obtain for the Lupus disks with that obtained from literature observations of disks in the Taurus-Auriga and Ophiuchus clouds, we find a significant difference that might be tentatively explained with viscous disk spreading. The present measurements do not allow us to rule out other mechanisms (e.g., radial drift) that could possibly explain such a discrepancy, and a complete survey of the Taurus-Auriga and Ophiuchus disks sample (at present probably observed in its high-mass end) is needed.

The paper is organized as follows. In Sect. 2 we present the sample selection, in Sect. 3 the modeling details and in Sect. 4 the results, with a comparison with previous results by Ansdell et al. (2016). The results of the analysis are discussed in Sect. 5 and in Sect. 6 we draw the conclusions. In Appendix A we report all the fit results for all the disks and in Appendix B the detailed results of Bayesian linear regressions.

2. Observations and sample selection

In this paper we analyze a sub-sample of the ALMA Cycle 2 continuum observations at 890 μm of the Lupus disks presented by Ansdell et al. (2016) (id: ADS/JAO.ALMA#2013.1.00220.S).

The Lupus disks in the Ansdell et al. (2016) sample have been observed with the an array configuration covering baselines between 21.4 and 783.5 m, corresponding to an average beam size (for the continuum) of 0.34′′ × 0.28′′, that is ~ 50 × 40 au at 150 pc. The bandwidth-weighted average frequency of the continuum observations is 335.8 GHz (890 μm) and the average rms between 0.25 and 0.41 mJy beam^-1. The flux calibration error of these observations is estimated to be 10%. For additional details on the observational setup and the data reduction we refer to Ansdell et al. (2016).

The observations in Ansdell et al. (2016) targeted a nearly complete (96% completeness) sample of the sources with Class II or Flat IR spectra in the Lupus star forming complex (I to IV clouds), for a total of 89 sources out of 93, detecting 61 of them in the continuum with >3σ significance. The sub-mm observations are complemented by the VLT/X-shooter spectroscopic survey by Alcalá et al. (2014, 2017) which derive fundamental stellar parameters for all the Class II sources of the region.

The sub-sample considered in this study has been selected from the total sample of 61 detected sources (Ansdell et al. 2016) by excluding: the edge-on disks (J16070854-3914075, Sz 118), the disks with clearly resolved gaps or holes (J16083070-3828268, RY Lup, Sz 111), the resolved and unresolved binaries (V856 Sco, Sz 74, Sz 123A), the sources for which we have no information on the stellar parameters (J15450634-3417378, J16011549-4152351) and two sources that exhibit an elongated irregular continuum emission (J16090141-3925119, J16070384-3911113, which could potentially be partially resolved binaries, although there is no spectroscopic confirmation, see Fig. 2 in Ansdell et al. 2016). We also exclude 14 sources whose Band 7 observations exhibit a very low signal-to-noise ratio at large uv-distances which makes the deprojected visibility profile compatible with being unresolved: these sources result to have an integrated flux F_cont< 4 mJy (Sz 106, J16002612-4153553, J16000060-4221567, J16085529-3848481, J16084940-3905393, V1192 Sco, Sz 104, Sz 112, J16073773-3921388, J16080017-3902595, J16085373-3914367, J16075475-3915446, J16092697-3836269, J16134410-3736462). This results in a sub-sample of 35 sources. In addition to these sources, we analyze ALMA 890 μm observations of Sz 82 (IM Lup) that have been taken by another observing program (PI: Cleeves, I.; id: ADS/JAO.ALMA#2013.1.00226.S). with comparable resolution and rms noise. After calibrating the IM Lup data set with the script provided by the ALMA Observatory, we performed self calibration with two rounds of phase calibration and one round of amplitude and phase.

Table 1 summarizes the properties of the 36 sources analyzed in this work, including their stellar properties and integrated 890 μm flux and rms. In Fig 1 we highlight the properties of the sub-sample in comparison to the complete sample from Ansdell et al. (2016). In the top panel, we show the distribution of stellar masses, ranging between 0.1 M_⊙ and 3 M_⊙: it is noteworthy that the sub-sample selected for this analysis includes all the stars in the 0.7–1M_⊙ mass bin except J16090141-3925119 that has an irregular shape: in this and in future plots the sources in this mass bin are identified with circled dots. In the bottom panel, we present the integrated continuum flux at 890 μm as a function of stellar mass, differentiating the sources whose analysis is presented in this paper (blue symbols) from the sources that we excluded (black symbols) according to the criteria listed above. The 14 red dots represent disks that were initially part of the 36 analyzed objects, but resulted in having a signal-to-noise ratio at long baselines too small to allow for a robust estimate of their disk structure, compatible with being unresolved. An exception among the red dots is J16102955-3922144 (marked in blue) on whose structure we have been able to obtain a marginal constraint.

3. Modeling

To study the structure of the disks we fit their continuum emission with a disk model that is based on the two-layer approximation (Chiang & Goldreich 1997; Dullemond et al. 2001). In the following we introduce the fundamental quantities that are needed to compute the disk emission and for more details we refer to Tazzari et al. (2016).

3.1. Disk model

Under the basic assumptions that, at each radius, the disk is vertically isothermal and in hydrostatic equilibrium, the two-layer approximation allows us to compute the disk continuum emission given the properties of the central star, a surface density profile and a dust grain size distribution. At mm wavelengths, most of the disk emission originates from the dust grains residing in the dense and geometrically thin disk midplane. Based on the conservation of the energy delivered by the star onto the disk surface and on its propagation to the disk midplane, the two-layer approximation is thus appropriate for reproducing the disk emission at mm wavelengths (Isella et al. 2009; Ricci et al. 2010; Trotta et al. 2013; Testi et al. 2016).

Following previous studies (Andrews et al. 2009; Trotta et al. 2013; Tazzari et al. 2016), we parametrize the gas surface density with a self-similar solution for an accretion disk (Lynden-Bell & Pringle 1974; Hartmann et al. 1998), using the parametrization: (1)where Σ₀ is a normalization factor, R₀ is a characteristic radius that we keep fixed to 10 au, γ is the power-law slope and R_c is the exponential cut-off radius. The dust surface density is given by: (2)where ζ is the dust-to-gas mass ratio, assumed to be constant and equal to the typical ISM value ζ = 0.01. The choice of the profile in Eq. (1)and of a constant dust-to-gas ratio are a clear simplification of reality, in which we expect ζ to change across the disk from both observational (de Gregorio-Monsalvo et al. 2013) and theoretical (Birnstiel & Andrews 2014) arguments. However, since in this study we are only analyzing the dust continuum emission, we cannot pose any constraint on the actual gas-to-dust variations in the disks. These choices are useful as they provide us with a simple parametrization of the dust distribution that we can directly compare to other studies. As a result, there are three free parameters describing the surface density: Σ₀, γ and R_c.

To compute the continuum emission we compute the dust opacity of the grain population using the Mie Theory, which allows us to compute the emissivity of a single spherical grain, and the Bruggeman mixing theory (Bruggeman 1935), which allows us to compute the effective dielectric constants for composite grains. For both the disk surface and midplane we use a MRN-like grain size distribution (Mathis et al. 1977), characterized by a number density n(a) ∝ a^{− q} for a_min ≤ a ≤ a_max and n(a) = 0 otherwise, where a is the grain radius. In the surface layer we use a_min = 10 nm and a_max = 1 μm with q = 3.5, typical of a population of small grains. This ensures that the surface layer is optically thick to the stellar radiation. In the disk midplane, where dust coagulation and settling are expected to occur (Dullemond & Dominik 2004), we use a_min = 10 nm and a_max = 1.023 cm with q = 3.0, which corresponds to a population of larger grains. The particular value of a_max has been chosen for continuity with previous analysis and reproduces the same opacity (per gram of dust) used by Ansdell et al. (2016)κ_{890 μm} = 3.37 cm² g^-1. Similarly to Trotta et al. (2013) and Tazzari et al. (2016), we assume spherical grains and we adopt the simplified volume fractional abundances found by Pollack et al. (1994): 20.6% carbonaceous materials, 5.4% astronomical silicates, 44% water ice and 30% vacuum, for an average grain density of 0.9 g cm^-3.

Finally, the disk appearance on sky is set by the disk inclination along the line of sight, defined as i = 0° for a face-on disk and i = 90° for an edge-on disk, and by the disk Position Angle, defined east-of-north from PA = 0° to PA = 180°.

3.2. Disk flaring

Computing a realistic dust temperature profile is key for a reliable estimate of their sub-mm continuum emission and therefore of their mass. The self-consistent fully-flared models based on the two-layer approximation (Chiang & Goldreich 1997; Dullemond et al. 2001) are typically too vertically thick and do not properly reproduce the spectral energy distribution in the far-IR. This is confirmed by theoretical and observational studies that require a reduced disk flaring (i.e., some degree of dust settling) in order to reconcile the far-IR and the sub-mm fluxes (e.g., Rodgers-Lee et al. 2014; Daemgen et al. 2016). These authors use the ratio between the fluxes in the far-IR and the J-band (a good proxy for the stellar photospheric emission) to estimate the disk flaring. Indeed, spectroscopic studies of disks around very low mass T Tauri stars and brown dwarfs have shown that the dust temperature strongly depends on the disk flaring (in addition to stellar luminosity), and to a lesser degree on the disk mass and other disk properties. In this work we use the spectral slope between the far-IR and the J-band to obtain a rough estimate of the disk flaring that characterizes the disks in our sample. Since the spectroscopic measurements are not available for all the sources in our sample, we derive an average disk flaring reduction factor that we use for all the sources that we have analyzed.

To this purpose, we use Herschel/PACS measurements at 70, 100 and 160 μm obtained from the Herschel Science Archive and performing the photometry as in van der Marel et al. (2016). We compute the spectral slope between the Herschel/PACS bands and the 2MASS J-band (1.235 μm) as (Adams et al. 1987; Daemgen et al. 2016): (3)which provides a model-independent estimate of the far-IR dust emission. In Fig. 2 we present the histograms of the computed spectral slopes.

On the same histograms the orange vertical lines represent the spectral slopes computed from our disk model for different values of the flaring reduction parameter f. By varying f, we can manually reduce the disk H_surf/R from a fully-flared profile (f = 1), which produces flatter spectral slopes (α ≥ −0.5), to progressively less flared models (f< 1), which produce progressively steeper spectral slopes (α ≤ −0.5). Given that Herschel measurements are available only for a subset of the sources in our sample and the approximate approach for deriving the flaring reduction factor, we decided to apply the same average value of f to all disks in the sample. By comparing synthetic and observed spectral slopes (Fig. 2) we find that a disk flaring reduced by a factor ≈3 (f = 0.3) gives a good representation of the slopes at all the three bands simultaneously. We adopt this value for all the fits that we perform in this study. To assess the impact of this choice on the constrained disk structure we have repeated the fits for different f values of 0.1, 0.3 and 1.0. We do not present these detailed checks here, but we observe that the average temperature (and therefore surface density normalization) change by a factor 1.5 at most with f varying between 0.1 and 1.0, while the surface density profile remains mostly unaltered with the cut-off radius R_c varying by 10% at most.

3.3. Modeling methodology

We perform the fits with a Bayesian approach, which produces probability distribution functions (PDFs) for the free parameters of the model by means of a Markov chain Monte Carlo (MCMC) algorithm. The free parameters are seven: Σ₀, γ and R_c which define the surface density, i and PA which define the disk appearance, Δα and Δδ which define the (RA, Dec) offset of the disk center with respect to the phase center of the observations. We consider these last two parameters as nuisance parameters: for each disk we fit the disk center to achieve a better matching between model and observations, however the information encoded in such offset is not relevant for the aims of this study¹.

Following the implementation developed by Tazzari et al. (2016), for a given set of values of the free parameters the model produces a synthetic image of the disk that is then Fourier-transformed and sampled in the same (u,v)-locations of the observed visibilities. We finally compute the χ² as: (4)where V^obs and V^mod are the observed and the synthetic visibilities, respectively, w_i is the weight associated with the observed visibilities at the (u_i,v_i) location and N is the total number of (u,v)-locations. As in Tazzari et al. (2016), the theoretical weights are calculated according to Wrobel & Walker (1999) and then re-scaled to ensure that ∑ _iw_i = 1. The posterior PDF is computed as exp(−χ²/ 2) within a rectangular domain in the parameter space (the region of interest), and zero outside such domain. The ranges defining the domain of exploration of the parameter space are detailed in Table 2.

The region of interest in the seven-dimensional parameter space is explored using an ensemble of Markov chains that evolve simultaneously according to the affine-invariant MCMC algorithm by (Goodman & Weare 2010). Two of the main advantages of using this algorithm are: first, several chains are initialized at random locations across the domain of interest, thus ensuring that the PDFs that we derive after they have converged do not depend on their initialization. Second, the algorithm enables a massive parallelization of the computation: in a reasonable amount of time it allows us to achieve a rather solid sampling of the posterior PDF out of which we can derive reliable PDFs for all the parameters.

In this study, we perform the fits using one hundred chains, which is a reasonable number for a seven-dimensional parameter space (approx. 10–20 chains per parameter). We initialize the chains in random positions in the domain of interest, making sure that they are not initialized too close to the borders in order to avoid computational issues. After initialization, we let the chains evolve for a burn-in phase of 1000–1500 steps (the actual number varies from source to source) and then we take 4000–5000 steps to achieve a good sampling of the posterior PDF. This results in approximately ~5−7 × 10⁵ evaluations of the posterior for the fit of each disk.

The product of a MCMC fit is the chain that results from collecting the locations of all the walkers throughout their evolution. Each element of the chain represents a sample of the posterior PDF. Therefore, to give an adequate representation of the results of the fit we always provide a plot of the chain, projected in the various dimensions. By marginalizing² the chain over one parameter we obtain an estimate of its PDF, out of which we derive its value as the median and uncertainty as the central interval (between 16% and 84% percentiles). By marginalizing the chain over two parameters of interest we obtain the 2D distribution of the samples from which we can study their correlation. To perform the fit we use the implementation of the MCMC algorithm provided by the Python package emcee (Foreman-Mackey et al. 2013), which allows us to exploit the massively parallel nature of the algorithm by running the fits on many cores simultaneously.

4. Derived disk properties

In this section we present the results of the fits. In order to present the analysis performed for each disk, here we illustrate the fit results for Sz 71, and in Appendix A we report the detailed plots for all the other sources.

The “staircase” plot in Fig. 3 shows the Markov chain resulting from the fit of Sz 71. On the main diagonal we show the histograms of the marginalized distribution, with vertical dashed lines indicating the 16th, 50th and 84th percentiles. The distributions are single-peaked, with profiles very close to Gaussians. For Sz 71 we find the following best fit parameters, corresponding to the model with maximum likelihood: γ = 0.27 ± 0.01, Σ₀ = 16.2 ± 0.25 gcm^-2, R_c = 85 ± 1 au, i = 40.8 ± 0.7°, PA = 37.5 ± 0.1°, Δα = −0.182′′ ± 0.002′′ and Δδ−0.559′′ ± 0.002′′. The red lines in Fig. 3 highlight the location of the best-fit model. The off-diagonal 2D plots represent the bivariate distributions between each pair of parameters which provide an immediate estimate of their correlation. The staircase plots for the other disks are presented in Appendix A.

A description of the physical structure of the best fit model for Sz 71 can be found in Fig. 4, where we plot the gas surface density and cumulative mass (left panel), the midplane temperature (middle panel) and the optical depth of the disk midplane at the observing wavelength (right panel) as a function of the distance from the star. In all plots, we highlight the radius containing the 95% of the mass as a vertical dotted red line. The disk model for Sz 71 has a very flat (γ ~ 0) surface density profile in the inner disk and a sharp exponential cut-off at R_c ~ 85 au, with 95% of the mass contained within 150 au. The disk is optically thin at 890 μm almost everywhere, except in the inner R < 2.6 au region which, anyway, gives a negligible contribution to the total mass. The midplane temperature decreases monotonically from 325 K in the innermost disk region (R ~ 0.1 au) to ~ 10 K at 100 au and then levels to 7 K, which is the minimum temperature allowed in the model, chosen to give a simple realization of the typical interstellar radiation field.

Figure 5 is a visual representation of how the distribution of models (i.e., the posterior samples shown in Fig. 3) compare with the observations. The top (bottom) panel shows the real (imaginary) part of the deprojected observation and model visibilities as a function of baseline length. The observation and model visibilities have been centered on the disk centroid (according to the inferred offsets Δα and Δδ) and then de-projected assuming the inferred values of i and PA. Of the model visibilities we show the posterior PDF (as the blue density indicator for each uv-bin), the median (black solid line), the 1σ central interval (the black dashed lines) and the 2σ central interval (the red dotted lines). In the case of Sz 71 the posterior PDF of the model visibilities has a very narrow peak, therefore these lines are very close to the median (cfr. with the broader model visibilities PDF derived for some disks in Appendix A). The real part of the observation and model visibilities match almost perfectly up to 300 kλ and are compatible within 2σ at higher uv-distances. The imaginary part of the observed visibilities is on average 0 (as it should be for a centered azimuthally symmetric surface brightness distribution) with a residual oscillating behavior at very low signal to noise level, probably due to some sort of asymmetry in the disk that cannot be described by our axisymmetric disk model.

In Fig. 6 we compare the observation to the best-fit model images. The three panels illustrate the images of the observations, of the model and of the residuals, derived from the respective visibilities with the CLEAN algorithm (Clark 1980) and a natural weighting scheme using the software CASA 4.5.0. The best-fit model (whose location in the parameter space is highlighted in Fig. 3) reproduces the observations extremely well (the residuals being lower than the 3σ level).

We note that our smooth model describes well all the disks (see Appendix A), with the exceptions of IM Lup and Sz 98 where we find systematic residuals. These residuals are, respectively, 14% and 4% of the total intensity and are probably related to the presence of radial inhomogeneities like rings or spirals (see, e.g., ALMA Partnership et al. 2015; Guidi et al. 2016; Isella et al. 2016; Pérez et al. 2016). These are minor discrepancies related to a very small fraction of disks in our sample and, while interesting to follow up in the future (van Terwisga et al., in prep.), do not affect the results we discuss in this paper.

In Table 3 we report the values of the parameters derived for Sz 71 and for all the other disks. For each disk, we provide estimates of the free parameters γ, Σ₀, R_c, i and PA, derived from the Markov chains as described in Sect. 3.3. For the 22 disks in the sample the Markov chains converged to single-peaked distributions with moderate to absent degeneracy. In all these cases, the best-fit model has usually a normalized chi-square of 1.0 ± 0.2 and the residuals are at very low signal-to-noise levels. From the chains, we compute some derived quantities such as the disk outer radius R_out, defined as the radius containing 95% of the model flux, and the total dust mass M_dust, computed by integrating the dust surface density: (5)where R_in = 0.1 au is the inner edge of the radial grid of the model (fixed for all the disks). The derived quantities are estimated as the median of their derived PDF, assigning an uncertainty that corresponds to the central interval between the 16th and 84th percentile.

4.1. Comparison with Ansdell et al. (2016)

We now compare our results with those by Ansdell et al. (2016), who derived dust masses and disk inclinations with a simplified method.

First, in Fig. 7 we compare the dust masses: the plot shows the ratio between the masses derived from our fits and those derived by Ansdell et al. (2016) by converting the spatially-integrated sub-mm continuum flux F_ν into dust mass: (6)where d is the distance, κ_{890 μm} = 3.37 cm² g^-1 (per gram of dust) is the dust opacity, B_ν(T) is the black-body brightness at the temperature T, and T_d = 20 K is the dust temperature.

The dust masses derived by Ansdell et al. (2016) with the simple conversion formula of Eq. (6)and a constant temperature of T_d = 20 K are accurate within a factor of 2 at a 1σ level for the majority of disks (15 out of 20).

From Fig. 7 it is also clear that our mass estimates are systematically larger than those by Ansdell et al. (2016). This discrepancy is not caused by a different dust opacity, as in our fits we assumed exactly the same opacity used by Ansdell et al. (2016). Rather, we interpret the discrepancy to originate from a different assumption on the dust temperature: while Ansdell et al. (2016) assumed a disk average temperature of 20 K for all the disks in the sample (regardless the spectral type, the mass and the luminosity of the central star), in our fits we use a physical model based on the two layer approximation (see Sect. 3.1) that takes into account the stellar properties and derives the radial profile of the midplane temperature T_mid(R) by solving the energy balance at each radius under the assumption of vertical hydrostatic equilibrium and some degree of dust settling set by the flaring parameter f. By checking the temperature profiles resulting from our fits, we find that in many cases (see Fig. 8) the disk-averaged temperature derived by our physical model (⟨ T_mid ⟩) is smaller than 20 K (but never smaller than 7 K by definition), thus explaining the tendency towards larger masses that characterizes our estimates. Moreover, since we do not observe a systematic trend of T_mid as a function of stellar mass and luminosity, we conclude that the assumption of a constant T_d (Ansdell et al. 2016) has not introduced a bias in the determination of M_d. Conversely the assumptions of Andrews et al. (2013) and van der Plas et al. (2016) on the dependence of the average T_d on the stellar parameters may introduce a spurious dependence on the derived M_d vs. M_⋆ relationship.

In the right panel of Fig. 7 we compare the inclinations derived for all disks. In all cases for which Ansdell et al. (2016) provide a measurement of inclination, their estimate is in very good agreement with ours. Moreover, in many cases we are able to put a more stringent constraint on the disk inclination, as shown by the smaller error bars. The improvement in the estimate of the disk inclination is likely due to the fit procedure, which in the case of Ansdell et al. (2016) is based on the CASA procedure uvmodelfit with the assumption of a gaussian bightness profile, while in our case benefits of a more extended exploration of the parameter space.

4.2. Distribution of surface density profiles and outer disk radii

In Fig. 9 we show the distribution of slopes γ (left panel), cut-off radii R_c (middle) and their correlation (right) obtained from our fits. The distribution of cut-off radii R_c has 6 disks with R_c< 25 au, 14 disks with 25 au ≤ R_c ≤ 100 au and two much larger disks with R_c ~ 200 au (Sz 98) and R_c ~ 430 au (IM Lup). The distribution of disk sizes (R_out, plot not shown) is similar to that of R_c: 4 disks are compatible with being unresolved (the inferred size is smaller than the spatial resolution), 14 disks have a size between 25 au and 100 au, and two largest disks Sz 98 and IM Lup have a size of ~ 280 au and ~ 500 au, respectively. Our inferred range of values of R_c, in line with the findings of Andrews et al. (2009, 2010) in Ophiuchus (14 au ≤ R_c ≤ 200 au) and Isella et al. (2009) in Taurus-Auriga (30 au ≤ Rc ≤ 230 au) who fit sub-mm observations of several bright disks using an exponentially-tapered power law surface density profile as we do here.

The distribution of γ (Fig. 9, left panel) is centrally peaked around γ = 0 and has a standard deviation of 0.6. In the right panel of Fig. 9 we show the distribution of γ as a function of R_c obtained for each disk, showing that there is not any particular trend between these two quantities. We notice that 10 disks are compatible with γ> 0 and other 10 disks with γ< 0. Interestingly, among the disks with negative γ, several of them are characterized by relatively small cut-off radii (40–70 au), which could be a signature of a partially resolved cavity.

These findings compare well with the trend that has been emerging in last years thanks to the significant improvement in the observational capabilities in the (sub-)mm window (Williams & Cieza 2011): while low angular resolution observations were usually compatible with large and positive γ values (e.g., ⟨ γ ⟩ ~ 0.9, Andrews et al. 2009, 2010; γ ~ 0.9, Hughes et al. 2008) mostly set by the fall-off of the outer disk (the only part clearly spatially resolved by such observations), higher angular resolution observations (capable of resolving the inner 100  au region) reveal that smaller γ values (⟨ γ ⟩ ~ 0.1, Isella et al. 2009; γ ~ 0.1, de Gregorio-Monsalvo et al. 2013) are required in order to reproduce the radial profile of the continuum brightness distribution consistently in the inner and outer disk.

It is worth noting that the surface density profiles have been determined assuming a dust opacity κ_λ that is constant with radius, choice that is justified by the current lack of information about the radial changes of the dust grain properties. Since a combination of grain growth and radial drift are expected to produce a size-sorting effect in the grain radial distribution (Birnstiel et al. 2010) with a consequent gradient in the dust opacity, we caution that the slope of Σ(R) derived here might be shallower than the real one (see the discussion in Banzatti et al. 2011; Trotta et al. 2013). As shown by Tazzari et al. (2016), such degeneracy can be broken by forward-modeling spatially resolved multi-wavelength observations.

4.3. Transition disks

Transition disks are protoplanetary disks that exhibit inner cavities or gaps in the distribution of their dust emission and, in some cases, also in that of the gas. In the observations by Ansdell et al. (2016) three disks already classified as TD (see references below Table 4) were detected with clearly resolved gaps in the continuum emission at ~ 50 au resolution (J16083070-3828268, RY Lup, Sz 111), three other sources (Sz 123A, Sz 100 and J16070854-3914075) showed marginal evidence for cavities with diameter of 0.4′′ and six sources (Sz 84, MY Lup, Sz 112, J16011549-4152351, J16102955-3922144, and J16081497-3857145) previously classified as TD did not exhibit any cavity or hole in the synthesized maps at 60 au resolution.

In Table 4 we report the findings of the fits of these disks. As explained in Sect. 2, we do not fit the three disks with clearly resolved gaps since the huge density depletion and the considerable size of such gaps is likely to have changed the heating of the dust in the remainder of the disk, which cannot be explained with a simple model based on the two-layer approximation. These disks will be analyzed in detail in van der Marel et al., in prep. Conversely, we apply our modeling to all the disks in the sample for which a gap is not observed in the synthesized maps: among these disks, five of them are found to have surface densities inwardly decreasing that are compatible with partially depleted gaps.

For Sz 100, one of the three sources with possible cavities, we obtain a robust fit with γ = −1.5 that confirms the presence of an inner hole with radius R_hole ≈ 46 au, where we defined R_hole as the radius where Σ_d(R) peaks. Unlike Sz 100, the other two sources with possible cavities Sz 123A and J16070854-3914075 were excluded from our sample, the former because it is a binary, the latter because it is edge-on. Finally, four out of the six disks classified as TD but with no evidence of cavities in the continuum maps were included in our sample: for Sz 84 and MY Lup we derive a surface density profile with a clear hole (γ = −1.0 and γ = −0.8) located respectively at R_hole ≈ 41 au and R_hole ≈ 34 au, for J16102955-3922144 the fit is more uncertain (uncertainty on γ is large), with a marginal evidence of a cavity at R_hole ≈ 31 au). The findings of our fits are confirmed from the fact that all these disks exhibit visibility profiles (cfr. deprojected visibility plots in Appendix A) whose real part tends toward (and in some cases reaches) negative values, compatible with the surface brightness of a disk with a central cavity. The other two disks classified as TD (Sz 112 and J16011549-4152351) were excluded from our sample due to a low integrated flux (below 4 mJy) and due to the lack of stellar parameters, respectively.

In addition to these disks, we also find evidence for the presence of holes in other two disks not classified as TD (J16000236-4222145 and Sz 129) for which we find robust estimates of negative γ values, respectively γ = −0.20 ± 0.02 and γ = −0.33 ± 0.02, and hole sizes of R_hole = 30 au and R_hole = 22 au (comparable with the spatial resolution of the observations). The surface density profiles corresponding to such γ values imply the presence of inner holes but the depletion factor inside R_hole is expected to be not as high as for γ ≤ −1. We thus conclude that for these two disks the evidence for an inner hole is tentative and to be confirmed with higher angular resolution observations. For a visual representation of the surface density profiles of these disks, see Fig. 10.

5. Discussion

5.1. Disk mass-size relation

During the evolution of a protoplanetary disk, the spatial distribution of its mass and angular momentum change dramatically. While the material gets accreted onto the star, the disk mass is a decreasing function of time, while the disk size may increase as well as remain constant, depending on the mechanism driving the angular momentum redistribution. While a viscosity-driven disk would increase in size as a result of the diffusive evolution (Lynden-Bell & Pringle 1974; Pringle 1981), if the disk angular momentum is lost through MHD winds, then partial suppression of the disk spreading can take place (Armitage et al. 2013; Bai & Stone 2013; Bai 2016; Suzuki & Inutsuka 2014). In the previous Section we have constrained the mass and the size of 22 Lupus disks and here we use such measurements to gain an insight on their evolutionary stage.

Recent (sub-)mm observations of protoplanetary disks at an angular resolution high enough to resolve their spatial structure (≲0.75′′ for nearby SFRs) seem to suggest that fainter disks are also more compact (Andrews et al. 2010; Piétu et al. 2014; Andrews 2015). Since at sub-mm wavelengths the disk emission is substantially optically thin, such trend can be interpreted in terms of a disk size – disk mass correlation. In Fig. 11 we show the 890 μm integrated flux (normalized to a common distance of 150 pc) as a function of cut-off radius for the sample of Taurus-Auriga and Ophiuchus disks (yellow points). The data have been collected from previous observations with angular resolution better than 0.75′′ from Andrews et al. (2009, 2010), Isella et al. (2010), Guilloteau et al. (2011), Piétu et al. (2014) and considers only “full” disks (binaries and transition disks are excluded). We note that the size measurements of this sample have been determined homogeneously by fitting the spatially resolved observations with an exponentially tapered power-law profile like the one in Eq. (1) and correspond to the exponential cut-off radius R_c. For some disks in Taurus-Auriga we extrapolated the flux at 890 μm from the observed 1.3 mm flux using an average spectral index α = 3.0 that corresponds to moderate grain growth. We note that this choice might have a minor effect in the inferred correlation, and future spatially resolved observations at several sub-mm/mm wavelengths will help removing this assumption.

To test the presence of a size-mass correlation in the Taurus-Auriga and Ophiuchus sample, we perform a linear regression using the Bayesian algorithm by Kelly (2007)³ which allows the uncertainties on both axes to be included in the computation. For the present case we assume uniform priors on the correlation coefficients and we take the medians of the marginalized posterior as their best-fit values (see full posterior in Appendix A). We obtain the following correlation: (7)with a correlation coefficient of 0.86 ± 0.06 and a standard deviation of 0.06 ± 0.02. Considering this sample alone, it is difficult to understand whether the observed correlation represents an evolutionary sequence or rather reflects the time evolution of disks with intrinsically different viscous timescales and initial disk conditions. Indeed, the size-mass measurements of the Taurus-Auriga and Ophiuchus disks in Fig. 11 give us a snapshot of their structure at a given moment of their evolution, and the large uncertainties in the relative ages of the stars in the sample do not allow us to distinguish between these three possible causes.

In order to shed light on the observed size-mass correlation, in Fig. 11 we add the results of this study, reporting the size measurements of 16 (out of 22 analyzed) Lupus disks. Since the Taurus-Auriga and Ophiuchus samples consider only full disks and exclude binaries and transition disks, we removed from the plotted Lupus sample the six disks with cavities described in the previous Section. In terms of stellar spectral types the Taurus-Auriga and Ophiuchus and the Lupus samples are not identical, but comparable: the former sample is mostly made of stars between K5 to M1 spectral types, while the latter extends to slightly later types, mostly between K7 to M5 types.

The resulting luminosity-size plot in Fig. 11 shows that, for a given integrated (sub-)mm flux, the Lupus disks tend to be slightly larger than the disks in Taurus-Auriga and Ophiuchus, with the exception of two disks, Sz 68 and Sz 83, which appear smaller than the average Lupus disks with a comparable flux. This evidence seems consistent throughout the range of cut-off radii between 10 and 200 au. We note that the Lupus disks in the 0.7–1.0M_⊙ mass range (circled blue dots) in which our sample is complete appear to be distributed randomly, with no signs of a particular correlation. If we apply to the Lupus sample the same linear regression used above, we find that also the Lupus disks exhibit a similar size-mass correlation: (8)with correlation coefficient of 0.73 ± 0.15 and a standard deviation of 0.12 ± 0.05. The slopes of the two linear correlations are similar within the uncertainties, both suggesting that larger disks appear brighter and smaller disks appear fainter. The Lupus correlation results systematically below the Taurus-Auriga and Ophiuchus one, confirming that the Lupus disks tend to be larger and fainter than the Taurus-Auriga and Ophiuchus ones. In this study we assumed a distance of 200 pc for the disks in the Lupus III cloud, however recent measurements of the Gaia space telescope suggest they could be closer, at a distance of ~ 150 pc, which would make them even fainter.

In order to test whether the Lupus disks are effectively less massive and more extended than the Taurus-Auriga and Ophiuchus ones, we need to check whether the two samples are distinguishable along a direction perpendicular to the inferred average correlation. To do that, we proceed as follows. First, we consider the two samples as a single sample and we perform a linear regression (same Bayesian algorithm used above), finding a correlation: (9)with correlation coefficient of 0.80 ± 0.06 and a standard deviation of 0.08 ± 0.02. As expected, the correlation results parallel to those derived for the two samples considered separately and with a vertical offset intermediate between them. Then, for each disk, we computed the distance from the fit δ_fit, which is defined to be positive for disks more massive and smaller than the linear relationship in Eq. (9), and negative for disks that are less massive and larger. In Fig. 12 (upper panels) , we show the distributions of (δ_fit) for the two samples of Lupus (blue histogram) and Taurus-Auriga and Ophiuchus (yellow histogram); the left plot is for the full sample, while the right one is for the subsample of stars between 0.7 and 1 M_⊙. In the bottom panels of Fig. 12 we show the empirical cumulative distribution functions (ECDF) of the two samples. We performed an Anderson-Darling (AD) two-sample test⁴ to check whether the two populations in Lupus and in Taurus-Auriga and Ophiuchus are consistent with being drawn from the same distribution of δ_fit. The AD test on the full sample gives a very low probability that the two samples are drawn by the same parent distribution (≤0.2%). When we restrict to the sample in the 0.7 ≤ M_⋆/M_⊙ ≤ 1 range, then the null hypothesis cannot be excluded (p ~ 8%).

5.2. Evidence for viscous evolution?

Considering that the typical disk dispersal time scale is 5–10 Myr and that the Lupus SFR is slightly older (by 1–2 Myr) than the Taurus-Auriga and Ophiuchus ones (Hernández et al. 2007; Fedele et al. 2010), viscous evolution could be a candidate mechanism to explain the systematic difference found in the two disk populations. Indeed, in the context of viscously evolving disks, while the inner disk material accretes onto the star, the outer disk radius spreads outwards. For the typical self-similar solutions adopted in this study (Eq. (1)), the disk mass and the exponential cut-off radius evolve such that, at any given time, , where M_disk is the total disk mass. As a result, the total disk mass decreases with time, while the disk size increases, compatible with the difference that we observe between the Lupus and the Taurus-Auriga and Ophiuchus disks in Fig. 11.

We note that while a viscously evolving disk model would explain the direction of the mass-size offset between the Lupus and the Taurus-Auriga and Ophiuchus populations in Fig. 11 (with older disks being larger and fainter), estimating whether the extent of such offset is compatible with the age difference remains problematic as it would require making assumptions on the viscosity and on the evolutionary stage (early or late phase w.r.t the viscous time scale) at least of one of the two populations.

5.2.1. Disk mm fluxes: caveats

Viscous disk evolution may not be the only mechanism that could produce a time evolution in the luminosity-size plot Fig. 11. Other processes such as disk photoevaporation, grain growth, radial drift and planet formation itself might induce time changes in the observed integrated sub-mm flux and dust outer edge. As an example, the generation of pressure maxima in the gas disk (e.g., induced by tidal interaction with forming planets) would substantially change the growth efficiency (and therefore the opacity) and migration rate of disk solids w.r.t. to a simpler scenario of a full disk with continuous mass distribution and a radially decreasing pressure profile. However, a detailed estimate of the signatures of these effects on the luminosity-size plot would require accurate numerical simulations of global disk evolution, which is beyond the scope of this paper. The current observations at 890 μm do not allow us to rule out these alternative processes. However, in the future, the application of a homogeneous analysis to observations of disks in other star forming regions with different mean ages will help in gaining a global view of disk evolution. In addition, the collection of observations at multiple sub-mm/mm wavelengths would reveal the grain growth level in these disks (Tazzari et al. 2016), providing important constraints on the dust emissivity that could be used to break the degeneracies between the different scenarios.

Theoretically the self-similar solutions (and their diffusive behavior) characterize the evolution of the gaseous component, while in this work we have constrained disk masses and sizes from the dust continuum emission. It is thus important to understand what are the potential biases in our estimates when comparing to the theoretical expectations. Regarding the disk masses, recent work by Manara et al. (2016) on the same Lupus data set found a correlation between the disk masses (inferred from the dust continuum emission) and the mass accretion rates onto the star (derived from uv excess measurements), in agreement with the theoretical expectations of viscous theory (Hartmann et al. 1998) for which M_disk/Ṁ_acc ≈ t_age. This broad correlation suggests that, despite the decoupling of the dust and gas in disks, the mm-continuum observations that we have discussed here provide a useful proxy of the total disk mass.

5.2.2. Outer radius: caveats

In this work we derived disk sizes assuming that the 890 μm dust continuum emission is tracing most of the disk material, namely that the dust is co-located to the gas, while in general this might not be the case. At this sub-mm wavelength the observations are mostly sensitive to the thermal emission of the large mm-sized grains and therefore might not be recovering the full spatial extent of the smaller dust particles. Indeed, while the small dust particles are tightly coupled to the gas and therefore are expected to be brought to large radii from the gas viscous spreading (in this respect, they would be good tracers of the gas distribution), as soon as they grow they become less coupled to the gas and start being subject to the inward pointing radial drift (Weidenschilling 1977; Brauer et al. 2008). The cut-off radius of the 890 μm dust continuum distribution is therefore a result of these two effects and it is not trivial to assess which one is dominating as they both depend on the grain size as well as on local gas properties (density, pressure, turbulence) that change with time.

Birnstiel & Andrews (2014) demonstrated that the combination of radial drift and viscous gas drag generates a sharp cut-off at the outer edge of the dust distribution already in the very early phases of disk evolution (before grain growth has had time to take place) and that such a feature is preserved throughout the disk viscous evolution. They also showed that in the late phases of disk evolution the dust outer edge follows the gas edge, which is roughly a factor 1.5× larger. This globally supports our finding of γ ~ 0 in many cases (i.e., disks with flat interior and a sharp cut-off) but is still unclear how the Birnstiel & Andrews (2014) results translate for a population of mm-sized grains. However, a detailed study of the time evolution of the outer edge of different grain size populations in viscously evolving disks is still missing.

5.2.3. Comparison with the minimum mass solar nebula

In Fig. 13 we collect the dust surface density profiles obtained for the disks that we fit. In the left plot we report all the inferred profiles, in the right plot only those of the disks orbiting stars with mass in the range 0.7 < M_⋆/M_⊙< 1. We compare these profiles with the surface density in solids of the minimum mass solar nebula (MMSN, Weidenschilling 1977) which estimates the dust surface density profile of the solar system primordial disk. In particular, we use the MMSN normalization by Hayashi (1981).

The surface density profiles of the Lupus disks appear generally less massive than the MMSN, only a few of them having a comparable or larger mass. We also note that the MMSN profile is smaller in size and has a steeper profile than the Lupus disks: considering that the Lupus disks are still young by planet formation time scales, this might suggest that the solar system disk was probably much more compact, or that the planets (or even the planetesimals) migrated very far inwards. Following the fact that the Lupus disks orbit stars with M_⋆ < M_⊙, it is also possible that the difference with the MMSN reflects an intrinsically lower initial disk mass distribution or a different – mass-dependent – time evolution of the disk surface density. For comparison, the surface density profiles of the disks in Taurus-Auriga and Ophiuchus (Andrews 2015) result more massive than the MMSN, in many cases by a factor between 4 and 10. This might in part be due to the bias towards higher masses of the Taurus-Auriga and Ophiuchus sample, or might be another signature that the Lupus disks are more evolved – and therefore dust depleted – than those in Taurus-Auriga and Ophiuchus.

6. Conclusions

In this paper we have analyzed the 890 μm continuum emission of 22 disks in the Lupus SFR which have been observed with ALMA at ~0.3′′ (~50 au) resolution (Ansdell et al. 2016).

• 1.

  We fit the spatially resolved disk continuum emission with aself-consistent disk model based on the two layer approximationand a realistic dust opacity computation. The fits are performeddirectly in the uv-plane.

• 2.

  For each disk, we derive the dust surface density profile Σ_d(R), the midplane temperature profile T_mid(R), and we constrain the disk inclination i and position angle PA. The disk (gas + dust) surface density profile is inferred with a constant gas-to-dust ratio of 100.

• 3.

  The disk masses are computed by integrating 100 × Σ_d(R) between an inner radius of 0.1 au and an outer radius R_out that contains 95% of the continuum flux. The masses are compatible within a factor of 2 with those derived by Ansdell et al. (2016) assuming a constant T_d = 20 K. Our masses tend to be slightly larger as our disk-averaged temperatures ⟨ T_mid ⟩ are generally smaller than 20 K. We observed no trend between T_mid and the stellar mass and luminosity.

• 4.

  The surface density profiles assumed self-similar solution of viscously evolving disks. The average radial slope for the sample of 22 disks is ⟨ γ ⟩ = 0 ± 0.6, calling for a flat disk interior and sharp outer disk edges. 18 out of 22 disks have cut-off radii R_c< 75 au. We observe no correlation between γ and R_c.

• 5.

  The spatially resolved continuum emission of almost all the disks is reproduced very well (residuals less than 3σ) by a smooth, monotonically decreasing, self-similar Σ_d(R), except for Sz 98 and IM Lup, for which we find a ring-like residual emission with flux ≲15% of the integrated flux.

• 6.

  We find a correlation between the sub-mm integrated flux and the cut-off radius inferred for the disks, which can be interpreted in terms of a disk mass-size correlation. A similar correlation was claimed for disks in the slightly younger regions of Taurus-Auriga and Ophiuchus. We observe that the Lupus disks appear generally larger and fainter (less massive) than those in Taurus-Auriga and Ophiuchus. The spreading of the disk material (gas and subsequently dust) induced by the viscous evolution can be a possible explanation of such difference but we could not rule out other processes that might be at play.

• 7.

  The Σ_d profiles inferred for the Lupus disks are generally less massive than the solar system protoplanetary disk. Also, the Lupus disks around stars with mass ~ 1 M_⊙ have a shallower Σ_d and a larger outer radius than the solar system protoplanetary disk, possibly suggesting that the planets (or already the planetesimals) in the solar nebula migrated very far inward during or after their formation.

---

Acknowledgments

M.T. and L.T. acknowledge support by the DFG cluster of excellence Origin and Structure of the Universe (www.universe-cluster.de). This work has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. M.T. and L.T. thank Eva Wirström, the Gothenburg Centre for Advanced Studies in Science and Technology and the participants of the workshop Origins of Habitable Planets hosted at Chalmers (Gothenburg) in June 2016 where considerable part of this work was fruitfully discussed. M.T. is grateful to Cathie Clarke for insightful discussions. The fits have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP) as part of the approved project “Dust evolution in protoplanetary disks”. M.T. acknowledges usage of the Max Planck Gesellschaft Hydra computing cluster for code testing, for which is grateful to Paola Caselli. Figures have been generated using the Python-based matplotlib package (Hunter 2007). Staircase plots of PDFs have been generated with a user-modified version of the Python-based triangle package (Foreman-Mackey et al. 2014). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This work was partly supported by the Italian Ministero dell’ Istruzione, Università e Ricerca through the grant Progetti Premiali 2012 – iALMA (CUP C52I13000140001). A.N. acknowledges funding from Science Foundation Ireland (Grant 13/ERC/I2907). C.F.M. acknowledges an ESA Research Fellowship. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00220.S, ADS/JAO.ALMA#2013.1.00226.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ

References

Appendix A: Fits of the individual sources

Here we report the results of the fits for the individual sources, following the order in Table 3. For each disk we report:

• a triangle plot showing the MCMC results, providing estimates ofthe parameters as median of the marginalized posterior;

• a deprojected uv-plot showing a comparison between observed and best-fit model visibilities;

• radial profiles of the total (dust + gas) surface density 100 × Σ_d(R), midplane temperature T_mid(R), optical depth τ_{890 μm};

• synthesized maps of the observations, of the best-fit model and of the residuals.

All the models fit well the observations with negligible residuals, except for Sz 98 and IM Lup where a smoothly decreasing surface density profile is not sufficient to explain the brightness distribution. In these cases ring-like residuals are left, respectively 14% and 4% of the total intensity, and are likely to be due to the presence of inhomogeneities like rings or spirals (ALMA Partnership et al. 2015; Guidi et al. 2016; Isella et al. 2016).

Appendix B: Results of Bayesian linear regressions

All Tables

All Figures