Small and Large Dust Cavities in Disks around Mid-M Stars in Taurus

We analyze disks observed in the ALMA project 2019.1.00566.S (PI: G. Herczeg), which aimed at studying the dependence of dust substructure properties on stellar mass by targeting 16 disks around M3–M5 stars selected from Taurus. The sample selection started from Luhman et al. (2017) with disk identification based on excess emission in Wide-field Infrared Survey Explorer (WISE) W2 and W3 bands. We excluded known binaries (e.g., Kraus et al. 2011; Daemgen et al. 2015) and sources with A_V > 3 mag or J-band brightness 1 mag fainter than the median to avoid edge-on disks and embedded objects. Furthermore, targets with decl. > 26° are excluded for visibility purposes. Like Long et al. (2019), the selection ignored the millimeter brightness of disks. Six out of the 16 proposed targets were observed with the C43-9/10 configuration and Band 6 receivers between 2021 August 27 and 2021 September 27.

For these six sources, their stellar properties (listed in Table 1) were reevaluated based on optical spectra when available. For 2M0450 and 2M0508, the spectral type and extinction are obtained from Luhman (2018), with uncertain accretion properties. We obtained flux-calibrated low-resolution optical spectra for 2M0412, 2M0434, and 2M0436 using UH88/SNIFS (3200–10000 Å; Lantz et al. 2004) and for 2M0507 (CIDA 12) using Palomar-Hale 5 m/DBSP (3200–8700 Å, with a gap from 5600 to 6250 Å; Oke & Gunn 1982). The SNIFS data reduction is described by Guo et al. (2018), while the DBSP data reduction is described by Herczeg & Hillenbrand (2014). For these four targets, the spectral type, extinction, and accretion rate are calculated from a simultaneous fit, following approaches described in Herczeg & Hillenbrand (2008, 2014) and using a plane-parallel slab model for accretion developed by Valenti et al. (1993). The properties of 2M0412 were presented in Long et al. (2023).

The spectral type for each star is converted to temperature using the relationship of Herczeg & Hillenbrand (2014). The luminosity is then calculated from the 2MASS J-band magnitude (Cutri et al. 2003) and assuming zero J-band veiling, J-band bolometric corrections from Pecaut & Mamajek (2013), and J-band extinction from the A_J /A_V ratio developed by Wang & Chen (2019). All properties are calculated from distances obtained from inverting the parallax in Gaia Collaboration et al. (2023). The masses are then estimated using the Somers et al. (2020) evolutionary tracks for pre-main-sequence stars with 51% spot coverage. These tracks lead to higher masses than standard evolutionary tracks (see discussion for 2M0412 in Long et al. 2023).

For our ALMA observations, the receivers were configured into four spectral windows, with two centered at 217.0 and 234.4 GHz for dust continuum with bandwidths of 1.875 GHz and the other two centered at 220.0 and 230.5 GHz targeting ¹²CO, ¹³CO, and C¹⁸O J = 2–1 lines. Line channel intervals are spaced at 0.244 MHz (∼0.3 km s⁻¹). Table 2 shows the detailed ALMA observation log.

The raw visibilities were pipeline calibrated using the specified CASA versions (6.1.1 for 2M0412, CIDA 12, and 2M0508; 6.2.1 for the rest) for each object that can be found in the QA2 report (CASA Team et al. 2022). We identified residual features of atmospheric absorption correction around channel 500 in the 234.4 GHz spectral window for all six targets and flagged corresponding channels 400–600. For dust continuum imaging, we flagged the channels within 25 km s⁻¹ of the rest frame of CO lines and then averaged the data set with a channel width of 125 MHz, which is the recommended maximum bandwidth to avoid bandwidth smearing for ALMA Band 6.

These observations were intended to provide snapshots of these disks. ¹⁷ Due to the short on-source time (∼15 minutes), the peak signal-to-noise ratio (S/N) for each source ranges from 5 to 20. We attempted one round of phase-only self-calibration on sources 2M0412, 2M0434, and 2M0450. For 2M0436 and 2M0508, the emission morphology was largely altered after self-calibration, so we did not apply the solutions. CIDA 12 has too low S/N to improve from self-calibration. The self-calibration solutions are not applied to CO line emission channels, since the improvement is negligible.

The continuum images were generated using the ''tclean'' task with multiscale imaging, with the Briggs ''robust'' weighting parameters being 0.5 for 2M0434, 2M0436, 2M0450, and 2M0508; 1 for 2M0412; and 2 for CIDA 12. The adopted weighting parameters are compromises between resolution and sensitivity. The uv coverage of observations produced elongated beams (aspect ratio ∼2). For better visualization, we performed uv-tapering to get more circular beams. Images with original beams and uv-tapered beams are shown in Figure 1. The beam sizes are 003–007, and the rms noise is 30–40 mJy beam⁻¹, with detailed values summarized in Table 2.

Zoom In Zoom Out Reset image size

For line images, after subtracting the continuum emission using the ''uvcontsub'' task, we applied uv-tapering to enhance the S/N of the line images. The beams are tapered to 02 for 2M0412 and 2M0434 and 01 for the other four targets. We detected (or marginally detected) ¹²CO for all six targets and ¹³CO for 2M0412 and 2M0434. The channel maps in velocity and moment maps of detected ¹²CO and ¹³CO are shown in Appendix A.

The peak fluxes of our disks range from 0.19 to 0.47 mJy beam⁻¹, which corresponds to peak S/Ns of about 6–12 from the non-uv-tapered images in Figure 1. Dust emission from these disks is faint, as expected for low-mass M stars. To quantify the dust emission morphology without biases introduced by the image reconstruction, we fit the brightness profiles in the visibility domain. The model profiles are selected based on eye identification of disk substructures. For 2M0450 and CIDA 12, the function to describe their deprojected brightness is a centrally peaked axisymmetric Gaussian profile:

