Molecules with ALMA at Planet-forming Scales
(MAPS). XIV. Revealing Disk Substructures in
Multiwavelength Continuum Emission

Anibal Sierra, Laura M. Pérez, Ke Zhang, Charles J. Law, Viviana V.

Guzmán, Chunhua Qi, Arthur D. Bosman, Karin I. Öberg, Sean M. Andrews,
Feng Long, et al.

► To cite this version:

Anibal Sierra, Laura M. Pérez, Ke Zhang, Charles J. Law, Viviana V. Guzmán, et al. Molecules with ALMA at Planet-forming Scales (MAPS). XIV. Revealing Disk Substructures in Multiwavelength Continuum Emission. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac1431. insu-03636625

HAL Id: insu-03636625

https://insu.hal.science/insu-03636625v1

Submitted on 19 Sep 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Molecules with ALMA at Planet-forming Scales (MAPS). XIV. Revealing Disk Substructures in Multiwavelength Continuum Emission

Anibal Sierra¹, Laura M. Pérez¹, Ke Zhang²,³,<fnref n="24" />, Charles J. Law⁴, Viviana V. Guzmán⁵, Chunhua Qi⁴, Arthur D. Bosman³, Karin I. Öberg⁴, Sean M. Andrews⁴, Feng Long⁴, Richard Teague⁴, Alice S. Booth⁶,⁷, Catherine Walsh⁷, David J. Wilner⁴, François Ménard⁸, Gianni Cataldi⁹,¹⁰, Ian Czekala¹¹,¹²,¹³,¹⁴,¹⁵,<fnref n="25" />, Jaehan Bae¹⁶,¹⁷,<fnref n="25" />, Jane Huang³,⁴,<fnref n="25" />, Jennifer B. Bergner¹⁸,<fnref n="25" />, John D. Ilee⁷, Myriam Benisty⁸,¹⁹, Romane Le Gal⁴,⁸,²⁰,²¹, Ryan A. Loomis²², Takashi Tsukagoshi⁹, Yao Liu²³, Yoshihide Yamato¹⁰, and Yuri Aikawa¹⁰

¹ Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile; asierra@das.uchile.cl

² Department of Astronomy, University of Wisconsin–Madison, 475 N. Charter St., Madison, WI 53706, USA

³ Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA

⁴ Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

⁵ Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile

⁶ Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands

⁷ School of Physics & Astronomy, University of Leeds, Leeds LS2 9JT, UK

⁸ Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁹ National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan

¹⁰ Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

¹¹ Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

¹² Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

¹³ Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

¹⁴ Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA

¹⁵ Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA

¹⁶ Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA

¹⁷ Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

¹⁸ University of Chicago Department of the Geophysical Sciences, Chicago, IL 60637, USA

¹⁹ Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UMI 3386, Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

²⁰ IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400 Toulouse, France

²¹ IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France

²² National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

²³ Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing 210023, People's Republic of China
*Received 2021 February 16; revised 2021 June 1; accepted 2021 June 2; published 2021 November 3*

## Abstract

Constraining dust properties of planet-forming disks via high-angular-resolution observations is fundamental to understanding how solids are trapped in substructures and how dust growth may be favored or accelerated therein. We use ALMA dust continuum observations of the Molecules with ALMA at Planet-forming Scales (MAPS) disks and explore a large parameter space to constrain the radial distribution of solid mass and maximum grain size in each disk, including or excluding dust scattering. In the nonscattering model, the dust surface density and maximum grain size profiles decrease from the inner disks to the outer disks, with local maxima at the bright ring locations, as expected from dust trapping models. The inferred maximum grain sizes from the inner to outer disks decrease from 1 cm to 1 mm. For IM Lup, HD 163296, and MWC 480 in the scattering model, two solutions are compatible with their observed inner disk emission: one solution corresponding to a maximum grain size of a few millimeters (similar to the nonscattering model), and the other corresponding to a size of a few hundred micrometers. Based on the estimated Toomre parameter, only IM Lup—which shows a prominent spiral morphology in millimeter dust—is found to be gravitationally unstable. The estimated maximum Stokes number in all the disks lies between 0.01 and 0.3, and the estimated turbulence parameters in the rings of AS 209 and HD 163296 are close to the threshold where dust growth is limited by turbulent fragmentation. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

*Unified Astronomy Thesaurus concepts: Circumstellar dust (236); Interstellar scattering (854); Radiative transfer (1335); Radio continuum emission (1340); Protoplanetary disks (1300)*

## 1. Introduction

High-angular-resolution observations of Class II disks have revealed that most of the large and bright disks have dust morphologies such as gaps and rings (e.g., Andrews et al. 2016; Isella et al. 2016; Fedele et al. 2017; Avenhaus et al. 2018;

Huang et al. 2018a; Cieza et al. 2021), vortices (e.g., Brown et al. 2009; Andrews et al. 2009; Casassus et al. 2013; Pérez et al. 2014; Robert et al. 2020), and spiral arms (e.g., Mouillet et al. 2001; Muto et al. 2012; Wagner et al. 2015; Pérez et al. 2016; Garufi et al. 2018). While the origins of these substructures are still under debate (see discussion in Andrews 2020), they provide a plausible mechanism for avoiding fast radial migration of dust grains toward the central star

<footnote n="24">²⁴ NASA Hubble Fellow.</footnote>

<footnote n="25">²⁵ NASA Hubble Fellowship Program Sagan Fellow.</footnote>
(Weidenschilling 1977). They also prevent the depletion of the raw dust material that can form planets. In addition, these structures may also facilitate planet formation as they concentrate dust grains into narrow spatial regions, which allows for more efficient growth.

Hydrodynamical simulations of protoplanetary disks with rings and gaps (e.g., Fouchet et al. 2007, 2010; Dipierro et al. 2015b), vortices (e.g., Richard et al. 2013; Surville & Barge 2015), spiral arms (e.g., Bae & Zhu 2018), and analytical models (e.g., Takeuchi & Lin 2002; Pinilla et al. 2012; Lyra & Lin 2013; Birnstiel et al. 2013) have shown that gas pressure maxima might act as dust traps, and they could be responsible for the observed substructures at millimeter wavelengths. In these models, the dust mass is enhanced at the pressure maxima and large grains are expected, due to dust size differential trapping (large grains of ~millimeter and centimeter sizes are more concentrated at the pressure maxima than small grains of ~100 µm; e.g., Ruge et al. 2016).

The width of the gas pressure bumps is an upper limit to the width of dust rings if the gas pressure maxima can effectively trap dust grains (e.g., Dullemond et al. 2018; Sierra et al. 2019). This scenario can be tested using high-angular-resolution images from molecular line emission and the millimeter dust continuum. The former is difficult, because the optically thin molecules that are needed to trace the underlying gas distributions at high spectral resolution are faint compared with the dust continuum emission. Furthermore, dust rings could be optically thick, complicating the extraction of a reliable gas column density profile. The width of the gas pressure bumps can also be inferred from the deviation of the gas angular velocity with respect to the Keplerian value (e.g., Rosotti et al. 2020). Even if the gas angular velocity profile can be estimated from observations, uncertainties in the mass of the central star, limitations in angular resolution, and the deprojection of the emission-line surface on the midplane can shift the position of the inferred pressure bump and the center of the dust continuum rings. Dust particles that are radially trapped in rings are also expected to be vertically settled; thus, the local scale height can be used (rather than the radial width of gas pressure bumps) to test dust trapping. Using this method, Dullemond et al. (2018) found that some of the rings in DSHARP disks (Andrews et al. 2018) are consistent with dust traps.

