© The Authors 2025

1 Introduction

The field of protoplanetary disks has been revolutionized in the past decade with the sensitivity and spatial resolution of the Atacama Large Millimeter Array (ALMA) telescope. Significant efforts have been dedicated to studying the demographics of protoplanetary disks in recent years with ALMA, providing a broader perspective on disk evolution and enabling direct comparisons and statistical analyses of exoplanets in conjunction with protoplanetary disks (Zhang et al. 2018; Lodato et al. 2019; van der Marel & Mulders 2021). Early ALMA survey programs yielded valuable insights into dust evolution in protoplanetary disks through various disk relations (e.g., Pinilla et al. 2020) despite being taken at relatively low resolutions (~0.25″; Ansdell et al. 2016; Barenfeld et al. 2016; Cieza et al. 2019). Subsequent high-resolution observations ($\sim 0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}04$) focused on the brightest disks, which are often large and harbor numerous extended substructures (Andrews et al. 2018a; Cieza et al. 2021), usually associated with pressure bumps halting the radial drift through dust traps (Pinilla et al. 2012). However, many disks remained unresolved in this initial reconnaissance, preventing a thorough investigation of the entire disk population.

The first discovery of a very small dust disk with highresolution ALMA observations was XZ Tau B, with a dust disk size of only 3.4 au and potentially an inner cavity, located within a binary system of 39 au separation (Osorio et al. 2016; Ichikawa et al. 2021); this is consistent with predictions from radial drift in binary systems (Zagaria et al. 2021). The first small disk in a single star system studied at high resolution was CX Tau, analyzed by Facchini et al. (2019). They found no substructure but measured a dust disk radius of 14 au, which was five times smaller than the CO extent, indicating efficient radial drift. In subsequent years, several additional studies (Long et al. 2019; Kurtovic et al. 2021; van der Marel et al. 2022; Miley et al. 2024; Shi et al. 2024) identified an increasing number of compact disks, defined as having dust disk radii smaller than 30 au for the purposes of this paper, which is also the size of Neptune’s orbit. We have identified 33 such compact disks in the literature to date, a few of which show centrally cleared cavities. However, the size distribution and structural classification of disks within a single star-forming region remain limited due to the limited number of high-resolution observations and accurate measurements of the disk radii for the smallest disks, which often remain unresolved.

One of the most important yet incomplete relationships for the most compact disks and within single star-forming regions is the size-luminosity relation (SLR). Assuming that the relation is not primordial, meaning that low-mass disks are born small while high-mass disks are born large, current literature identifies three distinct slopes that define the behavior and evolution of disks within this relation. Dust grain growth in disks is limited by two main barriers (Birnstiel et al. 2012): the drift barrier, largely driven by radial drift (Weidenschilling 1977), and the fragmentation barrier, shaped by turbulence which causes grain collisions and destruction (Voelk et al. 1980). These barriers determine how disks evolve along the SLR. Rosotti et al. (2019) studied the SLR in the context of these two primary dust growth barriers, introducing two distinct slopes in this relation. The first is the drift-dominated slope, where a disk’s position along the SLR is mainly influenced by its dust mass, which affects both its luminosity and radius. For disks with a smooth density profile, this slope follows $F_{mm}\propto R_{eff}^2$, where F_mm is flux and R_{e f f} its disk size, which is typically defined as the radius enclosing a specified fraction (often 68% or 90%) of the total flux from the disk. However, in fragmentation-dominated disks, the slope becomes steeper, following F_mm ∝ R_{e f f}, as fragmentation increases the mass retained in the disk, increasing the flux for the same radius compared to the drift-dominated scenario, thereby altering the disk’s position on the SLR. The third slope, described by Zormpas et al. (2022), is the trap-dominated slope, which implies that disks have strong dust traps. In this case, the SLR behaves differently compared to smooth, non-trap disks. This slope assumes that a disk begins its evolution with a planet forming within it, which shifts its evolution along the SLR. The relation for these disks follows $F_{mm}\propto R_{eff}^{5/4}$, which is steeper than the drift-dominated slope but less steep than the fragmentation-dominated one. More recently, Delussu et al. (2024) continued this exploration, suggesting that even smooth disks might have hidden substructures. They were able to reproduce the observed $F_{mm}\propto R_{eff}^2$ slope by having optically thick high-flux and optically thin low-flux disks with substructures.

Several factors contribute to deviations from the SLR, as noted by Zormpas et al. (2022). Dust properties, such as variations in opacity or porosity, can shift a disk’s position along the SLR by affecting its luminosity without changing its size. Additionally, the turbulence parameter α plays a significant role in modifying a disk’s location on the SLR. Both factors are essential for the drift and fragmentation barriers, which in turn are key to determining dust evolution and growth. This highlights the importance of studying the SLR in greater detail and at the highest possible resolution, as recent research has done.

Furthermore, constraining disk parameters and extrapolating the SLR to the smallest disks results in disk sizes of only a few au which has significant implications for understanding dust substructures and the comparison with exoplanet populations. Recently, there have been several efforts to directly link exoplanet observations to the properties of their birthplace, the protoplanetary disks. Key parameters have been studied, such as the available bulk mass in protoplanetary disks to form giant and terrestrial exoplanets (Manara et al. 2018; Mulders et al. 2021b), as well as the connection between disk substructures and exoplanet demographics (Lodato et al. 2019; van der Marel & Mulders 2021; Zhang et al. 2023). Specifically, van der Marel & Mulders (2021) aimed to understand whether disk substructures influence exoplanet formation scenarios or are linked to the observed disk dichotomy. Their work suggests that the majority of close-in rocky exoplanets around M stars likely formed in the more abundant smooth, compact disks in the absence of giant planets at wide orbits, which would have prevented radial drift (Mulders et al. 2021a). These compact disks can form planets, particularly super-Earths through pebble accretion, under the influence of substantial radial drift, concentrating sufficient material in the inner regions (Sanchez et al. 2024) under the assumption that the bulk of the initial dust mass in embedded disks (Class 0/I) decreases rapidly to the mass in the protoplanetary disk stage (Class II) through radial drift (Appelgren et al. 2025). However, due to the relatively low resolution of observations of protoplanetary disks, compared to the regions where most exoplanets are found, the relation between such disks and rocky exoplanets, if any, remains inconclusive.

In this work, we present a near-complete high-resolution survey of 73 Class II protoplanetary disks in the young, nearby Lupus star-forming region. The Lupus region, part of the Scorpius-Centaurus association, consists of several subgroups (Lupus 1–9), each associated with distinct molecular clouds. It is a young (1–2 Myr) and nearby (150–200 pc) region, comparable in proximity and age to Taurus (Comerón 2008). This study is based on new ALMA continuum observations as well as ALMA archival data. We measure dust disk radii down to 0.6 au, fit radial profiles to determine substructures and study the relation between the disk flux density and its size. We discuss the observations and data reduction in Sect. 2, analyze the continuum visibilities to measure disk sizes and substructure in Sect. 3, and discuss the implication of our results and a comparison with a model grid in Sect. 4. We summarize our findings in Sect. 5.

2 Observations