$$\begin{eqnarray}&&{I}_{\mathrm{GP}}\,(r)={10}^{f}\;\exp \;\left(-\displaystyle \frac{{r}^{2}}{2{\sigma }^{2}}\right)\end{eqnarray} \tag{ 1 }$$

where I_GP(r) is the Gaussian brightness profile as a function of disk radius r, with 10^f and σ being the amplitude and Gaussian width, respectively. For the other four disks showing cavities and rings, we modeled their morphologies with a radially asymmetric Gaussian ring whose inner and outer width can differ:

$$\begin{eqnarray}&&{I}_{\mathrm{AGR}}=\left\{\begin{array}{l}{10}^{f}\;\exp \;\left(-\displaystyle \frac{{\left(r-{r}_{\mathrm{peak}}\,\right)}^{2}}{2{\sigma }_{1}^{2}}\right)\ (r\leqslant {r}_{\mathrm{peak}})\\ {10}^{f}\;\exp \;\left(-\displaystyle \frac{{\left(r-{r}_{\mathrm{peak}}\,\right)}^{2}}{2{\sigma }_{2}^{2}}\right)\ (r\gt {r}_{\mathrm{peak}})\end{array}\right.,\end{eqnarray} \tag{ 2 }$$

where 10^f , r_peak, and σ_i are the amplitude, peak location, and ring width, respectively, for the inner and outer side. A radially asymmetric Gaussian ring has been previously used to describe the rings in transition disks (e.g., Pinilla et al. 2017a, 2018; Huang et al. 2020; Kurtovic et al. 2021). This model fitting is motivated by the accumulation of dust particles trapped in radial pressure bumps (Pinilla et al. 2017a). Dust particles in the outer disk are expected to take a longer time to grow and drift radially toward the pressure bump. Hence, the external width of the ring is expected to be larger than the internal width (e.g., Figure 4 in Pinilla et al. 2015).

Using the adopted radial profiles, we first generate corresponding face-on images. Each model image is projected with an inclination (i) and a position angle (PA) and combined with a spatial offset (δ_R.A., δ_decl.), which adds another four parameters. We generate each model image with a pixel size of 3 mas, which is far smaller than our synthesized beam (∼40 mas). We then use galario (Tazzari et al. 2018) to transform each model image into complex visibilities sampled at the same (u,v) points as those in the observations and calculate the ${\chi }^{2}={\sum }_{k}\,{w}_{k}| {V}_{\mathrm{obs}}\,({u}_{k},{v}_{k})-{V}_{\mathrm{mod}}\,({u}_{k},{v}_{k})\,{| }^{2}$, where w_k is the observed visibility weight; the likelihood function is then calculated as $L\propto \exp \,(-{\chi }^{2}/2)$. We assumed a uniform prior probability distribution for the parameters above (for inclination the probability $\propto \sin \;i$ to get a uniform disk orientation).

We sample the posterior probability distribution using a Markov Chain Monte Carlo (MCMC) method (emcee; Foreman-Mackey et al. 2013) with 50 walkers and 10,000 steps. After running MCMC, we get posterior samples after 10 times the integrated autocorrelation time to ensure stationary posterior samples. Results for each model parameter are shown in Table 3, with the reference values adopted as the median values of posterior samples and their uncertainties estimated from the 16th and 84th percentiles. The comparisons between observations and models in the visibility plane, the image plane, and deprojected azimuthally averaged radial profiles are shown in Figure 2. The brightness profiles of each model are shown in Figure 3. All best-fit model images show reasonable matches with observations, with no residuals above 5σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Dust Disk Model Parameters from uv Modeling

Source Name	2M0412	2M0434	2M0436	2M0450	CIDA 12	2M0508	Unit
Model	2AGR	GP&AGR	AGR	GP	GP	AGR
f1	${8.63}_{-0.03}^{+0.03}$	${9.83}_{-0.01}^{+0.01}$	${9.98}_{-0.03}^{+0.03}$	${10.45}_{-0.05}^{+0.05}$	${9.54}_{-0.08}^{+0.09}$	${9.93}_{-0.03}^{+0.03}$	log10 Jy sr−1
rpeak,1	${415.10}_{-6.41}^{+8.71}$	⋯	${84.68}_{-7.35}^{+6.05}$	⋯	⋯	${39.80}_{-4.54}^{+6.38}$	mas
σ11	${5.93}_{-4.21}^{+6.51}$	${268.38}_{-3.02}^{+3.04}$	${32.90}_{-6.23}^{+5.84}$	${16.89}_{-1.28}^{+1.44}$	${54.06}_{-11.15}^{+11.11}$	${4.69}_{-3.47}^{+5.94}$	mas
σ12	${117.71}_{-11.45}^{+13.71}$	⋯	${6.88}_{-4.63}^{+5.34}$	⋯	⋯	${41.91}_{-4.01}^{+3.61}$	mas
f2	${8.96}_{-0.02}^{+0.02}$	${9.95}_{-0.03}^{+0.03}$	⋯	⋯	⋯	⋯	log10 Jy sr−1
rpeak,2	${766.27}_{-14.94}^{+7.44}$	${67.41}_{-7.60}^{+8.91}$	⋯	⋯	⋯	⋯	mas
σ21	${37.05}_{-12.02}^{+7.31}$	${10.15}_{-6.71}^{+7.94}$	⋯	⋯	⋯	⋯	mas
σ22	${63.74}_{-5.68}^{+8.73}$	${56.22}_{-5.52}^{+4.81}$	⋯	⋯	⋯	⋯	mas
Inc	${15.90}_{-0.82}^{+0.70}$	${68.54}_{-0.23}^{+0.23}$	${53.46}_{-1.49}^{+1.41}$	${43.26}_{-10.37}^{+7.84}$	${65.53}_{-12.65}^{+7.43}$	${54.34}_{-1.78}^{+1.67}$	deg
PA	${124.14}_{-2.06}^{+2.81}$	${96.67}_{-0.24}^{+0.24}$	${46.96}_{-2.05}^{+2.03}$	${63.01}_{-12.87}^{+12.43}$	${169.94}_{-6.59}^{+5.22}$	${98.90}_{-2.10}^{+2.10}$	deg
ΔR.A.	${2.77}_{-2.06}^{+1.86}$	${10.89}_{-1.05}^{+1.03}$	${7.86}_{-1.15}^{+1.12}$	$-{3.13}_{-0.97}^{+0.96}$	${9.05}_{-3.30}^{+3.21}$	${12.36}_{-1.40}^{+1.38}$	mas
Δdecl.	$-{0.69}_{-2.12}^{+2.11}$	${10.37}_{-0.49}^{+0.50}$	${5.76}_{-1.11}^{+1.11}$	$-{8.98}_{-0.92}^{+0.91}$	$-{2.15}_{-6.18}^{+6.02}$	${5.77}_{-1.12}^{+1.10}$	mas
F1.3 mm	${17.58}_{-0.32}^{+0.33}$	${30.73}_{-0.24}^{+0.26}$	${2.71}_{-0.05}^{+0.05}$	${0.85}_{-0.03}^{+0.03}$	${0.60}_{-0.06}^{+0.06}$	${3.02}_{-0.09}^{+0.09}$	mJy
Mdust	${11.50}_{-1.15}^{+1.15}$	${19.31}_{-1.93}^{+1.93}$	${1.68}_{-0.17}^{+0.17}$	${0.53}_{-0.05}^{+0.05}$	${0.51}_{-0.05}^{+0.05}$	${2.61}_{-0.26}^{+0.26}$	M⊕
R68%	${118.0}_{-0.3}^{+0.5}$	${55.1}_{-0.5}^{+0.5}$	${11.8}_{-0.2}^{+0.2}$	${3.6}_{-0.3}^{+0.3}$	${13.8}_{-2.9}^{+2.8}$	${15.6}_{-0.5}^{+0.5}$	au
R90%	${126.0}_{-0.6}^{+0.9}$	${81.2}_{-0.8}^{+0.8}$	${13.2}_{-0.3}^{+0.4}$	${5.1}_{-0.4}^{+0.5}$	${19.6}_{-4.1}^{+4.0}$	${20.4}_{-0.8}^{+0.8}$	au

Note. Model row: AGR for asymmetric Gaussian ring, GP for Gaussian profile.

Download table as:  ASCIITypeset image

3.2.1. Disks with Substructures in Images

3.2.2. Potential Cavities around 2M0450 and CIDA 12

Gaussian profiles (Equation (1)) can only describe continuous disks with a monotonically decreasing radial profile. The Nuker profile, which was first introduced to describe the brightness profiles of elliptical galaxies (Lauer et al. 1995), can also reproduce ringlike emission with a central depression. The Nuker profile is formulated as

$$\begin{eqnarray}&&{I}_{\nu }\,(r)={10}^{f}{\,\left(\displaystyle \frac{r}{{r}_{t}}\right)}^{-\gamma }{\,\left[1+{\left(\displaystyle \frac{r}{{r}_{t}}\right)}^{\alpha }\right]}^{(\gamma -\beta )/\alpha },\end{eqnarray} \tag{ 3 }$$

where r_t is the transition radius, α is the transition index, β is the outer disk index, and γ is the inner disk index. The Nuker profile behaves as r^−γ at r ≪ r_t and r^−β at r ≫ r_t . The transition index α determines how smooth/sharp the transition is between these two asymptotic behaviors (see Figure 2 in Tripathi et al. 2017). With β > 0 and γ < 0 the Nuker profile can describe a ringlike morphology.

Appendix C presents the detailed modeling with Nuker profiles on 2M0450 and CIDA 12. Interestingly, the best-fit Nuker profiles for the two smooth-appearing disks show depleted inner emission with emission peak at ∼1.7 au for 2M0450 and 5.7 au for CIDA 12 (see Figure C1).

However, the Nuker profile has more parameters than the Gaussian profile, which results in a larger Bayesian Information Criterion (BIC) ¹⁸ value. Hence, we still treat 2M0450 and CIDA 12 as nonstructured disks in this work. Future higher angular resolution observations are needed to check their potential small cavities.

2M0436 shows a ring with a wider inner width, which is contrary to models of dust trapping in pressure bumps. One possibility is that an external object has truncated the outer disk, which we expect to be substellar since no stellar companions have been detected (Kraus et al. 2011; Kraus & Hillenbrand 2012). Constraining the mass of the substellar object given current knowledge would be highly uncertain and beyond the scope of this work. Other possibilities are the presence of central compact emission blending with the ring emission and the blending of multiple unresolved rings. Compact emission at the center is also hinted at because the model intensity of 2M0436 slightly underpredicts the emission of the very inner disk (right column in Figure 2), though this difference is not statistically significant. An example is CIDA 1, where a model considering only one ring results in a ring with wider inner width (Kurtovic et al. 2021), while higher angular resolution observations reveal a compact inner disk and the inclusion of an inner disk recovers a ring with wider outer width (Pinilla et al. 2021).

In Appendix B, we show the results of modeling 2M0436 with a ring and a central point. The flux of the added central point emission is ${0.04}_{-0.03}^{+0.03}$ mJy (corner plot statistics in Figure B2). The new radial profile increases the inner emission and matches the observations better than a ring-only model, while the ring still displays a wider inner width after including a point emission. However, the result of modeling with central point emission could be limited by current observations; deeper sensitivity and higher angular resolution are necessary for future observations.

The millimeter fluxes, masses, and sizes of the dust disk are not explicit parameters in the models and are inferred from fitting results. Since our models are axisymmetric, the total fluxes are obtained through F_ν = f_ν ( ∞ ), where f_ν (r) is the flux within radius r:

$$\begin{eqnarray}&&{f}_{\nu }\,(r)=2\pi {\int }_{0}^{r}\,{I}_{\nu }\,(r^{\prime} )\,r^{\prime} {dr}^{\prime} .\end{eqnarray} \tag{ 4 }$$

The dust disk sizes are defined to be a radius that encloses 68% and 90% of total fluxes, consistent with choices from previous studies (see, e.g., the review by Miotello et al. 2023). The continuum fluxes, sizes, and their uncertainties are then computed from the last 5000 steps of the chains, as shown in Table 3. Most adopted models reproduce lower fluxes than those seen at short baselines (left column in Figure 2), likely due to the fact that either large-scale emission is resolved out by high angular resolution beams or faint emission is buried by the noise, so that we lack extra components in our modeling that could account for the flux difference between models and observations. Future observations with deeper sensitivity and shorter baselines are needed to fix this issue.

For CIDA 12, Nuker profiles reproduce millimeter flux better than that from Gaussian profiles. The visibility near zero-spacing baselines indicates that its flux is above 1 mJy, and it is consistent with ${1.16}_{-0.09}^{+0.09}$ mJy reported in Akeson et al. (2019), specifically the Gaussian profile gives ${0.60}_{-0.06}^{+0.06}$ mJy, while the Nuker profile recovers a flux of ${1.20}_{-0.51}^{+0.98}$ mJy.

Under the assumption that the dust disk is optically thin at millimeter wavelengths, the dust disk mass M_dust is estimated following Hildebrand (1983) as

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}\approx \displaystyle \frac{{D}^{2}{F}_{\nu }}{{\kappa }_{\nu }\,{B}_{\nu }\,({T}_{\mathrm{dust}})},\end{eqnarray} \tag{ 5 }$$