Dust trapping can also be tested using multiwavelength dust continuum observations. For example, the dust size differential trapping predicts narrower rings in the dust continuum emission at lower frequencies (e.g., Pinilla et al. 2015), or narrower distributions around vortex centers (e.g., Lyra & Lin 2013; Sierra et al. 2017). This behavior occurs because the grain size *a* traced at a wavelength *λ* follows *a* = λ/2π (e.g., D'Alessio et al. 2001). Thus, if the ring width or vortex area is measured at different frequencies, the presence of dust trapping can be determined. For instance, Cazzoletti et al. (2018) found a decrease in the vortex area with wavelength, which suggests that dust trapping is active in HD 135344B. Moreover, multiwavelength millimeter observations enable a direct determination of dust grain size by fitting the continuum spectrum. In Carrasco-González et al. (2019), the millimeter continuum spectrum of the HL Tau disk was fitted using ALMA and VLA observations, finding dust grains of millimeter sizes, which tended to be larger in the bright rings and smaller in the gaps.

Fitting the continuum spectrum also allows estimates of the dust temperature and surface density when sufficient frequency coverage is available. At least four wavelengths are needed to fit the dust surface density, maximum grain size, and dust temperature and compute their covariance matrix. However, such fitting can be performed with a smaller number of wavelengths if additional constraints on the physical properties are available, e.g., dust temperature. The surface density of solids computed from multiwavelength dust emission is fundamentally important to constrain the dust mass available in disks to form rocky planets. Previous dust mass estimations from many disks seem insufficient to explain the large incidence of massive exoplanets (Greaves & Rice 2010; Manara et al. 2018) and the core of giant planets, unless a large fraction of mass is already in large bodies that the millimeter observations are not fully tracing, as proposed by Najita & Kenyon (2014) in the disks of the Taurus-Auriga region.

Zhu et al. (2019) provide another possible solution to this problem. They found that optically thick disks can mimic optically thin disks if scattering is properly included in the radiative transfer equation, resulting in the underestimation of the dust mass by a factor from 3 to 30. The emission from an optically thick region can be a factor of 4 smaller than that without scattering if albedo is larger than ≳0.6 (Sierra & Lizano 2020). However, scattering can be neglected (even if albedo is high) if the emission is optically thin.

Over the past few decades, the albedo was assumed to be small at millimeter wavelengths, such that scattering effects could be neglected. However, Miyake & Nakagawa (1993) pointed out that scattering opacity could be dominating the dust opacity in disks at millimeter wavelengths. Since then, many models have included scattering using high albedo (e.g., D'Alessio et al. 2001). However, nonresolved observations from the optically thick inner disks were not able to observe the effects of scattering. Since ALMA, where the angular resolution was high enough to resolve the inner few tens of astronomical units in disks, many models have used scattering properties to model the inner disks (e.g., Soon et al. 2017; Liu 2019; Carrasco-González et al. 2019; Ueda et al. 2020).

In this paper, we determine the maximum grain size and dust surface density of five protoplanetary disks, by fitting spatially resolved multifrequency ALMA continuum observations profiles, and we investigate any evidence of dust trapping in the rings. We compute and compare the results from models with and without scattering in the radiative transfer equation.

This paper is organized as follows. Section 2 summarizes the dust continuum observations. Section 3 explains how we compute the dust properties, and the results are shown in Section 4. The optical depth, Toomre parameter, and maximum Stokes number are also estimated and presented. In Section 5, we discuss our results and how they depend on our assumptions and modeling, as well as highlight notable results for each disk. We present our conclusions in Section 6.

## 2. Observations

As part of the ALMA Large Program Molecules with ALMA at Planet-forming Scales (MAPS), we obtained continuum observations for the disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480 at Bands 3 and 6. Additional continuum ALMA observations in Bands 4 and 7 were compiled from archival data for AS 209 (Qi et al. 2019)
<table><caption>**Table 1**

Properties of the Millimeter Continuum Images</caption><thead><tr><th rowspan="2">Source</th><th rowspan="2">ALMA Band</th><th rowspan="2">ν (GHz)</th><th colspan="2">Synthesized Beam</th><th rowspan="2">rms Noise (μJy beam<sup>-1</sup>)</th><th rowspan="2">Peak Intensity (mJy beam<sup>-1</sup>)</th><th rowspan="2">Peak S/N</th><th rowspan="2">ALMA Project Code</th></tr><tr><th>(mas×mas; deg)</th></tr></thead><tbody><tr><td rowspan="4">IM Lup</td><td>3</td><td>94</td><td>243 × 194;</td><td>-81.5</td><td>14.1</td><td>6.44</td><td>457</td><td>2018.1.01055.L</td></tr><tr><td>3</td><td>106</td><td>227 × 182;</td><td>-77.1</td><td>14.1</td><td>6.90</td><td>489</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>226</td><td>116 × 81;</td><td>-80.7</td><td>19.3</td><td>15.0</td><td>777</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>257</td><td>91 × 84;</td><td>-87.2</td><td>19.2</td><td>16.1</td><td>838</td><td>2018.1.01055.L</td></tr><tr><td rowspan="6">GM Aur</td><td>3</td><td>94</td><td>341 × 209;</td><td>-28.0</td><td>14.3</td><td>2.15</td><td>150</td><td>2018.1.01055.L</td></tr><tr><td>3</td><td>106</td><td>295 × 210;</td><td>5.7</td><td>12.4</td><td>2.29</td><td>185</td><td>2018.1.01055.L</td></tr><tr><td>4</td><td>145</td><td>57 × 35;</td><td>-13.8</td><td>11.3</td><td>0.48</td><td>42</td><td>2017.1.01151.S</td></tr><tr><td>6</td><td>226</td><td>130 × 94;</td><td>9.1</td><td>16.7</td><td>5.83</td><td>349</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>257</td><td>117 × 83;</td><td>-1.5</td><td>21.9</td><td>5.74</td><td>262</td><td>2018.1.01055.L</td></tr><tr><td>7</td><td>283</td><td>270 × 162;</td><td>0.1</td><td>57.6</td><td>22.5</td><td>390</td><td>2015.1.00678.S</td></tr><tr><td rowspan="5">AS 209</td><td>3</td><td>94</td><td>263 × 189;</td><td>-78.3</td><td>11.6</td><td>3.94</td><td>339</td><td>2018.1.01055.L</td></tr><tr><td>3</td><td>106</td><td>253 × 193;</td><td>-74.2</td><td>13.2</td><td>3.86</td><td>294</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>226</td><td>109 × 77;</td><td>-86.7</td><td>15.6</td><td>4.91</td><td>314</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>257</td><td>99 × 78;</td><td>-83.1</td><td>21.9</td><td>6.12</td><td>279</td><td>2018.1.01055.L</td></tr><tr><td>7</td><td>283</td><td>262 × 180;</td><td>-73.8</td><td>54.24</td><td>32.1</td><td>592</td><td>2015.1.00678.S</td></tr><tr><td rowspan="4">HD 163296</td><td>3</td><td>94</td><td>251 × 184;</td><td>-88.3</td><td>11.2</td><td>10.1</td><td>902</td><td>2018.1.01055.L</td></tr><tr><td>3</td><td>106</td><td>248 × 184;</td><td>-88.1</td><td>13.2</td><td>10.7</td><td>811</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>226</td><td>110 × 80;</td><td>-80.8</td><td>17.7</td><td>14.2</td><td>802</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>257</td><td>119 × 76;</td><td>-81.5</td><td>20.3</td><td>17.6</td><td>867</td><td>2018.1.01055.L</td></tr><tr><td rowspan="4">MWC 480</td><td>3</td><td>94</td><td>325 × 221;</td><td>-24.5</td><td>14.6</td><td>11.3</td><td>774</td><td>2018.1.01055.L</td></tr><tr><td>3</td><td>106</td><td>295 × 212;</td><td>7.5</td><td>12.4</td><td>11.8</td><td>952</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>226</td><td>131 × 93;</td><td>7.5</td><td>16.9</td><td>24.5</td><td>1450</td><td>2018.1.01055.L</td></tr><tr><td>6</td><td>257</td><td>121 × 85;</td><td>-1.2</td><td>22.0</td><td>24.4</td><td>1110</td><td>2018.1.01055.L</td></tr></tbody></table>

and GM Aur (Huang et al. 2020). Details on the acquisition of these continuum observations, the calibration of interferometric visibilities, and the self-calibration procedure applied to these data and imaging strategies can be found in Öberg et al. (2021) and Czekala et al. (2021).

For all MAPS observations, we used the self-calibrated visibilities to create one image for each of the four correlator setups for Bands 3 and 6 described in Table 4 of Öberg et al. (2021). For each correlator setup, the continuum-only spectral window was combined with all other spectral windows where line emission was removed, by discarding channels away from the line center by ~15 km s⁻¹, and those near the disk wind in HD 163296 (for details see Booth et al. 2021). This results in four independent intraband continuum images, corresponding to a central frequency of 94, 106, 226, and 257 GHz or a wavelength of 3.20, 2.84, 1.33, and 1.17 mm, respectively. Additionally, the Band 4 and 7 continuum observations for AS 209 and GM Aur were also imaged starting from the self-calibrated visibilities, resulting in two independent continuum images: Bands 4 and 7, corresponding to a central frequency of 145 and 284 GHz or a wavelength of 2.07 and 1.06 mm, respectively. The continuum images we present here were generated using CASA version 6.1.0 (McMullin et al. 2007) and deconvolved with the TCLEAN task and Briggs weighting with a robust parameter of 0, which achieves a good balance between angular resolution and signal-to-noise ratio (S/N). For all images, the image scales were set to [0, 10, 20, 30, 60] pixels with a cell size of 0.01 pixel⁻¹. The units of the imaging residuals obtained from CASA are Jy {dirty beam}⁻¹, while the units of the convolved CLEAN models are Jy {CLEAN beam}⁻¹; thus, the final CLEANed images are computed after the beam units are corrected. This correction is calculated using the ratio between the volume of the dirty and CLEAN beams,

which spans ~0.7–1.1 for the dust continuum emission of the disks studied. This correction does not modify the total flux of the dust continuum emission because the residuals are small compared with the bright model image. A full description of this correction (which is called the “JvM correction” throughout the MAPS papers) can be found in Czekala et al. (2021), along with the general imaging strategies.

Table 1 summarizes all the relevant imaging properties for each target and each wavelength imaged. The rms noise and peak intensities are directly computed from each map. The rms noise is defined as the standard deviation of the intensity in a box far from the disk, while the peak intensity is the largest intensity within a region that encloses the disk. The peak S/N is the ratio between these two quantities.

The new images are presented in Figure 1, while the Band 4 and 7 images are consistent with those published in Qi et al. (2019) and Huang et al. (2020) and thus are not presented. Figure 1 shows the continuum brightness temperature at each observed wavelength for each disk. The brightness temperature was computed using the full Planck law instead of the Rayleigh-Jeans approximation. Four disks have concentric rings and gaps (AS 209, GM Aur, HD 163296, MWC 480), and one disk (IM Lup) has spiral arms and a gap, only clear in the Band 6 observations. The substructures of these disks in continuum emission are consistent with those observed at similar or higher angular resolution (Long et al. 2018; Huang et al. 2018a; Guzmán et al. 2018; Isella et al. 2018; Huang et al. 2020). A full analysis of the radial profiles from the continuum and line emission for each disk can be found in Law et al. (2021). In that work, the radial locations of all dust substructures are consistent with previous studies, and in addition, a new ring and gap are also reported in the continuum images: one in MWC 480 and another in IM Lup.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Brightness temperature of the continuum emission from the disks in IM Lup, GM Aur, AS 209, HD 163296, and MWC 480 (ordered by increasing stellar mass; see Table 1 in Öberg et al. 2021), from top to bottom, and four ALMA Bands at central frequencies of 94, 106, 226, and 257 GHz, from left to right. The beam is shown in the lower left corner of each panel. The color bars are stretched such that faint structures in the outer disk can be highlighted. Details about each image are described in Table 1.</figcaption></figure>

The Band 3 angular resolution and sensitivity are lower than at all other bands, and this is the main limitation of the multiwavelength analysis of these observations (Section 3). For this reason, all Band 3 visibilities are combined using the CONCAT task in CASA, and then we obtain a Band 3 image with a central frequency of 100 GHz ($\lambda$ = 3.0 mm) for further analysis. The Band 6 data sets are not merged because their angular resolution and sensitivity are higher than the Band 3 images, and they are used to fit the spectral energy distribution (SED) of the disks (Section 3).

To ensure the highest possible resolution from the data, the radial profiles were independently derived by two methods: first, we reimage all the data using a robust = -2 (uniform

weighting), in order to increase the angular resolution on the Band 3 image, which sets the limit on the smallest resolution element for this multiwavelength study. The RESTORINGBEAM task in CASA is used to obtain a circular beam. The radial profiles for Bands 4, 6, and 7 are then obtained from convolving all images to the Band 3 beam size using IMSMOOTH in CASA<fnref n="26" /> and averaging the emission in concentric ellipses with the inclination and position angle summarized in Table 1 in Öberg et al. (2021). The uncertainty of each intensity profile is given by $\Delta I_{\nu} = (\sigma_{I} + \text{rms}_{\nu})/\sqrt{n}$,

<footnote n="26">²⁶ IMSMOOTH takes into account the change in the units of the output image by rescaling it by the ratio of the input and output beams</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Azimuthally averaged intensity profiles at 100 GHz (red), 145 GHz (orange), 226 GHz (green), 257 GHz (blue), and 284 GHz (magenta). The shaded areas are the uncertainties ΔI<sub>ν</sub> of each profile, and the dashed lines are the absolute flux calibration uncertainty. Vertical dashed lines mark the position of bright rings in each disk. The horizontal black line in the lower left corner is the beam size of each disk.</figcaption></figure>

where σ<sub>I</sub> is the standard deviation around the mean intensity, rms<sub>ν</sub> is the rms noise in each map, and n = Ω<sub>ring</sub>/Ω<sub>beam</sub> is the number of beams within a ring of the disk, where Ω<sub>ring</sub> and Ω<sub>beam</sub> are the solid angles of the ring and beam, respectively. Second, we confirm that employing such a robust parameter does not introduce false structures in the radial profiles, as they are consistent with those obtained by fitting the visibilities with a nonparametric model using the software FRANKENSTEIN (Jennings et al. 2020). The final image products at the common highest resolution obtained from CLEAN imaging and from FRANKENSTEIN are shown in Appendix A.

Figure 2 presents the azimuthally averaged intensity profiles (I<sub>ν</sub>), which were computed using 10 radial bins per resolution element. The solid lines are the azimuthally averaged intensity for different frequencies (see color legend), and the shaded area is the uncertainty ΔI<sub>ν</sub> of each profile. The dashed lines are the absolute flux calibration uncertainties associated with each band. We set nominal values of 10% for Band 7 and 5% for Bands 6, 4, and 3. Note that this error has a magnitude close to ΔI<sub>ν</sub> in the inner disks, while it can be neglected in the outer disks. The beam size is indicated in the lower left corner of each panel. The vertical dashed lines in each panel are the positions of bright rings identified in previous studies (see references in Section 5 and the radial structures found by Law et al. 2021) and that can also be distinguished in our data at this resolution.

Figure 3 shows the spectral indices (α = d log(I<sub>ν</sub>)/d log(ν)) of the MAPS observations for each disk between 226 GHz and 100 GHz and between 257 GHz and 100 GHz. The spectral indices were directly computed from the slope of the spectrum at each radial point in the intensity profiles. The shaded areas are the uncertainties in the spectral indices associated with ΔI<sub>ν</sub>, while the absolute flux calibration uncertainty is constant at all

radii and is represented by two vertical lines in the blank panel. In general, spectral indices increase with disk radius, from α ~ 2 in the inner disk (optically thick emission) to α ~ 4 in the outer disk (optically thin), with local minima at the position of the bright rings. This behavior is consistent with previous dust continuum observations, e.g., in the disk around the HL Tau disk (Carrasco-González et al. 2019).

### 3. Methodology

The dust thermal emission of a plane-parallel slab of the disk is given by

$$I_{\nu} = B_{\nu}(T_{d})[1 - \exp(-\tau_{0}(\nu/\nu_{0})^{\beta})], \quad (1)$$

where B<sub>ν</sub>(T<sub>d</sub>) is the Planck function, T<sub>d</sub> is the dust temperature, τ<sub>0</sub> is the optical depth of dust at frequency ν<sub>0</sub>, and β is the dust opacity spectral index. The temperature, optical depth, and opacity spectral index are parameters that can be fitted to model the multiwavelength dust continuum emission of each disk. However, the optical depth τ<sub>ν</sub> = Σ<sub>d</sub>κ<sub>ν</sub> at frequency ν depends on the dust surface density (Σ<sub>d</sub>) and the opacity spectrum κ<sub>ν</sub> = κ<sub>0</sub>(ν/ν<sub>0</sub>)<sup>β</sup>, where the opacity coefficient κ<sub>0</sub> and β depend on dust properties such as dust composition, dust size distribution, maximum grain size, or porosity. In this work we fix the dust composition and the particle size distribution and only consider compact dust spheres (more details about the opacity properties below); thus, the magnitudes of the opacity coefficient and the opacity spectral index are only determined by the maximum grain size (a<sub>max</sub>).

Hence, the dust thermal emission at each radius can be fitted using the dust temperature, dust surface density, and maximum grain size. If scattering is taken into account, Equation (1) is no
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Spectral indices between 100 GHz and 226 GHz (blue curves), and between 100 GHz and 257 GHz (orange curves). The shaded areas are the uncertainties of each profile, while the vertical error bars represent the ALMA absolute flux calibration uncertainty. Vertical dashed lines mark the position of bright rings in each disk. The horizontal lines indicate $\alpha = 3.7$ (typical value in the ISM) and $\alpha = 2.0$ (optically thick emission in the Rayleigh–Jeans limit).</figcaption></figure>

longer valid. In that case, the solution to the radiative transfer equation for a vertically isothermal slab<fnref n="27" /> can be written as (Sierra et al. 2019)

$$I_{\nu} = B_{\nu}(T_{\mathrm{d}})[1 - \exp(-\tau_{\nu}/\mu) + \omega_{\nu}F(\tau_{\nu}, \omega_{\nu})], \quad (2)$$

where

$$F(\tau_{\nu}, \omega_{\nu}) = \frac{1}{\exp(-\sqrt{3}\epsilon_{\nu}\tau_{\nu})(\epsilon_{\nu} - 1) - (\epsilon_{\nu} + 1)} \times \left[ \frac{1 - \exp(-\frac{\sqrt{3}\epsilon_{\nu} + 1/\mu}{\tau_{\nu}})\tau_{\nu}}{\sqrt{3}\epsilon_{\nu}\mu + 1} + \frac{\exp(-\tau_{\nu}/\mu) - \exp(\sqrt{3}\epsilon_{\nu}\tau_{\nu})}{\sqrt{3}\epsilon_{\nu}\mu - 1} \right] \quad (3)$$

is the factor that modifies the emergent intensity from the typical nonscattering case, $\epsilon_{\nu} = \sqrt{1 - \omega_{\nu}}$, $\mu = \cos(i)$ is the cosine of the disk inclination, and $\tau_{\nu}$ is the total optical depth (absorption + scattering), which can be written in terms of the opacity coefficient and albedo ($\omega_{\nu}$) as $\tau_{\nu} = \sum_{d}\kappa_{\nu}/(1 - \omega_{\nu})$. In this work, the albedo is given by the effective albedo (Henyey & Greenstein 1941), which takes into account the nonisotropic scattering effects via the asymmetry parameter (Birstiel et al. 2018). The effective albedo also depends on the maximum grain size. Equation (2) reduces to Equation (1) when $\omega_{\nu} = 0$.

Equations (2)–(3) were found by direct integration of the analytical solution of Miyake & Nakagawa (1993) for a vertically isothermal slab, and a similar solution was derived by Birnstiel et al. (2018) using the Eddington approximation. Both

models predict changes in the spectral indices (with respect to the nonscattering case) when scattering is considered. Particularly, scattering can explain spectral indices smaller than the typical value of 2 in the optically thick regime (Liu 2019; Zhu et al. 2019; Sierra & Lizano 2020).

From Equation (2), we can constrain the dust properties of each disk by fitting their spatially resolved continuum spectrum from the Band 3–Band 7 intensity profiles, as was recently done for the HL Tau disk over a similar range of wavelengths (Carrasco-González et al. 2019). We will consider the cases with and without scattering. Even when the former is a more realistic approximation to the radiative transfer effects taking place in protoplanetary disks, there are many previous analyses in disks without scattering, from which we can compare our results. In order to avoid degeneracy between the dust temperature and dust surface density (see Appendix B), the dust temperature of each disk is fixed by the midplane temperature computed in Zhang et al. (2021), convolved to the angular resolution of each disk (see Appendix C). While we assume the same temperature profiles in the scattering and nonscattering cases, in reality they may differ if the albedo is high owing to differing amounts of disk heating (e.g., Dullemond & Natta 2003).

Additionally, throughout this work, we adopt the dust opacity properties given by Birnstiel et al. (2018) from the DSHARP collaboration. The particle size distribution is assumed to follow $n(a)da \propto a^{-p}da$, with a minimum grain size much smaller (~0.05 µm) than the maximum grain size (~mm, cm). The effects of the minimum grain size on the opacity properties can be neglected if $p < 4$ (Draine 2006). Typically, a value of $p = 3.5$ is assumed, which was found by Mathis et al. (1977) in the interstellar medium (ISM). However, the slope of the particle size distribution tends to decrease in

<footnote n="27">²⁷ The vertical isothermal slab is a good approximation, as most of the millimeter continuum emission is coming from a thin layer close to the midplane (where the dust scale height is much smaller than the gas scale height), with no abrupt changes in the vertical temperature distribution (e.g., Villenave et al. 2020).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left panel: optically thin spectral index $\alpha^{\text{thin}}$ as a function of the maximum grain size. The curves represent different slopes of the particle size distribution (see color legend). The dashed horizontal line is the reference spectral index in our disk sample at $\sim 100$ au. Right panel: maximum grain size as a function of the slope of the particle size distribution when the optically thin spectral index is $\alpha = 3$. The gray shaded area is the region where no solution is found.</figcaption></figure>

protoplanetary disks owing to dust growth (Drążkowska et al. 2019). Birnstiel et al. (2012) found that this slope changes from $p \sim 3.5$ when the maximum grain size is regulated by fragmentation to $p \sim 2.5$ when it is regulated by drift. The magnitude of $p$ determines the opacity properties, in particular, the optically thin spectral index. In our disk sample (Figure 3), the spectral indices are close to $\alpha \sim 3$ at $\sim 100$ au, where one expects optically thin emission.

The left panel of Figure 4 shows the optically thin spectral index $\alpha^{\text{thin}}$ as a function of the maximum grain size for different slopes $p$. The dashed horizontal line is the reference value of $\alpha = 3$ in the outer disks for our sample. We note that values of $p \gtrsim 3.4$ are not able to explain the reference spectral index with dust grains smaller than 1 m. To visualize this, the right panel of Figure 4 shows the maximum grain size, when the optically thin spectral index is 3, as a function of $p$. This particular maximum grain size rapidly increases for $p \gtrsim 3.2$. Then, as the magnitude of $p$ cannot be easily constrained, we assume a value of $p = 2.5$, which gives a lower limit for the maximum grain size in the disk.

The set of parameters (dust surface density and maximum grain size) that best describes the observations is computed as follows. At each disk radius, we define the probability of observing the intensities $I_{\nu_1}, I_{\nu_2}, \dots, I_{\nu_n}$ given an intensity model with parameters $a_{\text{max}}, \Sigma_{\text{d}}$ by the likelihood function

$$p(I_{\nu_1}, I_{\nu_2}, \dots, I_{\nu_n} | a_{\text{max}}, \Sigma_{\text{d}}) \propto \exp(-\chi^2/2), \quad (4)$$

where

$$\chi^2 = \sum_n w_{\nu_n} \times \left( \frac{I_{\nu_n} - I_{\nu_n}^{\text{model}}}{\epsilon_{\nu_n}} \right)^2 \quad (5)$$

is the chi-square statistic and $\epsilon_{\nu_n}$ is the total uncertainty of the intensity profiles, given by $\epsilon_{\nu_n}^2 = \Delta I_{\nu_n}^2 + (\delta I_{\nu_n})^2$, where $\Delta I_{\nu_n}$ is the uncertainty from the intensity profiles (shaved vertical region in Figure 2) and $\delta$ is the coefficient associated with the uncertainty of the flux calibration, which were described in

Section 2. The relative weight of each wavelength is given by $w_{\nu_n}$. We choose a relative weight of 2 for Band 3 because two independent data sets were merged to obtain this profile. All the other individual bands have a relative weight of 1.

The explored parameter space is a logarithmically spaced grid between $10 \ \mu\text{m} \le a_{\text{max}} \le 3 \ \text{cm}$ and $0.001 \ \text{g cm}^{-2} < \Sigma_{\text{d}} < 10 \ \text{g cm}^{-2}$ (these are sensible ranges in disks), each divided into 200 bins. The best values for the maximum grain size ($a_{\text{best}}$) and the dust surface density ($\Sigma_{\text{best}}$) are chosen as those that maximize the likelihood at each radius. In order to constrain the maximum grain size and dust surface density, we marginalize the joint probability $p(a_{\text{max}}, \Sigma_{\text{d}})$ and obtain the probability distributions of these two parameters, $P_{\text{amax}}$ and $P_{\Sigma}$.

As an example of our procedure, Figure 5 shows the fitting process for HD 163296 at a radius of 85 au in the nonscattering case. Shaded areas in the left panel are the regions in the parameter space able to reproduce the observed intensity (within the uncertainty) at each wavelength (red, green, and blue for frequencies 257, 226, and 100 GHz, respectively). To fit all observations simultaneously, we need these regions to intersect. The color scale in the right panel is the normalized joint probability distribution $p(a_{\text{max}}, \Sigma_{\text{d}})$, and the white dashed line is the $1\sigma$ contour around the maximum probability. Marginal normalized probabilities, $P_{\text{amax}}$ and $P_{\Sigma}$, are shown on the top and right of this panel, with orange dashed lines marking the best value and gray dashed lines marking the $1\sigma$ range. Repeating this procedure at all sampled radii allows us to infer $a_{\text{max}}$ and $\Sigma_{\text{d}}$ for the disks studied.

As an example of our procedure for the model that includes scattering, Figure 6 shows the fitting process at 20 au for HD 163296 in the scattering model. In this case, two different solutions agree with the observations, with $a_{\text{max}}$ being either millimeter-sized grains or few-hundred-micron grain sizes. We note that no such degeneracy is present if scattering is neglected. In the scattering case, this degeneracy appears in the inner disks (the radial extent where this solution is valid for
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Fitting process for HD 163296 at a radius of 85 au in the nonscattering model. Left panel: shaded areas are the regions in parameter space that are able to reproduce the observed emission at each wavelength (see color legend in the left top panel). Right panel: normalized joint probability distribution (color scale); the white dashed line is the $1\sigma$ contour. Blue lines in the top and right subpanels are the marginal probabilities for $a_{\text{max}}$ and $\Sigma_d$. Orange dashed lines mark the best-fit values, and gray dashed lines mark the $1\sigma$ constraints.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Fitting process for HD 163296 at a radius of 20 au in the scattering model. Left panel: shaded areas are the regions in parameter space that are able to reproduce the observed emission at each wavelength (see color legend in the lower right corner). Right panel: normalized joint probability distribution (color scale). Blue lines in the top and right subpanels are the marginal probabilities for $a_{\text{max}}$ and $\Sigma_d$.</figcaption></figure>

each disk is discussed below) because scattering only modifies the spectral indices when the absorption opacity is optically thick (e.g., Zhu et al. 2019; Sierra & Lizano 2020). In Figure 6, the left panel shows the regions that are consistent with the emergent intensities of each wavelength within the

uncertainties. Dust grains close to 100 µm match with the three wavelengths. Additionally, there is a region of large grains (>$1$ mm) that can also fit the data. The right panel is the normalized probability, while the subpanels around this panel are the normalized marginal probabilities ($P_{\text{max}}$, $P_{\Sigma}$).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Probability distribution of the dust surface density for all disks in our sample with nonscattering properties. The white solid line is the best-fit dust surface density ($\Sigma_{best}$). The color scale is the marginal probability $P_{\Sigma}$. Dashed black lines are the $1\sigma$ uncertainties, and the horizontal white line in the lower left corner is the beam size of each disk. Vertical dashed lines mark the position of bright rings. The dust mass in each panel corresponds to the integrated dust surface density within a radius where the enclosed flux at Band 3 is 95% of the total flux at this band.</figcaption></figure>

The small-grain and large-grain solutions provide a good fit to the observations, and we are not able to distinguish which one is better. Observations at longer wavelengths (e.g., $\lambda = 7$ mm) are needed to dismiss one of the solutions. For that reason, in the scattering case we split the probabilities into two regimes: small grains ($a_{max} < 300$ µm) and large grains ($a_{max} > 300$ µm), and we present both results computing the $1\sigma$ constraints as in the nonscattering case.

## 4. Results

Using the radial profiles in Figure 2 and the methodology described in the previous section, the dust surface density and maximum grain size profiles for all the disks are computed and presented here. We constrain the radial dust properties for each disk from half of the beam size to the radial position where the enclosed flux at Band 3 (100 GHz) is 95% of the total flux at this frequency. This percentage is computed using the radial profiles from the CLEANed images in Figure 2. The influence of this outer radius on the inferred dust masses is discussed in Section 5.

The results are presented as follows: First, we focus the results in the nonscattering model (Section 4.1), where Equation (1) is used to fit the SED. Then, scattering effects are taken into account (Section 4.2), using the general solution in Equations (2)–(3). In both cases, the optical depths at all wavelengths are computed, and dynamical parameters such as the Toomre parameter and the Stokes number of the maximum grain size are estimated. As mentioned at the end of the previous section, in the scattering case we present two solutions in the inner disk, one of them resulting in millimeter grain sizes and the second one corresponding to some-hundred-

micrometer grain sizes. All the radial profiles of each disk are discussed in Section 5.

### 4.1. Nonscattering Model

Figures 7 and 8 show the probability distributions for the dust surface density and maximum grain size, respectively. In both figures, the white solid lines are the best values ($\Sigma_{best}$, $a_{best}$), while the color scales are the marginal probabilities $P_{\Sigma}$ and $P_{a_{max}}$. The vertical dashed lines mark the positions of the bright rings in each disk that we are able to distinguish in our data (Law et al. 2021). The spatial resolution (beam size) is shown in the lower left corner of each panel. The black dashed lines in Figures 7 and 8 are the $1\sigma$ levels. The total dust mass for each disk is shown in the upper right legend of each panel. The dust mass uncertainty is computed from the $1\sigma$ level of each curve. In all cases, the dust surface density and maximum grain size tend to decrease with the disk radius, with local maxima in the rings. The properties of each disk are discussed in Sections 5.2.1–5.2.5.

In Figure 8, the green dashed line is a deconvolved power-law fit to the maximum grain size, i.e., we find a power law that, after convolution, fits the maximum grain size radial profiles ($a_{max} = \text{Convolution}[a_{10}(\text{Radius}/10\,\text{au})^{-b}]$). The free parameters of this power law are the maximum grain size at 10 au ($a_{10}$) and the slope $b$ of the power law. These values and the dust masses of each disk are summarized in Table 2.

The optical depths at different wavelengths are computed as $\tau_{\nu} = \Sigma_{d}\kappa_{\nu}$, where the dust opacity coefficient is computed adopting the inferred maximum grain size. Figure 9 shows the optical depths of the five disks at different wavelengths. Three of the disks (IM Lup, HD 163296, MWC 480) have an optical
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Probability distribution of the maximum grain size for all disks in our sample with nonscattering properties. The white solid line is the best-fit maximum grain size ($a_{best}$). The color scale is the marginal probability $P_{amax}$. Dashed black lines are the $1\sigma$ uncertainties, and the horizontal white line in the lower left corner is the beam size of each disk. The green dashed line is a power-law fit to the maximum grain size with radius; vertical dashed lines mark the position of the bright rings.</figcaption></figure>

<table><caption>**Table 2**
Inferred Dust Properties in the Nonscattering and Scattering Models with Large Grains ($a_{max} > 300~\mu\text{m}$)</caption><thead><tr><th rowspan="2">Disk</th><th colspan="3">Nonscattering Model</th><th colspan="3">Scattering Model</th></tr><tr><th>Dust Mass<br>(× 10<sup>-3</sup>M<sub>☉</sub>)</th><th>a<sub>10</sub><br>(cm)</th><th>b</th><th>Dust Mass<br>(× 10<sup>-3</sup>M<sub>☉</sub>)</th><th>a<sub>10</sub><br>(cm)</th><th>b</th></tr></thead><tbody><tr><td>IM Lup</td><td>1.92<sup>+1.91</sup><sub>-0.46</sub></td><td>1.03 ± 0.88</td><td>0.99 ± 0.04</td><td>3.64<sup>+1.99</sup><sub>-1.41</sub></td><td>1.15 ± 0.08</td><td>1.06 ± 0.04</td></tr><tr><td>GM Aur</td><td>0.74<sup>+0.34</sup><sub>-0.13</sub></td><td>0.82 ± 0.05</td><td>0.70 ± 0.03</td><td>0.73<sup>+0.32</sup><sub>-0.13</sub></td><td>0.82 ± 0.05</td><td>0.70 ± 0.03</td></tr><tr><td>AS 209</td><td>0.91<sup>+0.76</sup><sub>-0.27</sub></td><td>0.83 ± 0.05</td><td>0.47 ± 0.03</td><td>0.75<sup>+0.72</sup><sub>-0.21</sub></td><td>0.79 ± 0.04</td><td>0.47 ± 0.03</td></tr><tr><td>HD 163296</td><td>0.90<sup>+0.79</sup><sub>-0.22</sub></td><td>0.73 ± 0.05</td><td>0.56 ± 0.04</td><td>0.83<sup>+1.35</sup><sub>-0.18</sub></td><td>0.57 ± 0.03</td><td>0.46 ± 0.03</td></tr><tr><td>MWC 480</td><td>0.97<sup>+0.16</sup><sub>-0.18</sub></td><td>0.49 ± 0.02</td><td>0.41 ± 0.03</td><td>1.19<sup>+1.77</sup><sub>-0.29</sub></td><td>0.31 ± 0.01</td><td>0.21 ± 0.03</td></tr></tbody></table>

depth close to or above 1 at Band 6 in the inner disks (≈20 au), but GM Aur and AS 209 have an optical depth <1 everywhere.

Additionally, the Toomre parameter ($Q$) and the Stokes number ($St(a_{max})$) of the maximum grain size can be estimated from the above results. The former quantifies the gravitational disk stability, in particular, if $1 < Q < 1.7$, the disk is stable to linear perturbations, but second-order perturbations could grow and make the disk gravitationally unstable (Toomre 1964). If $Q < 1$, the disk is also unstable to linear perturbations.

The Stokes number describes the dust dynamics, in particular, it quantifies the coupling between dust and gas (Whipple & Planet 1972). Dust grains with a large Stokes number ($St \gg 1$) are decoupled from the gas, i.e., they are not strongly affected by the gas dynamics, while dust grains with a small Stokes number ($St \ll 1$) are coupled with the gas. These parameters are defined as

$$Q = \frac{c_s \Omega}{\pi G \Sigma_g}, \qquad (6)$$

$$St = \frac{\pi \rho_m a}{2 \Sigma_g}, \qquad (7)$$

where $c_s$ is the sound speed; $\Omega$ is the disk angular rotation frequency (assumed to be Keplerian); $G$ is the gravitational constant; $\Sigma_g$ is the gas surface density; $\rho_m = 1.675 \text{ g cm}^{-3}$ (Birnstiel et al. 2018) is the dust bulk density, computed from the dust composition mix of solids (water ice, silicates, troilite, and refractory organics); and $a$ is the grain size. In particular, we compute the Stokes number associated with the maximum grain size $St(a_{max})$. We subsequently refer to this parameter as the maximum Stokes number. The sound speed is computed from the midplane temperature as

$$c_s^2 = \frac{k_B T}{\mu_m m_H}, \qquad (8)$$

with a mean molecular weight $\mu_m = 2.34$, and where $k_B$ and $m_H$ are the Boltzmann constant and the hydrogen mass, respectively.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Dust optical depths at 100 GHz (red), 145 GHz (orange), 226 GHz (green), 257 GHz (blue), and 284 GHz (magenta) for the five disks studied in the nonscattering model. The shaded areas are the uncertainties of each profile. The horizontal black line in the lower left corner is the beam size of each disk. Vertical dashed lines mark the position of the bright rings, and the horizontal dotted line is a reference value where the optical depth is 1.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Toomre Q parameter (left) and maximum Stokes number (right) for the five disks (see legend on the left panel) assuming a dust-to-gas mass ratio of 1/100. The horizontal dashed lines in the left panel ($Q<1, 1.7$) delimit the regions where the disks become gravitationally unstable to linear and second-order perturbations, respectively. The shadowed areas are the uncertainties of each curve, which were computed from the uncertainty of the dust surface density and maximum grain size.</figcaption></figure>

The gas surface density is assumed to be 100 times the dust surface density. The influences of these assumptions on the results are discussed in Section 5. Figure 10 shows the Toomre parameter (left panel) and the Stokes number (right panel) of the five disks (see color legend). Most disks appear to be gravitationally stable, except IM Lup, while the Stokes number is below 1 across the sample at all radii (the properties of each disk are discussed in Sections 5.2.1–5.2.5).

## 4.2. Scattering Model

We now consider the case where scattering is taken into account in the radiative transfer calculations (Equations (2)–(3)) and recompute the probability distributions. As mentioned in Section 3, the inner disks have two possible solutions, one

with large grains ($a_{max} > 300$ µm) and a second one with small grains ($a_{max} < 300$ µm).

Figures 11 and 12 show the dust surface density and maximum grain size, respectively, from the scattering case and large grains. The white solid lines are the best values ($\Sigma_{best}$, $a_{best}$), and the color scales are the marginal probabilities ($P_{\Sigma}$, $P_{amax}$). The dust masses and the free parameter in the power-law fit to $a_{best}$ are also summarized in Table 2.

Figure 13 shows the optical depths associated with the absorption coefficient ($\tau_{\nu}^{abs} = \Sigma_{d}\kappa_{\nu}$) and the total optical depth, where scattering opacity is also included ($\tau_{\nu}^{sca} = \Sigma_{d}\kappa_{\nu}/(1 - \omega_{\nu})$), for the five disks at different wavelengths. As the albedo is large, total opacities increase by one order of
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Probability distribution of the dust surface density for all disks in our sample when scattering is taken into account in the radiative transfer model and the large-grain solution ($a_{max} > 300$ µm) is considered for the inner disk. The white solid line is the best-fit dust surface density ($\Sigma_{best}$). The color scale is the marginal probability $P_{\Sigma}$. Dashed black lines are the $1\sigma$ uncertainties, and the horizontal white line in the lower left corner is the beam size of each disk. Vertical dashed lines mark the position of the bright rings. The dust mass in each panel corresponds to the integrated dust surface density within a radius where the enclosed flux at Band 3 is 95% of the total flux at this band.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Probability distribution of the maximum grain size for all disks in our sample when scattering is taken into account in the radiative transfer model and the large-grain solution ($a_{max} > 300$ µm) is considered for the inner disk. The white solid line is the best-fit maximum grain size ($a_{best}$). The color scale is the marginal probability distribution $P_{a_{max}}$. Dashed black lines are the $1\sigma$ uncertainties, and the horizontal white line in the lower left corner is the beam size of each disk. The green dashed line is a power-law fit to the maximum grain size with radius; vertical dashed lines mark the position of the bright rings.</figcaption></figure>

magnitude compared to the nonscattering opacities. The Toomre Q parameter and the maximum Stokes number for the scattering case are quite similar to those derived when

scattering is neglected (Figure 10), since the maximum grain size and dust density have similar behaviors and magnitudes; thus, they are not presented again.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Dust optical depths at 100 GHz (red), 145 GHz (orange), 226 GHz (green), 257 GHz (blue), and 284 GHz (magenta) of the five disks (see legend in the upper right corner of each panel) in the scattering models and large grains ($a_{max} > 300$ µm). The colored dashed and solid lines are the optical depths associated with absorption and scattering + absorption, respectively. The shaded areas are the uncertainties of each profile. The horizontal black line in the lower left corner is the beam size of each disk. Vertical dashed lines mark the position of the bright rings, and the horizontal dotted line is a reference value where the optical depth is 1.</figcaption></figure>

On the other hand, the solution for small grains ($a_{max} < 300$ µm) is particularly interesting because if the maximum grain size is a few hundred micrometers, the disks could be consistent with the polarization pattern observed in many disks (Kataoka et al. 2015).

Figure 14 shows the small-grain solution for three of the disks, as GM Aur and AS 209 have absorption optical depth smaller than 1 at all radii (Figure 13) and thus are not consistent with this small-grain solution. We emphasize that scattering has important effects on the emergent intensity only when the optical depth associated with the absorption component is thick (Zhu et al. 2019; Sierra & Lizano 2020). The top panels of Figure 14 are the dust surface densities, the middle panels are the maximum grain sizes, and the bottom panels are the optical depths considering scattering and absorption opacity. The radial extent of each disk is shown according to the region where a solution with some-hundred-micrometer grain sizes is found. This small-grain solution requires a large amount of solids. This occurs because the dust opacity of ~100 µm grains at millimeter wavelengths is smaller than the dust opacity of millimeter grains. For example, at $\lambda = 1.3$ mm, the dust opacity of 1 mm grains is a factor of ~4.5 larger than the opacity of 100 µm grains. Then, the small grains need to have a larger amount of mass in order to reach the same intensity of millimeter grain sizes. The dust mass (only within the radius shown in Figure 14) is shown in the lower left corner of each panel, which is a factor of 2–5 larger than the whole mass associated with millimeter grains (Table 1).

Figure 15 shows the Toomre parameter (left panel) and the maximum Stokes number (right panel) for the small-grain solution and the three disks. The Toomre parameter and the maximum Stokes number are small compared with those in Figure 10, since the maximum grain sizes are smaller and thus

the dust surface densities are higher compared with the results from the large-grain solution.

## 5. Discussion

Radial profiles of the dust surface density, maximum grain size, optical depths, Toomre *Q* parameter, and maximum Stokes number of the five disks were computed in Section 4. Here we discuss the results from the nonscattering (Section 4.1) and scattering (Section 4.2) models, discuss the main properties of each disk, and compare our results with constraints from previous works. However, before discussing our results, it is important to address how they depend on the assumptions and modeling.

First, all the inferred properties depend on angular resolution, which in our sample varies from 19 to 39 au. The contrasts between rings and gaps are expected to be higher than those presented in this work if the rings and gaps are not resolved at all wavelengths. Beam averaging from optically thick regions (with low spectral index) and optically thin regions (with large spectral index) could lead to overestimation of the spectral indices at the location of rings and simultaneously lead to underestimation of $\alpha$ in gaps (e.g., see an example and discussion in Lin et al. 2020). Multiwavelength observations that could spatially resolve the dust rings at all wavelengths could better resolve these different regions and better infer their dust properties.

The dust surface density and maximum grain size are the main properties from which all other parameters are computed. The maximum grain size is assumed to be constant above the midplane, i.e., we do not consider settling. The magnitude of the maximum grain size depends on the assumed slope of the particle size distribution ($p$). In particular, we choose $p = 2.5$,
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Probability distribution of the dust surface density (top panels), maximum grain size (middle panels), and optical depths (bottom panels) for three disks in our sample when scattering is taken into account in the radiative transfer model and the small-grain solution ($a_{max} < 300$ µm) is considered for the inner disk. The color scales in the top and middle panels are the marginal probabilities for each parameter, and the black dashed lines are the $1\sigma$ uncertainties. The dust mass (only within the radial extent where a solution is found) of each disk is shown in the legend of the top panels. The optical depths in the bottom panels consider scattering and absorption. The horizontal dashed line in the bottom panels is a reference value where the optical depth is 1. The shaded areas are the uncertainties of each profile. Different colors correspond to different frequencies (see colored legends).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Toomre Q parameter (left) and maximum Stokes number (right) for three of the disks (see legend on the left panel) in the scattering model with small grains ($a_{max} < 300$ µm). The horizontal dashed lines in the left panels ($Q < 1, 1.7$) delimit the regions where the disk becomes gravitationally unstable to linear and second-order perturbations, respectively. The shaded areas are the uncertainties of each curve.</figcaption></figure>

which is a typical value for the slope when the maximum grain size is limited by radial drift (Birnstiel et al. 2012) and gives a lower limit to the inferred maximum grain size, as we show in Section 3. Recently, Macías et al. (2021) fitted the dust properties of the disk around TW Hya. In their modeling, the value of *p* was initially fixed and then considered as a free parameter. In the former case, the maximum grain size distribution has local maxima and minima that correlate with

the dust rings and gaps. In the latter case, most of the radial substructure of the maximum grain size disappears, but the radial profile for *p* exhibits local minima and maxima at rings and gaps, respectively. The maximum grain size and the slope *p* are highly degenerate. A more realistic model (e.g., dynamical dust simulations) should be able to account for changes in *p* in gaps and rings. However, this is beyond the scope of this paper.
On the other hand, as discussed in Section 3, the dust temperature model was fixed for each disk using the midplane temperature from Zhang et al. (2021). This avoids degeneracy, especially in the optically thin regime, where the emergent intensity scales as $I_{\nu} \propto \Sigma_{d}T_{d}$, and one can only constrain the product between the temperature and dust surface density. However, uncertainties in the emergent intensities are propagated to the dust surface density and maximum grain size, and in some cases the observed multiwavelength emission can be reproduced using a large range of dust surface densities and maximum grain sizes, even if these two are the only free parameters. More details about this degeneracy are described in Appendix B.

The optical depths are directly computed from the dust surface density and from the maximum grain size, which sets the opacity coefficient and albedo in the scattering case. Thus, there are no additional assumptions when computing the optical depth of the dust continuum emission.

However, when computing the Toomre Q parameter and the maximum Stokes number, additional assumptions about the gas properties are needed. First, the gas surface density is assumed to be 100 times the dust surface density. This factor comes from the typical value in the ISM. However, Ansdell et al. (2016) computed the dust mass from 890 µm continuum observations and the gas mass from CO isotopologue in a survey of disks in Lupus and found that most disks in this region have a dust-to-gas mass ratio ($\epsilon$) larger than 1/100, even when specific examples with smaller dust-to-gas mass ratio have also been inferred (e.g., in the TW Hydra disk; Zhang et al. 2017).

Additionally, the disk sizes obtained from the millimeter dust continuum emission (see a complete analysis of the MAPS disks sizes in Law et al. 2021) tend to be smaller than that obtained from the gas molecular line emission (e.g., Piétu et al. 2005; Isella et al. 2007). Then, it is expected that most of the continuum emission at these wavelengths comes from regions where the dust-to-gas mass ratio has been enhanced with large grains by radial drift.

Numerical simulations of dust trapping in rings (e.g., Gonzalez et al. 2017; Pinilla et al. 2020) or in spiral arms (e.g., Dipierro et al. 2015a) have also found changes in the dust-to-gas mass ratio as dust grains migrate toward pressure maxima. However, the magnitude of $\epsilon$ depends on factors such as the amount of solids in the simulation, turbulence, fragmentation velocity, the evolutionary stage (all poorly constrained), or whether the dust back reaction is taken into account or not. In these works, $\epsilon$ varies from ~0.003 to ~1.

As the Toomre parameter and Stokes number are inversely proportional to the gas surface density ($Q$, St $\propto \Sigma_{g}^{-1} = \epsilon \Sigma_{d}^{-1}$), they are only vertically shifted with respect to the results in Figure 10 (where we assume $\epsilon = 1/100$) if the dust-to-gas mass ratio is modified. For example, in Dullemond et al. (2018), they assume that the upper limit of the gas surface density is the dust surface density (i.e., $\epsilon = 1$). In that upper limit, the Stokes numbers in Figure 10 lie within 1 and 10, and all the disks would be gravitationally stable. We also compute the Stokes number using the gas surface density derived from the CO column densities by Zhang et al. (2021), and we find that the Stokes numbers for the five disks lie between $10^{-2}$ and $10^{-1}$, as in Figure 10.

Thus, scaling the dust mass by a factor of 100 is only a first approximation to the gas mass. Many physical processes can

<table><caption>**Table 3**
Percentage of the Enclosed Flux within the Fitting Radius</caption><thead><tr><th rowspan="2">Frequency<br>(GHz)</th><th colspan="5">Percentage</th></tr><tr><th>IM Lup</th><th>GM Aur</th><th>AS 209</th><th>HD<br>163296</th><th>MWC 480</th></tr></thead><tbody><tr><td>100</td><td>95%</td><td>95%</td><td>95%</td><td>95%</td><td>95%</td></tr><tr><td>145</td><td>...</td><td>93%</td><td>...</td><td>...</td><td>...</td></tr><tr><td>226</td><td>78%</td><td>91%</td><td>91%</td><td>84%</td><td>88%</td></tr><tr><td>257</td><td>78%</td><td>91%</td><td>93%</td><td>85%</td><td>88%</td></tr><tr><td>284</td><td>...</td><td>91%</td><td>93%</td><td>...</td><td>...</td></tr><tr><td>Fitting radius</td><td>160 au</td><td>188 au</td><td>136 au</td><td>122 au</td><td>113 au</td></tr></tbody></table>

**Note.** The fitting radius is defined as the radius where the enclosed flux at Band 3 is 95% of the total flux. These percentages were computed using the radial profiles from the CLEANed images at the same angular resolution (Figure 2).

change this value and modify the Toomre parameter and the maximum Stokes number in Figure 10. In addition, the factor of 100 is assumed constant thought the disk, i.e., the shape of the gas and dust surface density are the same; however, this approximation is valid for disks with grains well coupled with the gas and/or a high turbulence (e.g., Birnstiel et al. 2013; Ruge et al. 2016), such that these solids can be used to trace both the gas and dust surface density. Smoother structures are expected in gas than in dust. On average, the widths of all gas features for the MAPS disks are much larger than the dust, and the chemical substructures have lower relative contrasts than the continuum substructures (Law et al. 2021). Then, one expects that the dust-to-gas mass ratio should be larger in rings (where dust grains are being trapped) than in gaps (where dust grains are depleted). Consequently, the assumption $\Sigma_{g} = \Sigma_{d}/\epsilon$ with a constant $\epsilon$ overestimates the gas surface density in rings and underestimates it in gaps. Thus, the Toomre parameter and maximum Stokes number shown in Figure 10 are overestimated in gaps and underestimated in rings.

On the other hand, the sound speed is computed from the dust temperature in the midplane, which is assumed to be the same as the gas temperature. However, the two components may not be in thermal equilibrium. Finally, we assume that the gas velocity is Keplerian, but the gas velocity is sub-Keplerian or super-Keplerian depending on the local pressure gradients (Takeuchi & Lin 2002). Even when a deviation from Keplerian angular velocity is smaller than 10% (e.g., Teague et al. 2018, 2019), it produces a sufficiently large observational signature. In fact, Rosotti et al. (2020) were able to estimate the width of the gas rings from the velocity curves in a disk with deviations of only ~3% in Keplerian velocity. Additional uncertainty is introduced owing to a lack of precise knowledge in the stellar masses or the assumption of equal dust composition and structure throughout the disks.

The dust masses estimated in this work are directly computed from integrating the dust surface densities. As mentioned in Section 3, the dust surface density is computed inside the radius where the enclosed flux at Band 3 is 95% of the total flux. However, the continuum emission at higher frequencies is more extended than this radius, but the mass beyond this radius is not considered in the total mass estimate. Table 3 shows the percentage of the total flux at different frequencies within the radius that encompasses 95% of the total flux in Band 3. At higher frequencies, the enclosed flux for GM Aur and AS 209 is larger than 90%, while for HD 163296 and
MWC 480 it is larger than 85%. For IM Lup, the most massive disk, the enclosed flux fraction is only 78%. Even when the dust masses inferred in this work could be interpreted as a lower limit to the total dust mass, the inferred dust surface densities decrease by more than one order of magnitude from the inner to the outer disk, and thus a nonsignificant contribution to the total mass is expected beyond the fitting radius.

In general, the dust masses agree or are within a factor of a few compared with previous estimates (each disk is discussed below). The reasons for the disagreement (in addition to the fitting radius) are diverse, but the main are the following: Typically, the dust mass is estimated from the millimeter flux as

$$M_{\text{dust}} = \frac{F_{\nu} d^2}{\kappa_{\nu} B_{\nu}(T_d)}, \qquad (9)$$

where $F_ν$ is the flux at some frequency $ν$, $d$ is the distance to the source, $κ_ν$ is the opacity coefficient at the same frequency, and $B_ν$ is the Planck function at temperature $T_d$. Equation (9) assumes that the disk is optically thin everywhere, which is not always valid in dense rings or the inner disk, especially at short millimeter wavelengths (e.g., at ALMA Band 7 or Band 6). The opacity coefficient depends on the dust composition assumed, but in general the dust opacities from Beckwith et al. (1990) or D’Alessio et al. (2001) have been most used in the literature. And recently, the DSHARP opacities (Birnstiel et al. 2018) have become the standard to compute opacity properties. In particular, the Beckwith et al. (1990) opacity prescription was found to be consistent with the dust opacities inferred at ALMA Bands 7, 6, and 3 in the disk around HH 212 (Lin et al. 2021). The absorption spectrum from Beckwith et al. (1990) is given by $κ_ν = 2.3(ν/230 \text{ GHz})^\beta \text{ cm}^2 \text{ g}^{-1}$, and a value of $β = 1$ is usually assumed to scale this opacity to other frequencies. For example, the opacity coefficient from Beckwith at 100 GHz is $1.0 \text{ cm}^2 \text{ g}^{-1}$. From the DSHARP opacities, the opacity coefficients for a grain size of 1 mm at 230 and 100 GHz are 1.89 and 0.31 $\text{cm}^2 \text{ g}^{-1}$, respectively. At 230 GHz, the opacity coefficient only differs by a factor of 1.2, but at 100 GHz, they differ by a factor of 3.2. The main reason for the disagreement is the assumed value of $β$, which depends on the grain size and can only be estimated given the spectral index between two wavelengths. In our modeling, the disks are not assumed to be optically thin, and the value of $β$ is adopted depending on the maximum grain size at each radius.

Thus, differences in dust composition, temperature, and nonresolved optically thick emission can be some of the reasons of disagreement between the inferred dust masses in our works and others.

## 5.1. Scattering versus Nonscattering

Scattering is expected to modify the spectral indices for optically thick emission (Liu 2019; Zhu et al. 2019; Sierra & Lizano 2020). In particular, spectral indices below 2 can be obtained in the optically thick regime if the albedo is high and increases with frequency. The latter condition is satisfied in ALMA Bands 7, 6, and 3 for spherical dust grains with sizes of about a hundred micrometers. For grains larger than 1 mm, scattering cannot reproduce spectral indices below ~2.1 (see

Figure 9 in Zhu et al. 2019 and Figure 3 in Sierra & Lizano 2020), while in the nonscattering case the optically thick spectral index is always 2.

These previous effects were computed in the Rayleigh-Jeans (RJ) regime. However, this regime is not valid for low dust temperatures. The RJ regime is only valid for $[ν/\text{GHz}] < < 20$ [$T/K$], which is not fully satisfied at Band 6 given the dust temperature models in our disk sample (Appendix C). When the RJ regime is not satisfied, the spectral indices (in both the scattering and nonscattering cases) decrease because the peak of the Planck function is displaced to the submillimeter range according to Wien’s law (e.g., see Figure A1 in Sierra & Lizano 2020).

Whether or not RJ is valid, the change of the optically thick spectral index with maximum grain size in the scattering case is responsible for the two possible solutions found in the inner disks. In particular, the optically thick spectral indices for few-hundred-micrometer and millimeter grain sizes are similar.

The solutions for the nonscattering case and the scattering case with large grains ($a_{max} > 300 \text{ µm}$) are the same in the outer disk. This occurs because the grains are optically thin with respect to absorption at these radii, and the radiative transfer solution with scattering (Equation (2)) reduces to the nonscattering solution (Equation (1)) in this limit, even if the albedo is high. By contrast, the solutions in the inner disks do depend on whether scattering is included or not. The best models are dominated by millimeter grain sizes in the scattering and nonscattering case, and the dust surface densities are also similar. However, the marginal distributions are modified, due to changes of the spectral index with grain size in this regime. From Table 1, one can see that the dust masses are similar in the scattering and nonscattering models, except for IM Lup, where the dust mass in the scattering model is a factor of 2 larger than the dust mass in the nonscattering model. This occurs because if scattering is neglected, the dust mass from the nonscattering model is underestimated in the optically thick regime (Zhu et al. 2019), and IM Lup is the most optically thick disk in our sample in Band 6.

The small-grain regime ($a_{max} < 300 \text{ µm}$) in the scattering case provides an alternative solution that is able to reproduce the observations in the inner disks. No solution is found in the outer disks in this case because the individual regions that can explain each wavelength do not overlap when $a_{max} < 300 \text{ µm}$. The solution in the inner disks corresponds to grain sizes between 100 and 300 µm, and the dust surface densities are higher than the large-grain solution. This occurs because the opacity coefficient of hundred-micrometer grains is smaller than the opacity of millimeter grains by a factor of ~5, i.e., more mass is required for small grains that have a smaller opacity coefficient, such that they can have the same amount of emission of larger grains with a larger opacity coefficient.

High dust surface densities and some-hundred-micrometer sizes were already suggested by Ueda et al. (2020) and Macías et al. (2021) in the inner region of the TW Hya disk. They fitted the SED using a scattering and a nonscattering model and found that the disk is consistent with 300 µm grains and large densities ($Σ_d = 10 \text{ g cm}^{-2}$) if scattering is taken into account. However, if they neglect the scattering effects, the SED could also be explained by optically thick emission and large grains.

Some of the disks in our sample have been studied with polarized observations. For example, Hull et al. (2018) found that polarized properties of the inner disk of IM Lup are
consistent with 61 µm grain sizes, but millimeter to centimeter
grain sizes are needed to reproduce the unpolarized spectrum
(consistent with this work). The polarization fraction of the
inner disk (0.5) of AS 209 at 870 µm is only 0.2% (Mori et al.
2019), a very low value compared with the expected polarized
fraction if the maximum grain size is around some hundred
micrometers (Kataoka et al. 2015). This supports the fact that
we find that some-hundred-micrometer sizes are not consistent
with the SED in the inner disk, and millimeter grains are
needed. Ohashi & Kataoka (2019) used the polarized
observations of HD 163296 (Dent et al. 2019) to model the
dust content within the disk. They found that the gaps should
have a large contribution of some-hundred-micrometer sizes to
explain the polarimetric observations; however, the rings could
have no contribution of these grains but much larger grains. No
polarized emission at low angular resolution (2.0–2.6) was
detected at Band 7 for GM Aur and MWC 480 (Hughes et al.
2013); however, Harrison et al. (2019) detected polarized
emission parallel to the disk minor axis in the inner disk of
MWC 480, with a polarization fraction of ~1%, which can be
explained by self-scattering of large grains (>1 mm). More
multiwavelength polarized observations are needed to model
and discern the maximum grain sizes in these disks.

In particular, the small-grain solution has two important
problems. First, it requires a large dust mass to produce sufficient
millimeter emission, and in that case, the optically thick inner
disks of IM Lup, HD 163296, and MWC 480 would be
gravitationally unstable if the gas-to-dust ratio is 100, as shown
in the left panel of Figure 15. The maximum Stokes numbers of
these grains (right panel of Figure 15) are ~2 orders of
magnitude smaller compared with the large-grain solution, and
they would be perfectly coupled with the gas, making it difficult
for the sharp ring structures observed in continuum emission to
exist. Second, this regime is only able to explain the inner disk.
Only millimeter-sized particles can explain the outer disk
emission, providing continuous dust profiles from the inner to
the outer disk, with no sharp transition.

Despite these physical problems associated with the small-
grain solution, it is consistent with the spherical grain sizes
expected from the polarization patterns observed in many disks
(e.g., Kataoka et al. 2015; Bacciotti et al. 2018; Hull et al.
2018; Mori et al. 2019), which are mainly tracing the
unresolved emission from the inner disks. However, since the
polarization degree of oblate dust grains is higher than that
from spheres (Kirchschlager & Bertrang 2020), the large-grain
solution may also be consistent with the observed polarization
properties.

Observations at larger wavelengths can also be used to
discern between the small- and large-grain solutions, even if the
observations are nonresolved. The expected 7 mm flux from the
large- and small-grain models differs by a factor of two. For
example, the estimated IM Lup flux at 7 mm for the large- and
small-grain models is ~1.9 and 1.1 mJy, respectively, and for
HD 163296, the fluxes are ~3.9 and 2.2 mJy, respectively.
Using the Australian Telescope Compact Array (ATCA),
Lommen et al. (2010) reported a total flux of 2.2 ± 0.16 mJy
at 6.8 mm for IM Lup, and using the Very Large Array (VLA),
Isella et al. (2007) reported a total flux of 4.5 ± 0.5 mJy for HD
163296. In the latter work, they also determined that the free-
free contribution at this wavelength is 1.2 mJy (~27%), such
that the dust continuum flux of HD 163296 at 7 mm is
3.3 ± 0.5 mJy. Figure 16 shows the SED for IM Lup and HD

163296 for the large- and small-grain solution and their
observed flux at 7 mm. For IM Lup, the free-free contribution
at 7 mm has not been estimated, but it should be ~14% or 50%
of the total flux in order to be consistent with the large- or
small-grain model, respectively. For HD 163296, the 7 mm
flux from the large-grain solution is slightly larger than the
observed continuum flux and differs from the small-grain
model by a factor of 1.5. Thus, the total flux at 7 mm provides
an additional constraint that suggests that the maximum grain
sizes around HD 163296 are millimeter sized.

Future polarized observations at long wavelengths (e.g.,
λ = 7 mm) will also help us to constrain grain sizes in the inner
disks, and more complex modeling should be formulated to
explain the whole structure of the inner disk with scattering, for
example, a two-layer model with small grains in the disk
surface (to explain the spectral indices) and large grains in the
midplane (dominating the dust opacity and mass).

Finally, note that the optical depths with scattering are larger
than the nonscattering case by one order of magnitude or more
(Figure 13). This large increase occurs because the albedo of
millimeter grains at radio frequencies is large (ω<sub>ν</sub> > 0.9;
Birnstiel et al. 2018), and the total optical depth increases by
a factor of 1/(1-ω<sub>ν</sub>) > 10. The disks tend to be close to the
optically thick regime (τ<sub>ν</sub> = 1) when scattering opacity is
turned on. In the small-grain solution (bottom panels in
Figure 14), the optical depths at Band 6 are optically thick,
while the optical depth at Band 3 is optically thin in all disks. In
the next sections we analyze the results from each disk,
focusing on the dust properties derived from the nonscattering
model and the scattering model with large grains.

5.2. Source Specific

5.2.1. IM Lup

IM Lup is a ~1.0 Myr old T Tauri star with a mass of
1.1 M<sub>☉</sub> (Andrews et al. 2018). It is the only disk in our sample
with spiral arms traced in the millimeter dust continuum.
A complete study of its morphology was done by Huang et al.
(2018b) using high-angular-resolution observations with ALMA
at Band 6 (44 mas or ~4 au resolution). They also found a ring at
134 au, which is plotted as a vertical dashed line in our figures as a
reference. In addition, Law et al. (2021) found a dark and a bright
ring at 209 and 220 au, respectively.

This disk has also been studied by multiwavelength
observations (optical, near-infrared, and millimeter emission)
by Pinte et al. (2008). They determined that dust grains are
millimeter sized in the disk midplane (consistent with our
results) and suggest that the disk may be gravitationally
unstable.

We found that the dust surface density (Figure 7)
monotonically decreases from ~3 g cm⁻² in the inner disk
(20 ~ au) to ~0.1 g cm⁻² in the outer disk (~140 au), while the
maximum grain size (Figure 8) decreases from ~1 cm in the
inner disk to several hundred micrometers in the outer disk.
The dust mass is 1.92<sup>+1.91</sup><sub>-0.46</sub> × 10<sup>-3</sup> M<sub>☉</sub> in the nonscattering
model and 3.64<sup>+1.99</sup><sub>-1.41</sub> × 10<sup>-3</sup> M<sub>☉</sub> in the scattering model with
large grains. The former is close to the dust mass of
0.7 × 10<sup>-3</sup> M<sub>☉</sub> in Pinte et al. (2008) and similar to 1.6 × 10<sup>-3</sup>
M<sub>☉</sub> in Cleeves et al. (2016).<sup>28</sup>

<footnote n="28"><sup>28</sup> The dust masses reported in previous works are updated using the distances summarized in Table 1 in Öberg et al. (2021).</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** SED for IM Lup (left) and HD 163296 (right). The fluxes are measured within the radius where the small-grain solution is found. The squares are the MAPS observations, while the circles are the predicted fluxes from the large (red) and small (blue) grain models. The green error bar is the nonresolved flux measured at long wavelength.</figcaption></figure>

The optical depth of the emission in IM Lup appears optically thick within ~30 au and optically thin beyond this radius (Figure 9). When scattering is taken into account (Figure 13), the total opacity is close to or above 1 for all wavelengths at all radii. The Toomre *Q* parameter (Figure 10) is close to or below 1.7 at all radii, indicative of a gravitationally unstable system. Such a small value is expected given the observed spiral arm structure in Band 6 (Figure 1). This is the only disk in our sample where the Stokes number decreases with radius beyond ~75 au, and its overall value is small ($1 \times 10^{-2}$ and $3 \times 10^{-2}$). This behavior occurs because of the high surface density, which tends to a constant value beyond ~100 au, while the maximum grain size monotonically decreases with radius. Such small values of the Stokes number are expected in dust trapping models by spiral arms when the fragmentation velocity is set to $10 \text{ m s}^{-1}$ (Dipierro et al. 2015a).

### 5.2.2. GM Aur

GM Aur is a T Tauri star with an age within ~3 and 10 Myr and a mass of $1.1 M_{\odot}$ (Macías et al. 2018). The GM Aur disk is a transition disk with a large inner cavity previously reported by many authors (e.g., Hughes et al. 2009; Andrews et al. 2011; Guilloteau et al. 2011; Hornbeck et al. 2016) and recently studied using multiwavelength observations by Macías et al. (2018) and Huang et al. (2020) at high angular resolution (~8 au in ALMA and ~26 au in VLA). They found two rings at 40 and 84 au, which are plotted as vertical dashed lines in our figures as a reference.

The spectral indices of this disk have local minima at the ring position (Figure 3). However, because of the low angular resolution in Band 3 and beam averaging of Band 6, it is not possible to clearly distinguish high contrasts between rings and gaps. When the spectral index is computed between Bands 4 and 6 using higher angular resolution profiles, the contrast between rings and gaps is clearer, as shown in Huang et al. (2020). The dust surface density (Figure 7) increases within the first 40 au, reaching a maximum value in the inner ring position, and then it decreases toward the outer disk, with a local maximum in the outer ring. The maximum grain size (Figure 8) monotonically decreases from the inner disk, where

$a_{\text{max}} \sim 8 \text{ mm}$, to the outer disk, where $a_{\text{max}} \sim 1 \text{ mm}$. GM Aur is optically thin for the entire disk (Figure 9). The largest optical depth is reached in the inner ring at Band 7, but it does not exceed the optically thick regime. When scattering opacity is included, the total opacity increases at all wavelengths, such that the optical depths are just below 1 at the ring positions. This disk is not consistent with a small-grain solution ($a_{\text{max}} < 300 \ \mu\text{m}$) in the scattering case owing to its small optical depth.

The dust mass is $0.74_{-0.13}^{+0.34} \times 10^{-3} M_{\odot}$ in the nonscattering model and the same for the scattering model with large grains, consistent with the dust mass in the best-fit model in Schwarz et al. (2021). This dust mass is a factor of ~3 smaller than the one derived in Macías et al. (2018) and a factor of 2 smaller than the estimation in McClure et al. (2016). The Toomre parameter *Q* (Figure 10) indicates that the disk is gravitationally stable, with a minimum close to the position of the outer ring at ~84 au. Based on the CO emission, Schwarz et al. (2021) found that the region between 70 and 120 au could be gravitationally unstable, due to the total gas mass being between 0.2 and 0.3 $M_{\odot}$. However, there is no evidence of this possible instability in the dust continuum maps. The maximum Stokes number is close to 0.1 in the inner disk and slowly decays to 0.05 at ~60 au. Then, it increases again in the outer disk.

### 5.2.3. AS 209

AS 209 is a ~1 Myr old T Tauri star with a mass of $1.2 M_{\odot}$ (Andrews et al. 2018). Its surrounding disk has been previously studied with ALMA Band 6 at a high angular resolution of 0.037 or ~4.5 au by Guzmán et al. (2018), revealing seven bright rings and a central component. From our radial profiles (Figure 2), we are able to distinguish the rings at 74.1 and 120.4 au, which are plotted together with the 15.1 au ring (Guzmán et al. 2018) as vertical dashed lines in our figures as a reference. Fedele et al. (2018) found that the deepest gap at ~95 au is consistent with the presence of a planet with a mass of 0.2 $M_{\text{Jup}}$, while the inner gap at ~60 au could be generated by a planet with an upper limit mass of 0.1 $M_{\text{Jup}}$. The latter is under debate, as it is not required to explain the observed structures (Zhang et al. 2018; Alarcón et al. 2021).
At both gaps, we found some evidence of local minima of dust mass and grain size. The dust surface density depletion is higher at 95 au than at 60 au gap location; at 95 au the dust surface density decreases by ~0.6 dex compared with the density in the outer ring at 120 au. This depletion may be larger based on the high contrasts observed in Band 6 (Guzmán et al. 2018).

The dust surface density and grain size have local maxima in the rings, which is consistent with dust trapping models (e.g., Lyra & Lin 2013; Ruge et al. 2016; Sierra et al. 2019), where larger grains are trapped and the dust mass is enhanced. Dullemond et al. (2018) found that the widths of both rings are smaller than the local scale height, which also supports the idea of dust trapping in these rings. The dust surface density decreases from ~1 g cm⁻² in the inner disk (~20 au) to ~0.03 g cm⁻² in the outer disk (~130 au), similar behavior to the dust surface density computed in Pérez et al. (2012) using particle size distributions with p = 3.0. The total dust mass we derive is 9.1±7.6 × 10⁻⁴ M☉ in the nonscattering model and 7.5±7.2±1.8 × 10⁻⁴ in the scattering model with large grains. These dust masses are a factor of ~2–4 larger than the dust mass of 2.6 × 10⁻⁴ M☉ estimated from the thermal continuum at 345 GHz (Andrews et al. 2009), the fiducial model of Fedele et al. (2018) with 3.2 × 10⁻⁴ M☉, or polarimetric data with SPHERE (Avenhaus et al. 2018) with a dust mass of 2.32 × 10⁻⁴ M☉.

The Toomre *Q* parameter for AS 209 reveals that this seems to be gravitationally stable. The estimated maximum Stokes number is between 2 × 10⁻² (at ~25 au) and 3 × 10⁻¹ (at ~100 au), with local minima in rings and local maxima in gaps (due to dust surface density morphology). The inferred Stokes numbers in the rings are consistent with the lower limits computed in Dullemond et al. (2018) (they assume a grain size of 200 µm) if they are scaled to the grain size derived in this work. Rosotti et al. (2020) estimated the ratio between the turbulent *α* parameter and the Stokes number in rings of AS 209 and HD 163296. In the case of AS 209, they found that *α*/St = 0.18 ± 0.04 at 74 au and *α*/St = 0.13 ± 0.02 at 120 au. This means that, using our estimated maximum Stokes number, the turbulent parameters at the bright ring location are *α* ~ 1.3 × 10⁻² at 74 au and ~9.1 × 10⁻³ at 120 au, which are close to the fragmentation limit computed in that work, i.e., where the dust growth is prevented owing to turbulent fragmentation.

### 5.2.4. HD 163296

HD 163296 is a Herbig Ae star with an age >6 Myr and a mass of 2.0 M☉ (Andrews et al. 2018). Its disk has been previously studied with ALMA Band 6 at high angular resolution (42 mas or 4.2 au; Isella et al. 2018), revealing four bright rings that were already identified by Isella et al. (2016). From our radial profiles, we are able to distinguish the rings at 14.4, 67.0, and 100.0 au, which are plotted as vertical dashed lines in our figures as a reference.

At the gap and ring positions, we obtained local minima and maxima in the dust surface density (Figure 7), respectively. However, the maximum grain size only has a maxima in the ring at 14.4 and 100 au. The ring at 66 au does not show strong evidence of dust trapping at this resolution. This can also be seen in the spectral index profiles (Figure 3), where the minimum at 100 au is deeper compared with that at 66 au. This is consistent with the results from Dullemond et al. (2018), where the resolved ring width at 66 au is larger than the local

scale height, disfavoring the dust trapping scenario. This is not the case for the ring at 99 au, which has a resolved ring with smaller than the local scale height. Then, the origin of the dust ring at 66 au seems to be different from that of the outer ring at 99 au. In particular, Zhang et al. (2021) found that the midplane CO snowline is located at 65 au. Thus, the origin of the dust ring at this radius seems to be a traffic jam caused by the CO depletion, where a continuum ring is expected owing to the sintering effect at a snowline (e.g., Okuzumi et al. 2016).

The estimated dust mass is 0.90±0.79 × 10⁻³ M☉ in the nonscattering model and 0.83±1.35±0.18 × 10⁻³ M☉ in the scattering model with large grains. Both dust masses are a factor of 2 smaller than that estimated using ALMA and SPHERE observations in Muro-Arena et al. (2018), but it is within the range (if the dust gas is 100 times smaller than the gas mass) estimated by the nonresolved observations of Isella et al. (2007) and Tilling et al. (2012), where the expected dust mass lies within 0.5 × 10⁻³ M☉ and 0.9 × 10⁻³ M☉. Booth et al. (2019) computed a total gas mass of 0.31 M☉, which implies a dust-to-gas mass ratio of 1/260 when compared with the dust mass from Isella et al. (2007). The emission of HD 163296 in all bands is optically thin, except the innermost ring at Band 6, where the optical depth is around 1 (Figure 9). If scattering is included, the total opacity is just below 1 in the outer rings at 66 and 100 au (Figure 13). Note that the disk is optically thin beyond ~40 au even when dust scattering opacity is included, which is consistent with the discussion of Zhu et al. (2019) based on the nontotal extinction of the CO emission in the rings.

The Toomre parameter *Q* (Figure 10) indicates that the entire disk is gravitationally stable with local minima at the ring position. The Toomre parameter from Booth et al. (2019) is ~6 at 110, close to the ~4 estimated in this work at the same position. The maximum Stokes number lies within ~10⁻² in the inner disk and reaches a value of 0.3 in the outer disks, with local minima at the rings. Similarly to AS 209, the Stokes numbers in the rings are consistent with the lower limits in Dullemond et al. (2018). According to the Rosotti et al. (2020) model, the ratio between the turbulent parameter and the maximum Stokes number at 67 and 100 au is 0.23 ± 0.03 and 0.04 ± 0.01, respectively. Using our estimates for the maximum Stokes number, this means that the turbulent parameters at these ring radii are *α* ~ 9.2 × 10⁻³ and ~2.8 × 10⁻³, respectively, consistent with constraints from line observations (e.g., Flaherty et al. 2017). Similar to AS 209, these values are close to the turbulent fragmentation limit, which could be preventing dust growth at these positions, where, according to the dust trapping models, the effects of radial drift become inefficient (e.g., Whipple & Planet 1972).

### 5.2.5. MWC 480

MWC 480 is a 7 Myr old Herbig Ae star with a mass of 2.1 M☉ (Simon et al. 2019). The MWC 480 disk was studied with ALMA at Band 6 with high angular resolution (0.′′14, ~22 au) by Long et al. (2018) and Liu et al. (2019), finding a gap centered at ~74 au and a ring at ~98 au. Additionally, Long et al. (2018) (using the residual/uv data) and Law et al. (2021) (studying the radial profiles from the images) inferred a dark and bright ring at 149 and 165 au, respectively. The bright rings at 98 and 166 au can be identified in our data, and vertical dashed lines are plotted in our figures as a reference.
Liu et al. (2019) find that the observed substructure is consistent with an embedded planet with a mass of 2.3 $M_{\text{Jup}}$ at a radius of 78 au. At this radius, we found a minimum in the dust surface density (Figure 7) and no clear evidence of a local minimum for the maximum grain size (Figure 8). However, this minimum cannot be dismissed, because our angular resolution is larger than the width of the gap resolved with higher angular resolution. Similarly, at the ring position we found a maximum in the dust surface density, but the maximum grain size does not have a local maximum. In general, the dust surface density decays from ~1 g cm⁻² in the inner disk (~20 au) to ~0.03 g cm⁻² in the outer disk (~120 au), while the maximum grain size is ~4 mm in the inner disk and ~2 mm in the outer disk.

The estimated dust mass is $0.97_{-0.18}^{+0.16} \times 10^{-3} M_{\odot}$ in the nonscattering model and $1.19_{-0.29}^{+1.77} \times 10^{-3} M_{\odot}$ in the scattering model with large grains. The dust mass in the nonscattering model is within a factor of ~1.8–2.6 smaller than that found in Piétu et al. (2006; ~3.2 $\times 10^{-3} M_{\odot}$), Guilloteau et al. (2011; ~2.3 $\times 10^{-3} M_{\odot}$), and Liu et al. (2019; $1.6_{-0.4}^{+0.5} \times 10^{-3} M_{\odot}$). The dust mass with scattering is within a factor of ~1.3–2 smaller compared with the same references. The emission in Band 3 appears to be optically thin, while for Band 6 the emission is optically thin outside of ~35 au. However, when scattering opacity is included, the total optical depth at Band 6 is optically thick within ≲60 au. The Toomre *Q* parameter is larger than 2 at all disk radii, such that no evidence of gravitational instability is found. Finally, the maximum Stokes number is ~9 $\times 10^{-3}$ in the inner disks and increases to ~0.1 in the outer disk.

## 6. Summary and Conclusions

In this paper we performed a multiwavelength analysis of the disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480 using ALMA observations at Bands 6 and 3 (and archival data for AS 209 in Band 7 and GM Aur in Bands 4 and 7). The Band 6 observations consist of two intraband data sets at a central frequency of 226 and 257 GHz (1.33 and 1.17 mm), while the two intraband data sets at Band 3 with central frequencies of 94 and 106 GHz (3.20 and 2.84 mm) are merged to obtain an image with a higher S/N at a central frequency of 100 GHz (3.0 mm). These observations provide at least three wavelengths that are used to fit the spatially resolved continuum spectrum from the intensity radial profiles. This fitting is performed by modeling the (sub)millimeter spectrum using the radiative transfer equation for two cases: when scattering is taken into account and when it is neglected. When scattering is considered, two possible solutions are consistent with observations in the inner disks of IM Lup (≪60 au), HD 163296 (≪40 au), and MWC 480 (≪60 au). The two solutions are characterized by a large-grain size regime ($a_{\text{max}} > 300$ µm) and small-grain size regime ($a_{\text{max}} < 300$ µm).

The dust temperature is fixed by the midplane temperature of Zhang et al. (2021), and a large parameter space of possible dust surface densities and grain sizes is explored. From this procedure, we constrain the radial profiles of dust surface density, maximum grain size, and optical depths at all observed wavelengths, which are able to explain the spectral properties of each disk. The dust opacity properties are computed using a slope of the particle size distribution of $p=2.5$, which represents the case where the maximum grain size is limited by drift (Birnstiel et al. 2012). This slope also gives a lower limit to the maximum grain size. Once these properties are

computed, radial profiles of the Toomre *Q* parameter and the maximum Stokes number are estimated for the five disks. Our main conclusions can be summarized as follows.

1. Millimeter grain sizes are inferred in the disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480 in the nonscattering model and scattering model in the large-grain size regime. The maximum grain size radial profiles have a negative slope from the inner disk (where the grain sizes are close to 1 cm) to the outer disk (where the grain sizes are slightly below 1 mm). The maximum grain size locally peaks in most of the known rings in these disks. In a few cases we do not find evidence for large grains in the ring locations, but this could be due to insufficient angular resolution. In these models, the dust surface density decreases with disk radius in all the disks in our sample, but with local maxima that coincide with the ring positions. Then, both the dust mass and the grain size are enhanced within most of the rings, which is consistent with models where pressure maxima act as dust size differential traps (large grains are more concentrated around the pressure maxima than small grains).

2. The inner disks (≪20 au) in IM Lup, HD 163296, and MWC 480 are found to be optically thick at Band 6, even if scattering opacity is not included. When scattering is taken into account, the total optical depths increase by a factor of 10; this occurs because the albedo of millimeter grain sizes is large at radio frequencies.

3. Grains of a few hundred micrometers in size are consistent with the emission from the inner disks of IM Lup, HD 163296, and MWC 480, when scattering is taken into account. This alternative solution occurs because scattering modifies the optically thick spectral indices and degeneracy arises. The disks around GM Aur and AS 209 are not consistent with some-hundred-micrometer sizes because they are optically thin, even in their inner disks. This particular solution has two problems: it requires a large amount of dust in small grains to reach the level of emission at all wavelengths, making the inner disks gravitationally unstable, and it is unable to reproduce the observed emission at all wavelengths for the outer disks. Additionally, the expected total flux of IM Lup and HD 163296 at $\lambda=7$ mm from the few-hundred-micrometer grain model is low compared with their observed flux at this wavelength, while it is consistent with the expected flux from millimeter grain sizes.

4. Our results strengthen the idea that IM Lup (which presents spiral arm structures) is a gravitationally unstable disk, as our estimated Toomre *Q* parameter is lower than 2 outside of ~15 au. In contrast, the Toomre *Q* parameter of GM Aur, AS 209, HD 163296, and MWC 480 is larger than 2 at all disk radii.

5. We estimate maximum Stokes numbers below 1 for all the disks. In general, the Stokes numbers in the four disks increases with the disk radius (except IM Lup), with local minima at the ring positions, where the dust surface density is a maximum. The estimated maximum Stokes numbers, combined with the constraints on the $\alpha$/St ratio of Rosotti et al. (2020), suggest that dust growth could be limited by turbulent fragmentation in the rings of AS 209 at 74 and 120 au and in the rings of HD 163296 at 67 au and 100 au.
6. Scattering is an important component of opacity in protoplanetary disks. The models where scattering is taken into account provide a more detailed approximation to the radiative transfer problem in disks. In many cases, the scattering model reduces to the nonscattering model because the emission is optically thin, where the scattering effects can be neglected. However, the main differences are found in the inner disks, where the optical depth can be larger than 1.

Observations at high angular resolution and sensitivity at multiple wavelengths are needed to better constrain the dust properties in protoplanetary disks. The degeneracy from the scattering models could be avoided if the (sub)millimeter SED is better sampled. For example, observations at longer wavelengths (e.g., 7 mm) can help to reach the midplane in the inner disk, where the ALMA wavelengths could not be tracing the dust properties. The solution that requires few hundreds of micrometer sizes in the inner disks could be rejected or strengthened with such observations. Also, higher angular resolution observations can be used to infer higher contrasts between the dust properties in rings and gaps, and this will help us to understand the origin of these substructures. Finally, further modeling that incorporates polarization data would help to simultaneously explain the polarization patterns and spectral index, such that the solutions with small and large grains in the scattering models could be consistent.

We are grateful to the anonymous referee for useful comments that helped to clarify some aspects of the paper.

This paper makes use of the following ALMA data: ADS/JAO. ALMA#2015.1.00678.S, ADS/JAO.ALMA#2018.1.01055.L, ADS/JAO.ALMA#2017.1.01151.S ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

A.S. and L.M.P. acknowledge support from ANID/CONICYT Programa de Astronomía Fondo ALMA-CONICYT 2018 31180052. A.S. thanks Enrique Macias, Carlos Carrasco, and Osmar Guerra for helpful discussions. L.M.P. acknowledges support from ANID project Basal AFB-170002 and from ANID FONDECYT Iniciación project No. 11181068. C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant DGE1745303. V.V. G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. K.Z. acknowledges the support of the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin–Madison with funding from the Wisconsin Alumni Research Foundation and the support of NASA through Hubble Fellowship grant HST-HF2-51401.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. A.S.B acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). A.D.B. acknowledges support from NSF AAG grant No. 1907653. C.W. acknowledges financial support from the University of Leeds, STFC,

and UKRI (grant Nos. ST/R000549/1, ST/T000287/1, MR/T040726/1). F.L. and R.T. acknowledge support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. F.M. acknowledges support from ANR of France under contract ANR-16-CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP “Origins of Life”). G.C. is supported by NAOJ ALMA Scientific Research grant code 2019-13B. I.C. was supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51405.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.B. acknowledges support by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. S.M.A. and J.H. acknowledges funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. J.H. acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51460.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.B.B. acknowledges support from NASA through the NASA Hubble Fellowship grant No. HST-HF2-51429.001-A, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.D.I. acknowledges support from the Science and Technology Facilities Council of the United Kingdom (STFC) under ST/T000287/1. K.I.Ö. acknowledges support from the Simons Foundation (SCOL No. 321183) and an NSF AAG grant (No. 1907653). R.L.G. acknowledges support from a CNES fellowship grant. T.T. is supported by JSPS KAKENHI grant Nos. JP17K14244 and JP20K04017. Y.Y. is supported by IGPEES, WINGS Program, the University of Tokyo. Y.A. acknowledges support by NAOJ ALMA Scientific Research grant code 2019-13B and Grant-in-Aid for Scientific Research 18H05222 and 20H05847. Y.L. acknowledges the financial support by the Natural Science Foundation of China (grant No. 11973090).

