© N. T. Kurtovic et al. 2021

1 Introduction

For every ten stars that are formed in the Milky Way, around two to five brown dwarfs (BDs) also form (e.g., Scholz et al. 2012; Mužić et al. 2019), and M-dwarfs represent about three-quarters of all the stars in our galaxy. Exoplanet discoveries show that short-period (< 50 days) sub-Neptune planets occur more frequently around M-dwarfs than around FGK stars (Mulders et al. 2015), but also a few giant planets have been discovered around BDs and very low mass stars (VLMS ≲ 0.1 M_⊙, e.g., Morales et al. 2019). This implies that planets of a large range of masses can form around these objects, although it remains an open question if these massive objects form as binary companions of the BDs and very low mass stars (VLMS), or as planets. Current models of planet formation through pebble or planetesimal accretion cannot explain the formation of giant planets around VLMS (Liu et al. 2020).

Observations of BDs and VLMS from the near-infrared to centimeter wavelength reveal the existence of circumstellar disks around these objects (e.g., Luhman 2006; Klein et al. 2003; Scholz et al. 2006, 2007; Ricci et al. 2012, 2017a,b; Daemgen et al. 2016; van der Plas et al. 2016; Sanchis et al. 2020), which are more compact and lower in dust mass when compared to disks around T-Tauri stars (e.g., Pinilla et al. 2017a; Ward-Duong et al. 2018; Hendler et al. 2017, 2020). The typical millimeter fluxes of such disks suggest that they have dust masses of few Earth masses, challenging the formation of giant planets through core or pebble accretion (Liu et al. 2020).

The core accretion scenario for planet formation assumes collisional growth from sub-μm-sized dust particles from the interstellar medium (ISM) to kilometer-sized bodies or planetesimals (e.g. Pollack et al. 1996). The collisions of particles and their dynamics within the disk are regulated by the interaction with the surrounding gas. Different physical processes lead to collisions of particles and their potential growth, such as Brownian motion, turbulence, dust settling, and radial drift (e.g., Brauer et al. 2008). All of these processes have a direct or indirect dependency on the properties of the hosting star, such as the temperature and mass. For instance, from theoretical calculations, settling and radial drift are expected to be more efficient in disks around VLMS and BDs, with BD disks being 15–20% flatter and with radial drift velocities being twice as high or even more in these disks compared to T-Tauri disks (Mulders & Dominik 2012; Pinilla et al. 2013).

With radial drift being a more pronounced problem in disks around BDs and VLMS, it is still unknown how this barrier of planet formation is overcome in these environments where the disks are more compact, colder, and have a lower mass. Millimeter-sized particles have been detected in BD and VLMS disks through measurements of the spectral index (Ricci et al. 2014; Pinilla et al. 2017a), which are only possible to explain when radial drift is significantly reduced by the presence of strong pressure bumps (Pinilla et al. 2013). The presence of pressure bumps produces substructures, such as rings, gaps, spiral arms, and lopsided asymmetries, with a different amplitude, contrasts, and locations depending on the origin of the pressure variations (e.g., Pinilla & Youdin 2017; Andrews 2020). Currently, due to sensitivity limitations, most of our observational knowledge about substructures comes from bright (and probably massive) disks, such as the DSHARP sample (Andrews et al. 2018a). A less biased sample of ALMA observations of disks in the star-formation region of Taurus has demonstrated that at least 33% of disks host substructures at a resolution of 0.1″, and the disks that do not have any substructures are compact (dust disk radii lower that ~ 50 au, Long et al. 2018, 2019). It remains an open question if compact disks are small because they lack pressure bumps or because current observations lack the resolution to detect rings and gaps in these disks (the scale of a radial pressure bump cannot be smaller than one local scale height, e.g., Dullemond et al. 2018).

Pinilla et al. (2018a) have thus far reported the lowest mass star with a resolved large dust cavity (radius ~ 20 au) in the disk around the M 4.5 star CIDA 1. In the context of planets creating such a cavity, it has been shown that a high planet-to-stellar mass ratio is needed to open a gap and trap particles in disks around VLMS, because in these cases, the disk scale height at a given location is higher than in moderate or high mass stars (Pinilla et al. 2017a; Sinclair et al. 2020). In a typical disk around a VLMS as CIDA 1, at least a Saturn-mass planet is needed to open a gap in the disk. This challenges our current understanding of substructures and the common idea that planets are responsible for their formation, since these disks around VLMS and BDs may not have enough mass to form such massive planets, although they may form from gravitational instability if the disks were much denser in their early stages (e.g., Mercer & Stamatellos 2020).

Based on the previous CIDA 1 observations, we selected a sample of five disks to observe with ALMA at a resolution of 0.1″ in the Taurus star-forming region, whose properties are similar to CIDA 1. Specifically, these disks around low mass stars are more massive compared to other disks with hosts in the same stellar regime. These observations included ¹² CO and ¹³CO and aim to estimate how common substructures are around VLMS. In addition, as radial drift is expected to be very efficient in these disks, they are excellent laboratories to search for the difference between the radial extent of the gas and dust in disks, which can be a direct signature of radial drift (e.g., Trapman et al. 2019).

This paper is organized as follows: Sect. 2 summarizes the target selection, the ALMA observations, and the calibration of the data. Section 3 presents the modeling of the data in the visibility plane for the continuum as well as in the image plane for the ¹²CO and ¹³CO emission. In Sect. 4, we discuss our results in the context of different origins for the seen substructures as well as the observed difference between the radial extent of the gas and dust. Section 5 summarizes the conclusions of this paper.

2 Target selection and observations

2.1 Target selection

For the new ALMA observations, we selected five disks around VLMS in the Taurus star-forming region. Given that this is the first small survey of high angular resolution observations of disks around VLMS, the sample was selected to optimize the chances of finding substructures. The criteria used to select the targets were as follows: (1) the target has been previously observed and detected by either SMA or ALMA in millimeter wavelengths (based on the observationsby Andrews et al. 2013 and Ward-Duong et al. 2018, respectively); (2) it has a stellar mass between ~ 0.1 and 0.2M_⊙; and (3) it has a high disk dust mass compared to the stellar mass. The last condition comes from observations of the M_dust−M_⋆ relation of transition disks and disks with substructures (Pinilla et al. 2020), which shows that those disks usually have higher dust masses compared to others with stellar hosts of a similar mass. From the list of targets that fulfilled the conditions, we selected the ones with the lowest optical extinction.

2.2 Observations

For sample completeness, we included the archival data of CIDA 1 in our list of targets, given that combining these previous archival observations results in a dataset with a similar sensitivity and angular resolution. For consistency with future works, the target properties used in this study (shown in Table 1) were taken from Herczeg & Hillenbrand (2014), with stellar masses derived following the method described in Pascucci et al. (2016), using distances inferred from Gaia DR2 (Gaia Collaboration 2018). The distances in parsecs was calculated as the inverse of the parallax. The spatial distribution of our sample in the Taurus optical extinction map is shown in Fig. A.1.

These targets were observed with ALMA at 0.87 mm (Band 7) as part of the project 2018.1.00310.S (PI: P. Pinilla), during Cycle 6, with the spectral setup configured to observe in four spectral windows centered at 331.3, 333.3, 344.0, and 345.8 GHz, with two of them centered to observe dust continuum emission, and two observing molecular line emission from ¹² CO (J = 3−2) and ¹³CO (J = 3−2). The frequency resolution for ¹²CO was 244.1 kHz per channel; while for ¹³CO and the continuum, it was 976.6 kHz. The most extended antenna configuration used was C43-8, providing an angular resolution of 0.08″ at best. Some of our sources had archival data from ALMA project 2012.1.00743.S (PI: G. van der Plas) and 2016.1.01511.S (PI: J. Patience), with observationsof ¹²CO and dust continuum, which were also included in the self-calibration and analysis. Archival Band 7 data of CIDA 1 was observed by the project 2015.1.00934.S (PI: L. Ricci) published in Pinilla et al. (2018a), and in this work we combined it with 2016.1.01511.S (PI: J. Patience). A summary of the observation details and the data considered for each target are shown in Table A.1.