where D is distance from Gaia DR3 (Gaia Collaboration et al. 2023), κ_ν is dust opacity for which we adopt a power law of κ_ν = 2.3 cm² g⁻¹(ν/230 GHz)^0.4 (Andrews et al. 2013), and B_ν is the Planck function where the T_dust is assumed at 20 K (e.g., Ansdell et al. 2016).

We assume 10% uncertainty in the source flux calibration (Diaz Trigo et al. 2019; see also Francis et al. 2020) to calculate the dust mass uncertainty. Table 3 summarizes the estimated dust masses. 2M0434 has the most massive disk in our sample, with a dust mass ${19.31}_{-1.93}^{+1.93}\,{M}_{\oplus }$, while CIDA 12 has the lowest dust mass with only ${0.51}_{-0.05}^{+0.05}\,{M}_{\oplus }$. We note that these dust masses serve as lower limits, since part of the disk could be optically thick in reality (see Section 4.1 for a brief calculation). The assumed dust temperature would also affect the estimated dust mass: for example, the dust mass would be higher/lower by ∼40%/50% if the actual averaged dust temperature is 15 K/30 K.

Gas disk observations and comparisons to dust disks are helpful to understand disk evolution as a whole. To compute the integrated fluxes of CO lines, extraction regions need to be determined. We first generate the cumulative flux distributions and determine the radii where CO emission no longer increases. Then, we draw circular regions as large as these radii and project the regions with the same inclination and PA as those of dust disks. CO fluxes are then retrieved from these elliptical regions. The uncertainties of the integrated flux are measured in nonemission channels with the same velocity range as those used to generate moment 0 images, as the standard deviation of integrated fluxes estimated from 1000 randomly distributed elliptical regions with the same shape as the above extraction regions. The fluxes and uncertainties of ¹²CO, ¹³CO, and C¹⁸O J = 2–1 lines are listed in Table 4, as well as the 3σ upper limits for the nondetections.

