EDP Sciences
Free Access
Issue
A&A
Volume 606, October 2017
Article Number A88
Number of page(s) 39
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201730890
Published online 17 October 2017

© 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) = Σ0(R/R0)γ, which can be either truncated at a radius Rout or – following the arguments of accretion theory (Lynden-Bell & Pringle 1974; Hartmann et al. 1998) – exponentially tapered with an e-folding radius Rc. 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).

Table 1

Source properties.

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 Rc 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 Fcont< 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.

thumbnail Fig. 1

Top: distribution of stellar masses for the full Lupus sample from Ansdell et al. (2016) (white bars) and for the sub-sample analyzed in this work (blue bars). Bottom: observed integrated continuum flux at 890 μm (Ansdell et al. 2016) as a function of stellar mass. Blue dots represent sources whose analysis is presented in this work. Red dots represent sources that appear unresolved and have been discarded in the rest of the analysis due low signal-to-noise ratio. In this and subsequent plots, the circled dots highlight sources in the 0.7–1 M mass bin. We also report transition disks (“×” symbols), sources with irregular shape (hexagons), binaries (diamonds), and the sources with integrated flux Fcont< 4 mJy that we excluded from our analysis (empty circles).

Open with DEXTER

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 Σ0 is a normalization factor, R0 is a characteristic radius that we keep fixed to 10 au, γ is the power-law slope and Rc 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: Σ0, γ and Rc.

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) ∝ aq for aminaamax and n(a) = 0 otherwise, where a is the grain radius. In the surface layer we use amin = 10 nm and amax = 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 amin = 10 nm and amax = 1.023 cm with q = 3.0, which corresponds to a population of larger grains. The particular value of amax 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 cm2 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.

thumbnail Fig. 2

Left: spectral indices for the sources in our sample measured between J-band (2MASS, 1.235 μm) and, respectively, 70 μm, 100 μm, 160 μm (Herschel/PACS) using the definition in Eq. (3). The vertical lines show the spectral indices obtained with our disk model for different values of the flaring reduction parameter f, from a fully-flared model (f = 1, dotted line) to less flared models (f = 0.5,0.3,01, respectively dash-dotted, dashed and solid line). Right: scale-height of the disk surface layer computed by the model.

Open with DEXTER

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 Hsurf/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 Rc 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: Σ0, γ and Rc 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 study1.

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 χ2 as: (4)where Vobs and Vmod are the observed and the synthetic visibilities, respectively, wi is the weight associated with the observed visibilities at the (ui,vi) 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 iwi = 1. The posterior PDF is computed as exp(−χ2/ 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.

Table 2

Parameter space explored by the Markov chains.

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 × 105 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 marginalizing2 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, Σ0 = 16.2 ± 0.25 gcm-2, Rc = 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.

thumbnail Fig. 3

Staircase plot of the chain resulting from the MCMC fit of Sz 71. The histograms on the main diagonal are the marginalized distributions of each parameter: from left to right, γ, Σ0, Rc, i and PA, Δα, Δδ. The vertical dashed lines indicate the 16%, 50% and 84% percentiles. The off-diagonal 2D plots show the correlation between each couple of parameter, with contour lines showing 0.5σ increments. The solid red lines highlight the coordinates of the best-fit model.

Open with DEXTER

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 Rc ~ 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.

thumbnail Fig. 4

Physical structure of the best-fit model describing Sz 71. Left: gas surface density profile computed as 100 × Σd(R), where 100 is the gas-to-dust ratio. The solid red line is the cumulative mass (right y axis) and the black filled circle highlights the cut-off radius Rc ≃ 85 au. In all three plots, the dotted red line indicates the radius containing 95% of mass and the vertical dashed black line half the synthesized beam size (25 ) au. Middle: dust temperature profile Tmid(R) derived self-consistently at each radius from the disk model. The dot-dashed horizontal line represents the mass-averaged dust temperature Tmid, which for Sz 71 is 14 K. Right: optical depth of the disk midplane at the observing wavelength: the sub-mm continuum emission is expected to be optically thick only within R < 2 au.

Open with DEXTER

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.

thumbnail Fig. 5

Comparison between model and observed visibilities of Sz 71 as a function of deprojected baseline length (uv-distance). The panel above (below) shows the real (imaginary) part of the visibilities. Black dots show the observed visibilities (binned every 30 kλ), the blue density indicators represent the PDF of the model visibilities in each uv-bin. The black solid line corresponds to the median of the model visibilities PDF, the black dashed lines the 16th and 84th percentiles, the red dotted lines the 2.7th and 97.7th percentiles.

Open with DEXTER

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).