*Software:* Astropy (Astropy Collaboration et al. 2013, 2018), CASA (McMullin et al. 2007), Frankenstein (Jennings et al. 2020), Matplotlib (Hunter 2007), Numpy (Harris et al. 2020)

## Appendix A
### Imaging and Radial Profiles

To analyze how the radial intensity profiles in Figure 2 depend on how we image the visibilities, we compute the radial profiles in three different ways and then compare them. Here we define $\theta_{B3}$ as the beam obtained from imaging the Band 3 data set with robust = -2 (uniform weighting), which always corresponds to the largest beam of our multiwavelength data set. The different ways one can derive the intensity profiles are as follows:

1. Imaging of all data sets with robust = -2 (uniform weighting), then convolving to a common resolution given by $\theta_{B3}$, and, finally, azimuthally averaging these common resolution images to obtain the intensity radial profiles.
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Convolved intensity profiles at B3 (red), B6(I) (green), and B6(II) (blue) of each disk (see legend in the upper right corner of each panel), and three different methods: fitting the visibilities (circles), robust = 0 and convolution (crosses), and robust = -2 and convolution (shaded area).</figcaption></figure>

2. Imaging the Band 6 data sets with robust = 0, then convolving to a common resolution given by $\theta_{B3}$, and, finally, azimuthally averaging these common resolution images to obtain the intensity radial profiles.