We detect (S/Ns of integrated fluxes above 3σ) ¹²CO in five disks and ¹³CO in 2M0412 and 2M0434. 2M0450 has a marginal detection of ¹²CO with an S/N of 2.7σ. The gas disk sizes are measured the same way as dust disk sizes, adopting the same inclination and PAs from the fits to the dust continuum emission. 2M0450 has the smallest gas disk, about 26 au, while 2M0412 has the largest gas disk, about 199 au. When compared to dust disk sizes, the ratio ranges from 1.6 (2M0412) to 5.1 (2M0450). Since the largest angular scale recovered by our observations is only ∼06, the gas emission at large radii is likely resolved out; cloud absorption could also result in underestimation of disk sizes (e.g., Long et al. 2022); by 10%–30% with a typical cloud emission line width of 1–2 km s⁻¹, we treat the gas radii R_gas,90% as lower limits. The measured gas disk sizes and size ratios compared to dust disks are summarized in Table 4.

Figure 3 shows the azimuthally averaged radial profiles of gas emission and the best-fit continuum model intensity profiles. CO emission around 2M0412 suffers from severe cloud absorption, so we show the profile extracted from the moment 8 map ¹⁹ instead. Among the six disks, 2M0412 shows a drop in emission inside the dust cavity, with peak emission located at ∼36 au for ¹²CO and 54 au for ¹³CO. The cloud absorption may affect these values. A tentative drop of ¹²CO is seen around 2M0508, with the peak location overlapping with the dust cavity size. No other clear evidence of gas emission depletion is found in the other four disks, under current ∼01 (02 for 2M0434) resolution for gas images.