Our data were obtained from Cycle 9 observations with ALMA in Band 6. These observations, with project code 2022.1.00154.S and PI Nienke van der Marel, were carried out in ten execution blocks between August 16, 2023, and September 11, 2023. The total observing time was 8.89 hours, with ~7.25 minutes spent on each source. In total, 33 sources were observed: Sz65, Sz66, Sz77, Sz72, Sz74A, Sz74B, J15592523-4235066, J16002612-4153553, Sz130, Sz131, Sz81A, J16073773-3921388, J16075475-3915446, J16080017-3902595, J16084940-3905393, J16085324-3914401, J16085373-3914367, J16092697-3836269, Sz102, Sz106, Sz108B, Sz110, Sz113, Sz117, Sz88A, Sz88B, Sz90, Sz95, Sz96, Sz97, V1192 Sco, V856 Sco A, and V856 Sco B. The final datasets for each source include six spectral windows with frequencies in the range 219–234 GHz and a total bandwidth of 2 GHz each for the continuum, as well as six spectral windows in the range 220–231 GHz with a total bandwidth of 1.875 GHz dedicated to ¹²CO and ¹³CO line emission, all with a maximum recoverable scale between 0.37″ and 0.715″. Additionally, we incorporated archival ALMA data to complete the sample of Lupus disks at high resolution (project code: 2018.1.01458.S, PI: Yen, Hsi-Wei). These observations were conducted between July 18, 2019, and July 19, 2019, with a total integration time of 1.1 hours and a total of eight sources: Sz71, Sz69, J15450887-3417333, Sz98, Sz123A, Sz100, J16083070-3828268, and Sz73. These observations include six spectral windows to observe line emission, each with a bandwidth of 58.6 MHz, and only one spectral window for the continuum with a total bandwidth of 2 GHz, between 232 and 234 GHz. Finally, a Band 4 dataset was used for EX Lup, as no high-resolution observations were available at any other band (project code: 2017.1.00388.S, PI: Liu, Hauyu Baobab). These observations were conducted on November 11, 2017, covering a frequency range of 145 to 161 GHz, with a total integration time of 5.1 minutes. We searched for line emission to flag but found none, indicating only the presence of the continuum. However, these data were not used for any part of the analysis and served only to provide a high-resolution image of the source.

All datasets were calibrated using the pipeline and scripts provided by the ALMA staff. Version 6.5.4.9 of the Common Astronomy Software Applications (CASA; McMullin et al. 2007) was utilized to analyze and process the data, as well as to clean and create the final images. We separated the line and continuum emission spectral windows for each source and averaged the channels in the continuum spectral windows. Self-calibration was attempted on the new datasets, but due to the limited observing time per source and the low signal-to-noise ratios (S/Ns), it was not possible to achieve satisfactory results. We used the TCLEAN task in CASA with the Multi-term Multi-Frequency Synthesis (MTMFS) deconvolver (Rau & Cornwell 2011), employing scales of 0, 1, 3, and 5 times the beam size to create the final images of each disk. Several weightings were tested depending on the source, and different weightings were selected from the datasets, as shown in Table A.1. Natural weighting emphasizes short baselines, improving sensitivity but at the cost of resolution. In contrast, uniform weighting gives more weight to longer baselines, maximizing resolution but reducing sensitivity and making the images noisier while Briggs weighting serves as a compromise between natural and uniform, balancing sensitivity and resolution. For each source, we selected the weighting that offered the highest resolution while maintaining sufficient sensitivity to capture and distinguish the entire disk structure effectively. For all sources, pixel sizes between 0.001 and 0.003 arcseconds were used to ensure the pixel size was approximately ten times smaller than the beam size for all images.

To complete the high-resolution Lupus disk sample, we supplemented our data with images obtained from several previous studies. These included Sz103, Sz76, Sz104, Sz112, J16011549-4152351, J16081497-3857145, J16000236-4222145, J16090141-3925119, and J16070384-3911113 from (van der Marel et al. 2022); HT Lup A, B, and C, GW Lup, IM Lup, RU Lup, Sz114, Sz129, and MY Lup from the Disk Substructures at High Angular Resolution Project (DSHARP; Andrews et al. 2018a; Kurtovic et al. 2018); GQ Lup from Wu et al. (2017); and Sz91 in Band 4 from Maucó et al. (2021). Additionally, we used RY Lup data from Francis & van der Marel (2020). For already resolved disks J16070854-3914075, Sz118, Sz84, as well as for two sources without high-resolution observations (J16102955-3922144 and J16000060-4221567), we used images from Ansdell et al. (2018). For V1094 Sco, we used the image presented in van Terwisga et al. (2018), and for Sz111, data from Rota et al., in prep. (project code 2018.1.00689.S, PI: Muto, Takayuki). Finally, RXJ1556.1-3655 and Sz133 data were sourced from Bosschaart et al. in prep (project code: 2022.1.01302.S, PI: Mulders, Gijs). The total sample consists of 73 disks in Lupus with all but 5 disks observed at very high angular resolution of ≲0.05″. An overview of the characteristics of the images created in this paper is presented in Table A.1, with all the images of the disks shown in Figure 1. The most compact disks with radii <0.15″ are displayed on a smaller scale in Figure 2.

3 Results

3.1 Dust continuum images

In this sample of Class II Lupus disks with continuum images at high angular resolution we find a range of morphologies across all sizes. Although most substructures are detected in the more massive and larger disks, several compact disks also exhibit sub-structure, such as gaps and inner cavities. Given our resolution, which limits detection to structures as small as approximately 4 au (~ 0.03 arcseconds), we may still miss small inner cavities. With this in mind, the smallest cavity detected in the image plane is in Sz72, with a radius of 4.1 au, right at the edge of our resolution limit. Most disks are well resolved, however, a few remained only marginally resolved. Notably, Sz104, J16075475-3915446, J15592523-4235066, and J16084940-3905393 are very poorly resolved even at $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}03$ resolution, which implies a radii of less than 2 au. J16080017-3902595 and J16000060-4221567 are also unresolved at $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}25$ but lack the high-resolution data of the bulk of the sample.

In total, we resolved 11 new disks with substructures: 10 cavities (J16083070, Sz100, Sz123A, EXLup, Sz108B, Sz90, Sz72, Sz131, J16092697, and Sz96) and 1 very faint ring and gap (Sz73). We have a total of six resolved binary systems with disk detections and within the field of view of the observations; Sz66 – Sz65, Sz74A – Sz74B, V56ScoA – V856ScoB, Sz81A – Sz81B, J16085324 – 3914401-J16085373-3914367, and the HT-Lup triple system. However, no disk detections were found in three additional well-known binaries: Sz88B, Sz108A, and Sz123B-C (Zurlo et al. 2021). The total flux for all targets is determined through aperture photometry, where the aperture is extended until the flux flattens, and inclination corrections are applied. The uncertainty is calculated as the standard deviation inside each aperture. Most flux measurements fall within an 10% uncertainty range compared to those reported at a low resolution by Ansdell et al. (2018). Some discrepancies could arise from the absence of short baselines in our observations, potentially leading to the loss of extended emission flux. Overall, the fluxes appear consistent with previous results, meaning there is no need for a detailed reevaluation and we proceed to the analysis.

3.2 Visibility modeling

To analyze the exact orientation and morphologies of these disks, we used the GALARIO code (Tazzari et al. 2018) to model the visibilities of the observations with the best fit possible. We did not fit any features that were not already visible in the continuum images. For all isolated disks, we computed the visibilities from an axisymmetric brightness profile using different models. For each model, we fitted the visibilities, V_mod, following the approach outlined in Tazzari et al. (2018). We explored the parameter space of each model using a Bayesian approach, employing the Markov chain Monte Carlo (MCMC) ensemble sampler provided by emcee (Foreman-Mackey et al. 2013). Each model generates a brightness profile, which is then transformed and compared to the observed visibilities. Posterior distributions are obtained by assuming a Gaussian likelihood. We employed three different models for the isolated disks, choosing the most simple structure we observed for each disk in the continuum images.