thumbnail Fig. 6

Comparison between the observed (left) and the best-fit model (middle) images for Sz 71. Residuals are shown in the right panel. The best-fit model has the following parameters γ = 0.25 ± 0.01, Σ0 = (16.2 ± 0.2)g/cm2,Rc = (85 ± 1)au. The three images have been produced by applying the CLEAN deconvolution algorithm with natural weighting to the observed, best-fit model and residual visibilities, respectively. Contour levels refer to –3 (dashed), 3, 6, 12, 24, 48, 50, 100, 150, etc. multiple of the rms, which here is σ = 0.3 mJy/beam. The FWHM of the synthesized beam is shown as a gray ellipse in the left panel.

Open with DEXTER

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.

Table 3

Fit results: free parameters and derived quantities.

thumbnail Fig. 7

Left: ratio between the dust masses derived in this work and those derived by Ansdell et al. (2016). Although derived with very different methods (for a detailed description, see the text), the masses are within a factor of 2 (gray area) for the vast majority of disks. Right: comparison between the disk inclinations i (along the line of sight) obtained in this work and those derived by Ansdell et al. (2016). In both plots, the disks are sorted from left to right by increasing stellar mass and the circled dots highlight sources with 0.7 ≤ M/M ≤ 1.

Open with DEXTER

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 γ, Σ0, Rc, 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 Rout, defined as the radius containing 95% of the model flux, and the total dust mass Mdust, computed by integrating the dust surface density: (5)where Rin = 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 cm2 g-1 (per gram of dust) is the dust opacity, Bν(T) is the black-body brightness at the temperature T, and Td = 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 Td = 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 Tmid(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 (Tmid) 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 Tmid as a function of stellar mass and luminosity, we conclude that the assumption of a constant Td (Ansdell et al. 2016) has not introduced a bias in the determination of Md. Conversely the assumptions of Andrews et al. (2013) and van der Plas et al. (2016) on the dependence of the average Td on the stellar parameters may introduce a spurious dependence on the derived Md vs. M relationship.

thumbnail Fig. 8

Mass-averaged midplane temperature as a function of stellar mass (left) and luminosity (right) for the disks for which we derive a reliable disk structure. Error bars on the y-axis reflect the distribution of models obtained from the fits. In more than half disks we find a disk temperature smaller than Td = 20 K used by Ansdell et al. (2016).

Open with DEXTER

thumbnail Fig. 9

Distribution of surface density slopes γ (left), exponential cut-off radii Rc(center) and the correlation between them (right) derived for the disks in the sample. The vertical dashed line in the top and middle panels gives a visual representation of the spatial resolution of the observations, estimated as half of the synthesized beam size (in our case ~ 25 au at a distance of 150 pc). There is no apparent correlation between γ and Rc.

Open with DEXTER

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 Rc (middle) and their correlation (right) obtained from our fits. The distribution of cut-off radii Rc has 6 disks with Rc< 25 au, 14 disks with 25 au ≤ Rc ≤ 100 au and two much larger disks with Rc ~ 200 au (Sz 98) and Rc ~ 430 au (IM Lup). The distribution of disk sizes (Rout, plot not shown) is similar to that of Rc: 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 Rc, in line with the findings of Andrews et al. (2009, 2010) in Ophiuchus (14 au ≤ Rc ≤ 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 Rc 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.

Table 4

Summary of fit results for disks with cavities.

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 Rhole ≈ 46 au, where we defined Rhole 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 Rhole ≈ 41 au and Rhole ≈ 34 au, for J16102955-3922144 the fit is more uncertain (uncertainty on γ is large), with a marginal evidence of a cavity at Rhole ≈ 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 Rhole = 30 au and Rhole = 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 Rhole 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.

thumbnail Fig. 10

Surface density profiles of the five disks having γ< 0 at a significance level larger than 3σ. The vertical dashed line gives a visual representation of half the beam size (25 au at 150 pc).

Open with DEXTER

thumbnail Fig. 11

890 μm integrated flux (normalized to a common distance of 150 pc) as a function of exponential cut-off radius, constrained by fitting spatially resolved observations with angular resolution better than 0.75′′. For some disks, the 890 μm flux was extrapolated from the measured 1.3 mm flux assuming an average spectral index α = 3.0. Yellow elements are Taurus-Auriga and Ophiuchus disks from the Andrews (2015) compilation; blue elements are Lupus disks analyzed in this paper. Gray lines represent the correlations found by the Bayesian linear regression for the Taurus-Auriga and Ophiuchus disks (dash-dotted), for the Lupus disks (dashed), and for both the samples taken together (solid).

Open with DEXTER

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 Rc. 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)3 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 test4 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%).