Putting disk properties together with the disk population can help constrain global disk evolution (Manara et al. 2023). Previous studies have shown that the millimeter luminosity of disks scales with stellar mass (e.g., Andrews et al. 2013; Ansdell et al. 2016, 2017; Barenfeld et al. 2016; Pascucci et al. 2016). The left panel of Figure 4 compares the six disks in this work with other Taurus disks with detected substructures in the L_mm–M_* plane (Long et al. 2018; Hashimoto et al. 2021; Kurtovic et al. 2021; Yamaguchi et al. 2021; Jennings et al. 2022; Zhang et al. 2023; review by Bae et al. 2023), as well as with other Taurus disks from Manara et al. (2023). The properties of the six disks follow the correlation for the whole Taurus population (Andrews et al. 2013). The four structured disks are brighter than the two smooth disks. It could be that the two smooth disks simply formed with less dust, while it could also be due to the lack of substructures to trap the dust and prevent dust from drifting toward the star. However, if the potential small cavities shown in Section 3.2.2 are indeed present in 2M0450 and CIDA 12, either their trapping efficiency is lower or the formation time of dust traps is later than the brighter disks.

Zoom In Zoom Out Reset image size

If disks are optically thick, the low L_mm might not correspond to a low dust mass. To have a simple estimate of the optical depth for the two smooth disks, we estimate the disk midplane temperature following Huang et al. (2018a):

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}\,(r)={\left(\displaystyle \frac{{\phi }_{\mathrm{fl}}\,{L}_{* }}{8\pi {r}^{2}{\sigma }_{\mathrm{SB}}}\right)}^{1/4},\end{eqnarray} \tag{ 6 }$$

where σ_SB is the Stefan–Boltzmann constant and ϕ_fl is the disk flaring angle. A conservative ϕ_fl = 0.02 is used (Huang et al. 2018a). Then, the optical depth τ_ν is calculated using

$$\begin{eqnarray}&&{I}_{\nu }\,(r)={B}_{\nu }\,({T}_{\mathrm{mid}}\,(r))\,(1-\exp \,(-{\tau }_{\nu }\,(r))),\end{eqnarray} \tag{ 7 }$$

where B_ν is Planck's law for blackbody radiation. The conservatively estimated peak τ is ∼0.5 for 2M0450 and 0.15 for CIDA 12, so that 2M0450 is partially optically thick and CIDA 12 is nearly optically thin. We expect the above arguments regarding dust trapping to hold for CIDA 12 at least.

The L_mm–M_* relation has been suggested to be flatter for disks with inner cavities (Pinilla et al. 2018, 2020; pink line in left panel of Figure 4). However, those studies lack data for the very low mass stars and disks with small cavities (<20 au). Furthermore, the flatter relation is not observed in any individual star-forming region. The two fainter disks with small inner cavities (2M0436 and 2M0508) deviate significantly from this flatter relation. Conservative estimates of their optical depth yield maximum values of 0.6 for 2M0436 and 0.5 for 2M0508, suggesting that the two disks are partially optically thick. Boulder formation could also lower the observable millimeter flux (Pinilla et al. 2020). A combination of optical thickness and boulder formation in these disks might explain their deviation.

In addition to stellar mass, the millimeter luminosity also scales with disk size (Tripathi et al. 2017; Andrews et al. 2018b). The right panel of Figure 4 compares the six disks in this work with other Taurus disks (Tripathi et al. 2017; Long et al. 2019; Kurtovic et al. 2021) in the R_eff − L_mm plane with the best-fit scaling relation ${R}_{\mathrm{eff}}\propto {L}_{\mathrm{mm}}^{0.53\pm 0.12}$ from Hendler et al. (2020). ²⁰ Numerical simulations have shown that disks fall along the observed relation ${L}_{\mathrm{mm}}\propto {R}_{\mathrm{eff}}^{2}$ if disks are in the radial-drift-dominated regime, while for disks with strong substructures they follow a relation of ${L}_{\mathrm{mm}}\propto {R}_{\mathrm{eff}}^{5/4}$ (Rosotti et al. 2019; Zormpas et al. 2022). The two relations intersect at the top end of the size–luminosity relation (i.e., for bright and large disks; Zormpas et al. 2022).

Most of our six disks follow the observed relation ${L}_{\mathrm{mm}}\propto {R}_{\mathrm{eff}}^{2}$, except 2M0412, whose large size stands out at its luminosity regime, suggesting that a pressure bump was built early at large disk radii of 2M0412 (Long et al. 2023). 2M0436, 2M0450, CIDA 12, and 2M0508 fall at the faint end of the ${L}_{\mathrm{mm}}\propto {R}_{\mathrm{eff}}^{2}$ relation and hence are likely drift dominated. For the structured 2M0436 and 2M0508, either their substructures are too weak to retain the dust in their rings or their rings formed late when the disks have already evolved to the faint end of the ${L}_{\mathrm{mm}}\propto {R}_{\mathrm{eff}}^{2}$ relation.

The gas disk is universally found to be more extended than the dust disk (e.g., Barenfeld et al. 2017; Ansdell et al. 2018; Long et al. 2022). In our six disks, we find that the two largest dust disks (2M0412 and 2M0434) have an R_90,gas/R_90,dust of ∼1.5 and 1.9, respectively. CIDA 12 and 2M0508 have a ratio of ∼2.7. The two smallest dust disks (2M0436 and 2M0450) have the most extreme ratios of around 5, which are also among the largest values for the disk population (Long et al. 2022). Early formation of pressure bumps in the two largest dust disks 2M0412 and 2M0434 may explain their lower size ratios. 2M0412 shows a large ring at large radii of 114 au, while 2M0434's ring peaks at about 10 au, which is far smaller than its size of 81 au, indicating that its extended emission outside the ring holds unresolved substructures. Trapman et al. (2019) suggest that disks with R_90,gas/R_90,dust > 4 can only be explained with dust evolution and radial drift. This is likely what 2M0436 and 2M0450 have undergone owing to their such extreme size ratios, which is also consistent with more efficient dust radial drift expected in lower-mass stars (Pinilla et al. 2013). These size ratios should be considered as lower limits since they are affected by cloud absorption and the observations lack short-spacing data (see the figures in Appendix A); hence, the gas emission is probably not recovered as well as dust emission.