3. Fitting the visibilities with a nonparametric model to derive a radial profile of the emission, which we convolve to the common $\theta_{B3}$ resolution.

The radial profiles used in this work correspond to the first method proposed. Due to the uncertainty in using a robust parameter of -2, and because the extended emission could be filtered out in Band 6, we adopt a second method: imaging the Band 6 data with a robust parameter of 0 (still higher angular resolution than $\theta_{B3}$), and then convolving these images to the same angular resolution given by $\theta_{B3}$. Azimuthally averaged intensity profiles are also computed over the image plane in this case. In the third case, we use the Frankenstein code (Jennings et al. 2020) to fit the visibility data (assuming that the disks are axisymmetric and using the standard hyperparameters $\alpha = 1.05$, $w_{smooth} = 10^{-4}$) and obtain deconvolved intensity profiles that did not go through the imaging process of the first or second method. These intensity profiles are then convolved to the same angular resolution $\theta_{B3}$.

Figure 17 shows the intensity profiles derived from the three methods for all disks. The shaded areas correspond to the radial profiles used in this work (method 1, robust = -2). The crosses correspond to the radial profile from method 2 with a robust parameter of 0 (the error area is not included because they overlap to those in method 1). The circles correspond to the radial profiles from fitting the visibilities (method 3). All these