The raw datasets were calibrated by applying the ALMA pipeline using the CASA version specified for each project(McMullin et al. 2007). Then, we used CASA v5.6.2 for the subsequent data handling and imaging. We extracted the dust continuum emission of every source by flagging the channels closer than 25 km s⁻¹ to the targeted molecular lines. To reduce the data volume, we averaged over time (6 s intervals) and channels (with a width of 125 MHz). Before combining all available observations for each source, the centroid spatial position of the emission was determined by fitting a Gaussian using the imfit task, and shifted using fixvis and fixplanets tasks to the centroid of the observations of extended baselines, shown in Table A.2. We also checked for a consistent flux calibration by comparing the amplitude of the emission in different executions. We found a discrepancy of 12% in the fourth compact observation of J0433 (2018-11-24), which we rescaled to match all the others.

In order to boost the signal-to-noise ratio (S/N) of each source, we performed self-calibration of the datasets in two steps: first we combined the compact configuration observations and performed self-calibration. Second, we combined those self-calibrated observations with the extended configuration observations and then self-calibrated again. Phase calibrations were applied until the improvement on the S/N was below ~ 5%. Only 1 amplitude calibration was applied in each step. The overall S/N improvement was between 1.5 ~ 4.0, depending on the source. The only source where self-calibration was not possible was J0415, because the initial S/N of 9 was too low for improvements. The final continuum images were generated using a Briggs robust parameter of 0.0 for CIDA 1 and 0.5 for the remaining sources. The image properties are summarized in Table A.2.

All the steps for the dust continuum emission calibrations, including centroid shifting, flux calibration, and self-calibration tables, were then applied to the molecular line emission channels. The continuum emission was subtracted using the uvcontsub task, and the images were generated using a robust parameter of 1. In MHO 6, J0420, and J0415, UV-tappering was applied in order to increase the S/N of the gas images. For MHO 6, we applied a UV-tappering of 0.13″ on the ¹² CO; while for CIDA 7 and J0415, we applied a UV-tappering of 0.1″ in both molecular line images. All our scripts of self-calibration and imaging are available online¹.

The final images of the dust continuum, moment 0, and moment 1 of the ¹² CO and ¹³CO (when detected) are shown in Fig. 1, and the velocity channel maps are in Appendix B. The details of the dust continuum and CO images can be found in the appendix as well, summarized in Tables A.2 and A.3, respectively.

3 Results

3.1 Morphology of the continuum emission

Our sources are spatially compact (radius of 0.1″ ~ 0.3″), and both their radial extent and substructures have sizes comparable to the synthesized beam shape. To avoid image reconstruction biases, we fit the deprojected brightness profile of the sources in the UV-plane. The continuum visibilities were extracted from the self-calibrated measurement set, and we used the central frequency of each channel to convert the UV-coordinates to wavelength units. We started by modeling every source with a central Gaussian profile, and then we increased the complexity of the profile if the residuals suggested it. We also guided our parametrization of the profiles based on the best fitting from frank (Jennings et al. 2020), which fits a nonparametric 1D radial brightness profile in the visibilities, using Gaussian processes. For CIDA 7, J0420, and J0415, the function that describes their brightness profile is a centrally peaked Gaussian profile, following an intensity given by: $$I_{\text{g}}(r)=f_1\exp \left(-\frac{r^2}{2{\sigma }_1^2}\right)\text{,}$$(1)

where I_g is the Gaussian intensity profile of the source as a function of the radius r.