Current studies of substructures mostly come from disks around solar-like stars (e.g., Andrews et al. 2018a; Long et al. 2018; Cieza et al. 2021). Van der Marel & Mulders (2021) found that substructures are less common around M stars, which was mostly based on studies with intermediate spatial resolution (∼25 au in radius) and biased toward larger disks, while detecting substructures around M stars needs higher resolution since they are smaller. Currently only a few disks around M stars (especially mid- to late type) have been imaged at high spatial resolution. In Taurus, Kurtovic et al. (2021) surveyed at a resolution of ∼01 (14 au) for a sample of six disks around M4–M5 stars and found two disks with a cavity and one ringed disk. The sample in Kurtovic et al. (2021) is biased to the brightest disks in the corresponding stellar regime and hence favorable for substructure detection. Hashimoto et al. (2021) reported an asymmetric dust ring around a cavity around the M4.5 star ZZ Tau IRS with resolution ∼02. Other studies toward the Lupus and Ophiuchus regions have also detected a few structured disks around mid- to late M stars with resolutions down to ∼4 au (González-Ruilova et al. 2020; Cieza et al. 2021; van der Marel et al. 2022).

Cavity+ring seems to be the most common substructure around low-mass stars. In our sample, substructures are detected in four (2M0412, 2M0434, 2M0436, and 2M0508) of six disks. The substructures we identify are all cavities surrounded by rings, with varying numbers in individual disks. Combined with structured disks around very low mass stars collected in Pinilla (2022), the cavity frequency seems higher than that around solar-like stars. Taking disks in Taurus as an example (Hashimoto et al. 2021; Kurtovic et al. 2021; this work), high spatial resolution imaging revealed 7 disks with cavities out of 13 disks around M3–M6 stars. In comparison to stars earlier than M3 in Taurus, only 4 out of 32 disks analyzed by Long et al. (2019) show cavities (5/32 considering the reanalysis in Zhang et al. 2023). However, the sample around mid- to late M stars used here might be more biased than that in Long et al. (2019). Based on SED modeling (see below and Figure 5), two out of the six disks in this work have cavities >1 au, which is higher than the transition disk fraction of ∼8% within the full Class II disks (van der Marel et al. 2016; similar fraction is obtained when limiting to M3–M6 stars). Furthermore, disks in this work have been imaged with ∼005 compared to 01 in Long et al. (2019). A complete sample around mid- to late M stars and higher-resolution imaging of the inner regions of disks around solar-like stars are needed to evaluate any difference in their cavity occurrence rate in millimeter wavelength.

Zoom In Zoom Out Reset image size

Gap+ring pairs are only detected in the large double ring disk 2M0412. 2M0434 has a ring located at about 10 au, while its dust disk size is about 80 au; the very extended emission in between may hold shallow rings and gaps inside. Rings around 2M0436 and 2M0508 have radial widths comparable to the beam size and hence are not resolved—at higher resolution they might be resolved into multiple narrower rings (e.g., Facchini et al. 2020; Pérez et al. 2020). It is also possible that fainter rings at large radii are not recovered with current sensitivity. Observations with a spatial resolution that is smaller than the pressure scale height and with deeper sensitivity are needed to have complete characterization of substructure type and occurrence rate in disks around mid- to late M stars.

Smooth disks 2M0450 and CIDA 12 are candidates to host small cavities (Section 3.2.2). Though their visibility profiles do not show clear null points (Figure 2), small cavities below 5 au could remain hidden at our current angular resolution, especially for CIDA 12, which has an inclination of ∼65° (see analysis of high-inclination disks by van der Marel et al. 2022).

ALMA revealed a central depression of large grains in four structured disks in our sample. The SED complements ALMA imaging by providing information on small grains, since deficits of near-infrared emission indicate an inner disk clearing of small dust grains. Figure 5 shows the collected photometry and the modeled SEDs with modeling processes given in Appendix D. In short, the modeling uses a continuous disk component with a cavity sized R_μm,in and an outer radius adopting R_90%,dust retrieved from millimeter images (Table 3). The goal is not to have a stringent constraint on disk parameters, since the available infrared photometry is limited, but to have a rough sense of the cavity sizes of small dust grains and how they compare to the millimeter-sized grains observed by ALMA. Since we lack knowledge of other parameters describing disk density profiles (e.g., power-law index of surface density and pressure scale height), the cavity size is degenerated with other parameters, including the power index for dust mass density and the pressure scale height. Through experiments, the constrained cavity sizes of small dust grains can differ by a few astronomical units under different combinations of other parameters. However, the uncertainty on cavity sizes does not affect the discussions below.