The first is a Gaussian like disk, $I(R)=I_0\cdot \exp (-\frac{1}{2}\cdot {\left(\frac{R}{r_c}\right)}^2),$(1)

where I₀ is the peak intensity, R the radial distance from the center, and r_c the width of the Gaussian.

The second is a ring with different slopes on each side, $$\begin{array}{rl}I(R)=& I_0((1-\mathrm{\Theta }(R-r_{\text{ring\hspace{0.17em}}}))\exp \left(\frac{-{(R-r_{\text{ring\hspace{0.17em}}})}^2}{2{\left(r_{w}a\right)}^2}\right)),\\ & +I_0(\mathrm{\Theta }(R-r_{\text{ring\hspace{0.17em}}})\exp \left(\frac{-{(r-r_{\text{ring\hspace{0.17em}}})}^2}{2{\left(r_{w}b\right)}^2}\right))\end{array}$$(2)

where I₀ is the peak intensity of the ring, r_ring the peak position of the ring, r_wa the width of the ring on the inner side, r_wb the width of the ring on the outer side, and Θ the step function: $$\mathrm{\Theta }(x)=\{\begin{array}{ll}1& \text{\hspace{0.17em}if\hspace{0.17em}}x\ge 0\\ 0& \text{\hspace{0.17em}if\hspace{0.17em}}x<0\end{array}.$$(3)

With the first term active when inside the position of the ring peak, and the second term when outside the position of the ring peak.

Regarding the third, only one disk, Sz 73, required a two-component model to describe an inner Gaussian and a faint ring. We modeled this disk as follows: $$\begin{array}{rl}I(R)=& I_0\cdot \exp (-\frac{1}{2}\cdot {\left(\frac{R}{r_c}\right)}^2)\\ & +I_{0}b\cdot \exp {(-\frac{1}{2}\cdot {\left(\frac{R-r_{\text{ring\hspace{0.17em}}}}{r_{\text{width\hspace{0.17em}}}}\right)}^2)}^{,}\end{array}$$(4)

where the left Gaussian describes a ring with I₀b as the intensity at the peak of the ring, r_ring as the radial position of the peak of the ring and r_width as the width of the ring.

For the binaries, we sampled a 2D image and computed the corresponding visibilities, rather than using a single radial profile. We adopted a method similar to that of Manara et al. (2019), where the visibilities of the two sources were summed and compared to a Gaussian likelihood to obtain the final visibilities. We used the previously described Gaussian model, where the total visibilities are expressed as the sum of the Fourier transforms of each of the Gaussian models. This is mathematically represented as follows: $$V_{modT}=V_{modGauss1}+V_{modGauss2}.$$(5)

For two of the three binaries in this study, we employed a model consisting of two Gaussian profiles, which were added together in the visibility plane. The model, V_modGauss1 is a function of the brightness profile parameters and includes an offset relative to the phase center of the image.

V856ScoA was a special case where we see significant flux even on the longest baselines indicating an unresolved point source in addition to a Gaussian disk. To account for this, we defined a combined model as follows: $$V_{modGauss1}=\mathcal{F}(I_0\cdot \exp (-\frac{1}{2}\cdot {\left(\frac{R}{r_c}\right)}^2))+\mathcal{F}(\delta (R)),$$(6)

where δ(R) is defined as a Gaussian with an extremely small width, effectively approximating a delta function. For all the single disks, we used UVMODELFIT within CASA to obtain the initial parameter estimates. For the MCMC sampling, we employed between 48 and 120 walkers, running the chains for 2000–30 000 steps to ensure convergence, with the binaries requiring more walkers and iterations to achieve satisfactory results. To verify convergence, we discarded the first 1000 steps (burn-in phase) and examined the corner plots. The final model parameters are represented by the median of the posterior probability distributions for each parameter, while the uncertainties are given by the 16th and 84th percentiles.

Subsequently, we derived the total flux and dust disk radius from the modeled images. The results for the 38 sources modeled are presented in Table B.1 with the corresponding distribution of disk radii shown in Figure 4. We excluded Sz65, Sz66, J16083070, and Sz98, from the GALARIO analysis as these have been previously analyzed in Miley et al. (2024), Villenave et al. (2019), and Gasman et al. (2023). The corresponding model images are shown in Figure 3 and the visibility plots in Appendix D.

For most disks, the visibility fitting worked very well, but there were a few cases that required specific adjustments. The Gaussian model failed to capture the flux at the shortest baselines for Sz130, Sz77, and J16002612. In the case of Sz130, this resulted in a lower flux (41% lower) compared to the literature values. For Sz77 and J16002612, the model appeared to overestimate the flux by 9, 7% and ~3%, respectively. It remains unclear whether these discrepancies are due to the low S/N of the observations or if a more complex emission pattern (potentially with a different orientation) is present that a simple Gaussian model cannot capture. For Sz130, the visibilities suggest the presence of some extended structure or emission that shows below 300 kλ, which we were unable to fit or determine its nature. Nevertheless, we opted to proceed with the model, as it represents only what we could directly observe.

Three disks were unresolved or only marginally resolved, J16075475, J16084940, and J15592523. These disks exhibit extremely small sizes and very low surface brightness, leading to a low S/N. Consequently, even with visibility fitting, the inclination and position angle (PA) could not be reliably determined. To improve the fitting of other parameters, the inclination and PA were fixed to face-on values. As a result, the radius measurements for these sources should be interpreted with caution. However, the overall fit remained satisfactory, and we retained these results for analysis and discussion in this paper. In the same case, we were unable to fit the PA of J15450887-341733, which best fit was between 0.0 or 180 degrees. To obtain a plausible result, we fixed the PA to 0.0.

The faint ring in Sz73 was not detected during the initial imaging process with CASA, but became apparent when analyzing the visibilities. However, we concluded that this is unlikely for many sources, as nothing was evident in the visibilities. Even the faint ring in Sz73, with surface brightness between 0.2 and 0.4 mJy, was detected above the 3-σ level. Naturally, this depends on the observation brightness and the rms, but this means that fainter rings below this threshold would remain undetected.

In the Herbig star V856ScoA disk, a point source of strong emission is observed. The origin of this emission remains unknown and further discussion is needed. This emission, was detected at first only in the visibility plots, where the emission at long baselines never approaches zero, but it is observed in the continuum image. A model combining a Gaussian with a very narrow Gaussian successfully reproduced this feature, indicating that the emission is quite intense. Possible explanations include free-free emission from the central star’s outflow or wind emission from the disk near the star.

The Sz74 binary system shows a small substructure in the disk around its companion (Sz74B), potentially indicative of a cavity accompanied by a ring or asymmetry. However, due to insufficient resolution, it was impossible to constrain the size of the cavity or determine the disk’s PA, complicating the fitting process. To achieve a better fit, we opted to model the system using two Gaussians, as the substructure could also result from binary interactions, such as the formation of a vortex or a dust pile in the disk. While this approach improved the fit for the binary, the asymmetry remained un-modeled, and the PA was still unconstrained, suggesting the source was very unresolved. Similarly, five additional disks (Sz113, Sz69, Sz103, and Sz117) also show signs of cavities but were only fitted with a Gaussian profile since they were facing the same issue as Sz74B. Higher-resolution observations of these systems are necessary to confirm the origin of asymmetries, to better constrain the disk’s orientations or the potential presence of a cavity.

Fig. 1 All protoplanetary disks in Lupus at high resolution from various projects and observations with ALMA (see Section 2), shown at the same spatial scale. The disks are arranged in descending order based on their total flux. The scale bar and Neptune’s orbit in the final panel assumes that each disk is at a distance of 160 pc. Several cavities and substructures are observed, though the smallest disks are barely visible.

Fig. 2 Zoomed-in view of the smallest disks in the Lupus molecular cloud. Some of these small disks exhibit substructure, although most appear featureless. The disks are arranged in descending order of their total flux, consistent with the disks shown previously in Fig. 1. The scale bar and Jupiter’s orbit in the final panel assumes that each disk is at a distance of 160 pc.

Fig. 3 Best disk models obtained from the GALARIO visibility fitting. Note that the scale changes from Sz88 onward to better highlight the most compact disks. The disks are arranged in descending order based on their total flux derived from the visibility fitting.

Fig. 4 Size distribution of the Lupus disks. Disks exhibiting substructures in the continuum images are shown in dark blue, while smooth disks without visible substructures are depicted in red and added on top of the blue histogram. The pie charts inside the plot show the percentages of disks with and without substructures for sizes above and below 30 au. The vertical dashed orange line indicates our resolution limit, which is set at 0.03″ (4 au). More than 67% of the sample belongs to the compact disks classification.

3.3 Disk size and substructure distribution

The GALARIO visibility fitting allows us to explore the size distribution of disks in the Lupus star-forming region at much smaller scales than previous works. We define disk sizes as the radius where the cumulative flux equals 68% of the total flux (Andrews et al. 2018b; Tripathi et al. 2017). The histogram in Figure 4 shows a broad range of disk radii up to ~160 au and a marked increase at sizes smaller than 20 au.

We find that two thirds of the disks are compact, with radii <30 au. On the other hand, the majority of disks with large radii exhibit substructures, specifically at radii ~>30 au, although they are relatively few in number. However, as we approach the resolution limit (indicated by the dashed orange line in Figure 4, the number of smooth disks increases significantly. This is expected, as we would not anticipate resolving any substructures at approximately three times the resolution limit, which is about 12 au. Consequently, while we can now better measure disk sizes and note an increase in the number of smaller disks, there remains insufficient resolution to discern whether a trend exists that correlates disk size with substructure occurrence. In fact, the better the resolution, the fewer smooth disks appear in observations, which points to a potential bias introduced by resolution wherein very small substructures may remain undetected in the smallest disks. For large disks with R_dust > 30 au, approximately four-fifths exhibit substructure detectable down to scales of 7 au. Conversely, in the compact disks, at most 72% (5 out of 7) are smooth, with substructures down to scales of 4 au. Considering all disks in the sample, the data suggest that about 3 out of 8 exhibit substructure, regardless of disk size.

We plot the disk radii versus their total flux density (scaled to 160 pc) in Figure 5. There is a clear relationship between the two observables and we apply a Bayesian linear regression method implemented in the LINMIX package (Joshua E. Meyers, Kelly 2007), $${\log }_{10}(R_{68})=\alpha +\beta {\log }_{10}\left(F_{\mathrm{m}\mathrm{m}}\right).$$(7)

The fitted values are an intercept (normalization) α = 0.66 ± 0.06 and slope (power law index) β = 0.61 ± 0.06, with an intrinsic dispersion of 0.116 ± 0.02.

This relation can indicate whether dust evolution is dominated by drift or by traps being present in the disk (Rosotti et al. 2019; Zormpas et al. 2022). The observed slope falls between the drift-dominated and the trap- and frag-dominated regimes, which is consistent with previous findings by Andrews et al. (2018b) and Hendler et al. (2020). Notably, a few disks lie significantly below the SLR of Lupus. Most of these disks are in binary systems (star markers in Fig. 5), and the interaction between companions likely truncates their radius, as discussed in Rota et al. (2022), affecting their position on the SLR. Another exception is RU Lup (square in Fig. 5), one of the most active T Tauri stars, which likely experiences increased flux caused by its strong accretion features (e.g., Gahm et al. 2008; Siwak et al. 2016), placing it below the SLR as well. Of the five disks that appear high above in the SLR, two (J16102955-3922144 and J16000060-4221567) have only been observed at low resolution (Ansdell et al. 2018); therefore, their radius might be overestimated. The other three disks have large cavities and may be more evolved. As mentioned by Zormpas et al. (2022), the presence and location of planets within a disk can alter the SLR by influencing the disk’s evolutionary path. Given that the SLR represents only a snapshot of a cluster, with a diverse mix of stellar masses and angular momenta, it remains challenging to draw definitive conclusions about the outliers.

Fig. 5 SLR for the Lupus star-forming region at high resolution. Disks with substructures are marked by purple circles outlined in dark blue, while smooth disks are represented by solid purple circles without outlines. Unresolved disks are indicated by downward-pointing triangles. Binary disks are marked with star symbols, and RU Lup is identified with a square marker for clarity. The drift-dominated slope, as described by Rosotti et al. (2019), is illustrated by the dashed red line, the trapdominated slope (Zormpas et al. 2022) in green, and the frag-dominated slope in black. Our fit, using the LINMIX package, is displayed in turquoise. All disks are normalized to a common distance of 160 pc for consistency.

3.4 Radiative transfer modeling

While disk dust masses are usually calculated using a linear relation with the millimeter-flux (e.g., Ansdell et al. 2018), this relation relies on the assumption that the dust continuum emission is optically thin. The existence of very small disks raises questions about how well we can determine their disk masses as the optical depth may be very high. We therefore created a large grid of models to calculate the expected total flux densities for different stellar parameters and disk masses spanning the range of observed radii. We performed radiative transfer at 1.3 mm for the generic protoplanetary disk model in the RADMC-3D software package (Dullemond et al. 2012) with a dust density distribution, $$\rho (r,z)=\frac{\mathrm{\Sigma }(r)}{\sqrt{2\pi }{H}_p}\exp (-\frac{z^2}{2H_p^2}),$$(8)

where r represents the radial distance from the star, H_p(r) is the disk vertical scale height, and Σ(r) is the dust surface density, $$\mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0{\left(\frac{r}{r_{\mathrm{o}\mathrm{u}\mathrm{t}}}\right)}^{-1},$$(9)

defined out to an outer radius r_out.

We produced a grid of 1728 models, varying three key parameters: 12 disk dust masses, 12 radius, 12 stellar luminosities (L_⊙). These parameters ranged from 4e-6 to 5.1e-3 M_⊙, 0.5 to 30 au, 0.0025 to 3 L_⊙, with values spaced in a log-uniform manner, respectively. Unsurprisingly, the star’s effective temperature has a minimal impact on the flux-derived dust mass of the models as the absorbed stellar radiation is fully reprocessed by the disk to an equilibrium that is determined only by the energy input and dust distribution; therefore, we did not include it in the parameter space. The parameter limits were based on Alcalá et al. (2017) for Lupus sources, later updated using Gaia DR2 distance corrections as described in Alcalá et al. (2019). The stellar parameters of the Lupus sources, where known, are provided in Tables C.1 and C.2.

To calculate the dust opacities needed inside RADMC-3D, we utilized the optool software (Dominik et al. 2021) with DSHARP opacities (Birnstiel et al. 2018). The dust opacity values were calculated using grain sizes from a_min = 0.050 μm to a_max = 3 mm, incorporating the full scattering matrix.

For each model, we created a ray-traced image assuming a distance of 160 pc and a face-on disk inclination and measured the total flux and R₆₈ radius in the same way as for the observations. We first compared the disk masses as derived from the simulated flux, M_flux, with the model input, M_model, in Figure 6. For the former, we applied the optically thin approximation for dust mass at 1.3 mm simplified in Ansdell et al. (2018), $$M_{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}}=\frac{F_{\nu }{d}^2}{{\kappa }_{\nu }{B}_{\nu }\left(T_{\text{dust}}\right)}\simeq 0.68{\left(\frac{d}{160\mathrm{p}\mathrm{c}}\right)}^2\left(\frac{F_{1.33\mathrm{m}\mathrm{m}}}{\mathrm{m}\mathrm{J}\mathrm{y}}\right)M_{\oplus }$$(10)

where the dust opacity, κ_ν = 2.3 cm² g⁻¹ (Beckwith et al. 1990), B_ν is the Planck function, and we assume a uniform dust temperature, T_dust = 20 K.

Figure 6 primarily highlights the saturation of the disk flux as the emission becomes optically thick, demonstrated by how the models deviate from the dashed gray line, which indicates equality of M_flux and M_dust. Indeed, the flux-derived mass scales in a systematic way with the stellar luminosity for each input mass. Small disks saturate at lower disk masses than larger disks due to their much higher average optical depths. Motivated by this behavior, we fit power laws to the model grid away from the saturated region, defined through the mean optical depth, $$\overline{\tau }=\kappa \overline{\mathrm{\Sigma }}=\frac{\kappa M_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}}}{\pi R_{\mathrm{o}\mathrm{u}\mathrm{t}}^2}\simeq 19.4\left(\frac{M_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}}}{M_{\oplus }}\right){\left(\frac{R_{\mathrm{o}\mathrm{u}\mathrm{t}}}{\mathrm{a}\mathrm{u}}\right)}^{-2},$$(11)