For CIDA 1 and MHO 6, we modeled the disk with a radially asymmetric Gaussian ring or a broken Gaussian from hereafter, that is to say the inner and outer width of the ring can differ. This profile is motivated by results of radially asymmetric accumulation of particles in pressure bumps (see e.g., Pinilla et al. 2015, 2017b). Such radially broken Gaussian profiles have been used to describe the morphology of different rings in transition disks and disks with substructures (e.g., Pinilla et al. 2018b; Huang et al. 2020), which is the same model used in Pinilla et al. (2018a) to model CIDA 1. The intensity profile is given by a ring as follows: $$I_{\text{bg}}(r)=\text{left{nn begin}aligned{f}_1\exp \left(-\frac{{(r-r_1)}^2}{2{\sigma }_1^2}\right)\text{for }r\le r_1\text{}{f}_1\exp \left(-\frac{{(r-r_1)}^2}{2{\sigma }_2^2}\right)\text{for }r>r_1\text{nn end}aligned\text{ right}\text{.,}$$(2)

where I_bg is the broken Gaussian intensity profile as a function of the radius, r₁ is the radial location of the ring peak intensity, and σ_1,2 are the Gaussian widths for the inner and outer sides of the ring, respectively.

Finally, for J0433, the profile is the sum of a centrally peaked Gaussian profile and two symmetric Gaussian rings, as suggested by frank. It is also the profile that creates the lowest amount of residuals from our experiments, such as the single Gaussian, Gaussian ring, and broken Gaussian ring. The intensity profile is $$I_{\text{J0433}}(r)={\sum }_{i=1}^3{f}_i\exp \left(-\frac{{(r-r_i)}^2}{2{\sigma }_i^2}\right)\text{,}$$(3)

where r_{i = 1} = 0, so the first Gaussian is peaked at the center.

The visibilities of each profile were computed by combining each model with a spatial offset (δ_RA, δ_Dec), inclination (inc), and positionangle (PA). Therefore, each model has four more free parameters in addition to those that describe the intensity profile. The Fourier transform and the χ² calculation were carried out with the galario package (Tazzari et al. 2018). The pixel size used in the models is 1 mas, which is several times smaller than the smallest resolvable scale of the observations. The χ² was scaled up by a factor of 2.667 since CASA does not account for the effective channel width, introduced by Hanning smoothing, when it averages the weights during data binning. We adopted a uniform prior probability distribution over a wide parameter range, such that walkers would only be initially restricted by geometric considerations (inc ∈ [0, 90], PA ∈ [0, 180], σ ≥ 0).

We used a Monte Carlo Markov chain (MCMC) routine based on the emcee package (Foreman-Mackey et al. 2013) to sample the posterior probability distribution of each parameter space. Furthermore, we ran more than 250 000 steps after converging to find the most likely set of parameters and the error bars, taken from the 16th and 84th percentile.

Our results for each parameter are shown in Table 2, while in Fig. 2 we show the models with and withoutconvolution (right most panel). The visibilities and radial profile were deprojected using the best inclination and position angle. The residual image was generated in CASA using the same parameters and procedure used for the observations, from a measurement set with its visibilities calculated by subtracting the best model from the data.

The radial profile recovered from the UV-modeling was used to measure the dust continuum emission radii (R_dust) that encloses 68 and 90% of the total flux (the dust R₆₈ and R₉₀, respectively). We used thedifferent sets of parameters sampled by the walkers to compute their radial profiles, and we found the 16th and 84th percentile on the R₆₈ and R₉₀ to calculate the error bars. The results are shown in Table 3 and in Fig. 3. The continuum emission in our sample extends up to 46 au in the biggest disks, MHO 6, and J0433, while we also have the smallest (13 au) and the dimmest (~1 mJy) Taurus disks ever resolved in 0.87 mm emission, CIDA 7, and J0415, respectively.

Fig. 1 ALMA observations of the disks around VLMS in Taurus. From left to right: dust continuum emission, 12 CO moments 0 and 1, and 13CO moments 0 and 1. All boxes are 3.0″ in size. The scale bar represents 20 au, and the white ellipses show the synthesized beam. The centers, beam sizes, measured fluxes, and sensitivity are detailed in Tables A.2 and A.3. Dashed lines show the region used to calculate the radial profiles and to measure R68 and R90. A zoomed-in version of the continuum images can be found in Fig. 2.

3.2 Rings and cavities in CIDA 1, MHO 6, and J0433

We find evidence of ring structures in three out of six disks in our sample. At the current resolution, CIDA 1 and MHO 6 show a single ring and a cavity, which are well described by a broken Gaussian profile. In CIDA 1, we find the inner side of the Gaussian to be more extended than the outer side, which is similar to the results obtained in Pinilla et al. (2018a). The ratio between the widths of the inner side and the outer side in our analysis is 1.2, and the peak of our ring model is located at 20.8 au.

For MHO 6, we find that a strong asymmetric ring peaked at 10 au, with its inner side having converged to σ = 0 au with a narrow error margin. To ensure that our result was not being affected by numerical biases related to the pixel size, we ran another MCMC with a pixel size of 0.4 mas (about 0.057 au), and we consistently recovered the same result for each parameter. Even though this steep transition could suggest an unresolved inner side of the ring, the best model is also driven by the non-axisymmetric emission of the disk, which is only noticeable when looking at the residuals (see Fig. 2). The contrast of these asymmetries is about 5% of the peak amplitude of the ring, meaning that most of the emission is still well described by a radially axisymmetric ring.

In J0433, our best model finds two rings located at 21 and 32 au, with gaps at 16 and 25 au. The brightness ratio between the first ring and first gap is about 4.4; whereas, the contrast is 1.2 between the second gap and ring. Since the brightness ratio fades at the 1σ error (as seen in Fig. 3), we cannot exclude that it is an extension of the outer side of the first ring.

3.3 Dust disk masses

Assuming the dust continuum emission is optically thin, we estimated the disk dust mass by following Hildebrand (1983): $$M_{\text{dust}}\approx \frac{d^2{F}_{\nu }}{{\kappa }_{\nu }{B}_{\nu }(T)}\text{,}$$(4)

where d is the distance in parsecs (given in Table 1, taken from Gaia DR2), κ_ν is the mass absorption coefficient, and B_ν(T) is the Planck function. For κ_ν, we assumed the opacity law ${\kappa }_{\nu }=2.3{\text{cm}}^2{\text{g}}^{-1}{(\nu /230\text{GHz})}^{0.4}$ (Andrews et al. 2013), while the temperature was assumed constant at T_dust = 20 K (e.g., Ansdell et al. 2016; Pinilla et al. 2018a). We estimated the uncertainty of the dust mass by taking the 10% uncertainty in the flux calibration. Our results for each source are compiled in Table A.2 and shown compared with the disk dust mass in the Taurus region in the Fig. 4. The lowest dust mass estimate was obtained for J0415, with only 0.22 ± 0.02 M_⊕ (about 2.1 ± 0.2 M_Mars), while the most massive dust disk was detected in J0433 with 14.32 ± 1.43 M_⊕.

3.4 ¹²CO and ¹³CO emission

3.4.1 Lines detection in each source

Both molecular lines were detected in all of our sources, except for J0415, where only ¹² CO was observed. As shown in Figs. 1 and 2, our angular resolution of ~ 90 mas was just enough to resolve the dust continuum and gas emission of this source, but the limited sensitivity did not allow us to constrain the geometric parameters as in the other systems.

For the two smooth sources, CIDA 7 and J0420, we found strong cloud contamination, so that the western side of the ¹² CO emission in J0420 is completely absorbed. The rotational pattern was, however, recovered in the ¹³CO (see Figs. 1 and B.9). On the other hand, in CIDA 7, we observed anextended asymmetric emission in the south region, with a velocity range of at least 0.8 km s⁻¹ (it was detected in two velocity channels). The contribution of this emission is not significant to the total flux of the gas emission, and it was not detected in ¹³CO. Although the S/N did not allow us to accurately recover geometric parameters from the gas emission, we can see that for both sources, the PA obtained from the dust continuum emission is in good agreement with the orientation of the major axis in the rotation pattern.

In the sources with detected substructures, both J0433 and CIDA 1 are affected by cloud contamination. For CIDA 1,the side most affected by extinction is the north side, while the same is true for the south side of J0433 (row 1 and 3 of Fig. 1). The least cloud-contaminated source is MHO 6 disk (also seen in Fig. A.1), which is the brightest and more extended disk in CO emission of our sample. Since the cloud contamination is weaker, the gas can be detected continuously in every channel from both high velocity ends. Given our good spatial detection of both lines, we used the package eddy (Teague 2019) to analyze the Keplerian rotation, which we further describe and discuss in Sects. 3.4.3 and 4.3.

Fig. 2 Visibility modeling versus observations of our sample. From left to right: (1) real part of the visibilities after centering and deprojecting the data versus the best fit model of the continuum data, (2) continuum emission of our sources where the scale bar represents a 10 au distance, (3) model image, (4) residual map (observations minus model), (5) and normalized, azimuthally-averaged radial profile calculated from the beam convolved images in comparison with the model without convolution (purple solid line) and after convolution (red solid line). In the right most plots, the gray scale shows the beam major axis.

Fig. 3 Normalized azimuthally averaged radial profiles of the dust continuum, 12 CO, and 13CO of each source. Each curve has a 68% error region. The orange and red solid curves correspond to the average intensity profile for the gas profiles of the deprojected images, while the blue solid line corresponds to the best χ2 solution from the visibility fit of the dust continuum. The dashed vertical line denotes the position and 1σ error of the R90 radius for the dust (blue) and the gas (red).

Fig. 4 Dust mass compared to stellar mass for the Taurus disks. We note that Md was computed using distances inferred from Gaia DR2, assuming a 20K disk temperature and using the most updated value of F0.87mm from Andrews et al. (2013), Tripathi et al. (2017), and Ward-Duong et al. (2018). Sources where F0.87mm was extrapolated from a F1.3mm measurement are plotted in faded colors (see Table 2 in Andrews et al. 2013), while direct measurements from SMA or ALMA are plotted in darker colors. Green triangles denotes upper limits, and the VLMS studied in this work are plotted as orange stars.

3.4.2 Radial gas profiles

The inclinations and position angles obtained from the continuum fitting were used to deproject the distances from the central star in the moment 0 images, and we calculated the azimuthally averaged radial profiles of the ¹² CO and ¹³CO emission. The vertical structure was neglected in this calculation. In CIDA 1, J0433, and J0420, the profiles were calculated from the side that is less affected by cloud contamination. In CIDA 7, we found that the contribution of the asymmetric emission in the south is negligible, but we nevertheless did not take it into account as we wanted to recover the disk emission profile (see masking in Fig. 1). The gas R₆₈ and R₉₀ radii that encloses 68 and 90% of the flux are shown in Table 3 and Fig. 3.

As we discuss in Sect. 4.4, the brightness in dust continuum and gas emission of J0415 are much lower than expected from previous SMA observations (Andrews et al. 2013). The low S/N of the detection prevented us from applying self-calibration, and so the sensitivity is two times worse than in the other VLMS disks. Therefore, we calculated the radial profile of the ¹²CO emission from an image generated by only considering the channels in the velocity range with line detection, and we did not apply clipping at 3σ.

3.4.3 Keplerian rotation of CO emission

Although we recovered the rotational pattern for all the disks in our sample in both CO isotopologues (with the exception of the ¹³CO in J0415), the strong cloud contamination and low S/N of our images prevented us from reliably recovering the dynamical mass from our sources without proper modeling. The only system where the cloud contamination does not completely make the central velocity channels extinct is MHO 6, where we have good spatial detection of both lines. As described in Sect. 2.2, we used CASA to generate the moment 1 image from the velocity maps; in addition, we generated velocity integrated images using the python package bettermoments (Teague & Foreman-Mackey 2018), with a quadratic and Gaussian method, and also by varying the root mean squre (RMS) clipping limit. The Keplerian rotation of these images was modeled using the package eddy (Teague 2019), with different models considering the stellar mass (M_⋆), the central velocity (VLSR), the flaring parameter (ψ), the position angle (PA), and the source center (x₀, y₀) as fixed or free parameters.

Depending on the free parameters used (e.g., fixing the PA or allowing the fitting of vertical structure) and also on the velocity integrated image used, the stellar mass recovered would vary in the range of 0.16~ 0.24 M_⊙, which is consistent with a stellar mass derived from evolutionary models. However, it is important to note that all the models and images would leave residuals, which span two times the channel velocity widths and are also strongly structured, as shown in Fig. B.2.

To have a referential value for the stellar host mass, we used CASA5.6 to get the position velocity diagrams (PV diagram from hereafter) along the major axis of each source, based on the position angle obtained from the continuum UV-modeling. The only exception is J0415, whose PV diagram was obtained along the east-west axis. The PV diagrams are shown in Figs. B.3 and B.4 for the sources with continuum substructures and smooth profiles, respectively.

4 Discussion

4.1 Detection of substructures

We only detected obvious substructures in the brightest disks of our sample (CIDA 1, MHO 6, and J0433), which are also the most radially extended disks in gas and dust emission. The existence of strong dust traps located at larger radial distances from the star is most likely the reason for this observational result, as the dust is allowed to stay for longer timescales in the outer disk, thereby increasing the emitting area. On the other hand, our limited angular resolution only allowed us to detect substructures in the most extended disks.

This direct detection of substructures in 50% of our sample does not represent the occurrence rate of substructures in all disks around VLMS, as our sample is biased towards the brightest disks and our spatial resolution is limited. Therefore, our observations only allowed us to directly detect deviations from a simple Gaussian profile in the disks where the extent of the emission is consistently larger than the synthesized beam size.

CIDA 7 and J0420 are a good example of the limitations that our datasets have when detecting substructures. Even though we are unable to confirm the existence of dust substructures in these systems, the UV-plot in Fig. 2 of CIDA 7 shows some structure which is not described by a Gaussian profile, while the residual image of J0420 suggests that there might be more substructured dust emission that is not described by a single centrally peaked Gaussian in these disks. We lack the resolution and sensitivity to confidently characterize them. Given that these two disks are not totally dust depleted, the expected efficient radial drift must have been counteracted by a dust trapping mechanism. The compactness of these disks suggest that the dust trap is located so close to the star that our resolution did not allow us to detect it. In theory, any pressure bumps cannot be smaller than the local pressure scale height, which implies that if they are located closer to their star, their radial extent is smaller than our current resolution. CIDA 7 and J0420 are good candidates to be targeted by deep observations with ALMA in the most extended antenna configuration, allowing us to detect substructures of ≤ 2 au in size at the distance of these targets, which is six and ten times smaller than the dust R₆₈ of those disks, respectively. Future observations will test if even the very small disks around VLMS are able to generate dust traps, as it is observed in the massive and extended ones.

For J0415, our UV-coverage and sensitivity resolves the emission, and the centrally peaked Gaussian model does not leave any significant residual. Higher sensitivity is needed in order to discern deviations from the Gaussian profile (see Fig. 2).

4.2 Origin of dust continuum rings and cavities

All the detected substructures in our sample resemble ring-like emission. MHO 6 is the only disk displaying what could be a hint of non-axisymmetric dust emission in the residuals at our current sensitivity. Rings structures are the most common type of substructure in moderate and massive stars, as shown by surveys such as DSHARP (Andrews et al. 2018a; Huang et al. 2018) and the Taurus survey at 1.3 mm (Long et al. 2018), and now a similar trend is found for the most extended disks around VLMS.

4.2.1 Detected structures coincide with possible CO iceline location

The iceline of each volatile marks the location in the disk where that volatile transitions from being mostly gas-phase to being frozen out on dust grains. It is possible that this phenomena could induce ring-like substructures in the dust continuum emission by changing the dust opacity and grain collisional fragmentation and growth rates (e.g., Okuzumi et al. 2016).

To investigate if any of the iceline locations of the major volatiles coincide with the location of our gaps, we calculated the midplane temperature of our disks by following Kenyon & Hartmann (1987), where, for an irradiated flared disk, the temperature was parametrized as: $$T(r)=T_{\star }{\left(\frac{R_{\star }}{r}\right)}^{1/2}{\varphi }_{\text{fl}}^{1/4}\text{,}$$(5)

with r is the distance from the star, T_⋆ and R_⋆ are the star temperature and radius, respectively, and ϕ_fl is the flaring angle, which we assume to be equal to 0.05 (Dullemond & Dominik 2004). The stellar radius R_⋆ was measured from the stellar luminosity L_⋆ by assuming black body emission and spherical symmetry ($L_{\star }=4\pi R_{\star }^2{\sigma }_{\text{sb}}{T}_{\star }^4$, with σ_sb the Stefan-Boltzmann constant). As a first approximation, we ignored the contribution of accretion in our calculations of the stellar luminosity (further discussed in Long et al. 2018). If the stellar radius is replaced in Eq. (5), we can obtain the distance r_i at which the temperature T_i is reached, given a star of fixed luminosity L_⋆, following: $$r_i=\frac{L_{\star }^{1/2}}{T_i^2}\cdot {\left(\frac{{\varphi }_{\text{fl}}}{4\pi {\sigma }_{\text{sb}}}\right)}^{1/2}\text{.}$$(6)

We considered the two coldest icelines presented in Zhang et al. (2015), which are CO with a sublimation temperature of ~ 20~ 28 K (we took the lower limit 20 K from Öberg et al. 2011 for direct comparison with Long et al. 2018), and the N₂ iceline which goes from ~12 ~ 15K. The iceline location was calculated for each star using the parameters shown in Table 1, and the results are shown in Fig. 5.

In the right panel of Fig. 5, the positions of the detected substructures are displayed as a function of the root of stellar luminosity. Most of the substructures detected lie within the region where the CO iceline is identified. However, the role that the COiceline plays in the morphology of our disk is not conclusive from our datasets. Even if we neglect the uncertainties introduced by the stellar parameters and the disk temperature, it is not clear if the structures that were generated by an iceline would result in a ring peaked at the iceline position or in a gap (Pinilla et al. 2017c; van der Marel et al. 2018).

Pinilla et al. (2017c) find that icelines do not strongly change the gas density around the icelines, nor can they carve a dust cavity as the ones we detected in MHO 6 and CIDA 1. Given that the peak of the rings detected in these disks are located very close to the position of the CO iceline, we cannot rule out the possibility that the iceline played a role in triggering the mechanism that is carving the cavity. In J0433, we find both peaks and gaps at the region (or close) where the CO iceline could be located. Given the radial compactness of the disks in our sample, the limited angular resolution, and the wide radial range covered by the CO iceline location, it is also possible that the substructures coincide with the iceline just by chance, in the same way that other surveys of substructures have not found a strong correlation between the position of substructures and iceline location (e.g., Long et al. 2018). Additional deeper observations of the CO isotopologues could allow a better constrained modeling of the temperature of radial and vertical profile of these disks, thus providing better evidence for the role the iceline plays on the substructures detected.

Fig. 5 Left panels: temperature profile of the midplane is obtained from Eq. (5) and plotted with a black curve. The blue curve shows the best radial profile obtained from the continuum UV-modeling, which was normalized to the peak emission. Vertical red and light blue shaded regions show the possible locations of the CO (light red) and N2 (light blue) icelines. Right panel: disk radius is shown versus the square root of stellar luminosity ($\text{}\sqrt{L_{\star }}$). The shaded regions show the iceline location for CO (light red) and N2 (light blue), obtained from Eq. (6). The peak location of the modeled continuum emission rings are shown in squares, while the gap locations of J0433 are shown with dots. The lower line in each target marks half of the theoretical resolution given by the longest baseline for each dataset, as measured from the center, and the upper line marks the continuum R90.

4.2.2 CIDA 1

Our modeling of the cavity in CIDA 1 gives similar results to those found in Pinilla et al. (2018a), thus we do not further explore new possibilities for its cavity origin. It is important to note that our finding of the inner side of the ring being wider than the outer side is opposite to what we expected from the physical motivation of using this model, which accounts for the timescale of grain growth from micrometer to millimeter particles in dust traps (Pinilla et al. 2017b, 2018a). One possibility is that there is unresolved emission inside the main ring, and the inner side of the Gaussian is blending with it. Higher angular resolution observations are needed to understand the true nature of the dust distribution of this system.

4.2.3 J0433

If the gaps are assumed to be generated by a planet-disk interaction, we can use the width of the gap to estimate the mass of this gap-carving planet, under the assumption that this single planet is located at the position of the gap, and the ring is peaked at the local pressure maxima of the gas. In this scenario, if the physical parameters of the disk are kept constant, a more massive planet would create a wider gap (e.g., Fung et al. 2014; Kanagawa et al. 2015; Rosotti et al. 2016).

For a crude estimate for the mass of the planet in the first J0433 gap, we followed a similar procedure as Long et al. (2018). In this approach, we assumed that the distance between the gap and ring scales with the Hill radius of the planet: $$R_{\text{Hill}}=r_{\text{p}}{(M_{\text{p}}/(3M_{\star }))}^{1/3}\text{,}$$(7)

where r_p is the location of the planet (the location of the gap coincides with the location of the planet, as the planet is carving the gap), M_p is the mass of the planet, and M_⋆ is the mass of the host star. If we consider the distance between the gap and ring to be 5 R_Hill (conservative upper limit for the gap width carved by a planet in Dodson-Robinson & Salyk 2011), and considering that our best model gives r_ring − r_gap = 5 au for the first ring, then the approximate mass of the planet would be M_p ~ 0.1 M_Jup.

This calculation has a large uncertainty depending on the type of simulation. For instance, Dodson-Robinson & Salyk (2011) estimated a maximum of4 R_Hill between theplanet location and the ring position, while Pinilla et al. (2012) estimated 7 R_Hill for planets as massive as Jupiter. Taking these two limits (4 R_Hill and 7 R_Hill), the estimated planet mass becomes $M_{\text{p}}={0.11}_{-0.07}^{+0.10}{M}_{\text{Jup}}$. In this calculation, we did not account for the physical conditions of the disk, such as turbulence or temperature, nor did we consider theminimum mass to open a gap in the disk. Therefore, it just gives us a crude estimate for the mass order of magnitude of a single planet carving the gap. A future analysis should consider higher angular resolution observations to better constrain the gap-ring morphology as well as dedicated hydro-simulations to estimate the planet candidate mass.

4.2.4 MHO 6

In MHO 6, the modeling indicates the existence of a central cavity that is slightly asymmetric. Several processes in protoplanetary disks could create this type of an inner cavity, such as photoevaporation due to stellar irradiation (e.g., Alexander & Armitage 2007; Owen et al. 2012; Owen & Kollmeier 2019), companions (e.g., Price et al. 2018), a planet-disk interaction (e.g., Rice et al. 2006; Zhu et al. 2011), and dead zones (e.g., Flock et al. 2015). In the following, we discuss each one of these scenarios for the MHO 6 cavity.

Photoevaporation

The 10 au dust cavity and the low accretion rate (5 × 10⁻¹¹ M_⊙ yr⁻¹ estimated from a UV excess in Herczeg & Hillenbrand 2008) are in agreement with the predicted evolution for a very low mass star (0.1 M_⊙) having its material stripped away by a photoevaporative flow (Owen et al. 2012). However, we do not see clear evidence of a gas depleted cavity in the ¹² CO emission, although the ¹³CO emission suggests a decrease in the gas density in the inner region. A follow-up with dedicated VLMS photoevaporation models and higher angular resolution of this target in different molecular lines would be able to characterize the impact of this mechanism in the cavity we detect. From the current observations, this mechanism is still an option for the observed cavity.

Companion

Binary companions produce cavities in circumbinary disks, with sizes in the range of 3− 5 times the binary semi-major axis, depending on the mass ratio, the eccentricity, and the disk viscosity (Artymowicz & Lubow 1994; Miranda et al. 2017; Ragusa et al. 2017). They also produce quasi-periodic variations in the accretion rate (Muñoz & Lai 2016). If MHO 6 was a circumbinary disk, the location of this companion could not be farther out than ~ 3 au (~ 20 mas) from the primary star. However, previous independent observations of MHO 6 have not found any evidence of multiplicity in this system (Briceño et al. 1998; Kraus & Hillenbrand 2007; Herczeg & Hillenbrand 2014, the latest reference identified companions in its survey as late-type as M0.0). A spectroscopic follow-up with radial velocities, given the high inclination of the system, could be useful to constrain the upper mass limit of a companion in the close inner region.

Embedded planets

The planet population around very low mass stars spans a wide range of masses, where even giant planets have been detected (e.g., Morales et al. 2019, a 0.46 M_Jup around a 0.1 M_⊙ VLMS), meaning that protoplanetary disks around VLMS probably have the potential to create such objects. If we make a simple calculation following the opening gap criterium introduced by Crida et al. (2006), the approximate mass needed for a single planet to open a gap and explain the cavity observed in MHO 6 is 1 or 0.6 M_Saturn (0.18 and 0.3 M_Jup), located at ~7 au for α = 10⁻³ and α = 10⁻⁴, respectively.If we assume a gas to dust ratio of 100, we estimate a disk mass of ~ 12.6 M_Saturn or ≈ 3.8 M_Jup given the dust mass obtained in Sect. 3.3. This means that the cavity opening planet is < 10% of the current estimated mass of the disk. This implies that this type of potential single planet could have formed within the disk of MHO 6, although we do not exclude the possibility of multiple planets being responsible of this central cavity.

Anothercharacteristic to take into consideration is the location of the peak brightness in the azimuthally averaged radial profile of ¹³CO and dust continuum emission. We find hints of the cavity in the ¹³CO line emission with its peak at 10 au when measured from the image. The same peak appears at 14 au when measured from the dust continuum image. Although these comparisons are strongly biased by the beam shape, it supports the idea that a planet may be responsible for the cavity formation since models predict such segregation between the outer edge of the gap in gas versus dust (de Juan Ovelar et al. 2013; Facchini et al. 2018). Future higher angular resolution observations could test this hypothesis by better characterizing the CO brightness profile.

Dead zone

A dead zone is a low-ionized region at the disk midplane, where the dense environment blocks the high energy radiation, suppressing the magneto-rotation instability and, therefore, the angular momentum transport. It has been shown that the presence of dead zones can open gaps and cavities by forming gas pressure bumps at the outer edge of the dead zones, where dust trapping is efficient (e.g., Regály et al. 2012; Flock et al. 2015). Pinilla et al. (2016) predict that cavities that formed by dead zones alone would have millimeter- and micrometer-sized particles concentrated at the peak of the gas density. As a result, the radial location of the ring in millimeter wavelengths and scatter light would be the same. If we neglect temperature effects, our images suggest that the peak of ¹³CO is closer to the star than the millimeter peak. Therefore, we expect for the micrometer-sized particles to peak closer to the star as well. In the dead zone scenario, Pinilla et al. (2016) also predict a strong gas depletion in the outer parts of the disk, further out from the dust trap, which is not currently seen in our observations of ¹³CO. If a magnetohydrodynamical wind is included in the models together with a dead zone, the observational diagnostics are very similar to planets. One step forward to try to disentangle between these models is to image this disk in scattered light and search for potential planets in the cavity. However, these disks are too faint to image given the limitations of current telescopes in the optical and near-infrared.

To summarize, the formation of the cavity observed in MHO 6 could be explained by one or a combination of the mechanisms we have discussed above. Several observational efforts can be done to disentangle these possibilities, such as a better estimate of the star accretion rate, a deeper search for planets or companions in its cavity in the optical and infrared wavelengths as well as imaging the disk in this regime, and very deep and higher angular resolution observations of molecular lines.

4.3 MHO 6 kinematics

Given that MHO 6 is the only disk with the gas emission detected across the whole velocity range, it was a good candidate for the dynamical mass measurement of its star. After several attempts to model MHO 6 rotation with different images, masks, and a combination of free parameters, we find our range of recovered stellar masses (between ~ 0.16 and 0.24 M_⊙) to be consistent with the values from pre-main sequence models of MHO 6 (0.09−0.20, M_⊙ Kraus & Hillenbrand 2009; Herczeg & Hillenbrand 2014; Ward-Duong et al. 2018), but all the kinematic models would leave strong structured residuals with an amplitude spanning two times the channel velocity width.

Although our observations have enough spatial resolution to resolve the vertical structure in some of our disks (as seen in Fig. 1 as a cone-like emission in the ¹² CO moment 0 of MHO 6 and J0433), the low S/N did not allow us to differentiate between the emission coming from the back side of the disk from that of the front side. The integrated velocity map contains the emission of both sides as if they were the same, and thus parameters such as M_⋆ or ψ from Z(r, ψ) ∝ r^ψ could not be recovered reliably.

This issue is not an exclusive problem in VLMS disks, but it applies to all the gas measurements with a high angular resolution and poor S/N. Future approaches to accurately recover M_⋆ should consider more robust methods, such as UV-modeling (similar to the DiskJockey code from Czekala et al. 2015).

MHO 6 is a good candidate to further study the kinematics, substructures, and the physics of planet formation in VLMS. Its brightness allows high S/N observations at high angular resolution with non-prohibitive integration times. Combining the datasets presented in this paper, plus an observation of ≥ 5 h of time on source with ALMA using long baseline configurations should have enough resolution and sensitivity to precisely characterize the kinematics of the CO isotopologues, as previously done with disks around T-Tauri and Herbig stars (e.g., Pinte et al. 2019; Teague et al. 2019).

4.4 Comparison of dust R₆₈ and L_mm with previous studies

All our sources have observed or estimated 0.87 mm fluxes, and a subset also had their dust continuum R₆₈ measured with previous UV-modeling. In the following, we compare and discuss our measurements with previous results, and we include the VLMS in the study of a size-luminosity relation.

Flux

J0415 was observed using the Submillimeter Array (SMA) at a wavelength of 1.3 mm, with a reported flux density of 12.6 ± 1.4 mJy in Andrews et al. (2013), where they extrapolated this measurement to 0.87 mm using F_ν ∝ ν^α, with α = 2.4 ± 0.5, thus obtaining F_0.87mm = 32.9 ± 15.2 mJy. Although for all our other five sources, the extrapolations from the SMA flux measurements are consistent within the error range, the flux at 0.87 mm received from J0415 is approximately 35 times dimmer than expected, as observed independently by two different ALMA projects. A possible explanation could be a sudden change in millimeter flux due to flares, as it has been observed in low mass stars (e.g., MacGregor et al. 2018). For now, the only ALMA observations of this disk are the ones presented in this work, which are both in Band 7. Future observations at 1.3 mm should solve this discrepancy with the SMA results.

Dust R₆₈

The radius enclosing 68% of the dust continuum emission (R₆₈) was previously estimated for MHO 6 and J0433, using UV-modeling on SMA 340 GHz observations. In Tripathi et al. (2017) and Andrews et al. (2018b), they estimated a value of $R_{68}={36.9}_{-5.7}^{+8.5}$au for MHO 6 ($0.2{6^{″}}_{-0.04}^{+0.06}$), and $R_{68}={58.9}_{-8.7}^{+6.9}$au for J0433 ($0.3{4^{″}}_{-0.05}^{+0.04}$), which are consistent with the values measured by this work. Both are slightly overestimated (1.6 times for J0433) compared to the values of this work due to the considerably lower angular resolution of SMA observations compared to ALMA (0.86″ × 0.80″ for MHO 6 and 0.61″ × 0.52″ for J0433).

R₆₈ – L_mm relation

We combined our VLMS measurements with the Taurus sample observed with the SMA at 340GHz from Tripathi et al. (2017) in order to compare this all with the recent analysis by Hendler et al. (2020), where they obtained a relation of $R_{68}\propto L_{\text{mm}}{}^{0.53}$ for the Taurus star-forming region. For comparison purposes with this study, we followed the same approach using the Bayesian linear regression described in Kelly (2007), which was implemented in the python package linmix (publicly available in github, see Meyers 2015), to fit a linear relation between R₆₈ and L_mm following: $${\log }_{10}(R_{68})=\alpha +\beta {\log }_{10}(L_{\text{mm}})\text{,}$$(8)

where α and β are the regression coefficients. Our best fit was calculated by using the median value of the last 200 000 steps after convergence. We find that including our VLMS sample does not statistically change the previous result. However, the inclusion of J0415 changes the steepness of this relation in about 1σ. If we consider J0415 as part of the fitting data, then the relation recovered is α = 2.09 ± 0.09, β = 0.47 ± 0.08. If we exclude J0415 from the fitting, the result is α = 2.19 ± 0.10, β = 0.59 ± 0.10. Both relations are close to the 1 σ limit of each other, and they also overlap with the previous calculation from Hendler et al. (2020), as shown in Fig. 6. If we check their relations’ intrinsic scatter from the linear regression, we obtain σ_scatter = 0.231 ± 0.041 and σ_scatter = 0.217 ± 0.040, respectively. Therefore, including J0415 does not significantly increase the scatter.

The difference in the steepness of the recovered relations is not statistically significant from what was previously found. However, it is not completely clear if the same single power law relation between the millimeter luminosity and the size of the disks holds along the whole luminosity range. Given that J0415 is the only source with its size measured in the ~ 1 mJy brightness range, it is unknown if disks have some mechanisms to remain extended even when they are low in dust content (thus flattening the relation between size and luminosity in the low luminosity regime), or if a J0415-extended dust size is part of the relation scatter that is also observed in bright disks. To understand if J0415 is an outlier and to test if the power law behavior of the R₆₈ − L_mm relation flattens or holds at the low brightness regime, we need more deep observations at a high angular resolution of disks with F_0.87mm < 10 mJy. This could be achieved by observing each source from several tens of minutes to a few hours in ALMA Band 7.

Fig. 6 Relation between log(R68) and log (F0.87mm). The points in blue are from Tripathi et al. (2017), with distances scaled to 140 pc using distances inferred from Gaia DR2, similar to Andrews et al. (2018b). The solid line in blue is the linear regression found by Hendler et al. (2020) (abbreviated as H+20), fitting the data from Tripathi et al. (2017). The VLMS from this work are shown in red. The bests linear regression fits and their 68% confidence intervals are plotted in red when J0415 is excluded from the fit, and in green when all points are considered.

4.5 J0415 dust radial extent

Although CIDA 7 has a dust content of at least ~27 times higher than J0415, it remains an open question as to how both can have a similar size (see Table 3). This result could be due to weaker dust traps in J0415 compared to those in CIDA 7, thus CIDA 7 can trap the dust more efficiently.

Alternatively, if the M_gas of the disk is very low, the millimeter grains would have a high Stokes number and radial drift would become negligible. In Pinilla et al. (2017a), they explore this scenario with a disk of 60 au in radius around a BD of 0.05 M_⊙. When M_gas = 2 × 10⁻² M_Jup, diffusion and drift still depleted the disk from millimeter particles; however, when M_gas = 2 × 10⁻³ M_Jup, the millimeter grains were decoupled from the gas, and they could remain in the disk for longer timescales. For J0415, if we assume a dust to gas ratio of 1/100, we obtain M_gas ≈ 4 × 10⁻² M_Jup, which does not seem to be low enough for a complete dust-gas decouple to occur. Additional observations are needed in order to characterize its gas radial density profile and to test this possibility.

Finally, it has been proposed that dust grains could be growing in a fractal manner, such that large aggregates would avoid radial drift by maintaining a low Stokes number (Kataoka et al. 2013). This scenario can be distinguished from the compact millimeter grains by measuring the opacity index β between the 1 and 3 mm emission (Kataoka et al. 2014). Those observations, however, would require a very high sensitivity with enough angular resolution to spatially resolve the radial profile of J0415.

4.6 R_gas/R_dust ratio

We find that four of our six disks show a ratio between the R_90,gas and R_90,dust that is very close to or above 4.0, which is similar to the ratio measured in CX Tau (Facchini et al. 2019) (ratio of 3.9 ± 0.5 at R₉₀), as compared in Fig. 7. Our values, however, are conservative measurements of radii ratios for the sources with strong cloud contamination (all the sources in our sample, except for MHO 6), and so they should be considered as a lower limit. As it can be seen in the velocity channel maps in Figs. B.5, B.7, B.8, B.9, and B.10, the channels that are less affected by cloud contamination are the highest velocity channels, which are closer to the star. As soon as the emission gets radially extended in the channels closer to the central velocity, the cloud extinction becomes so high that we cannot see the disk emission anymore. Our integrated flux moment 0 images and the gas emission are then biased towards the compact emission located close to the star, underestimating the R_gas measurement. This effect is particularly strong in J0420 and J0415, as they are the dimmest sources in this work, and J0420 also has the highest extinction (see Fig. A.1).

Our observations show that it is common for bright disks around VLMS to have a gas radial extension of > 3 times the millimeter dust radial extension, which is consistent with the efficient radial drift expected for millimeter-sized grains in VLMS disks (Pinilla et al. 2013; Zhu et al. 2018). In fact, recent thermochemical modeling including dust evolution by Trapman et al. (2019, 2020) shows that ratios of R_gas∕R_dust > 4 cannot be solely explained by effects of optical depths; additionally, in those cases, radial drift is required to explain the gas and dust size difference. Despite the limited sensitivity and cloud contamination, our disks are still very close and even above that limit. Although similar ratios have been observed in sources with a moderate stellar mass (~ 0.5 M_⊙, see Ansdell et al. 2018), as shown in the upper panel of Fig. 7, the direct comparison between our works ishindered by the data differences, such as in the sensitivity and angular resolution, and also by the different approaches used to obtain the gas and dust radii: our continuum radial extension was calculated from the visibilities rather than the images, and we did not apply Keplerian masking in the flux integrated gas images. A more complete sample, which includes more observations of disks around very low and moderate mass stars, and also uniform data sensitivity and analysis are required to confirm if there is a general trend where R_gas∕R_dust increases towards lower stellar masses.

According to Trapman et al. (2019, 2020), they expect more massive disks to have a higher R_gas∕R_dust ratio, driven by a larger observed R_gas due to the greater total CO content, and also because the higher dust content would produce a more efficient grain growth and inward radial drift. However, we do not observe this trend in the lower panel of Fig. 7. Apparently, a decreasing radii ratio is obtained towards higher dust mass disks, which could be the result of an efficient radial drift in disks around VLMS, and the linear relation between log (M_⋆) − logM_dust. This supports the idea that smaller disks result from a fast radial drift (Long et al. 2019) due to their inability to trap dust in the outer regions, while more massive disks are more capable of creating dust traps father away from the central star.

CIDA 7stands out in our sample as having the most extreme R_gas∕R_dust ratio observed so far, with the gas being six times more extended than the dust, which is well beyond the ratio of four limit from Trapman et al. (2019). This confirms that radial drift is responsible for the compact size of this source. However, it is not completely dust depleted, so radial drift has been counteracted by another mechanism, or a combination of them. The southern non-Keplerian emission detected in this system (in the velocity channels 4.4 and 4.8 km s⁻¹ from Figs. B.8 and 1) was masked when measuring the R_90,gas, so the ratio of six is between the disk rotating gas and the dust size recovered from the model. We were unable to determine the origin of this extended emission in the south, as it is only detected in two different channels and it does not appears to be axisymmetric. A multiwavelength follow-up, with high sensitivity and angular resolution, might be required to understand the nature of this emission, as it could be explained by several different mechanisms, such as winds, outflows, interactions with external companions, an interaction with the surrounding cloud and envelope, among others.

The lowest gas-to-dust size ratio in our sample, measured in J0415, is likely due to the combined effects of lower than expected brightness, a low sensitivity due to our inability to apply self-calibration, and the extinction due to cloud contamination. Although we were unable to confidently recover the gas radius, our observations set a lower limit for its radial extension, and we confirm that in this very low disk mass regime, the gas emission is still more extended than dust emission. A more precise measurement of the gas radius requires a combination of deeper observations of lines less affected by the surrounding cloud and envelope as well as line modeling.

Fig. 7 Upper panel: Rgas∕Rdust as a function of the stellar mass of the disks in the Lupus star-forming region (SFR) reported in Ansdell et al. (2018) (abbreviated as A+18), CX Tau (Facchini et al. 2019), and the targets reported in this paper. For CX Tau, we took the R90 radii of the dust and gas. The VLMS of this work are shown in red for the sources with a substructure, and they are in yellow for the smooth sources, with the values from Table 3. Ratios from the Lupus SFR by Ansdell et al. (2018) were calculated from R90 and a different color was used for disks identified as transition disks (TD). Lower panel: Rgas ∕Rdust of the same targets, but compared to the dust mass of the disks. The dust mass was calculated from Eq. (4) under the assumption that T = 20 K for the whole disk midplane, with distances from Gaia DR2.

5 Conclusions

To understand the process of planet formation in VLMS and compare it with our current knowledge of planet formation in solar type stars, we observed and studied a sample of the brightest disks around VLMS in Taurus, at a 0.1″ resolution and at a 0.87 mm wavelength. This sample is composed of CIDA 1, MHO 6, J0433, CIDA 7, J0420, and J0415 (2MASS names in Sect. 2). Here we summarize our main conclusions as follows.

• Detection rate of substructures: millimeter dust substructures were directly detected in only 50% of the targets in our sample. Our results suggest that the detection of substructures in disks around VLMS is limited by angular resolution and sensitivity, since the dust radial extent is very small and these disks are also very faint. Deep, high angular resolution observations over a non-brightness biased sample of VLMS should confirm the ubiquity of substructures in these disks.

• Substructured disks: substructures were detected in CIDA 1, MHO 6, and J0433; with the latest two being new detections. These three disks are the brightest and largest in our sample. They all have axisymmetric ring-like substructures, and only MHO 6 shows a weak asymmetry of amplitude less than 5% of the peak brightness. Both CIDA 1 and MHO 6 show central cavities in their emission. If we assume that a planet-disk interaction is the origin of the MHO 6 cavity, then a Saturn-mass planet (0.3 M_Jup) is needed (as in the case of CIDA 1 Pinilla et al. 2018a). This planet should be located around 7 au. However, we cannot exclude other mechanisms that can also explain the origin of this cavity, such as multiple planets, a dead zone, a binary companion, or photoevaporation. Our UV-modeling of J0433 suggests that this disk could have two rings located at 21 and 32 au. However, we cannot confirm the separation between both at the 1 σ errorbar.We estimate that a planet of ~0.1 M_Jup in mass could explain the first gap-ring. The substructures were detected within the region where the CO iceline could be located. Our current datasets lack the necessary S/N and resolution to properly characterize the vertical and radial temperature profile of the CO isotopologues, and so future deeper observations will be needed to determine if the iceline played any role in triggering or maintaining the substructures observed.

• Smooth disks: the dust disks in CIDA 7, J0420, and J0415 are the less radially extended, less massive disks of our sample. With an angular resolution of 0.1″, these disks are well described by a single Gaussian radial profile, which we used to measure their sizes. CIDA 7 is the most compact of them, with an R₉₀ = 13.16 au, which is similar to the 15.46 au from J0415. Yet, the dust mass estimate suggests that CIDA 7 is ~27 times more massive. In J0420, the residual continuum image shows some structured non-axisymmetric emission with 5σ peaks. However, this emission is very low in contrast to the smooth emission, which is over 300σ at its peak. Higher angular resolution observations are needed to describe the potential substructures in these disks.

• Size-luminosity relation: the disks in our sample follow a similar relation between L_mm − R₆₈ as the one found for bright disks in the same star-forming region (see Hendler et al. 2020). However, our single measurement of a disk size in the low luminosity regime (J0415) needs to be complemented with deeper additional observations of other sources with a low stellar mass and low disk brightness. These measurements will help us understand the behavior of the size-luminosity relation across the whole range of disk sizes, enabling us to test if a single power law describes it.

• Evidence of efficient radial drift: when considering the dust and gas radii as the location where 90% of the emission is enclosed, four out of six disks in our sample show a ratio between R_gas∕R_dust above 3.5. This is expected for disks where radial drift is depleting the dust, suggesting that radial drift is more efficient in VLMS than in moderate or high mass stars. However, the analysis of the sizes of disks in Lupus by Ansdell et al. (2018) also allowed them to find moderate stellar mass sources with similar disk size ratios as those observed in VLMS, which suggest that we need more observations to confirm this trend. Our comparison with Ansdell et al. (2018) does not account for the differences in data acquisition and R₉₀ calculation, which should be considered in future works, aiming for a uniform analysis of the extension of disks between different star-forming regions. The most extreme case of high R_gas∕R_dust in our VLMS is observed in CIDA 7, with a value of six. This very high R_gas∕R_dust ratio suggests that strong radial drift is at play (Trapman et al. 2019), raising the question about how this disk remains massive in dust.

Our observations do not exclude giant planet formation as an explanation for the substructures detected. Disks around VLMS follow similar trends as those that have been observed in disks around higher mass stars, based on our sample of bright disks. A future confirmation of a deviation from current correlations of physical parameters will require the recalculation of fluxes and sizes of the Taurus disks by using more sensitive and higher angular resolution observations from ALMA in a larger and more complete sample.

Acknowledgements

We thank Timmy Delage and Reema Joshi for their very helpful comments. We also thank Sean M. Andrews for the insightful discussions, which significantly improved this manuscript. N.K. and P.P. acknowledge support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. This Project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 823823 (DUSTBUSTERS). This work was partly supported by the Deutsche Forschungs-Gemeinschaft (DFG, German Research Foundation) - Ref no. FOR 2634/1 TE 1024/1-1. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00743.S, ADS/JAO.ALMA#2015.1.00934.S, ADS/JAO.ALMA#2016.1.01511.S, and ADS/JAO.ALMA#2018.1.00310.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Appendix A Summary of ALMA observations

A summary of the observation log is provided in Table A.1. The difference for configurations “Compact” and “Extended” denotes the stage at which we started its self-calibration based on the spatial extension of the antenna array, see Sect. 2.2. The properties of the dust continuum emission images and gas emission datacube are shown in Tables A.2 and A.3, respectively. In Fig. A.1 we show the reddening map from Schlafly et al. (2014) of the Taurus star-forming region, with the location of our sources. The reddening values for the sources are 0.61 for CIDA 1, 0.48 for MHO 6, 0.51 for J0433, 0.50 for CIDA 7, 1.38 for J0420, and 0.44 for J0415.

Fig. A.1 Spatial distribution of the Taurus sample studied in this work. The background is an extinction map compiled by Schlafly et al. (2014).

Appendix B Spectralprofiles, eddy fitting, and channel maps

The spectralprofiles of ¹²CO and ¹³CO for each source, shown in Fig. B.1, were obtained with the CASA 5.6.2 software by placing an elliptical mask in the region where emission was detected. The errorbars represent the standard deviation of the noise, measured from the channels without emission.

Fig. B.1 Spectral profile of the 12CO and 13CO in our sample. The gray line denotes the 0 Jy level.

The moment images shown in Fig. B.2 were also computed with CASA 5.6.2 by clipping the emission in each channel at the 3 σ level. For the Keplerian fitting with eddy (Teague 2019), we fixed the inclination to be one recovered from the continuum UV-modeling. The fitting was performed inside an elliptical masked region, with the same inclination of the source, extending up to 0.9″ in the ¹² CO emission, and up to 0.6″ in the ¹³CO. Each image was modeled separately.

The models and residuals shown in Fig. B.2 allowed the center of the source (x₀, y₀), the position angle (PA), the mass of the star (M_⋆), the central velocity (VLSR), and the vertical structure (Z(r) = z₀ ⋅ r^ϕ) to vary. Thecentral velocity was consistent between the two lines, and PA was also consistent with continuum. For the stellar mass, ¹² CO would give a value of M_⋆ = 0.195 ± 0.02 M_⊙, and M_⋆ = 0.19 ± 0.1 M_⊙ from ¹³CO. Fixing someof these parameters to simplify the model, or using images obtained with bettermoments, would change M_⋆ between 0.16~0.24 M_⊙, but the structure of residuals would remain the same.

Fig. B.2 Keplerian fitting of 12CO and 13CO. The grid shows the best surface recovered, extending up to the distance of the mask for the fitting.

Fig. B.3 Position velocity diagrams for the sources where substructures are detected in the continuum emission. Lines follow the Keplerian velocities for central mass objects of 0.1, 0.2, and 0.3 M⊙ (solid, dashed, and dotted, respectively).

Fig. B.4 Same as Fig. B.3, but for the sources without a substructure.

Fig. B.5 Images of CIDA 1 12CO channel maps are shown in the top, and 13CO channel maps are shown in the bottom. Each square is 3.0″ in width andheight, centered at the samed 13CO channel maps are shown in the bottom position as the dust continuum image. The contour level traces the 5σ emissionin the continuum image. The scale bar in the lower left panel is 20 au in size, and the beam size is found in the lower left corner of the same panel. Central velocities in km s−1 of each channel are given in the upper left corner of each panel.

Fig. B.6 Same as Fig. B.5, but for MHO 6.

Fig. B.7 Same as Fig. B.5, but for J0433.

Fig. B.8 Same as Fig. B.5, but for CIDA 7.

Fig. B.9 Same as Fig. B.5, but for J0420.

Fig. B.10 Same as Fig. B.5, but for J0415.

References

All Tables

All Figures

Fig. 1 ALMA observations of the disks around VLMS in Taurus. From left to right: dust continuum emission, 12 CO moments 0 and 1, and 13CO moments 0 and 1. All boxes are 3.0″ in size. The scale bar represents 20 au, and the white ellipses show the synthesized beam. The centers, beam sizes, measured fluxes, and sensitivity are detailed in Tables A.2 and A.3. Dashed lines show the region used to calculate the radial profiles and to measure R68 and R90. A zoomed-in version of the continuum images can be found in Fig. 2.

In the text

Fig. 2 Visibility modeling versus observations of our sample. From left to right: (1) real part of the visibilities after centering and deprojecting the data versus the best fit model of the continuum data, (2) continuum emission of our sources where the scale bar represents a 10 au distance, (3) model image, (4) residual map (observations minus model), (5) and normalized, azimuthally-averaged radial profile calculated from the beam convolved images in comparison with the model without convolution (purple solid line) and after convolution (red solid line). In the right most plots, the gray scale shows the beam major axis.

In the text

Fig. 3 Normalized azimuthally averaged radial profiles of the dust continuum, 12 CO, and 13CO of each source. Each curve has a 68% error region. The orange and red solid curves correspond to the average intensity profile for the gas profiles of the deprojected images, while the blue solid line corresponds to the best χ2 solution from the visibility fit of the dust continuum. The dashed vertical line denotes the position and 1σ error of the R90 radius for the dust (blue) and the gas (red).

In the text

Fig. 4 Dust mass compared to stellar mass for the Taurus disks. We note that Md was computed using distances inferred from Gaia DR2, assuming a 20K disk temperature and using the most updated value of F0.87mm from Andrews et al. (2013), Tripathi et al. (2017), and Ward-Duong et al. (2018). Sources where F0.87mm was extrapolated from a F1.3mm measurement are plotted in faded colors (see Table 2 in Andrews et al. 2013), while direct measurements from SMA or ALMA are plotted in darker colors. Green triangles denotes upper limits, and the VLMS studied in this work are plotted as orange stars.

In the text

Fig. 5 Left panels: temperature profile of the midplane is obtained from Eq. (5) and plotted with a black curve. The blue curve shows the best radial profile obtained from the continuum UV-modeling, which was normalized to the peak emission. Vertical red and light blue shaded regions show the possible locations of the CO (light red) and N2 (light blue) icelines. Right panel: disk radius is shown versus the square root of stellar luminosity ($\text{}\sqrt{L_{\star }}$). The shaded regions show the iceline location for CO (light red) and N2 (light blue), obtained from Eq. (6). The peak location of the modeled continuum emission rings are shown in squares, while the gap locations of J0433 are shown with dots. The lower line in each target marks half of the theoretical resolution given by the longest baseline for each dataset, as measured from the center, and the upper line marks the continuum R90.

In the text

Fig. 6 Relation between log(R68) and log (F0.87mm). The points in blue are from Tripathi et al. (2017), with distances scaled to 140 pc using distances inferred from Gaia DR2, similar to Andrews et al. (2018b). The solid line in blue is the linear regression found by Hendler et al. (2020) (abbreviated as H+20), fitting the data from Tripathi et al. (2017). The VLMS from this work are shown in red. The bests linear regression fits and their 68% confidence intervals are plotted in red when J0415 is excluded from the fit, and in green when all points are considered.

In the text

Fig. 7 Upper panel: Rgas∕Rdust as a function of the stellar mass of the disks in the Lupus star-forming region (SFR) reported in Ansdell et al. (2018) (abbreviated as A+18), CX Tau (Facchini et al. 2019), and the targets reported in this paper. For CX Tau, we took the R90 radii of the dust and gas. The VLMS of this work are shown in red for the sources with a substructure, and they are in yellow for the smooth sources, with the values from Table 3. Ratios from the Lupus SFR by Ansdell et al. (2018) were calculated from R90 and a different color was used for disks identified as transition disks (TD). Lower panel: Rgas ∕Rdust of the same targets, but compared to the dust mass of the disks. The dust mass was calculated from Eq. (4) under the assumption that T = 20 K for the whole disk midplane, with distances from Gaia DR2.

In the text

	Fig. A.1 Spatial distribution of the Taurus sample studied in this work. The background is an extinction map compiled by Schlafly et al. (2014).
In the text

	Fig. B.1 Spectral profile of the 12CO and 13CO in our sample. The gray line denotes the 0 Jy level.
In the text

	Fig. B.2 Keplerian fitting of 12CO and 13CO. The grid shows the best surface recovered, extending up to the distance of the mask for the fitting.
In the text

	Fig. B.3 Position velocity diagrams for the sources where substructures are detected in the continuum emission. Lines follow the Keplerian velocities for central mass objects of 0.1, 0.2, and 0.3 M⊙ (solid, dashed, and dotted, respectively).
In the text

	Fig. B.4 Same as Fig. B.3, but for the sources without a substructure.
In the text

Fig. B.5 Images of CIDA 1 12CO channel maps are shown in the top, and 13CO channel maps are shown in the bottom. Each square is 3.0″ in width andheight, centered at the samed 13CO channel maps are shown in the bottom position as the dust continuum image. The contour level traces the 5σ emissionin the continuum image. The scale bar in the lower left panel is 20 au in size, and the beam size is found in the lower left corner of the same panel. Central velocities in km s−1 of each channel are given in the upper left corner of each panel.

In the text

	Fig. B.6 Same as Fig. B.5, but for MHO 6.
In the text

	Fig. B.7 Same as Fig. B.5, but for J0433.
In the text

	Fig. B.8 Same as Fig. B.5, but for CIDA 7.
In the text

	Fig. B.9 Same as Fig. B.5, but for J0420.
In the text

	Fig. B.10 Same as Fig. B.5, but for J0415.
In the text