For the four disks with cavities in millimeter images, three disks show evidence of inner dust clearing from their SEDs. 2M0412's SED is consistent with R_μm,in ∼ 0.9 au for small dust grains, which is far smaller than its ∼60 au millimeter cavity. 2M0434 and 2M0508 are consistent with R_μm,in ∼ 4 au cavities. This test of simple SED fitting suggests that micron-sized dust grains are not as depleted as millimeter-sized grains in the cavity, consistent with a morphology that is also seen around larger and more massive disks with both scattered light and ALMA images (Villenave et al. 2019). R_μm,in for 2M0436 is around 0.03 au, consistent with dust sublimation radius, so there is no evidence for clearing of small grains. There is evidence of an inner disk around 2M0436 from both its brightness profile (Figure 4) and its visibility fitting result (Section 3.1). Models with the same inclination as that from 2M0436's ALMA image have a bit of difficulty in reproducing the high flux at near-infrared. These suggest that 2M0436 may host a puffed up or misaligned inner disk (e.g., Dullemond et al. 2001; van der Marel et al. 2018, 2022). For the two smooth disks, 2M0450 shows no clear small grain cavity with R_μm,in consistent with dust sublimation radius, while CIDA 12 shows evidence of a 0.3 au cavity. Further evaluation of these potential cavities will require follow-up at much higher spatial resolution and sensitivity. Scattered light imaging would likely require the next generation of ground-based telescopes, since these disks are too faint to observe with current ground-based telescopes and the cavities are too small for JWST imaging.

Rings are the most common substructures in our sample, similar to surveys on disks around early-type stars such as DSHARP (Andrews et al. 2018a; Huang et al. 2018a) and the Taurus survey (Long et al. 2018). With the current resolution and sensitivity of our observations, the detected rings all appear axisymmetric, without significant azimuthally asymmetric features. Most rings encircle cavities, except for the outer ring of 2M0412 (see Long et al. 2023 for a detailed analysis of the origins of substructures of the 2M0412 disk). In short, planet–disk interaction is preferred over other mechanisms in 2M0412, with Saturn-mass planets capable of shaping the observed disk morphology. However, a combination of dead zones and photoevaporation cannot be ruled out. Below we focus on discussion of possible origins ²¹ for substructures of 2M0434, 2M0436, and 2M0508, which show small cavities around 10 au.

Condensation Fronts. Major volatiles in protoplanetary disks freeze onto dust grains from the gas phase as the disk temperature decreases outward. Across these regions, which are referred to as condensation fronts or ice lines, the dust opacity and critical fragmentation velocity are expected to change and further ring/gap substructures can be produced (e.g., Zhang et al. 2015; Okuzumi et al. 2016).

We consider an irradiated flared disk with a disk midplane temperature formulated as Equation (6), where the flaring angle ϕ_fl is assumed to be 0.02 (e.g., Dullemond et al. 2018; Huang et al. 2018a). The ice-line location of a volatile can be calculated from a stellar luminosity and the freezing point of the molecule. Since ice lines of abundant species like H₂O and NH₃ are close to the host star (<2 au) and are not resolvable by our observations, we focus on species with lower condensation temperatures in Zhang et al. (2015), which are clathrate-hydrated CO and N₂ (41–46 K), CO (23–28 K), and N₂ (12–15 K).

The comparison between ice-line locations and brightness profiles from best-fit models is shown in Figure 6. Interestingly, the peak emission radii are all close to the CO ice lines. However, inferring the ice-line radius is a challenging task. Besides the uncertainties in the estimation of disk midplane temperature (e.g., Liu 2021), ice lines can be thermally unstable and dynamically evolve on timescales from 1000 to 10,000 yr (Owen 2020). In simulations, ice lines have not been found to carve out cavities (see Pinilla et al. 2017b for a simulation around a Herbig AeBe star). Deeper observations of CO isotopologues are needed to better constrain disk midplane temperature, and future simulations tailored to M stars could together help determine the contribution of ice lines to the cavities/rings detected.

Zoom In Zoom Out Reset image size

Photoevaporation and MHD Wind. When the mass-loss rate by photoevaporation surpasses the accretion rate through the disk, a gas and dust gap will open in the disk. Recent photoevaporation models by Picogna et al. (2019) predict transition disks with accretion rates ≤10⁻⁹ M_⊙ yr⁻¹ and maximum cavity sizes of about 30 au when surrounding solar-type stars. The lower mass-loss rate at large radii around very low mass stars (Picogna et al. 2021) will likely lead to a smaller maximum cavity size, which we expect to be above 10 au. This is because the old photoevaporation models around 0.1 M_⊙ stars by Owen et al. (2012) have cavities ≤10 au and new mass-loss profiles by Picogna et al. (2019) are more efficient at removing material at larger disk radii, which will lead to larger cavities. We therefore expect that the maximum cavity size around very low mass stars will fall between 10 and 30 au. 2M0434's 10 au cavity, 2M0436's 12 au cavity, and their accretion rates of about a few × 10⁻¹⁰ M_⊙ yr⁻¹ fall in the predicted parameter space.

Further evaluating the role of photoevaporation in the creation of cavities will require additional observations. Possible discriminators include high angular resolution imaging of the gas—in this paper, we have limited sensitivity that allows us to only tentatively identify a cavity around one source, 2M0508. In addition, empirical measurements of photoevaporation rates through lines of O i and Ne ii (Pascucci et al. 2011) but the results are challenging to interpret (e.g., Banzatti et al. 2019; Fang et al. 2023b; Rab et al. 2023), given the possibility that these lines may form in a magnetothermal wind (e.g., Wang et al. 2019).

Furthermore, accreting transition disks could be sustained by other mechanisms. Gárate et al. (2021) showed that a disk with X-ray photoevaporation in combination with dead zones can reproduce both the accretion rates and gap sizes observed in transition disks. Their models predict a long-lived compact inner disk, which is not seen in our ALMA images under current resolution and sensitivity. MHD wind as another scenario, if the large-scale magnetic field distribution is not significantly modified, could generate a positive radial slope of the gas surface density (e.g., Suzuki et al. 2016), which inhibits the inward drift of pebbles, and sustain accretion rates similar to those of classical T Tauri stars (e.g., Lesur 2021; Martel & Lesur 2022).