which holds independently of the surface density profile and the normalization is for our mean observing wavelength, 1.33 mm. For very optically thin emission, $\overline{\tau }<0.25$, $$M_{\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}}\simeq 7{\left(\frac{R_{\mathrm{o}\mathrm{u}\mathrm{t}}}{1\mathrm{a}\mathrm{u}}\right)}^{-0.5}{\left(\frac{L_{\ast }}{1L_{\odot }}\right)}^{0.3}{M}_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}},$$(12)

to within 10%. As $\overline{\tau }$ increases, the behavior begins to deviate from a linear dependence on M_dust and is a more complicated function of radius and luminosity. For optically thick emission, $\overline{\tau }$ > 1, the mass dependence is so weak that the flux is no longer a reliable measure of dust mass. This presents a chicken and egg problem because we do not know how accurately we can determine the disk mass without knowing the optical depth and vice versa. However, the SLR provides a way forward.

Fig. 6 Comparison of the disk dust mass from a grid of 1728 models versus the mass derived from flux for the same models. In the top panel, stellar luminosity and effective temperature are fixed, illustrating how disk radius varies with flux and dust mass (see the colorbar with the radius values). The bottom panel shows the effect of changing stellar luminosity while keeping the disk radius and effective temperature constant. The dotted gray line highlights the masses from the flux that equal the disk dust mass introduced in the models.