profiles agree, except IM Lup Band 3, where we were not able to use Frankenstein to model the visibilities, and for that reason it is not shown.

Figures 18 and 19 show the brightness temperature images of the final products from CLEAN (method 1) and Frankenstein (method 3), respectively. The maps from CLEAN are noisier than those shown in Figure 1, since they are computed using a smaller robust parameter. A robust = -2.0 improves the angular resolution but leads to reduced S/N. The latter is improved by azimuthally averaging the data, as the uncertainty in the radial profiles decreases by a factor of $\sqrt{n}$, where n is the number of beams within a ring of the disk.

The images from Frankenstein are similar to those obtained from CLEAN, and the radial profiles match (Figure 17), so then we confirm that false structures are not introduced when imaging the data using a robust parameter = -2.

Recently, the DSHARP observations of IM Lup, AS 209, and HD 163296 were modeled in Jennings et al. (2021) using Frankenstein. They found deeper gaps and brighter rings in the inner region (~10 au) of AS 209 and HD 163296 compared with the previous constraints in Huang et al. (2018a). The DSHARP observations are able to resolve these small substructures (compared with this work), as they probe baselines upon 10 Mλ, while the observations presented here cover baselines <~2 Mλ. Our goal, rather than finding small-scale substructure (with all the continuum substructures presented in Law et al. 2021), is to double-check that the radial profiles from CLEAN can also be recovered by an alternative methodology.
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Brightness temperature images of the disks obtained from CLEAN with a robust parameter = -2. The top, middle, and bottom panels are observations at 100, 226, and 257 GHz, respectively.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Brightness temperature images of the disk obtained from Frankenstein. The top, middle, and bottom panels are the observations at 100, 226, and 257 GHz, respectively.</figcaption></figure>
## Appendix B