thumbnail Fig. 12

Top panels: distribution of the parameter δfit (see text) for the Lupus (blue) and Taurus-Auriga and Ophiuchus (yellow) samples. In the left panel we show the full samples, while in the right one we restrict to the range 0.7 ≤ M/M ≤ 1. Bottom panels: empirical cumulative distribution functions for the values of δfit for the same samples as in the upper panels.

Open with DEXTER

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 Mdisk 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 Mdisk/acctage. 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.

thumbnail Fig. 13

Left: dust surface density profiles inferred for all the disks in our Lupus sample (gray curves) compared to the Hayashi (1981) dust surface density model for the solar system (black thick line). The inferred surface density profiles are dashed where they are not spatially resolved, that is for R< 25 au. Right: dust surface density profiles for the sources in the 0.7–1M stellar mass bin.

Open with DEXTER

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 Tmid(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 Rout 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 Td = 20 K. Our masses tend to be slightly larger as our disk-averaged temperatures Tmid are generally smaller than 20 K. We observed no trend between Tmid 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 Rc< 75 au. We observe no correlation between γ and Rc.

  • 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.


1

In principle it is possible to relate the offsets inferred from the fit to the proper motion of each single object as shown in Tazzari et al. (2016), but this is beyond the scope of this paper.

2

Integrating over all but the one parameter of interest.

3

We use the linmix Python package which implements Kelly (2007) and is available here: https://github.com/jmeyers314/linmix

4

We used the implementation of the test as provided in scipy.stats.anderson_ksamp, which is based on Scholz & Stephens (1987).

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 Tmid(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).

thumbnail Fig. A.1

Fit results for Sz 65. Top row: on the left, the staircase plot showing the MCMC chains as in Fig. 3; on the right, the comparison of the observations and model deprojected visibilities as a function of uv-distance as in Fig. 5. Middle row: plots showing the physical structure of the disk as in Fig. 4. Bottom row: synthesized images of observations, model and residual visibilities as in Fig. 6. In the images σ = 0.27 mJy/beam.

Open with DEXTER

thumbnail Fig. A.2

Fit results for J15450887-3417333, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER

thumbnail Fig. A.3

Fit results for Sz 68, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER

thumbnail Fig. A.4

Fit results for Sz 69, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER

thumbnail Fig. A.5

Fit results for Sz 71, presented as in Fig. A.1. In the images σ = 0.23 mJy/beam.

Open with DEXTER

thumbnail Fig. A.6

Fit results for Sz 73, presented as in Fig. A.1. In the images σ = 0.33 mJy/beam.

Open with DEXTER

thumbnail Fig. A.7

Fit results for IM Lup, presented as in Fig. A.1. In the images σ = 0.6 mJy/beam.

Open with DEXTER

thumbnail Fig. A.8

Fit results for Sz 83, presented as in Fig. A.1. In the images σ = 0.4 mJy/beam.

Open with DEXTER

thumbnail Fig. A.9

Fit results for Sz 84, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER

thumbnail Fig. A.10

Fit results for Sz 129, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER

thumbnail Fig. A.11

Fit results for J16000236-4222145, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER

thumbnail Fig. A.12

Fit results for MY Lup, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER

thumbnail Fig. A.13

Fit results for Sz 133, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER

thumbnail Fig. A.14

Fit results for Sz 90, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER

thumbnail Fig. A.15

Fit results for Sz 98, presented as in Fig. A.1. In the images σ = 0.33 mJy/beam.

Open with DEXTER

thumbnail Fig. A.16

Fit results for Sz 100, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER

thumbnail Fig. A.17

Fit results for Sz 108B, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER

thumbnail Fig. A.18

Fit results for J16085324-3914401, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER

thumbnail Fig. A.19

Fit results for Sz 113, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER

thumbnail Fig. A.20

Fit results for Sz 114, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER

thumbnail Fig. A.21

Fit results for J16102955-3922144, presented as in Fig. A.1. In the maps, σ = 0.25 mJy/beam.

Open with DEXTER

thumbnail Fig. A.22

Fit results for J16124373-3815031, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER

Appendix B: Results of Bayesian linear regressions

thumbnail Fig. B.1

Staircase plots showing the posterior PDF resulting from the application of the Bayesian regression method by Kelly (2007) to infer a law log (F890 μm) = βlog (Rc) + α for the Lupus sample (top left), for the Taurus-Auriga and Ophiuchus sample (top right) and for both the samples together (bottom). Each plot shows the PDF of the intercept (α), the slope (β), the intrinsic scatter of data points (σ) and the correlation degree (corr). The red lines highlight the posterior peak.

Open with DEXTER

All Tables

Table 1

Source properties.

Table 2

Parameter space explored by the Markov chains.

Table 3

Fit results: free parameters and derived quantities.

Table 4

Summary of fit results for disks with cavities.

All Figures

thumbnail Fig. 1

Top: distribution of stellar masses for the full Lupus sample from Ansdell et al. (2016) (white bars) and for the sub-sample analyzed in this work (blue bars). Bottom: observed integrated continuum flux at 890 μm (Ansdell et al. 2016) as a function of stellar mass. Blue dots represent sources whose analysis is presented in this work. Red dots represent sources that appear unresolved and have been discarded in the rest of the analysis due low signal-to-noise ratio. In this and subsequent plots, the circled dots highlight sources in the 0.7–1 M mass bin. We also report transition disks (“×” symbols), sources with irregular shape (hexagons), binaries (diamonds), and the sources with integrated flux Fcont< 4 mJy that we excluded from our analysis (empty circles).

Open with DEXTER
In the text
thumbnail Fig. 2

Left: spectral indices for the sources in our sample measured between J-band (2MASS, 1.235 μm) and, respectively, 70 μm, 100 μm, 160 μm (Herschel/PACS) using the definition in Eq. (3). The vertical lines show the spectral indices obtained with our disk model for different values of the flaring reduction parameter f, from a fully-flared model (f = 1, dotted line) to less flared models (f = 0.5,0.3,01, respectively dash-dotted, dashed and solid line). Right: scale-height of the disk surface layer computed by the model.

Open with DEXTER
In the text
thumbnail Fig. 3

Staircase plot of the chain resulting from the MCMC fit of Sz 71. The histograms on the main diagonal are the marginalized distributions of each parameter: from left to right, γ, Σ0, Rc, i and PA, Δα, Δδ. The vertical dashed lines indicate the 16%, 50% and 84% percentiles. The off-diagonal 2D plots show the correlation between each couple of parameter, with contour lines showing 0.5σ increments. The solid red lines highlight the coordinates of the best-fit model.

Open with DEXTER
In the text
thumbnail Fig. 4

Physical structure of the best-fit model describing Sz 71. Left: gas surface density profile computed as 100 × Σd(R), where 100 is the gas-to-dust ratio. The solid red line is the cumulative mass (right y axis) and the black filled circle highlights the cut-off radius Rc ≃ 85 au. In all three plots, the dotted red line indicates the radius containing 95% of mass and the vertical dashed black line half the synthesized beam size (25 ) au. Middle: dust temperature profile Tmid(R) derived self-consistently at each radius from the disk model. The dot-dashed horizontal line represents the mass-averaged dust temperature Tmid, which for Sz 71 is 14 K. Right: optical depth of the disk midplane at the observing wavelength: the sub-mm continuum emission is expected to be optically thick only within R < 2 au.

Open with DEXTER
In the text
thumbnail Fig. 5

Comparison between model and observed visibilities of Sz 71 as a function of deprojected baseline length (uv-distance). The panel above (below) shows the real (imaginary) part of the visibilities. Black dots show the observed visibilities (binned every 30 kλ), the blue density indicators represent the PDF of the model visibilities in each uv-bin. The black solid line corresponds to the median of the model visibilities PDF, the black dashed lines the 16th and 84th percentiles, the red dotted lines the 2.7th and 97.7th percentiles.

Open with DEXTER
In the text
thumbnail Fig. 6

Comparison between the observed (left) and the best-fit model (middle) images for Sz 71. Residuals are shown in the right panel. The best-fit model has the following parameters γ = 0.25 ± 0.01, Σ0 = (16.2 ± 0.2)g/cm2,Rc = (85 ± 1)au. The three images have been produced by applying the CLEAN deconvolution algorithm with natural weighting to the observed, best-fit model and residual visibilities, respectively. Contour levels refer to –3 (dashed), 3, 6, 12, 24, 48, 50, 100, 150, etc. multiple of the rms, which here is σ = 0.3 mJy/beam. The FWHM of the synthesized beam is shown as a gray ellipse in the left panel.

Open with DEXTER
In the text
thumbnail Fig. 7

Left: ratio between the dust masses derived in this work and those derived by Ansdell et al. (2016). Although derived with very different methods (for a detailed description, see the text), the masses are within a factor of 2 (gray area) for the vast majority of disks. Right: comparison between the disk inclinations i (along the line of sight) obtained in this work and those derived by Ansdell et al. (2016). In both plots, the disks are sorted from left to right by increasing stellar mass and the circled dots highlight sources with 0.7 ≤ M/M ≤ 1.

Open with DEXTER
In the text
thumbnail Fig. 8

Mass-averaged midplane temperature as a function of stellar mass (left) and luminosity (right) for the disks for which we derive a reliable disk structure. Error bars on the y-axis reflect the distribution of models obtained from the fits. In more than half disks we find a disk temperature smaller than Td = 20 K used by Ansdell et al. (2016).

Open with DEXTER
In the text
thumbnail Fig. 9

Distribution of surface density slopes γ (left), exponential cut-off radii Rc(center) and the correlation between them (right) derived for the disks in the sample. The vertical dashed line in the top and middle panels gives a visual representation of the spatial resolution of the observations, estimated as half of the synthesized beam size (in our case ~ 25 au at a distance of 150 pc). There is no apparent correlation between γ and Rc.

Open with DEXTER
In the text
thumbnail Fig. 10

Surface density profiles of the five disks having γ< 0 at a significance level larger than 3σ. The vertical dashed line gives a visual representation of half the beam size (25 au at 150 pc).

Open with DEXTER
In the text
thumbnail Fig. 11

890 μm integrated flux (normalized to a common distance of 150 pc) as a function of exponential cut-off radius, constrained by fitting spatially resolved observations with angular resolution better than 0.75′′. For some disks, the 890 μm flux was extrapolated from the measured 1.3 mm flux assuming an average spectral index α = 3.0. Yellow elements are Taurus-Auriga and Ophiuchus disks from the Andrews (2015) compilation; blue elements are Lupus disks analyzed in this paper. Gray lines represent the correlations found by the Bayesian linear regression for the Taurus-Auriga and Ophiuchus disks (dash-dotted), for the Lupus disks (dashed), and for both the samples taken together (solid).

Open with DEXTER
In the text
thumbnail Fig. 12

Top panels: distribution of the parameter δfit (see text) for the Lupus (blue) and Taurus-Auriga and Ophiuchus (yellow) samples. In the left panel we show the full samples, while in the right one we restrict to the range 0.7 ≤ M/M ≤ 1. Bottom panels: empirical cumulative distribution functions for the values of δfit for the same samples as in the upper panels.

Open with DEXTER
In the text
thumbnail Fig. 13

Left: dust surface density profiles inferred for all the disks in our Lupus sample (gray curves) compared to the Hayashi (1981) dust surface density model for the solar system (black thick line). The inferred surface density profiles are dashed where they are not spatially resolved, that is for R< 25 au. Right: dust surface density profiles for the sources in the 0.7–1M stellar mass bin.

Open with DEXTER
In the text
thumbnail Fig. A.1

Fit results for Sz 65. Top row: on the left, the staircase plot showing the MCMC chains as in Fig. 3; on the right, the comparison of the observations and model deprojected visibilities as a function of uv-distance as in Fig. 5. Middle row: plots showing the physical structure of the disk as in Fig. 4. Bottom row: synthesized images of observations, model and residual visibilities as in Fig. 6. In the images σ = 0.27 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.2

Fit results for J15450887-3417333, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.3

Fit results for Sz 68, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.4

Fit results for Sz 69, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.5

Fit results for Sz 71, presented as in Fig. A.1. In the images σ = 0.23 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.6

Fit results for Sz 73, presented as in Fig. A.1. In the images σ = 0.33 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.7

Fit results for IM Lup, presented as in Fig. A.1. In the images σ = 0.6 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.8

Fit results for Sz 83, presented as in Fig. A.1. In the images σ = 0.4 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.9

Fit results for Sz 84, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.10

Fit results for Sz 129, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.11

Fit results for J16000236-4222145, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.12

Fit results for MY Lup, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.13

Fit results for Sz 133, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.14

Fit results for Sz 90, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.15

Fit results for Sz 98, presented as in Fig. A.1. In the images σ = 0.33 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.16

Fit results for Sz 100, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.17

Fit results for Sz 108B, presented as in Fig. A.1. In the images σ = 0.25 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.18

Fit results for J16085324-3914401, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.19

Fit results for Sz 113, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.20

Fit results for Sz 114, presented as in Fig. A.1. In the images σ = 0.2 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.21

Fit results for J16102955-3922144, presented as in Fig. A.1. In the maps, σ = 0.25 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. A.22

Fit results for J16124373-3815031, presented as in Fig. A.1. In the images σ = 0.3 mJy/beam.

Open with DEXTER
In the text
thumbnail Fig. B.1

Staircase plots showing the posterior PDF resulting from the application of the Bayesian regression method by Kelly (2007) to infer a law log (F890 μm) = βlog (Rc) + α for the Lupus sample (top left), for the Taurus-Auriga and Ophiuchus sample (top right) and for both the samples together (bottom). Each plot shows the PDF of the intercept (α), the slope (β), the intrinsic scatter of data points (σ) and the correlation degree (corr). The red lines highlight the posterior peak.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.