3.5 Dust mass estimates

As shown in Sect. 3.3, our survey of compact disks extends the SLR to much lower disk radii than previous work. For a given stellar luminosity, our model grid calculates the flux for a given disk mass and radius. We then invert this relationship to map the disk mass on the flux-radius plane in Figure 7, which allows a direct comparison with the two principal observables in our survey.

The locus of the model shifts horizontally with the stellar luminosity but in principle, this allows the disk mass to be estimated from our model. Relatively large, faint disks have low $\overline{\tau }$ and the millimeter flux scales with the mass. However, disks that are relatively bright for their size may be optically thick and very massive. A handful of outliers that lie on or beyond the yellow region of this plot are binaries with relatively low radii for their flux due to tidal truncation and are not well represented by our models. However, most disks lie in the regime of intermediate optical depth, 0.25 < $\overline{\tau }$ < 1, where radiative transfer modeling is necessary to determine their dust mass.

The observations extend in the same direction as the $\overline{\tau }$ contours. For small optical depths, the flux is proportional to mass so, based on Equation (11), this is similar to the relation found by Andrews et al. (2018b) who found that the millimeter luminosity scales with the disk surface area. They suggested a possible explanation in which the emission comes from optically thick substructures with a filling fraction of ~0.3.

Using the stellar luminosities tabulated in Appendix C, we can create a bespoke radius-flux grid for each disk and thereby determine its dust mass. The uncertainties in each observable are readily propagated through. Several disks lie at the edge of their grids (Sz74A, Sz113, Sz102, J16084940-3905393, J15450887-3417333, Sz81A, Sz112, Sz88, Sz69, Sz130, HT LupA, Sz104, J16002612-4153553, and J16085373-3914367) due to being part of a binary system, having a poorly constrained luminosity, or being highly optically thick with $\overline{\tau }$ ≫ 5, resulting in inferred masses that were either indeterminate or exceptionally large with substantial uncertainties. The distribution of dust masses for the 28 disks with well determined values and known luminosities and the cumulative distribution of all Lupus disks are presented in Figure 8 and Table 1. Additionally, the τ values are plotted in Figure 9.

Fig. 7 Comparison of the model grid with observed disk radii and total flux (normalized to 160 pc). The image represents the model dust mass on a log scale (also in brown contours), and the black contours show the mean optical depth, $\overline{\tau }$, defined in Equation (11), with values of 0.25 and 1. The dark blue edges indicate a visible substructure in this source and a star-shaped marker distinguishes binary disks.

3.6 The exoplanet population in the Lupus substructures

To gain insight into the potential exoplanet population within the Lupus protoplanetary disk sample, we estimated the planetary masses that could dynamically cause each of the observed gaps. We employed DBNets, a deep-learning tool that utilizes convolutional neural networks to analyze observations of dust continuum emission and predict the mass of the gap-opening planets, as developed by Ruzza et al. (2024). DBNets takes several inputs: the continuum emission image of the disk (in the form of a FITS file), the center of the image (in pixels), the disk’s orientation (inclination and PA), the distance to the source, and the position of the planet (in astronomical units), which corresponds to the center of the gap being studied. In total, we analyzed 25 substructures within our sample. For well-studied disks such as RULup, IMLup, Sz129, Sz114, and GWLup, we used the planet masses already estimated in Ruzza et al. (2024).

For the remaining 16 disks and 17 substructures, we ran the code using the continuum FITS files provided in Section 2 (11 from the new images and 5 from other projects). For cavity positions we used $\frac{r_{\text{ring\hspace{0.17em}}}-r_{w}a}{2}$ while for the only gap, we used $\frac{r_{\text{ring}}-r_{\text{widh}}+r_c}{2}$. Inclinations and PAs, which were not obtained in Section 3 of this paper, were sourced from previous studies. Specifically, for Sz98, we adopted values from Gasman et al. (2023), for J16083070-3828268 from Villenave et al. (2019), for Sz91 from Maucó et al. (2021), and for Sz76 and J16090141-3925119 from van der Marel et al. (2022). For RYLup and Sz111, we referred to van der Marel et al. (2018). The results, obtained for a range of α-turbulence values between 10⁻²−10⁻⁴, including predicted planet masses, their radial positions, and whether the uncertainties exceed DBNet’s rejection threshold, are summarized in Table 2. It is important to note that this rejection threshold defines the significance of the results, meaning these values should be interpreted with caution. Additionally, in Figure 10 we show the full population of known exoplanets around M and K stars (sourced from the NASA Exoplanet Archive in November 2024) alongside our predicted exoplanet population in the Lupus star-forming region. Most exoplanets in the archive were discovered via radial velocity (RV), providing minimum mass values (Msini), and transit timing variations, which yield more precise masses when combined with RV data. The disk dust masses obtained in the previous section, where applicable, and the dust masses derived from the flux for the larger disks are plotted at the bottom of Figure 11 together with the planet masses for comparison. Moreover, we show the ring peak position beyond each gap and gap widths (this work and the DSHARP disks; Zhang et al. 2018) at the top of Figure 11.

