Size and structures of disks around very low mass stars in the Taurus star-forming region

We aim to estimate if structures, such as cavities, rings, and gaps, are common in disks around VLMS and to test models of structure formation in these disks. We also aim to compare the radial extent of the gas and dust emission in disks around VLMS, which can give us insight about radial drift. We studied six disks around VLMS in the Taurus star-forming region using ALMA Band 7 ($\sim 340\,$GHz) at a resolution of $\sim0.1''$. The targets were selected because of their high disk dust content in their stellar mass regime. Our observations resolve the disk dust continuum in all disks. In addition, we detect the $^{12}$CO ($J=3-2$) emission line in all targets and $^{13}$CO ($J=3-2$) in five of the six sources. The angular resolution allows the detection of dust substructures in three out of the six disks, which we studied by using UV-modeling. Central cavities are observed in the disks around stars MHO\,6 (M5.0) and CIDA\,1 (M4.5), while we have a tentative detection of a multi-ringed disk around J0433. Single planets of masses $0.1\sim0.4\,M_{\rm{Jup}}$ would be required. The other three disks with no observed structures are the most compact and faintest in our sample. The emission of $^{12}$CO and $^{13}$CO is more extended than the dust continuum emission in all disks of our sample. When using the $^{12}$CO emission to determine the gas disk extension $R_{\rm{gas}}$, the ratio of $R_{\rm{gas}}/R_{\rm{dust}}$ in our sample varies from 2.3 to 6.0, which is consistent with models of radial drift being very efficient around VLMS in the absence of substructures. Our observations do not exclude giant planet formation on the substructures observed. A comparison of the size and luminosity of VLMS disks with their counterparts around higher mass stars shows that they follow a similar relation.


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 ( 0.1 M , VLMS, 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. 2006Scholz et al. , 2007Ricci et al. 2012;Daemgen et al. 2016;van der Plas et al. 2016;Ricci et al. 2017a,b;Sanchis et al. 2020), which are more compact and lower in dust mass when compared to disks around T-Tauri stars Article number, page 1 of 26 arXiv:2012.02225v1 [astro-ph.EP] 3 Dec 2020 (e.g., Pinilla et al. 2017c;Ward-Duong et al. 2018;Hendler et al. 2017Hendler et al. , 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. 2017c), 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. 2018b). 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. 2018Long et al. , 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. (2018b) have thus far reported the lowest mass star with a resolved large dust cavity (radius ∼20 au) in the disk around the M4.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. 2017c;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 12 CO and 13 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: Section 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 12 CO and 13 CO emission. In Section 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.

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 observations by Andrews et al. (2013) and Ward-Duong et al. (2018), respectively); (2) it has a stellar mass between ∼ 0.1 − 0.2 M ; 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 , which shows that those disks usually have higher dust masses compared to others with stellar hosts of similar mass. From the list of targets that fulfilled the conditions, we selected the ones with the lowest optical extinction.

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 et al. 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 Figure A Table A.1 in the appendix. 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 −1 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 selfcalibration 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 12 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 1 .
The final images of the dust continuum, moment 0, and moment 1 of the 12 CO and 13 CO (when detected) are shown in Figure 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 Table A.2 and A.3, respectively.

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:  (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, azimuthallyaveraged 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.
Article number, page 5 of 26 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. 2015Pinilla et al. , 2017a. 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. 2018;Huang et al. 2020), which is the same model used in Pinilla et al. (2018b) to model CIDA 1. The intensity profile is given by a ring as follows: where I bg is the broken Gaussian intensity profile as a function of the radius, r 1 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 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 position angle (PA). Therefore, each model has four more free parameters in addition to those that describe the intensity profile. The Fourier transform and the χ 2 calculation were carried out with the galario package . The pixel size used in the models is 1 mas, which is several times smaller than the smallest resolvable scale of the observations. The χ 2 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 250000 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 without convolution (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 68 and R 90 , respectively). We used the different sets of parameters sampled by the walkers to compute their radial profiles, and we found the 16th and 84th percentile on the R 68 and R 90 to calculate the error bars. The results are shown in Table 3 and in Figure 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.

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. (2018b). 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 Figure  3), we cannot exclude that it is an extension of the outer side of the first ring.

Dust disk masses
Assuming the dust continuum emission is optically thin, we estimated the disk dust mass by following Hildebrand (1983): 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 κ ν = 2.3 cm 2 g −1 (ν/230 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. 2018b). We estimated the uncertainty of the dust mass by Article number, page 6 of 26 Table 2. Best parameters from UV-modeling, following equations (1), (2), and (3). "mas" stands for milliarcsecond. The resulting F0.87mm of each model is given in the last row (the measured F0.87mm from the data is in  Normalized azimuthally averaged radial profiles of the dust continuum, 12 CO, and 13 CO 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).
taking the 10% uncertainty in the flux calibration. Our results for each source are compiled in Table A