Embedded Planets. Giant planets can open gaps in the gas disks and form gas pressure bumps that trap dust particles outside its orbit. Dust cavities will then form after the inner dust disk materials are accreted onto the star. By applying the analytical criterion for opening a gap in the gas (Crida et al. 2006), the estimated masses of planets needed to open the small cavities in our disk sample are as follows: for viscous α ∼ 10⁻⁴ and 10⁻³, 0.1M_Jup and 0.3M_Jup at ∼8 au for 2M0434, 0.1M_Jup and 0.3M_Jup at ∼10 au for 2M0436, and 0.1M_Jup and 0.2M_Jup at ∼6 au for 2M0508. The locations of planets are assumed to be five times the Hill radius ${r}_{{\rm{H}}}={r}_{p}{\,({M}_{p}/3{M}_{* })}^{1/3}$ away from the ring peak locations (e.g., Dodson-Robinson & Salyk 2011). Since only 2M0508 shows a tentative gas cavity, these planet masses should be treated as upper limits. Moreover, if the disk mass and angular momentum transport are dominated by MHD wind rather than turbulent viscosity, planets can open wider and deeper gaps more easily (e.g., Elbakyan et al. 2022; Aoyama & Bai 2023; Wafflard-Fernandez & Lesur 2023). As a result, the masses of embedded planets could be a factor of a few to 10 smaller than estimated above (Elbakyan et al. 2022).

Simulations of planet–disk interaction predict the segregation between millimeter-sized dust continuum emission peak and gas emission peak (or micron-sized dust; e.g., de Juan Ovelar et al. 2013; Facchini et al. 2018). We find hints of such segregation between micron-sized and millimeter-sized dust grains from SEDs (Section 4.2). For gas emission, only 2M0508 shows tentative inner gas depletion, and the gas peak location overlaps with the ring location. Deeper sensitivity and higher angular resolution of molecular line mapping are needed to resolve the gas depletion in the inner disk if present.

Taking disk dust masses in Section 3.4 and assuming a gas-to-dust mass ratio of 100, the estimated current disk masses are 6.1M_Jup for 2M0434, 0.5M_Jup for 2M0436, and 0.8M_Jup for 2M0508. If we assume a constant accretion rate after disk formation, then the total mass accreted onto the star from the disk is 1.5M_Jup for 2M0434 and 4.5M_Jup for 2M0436 (using stellar ages with 50% spot coverage). The planet masses needed to open cavities for 2M0434 and 2M0436 are within 10% of their estimated total disk masses at initial stages. If 2M0508 has an accretion rate ∼3 × 10⁻¹⁰ M_⊙ yr⁻¹, similar to those in our sample (Table 1), and that accretion rate has remained constant (most likely, the accretion rate at earlier stages of formation was higher; e.g., Fischer et al. 2017; Fang et al. 2023a), its potential planet mass is also within 10% of disk mass. Hence, from the perspective of disk mass budget, the three small cavity disks had enough materials to form those Saturn-mass planets (e.g., Boss 2006; Lin et al. 2018).

If the observed cavities are carved out by giant planets, how does the current detection rate of cavities around mid- to late M stars compare to the occurrence rate of exoplanets? The detection rate of cavities around mid- to late M stars in Taurus is around 24%, given the fraction of 7/13 for cavities in disks (Section 4.2) and accounting for the disk fraction around 45% for M3–M6 stars in Taurus (Esplin & Luhman 2019; the latter correction applies only if the age spread of the sample is much less than the disk survival timescale). However, this is likely higher than the real cavity rate, since the sample of 13 disks is slightly biased to the brighter end of disks around M3–M6 stars. As for exoplanets around M dwarfs, the CARMENES survey obtained occurrence rates for giant planets with ${M}_{p}\;\sin \;i\gt$ 100 M_⊕ and periods of 1–1000 days of ${0.021}_{-0.011}^{+0.018}$ planets per star and ${0.045}_{-0.016}^{+0.021}$ planets per star for stars less and more massive than 0.337 M_⊙, respectively (Ribas et al. 2023). Wide-orbit (>10 au) giant planets more massive than Jupiter are found to be rare around M dwarfs (e.g., Bowler et al. 2015). Currently, the constraints on sub-Jovian giants on orbits of a few astronomical units around M dwarfs are limited. Measurements around Sun-like stars show that giant planet occurrence is enhanced by a factor of four beyond 1 au compared to within 1 au (e.g., Fernandes et al. 2019; Fulton et al. 2021), while it is unclear whether the result could be extrapolated to M dwarfs. Samples from both disks and exoplanets need to be developed to make a robust comparison.

Binary Companions. Binary stars could also carve out inner cavities in circumbinary disks, with the cavity sizes expected to be 3–5 times the binary semimajor axis (e.g., Artymowicz & Lubow 1994; Miranda et al. 2017). Hence, stellar companions within 4 au could be responsible for the 7–12 au cavities in the three disks. To our knowledge, no evidence of stellar multiplicity has been found in our sample (Kraus et al. 2011; Kraus & Hillenbrand 2012), although the scales of a few au have yet to be explored. If the companions are sufficiently massive (mass ratios > 0.1–0.2), the binary would produce eccentric cavities (0.05–0.35 for circular binaries; Miranda et al. 2017; Ragusa et al. 2020) and cause a shift between the stellar mass center and the geometric center of the cavity. For the ∼10 au cavities in our sample, resolutions better than the binary separation are needed to detect the potential shift.

Among the mechanisms described above, photoevaporation and planets are more preferred origins of small cavities (around 7 and 10 au) in our sample. To distinguish between these possibilities, more deep and higher angular resolution observations of molecular lines are needed. Furthermore, multiwavelength analysis of the spectral index at the cavity edge might help to distinguish photoevaporation and planet scenarios (Picogna et al. 2023), where a cavity created by photoevaporation is expected to have a higher spectral index owing to its lower dust filtering efficiency.