From Table 2 we see that the range of inferred planet masses in Lupus, ~20–2000 M_⊕, is comparable to similarly derived estimates from observations of disks in Taurus (Zhang et al. 2023) and the σ Orionis cluster (Huang et al. 2024a). In addition, our high-resolution images reveal gaps at such small radii that the inferred planets overlap with the cluster of RV-detected exoplanets.

Fig. 8 Top panel: masses of the compact disks in Lupus (R < 30 au). These masses were obtained through the RADMC-3D model grid for τ < 5. The dashed black line indicates where the RADMC-3D dust mass equals the optically thin dust mass calculation (see Eq. 10). Top-right panel: Kernel density estimation of the disk dust masses. Bottom panel: cumulative distribution functions (CDFs) of all the disk masses in Lupus, using two approaches: in teal, all masses are estimated with the optically thin approximation based on flux values (Eq. (10)), and in pink, we substitute the masses with values derived from our model grid where available.

Fig. 9 Optical depths of the Lupus disks using a kernel density estimation for points with radii <30 au and obtained through Eq. (11) and the grid measurements.

Fig. 10 Lupus exoplanet masses obtained using DBNets plotted against the semimajor axis, assuming that the gaps are carved by planets. In blue, we show all observed exoplanets around M and K stars from the NASA Exoplanet Archive catalog, while the orange markers represent estimates for the Lupus sample. Brown outlines indicate cavities as the type of substructure, while green outlines denote gaps and rings. Additionally, planet mass estimates from Zhang et al. (2023), Huang et al. (2024a), and Ruzza et al. (2024) are highlighted in crimson circle markers.

Fig. 11 Top panel: gap widths and locations of the ring peaks beyond each gap in our disks, as well as the DSHARP disks (Zhang et al. 2018) in Lupus. Gaps are shown in teal, and cavities in pink. Bottom panel: disk dust masses versus planet masses from the DBNets analysis, using the same colors to denote the substructure type.

4 Discussion

4.1 New insights into the substructures and origin of compact disks

Our high-resolution compilation of the Lupus disk population shows a large number of compact disks with a relatively small fraction of substructures, whereas substructures are common in large disks with radii >30 au, with 79% of them displaying rings or gaps (19 out of 24). However, since the majority of Lupus disks are compact, this fraction may not be representative of the entire disk population. Many previous studies have shown that initially featureless disks reveal greater complexity when the resolution is increased (Long et al. 2019; Kurtovic et al. 2021; Shi et al. 2024; Yamaguchi et al. 2024). Here, we have detected new substructures, mainly cavities down to 4 au in radius, but only in 24% of the compact disks, excluding those that remain unresolved even in the visibility analysis.

Most of the structured compact disks show cleared cavities rather than rings, with the notable exception of Sz73. At a lower resolution, van der Marel et al. (2018) found that 15%¹ of the Lupus disk population show large cavities of >20 au radius. Our high-resolution survey increases the resolved transition disk fraction to 24.6% (18 out of 73), including cavities down to 4 au radius.

Sz73 is the sole structured compact disk with a ring. This is faint, between 0.2 and 0.4 mJy, but should be detectable in most of the other disks in the survey. Furthermore, by utilizing GALARIO, we reduced analysis biases and are therefore reasonably confident that there are few, if any, other Sz73-like objects in Lupus.

Our results thus confirm that many of the unresolved disks from Ansdell et al. (2018) are in fact very compact, consistent with the hypothesis from van der Marel & Mulders (2021). They speculate that disks that did not form giant planets early on, and thereby did not halt the pebble drift, will rapidly become compact, following dust evolution models (Pinilla et al. 2020; Appelgren et al. 2025). Subsequently, Sanchez et al. (2024) demonstrated that such drift-dominated disks are capable of producing the close-in Super-Earth population around M dwarfs through pebble accretion.

The observed properties of the compact disks in Lupus are consistent with this scenario. Most of our compact disks (73%) are found around M-type stars, with only 14% around K-type stars (the rest do not have a determined Spectral type). However, this distribution reverses in larger disks: 50% are around M-type stars, but already 41% hosted by K-type stars. Furthermore, Sanchez et al. (2024) showed that the disk dust masses in their models are a few M_⊕ and a few au in size after 1–2 Myr due to the inward drift in combination with inner planet formation through pebble accretion. The observed dust masses are in that case the remnant of planet formation and not representative for the solid mass budget for planet formation. A fraction of the formed planets in their simulations reach pebble isolation mass, creating small dust traps, which may be the explanation for some of the observed small cavities in compact disks. While forming wide gaps in the first 10 au of the disk is challenging, Kuwahara et al. (2024) proposed that low-mass planets could carve such gaps in the inner disk (<10 au), with widths as large as their position. For compact disks without observed substructure, either planets have not reached pebble isolation mass or they are located so close in that the dust traps remained unresolved. Predictions of the exact scenarios for individual compact disks will be explored in a follow-up study (Guerra-Alvarado et al., in prep.).

Also, in the sample of large, well-resolved disks in Lupus, there are several disks that do not show substructure (~25% of the disks with radius >15 au) (e.g., J16011549, J16000236, Sz133, Sz65, RXJ1556.1, GQ Lup, J16081497, J15450887, and HTLupA). Ribas et al. (2023) identified a disk around MP Mus without any substructure at 4 au resolution and proposed that in such evolved systems, large grains should have drifted onto the star unless some mechanism was preventing this. This means that substructures could remain undetected due to high optical depth at 1.3 mm or the substructures could be smaller than the current resolution limits of the observations. This could also be the case for some of the disks in the Lupus sample. An alternative possibility for the lack of substructure could be that the gas surface densities are high enough to drag the millimeter-emitting dust grains along such that there is little radial drift. If this were the case, the total disk masses would be ~10–20% of the stellar mass, near the limit of gravitational stability (Williams et al. 2024) and the mass of millimeter-emitting dust grains is a small fraction, ~5–10%, of the total solid mass. Testing this scenario requires more sensitive CO observations to measure the gas disk radius.

For compact disks specifically, Tong & Alexander (2025) found that such disks can form in magnetorotational-instability-active regions beyond the dead zone. In these regions, and particularly below ~6 au, dust traps are ineffective, and the optical depth of the dust emission obscures disk substructures.

4.2 Size–luminosity relation

The millimeter luminosity in Lupus, which increases with radius as $$R_{68}^{1/\beta }$$ where β = 0.61 ± 0.06, is consistent with the findings of Andrews et al. (2018b) and Tripathi et al. (2017) for a compilation of disks, as well as with Hendler et al. (2020) for Lupus. These results align with the drift-dominated scenario described in Rosotti et al. (2019). This also agrees with the super-Earth formation scenario from Sanchez et al. (2024), which requires significant radial drift. There is no clear distinction in the SLR for disks with or without observed substructures, as their radii and millimeter fluxes appear scattered. When fitted independently, we found that the slope of disks with substructures (β = 0.47 ± 0.09) is the same as for disks without substructures (β = 0.55 ± 0.08) within the error bars, and the latter group is not favored more than the former for the expected drift-dominated slope.