Lines detection in each source
Both molecular lines were detected in all of our sources, except for J0415, where only 12 CO was observed. As shown in Figures 1 and 2, our angular resolution of ∼ 90 mas was just enough to resolve the dust continuum and gas emission Article number, page 7 of 26  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.
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 12 CO emission in J0420 is completely absorbed. The rotational pattern was, however, recovered in the 13 CO (see Figures 1 and B.9). On the other hand, in CIDA 7, we observed an extended asymmetric emission in the south region, with a velocity range of at least 0.8 km s −1 (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 13 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 Figure 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 Sections 3.4.3 and 4.3.

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 12 CO and 13 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 Figure 1). The gas R 68 and R 90 radii that encloses 68% and 90% of the flux are shown in Table 3  As we discuss in Section 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 12 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 σ .

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 13 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 Section 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 0 , y 0 ) 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 Figure 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 Figures B.3 and B.4 for the sources with continuum substructures and smooth profiles, respectively.

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 res-olution 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 68 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 Figure 2).

Origin of dust continuum rings and cavities
All the detected substructures in our sample resemble ringlike 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. 2018b;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.

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 Article number, page 9 of 26 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: 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 = 4πR 2 σ sb T 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 Equation (5), we can obtain the distance r i at which the temperature T i is reached, given a star of fixed luminosity L , following: 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 2 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 Figure 5.
In the right panel of Figure 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 CO iceline 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 van der Marel et al. 2018). Pinilla et al. (2017b) 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.

CIDA 1
Our modeling of the cavity in CIDA 1 gives similar results to those found in Pinilla et al. (2018b), 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. 2017a(Pinilla et al. , 2018b. 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.

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: 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 of 4 R Hill between the planet 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 p = 0.11 +0.10 −0.07 M Jup . In this calculation, we did not account for the physical conditions of the disk, such as turbulence or temperature, nor did we consider the minimum 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.

MHO 6
In MHO 6, the modeling indicates the existence of a central cavity that is slightly asymmetric. Several processes  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 −11 M yr −1 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 12 CO emission, although the 13 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. , located at ∼7 au for α = 10 −3 and α = 10 −4 , 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 Section 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.
Another characteristic to take into consideration is the location of the peak brightness in the azimuthally averaged radial profile of 13 CO and dust continuum emission. We find hints of the cavity in the 13 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 Article number, page 11 of 26 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 13 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 13 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.

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 − 0.24 M ) to be consistent with the values from pre-main sequence models of MHO 6 (0.09 − 0. 20M 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 Figure 1 as a cone-like emission in the 12 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 hr 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).

Comparison of dust R 68 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 68 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 68 : The radius enclosing 68% of the dust continuum emission (R 68 ) 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. (2018a), they estimated a value of R 68 = 36.9 +8.5 −5.7 au for MHO 6 (0.26 +0.06 −0.04 ), and R 68 = 58.9 +6.9 −8.7 au for J0433 (0.34 +0.04 −0.05 ), 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 68 − 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 ∝ L 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 68 and L mm following: log 10 (R 68 ) = α + β log 10 (L mm ), where α and β are the regression coefficients. Our best fit was calculated by using the median value of the last 200000 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 Figure 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 68 − 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.

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. (2017c), they explore this scenario with a disk of 60 au in radius around a BD of 0.05M . When M gas = 2 · 10 −2 M Jup , diffusion and drift still depleted the disk from millimeter particles; however, when M gas = 2·10 −3 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 −2 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 mm 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.

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 ) (ratio of 3.9 ± 0.5 at R 90 ), as compared in Figure 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 Figures 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. (2019Trapman et al. ( , 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 Figure 7, the direct comparison be-  , 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/R dust of the same targets, but compared to the dust mass of the disks. The dust mass was calculated from Equation (4) under the assumption that T = 20 K for the whole disk midplane, with distances from GAIA DR2.
tween our works is hindered 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) and Trapman et al. (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 Figure 7. Apparently, a de-creasing 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 ) − log M 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 7 stands 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 −1 from Fig. 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.

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 Section 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 Article number, page 14 of 26 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 planetdisk 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. 2018b). 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 90 = 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 68 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 ob-servations to confirm this trend. Our comparison with Ansdell et al. (2018) does not account for the differences in data acquisition and R 90 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.