Millimeter (mm) and submillimeter (submm) wavelength observations are particularly useful to study dust properties in protoplanetary disks because the dust thermal emission of the outer disk, where the bulk of the dust mass resides, can be probed at these wavelengths (Testi et al. 2014;Andrews 2015). Demographic studies based on mm and submm wavelength surveys of the Class II population from nearby star-forming regions (Ansdell et al. 2016;Barenfeld et al. 2016;Pascucci et al. 2016;Cox et al. 2017;Cieza et al. 2018;Cazzoletti et al. 2019) have found positive correlations between various disk properties: disk mass with stellar mass (M disk -M , Andrews et al. 2013;Ansdell et al. 2016;Pascucci et al. 2016), disk size with luminosity (Andrews et al. 2010(Andrews et al. , 2018;;Tazzari et al. 2017;Tripathi et al. 2017), and mass accretion rate onto the central star with disk mass ( Ṁacc -M disk , Manara et al. 2016b;Mulders et al. 2017;Rosotti et al. 2017).
These relations are poorly constrained in the brown dwarf (BD) and very low-mass (VLM) star regimes because these surveys focused primarily on disks around more massive stars. Therefore, observations at mm and submm wavelengths targeting BDs and VLM stars are necessary in order to extend these demographic studies and to investigate their formation mechanisms and ongoing physical processes in their disks.
General interest in BDs and VLM stars has increased substantially thanks to the recent exoplanet discoveries around VLM objects. The most thrilling cases are Trappist-1 (Gillon et al. 2017), a ∼0.085 M VLM star that hosts seven rocky planets in Article published by EDP Sciences A114, page 1 of 23 A&A 633, A114 (2020) a packed orbital configuration, and Proxima B (Anglada-Escudé et al. 2016), an Earth-like planet orbiting our closest neighboring star (M = 0.12 M ), at a distance of only ∼1.3 pc from Earth. These and other discoveries (e.g., 2M1207b and 2M J044144b, Chauvin et al. 2004;Todorov et al. 2010) suggest that planets orbiting BDs and VLMs may be a common outcome of their formation.
The study of the early stages of BDs and VLM stars is crucial to understanding the viability of planet formation around these objects and to determine the properties of the potential planetary systems that may form. In Klein et al. (2003), millimeter emission of dust from disks around BDs was detected for the first time. Like stars, BDs are often found surrounded by a protoplanetary disk in their early stages (Comeron et al. 1998;Natta & Testi 2001;Scholz 2008), where planet formation is expected to take place. The disk fraction for stellar and BD populations is found to be similar (Luhman 2012). Disk accretion (Jayawardhana et al. 2003;Scholz & Eislöffel 2004;Muzerolle et al. 2003Muzerolle et al. , 2005) ) and outflows (Natta et al. 2004;Whelan et al. 2005) also occur in the early stages of BDs, analogous to those around more massive stars.
In this study we conducted a systematic survey of BD disks in the Lupus star-forming region, observing the full known sample of BD disks from a single region with the Atacama Large Millimeter/submillimeter Array (ALMA) in the same band for the first time. Previous ALMA observations of BD disks studied incomplete samples of the known BD population of other regions (Testi et al. 2016;van der Plas et al. 2016;Ward-Duong et al. 2018).
Dust disk masses, dust emission distribution profiles, and dust disk characteristic sizes are determined from these observations. The last two properties are inferred from interferometric modeling of the dust disk emission. The characteristic size of the dusty disks is crucial to constraining the ongoing disk evolution processes (e.g., radial drift, grain growth). However, its determination is not straightforward. Firstly, the disk emission needs to be sufficiently resolved. For disks around BDs and VLM stars, this is only possible using state-of-art facilities, like ALMA, that provide the high resolution and sensitivity required at these wavelengths. In addition, a general size definition is needed for a reliable comparison between observations and theoretical models. In this work we use the radius enclosing 68% of the object's emission distribution; this definition is representative of the physical size of the object (Tripathi et al. 2017), and is independent of the model used to fit the observations. Another important disk property that can be derived from submm observations is the disk mass. For the formation of rocky planets, the dust mass in disks should be larger than the mass of the resulting planets. However, comparing the mass derived from disk emission with the results from exoplanetary surveys, there is an apparent lack of material to produce the known planetary systems (Greaves & Rice 2010;Williams 2012;Najita & Kenyon 2014;Mulders et al. 2015;Pascucci et al. 2016;Testi et al. 2016;Manara et al. 2018).
The inferred disk properties of the young BD population are compared to the properties of disks around stars in the same region, with the aim of testing whether the known relations for stars hold for disks around BDs. In addition to the BD observations, a further seven T Tauri star (TTS) disks are characterized and modeled here for the first time.
The study is organized as follows: the target selection is described in Sect. 2. A summary of the observations and the data processing can be found in Sect. 3. Section 4 provides a description of the modeling employed for the determination of disk properties, together with the modeling results. The demographic comparison of the inferred properties between BD and stellar disks, and the planet formation implications from the measured dust masses of the BD disks are discussed in Sect. 5, and the main conclusions of this work are presented in Sect. 6.
The list of selected targets of the Lupus BD disks survey (Cycle 5; PI: L. Testi, Project ID: 2017.1.01243.S) encompasses all the known BDs in the Lupus region I-IV that were not observed previously with ALMA Band 7. Our population of BD disks in the Lupus star-forming region consists of all the known objects from the region census (Merín et al. 2008;Mužić et al. 2014Mužić et al. , 2015) ) that show excess emission in at least two mid-infrared bands (Spitzer IRAC/MIPS). These objects have been spectroscopically classified as spectral type (SpT) M6 or later (excluding L type or later candidates), and with estimated masses of ≤0.09 M . L and later type candidates have been excluded. Eleven sources in Lupus satisfied these selection criteria; seven were the targets for the new observations and the remaining four had already been observed in the Lupus disks survey (Cycle 2; PI: J. Williams, Project ID: 2013.1.00220.S). Based on radial velocity analysis and X-shooter spectra, two sources in the sample, IRAS 15567-4141 and SSTc2d J160034.4-422540, were recently excluded from being Lupus members, and are likely background giants (Frasca et al. 2017;Alcalá et al. 2017), in agreement with the poorly constrained parallaxes from Gaia DR2 (Gaia Collaboration 2018). Therefore, these will not be discussed further in this paper; only five new targets are discussed. All the studied BDs are isolated systems, except SONYC-Lup3-7, which might form a very wide (∼7 ) binary system with SONYC-Lup3-6, but this last object has no confirmed membership to the region (Mužić et al. 2014).
In Table 1, we list all the sources analyzed in this work: the known BD population (five objects from the new observations and four from previous observations), together with the seven disks around stars observed in the Lupus completion survey (Cycle 5; PI: S.E. van Terwisga, Project ID: 2016.1.01239.S). The last two objects in the table are those that were observed but later excluded from the Lupus census (Frasca et al. 2017;Alcalá et al. 2017). The names, sky positions, and main stellar properties of the central stars are included in the Table 1. The stellar properties shown in the table (SpT, effective temperature T eff , extinction in V-band A V , stellar luminosity L and M ) were reported by Alcalá et al. (2014Alcalá et al. ( , 2017) ) and Mužić et al. (2014). The methodology for the stellar luminosity derivation between these studies differs, but nevertheless the agreement between the two methods is very good, as shown in Manara et al. (2016a). Values for L have been adjusted accounting for updated distances by Gaia DR2 (distance estimated as the inverse of the parallax, Gaia Collaboration 2018). The stellar mass is derived from the pre-main sequence (MS) evolutionary models of Baraffe et al. (2015), estimated from the position in the Hertzsprung-Russell (HR) diagram. For objects with estimated mass >1.4 M and objects laying above the 1 Myr isochrone, the tracks from Siess et al. (2000) are used instead. Stellar mass uncertainties are computed with a Monte Carlo approach (as described in Alcalá et al. 2017), which takes into account the associated uncertainties of the stellar properties L and T eff used to infer the mass.
The BD disk population is compared to the sample of young stellar objects (YSOs) in Lupus that have a protoplanetary disk and estimated stellar mass >0.09 M . Thanks to the inclusion
BDs from this survey J154518. 5-342125 15:45:18.53 -34:21:24 Notes.
(1) Stellar properties from Alcalá et al. (2014Alcalá et al. ( , 2017) ) and adjusted to the new Gaia DR2 parallaxes.
(2) Stellar properties from Mužić et al. (2014), and adjusted to the Gaia DR2 parallaxes.
(3) Gaia parallax unknown, mean distance of Lupus region considered. (4) Two sets of stellar properties (Alcalá et al. 2017;Mužić et al. 2014). ( 5) Subluminous object (see Alcalá et al. 2014Alcalá et al. , 2017)). ( 6) No estimation of the stellar mass since it falls below the zero-age main sequence in the HR diagram (see Alcalá et al. 2017). ( 7) Categorized as background sources (Frasca et al. 2017;Alcalá et al. 2017).
of the seven stellar disks from the Lupus completion survey, we have the largest sample of stellar disks in the Lupus clouds (I-IV) observed with ALMA in Band 7. The stellar disk population is assembled from different census of the region from 2MASS, Spitzer, and Herschel surveys (Hughes et al. 1994;Comerón 2008;Merín et al. 2008;Mortier et al. 2011;Dunham et al. 2015;Bustamante et al. 2015). Of these disks 82 were observed with ALMA in Band 7 in the original Lupus disks survey (ALMA Cycle 2; PI: J. Therefore, the BD and stellar disk samples of our demographic study consist of 9 and 91 sources respectively. The HR-diagram for the studied disk population is shown in Fig. 1 using the stellar properties from previous studies as described before. The BD disks are marked in red, and the stellar disk population in blue. Subluminous sources, those with luminosities lower than those expected for YSOs with an age of ∼3 Myr (likely due to gray obscuration Alcalá et al. 2014Alcalá et al. , 2017)), are illustrated with square symbols. From the re-adjusted L using the more accurate Gaia DR2 parallaxes, two objects (J16081497-3857145 and J16085373-3914367) are now added to the list of subluminous objects of the region. Baraffe et al. (2015) are overlaid in the figure. Objects with luminosities that would correspond to older ages than expected are considered subluminous and marked as squares. A number of points lay on top of each other (e.g., 4 BDs near the 1 Myr and the 0.1 M lines.) Fig. 2. Dust continuum images at 890 µm of the Lupus BDs disks survey from ALMA Band 7 observations. The beam size FWHM is 0.27 × 0.24 for the J154518.5-342125 map (robust parameter of -1), and 0.36 × 0.33 for the rest of the maps (robustness = +0.5). The average beam position angle is PA = 28 • . The contours are drawn at increasing (or decreasing) 3σ intervals as solid (dashed) lines.
ALMA observed our targets on 2018 April 1 and April 2 with 44 and 42 antennas respectively, each of 12 m in diameter. The baselines ranged between 15.1 and 704.1 m for the array configuration of the first day, and between 15.1 and 629.2 m for the configuration of the second day. Four spectral windows were set for the continuum observation of the targets, centered at 334. 432, 336.321, 345.889 and 347.821 GHz and bandwidths of 2, 1.875, 1.875 and 0.938 GHz respectively (total receivers bandwidth of ∼6.688 GHz). The calibrators for the observations were J1517-2422 for flux and passband, and J1610-3958 for the complex gain calibration, the same in both executions. The flux density scale accuracy is expected to be of 10% for observations of the Lupus BD disks survey. Twelve scans of 60-62 seconds duration each were performed for every target, for a total integration time of more than 12 min per source.
The CASA 5.3.0 software has been used for the interferometric visibilities calibration and imaging. The continuum maps are produced using the channels free from spectral line emission, with Briggs weighting of the visibilities (-1.0 and +0.5 robustness for resolved and unresolved objects respectively). None of the sources are bright enough to perform self-calibration. The full width at half maximum (FWHM) of the synthesized beam is 0.27 × 0.24 for robust parameter = -1, and 0.36 × 0.33 for robustness = +0.5, with average position angle (PA) of 28 • . The continuum maps of the five BD disks observed are shown in Fig. 2. The sensitivity for the BD disks survey is improved by a factor of about three with respect to the previous Lupus disk surveys, allowing us to detect fainter emission. Three BD disks are detected, J154518.5-342125, SONYC-Lup3-7, and Lup706, with respective signal-to-noise ratio (S/N) of 42, 8, and 7. Emission is not detected from the two other disks (AKC2006-18 and SONYC-Lup3-10) or from the two background objects (IRAS 15567-4141 and J160034.4-422540).
The main results of the observations are reported in Table 2. This table includes the total disk flux, peak intensity, and the rms of the image. These values from the observations are obtained using identical methodology to the results presented in Ansdell et al. (2016) for the original Lupus disks survey. The continuum flux is inferred from the uvmodelfit task in CASA: emission is fitted with an elliptical Gaussian in cases where the resulting FWHM along the major axis from the fit is at least five times its uncertainty, otherwise the emission is fitted as a point source. For sources with resolved structure, the flux is obtained from a curve of growth method with increasing circular apertures centered at the peak emission of the object. The rms is computed from an annulus of 4-9 radius centered on the detected emission, or on the expected source position if no disk emission is detected. The flux upper limits of the nondetected BD disks are displayed in the F cont column of the table. In this work, upper limits of nondetected BD and stellar disks are computed as three times the rms level above the measured flux within the beam size (centered at the expected source position),which corresponds to a 99.87% confidence level. This differs slightly from the upper limits in Ansdell et al. (2016), considered to be 3σ.
We also derived the flux densities and image rms for the sources observed with ALMA in Band 7 in the Lupus completion survey (see also van Terwisga et al. 2018van Terwisga et al. , 2019)). For these sources we followed the same procedure as for the BDs disks. As reported by van Terwisga et al. ( 2018 a fluxes ratio of 1.3 ± 0.009 between both observations. Since GQ Lup and the rest of the disks of the Lupus completion survey were observed on the same day, we applied that factor to the measured fluxes of all these sources. The results are included in Table 2. For the observation results from the Lupus disks survey (ALMA Project ID: 2013.1.00220.S), we refer to Ansdell et al. (2016).
Previous work characterizing interferometric observations of protoplanetary disks modeled the continuum emission with either physical or empirical models. A physical model commonly used is the two-layer approximation model (Ricci et al. 2014;Testi et al. 2016;Tazzari et al. 2017). Although this model can successfully describe the spectral energy distribution (SED) of TTS disks (Chiang & Goldreich 1997;Dullemond et al. 2001), the model assumes a simplified physical structure of the disk.
In order to provide an observational characterization of the emission, we prefer to fit empirical analytical functions to the emission profile to allow a more straightforward comparison of the disk properties.
We model the extended disks around stars with the Nuker profile used in Andrews et al. (2018) to characterize the Lupus disks observed with ALMA in Band 7 from Ansdell et al. (2016). By using the Nuker profile we ensure the homogeneity on the characterization of the Lupus disks sample, which is a key aspect of the demographic study discussed in the following section. We follow the modeling described in Tripathi et al. (2017), using the Lauer et al. (1995) formulation:
where ρ t is the transition radius that sets the boundary between the inner and outer regimes of the radial profile, γ and β are the inner and outer disk slopes and α is a factor that determines the smoothness of the transition between both regimes. The disk is assumed to be azimuthally symmetric. The total number of parameters used to model the extended disks with the Nuker profile are 9: ρ t , γ, β, α, the total disk flux density F tot , and four additional geometrical parameters connected to the observation: inclination of the disk in the sky towards the observer (i, 0 • face-on disk, 90 • edge-on), the position angle in the sky plane (PA, defined east of north), and right ascension and declination off-sets to the phase center of the observations (∆RA and ∆Dec).
A simple parametrized Gaussian function has been used to fit the BD disks to reduce the number of free parameters because the emission maps of these objects are extremely compact; only one BD disk is marginally resolved. This function can be used to model moderate-resolution observations to characterize the disk brightness profile and its size, as shown in
where ρ is the projected radius in the sky in arcsec, I 0 is a normalization factor, and σ is the standard deviation of the Gaussian profile. The disk is assumed to be azimuthally symmetric. Six free parameters were used to model the BD disks: two from the Gaussian model (I 0 and σ), together with the observational parameters analogous to the Nuker model (i, PA, ∆RA and ∆Dec).
To perform the fits, the Galario (Tazzari et al. 2018) and emcee (Foreman-Mackey et al. 2013) packages were used. For a detailed description of the methodology we refer to Tazzari et al. (2016Tazzari et al. ( , 2017)). To run the affine invariant Markov chain Monte Carlo (MCMC) from Goodman & Weare (2010), 200 walkers were used in order to investigate the parameter space for each disk ( 20-30 walkers for each parameter); the computation ran for 20 000 steps per walker, which guaranteed convergence in all fitted disks.
An appropriate definition of the disk size is necessary for a proper characterization of disks from observations, and for comparison to theoretical models. An approach commonly used is to extract the radius from the models used to fit the emission. The transition radius ρ t from the Nuker profile definition (Eq. ( 1)) provides misleading information on the disk size, since systems with very different architectures and extents might have similar values of ρ t . When fitting to a power-law with an exponential cut-off, an analogous problem arises if the cut-off radius ρ c from Guilloteau et al. (2011) parametrization is used. The outer radius R out was used in Ricci et al. (2014) and in Testi et al. (2016) to fit BD disk observations, defined as the outermost radius of their modeled surface density. All these definitions may be useful in specific studies, but are not suitable for a general characterization of disk sizes from their emission. A more convenient size definition is the radius enclosing a certain fraction of the total disk emission. This definition with a fraction of 68% of the total disk emission has been used in many recent works (e.g., Tripathi et al. 2017;Andrews et al. 2018;Facchini et al. 2019;Long et al. 2019;Manara et al. 2019). In Tazzari et al. (2017), 95% of the total disk emission was used. To avoid confusion with the different terminology used in the literature, we simply refer to them as 68% (R 68% ) and 95% (R 95% ) flux radii. In the nomenclature used throughout this work, R refers to the radius in the system reference frame (typically in au), while ρ stands for the projected radius in the sky plane in arcsec.
We tested both R 68% and R 95% to determine the quality of each radius as the characteristic size for the disk emission (details in Appendix A). For this test we fitted the same disk to various models and inferred R 68% and R 95% for each model. This test shows that the dispersion on R 68% is much smaller than for R 95% . Thus, we consider R 68% as the most reliable size definition for our sample. Nevertheless, in the modeling results of the disks (Table 3), we include both R 68% and R 95% for completeness.
Additionally, we fitted several disks that were previously modeled in Andrews et al. (2018) in order to test the proper functioning of our modeling tool. This additional test is also included in Appendix A; the resulting ρ 68% from this work and from Andrews et al. (2018) are in very good agreement.
To ease the comparison with the existing surveys of disks around stars in the Lupus clouds (Ansdell et al. 2016), we provide an estimate of the disk dust masses using the simplifying assumption of optically thin emission and using an average temperature of 20 K. We note that these assumptions may lead to underestimation of the disk mass in cases where the emission is optically thick or the average temperature is lower than the assumed value.
For the dust mass determination, assumptions on the dust temperature and opacity are needed. A dependence of the dust temperature with the stellar luminosity was first proposed on theoretical grounds (e.g., Yorke et al. 1993;Sonnhalter et al. 1995). More recently, from radiative transfer modeling of mm observations, Andrews et al. (2013) proposed a single-value mean temperature for each disk that could be used to estimate the disk mass for objects with L ∈ [0.1, 100] L . Soon after, van der Plas et al. ( 2016) suggested a more flattened relation for VLM objects. In Ballering & Eisner (2019), the correlation of disk temperature with the stellar luminosity was derived using simplistic radiative transfer models from SED fitting of Taurus disks. Daemgen et al. (2016) and Tazzari et al. (2017) showed that depending on assumptions on disk size and vertical structure, similar T dust even with very different luminosities are compatible with the data. From the different studies, it is unclear whether or not there is a simple and general relation between dust temperature and stellar properties. Using additional relations of the temperature with other stellar properties as proposed by Andrews et al. (2013), and van der Plas et al. ( 2016) might introduce spurious results in our analysis, or even erase possible relations between different disk properties.
We therefore compute the dust mass of each BD and stellar disk assuming a constant dust opacity of κ 890µm = 2 cm 2 g -1 , following previous ALMA Band 7 observations for VLMs and BDs (Ricci et al. 2014;Testi et al. 2016), and an averaged dust temperature of T dust = 20 K, as in Pascucci et al. (2016), and
The results of the modeling and the derived properties of dust mass and disk size are presented in this section. The results for J154518.5-342125 are shown in Figs. 3-5 and7, while the results for the remaining disks whose fits converged are in Appendix C.
In Fig. 3 we show the real and imaginary part of the observed and modeled visibilities as a function of baseline. The visibilities were first centered using ∆RA and ∆Dec from the model with lowest χ 2 , and were then de-projected taking i and PA (for a detailed description, see Tazzari et al. 2017Tazzari et al. , 2018)).
The posterior probability distribution functions (PDFs) of the free parameters (of the Gaussian model, or from the Nuker profile), are shown in the top panels of plots of the fits in Appendix C). This limitation does not affect the characterization of their disk sizes.
The modeled emission distribution of the disk is shown in the top panel of Fig. 5. The bottom panel shows the normalized cumulative flux f cumul derived from Eq. (A.1) for the respective models of the top panel. In both plots, the inferred values of R 68% and R 95% radii are included as vertical dashed and dotted lines. These radii are computed for each model of the MCMC; the final values of R 68% and R 95% are the median of their respective PDFs, with upper and lower errors as the median ±1σ (example of a R 68% PDF in Fig. 6).
In Fig. 7, we show the reconstructed (observed, modeled, and residuals) continuum emission of the source in the sky plane.
Residuals are at noise level on all fitted disks, indicating that the model represents the observation faithfully.
The results of the free parameters from the fits and the derived disk properties can be found in Table 3. The values shown are the median of their respective PDF. Lower (upper) uncertainties are obtained from the range between median and 16th (84th) percentiles of their posterior distribution. The missing values are for the cases in which the fit did not find a clear convergence. For source Sz 102, from our fit results, gray obscuration due to its inclination (∼60 • ) would not explain its subluminous nature. Strong episodic accretion as suggested by Baraffe & Chabrier (2010) could explain its luminosity: this effect reduces the radius of the star, increasing its temperature and resulting in different pre-MS path and a lower luminosity. Another viable explanation would be a misalignment of the inner and the outer disk.
The radii enclosing 68 and 95% of the total flux (R 68% and R 95% ) are specified in Table 3. The values are inferred from their respective PDFs (example in Fig. 6), as derived from the model parameter results. The radii of the detected BDs disks are unfortunately poorly determined due to the compactness of the sources combined with the low S/N of their continuum emission at this waveband. Only for J154518.5-342125 can we properly quantify its size, since its continuum emission is detected with enough S/N and is marginally resolved. We consider a disk to be marginally resolved if the disk emission is of similar spatial scale to the beam size in the image plane and the observed visibilities can be fitted by a Gaussian function with well constrained σ. For all the other BD sources we provide upper limits of their sizes as 95% confidence level, inferred from the PDF of R 68% and R 95% . On the other hand, we determined the emission distribution size A114, page 7 of 23 A&A 633, A114 (2020) Fig. 7. Observed (left), model (center), and respective (right panel) residuals for the continuum emission of J154518.5-342125 observed with ALMA. The modeled emission map is reconstructed from the synthetic visibilities with the lowest χ 2 from the interferometric modeling. The contours are drawn at increasing (or decreasing) 3σ intervals as solid (dashed) lines.
Table 3. Results from the modeling of the studied disks, together with inferred sizes (R 68% and R 95% , in au), and total disk dust mass (M dust ).
Object log 10 (I 0 )
BDs from this survey J154518.5-342125 10.56 14.4 +3.9 -3.5
1.4 +0.4 -0.3 50.0 +0.9 -0.9
30.8 +1.5
-1.6
69.2 +2.9 -2.9 -0.032 +0.001 19.1 ± 1.9
Notes. The first 9 objects are the full list of known BDs in Lupus, fitted to a Gaussian model. The six free parameters are: normalization factor of the emission profile log I 0 , standard deviation of the Gaussian profile σ, inclination i, position angle PA, and right ascension and declination off-sets to the phase center of the observations ∆RA and ∆Dec. The last 7 disks were fitted to a Nuker profile. The nine free parameters of these fits are: transition radius ρ t , inner and outer slopes γ and β, smoothness parameter α, total disk flux density F tot , and the geometrical parameters of the observation i, PA, ∆RA and ∆Dec.
for six out of seven disks of the Lupus completion survey; for the remaining one (Sz77) we provide an upper limit. Previous ALMA observations of BD disks in other regions showed that most of the objects were too compact to be resolved (van der Plas et al. 2016;Testi et al. 2016). Using different methodology to define and derive disk radii, Ricci et al. (2014) and Testi et al. (2016) showed that some BD disk radii (R) in Taurus may extend beyond R 80 au, while in ρ-Oph BD disks seem to all have R 25 au.
Nevertheless, we should bear in mind that the few BD disks with well determined sizes are among the brightest and most massive of the BD population of their respective regions (Lupus,
Constant (20 K) 1.69 ± 0.19 -3.89 ± 0.13 0.60 ± 0.05 Andrews et al. (2013) 0.95 ± 0.18 -4.08 ± 0.12 0.57 ± 0.05 van der Plas et al. ( 2016) 1.25 ± 0.18 -3.98 ± 0.12 0.58 ± 0.05
F 890 µm [mJy] -1.31 ± 0.17 -0.41 ± 0.27 0.41 ± 0.05 100
Ṁacc [M yr -foot_0 ] Constant (20 K) 0.69 ± 0.14 -7.26 ± 0.36 0.56 ± 0.08 Andrews et al. (2013) 0.63 ± 0.17 Taurus and Ophiuchus); they are likely not representative of the BD population.
The total disk dust mass (last column in Table 3) is computed from the assumptions detailed in Sect. 4.2. As all detections have good S/N (Table 2), the main uncertainty when comparing samples observed at different times is the flux calibrator uncertainty (∼10%, see Sect. 3). Dust mass upper limits for nondetected disks around BDs and stars are obtained from the respective continuum flux upper limits as described in Sect. 3.
The total dust mass for the detected BD disks range between 0.2 and 3.2 M . This means that our sources are within the lightest protoplanetary disks known to date. In particular, SONYC-Lup3-7 is the BD disk with the lowest dust mass estimate, independent of the prescription used for the dust disk mass determination. Comparing our dust mass results of BD disks in Lupus to the results of BD disks in other regions, our results are found to be similar. In Testi et al. (2016), a sample of 17 BD disks in the ρ Ophiuchus region were observed and their dust mass estimates are within 0.5 and 6.3 M , with the same assumptions of dust temperature and opacity as for our results. The dust masses of Taurus disks around BD and VLMs (Ward-Duong et al. 2018) range between ∼0.25 and ∼16.7 M , using the same temperature and opacity values as in this work.
We performed a demographic comparison between the BD and stellar disk populations of Lupus. For this analysis we use the derived disk properties to test whether the known relations for disks around stars are also relevant for BD disks. The observational datasets of both populations were obtained at the same facility (ALMA, Band 7), and the derivation of the properties has been conducted with homogeneous methodology for the entire disk population. In this manner we eradicate systematic errors due to the mixing of diverse datasets handled with different methods.
In addition, we updated the relations between disk properties in Lupus with the largest sample of disks in the region observed with ALMA in Band 7 thanks to the incorporation of the seven stellar disks from the Lupus completion survey to the stellar population. In Table 4 we summarize all the correlations discussed throughout this section.
As a preliminary step, we show in Fig. 8 the relation of the respective observables of M and M dust , that is, the stellar luminosities L and the fluxes at 890 µm wavelength (scaled to a distance of 158.5 pc). From the figure, there is a continuity of the correlation for any range of L , and it holds for the BD population. The linear regression shown in the figure is for the entire population (stars and BDs), obtained following the Bayesian method described in Kelly (2007) 1 . Uncertainties and upper limits are taken into account, while subluminous objects are excluded for the fit. Testi et al. (2016) found potential evidence of BD disks being less massive than stellar disks, based on the analysis of an incomplete sample of BD disks in Ophiuchus. Our Lupus sample allows us to check whether similar results hold in this star forming region. In Fig. 8, there is no obvious trend for BDs to have very significantly smaller 890 µm fluxes than stars with similar luminosities. To quantify this comparison, we followed a similar procedure as in Testi et al. (2016) based on a statistical comparison of the two populations, and analyzed whether the distribution of the M dust /M ratios in the sample of BD disks is consistent with being drawn from the same distribution as for the stars.
Figure 9 shows the cumulative distributions and the histograms of the values of the M dust /M ratios for the Lupus samples. Dust mass of each object is derived following Sect. 4.2, M as described in Sect. 2. The histogram shows that the values of the BD ratios are similar to the stellar population ratios, unlike the Ophiuchus sample in Testi et al. (2016). We performed the Anderson-Darling testfoot_1 to study the null hypothesis that the two samples are drawn from the same underlying population, obtaining a probability of 6% that the BD and stellar disk populations are drawn from the same distribution. Although it is a low percentage, it is below 2σ significance. Moreover, the likelihood increases to ∼80-90% if we use the dependence of T dust with L from Andrews et al. (2013), and van der Plas et al. ( 2016). Thus, the data are consistent with the null hypothesis to be correct. Our analysis of the Lupus sample does not show a statistically different fraction of dust mass around BDs as compared to stars. We caution that our sample of VLM stars and BDs in Lupus is very limited and that the results of Testi et al. (2016) were based on highly incomplete and inhomogeneous samples. Further studies with larger and/or unbiased samples are needed to make firm conclusions on this matter.
The result of the previous analysis is also confirmed by inspecting the dependence of M dust on M (Fig. 10). In Appendix B, this dependence is shown for the other T dust prescriptions. The linear regression result is consistent with those of Ansdell et al. (2016), andPascucci et al. (2016) when using the same assumptions of dust opacity and temperature. The slope (α) and intercept (β) for the stellar population are α = 1.73 ± 0.25 and β = -3.88 ± 0.14 respectively (inferred using 1. linmix package, and including upper limits of ALMA nondetections). As consequence of incorporating the BD population into the fit, there is a substantial reduction of the uncertainty of α and β thanks to the extension of the mass range over one order of magnitude: the slope and intercept become 1.69 ± 0.19 and -3.89 ± 0.13. If we compute a linear regression taking into account only the BD sample, we obtain a slope and intercept that is in agreement with the stellar fit, although the uncertainties in this case are large due to the short range in both axes of the BD population.
The existence of a correlation between the disk luminosity and its size was first shown using pre-ALMA observations (Andrews et al. 2010;Piétu et al. 2014). For the Lupus disk population, this dependence was confirmed in Tazzari et al. (2017) and Andrews et al. (2018). In Fig. 11, we show the updated relation for the Lupus disk population, including the seven new measurements from this paper (one BD disk and six disks around stars). The linear regression shown in the figure is obtained for the stellar disk population, excluding the upper limits of the disks with poorly constrained sizes.
The BD disk with well constrained R 68% (J154518.5-342125) is in very good agreement with the relation for stars. The result suggests that this BD disk is a scaled-down equivalent of the very extended disks around stars that show substructure. Nevertheless, since its central object mass is near the BD/VLM boundary, this result might not be representative of the full BD population. Higher-angular-resolution observations of the BD population are needed in order to obtain reliable estimates of their sizes. The estimated size upper limits of the other BD disks provide limited constraints on the relation. The compactness of the BD disks might be indicative that BD disks follow the size-luminosity relation of stars, as suggested by Hendler et al. (2017) from SED fitting of disk observations, and also from the results for ρ-Ophiuchus of Testi et al. (2016). If these objects were to follow the same relation as stars, their R 68% would range between 1 and 10 au. There is now evidence of optically thick emission in the inner ( 50 au) regions of disks around stars (Huang et al. 2018;Liu 2019;Zhu et al. 2019). Likewise, BD disks in Lupus might be optically thick, as suggested by their compact continuum emission. In Fig. 11, we show two optically thick (optical depth τ >> 1) fiducial models, the first one (green) assuming a constant T dust of 20 K, and a second model (purple) with radial dependence of T dust [K] ≈ 30 × ( L L ) 0.25 × ( R 10 ) -0.5 (Andrews et al. 2013). The emission of these models is described by I ν (R) = F B ν (T dust ), where F is a filling factor that describes the fraction of the disk emission distribution that is optically thin: F = 1 if the disk emission is optically thick, 0 < F < 1 for a partially optically thick disk (analogous to Tripathi et al. 2017;Andrews et al. 2018). The optically thick curves in the figure are built considering a of disks with increasing outer size. Objects laying on the line are compatible with being fully optically thick. Additionally, in optically thick disks, the inferred R 68% trace the location of large grains rather than the physical outer radius of the disk (Rosotti et al. 2019).
The only BD with determined dust disk size (J154518.5-342125) lays below both fiducial models. Its dust emission can be understood as optically thin with a fraction of the disk emission distribution being optically thick. If its disk emission is partially optically thin, a portion of dust is not observed, thus the inferred dust mass is underestimated. The upper limits of the remaining BD disks are far below the optically thick models, although their exact positions in the R 68% -F 890µm are unknown.
A linear correlation between mass accretion rate onto the central object ( Ṁacc ) and the disk mass is expected if disks evolve viscously (e.g., Dullemond et al. 2006;Natta et al. 2007;Lodato et al. 2017). Observational evidence for this correlation was first reported by Manara et al. (2016b) for the Lupus disks, and Mulders et al. (2017) in the Chamaeleon I region. The Ṁacc -M dust relation for Lupus disks is shown in Fig. 12. The x-axis of the figure is an estimate of the total disk mass based on our derivation of the disk dust mass and assuming a gas-todust ratio of 100. The Ṁacc and its uncertainty is taken from the X-shooter observations presented by Alcalá et al. (2014Alcalá et al. ( , 2017)). The accretion rate values have been recomputed with the new accretion luminosities that correspond to the parallaxes from Gaia DR2. One BD disk (SONYC-Lup3-10) was not characterized from X-shooter observations. Although the H α emission line is known (Mužić et al. 2014), we have excluded this BD from the analysis in order to have a fully homogeneous sample for our statistical study.
The linear regression for the stellar population in Fig. 12 has been obtained excluding BDs (due to the different BD disks behavior compared to stellar disks, demonstration below), nondetections from ALMA, upper limits of Ṁacc , and subluminous sources. The resulting slope is α = 0.69 ± 0.14, while the intercept is β = -7.26 ± 0.36. When using the same assumptions of dust temperatures and opacities, the linear regression is consistent with the results presented in Manara et al. (2016b). Brown dwarfs have systematically lower accretion rates than stars for the same disk mass. This is also seen in the relation between the more directly observed properties (in Appendix D), and is independent of the considered prescription of the dust temperature (results using other prescriptions in Appendix B).
We inspected the M disk / Ṁacc ratio for the Lupus disk population to confirm or deny this trend. This ratio can be understood as the accretion depletion timescale (or disk age, as in Jones et al. 2012, see also Rosotti et al. 2017), and provides an estimate of the survival time of the disk, assuming that accretion onto the central object remains constant and that accretion is the dominant mechanism for the depletion of the disk. In Fig. 12, we plot different lines indicating accretion depletion timescales of 0.1, 1, and 10 Myr.
As in Sect. 5.1.1, we conducted a statistical analysis of the two populations in order to confirm the behavior of the BD disks. We compare the BD population with the subsample of stars with disk masses within the range of the BD disk masses (in other words, all disks with log 10 (M disk ) < -3 in Fig. 12). This is done in order to remove the more massive disks from the stellar sample, which may accentuate the difference between populations. Subluminous objects and those with upper limits of the Ṁacc BD disk population has a significantly larger accretion depletion timescale with respect to the stellar population, with only a 0.6% probability that the BD and stellar disk populations are drawn from the same original distribution. The median value of this timescale for the BD sample is 9.5 Myr, while this timescale is 1.4 Myr for the stellar subsample considered. The results hold when using other prescriptions of the disk dust mass for this test, with even lower probabilities (∼0.05%).
The result of the accretion depletion timescale is obtained using the total disk mass, which is estimated assuming that the emission of the dust is optically thin, and that dust mass traces the total disk mass. If submm emission of BD disks is optically thick, the disk masses are underestimated, and consequently the accretion depletion timescale is larger than the estimated values. Thus, the difference on the accretion depletion timescale would be even more pronounced if BD disks were optically thick at these wavelengths.
A larger accretion depletion timescale may reflect a difference in the accretion process of BDs with respect to stars. If viscous evolution models are invoked to explain accretion onto the central star, a weaker accretion would imply a lower α parameter (Shakura & Sunyaev 1973) in BD disks. Since the turbulence in viscous disks depends on α, a lower Ṁacc implies a less turbulent disk. Thus, in viscous disks, our result suggests that the α parameter of disks around BDs is lower than around stars (in contrast with the results of Mulders & Dominik 2012), and consequently BD disks would be less turbulent. A bi-modal behavior of accretion has been suggested observationally in Alcalá et al. (2017), andManara et al. (2017), and predicted from theoretical modeling by Vorobyov & Basu (2009); in those studies the two suggested modes were between VLMs (M < 0.2 M ) and more massive stars. When performing a statistical comparison in our sample between the VLM population (0.1 M < M < 0.2 M ) and more massive stars, we obtain a likelihood in the Anderson-Darling test of ∼9% (averaged over the three different prescriptions of M dust used in this work). Thus, VLM stars might show this behavior as well, but less pronounced, and with much lower statistical significance.
A lower viscosity in disks around BDs/VLM stars compared to disks around more massive stars could be explained with a globally lower ionization rate (e.g., see Mohanty et al. 2005). As a consequence of the general correlation between X-ray and bolometric luminosities, young BDs/VLM stars have slightly lower X-ray luminosities than more massive young low-mass stars (see, e.g., Gregory et al. 2016); this might lead to slightly lower ionization rates in the disks of BDs and VLM objects. Some BD/VLM disks being flatter than disks around more massive stars, which would decrease the irradiation cross-section, could also explain a lower ionization rate (as suggested from SED models and observations by Ercolano et al. 2009;Pascucci et al. 2003;Apai et al. 2004;Allers et al. 2006). However, there is evidence of disks around BDs being flared (e.g., Natta & Testi 2001;Natta et al. 2002;Furlan et al. 2011), and therefore this last possibility would need more detailed and extended investigation.
The exoplanetary systems recently discovered around BDs and VLMs can be used to study the ability of BDs to form planets. Since the planets hosted by Trappist-1 and Proxima Centauri are most likely of rocky composition, the total planetary mass of these systems can be compared to the estimates of the disk dust mass of the BD disk population in Lupus. The total mass of the seven known planets (Gillon et al. 2017) in Trappist is 4 M (Wang et al. 2017). Proxima B, the planet hosted by our closest neighbor (Anglada-Escudé et al. 2016), has a minimum mass of 1.3 M . In (Bixel & Apai 2017) the planet mass was predicted to be 1.63 +1.66 -0.72 M with 95% confidence level. The remaining known exoplanets around BDs have estimated masses of at least several times that of Earth. A considerable fraction of them have been detected via microlensing (e.g., Jung et al. 2018), with typically much higher estimated masses. Thus the picture for these planets is analogous to the Trappist-1 planets and Proxima B.
From theoretical predictions of planet formation around BDs via core accretion (Payne & Lodato 2007), disk masses (gas and dust) on the order of a few Jupiter masses are required in order to form Earth-like planets around BDs. Not only do none of the BD disks in Lupus (this paper) and ρ-Oph (Testi et al. 2016) have enough mass available at their current stage to form a planetary system, but even the available mass in solids is smaller than the total planetary mass in Trappist-1. The efficiency with which the available mass is converted into the final planetary rocky cores might be boosted by internal recycling of the disk material, but it is unlikely to reach an efficiency close to unity (Manara et al. 2018, and references therein). On the other hand, the tentative result from Sect. 5.1.3 of lower viscosity and ionization rates on BD disks might contribute to the presence of an extended dead zone in the disk, which would boost the planet formation process.
A plausible explanation to alleviate this divergence is that the determination of dust mass from continuum emission flux might be underestimated, as pointed out in Ballering & Eisner (2019). This might be the case if the emission at this wavelength is optically thick; consequently the inferred M dust provides only a lower limit of the disk dust mass. In Fig. 11, disks laying on the τ >> 1 fiducial models are consistent with their emission being fully optically thick. The only BD disk with a well-determined size (J154518.5-342125) is below these models. This suggests that the emission of this BD disk is optically thin with small regions of the disk being optically thick. The inferred dust mass of this disk is underestimated by an unknown fraction. This can help to explain the mass difference with exoplanetary systems. Nevertheless, it seems unlikely that partial optically thin emission alone can account for this large difference in solid masses.
A likely possibility is that planets might have already formed at this stage of disk evolution (Greaves & Rice 2010;Najita & Kenyon 2014;Manara et al. 2018;Dodds et al. 2015). If this is indeed the case, the formation of planetary rocky cores would have occurred within the first million years (considering the estimated ages of Lupus and ρ-Ophiuchus). While direct detection of planets embedded in protoplanetary disks is extremely difficult (see Sanchis et al. 2020;Johns-Krull et al. 2016), the presence of circumplanetary disks, if confirmed, sets a strong indirect evidence of young planets in these disks (Keppler et al. 2018;Isella et al. 2019;Pérez et al. 2019). Analysis of the gas kinematics can also be used as an indirect method to study embedded planets (Teague et al. 2018;Pinte et al. 2018). Other indirect indications, such as the existence of gaps, spirals, asymmetries, and dust processing are observed frequently, and suggest that planets might already have formed (e.g., ALMA Partnership et al. 2015;Zhang et al. 2018;Lodato et al. 2019;Pinilla et al. 2018).
In this work we presented new Band 7 ALMA observations of five protoplanetary disks around BDs in Lupus. Combined with previous observations, we analyzed the submm disk properties of the known population of BDs and VLM objects with infrared excess. From the continuum fluxes and modeling the visibilities, we inferred total dust disk masses and characteristic sizes of the disk population. Due to the extremely compact emission of the BD disks in Lupus, the size determination was only possible on one BD disk, while for the other detected disks we present upper limits on the size.
We updated the relations of M -M dust , size-luminosity, and M dust -Ṁacc relations extending them down to the substellar regime. Brown dwarf disks in Lupus follow the relation for stars between stellar mass and dust disk mass. They show no statistical difference from the stellar disk population on the disk mass fraction, however we note the apparent lack of massive BD disks. On the other hand, the accretion depletion timescale (inferred assuming that dust mm continuum emission is a reliable proxy of the total disk mass) of the BD population is significantly longer than for stars (9.5 Myr vs. 1.4 Myr), which in viscously evolving disks may imply a lower α value, possibly linked to a globally lower ionization rate. Lastly, we inspected the ability of these objects to form planets, comparing the estimated disk dust masses with the rocky planetary masses in known exoplanetary systems. The estimated disk dust masses around brown dwarfs are very low, suggesting that either these systems are unable to form planets, or more likely, that rocky planetary cores have already formed within the first million years. Optically thick emission in BD disks can alleviate this mass discrepancy.
Implemented with the linmix Python package, https://linmix. readthedocs.io/en/latest/index.html
Using scipy.stats Python module, https://docs.scipy.org/ doc/scipy/reference/stats.html