While our slope aligns with the drift-dominated scenario, there are still some outliers that fall below the drift-dominated slope, which we did not address in Section 3. Rosotti et al. (2019) predicted a population of fragmentation-dominated disks (lying closer to the fragmentation-dominated slope), which could explain some of the disks outliers. However, various factors may also contribute to these deviations in the SLR (Zormpas et al. 2022). For instance, dust properties, such as variations in opacity or porosity, and changes to the turbulence parameter, α, can shift a disk’s position along the SLR and potentially account for some of the observed outliers. Specifically, higher α values increase the luminosity, explaining disks positioned higher on the SLR. Variations in opacity not only affect luminosity but also influence disk size as they evolve. Lastly, extremely high porosity can lead to disks with low luminosity while causing minimal changes in disk size, which could explain the disks located at the lower end of the SLR.

More recently, for Taurus, Long et al. (2019) and Shi et al. (2024) also found an SLR with a slope consistent with drift-dominated disks by including more high-resolution observations. Nevertheless, they also found two disks, below 5 au, that fall below this slope.

4.3 Dust masses and optical depths

We create a grid of radiative transfer models that predict the flux for a range of disk radii, dust masses, and stellar luminosities. Comparing our observations of size and flux, we then estimate disk masses, finding a range from 0.3 to 26.3 M_⊕. The average optical depth for these disks shows a narrower variation and clusters around 0.7 for most disks. Most of the masses are consistently higher than those from the optically thin approximation based on flux, but they remain relatively comparable and low. These new masses remain too biased to observe any trends or changes in disk relations, such as M_dust − R_dust and M_dust − M_*. When considering all disks, these relations remain unchanged. Additionally, comparing the optical depths of these disks with other parameters does not reveal any clear trends. This likely indicates that relations like M_dust−R_dust and M_dust−M_* are not significantly impacted by associated optical depth effects. While we might have expected many compact disks to be optically thick, our findings indicate otherwise. A similar pattern was observed in the rings of the DSHARP survey (Dullemond et al. 2018), where optical depths clustered between 0.2 and 0.5, rather than being optically thick. Stammler et al. (2019) showed that if the dust density is regulated by planetesimal formation (small grains turning into large objects), typical millimeter continuum optical depths are ~0.5. This means that the dust mass distribution in Lupus is likely not representative for the solid mass budget for planet formation, indicating that the bulk of the pebbles have already been sequestered in larger bodies, perhaps even planets, especially when considering efficient pebble accretion.

4.4 Impact of compact disks in the exoplanet population

Under the assumption that the observed rings and gaps are produced by planets, we estimate their masses using DBNets. We are in agreement with the planet masses found in previous studies, including Zhang et al. (2023), Huang et al. (2024a), Ruzza et al. (2024), and Kurtovic et al. (2021). However, we extend their findings to include smaller substructures, with gaps as small as 4 au. As previously mentioned, the planet masses and semimajor axes for the smallest disks in our sample align closely with those of exoplanets primarily detected through RV methods and where the occurrence rates of giant exoplanets peak between 1 and 10 au (Fernandes et al. 2019; Fulton et al. 2021). This strengthens the case for planets forming within the observed gaps in several protoplanetary disks. Furthermore, the overlap between these gap radii and the peak in the exoplanet occurrence suggests that at least some exoplanets may form in situ in compact disks rather than migrating to their observed location or that they have already migrated to at the time the compact disks are observed.

In Fig. 10, we observe an intriguing lack of planets and/or gap centers in Lupus disks between 5–15 au, this gap remains between 5–8 au even with additional data. Although this needs to be properly quantified, since there is additional uncertainty in the case of cavities because the planet’s location can be any value within the cavity and we need more observations of disks at high angular resolution to conclude that this “valley” is significant, we can still speculate about the possible causes if this gap is real. For this to occur, a break in the disk properties is necessary in these two regions. While there is no straightforward explanation at the moment, if the snow lines (H₂O and CO) in Lupus are located around ~5 au and ~15 au, this could give rise to a bimodal population of planets and gaps if planet formation is enhanced at these snow lines.

Alternatively, though less likely, we could be observing a different regime where the smallest gaps in our sample are caused by internal photoevaporation rather than planet formation (Owen & Clarke 2012). Disks with cavities under 10 au, such as Sz108B, J16092697-3836269, Sz72, Sz90, Sz96, Sz131, and Sz76, have accretion rates of log₁₀[M_⊙· yr⁻¹] = −9.5, −8.1, −8.6, −8.9, −9.3, −9.1, and −9.1, respectively. Most of these align with the photoevaporation models by Ercolano et al. (2018) and Picogna et al. (2019) and could potentially be explained by this mechanism. However, J16092697-3836269 and Sz72, with larger accretion rates, fall outside the regime predicted by current photoevaporation models, suggesting that other mechanisms, like planets, may be causing these substructures. Generally it is difficult to prove whether cavities are caused by planets, photoevaporation, or a combination of the two (Gárate et al. 2023) without direct observations or limits on the planets itself. Finally, other physical mechanisms, such as dead zones, can also explain the formation of small cavities without invoking planets (Flock et al. 2015; Pinilla et al. 2016; Gárate et al. 2021).

For the remaining planets predicted by DBNets, Fig. 10 shows that they occupy a region where only a few exoplanets have been discovered (between 14–117 au and 20–2008 M_⊕). This could be explained by three possibilities: (1) it is challenging to detect exoplanets in this region due to observation limitations that come with low planet masses or large distances from the star; (2) planets may have formed in these regions but have since migrated inward, closer to their stars, where we now observe them (Kley & Nelson 2012; Lodato et al. 2019); or (3) the gaps observed in protoplanetary disks may not be the result of planet-disk interactions, indicating other formation processes at work, for example snow lines (Zhang et al. 2015) or zonal flows (Johansen et al. 2009).

In Fig. 11, there appears to be a relation between gap widths and the ring peak positions beyond each gap, with gaps and cavities following distinct trends. Additionally, this figure also reveals a lack of substructures (ring peaks) between ~13–27 AU, which is closely related to the cavity gap mentioned before in Fig. 10. The absence of substructures, either cavities or rings, in this region of the disks appears to persist even when considering a more direct observable. There does not appear to be a straightforward relationship between the disk dust mass and the mass of the planet carving the gap. However, the derived planet masses are much larger than the dust masses, even though above 10 M_⊕ most of the mass is probably in the gas envelope, this could imply that the millimeter emission we observe is coming from left-over dust particles that have not turned into planetesimals yet, meaning that the Lupus disks could have already finished most of their planet formation. On the other hand, the observed millimeter dust emission could be just a tracer of an unseen, larger population of solids.

Assuming that for other star-forming regions compact disks around M dwarfs are also common (Ansdell et al. 2016; Mulders et al. 2018), understanding the evolutionary paths of these types of disks is crucial for planet formation and dust evolution. The fact that we observe numerous compact disks in the Lupus region, both with and without substructures, aligns with the scenario from van der Marel & Mulders (2021), in which such disks, undergoing significant radial drift, can supply enough dust material to form multiple super-Earths (Sanchez et al. 2024). This offers a potential explanation for the origin of the exoplanet populations observed around M stars. In addition, the fact that we are still observing these disks indicates that another mechanism is hindering dust drift and trapping particles. Instead of forming a single Jupiter or Saturn-mass planet, several super-Earths or smaller planets may be forming in the inner regions of these disks, collectively stopping the drift (Huang et al. 2024b).