### Degeneracy of the Inferred Dust Properties

Here we explore the degeneracy in the inferred dust parameters in the scattering and nonscattering cases. We assume that the dust temperature is known, as we did in our analysis. Then, given a maximum grain size ($a_{max}$) and dust surface density ($\Sigma_d$), the intensities ($I_\nu$) at many wavelengths can be computed (Equations (1) or (2)),

$$ (a_{\max}, \Sigma_d) \Rightarrow (I_{\nu_1}, I_{\nu_2}, \dots, I_{\nu_n}). \qquad (B1) $$

We want to know whether there is a different set of parameters ($a'_{\max}$, $\Sigma'_d$) such that they can have a similar continuum spectrum, i.e.,

$$ (a'_{\max}, \Sigma'_d) \Rightarrow (I'_{\nu_1}, I'_{\nu_2}, \dots, I'_{\nu_n}) \approx (I_{\nu_1}, I_{\nu_2}, \dots, I_{\nu_n}). \qquad (B2) $$

We define the probability of inferring ($a'_{\max}$, $\Sigma'_d$) given the “real” model ($a^r_{\max}$, $\Sigma^r_d$) as

$$ p(a'_{\max}, \Sigma'_d | a^r_{\max}, \Sigma^r_d) = \exp(-\chi^2/2), \qquad (B3) $$

where the chi-squared statistic is given by

$$ \chi^2 = \sum_{\nu} \left( \frac{I_{\nu} - I'_{\nu}}{\epsilon_{\nu}} \right)^2, \qquad (B4) $$

where $\epsilon_\nu = 0.1 I_\nu$ is the error in the intensities, which is assumed to be the same as the uncertainty from the flux calibration. We assume that the observed frequencies are the same as in this work: $\nu = 257, 226$, and 100 GHz.

To explore possible degeneracies in our model, we include random fluctuations in the intensities and look for parameters that result in an equivalent continuum spectrum. For this, given a set of “real” parameters ($a^r_{\max}$, $\Sigma^r_d$), the intensities $I_\nu$ over our wavelengths are computed. Then, a random flux factor is applied independently to each wavelength, and one can find a set of inferred parameters $a'_{\max}$, $\Sigma'_d$ that have a similar continuum spectrum given their probability (Equation (B3)).

This process is repeated many times using a random flux factor that follows a Gaussian distribution centered at one with a width of 0.1 (the uncertainty in the flux calibration). When a set of parameters ($a'_{\max}$, $\Sigma'_d$) with a high probability is found (we choose $p > 0.9$), these properties are plotted against the real parameters ($a^r_{\max}$, $\Sigma^r_d$).

Figure 20 shows the results for this process in a model with scattering properties. The left, middle, and right panels correspond to a model where the real dust surface density is $\Sigma^r_d = 0.1, 1.1, 11.0 \text{ g cm}^{-2}$, respectively. These values were chosen to cover different optical depth regimes. The top panels show the inferred grain size as a function of the real grain size, and bottom panels show the inferred dust surface density as a function of the real grain size. Different colors are used to represent inferred parameters that overestimate (blue), underestimate (red), or recover (green) the real dust surface density.

In principle, the realizations shown in the top panels of Figure 20 should follow a line with a slope of 1, while the

bottom panels should follow a horizontal line at the real dust surface density value. Although most realizations recover the original dust properties, close to this line, there are some where the inferred dust surface density and grain size differ from the real ones.

For example, in the case of $\Sigma^r_d = 0.1 \text{ g cm}^{-2}$ (optically thin regime), the continuum spectrum of dust grains with a real size of $\lesssim 100 \text{ µm}$ is similar to that of grains with a size of $\sim 1 \text{ mm}$ (red horizontal branch), but with a smaller dust density ($\Sigma'_d \sim 10^{-2} \text{ g cm}^{-2}$). This is due to the degeneracy of the opacity spectral index (middle panel of Figure 4 in Birnstiel et al. 2018). In the same way, the continuum spectrum of millimeter grains is similar to that of small grains (blue vertical branch), but with a larger inferred dust density. The continuum spectrum of grains larger than $\sim 1 \text{ mm}$ is also consistent with models that overestimate the dust mass and grain size, or with models that underestimate the dust mass and grain size.

In the middle panels ($\Sigma^r_d = 1.1 \text{ g cm}^{-2}$), the disks are in a transition between the optically thin and optically thick regime. The opacity of small grains ($a^r_{\max} < 100 \text{ µm}$) is very small at millimeter wavelengths, such that the emission is still in the optically thin regime, and the degeneracy is equivalent to the previous case. The continuum spectrum of millimeter grains is similar to that of grain sizes around $\sim 200-300 \text{ µm}$ that overestimates the dust density (blue horizontal branch). The degeneracy of grains larger than some millimeters is similar to the optically thin case, but with a significant increase of the number of models that can fit the continuum spectrum.

In the right panels ($\Sigma^r_d = 11 \text{ g cm}^{-2}$), where the disks are optically thick, the degeneracy between the small grains ($a^r_{\max} < 100 \text{ µm}$) and millimeter grains disappears. However, the degeneracy between millimeter grain sizes and hundred-micrometer grains remains. This degeneracy extends to centimeter grain sizes (horizontal red branch).

In the optically thin regime (bottom left panel of Figure 20), the dust surface density could be underestimated for $\lesssim 100 \text{ µm}$ and for centimeter grains. For $\gtrsim 3 \text{ mm}$ grains, the number of models that tends to overestimate the dust mass is larger than those that underestimate it. In the intermediate regime (bottom middle panel of Figure 20), the mass is highly overestimated for millimeter grains, while for micrometer and centimeter grains they are equivalent to the previous optically thin case. Finally, in the optically thick regime (bottom right panel of Figure 20) the mass is well constrained for grains smaller than $\lesssim 100 \text{ µm}$, but it could be overestimated or underestimated for millimeter or larger grains. These results are for the wavelengths considered in this work; the addition of more observations at different wavelengths (submillimeter and centimeter) will help better constrain the dust properties in these objects. In models where the temperature is also unknown, the degeneracy between the models will increase, unless observations at more wavelengths can help reduce the number of equivalent models.
<figure><img src="image_22.png" /><figcaption>Figure 20. Maximum grain size and dust surface density degeneracy for the particular set of wavelengths observed. The left, middle, and right panels correspond to a model with a dust surface density of 0.1, 1.1, 11.0 g cm⁻², respectively. Top and bottom panels are the inferred maximum grain size and dust surface density, respectively, that have a similar SED at the three observed wavelengths (B3, B6(I), B6(II)). The blue, red, and green colors represent models where the dust surface density is overestimated, underestimated, and recovered, respectively.</figcaption></figure>

## Appendix C
Dust Temperature Profiles

In order to avoid degeneracy, especially in the optically thin regime, we fix the midplane temperature profiles using the models from Zhang et al. (2021). In their model, they iterate the density properties of a population of small grains (with sizes of 1 µm), large grains (1 mm grains), and gas and generate dust continuum images at Band 6, which are then compared with the observed images at high resolution (Figure 1). From this procedure, they obtain a bidimensional structure (density and temperature) for gas and both dust populations.

We define the dust midplane temperature profiles as

$$T_{\text{mid}}(r) = \frac{\int \rho_{\text{d}}(r, z) T(r, z) dz}{\int \rho_{\text{d}}(r, z) dz}, \qquad (\text{C5})$$

where $\rho_{\text{d}}(r, z)$ is the density from the 1 mm grains and $T(r, z)$ is the dust temperature at the cylindrical coordinates $(r, z)$. The integral is evaluated within an aspect ratio of $z/r = 1/10$. We tested changing this aspect ratio, but the final temperature profiles are only slightly modified. This occurs because the temperature from the denser regions close to the midplane dominates the integrals.

Figure 21 shows the midplane temperature for the five disks (blue lines). Additionally, for reference we fitted the radial profiles to a power-law function of the radius (dashed orange line)

$$T = T_{50}(r/50 \text{ au})^{-q}. \qquad (\text{C6})$$

The magnitude of the temperature at 10 au and the slope of the power law from the fit are summarized in Table 4.

In all cases, the midplane profiles are convolved to the angular resolution of each disk and then used to fit the dust surface density and maximum grain size (Section 4).
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Dust temperature profiles for the five disks. The blue line is the midplane temperature defined by Equation (C6), and the orange line is a power-law fit to the midplane temperature. Vertical dashed lines mark the position of bright rings in each disk.</figcaption></figure>

<table><caption>**Table 4**
Power-law Fit of the Midplane Dust Temperature</caption><thead><tr><th>Source</th><th><i>T</i><sub>50</sub> (K)</th><th><i>q</i></th></tr></thead><tbody><tr><td>IM Lup</td><td>12.7 ± 0.1</td><td>0.52 ± 0.01</td></tr><tr><td>GM Aur</td><td>22.3 ± 0.4</td><td>0.51 ± 0.02</td></tr><tr><td>AS 209</td><td>19.3 ± 0.2</td><td>0.37 ± 0.01</td></tr><tr><td>HD 163296</td><td>27.7 ± 1.3</td><td>0.46 ± 0.04</td></tr><tr><td>MWC 480</td><td>25.0 ± 0.4</td><td>0.45 ± 0.01</td></tr></tbody></table>

## ORCID iDs

Anibal Sierra https://orcid.org/0000-0002-5991-8073

Laura M. Pérez https://orcid.org/0000-0002-1199-9564

Ke Zhang https://orcid.org/0000-0002-0661-7517

Charles J. Law https://orcid.org/0000-0003-1413-1776

Viviana V. Guzmán https://orcid.org/0000-0003-4784-3040

Chunhua Qi https://orcid.org/0000-0001-8642-1786

Arthur D. Bosman https://orcid.org/0000-0003-4001-3589

Karin I. Öberg https://orcid.org/0000-0001-8798-1347

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Feng Long https://orcid.org/0000-0002-7607-719X

Richard Teague https://orcid.org/0000-0003-1534-5186

Alice S. Booth https://orcid.org/0000-0003-2014-2121

Catherine Walsh https://orcid.org/0000-0001-6078-786X

David J. Wilner https://orcid.org/0000-0003-1526-7587

François Ménard https://orcid.org/0000-0002-1637-7393

Gianni Cataldi https://orcid.org/0000-0002-2700-9676

Ian Czekala https://orcid.org/0000-0002-1483-8811

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Jane Huang https://orcid.org/0000-0001-6947-6072

Jennifer B. Bergner https://orcid.org/0000-0002-8716-0482

John D. Ilee https://orcid.org/0000-0003-1008-1142

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Romane Le Gal https://orcid.org/0000-0003-1837-3772

Ryan A. Loomis https://orcid.org/0000-0002-8932-1219

Takashi Tsukagoshi https://orcid.org/0000-0002-6034-2892

Yao Liu https://orcid.org/0000-0002-7616-666X

Yoshihide Yamato https://orcid.org/0000-0003-4099-6941

Yuri Aikawa https://orcid.org/0000-0003-3283-6884

## References

Alarcón, F., Bosman, A. D., Bergin, E. A., et al. 2021, ApJS, 257, 8

Andrews, S. M. 2020, ARA&A, 58, 483

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJL, 820, L40

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

Bacciotti, F., Girart, J. M., Padovani, M., et al. 2018, ApJL, 865, L12

Bae, J., & Zhu, Z. 2018, ApJ, 859, 118

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924

Birnstiel, T., Dullemond, C. P., & Pinilla, P. 2013, A&A, 550, L8

Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJL, 869, L45

Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

Booth, A. S., Tabone, B., Ilee, J. D., et al. 2021, ApJS, 257, 16

Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJL, 882, L31

Brown, J. M., Blake, G. A., Qi, C., et al. 2009, ApJ, 704, 496

Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71

Casassus, S., van der Plas, G. M., Perez, S., et al. 2013, Natur, 493, 191

Cazzoletti, P., van Dishoeck, E. F., Pinilla, P., et al. 2018, A&A, 619, A161

Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110

Czekala, I., Loomis, R. A., & Teague, R. 2021, ApJS, 257, 2

D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321

Dent, W. R. F., Pinte, C., Cortes, P. C., et al. 2019, MNRAS, 482, L29
Dipierro, G., Pinilla, P., Lodato, G., & Testi, L. 2015a, MNRAS, 451, 974

Dipierro, G., Price, D., Laibe, G., et al. 2015b, MNRAS, 453, L73

Draine, B. T. 2006, ApJ, 636, 1114

Drażkowska, J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H. 2019, ApJ, 885, 91

Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJL, 869, L46

Dullemond, C. P., & Natta, A. 2003, A&A, 408, 161

Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72

Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Fouchet, L., Gonzalez, J. F., & Maddison, S. T. 2010, A&A, 518, A16

Fouchet, T., Maddison, S. T., Gonzalez, J. F., & Murray, J. R. 2007, A&A, 474, 1037

Garufi, A., Benisty, M., Pinilla, P., et al. 2018, A&A, 620, A94

Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Greaves, J. S., & Rice, W. K. M. 2010, MNRAS, 407, 1981

Guilloteau, S., Dutrey, A., Piétu, V., & Boehler, Y. 2011, A&A, 529, A105

Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L48

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Harrison, R. E., Looney, L. W., Stephens, I. W., et al. 2019, ApJL, 877, L2

Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70

Hornbeck, J. B., Swearingen, J. R., Grady, C. A., et al. 2016, ApJ, 829, 65

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018a, ApJL, 869, L42

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020, ApJ, 891, 48

Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018b, ApJL, 869, L43

Hughes, A. M., Andrews, S. M., Espaillat, C., et al. 2009, ApJ, 698, 131

Hughes, A. M., Hull, C. L. H., Wilner, D. J., & Plambeck, R. L. 2013, AJ, 145, 115

Hull, C. L. H., Yang, H., Li, Z.-Y., et al. 2018, ApJ, 860, 82

Hunter, J. D. 2007, CSE, 9, 90

Isella, A., Guidi, G., Testi, L., et al. 2016, PhRvL, 117, 251101

Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L49

Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213

Jennings, J., Booth, R. A., Tazzari, M., Clarke, C. J., & Rosotti, G. P. 2021, arXiv:2103.02392

Jennings, J., Booth, R. A., Tazzari, M., Rosotti, G. P., & Clarke, C. J. 2020, MNRAS, 495, 3209

Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78

Kirchschlager, F., & Bertrang, G. H. M. 2020, A&A, 638, A116

Law, C., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 3

Lin, Z.-Y. D., Lee, C.-F., Li, Z.-Y., Tobin, J. J., & Turner, N. J. 2021, MNRAS, 501, 1316

Lin, Z.-Y. D., Li, Z.-Y., Yang, H., et al. 2020, MNRAS, 496, 169

Liu, H. B. 2019, ApJL, 877, L22

Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A&A, 622, A75

Lommen, D. J. P., van Dishoeck, E. F., Wright, C. M., et al. 2010, A&A, 515, A77

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Lyra, W., & Lin, M.-K. 2013, ApJ, 775, 17

Macías, E., Espaillat, C. C., Ribas, Á, et al. 2018, ApJ, 865, 37

Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33

Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser., 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Miyake, K., & Nakagawa, Y. 1993, Icar, 106, 20

Mori, T., Kataoka, A., Ohashi, S., et al. 2019, ApJ, 883, 16

Mouillet, D., Lagrange, A. M., Augereau, J. C., & Ménard, F. 2001, A&A, 372, L61

Muro-Arena, G. A., Dominik, C., Waters, L. B. F. M., et al. 2018, A&A, 614, A24

Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJL, 748, L22

Najita, J. R., & Kenyon, S. J. 2014, MNRAS, 445, 3315

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Ohashi, S., & Kataoka, A. 2019, ApJ, 886, 103

Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82

Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Sci, 353, 1519

Pérez, L. M., Carpenter, J. M., Chandler, C. J., et al. 2012, ApJL, 760, L17

Pérez, L. M., Isella, A., Carpenter, J. M., & Chandler, C. J. 2014, ApJL, 783, L13

Piétu, V., Dutrey, A., Guilloteau, S., Chapillon, E., & Pety, J. 2006, A&A, 460, L43

Piétu, V., Guilloteau, S., & Dutrey, A. 2005, A&A, 443, 945

Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81

Pinilla, P., Pascucci, I., & Marino, S. 2020, A&A, 635, A105

Pinilla, P., van der Marel, N., Pérez, L. M., et al. 2015, A&A, 584, A16

Pinte, C., Padgett, D. L., Ménard, F., et al. 2008, A&A, 489, 633

Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160

Richard, S., Barge, P., & Le Dizès, S. 2013, A&A, 559, A30

Robert, C. M. T., Méheut, H., & Ménard, F. 2020, A&A, 641, A128

Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173

Ruge, J. P., Flock, M., Wolf, S., et al. 2016, A&A, 590, A17

Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, ApJS, 257, 20

Sierra, A., & Lizano, S. 2020, ApJ, 892, 136

Sierra, A., Lizano, S., & Barge, P. 2017, ApJ, 850, 115

Sierra, A., Lizano, S., Macías, E., et al. 2019, ApJ, 876, 7

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42

Soon, K.-L., Hanawa, T., Muto, T., Tsukagoshi, T., & Momose, M. 2017, PASJ, 69, 34

Surville, C., & Barge, P. 2015, A&A, 579, A100

Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344

Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018, ApJ, 868, L13

Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019, ApJL, 884, L56

Tilling, I., Woitke, P., Meeus, G., et al. 2012, A&A, 538, A20

Toomre, A. 1964, ApJ, 139, 1217

Ueda, T., Kataoka, A., & Tsukagoshi, T. 2020, ApJ, 893, 125

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164

Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJL, 813, L2

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Whipple, F. L. 1972, in Proc. 21st Nobel Symp., From Plasma to Planet, ed. A. Elvius (New York: Wiley-Interscience), 211

Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., & Schwarz, K. R. 2017, NatAs, 1, 0130

Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18