5 Summary and conclusions

We have analyzed a complete sample of protoplanetary disks in Lupus using the highest-resolution and highest-sensitivity data to date. The key findings of this work are as follows:

• We gathered high-resolution images of Lupus disks and used new observations of 33 faint disks to complete the Lupus sample of Class II disks at high resolution. Our findings reveal that over 67% of the disks in Lupus have dust radii smaller than approximately 30 au. Additionally, we discovered 11 new disk cavities, including one of the smallest cavities measured to date, with a radius of 4.1 au;

• Through visibility modeling, we measured the dust disk radii of several compact disks, finding sizes as small as 0.6 au;

• The observed properties of the compact disks are consistent with planet formation models which predict efficient formation of close-in super-Earths in drift-dominated disks around M dwarfs (Sanchez et al. 2024);

• We revisited the SLR, finding good agreement with previous works, which found $F_{mm}\propto R_{eff}^2$. This is consistent with drift-dominated disks. We extended the SLR by including the smallest disk sizes;

• By comparing the total flux and radius for each disk with radiative transfer models tailored to the stellar luminosity of each source, we estimated the dust mass and average optical depths of the compact disks. Both are generally low, with median values of ~3.38 M_⊕ and ~0.7, respectively. This implies that the observed dust masses in Lupus are not representative of the total solid mass budget for exoplanet formation and that the bulk of the pebbles may have already been converted into boulders or even planets or have already drifted inward and sublimated close to the host star;

• We estimated planet masses for each substructure in our Lupus sample using DBNets. The smaller cavities in our study align well with the population of RV-detected exoplanets, while the larger ones are consistent with previous findings and lie in a region of the parameter space where no exoplanets are currently detected. However, no clear correlation was observed between planet masses and disk dust masses in our sample.

Acknowledgements

We thank the anonymous referee for the insightful comments and suggestions, which helped improve the quality of this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2022.1.00154.S (PI: van der Marel, Nienke), ADS/JAO.ALMA#2018.1.01458.S (PI: Yen, Hsi-Wei), ADS/JAO.ALMA#2017.1.00388.S (PI: Liu, Hauyu Baobab), ADS/JAO.ALMA#2018.1.00689.S (PI: Muto, Takayuki) and ADS/JAO.ALMA#2022.1.01302.S (PI: Mulders, Gijs). ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. G.D.M. acknowledges support from FONDECYT project 11221206 and the ANID BASAL project FB210003.

Appendix A ALMA image characteristics

Appendix B GALARIO fitting results

Appendix C Lupus stellar parameters

Appendix D Visibility models

Fig. D.1 Visibility plots comparing the observed and modeled visibilities.

References

All Tables

All Figures

Fig. 1 All protoplanetary disks in Lupus at high resolution from various projects and observations with ALMA (see Section 2), shown at the same spatial scale. The disks are arranged in descending order based on their total flux. The scale bar and Neptune’s orbit in the final panel assumes that each disk is at a distance of 160 pc. Several cavities and substructures are observed, though the smallest disks are barely visible.

In the text

	Fig. 2 Zoomed-in view of the smallest disks in the Lupus molecular cloud. Some of these small disks exhibit substructure, although most appear featureless. The disks are arranged in descending order of their total flux, consistent with the disks shown previously in Fig. 1. The scale bar and Jupiter’s orbit in the final panel assumes that each disk is at a distance of 160 pc.
In the text

	Fig. 3 Best disk models obtained from the GALARIO visibility fitting. Note that the scale changes from Sz88 onward to better highlight the most compact disks. The disks are arranged in descending order based on their total flux derived from the visibility fitting.
In the text

Fig. 4 Size distribution of the Lupus disks. Disks exhibiting substructures in the continuum images are shown in dark blue, while smooth disks without visible substructures are depicted in red and added on top of the blue histogram. The pie charts inside the plot show the percentages of disks with and without substructures for sizes above and below 30 au. The vertical dashed orange line indicates our resolution limit, which is set at 0.03″ (4 au). More than 67% of the sample belongs to the compact disks classification.

In the text

Fig. 5 SLR for the Lupus star-forming region at high resolution. Disks with substructures are marked by purple circles outlined in dark blue, while smooth disks are represented by solid purple circles without outlines. Unresolved disks are indicated by downward-pointing triangles. Binary disks are marked with star symbols, and RU Lup is identified with a square marker for clarity. The drift-dominated slope, as described by Rosotti et al. (2019), is illustrated by the dashed red line, the trapdominated slope (Zormpas et al. 2022) in green, and the frag-dominated slope in black. Our fit, using the LINMIX package, is displayed in turquoise. All disks are normalized to a common distance of 160 pc for consistency.

In the text

Fig. 6 Comparison of the disk dust mass from a grid of 1728 models versus the mass derived from flux for the same models. In the top panel, stellar luminosity and effective temperature are fixed, illustrating how disk radius varies with flux and dust mass (see the colorbar with the radius values). The bottom panel shows the effect of changing stellar luminosity while keeping the disk radius and effective temperature constant. The dotted gray line highlights the masses from the flux that equal the disk dust mass introduced in the models.

In the text

Fig. 7 Comparison of the model grid with observed disk radii and total flux (normalized to 160 pc). The image represents the model dust mass on a log scale (also in brown contours), and the black contours show the mean optical depth, $\overline{\tau }$, defined in Equation (11), with values of 0.25 and 1. The dark blue edges indicate a visible substructure in this source and a star-shaped marker distinguishes binary disks.

In the text

Fig. 8 Top panel: masses of the compact disks in Lupus (R < 30 au). These masses were obtained through the RADMC-3D model grid for τ < 5. The dashed black line indicates where the RADMC-3D dust mass equals the optically thin dust mass calculation (see Eq. 10). Top-right panel: Kernel density estimation of the disk dust masses. Bottom panel: cumulative distribution functions (CDFs) of all the disk masses in Lupus, using two approaches: in teal, all masses are estimated with the optically thin approximation based on flux values (Eq. (10)), and in pink, we substitute the masses with values derived from our model grid where available.

In the text

	Fig. 9 Optical depths of the Lupus disks using a kernel density estimation for points with radii <30 au and obtained through Eq. (11) and the grid measurements.
In the text

Fig. 10 Lupus exoplanet masses obtained using DBNets plotted against the semimajor axis, assuming that the gaps are carved by planets. In blue, we show all observed exoplanets around M and K stars from the NASA Exoplanet Archive catalog, while the orange markers represent estimates for the Lupus sample. Brown outlines indicate cavities as the type of substructure, while green outlines denote gaps and rings. Additionally, planet mass estimates from Zhang et al. (2023), Huang et al. (2024a), and Ruzza et al. (2024) are highlighted in crimson circle markers.

In the text

	Fig. 11 Top panel: gap widths and locations of the ring peaks beyond each gap in our disks, as well as the DSHARP disks (Zhang et al. 2018) in Lupus. Gaps are shown in teal, and cavities in pink. Bottom panel: disk dust masses versus planet masses from the DBNets analysis, using the same colors to denote the substructure type.
In the text

	Fig. D.1 Visibility plots comparing the observed and modeled visibilities.
In the text