# Modest Dust Settling in the IRAS04302+2247 Class I Protoplanetary Disk

M. Villenave, L. Podio, G. Duchêne, K R Stapelfeldt, C. Melis, C.

Carrasco-Gonzalez, V. Le Gouellec, F. Ménard, M. de Simone, C. Chandler, et al.

► To cite this version:

M. Villenave, L. Podio, G. Duchêne, K R Stapelfeldt, C. Melis, et al. Modest Dust Settling in the IRAS04302+2247 Class I Protoplanetary Disk. The Astrophysical Journal, 2023, 946 (2), pp.70. 10.3847/1538-4357/acb92e. hal-04956210

HAL Id: hal-04956210

https://hal.science/hal-04956210v1

Submitted on 19 Feb 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>

# Modest Dust Settling in the IRAS04302+2247 Class I Protoplanetary Disk

M. Villenave¹, L. Podio², G. Duchêne³,⁴, K. R. Stapelfeldt¹, C. Melis⁵, C. Carrasco-Gonzalez⁶, V. J. M. Le Gouellec⁷, F. Ménard⁴, M. de Simone⁸, C. Chandler⁹, A. Garufi², C. Pinte⁴,¹⁰, E. Bianchi¹¹, and C. Codella²,⁴

¹ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA; marion.f.villenave@jpl.nasa.gov

² INAF—Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy

³ Astronomy Department, University of California, Berkeley, CA 94720, USA

⁴ Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁵ Center for Astrophysics and Space Sciences, University of California, San Diego, La Jolla, CA 92093-0424, USA

⁶ Instituto de Radioastronomía y Astrofísica (IRyA-UNAM), Morelia, Michoacán 58089, Mexico

⁷ SOFIA Science Center, USRA, NASA Ames Research Center, Moffett Field, CA 94045, USA

⁸ European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany

⁹ National Radio Astronomy Observatory, PO Box O, Socorro, NM 87801, USA

¹⁰ School of Physics and Astronomy, Monash University, Clayton Vic 3800, Australia

¹¹ Excellence Cluster ORIGINS, Boltzmannstraße 2, D-85748 Garching bei München, Germany

Received 2022 December 5; revised 2023 January 16; accepted 2023 February 2; published 2023 March 31

## Abstract

We present new Very Large Array observations, between 6.8 and 66 mm, of the edge-on Class I disk IRAS04302+2247. Observations at 6.8 mm and 9.2 mm lead to the detection of thermal emission from the disk, while shallow observations at the other wavelengths are used to correct for emission from other processes. The disk radial brightness profile transitions from broadly extended in previous Atacama Large Millimeter/submillimeter Array 0.9 mm and 2.1 mm observations to much more centrally brightened at 6.8 mm and 9.2 mm, which can be explained by optical depth effects. The radiative transfer modeling of the 0.9 mm, 2.1 mm, and 9.2 mm data suggests that the grains are smaller than 1 cm in the outer regions of the disk, allowing us to obtain the first lower limit for the scale height of grains emitting at millimeter wavelengths in a protoplanetary disk. We find that the millimeter dust scale height is between 1 au and 6 au at a radius 100 au from the central star, while the gas scale height is estimated to be about 7 au, indicating a modest level of settling. The estimated dust height is intermediate between less evolved Class 0 sources, which are found to be vertically thick, and more evolved Class II sources, which show a significant level of settling. This suggests that we are witnessing an intermediate stage of dust settling.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); Radiative transfer (1335); Dust continuum emission (412)*

## 1. Introduction

As the birthplaces of planets, protoplanetary disks are key to understanding the diversity of the observed exoplanet population. Within disks, submicron-sized particles grow to large pebbles sizes, which eventually aggregate to form planetesimals or planetary cores (Drazkowska et al. 2022). This process can be accelerated in high-dust density regions, such as radial substructures or vertically thin dust layers. While radial substructures have been detected in numerous disks (Andrews et al. 2018), fewer measurements of the vertical extent of millimeter dust in protoplanetary disks have been performed (e.g., Pinte et al. 2016; Doi & Kataoka 2021; Villenave et al. 2022).

The vertical thickness of dust in a disk is set by the efficiency of vertical settling. This mechanism allows dust particles to concentrate into the midplane, with concentrations that depend on their interaction with the gas (Weidenschilling 1977). Larger particles are expected to be more decoupled from the gas and more affected by vertical settling. These large dust grains (e.g., millimeter-sized) will thus be more concentrated into the midplane than smaller particles (e.g., micron-sized), which remain well mixed with the gas, and up to high

altitudes (Barrière-Fouchet et al. 2005). In addition, the vertical settling efficiency depends on the turbulence level of the disk and on its evolutionary stage.

Star formation is divided into several classes: Class 0 corresponds to embedded protostars, Class I objects present both a disk and a prominent envelope around the central star, and Class II systems only have the disk left (e.g., Andre et al. 2000). Direct measurements of the millimeter dust scale height in several Class I/II and Class II disks have found that the outer regions are very settled, with a typical scale height of less than 1 au at a radius of 100 au (Pinte et al. 2016; Villenave et al. 2020; Doi & Kataoka 2021; Liu et al. 2022; Villenave et al. 2022). In contrast, observations of younger systems, such as HH212, VLA 1623 West, and L1527, at the Class 0 and 0/I stages, have revealed much thicker disks, which are possibly not affected by vertical settling (Sakai et al. 2017; Lin et al. 2021; Lee et al. 2022; Michel et al. 2022; Ohashi et al. 2022; Sheehan et al. 2022). Further constraints on the evolution of the vertical extent of protoplanetary disks at different evolutionary stages, specifically in the Class I stage, are important for determining the efficiency of this mechanism with age.

In this work, we focus on the IRAS04302+2247 protoplanetary disk system (hereafter, IRAS04302), located in the L1536 cloud in the Taurus star-forming region (d = 161 ± 3 pc; Galli et al. 2019). This disk, classified as Class I by Kenyon & Hartmann (1995), has been the subject of a number of studies at

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<table><caption>**Table 1**
VLA Observation Log Table</caption><thead><tr><th>B</th><th>λ<br>(mm)</th><th>ν<br>(GHz)</th><th>Beam Size<br>(")</th><th>Beam PA<br>(°)</th><th>rms<br>(μJy/beam)</th><th>TToS<br>(min)</th><th>Project</th><th>Observation<br>Date</th><th>Flux<br>Calibrator</th><th>Phase<br>Calibrator</th></tr></thead><tbody><tr><td>Q</td><td>6.8</td><td>40–48</td><td>0.19 × 0.15</td><td>−6</td><td>15.8</td><td>45</td><td>21B-100</td><td>2021/12/18+21</td><td>3C48</td><td>J0431+2037</td></tr><tr><td rowspan="2">Ka</td><td rowspan="2">9.2</td><td>28 – 29</td><td rowspan="2">0.19 × 0.18</td><td rowspan="2">87</td><td rowspan="2">4.3</td><td rowspan="2">510</td><td rowspan="2">21B-183</td><td>2021/09/30</td><td>3C147</td><td>J0426+2327</td></tr><tr><td>+36 – 37</td><td>2021/10/04</td><td>3C286</td><td>J0426+2327</td></tr><tr><td rowspan="2">K</td><td rowspan="2">14</td><td rowspan="2">18 – 23</td><td rowspan="2">0.37 × 0.29</td><td rowspan="2">−14</td><td rowspan="2">8.3</td><td rowspan="2">10</td><td rowspan="2">21B-100</td><td rowspan="2">2021/12/18+21</td><td rowspan="2">3C48</td><td rowspan="2">J0431+2037</td></tr><tr><td>2021/12/18+21</td><td>3C48</td><td>J0431+2037</td></tr><tr><td>C1</td><td>40</td><td>7 – 8</td><td>1.07 × 0.88</td><td>0</td><td>9.4</td><td>8</td><td>21B-100</td><td>2021/12/18+21</td><td>3C48</td><td>J0431+2037</td></tr><tr><td>C2</td><td>66</td><td>4 – 5</td><td>1.70 × 1.42</td><td>−2</td><td>10.6</td><td>8</td><td>21B-100</td><td>2021/12/18+21</td><td>3C48</td><td>J0431+2037</td></tr></tbody></table>

**Note.** “B”: observing band. “λ”: wavelength corresponding to the average frequency representative of the image. “ν”: frequency range of the observations. “TToS”: total time on source.

various wavelengths, providing useful insights into its structure. The scattered light images (Lucas & Roche 1997; Padgett et al. 1999) show a clear dark lane, indicating the presence of a disk that is seen almost edge on (90° ± 3°; Wolf et al. 2003), and a bipolar nebula, which is dominated by a prominent envelope structure. The source has also been observed at millimeter wavelengths, both in continuum and in different molecular lines (Wolf et al. 2008; Podio et al. 2020; van't Hoff et al. 2020; Villenave et al. 2020). The Atacama Large Millimeter/ submillimeter Array (ALMA) observations at 2.1 mm obtained by Villenave et al. (2020) resolved the minor axis size of the disk (beam size 0.′′09 × 0.′′04), showing that the disk is flared—i.e., with a minor axis width increasing with distance from the star. Here, we present new Karl G. Jansky Very Large Array (VLA) observations at wavelengths between 6.8 and 66 mm. These observations are expected to probe grains that are up to one order of magnitude larger than those that have been probed by previous ALMA observations, which are predicted to be more affected by vertical settling. We present the observations in Section 2 and the observational results in Section 3. Then, in Section 4, we compare 0.9 mm, 2.1 mm, and 9.2 mm observations to a grid of radiative transfer models, with the aim of obtaining constraints on the vertical and radial distributions of the largest grains in the disk. The results are discussed in Section 5, and we conclude this work in Section 6.

## 2. Observations and Data Reduction

In this section, we present the new VLA observations of IRAS04302 at 6.8 mm, 9.2 mm, 14 mm, 40 mm, and 66 mm. In addition, we use previously published photometric data points from Gräfe et al. (2013) as well as ALMA 0.9 mm and 2.1 mm observations from Villenave et al. (2020). We refer the reader to Villenave et al. (2020) for details about the calibration of the ALMA observations used in this work.

### 2.1. VLA Ka Data

IRAS04302 was observed in the Ka band (9.2 mm; Project: 21B-183; PI: Podio) in 2021 September and October, with the VLA B configuration, for a total of 8.5 hr on source. In Table 1, we report the observing dates, frequency ranges, and the flux and phase calibrators that were used for these observations.

We calibrated the raw data using the CASA VLA data reduction pipeline, version 6.2.1.7 (CASA Team et al. 2022). Before producing the final images, we identified a mismatch in the flux density of a factor of 1.07 between the observations

taken on September 30/October 21 and those taken on October 4/October 21. The observations made on September 30 and October 21 used the flux calibrator 3C147, which is known to be variable,<fnref n="12" /> contrary to 3C286, which was used for the other two observations in the Ka band. We thus adjusted the fluxes of the observations from September 30 and October 21 to match those of October 4 and 22, using the *genecal* and *applycal* CASA tasks. To maximize the dynamic range of the image, we then performed phase self-calibration on the observations. Finally, we produced the final continuum image using the *tclean* task, with a Briggs weighting (robust 0.5). The resulting continuum beam size and rms are reported in Table 1.

The 9.2 mm image resulting directly from the data is presented in Figure 1. It is likely that the central beam (~0.′′2) includes contributions from processes other than dust thermal emission, such as free-free emission, from ionized jets, or disk winds and gyrosynchrotron emission, from coronal processes (see Melis et al. 2011, and references therein). Throughout the text, we also refer to these processes under the terms “non-dust emission” or “wind/coronal emission.” To subtract such non-dust emission, we use other available VLA data (see Section 2.2), and follow the method previously implemented by Carrasco-González et al. (2019). We describe the methodology in detail in Appendix A. Ultimately, in this paper, when modeling the disk (Section 4), we consider only the regions outside 1 beam size, to limit the effects of the potential variability of processes other than continuum dust emission.

### 2.2. VLA Q, K, and C Data

IRAS04302 was also observed in the Q, K, and C bands (Project: 21B-100; PI: Villenave), corresponding to wavelengths between 6.8 and 66 mm, in 2021 December, with the VLA in the B configuration. Table 1 presents the frequency ranges, times on source, observing dates, and calibrators that were used for these observations. We note that the C-band observations were performed with two widely separated basebands of 1 GHz bandwidth each, which we separated into C1 (40 mm) and C2 (66 mm), to increase the wavelength range probed by these observations. We calibrated the raw data using the CASA pipeline, version 6.2.1.7. To increase the signal-to-noise ratio (S/N) of the emission, we produced the final images using the *tclean* task with a natural weighting. We report the

<footnote n="12">¹² See the VLA manual: science.nrao.edu/facilities/vla/docs/manuals/oss/performance/fdscale.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** VLA observations of IRAS04302. The scale bar represents the flux level; its values range between the 1σ level (rms; Table 1) and the peak flux of each map. This corresponds to brightness ratios (max/1σ) of 7, 15, 34, 71, and 22, respectively, for the 66 mm, 40 mm, 14 mm, 9.2 mm, and 6.8 mm data. We show the beam size (ellipse) in the bottom left corner of each panel.</figcaption></figure>

beam sizes and rms values in Table 1, and show the final images in Figure 1. Similar to the observations at 9.2 mm, the central regions of the 6.8 mm observations are likely affected by processes other than dust thermal emission, which we correct by following the methodology described in Appendix A.

# 3. Results

## 3.1. Continuum Emission, Fluxes, and Sizes

In Figure 1, we present the new 66 mm, 40 mm, 14 mm, 9.2 mm, and 6.8 mm observations, without any correction for coronal/wind emission. We find that the 66 mm, 40 mm, and 14 mm observations are unresolved. The spectral indices obtained for these three wavelengths (see Section 3.2) are inconsistent with dust emission, so we conclude that they are dominated by coronal/wind emission. On the other hand, the 9.2 mm and 6.8 mm observations are strongly centrally peaked, with detection of extended structure along the disk's major axis. The extended component corresponds to thermal emission from the disk, while the central beam is dominated by emission from other processes. In Figure 2, we show the ALMA 0.9 mm and 2.1 mm observations, along with the VLA 6.8 mm and 9.2 mm observations, after correction for processes other than dust continuum emission (see Appendix A). Even after correction, we find that the emission at 6.8 mm and 9.2 mm is centrally peaked (see also Figure 3). This could suggest a higher dust concentration or temperature within the inner region of the disk, although caution is needed, due to the potential variability of free-free and gyrosynchrotron emission.

We present the flux densities of the source at the different wavelengths in Table 2. For the observations where the source is unresolved—namely, the 14 mm, 40 mm, and 66 mm observations—we fit the visibilities by a point source, using the uvmodelfit CASA task. On the other hand, we estimate the flux densities of the extended 6.8 mm and 9.2 mm observations in the images by using an aperture of 0.′′8 × 4′′ (129 × 644 au). We also re-estimate the 2.1 mm and 0.9 mm fluxes from the ALMA observations of Villenave et al. (2020), using the same aperture, which encompasses all the emission at the four bands. For all wavelengths, we evaluate the uncertainty by taking the quadratic sum of the rms (Table 1) and the systematic uncertainty of the observatory (10% for $\lambda \le 13\mu$m and 5% otherwise; see the VLA manual).

We estimate the disk size using the cumulative flux technique along the major axis. Previous literature works

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Side-by-side comparison of the ALMA 0.9 mm and 2.1 mm and the VLA 6.8 mm and 9.2 mm observations, after correction for emission from processes other than dust continuum emission (Appendix A). The contours indicate the emission at the 3σ, 5σ, and 20σ levels of each map. The signal-to-noise (S/N) of each map (max/1σ) is 83, 17, 11, and 34, respectively for the ALMA 0.9 mm and 2.1 mm and the VLA 6.8 mm and 9.2 mm data. Because of the correction for non-dust emission, the S/Ns of the 6.8 mm and 9.2 mm data are decreased, compared to those reported in Figure 1.</figcaption></figure>

(e.g., Ansdell et al. 2018; Long et al. 2018), focusing mostly on low- to mid-inclination systems, generally estimated the cumulative flux using elliptical apertures (related to the disk inclination) of increasing radii. However, IRAS04302 is highly inclined and typically unresolved in the minor axis direction, so using an elliptical aperture to measure the disk size could bias the results. A better approach is to estimate the radius from the cumulative flux that is obtained from the major axis cut directly (see also Rota et al. 2022, who made a similar choice when measuring the gas size of the edge-on disk around HK Tau B). We calculate the cumulative flux at increasingly larger radii along the major axis cut, where the major axis cut corresponds to a 1 pixel line along the position angle (PA) of the disk, passing by the continuum peak. We then estimate the disk radius $R_{95\%}$ as the radius containing 95% of the total flux, and the effective disk radius as that containing 68% of the total flux ($R_{68\%}$). We note that for 6.8 mm and 9.2 mm, we use the maps corrected for coronal/wind emission to estimate the disk size. However, we checked that using the non-corrected maps only led to minimal changes in the apparent size, of order $\lesssim 0.′′1$. To obtain a consistent comparison between the new centimeter-range and the previous millimeter-range observations, we also apply this technique to the previously published ALMA
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Top: major axis profiles of the 0.9 mm, 2.1 mm, 6.8 mm, and 9.2 mm observations (Figure 2), normalized to their peak. The vertical lines illustrate the locations of R<sub>95%</sub> (Table 2). Bottom: major axis profiles of the spectral index maps (Figure 4). The gray regions in both panels indicate the regions that are possibly still impacted by coronal/wind emission. The region is larger in the bottom plot because the 6.8 mm and 9.2 mm data were generated at a larger resolution than the 0.9 mm image.</figcaption></figure>

<table><caption>**Table 2**
Fluxes and Major Axis Sizes of the VLA and ALMA Observations</caption><thead><tr><th>λ (mm)</th><th>Flux (μJy)</th><th>R<sub>68%</sub> (")</th><th>R<sub>95%</sub> (")</th></tr></thead><tbody><tr><td>0.9</td><td>267500 ± 26752</td><td>0.90</td><td>1.52</td></tr><tr><td>2.1</td><td>44500 ± 4450<sup>a</sup></td><td>0.72</td><td>1.51</td></tr><tr><td>6.8</td><td>1850 ± 186<sup>b,c</sup></td><td>0.52</td><td>1.37</td></tr><tr><td>9.2</td><td>820 ± 82<sup>b</sup></td><td>0.39</td><td>1.01</td></tr><tr><td>14</td><td>294 ± 17</td><td>...</td><td>...</td></tr><tr><td>40</td><td>139 ± 11</td><td>...</td><td>...</td></tr><tr><td>66</td><td>62 ± 11</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.** The major axis sizes at 6.8 mm and 9.2 mm were measured on the images that were corrected for emission from processes other than dust continuum emission.

<sup>a</sup> The 2.1mm flux reported here is ~20% higher than the flux reported in Villenave et al. (2020), possibly due to differences in aperture sizes.

<sup>b</sup> After subtracting for non-dust emission, the total fluxes at 6.8 mm and 9.2 mm are 1690 ± 170 μJy and 654 ± 65 μJy, respectively.

<sup>c</sup> The flux density of the low-S/N 6.8 mm data strongly depends on the aperture size, likely due to the presence of large-scale background structure. It can be decreased by a factor of nearly 2, if the aperture size is smaller.

2.1 mm and 0.9 mm observations (Villenave et al. 2020). We report the results in Table 2. We find that both the 68% and 95% flux radii decrease with wavelength. This is also visible in the top panel of Figure 3, which shows the major axis profiles of the 0.9 mm, 2.1 mm, 6.8 mm, and 9.2 mm observations.

## 3.2. Spectral Indices

The spectral index of millimeter–centimeter emission, defined as $\alpha$ in $F_{\nu} \propto \nu^{\alpha}$, can be used to study the dust and optical depth properties within a disk. From our new VLA fluxes at 66 mm, 40 mm, and 14 mm (see Table 2), we estimate

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Centimeter and millimeter spectral index maps.</figcaption></figure>

a spectral index of $\alpha_{14–66 mm} = 0.97 \pm 0.24$. Using the integrated fluxes reported in Table 2 and from Gräfe et al. (2013), we also estimate a spectral index between 0.9 mm and 9.2 mm of $\alpha_{0.9–9.2 mm} = 2.54 \pm 0.06$. At longer wavelengths, the spectral index is dominated by coronal/wind emission, while dust dominates at shorter wavelengths.

We also compute the spectral index maps of the centimeter–millimeter emission of IRAS04302. To do so, we first generate 6.8 mm and 9.2 mm images at the same resolution as the “restored” resolution used by Villenave et al. (2020; 0"'30 × 0"'24; see their Table A.1), with the imsmooth CASA task. We then generate spectral index maps using the immath CASA task on the images at the same resolution. We present the resulting maps in Figure 4 and the major axis cuts in Figure 3.

We find that the spectral indices increase at larger radii and also at longer wavelengths. This behavior is consistent with the expected decrease in optical depth at outer radii. At the center of the disk, the emission is presumably very optically thick, and we obtain a spectral index of 2 at nearly all wavelengths. This suggests that the subtraction of the emission from processes other than dust continuum emission (Appendix A) was reasonable, as wind/coronal emission is generally associated with lower spectral index (e.g., Rodmann et al. 2006). When using the longer wavelengths (9.2 mm), we see a steeper decrease in the spectral index within 0"'5 of the center of the disk. This would also be expected in the case of a concentration of larger particles at the center of the disk. However, given the high optical depths expected at these radii, this cannot be clearly concluded from this simple analysis.

# 4. Radiative Transfer Modeling

## 4.1. Methodology

### 4.1.1. Model Description

To characterize the distribution of the large dust in IRAS04302, we model the new 9.2 mm VLA and previous 0.9 mm and 2.1 mm ALMA observations of the system with the radiative transfer code mcfost (Pinte et al. 2006, 2009). Mcfost solves the temperature structure using Monte Carlo methods that are based on the disk structure and dust properties, properly accounting for scattering effects. We do not consider the 6.8 mm data in this analysis, because they have
<table><caption>**Table 3**
Overview of the Parameter Ranges of the Grid</caption><thead><tr><th>H<sub>ld,100 au</sub></th><th>(au)</th><th>[1, 3, 4, 5, 6]</th></tr></thead><tbody><tr><td>i</td><td>(°)</td><td>[88, 89, 90]</td></tr><tr><td>a<sub>max</sub></td><td>(μm)</td><td>[100, 1000, 10000]</td></tr><tr><td>p</td><td>(–)</td><td>[0.5, 1, 1.5]</td></tr><tr><td>β</td><td>(–)</td><td>[1.06, 1.14, 1.22, 1.3]</td></tr><tr><td>log<sub>10</sub>(M<sub>dust</sub>/M<sub>☉</sub>)</td><td>(–)</td><td>[-5.0, -4.75, -4.5, -4.25, -4.0<br>-3.75, -3.5, -3.25, -3.0, -2.75<br>-2.5, -2.25, -2.0]</td></tr></tbody></table>

**Note.** H<sub>ld,100 au</sub>, i, a<sub>max</sub>, β, and M<sub>dust</sub> are, respectively, the large dust scale height at 100 au, the inclination, the maximum grain size, the surface density exponent, the flaring parameter, and the dust mass in the large grain layer.

a lower S/N and the same angular resolution as the 9.2 mm data (see Table 1).

We assume that the disk is axisymmetric, smooth, and that the surface density follows a power-law distribution: Σ(r) ∝ r⁻$^{p}$, for R<sub>in</sub> < r < R<sub>out</sub>. The vertical extent of the grains is parameterized, such that H(r) = H<sub>100 au</sub>(r/100 au)<sup>β</sup>. We assume that the grain size follows a power-law distribution, such that n(a)da ∝ a⁻³.⁵da.

In addition, as described in Appendix B, we implement a simplified parameterization of vertical settling, with a disk represented by a layer of small grains that are located up to high altitudes in the disk (characterized by H<sub>sd,100 au</sub>) and a layer of larger dust that is more concentrated toward the midplane (characterized by H<sub>ld,100 au</sub>; see, e.g., Villenave et al. 2022). The complete model, described in Appendix B, also includes an envelope to reproduce the thermal part of the spectral energy distribution (SED).

### 4.1.2. Fixed and Varied Parameters

We use a stellar effective temperature of T<sub>eff</sub> = 4500 K and a stellar radius of 3.7R<sub>☉</sub> (Gräfe et al. 2013). We fix the inner disk radius to R<sub>in</sub> = 0.1 au and the outer radius to R<sub>out</sub> = 300 au, based on the apparent size of the 0.9 mm observations. For the layer of big grains, we fix the minimum grain size to a<sub>min</sub> = 10 µm. Following Gräfe et al. (2013), we also assume that the dust grains are composed of a mixture: 62.5% of astronomical silicates and 37.5% of graphite.

We then construct a grid of models by varying several parameters: the dust scale height H<sub>ld,100 au</sub>, the inclination i, the maximum grain size a<sub>max</sub>, the surface density exponent p, the flaring exponent β, and the dust mass M<sub>dust</sub>. We report the range of explored parameters in Table 3. We note that we only consider inclinations between 88° and 90°. These values are chosen based on previous modeling of scattered light observations, which found an inclination of 90° ± 3° (Wolf et al. 2003; Lucas & Roche 1997), and based on the shape of the 2.1 mm image, which implies an inclination >84° (Villenave et al. 2020). Our initial modeling, which tried models between 84° and 88°, produced only very bad fits, and thus we limited our grid to a minimum inclination of 88°. Our grid consists of 5 × 3 × 3 × 3 × 4 × 13 = 7020 models. For each set of parameters, we compute the 0.89 mm, 2.1 mm, and 9.2 mm images, and convolve the resulting images by the ALMA or VLA 2D Gaussian beam.

### 4.1.3. Fitting Procedure

The goal of this modeling is to reproduce the shapes of the millimeter and centimeter emission in order to obtain constraints on the vertical extent of the larger dust particles. Thus, we first build our models by including only a disk region with large grains (a<sub>min</sub> = 10 µm, with varied a<sub>max</sub>; see Table 3). While adding a layer of smaller grains and envelope allows for a more realistic disk temperature structure, calculating such models is significantly more computationally expensive than when only including the layer of big grains. In Appendix B, we check that adding such layers does not significantly affect the shape of the millimeter/centimeter emission, while allowing a good match in the thermal part of the SED between the model and the data. Our modeling is thus based on a layer of big grains only.

Our modeling strategy follows two key steps, as summarized hereafter:

1. For each set of parameters (H<sub>ld,100 au</sub>, i, a<sub>max</sub>, p, and β), we first determine the dust mass M<sub>dust</sub> that offers the best match with the observations along the major axis only. This allows us to obtain optical depths of the best-mass models that are consistent with the observations.

2. We then evaluate the agreement of each best-mass model with the data, by considering both the major and minor axis profiles at the three different wavelengths.

For both steps, we determine the best models by considering the observations at 0.89 mm, 2.1 mm, and 9.2 mm. For each wavelength, we evaluate the model's accuracy in reproducing either the major or minor axis profiles of the data by estimating a χ², defined by:

$$ \chi^2 = \sum \frac{1}{n} \left( \frac{F_d - F_m}{\sigma} \right)^2, \qquad (1) $$

where F<sub>d</sub> and F<sub>m</sub> are the respective fluxes of the data and the model along the major or minor axis, normalized so that their peak intensity is equal to 1, σ is the normalized rms of the data, and n is the number of pixels along the cut. We note that for the 9.2 mm observations, the region within one beam from the central star is not considered, and thus we normalize the data and model major axis profiles to 1 at the distance of one beam from the central star. The accuracy of each model along the major axis is then given by χ<sup>2</sup><sub>maj</sub>, taken as the mean χ<sup>2</sup><sub>maj</sub> between 0.89 mm, 2.1 mm, and 9.2 mm. For each set of parameters (H<sub>ld,100 au</sub>, i, a<sub>max</sub>, p, and β), the first step determines the best-mass model (best M<sub>dust</sub>), which corresponds to that with the lowest χ<sup>2</sup><sub>maj</sub>.

Then, for the second step, we evaluate the agreement of each best-mass model with the data by considering both the major and minor axis profiles at the three different wavelengths. In addition to χ<sup>2</sup><sub>maj</sub>, we evaluate χ<sup>2</sup><sub>min</sub> along the minor axis direction. We generate minor axis profiles by taking the mean of the cuts at all distances along the major axis (see Villenave et al. 2020), excluding the central region for the 9.2 mm comparisons. For each wavelength, we then estimate χ<sup>2</sup><sub>min</sub> using Equation (1) along the minor axis profiles, and obtain χ<sup>2</sup><sub>min</sub> as the average between 0.89 mm, 2.1 mm, and 9.2 mm. Finally, the global agreement between our model and the data used in the second step is obtained by $\phi^2 = \frac{\bar{\chi}_{maj}^2 + \bar{\chi}_{min}^2}{2}$.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Best dust masses $M_{\text{dust}}$ for each combination of ($H_{\text{ld,100 au}}$, $i$, $a_{\text{max}}$, $p$, and $\beta$) in the grid, based on the first modeling step (using $\bar{\chi}_{\text{maj}}^2$; blue bars). For comparison, the red bars show the dust masses of the top 10% of best models, based on $\phi^2$. In our initial grid (i.e., before the first step), there are $5 \times 3 \times 3 \times 3 \times 4 = 540$ models in each mass bin (Table 3).</figcaption></figure>

## 4.2. Results

### 4.2.1. Grid Best Models

We show the distributions of the dust masses for each combination of ($H_{\text{ld,100 au}}$, $i$, $a_{\text{max}}$, $p$, and $\beta$), as obtained during the first step, in Figure 5. The median of the distributions is $M_{\text{dust}} = 3.2 \times 10^{-4} M_{\odot}$, and the 25% and 75% quartiles reach values of $1.8 \times 10^{-4} M_{\odot}$ and $1.8 \times 10^{-3} M_{\odot}$, respectively. For comparison, we also indicate the dust masses of the top 10% of best models (based on $\phi^2$) in this figure. They fall within the range of the best dust masses of the distribution from the first step.

In Figure 5, we can also see that 24 models saturate to the high mass limit of our grid, namely at $M_{\text{dust}} = 10^{-2} M_{\odot}$. These are exclusively models with $H_{\text{ld,100 au}} = 1$ au and $i = 88^{\circ}$, which are typically more centrally peaked than models at higher inclination or with a larger scale height. They require a high dust mass to become sufficiently optically thick to reproduce the observations. In fact, nearly all models with this combination of scale height and inclination are extremely high mass. However, it is very unlikely that the dust mass in IRAS04302 is as high or higher than $M_{\text{d}} = 10^{-2} M_{\odot}$, because, assuming a gas-to-dust ratio of 100, this would imply a star-to-dust mass ratio close to or higher than 1 ($M_{\star} \sim 1.6 M_{\odot}$; Z.-Y. D. Lin et al., in preparation). Thus, this first modeling step indicates that IRAS04302 likely does not have the combination of $H_{\text{ld,100 au}} = 1$ au and $i = 88^{\circ}$.

After performing both fitting steps, we find that the parameters of our best model are $H_{\text{ld,100 au}} = 3$ au, $i = 89^{\circ}$, $a_{\text{max}} = 1$ mm, $\beta = 1.14$, $p = 1.5$, and $M_{\text{dust}} = 3 \times 10^{-4} M_{\odot}$. This model is associated with the following metrics: $\phi^2 = 8.1$, $\chi_{\text{maj}}^2 = 5.5$, and $\chi_{\text{min}}^2 = 10.6$. We show the major and minor axis profiles of this best model in Appendix B, and the maps in Appendix C.

In addition, we also compute the mean and median $\phi^2$, $\chi_{\text{min}}^2$, and $\chi_{\text{maj}}^2$ values for each $H_{\text{ld,100 au}}$, $i$, $a_{\text{max}}$, $p$, and $\beta$ value of the grid, showing them in Figure 6. We find that $H_{\text{ld,100 au}}$ is relatively well constrained, while there are trends for the values of $a_{\text{max}}$ and $i$. Specifically, the models indicate that $1 \text{ au} < H_{\text{ld,100 au}} < 6 \text{ au}$ and suggest that $a_{\text{max}} < 1$ cm and

$i > 88^{\circ}$. For the scale height, we note that $\bar{\chi}_{\text{min}}^2$ and $\bar{\chi}_{\text{maj}}^2$ favor opposite values, such that the constraint that we obtain is a compromise between the two. On the other hand, the profiles for $\beta$ and $p$ appear relatively flat, and these parameters are thus not constrained.

### 4.2.2. Radial Extent and Maximum Grain Size

As illustrated in Figure 6, the best maximum grain size in the disk is mostly determined by the major axis profiles. To better understand the effect of grain size on the major axis profiles at our three wavelengths of interest, we present some models in Figure 7. We choose to show models where only the maximum grain size varies, while the other parameters are fixed to $H_{\text{ld,100 au}} = 3$au, $i = 90^{\circ}$, $p = -1$, and $\beta = 1.14$ (the best values from Figure 6). The corresponding dust masses of each model (as determined in our first modeling step) are $M_{\text{dust}} = 10^{-4} M_{\odot}$ for the models with $a_{\text{max}} = 100$ µm and $a_{\text{max}} = 1$ mm, and $M_{\text{dust}} = 3 \times 10^{-4} M_{\odot}$ for the model with $a_{\text{max}} = 1$ cm.

In Figure 7, we find that the models with $a_{\text{max}} = 100$ µm and $a_{\text{max}} = 1$ mm can simultaneously reproduce the radial extent of the data at 0.9 mm, 2.1 mm, and 9.2 mm, even though radial drift is not included in the modeling. On the contrary, the model with $a_{\text{max}} = 1$ cm is unable to reproduce the brightness of the three bands simultaneously. Because of the presence of the large grains, the radial profile of that model at 9.2 mm is too radially extended to reproduce the data. Most combinations of ($H_{\text{ld,100 au}}$, $i$, $p$, and $\beta$) also exclude a maximum grain size of $a_{\text{max}} = 1$ cm, as can be seen in Figure 6. In Appendix C, we show a similar result, when considering the $\bar{\chi}_{\text{maj}}^2$ values for different pairs of ($a_{\text{max}}$, $p$) and ($a_{\text{max}}$, $i$).

This result indicates that if there are some grains of 1 cm in the disk, they cannot be well mixed with the gas and extend to the outer disk regions, because the emission would then appear too radially extended at 9.2 mm. However, such grains of 1 cm might be more concentrated radially than smaller grains, which cannot be tested with the current data and modeling. Besides, the models show that if the maximum grain size is between $a_{\text{max}} = 100$ µm and $a_{\text{max}} = 1$ mm, optical depths effects alone can be the origin of the apparent radial segregation between millimeter and centimeter observations. This contrasts with the previous radiative transfer modeling results for another edge-on Class I source, CB 26, where a maximum grain size of 5 cm was found throughout the disk (Zhang et al. 2021).

### 4.2.3. Vertical Extent of Large Dust Particles

Our grid models were performed for five different scale heights of the large dust grains. In Figure 8, we present the minor axis profiles of the data and the models, integrated over the major axis extent of the disk (see Section 4.1.3). We represent models with $a_{\text{max}} = 100$ µm, $p = -1$, $\beta = 1.14$, $i = 90^{\circ}$, and varying $H_{\text{ld,100 au}}$. The models with $H_{\text{ld,100 au}} = 1, 3$, and $6$ au have dust masses of $M_{\text{dust}} = 3 \times 10^{-5}, 10^{-4}$, and $3 \times 10^{-4} M_{\odot}$, respectively. Due to the high inclination of the models, these cuts are dominated by the vertical extent of the disk.

The ALMA 2.1 mm image is the most resolved in the vertical direction (Villenave et al. 2020); it is thus the wavelength that is expected to be able to provide the most constraints on the vertical extent of the disk. Indeed, we find that while the 1 au and 3 au models correctly match the marginally resolved and unresolved minor axis profiles of the
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Integrated median $\phi^2$, $\bar{\chi}_{\min}^2$, and $\bar{\chi}_{\max}^2$ for all the best-mass models in our grid. The 40% and 60% percentiles are represented by the shaded contours. Models with $\phi^2 > 25$ are typically bad fits to the data.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Effect of grain size on the major axis profiles, for models with $H_d = 3$ au, $i = 90^\circ$, $a_{max}$, $p = -1$, and $\beta = 1.14$. The shaded contours show three times the rms of the observations. The gray region in the 9.2 mm panel indicates the zone that is potentially contaminated by coronal/wind emission, and not considered in this analysis.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Effect of scale height on the minor axis profiles, for models with $a_{max} = 100$ µm, $p = -1$, $\beta = 1.14$, and $i = 90^\circ$.</figcaption></figure>

ALMA 0.9 mm and VLA 9.2 mm images, they appear to be too vertically thin to reproduce the 2.1 mm images. In contrast, we find that the ALMA 0.9 mm minor axis profile is not well reproduced by the models with a scale height of 6 au, as they appear too vertically extended to match the profiles. As shown in Figure 6, the best compromise between all bands is achieved

for a scale height of 3 au, while models with scale heights of 1 au or 6 au can be excluded, as they are unable to reproduce either the 2.1 mm or the 0.9 mm minor extent. In Appendix C, we show that these constraints on the large dust vertical extent are similar for the different inclinations and flaring exponents in our grid. In particular, we emphasize that while the models
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Top: minor axis size as a function of radius of the 2.1 mm data and models for $\beta = 1.14$ and varying $H_{\text{ld,100 au}}$. Bottom: minor axis size as a function of radius of the 2.1 mm data and models with $H_{\text{ld,100 au}} = 6$ au and varying $\beta$. The shaded regions show the beam size in the direction of the minor axis cuts.</figcaption></figure>

shown in Figure 8 have an inclination of 90°, our grid also explores less inclined disks, allowing us to exclude the possibility of IRAS04302 being both less inclined (88° of inclination) and very settled (see Section 4.2.1).

Up to this point, we have focused on the minor axis cut averaged over the full disk. However, Villenave et al. (2020) have previously identified that the minor axis size of IRAS04302, measured at 2.1 mm, increases with radius. Now, we present the impacts of the different parameters on the variation of the minor axis extent at different distances from the central star. Following Villenave et al. (2020), we fit a Gaussian to the minor axis profiles obtained at different locations. We plot the resulting FWHM in Figure 9 for two sets of models, varying either $H_{\text{ld,100 au}}$ or $\beta$. In the two figures, we fix $i = 90°$, $a_{\text{max}} = 100$ µm, and $p = -1$. The top figure shows the impact of the dust scale height on the minor axis size variation, for $\beta = 1.14$, and the bottom figure shows the effect of the flaring on the shape of the minor axis size variation, for $H_{\text{ld,100 au}} = 6$ au$^{13}$.

We find that while the $H_{\text{ld,100 au}} = 1$ au model does not show any clear variation in the minor axis size with respect to the distance to the star, every other model does. The model with a dust scale height of 6 au provides the best match to the profile and absolute value. On the other hand, the bottom part of Figure 9 shows that the impact of the flaring exponent is less important than the impact of the scale height. For the range of

flaring values explored within our grid, and a dust scale height of 6 au, all models are consistent with the data within the uncertainties, which does not permit the constraining of the flaring of the large dust particles in this system.

To summarize, using a grid of radiative transfer models, we are able to constrain the scale height of the large grains in the Class I disk IRAS04302. We find that it is of a few astronomical units at a radius of 100 au, where values lower than 1 au and higher than 6 au are excluded from our analysis. We obtain an upper limit that is consistent with the results of Z.-Y. D. Lin et al. 2023 (in preparation), who model high-angular resolution 1.3 mm observations of IRAS04302, and that is also consistent with the previous modeling results from Gräfe et al. (2013), using vertically unresolved millimeter data. However, the flaring exponent $\beta$ of the millimeter disk could not be constrained, as its effects are minimal on the minor axis profile radial variation.

We note that no model is simultaneously able to reproduce the 0.9 mm and 2.1 mm data (see Figure 8). We suggest that this might be due to the simplifying assumptions that we made for our modeling—namely, that we consider fully mixed grains without substructures, that we use a unique maximum grain size, and/or that we determine the dust mass by simultaneously reproducing the major axis profiles of the 0.9 mm, 2.1 mm, and 9.2 mm observations. In particular, our results might suggest that the optical depth change between 0.9 mm and 2.1 mm is less than what is assumed in the models, such that the 0.9 mm image could appear less vertically extended. Besides, the 0.9 mm image is also only marginally resolved in the vertical direction, such that further observations at higher angular resolution would allow the further constraining of the vertical extent of the large dust in IRAS04302. We note that Z.-Y. D. Lin et al. 2023 (in preparation) estimate a dust scale height of about 6 au at 100 au, when working on independent observations of IRAS04302 at 1.3 mm with high angular resolution. Their constraints indeed appear to be more consistent with our results based on our most resolved observation (2.1 mm), rather than the results based on the less resolved 0.9 mm image.

# 5. Discussion

## 5.1. Dust and Gas Vertical Extent in IRAS04302

In Section 4.2.3, we find that the large dust scale height is constrained within $1\text{au} < H_{\text{ld,100 au}} < 6\text{au}$. This result can be compared to the gas scale height that we estimate in Appendix B. We use the midplane temperature of our best model, to which we add an envelope and a disk layer of small grains, to obtain a better representation of the disk temperature structure. We obtain a midplane temperature of 21 K at 100 au from the central star,$^{14}$ which translates into a gas scale height of $H_{\text{g,100au}} \sim 7$ au for a stellar mass of $M_* = 1.6 M_{\odot}$ (Z.-Y. D. Lin et al. 2023, in preparation). The gas scale height appears to be slightly larger than our upper limit for the large dust scale height, suggesting that at least a modest amount of vertical settling occurs in the Class I disk IRAS04302, in agreement with the previous conclusions of Gräfe et al. (2013).

From the difference in the scale height between the gas and the millimeter dust grains, we can also estimate the degree of coupling of dust with gas, parameterized by the $\alpha/\text{St}$ ratio,

<footnote n="13">13 The models with varying flaring have dust masses of $M_{\text{dust}} = 3 \times 10^{-4}$ for $\beta = 1.06$ and 1.14 and of $M_{\text{dust}} = 3 \times 10^{-4} M_{\odot}$ for $\beta = 1.22$ and 1.3, respectively. Those with varying scale heights are the same as in Figure 8, and the model with $H_{\text{ld,100 au}} = 4$ au has a dust mass of $M_{\text{dust}} = 2 \times 10^{-5} M_{\odot}$.</footnote>

<footnote n="14">14 This temperature is in perfect agreement with the location of the CO snowline from van't Hoff et al. (2020), based on the disappearance of CO from the midplane at 100 au (see also Podio et al. 2020).</footnote>
<table><caption>**Table 4**
Stellar Mass, Outer Radius, Bolometric Luminosity, Accretion Rate, Evolutionary Stage, and [α/St]$_{100au}$ for Seven Disks with Constraints on Their Large Dust Height from ALMA Observations</caption><thead><tr><th></th><th><i>M</i><sub>*</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>R</i><sub>out,mm</sub><br>(au)</th><th><i>L</i><sub>bol</sub><br>(<i>L</i><sub>☉</sub>)</th><th>Accretion<br>(<i>M</i><sub>☉</sub> yr<sup>-1</sup>)</th><th>Class</th><th>[α/St]<sub>100 au</sub></th><th>Dust Height</th><th>References</th></tr></thead><tbody><tr><td>HH212</td><td>0.2–0.3</td><td>60</td><td>9.0</td><td>2–4 × 10<sup>-6</sup></td><td>0</td><td>...</td><td>Vertically thick</td><td>1, 2, 3, 4</td></tr><tr><td>L1527</td><td>0.2–0.4</td><td>54</td><td>2.75</td><td>5–11 × 10<sup>-7</sup></td><td>0/I</td><td>...</td><td>Vertically thick</td><td>5, 6, 7, 8</td></tr><tr><td>VLA 1623 W</td><td>0.45</td><td>50</td><td>0.13</td><td>2.5 × 10<sup>-8</sup></td><td>0/I</td><td>...</td><td>Vertically thick</td><td>9, 10, 11, 12</td></tr><tr><td>IRAS04302</td><td>1.6</td><td>300</td><td>0.15–0.67</td><td>0.7–3 × 10<sup>-8</sup></td><td>I</td><td>>2 × 10<sup>-2</sup></td><td>1 au &lt; <i>H</i><sub>ld,100 au</sub> &lt; 6 au</td><td>This work, 13, 14</td></tr><tr><td>HL Tau</td><td>1.7</td><td>150</td><td>3.5–15</td><td>8.7 × 10<sup>-8</sup></td><td>I/II</td><td>...</td><td><i>H</i><sub>ld,100 au</sub> &lt; 1 au</td><td>14, 15, 16</td></tr><tr><td>HD163296</td><td>2.0</td><td>240</td><td>...</td><td>4.5 × 10<sup>-7</sup></td><td>II</td><td>&lt;1 × 10<sup>-2</sup></td><td><i>H</i><sub>ld,100 au</sub> &lt; 1 au</td><td>17, 18, 19, 20, 21</td></tr><tr><td>Oph163131</td><td>1.2</td><td>150</td><td>...</td><td>...</td><td>II</td><td>&lt;5 × 10<sup>-3</sup></td><td><i>H</i><sub>ld,100 au</sub> &lt; 0.5 au</td><td>22, 23</td></tr></tbody></table>

**Note.** For HH212, L1527, VLA 1623 W, and IRAS04302, we recalculate the accretion rates based on their bolometric luminosity (see the main text). For HD163296 and HL Tau, we report estimates based on the UV excess or HI recombination line luminosity. The outer radii and dust heights correspond to the radial extents or the dust heights inferred for the dust detected at millimeter wavelengths with ALMA.

**References:** (1) Lee et al. (2014); (2) Codella et al. (2014); (3) Lee et al. (2017); (4) Lin et al. (2021); (5) Maret et al. (2020); (6) Tobin et al. (2008); (7) Sheehan et al. (2022); (8) Ohashi et al. (2022); (9) S. Mercimek et al. 2023, in preparation; (10) Harris et al. (2018); (11) Murillo et al. (2018); (12) Michel et al. (2022); (13) Z.-Y. D. Lin et al. 2023, in preparation; (14) Robitaille et al. (2007); (15) Pinte et al. (2016); (16) Beck et al. (2010); (17) Teague et al. (2019); (18) Muro-Arena et al. (2018); (19) Mendigutía et al. (2013); (20) Doi & Kataoka (2021); (21) Liu et al. (2022); (22) Flores et al. (2021); (23) Villenave et al. (2022).

where St is the Stokes number and α is the turbulence parameter. Indeed, assuming that turbulence balances gravity, for grains in the Epstein regime (St ≪ 1), and under the assumption of z ≪ H_g, the dust scale height can be written as:

$$H_d = H_g \left( 1 + \frac{\text{StSc}}{\alpha} \right)^{-1/2}, \qquad (2)$$

with Sc being the Schmidt number, the ratio between turbulent viscosity and turbulent diffusivity. Following previous studies (Dullemond et al. 2018; Rosotti et al. 2020; Villenave et al. 2022), we assume Sc = 1, which is valid for particles with St ≪ 1 (Youdin & Lithwick 2007; see also Johansen & Klahr 2005). For H<sub>ld,100 au</sub> > 1 au and H_g ~ 7 au, we find that [α/St]<sub>100au</sub> > 2 × 10<sup>-2</sup>. This value is higher than previous estimates for Class II systems (Table 4), which we discuss in Section 5.2.

## 5.2. Comparison with Other Systems

We now aim to look for any trends in the variation of the settling efficiency with different parameters, such as the stellar mass, disk outer radius, and time. In Table 4, we compile various parameters for seven disks at different evolutionary stages, with constraints on the large dust vertical height (HH212, L1527, VLA 1623 W, IRAS04302, HL Tau, HD163296, and Oph163131; Pinte et al. 2016; Nakatani et al. 2020; Doi & Kataoka 2021; Lin et al. 2021; Liu et al. 2022; Michel et al. 2022; Ohashi et al. 2022; Sheehan et al. 2022; Villenave et al. 2022).

The stellar masses that are gathered in Table 4 are all dynamical estimates, based on the Keplerian rotation of the gas. For the Class I and II disks, the abundant <sup>12</sup>CO lines are used, while rarer species such as C<sup>17</sup>O, HCO<sup>+</sup>, or SO are considered for the more embedded Class 0 protostars (see the references in Table 4). We also recalculate the mass accretion rates of the least evolved systems (HH212, L1527, VLA 1623 W, and IRAS04302) based on their bolometric luminosity, using M<sub>acc</sub> = (L<sub>bol</sub>R<sub>*</sub>)/(GM<sub>*</sub>(1 - R<sub>*</sub>/R<sub>in</sub>)) (Gullbring et al. 1998), assuming that the bolometric luminosity is dominated by the accretion luminosity. We assume a stellar radius of R<sub>*</sub> = 2R<sub>☉</sub> (Lee 2020) and the inner radius of the accretion disk as

$R_{in} = 5R_{*}$ (Gullbring et al. 1998). For HD163296 and HL Tau, we report estimates based on the UV excess or HI recombination line luminosity, and we note that no clear signs of accretion are found in the spectra of Oph163131 (Flores et al. 2021). Finally, the last column of Table 4 reports the dust heights of the different systems. For the most evolved sources, specific radiative transfer modeling were performed, allowing the authors to quantify the vertical extents of the large dust particles and to compare them with those of the gas (Pinte et al. 2016; Doi & Kataoka 2021; Wolff et al. 2021; Liu et al. 2022; Villenave et al. 2022). On the other hand, for the younger Class 0 systems, determining the gas height is complex, due to the presence of a dense envelope, and the different papers do not provide quantitative measurements of the gas and dust scale heights. Nevertheless, by using radiative transfer or geometrical analysis in the uv-plane, or by measuring the apparent size in the image plane, the authors have argued that the young disks of HH212, L1527, and VLA 1623 W were geometrically thick, and possibly not affected by vertical settling (Lin et al. 2021; Michel et al. 2022; Ohashi et al. 2022; Sheehan et al. 2022).

From Table 4, we find that the millimeter dust height appears to correlate with the evolutionary class. In particular, the Class I IRAS04302 appears to be less settled than the Class II and I/II disks (HD163296, Oph163131, and HL Tau). In addition, while we have identified at least a moderate level of settling in IRAS04302, previous studies of younger Class 0 and 0/I systems (HH212, L1527, and VLA1623W) have suggested that settling did not occur at all in these young sources. Our results thus seem to support the evolution of the settling efficiency with time.

In Table 4, we also find that the younger sources have the least massive stars and the least extended disks, which might impact the settling efficiency. To investigate the differences in the dust height between the disks, we estimate the timescale for settling under simple assumptions. When only the thermal speed of molecules and the vertical settling are considered—in other words, in the absence of vertical turbulence and infalling material—the timescale for settling is given by Armitage (2015):

$$t_{\text{settle}} = \frac{\rho}{\rho_m} \frac{v_{\text{th}}}{a} \frac{1}{\Omega_K^2} \propto \rho \frac{\sqrt{T}}{a} \frac{R^3}{M_*}, \qquad (3)$$
where $\rho_m$ is the dust material density, $\rho$ is the gas density, $a$ is the particle size, $v_{th} = \sqrt{8k_B T/\pi\mu m_H}$ is the thermal speed of the molecules, and $\Omega_K = \sqrt{GM_*/R^3}$ is the disk rotation frequency, which depends on the stellar mass ($M_*$) and the local radius ($R$). The settling timescale strongly depends on the radius, particle size, stellar mass, and gas density. Specifically, settling is expected to be faster for larger particles, larger stellar mass, and at the inner regions of the disk. If they were in the same evolutionary stage, smaller disks—such as HH212, L1527, and VLA 1623 W—should thus appear more settled than larger disks—such as IRAS04302, HL Tau, HD163296, and Oph163131—which we do not observe.

Using the generic numerical values given by Sheehan et al. (2022)—$\rho = 10^{-12} \text{ g cm}^{-3}$, $\rho_m = 3 \text{ g cm}^{-3}$, and $T = 50 \text{ K}$—we find that the settling timescale for grains of 1 mm is $t_{settle,100\text{au}} \sim 0.4-0.09 \text{ Myr}$ at 100 au, and $t_{settle,50\text{au}} \sim 0.05-0.01 \text{ Myr}$ at 50 au, for $M_* = 0.45-2M_\odot$. At 100 au from the star, these estimates are comparable to the lifetimes of Class 0 systems (0.13–0.26 Myr; Dunham et al. 2015; see also Kristensen & Dunham 2018). This suggests that, if they are located at 100 au from the star, 1 mm grains do not have time to settle at the Class 0 stage, while at 50 au they are expected to settle before the end of the Class 0 stage. On the other hand, Class I objects have longer lifetimes (0.27–0.52 Myr; Dunham et al. 2015), and millimeter-sized particles, if present, should already have totally settled at both 50 au and 100 au from the star.

Yet, we identify only a modest level of settling in IRAS04302, the only Class I disk in Table 4, and there is no evidence for settling in the small Class 0 and 0/I sources HH212, L1527, and VLA 1623 W. This suggests that the turbulence could be stronger in the earliest stages of star formation, preventing the dust grains from settling within the predicted timescale without turbulence. Interestingly, the values of $\alpha/St$ presented in Table 4 also indicate that, at 100 au from the star, turbulence is stronger in the Class I disk IRAS04302 than in the Class II disks Oph163131 and HD163296. In addition, more evolved systems tend to have lower accretion rates (see Table 4), which would also be consistent with the evolution of turbulence with time, even though the accretion rate does not necessarily relate to the turbulence at 100 au.

Alternatively, the difference in dust heights between the evolutionary classes might also be related to the differences in the environments external to the disk. Specifically, with an infalling envelope at the earliest stages of star formation, fresh grains are constantly being resupplied to the disk, both radially and vertically. So while the first grains to fall in may have had time to settle in the disk, recently added grains would not. Yet, recent studies have shown that the grains falling onto the disk that are coming from the infalling envelope cannot have grown larger than a few tens of microns, for timescale and density reasons (e.g., Silsbee et al. 2022), such that it is not clear whether this effect would be sufficient to explain the thickness of young disks when they are observed at millimeter wavelengths. Further studies are needed in order to determine which mechanism is dominant in setting the apparent heights of disks at millimeter wavelengths for different evolutionary stages.

## 6. Conclusions

We present new VLA observations of the Class I disk IRAS04302, at five different wavelengths between 6.8 mm and 66 mm. The disk is detected at both 6.8 mm and 9.2 mm, while processes other than dust continuum emission dominate the observations at 14 mm, 40 mm, and 66 mm. We compare the disk size at centimeter wavelengths with its extent as seen in the millimeter ALMA observations at 0.9 mm and 2.1 mm, and we find that the disk is about 70% smaller at 9.2 mm than at 0.9 mm. The spectral index maps that are obtained with the different millimeter and centimeter observations allow us to identify an increase in the spectral index with radius, from ~2 in the center to ~3 in the outer disk, consistent with the size differences observed in the images. We also find that the spectral index obtained when including the centimeter observations is consistently higher than the spectral index obtained using only the 0.9 mm and 2.1 mm images (by about 0.5 dex, except at the center of the disk), indicating that the observations at 9.2 mm and 6.8 mm are more optically thin than the ALMA observations.

With the goal of characterizing the distribution of the large dust in IRAS04302, we produced a grid of radiative transfer models, aiming to reproduce the 0.9 mm, 2.1 mm, and 9.2 mm data. The models implement one layer of large grains without substructures, and we determine the dust mass for each combination of parameters, based on the major axis profiles at the three wavelengths previously mentioned. We show in Appendix B that adding a small grain disk layer and envelope does not affect the shape of the millimeter/centimeter emission, but does allow the recovery of a more realistic temperature structure.

Our results indicate that if there are some grains of 1 cm in the disk, they cannot be well mixed with smaller dust particles and extend to the outer disk regions, because the emission would then appear too radially extended at 9.2 mm. However, such grains of 1 cm might be more concentrated radially than smaller grains, which cannot be tested with the current data and modeling. Moreover, models with maximum grain sizes of 100 µm or 1 mm are able to reproduce the 0.9 mm, 2.1 mm, and 9.2 mm major axis profiles simultaneously, suggesting that the apparent size difference between the millimeter and centimeter observations could mostly be due to optical depth effects. The models also favor a disk inclination that is very close to edge on ($i > 88^\circ$). However, no clear trends were found for the different values of the surface density exponent $p$ and the flaring exponent $\beta$ as tested within our grid.

Our modeling strategy also allows us to determine a plausible interval for the large dust scale height of 1 au $< H_{\text{ld,100 au}} < 6$ au. When compared to our estimate of the gas scale height $H_{\text{g,100au}} \sim 7$ au, based on the midplane temperature profile of our best model, we identify that the millimeter dust in the disk is subject to at least moderate vertical settling. The level of settling in this Class I disk contrasts with previous results for Class 0 systems, where no settling occurs, and for Class II disks, which are found to be thinner. By estimating the timescale of the vertical settling for millimeter-sized particles and comparing it to the lifetimes of Class 0 and Class I systems, we suggest that this variation of the settling efficiency with time is linked to some variation of the turbulence or of external conditions with time.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Resulting 6.8 mm and 9.2 mm images before (left) and after (middle) correction for non-disk emission, at the same angular and brightness scale. The right panel shows the major axis cut to the same brightness scale, illustrating that this correction only affects the central beam of the emission (marked by the filled gray region).</figcaption></figure>

The research of M.V. was supported by an appointment to the NASA Postdoctoral Program at the NASA Jet Propulsion Laboratory, administered by Oak Ridge Associated Universities, under contract with NASA. E.B. acknowledges the Deutsche Forschungsgemeinschaft (DFG; the German Research Foundation) under Germany's Excellence Strategy —EXC 2094—390783311. C.C.-G. acknowledges support by UNAM DGAPA-PAPIIT grant IG101321 and CONACYT Ciencia de Frontera grant No. 86372. G.D. acknowledges support from NASA grant 80NSSC18K0442 as well as NNX15AC89G and NNX15AD95G/NExSS. C.C. and L.P. acknowledge the EC H2020 research and innovation program for the project "Astro-Chemical Origins" (ACO; No. 811312) and the PRIN-MUR 2020 MUR BEYOND-2p (Astrochemistry beyond the Second Period Elements; Prot. 2020AFB3FX). This paper makes use of the following ALMA data: ADS/JAO. ALMA\#2016.1.00460.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.

**Facilities:** ALMA, VLA.

**Software:** CASA (CASA Team et al. 2022), `mcfost` (Pinte et al. 2006, 2009), `Matplotlib` (Hunter 2007), `Numpy` (Harris et al. 2020).

# Appendix A

## Correction for Processes Other than Dust Emission

The variety of the VLA observations allows us to perform a detailed analysis to remove emission from processes other than dust thermal emission (also referred to as non-dust or wind/coronal emission throughout the text) at the center of the disk. To do so, we follow the approach highlighted in Carrasco-González et al. (2019). We first align the observations at the different wavelengths, using the CASA task `fixplanets`. We then produce an image combining the 66 mm, 40 mm, 14 mm, and 6.8 mm data together, aiming to model the non-dust emission of the disk. We use the `tclean` task with `nterm = 2`, which allows us to obtain both the flux and the spectral index of the emission. We also assume a Briggs robust parameter of 0 in order to be sensitive only to the brightest emission of the maps.

We note that at 6.8 mm and 9.2 mm, it is possible that both dust and non-dust emission are present in the inner regions. To determine whether to use the 6.8 mm or 9.2 mm observations to model the wind/coronal emission, we look at the shape of the emission when imaged with a robust Briggs parameter of 0. First, we note that this parameter value is chosen to decrease the sensitivity of the imaging, such that only the brightest emission (expected to be from processes other than dust thermal emission) will be visible in the final image. In addition, we expect the wind/coronal emission to be concentrated within one beam, while the dust emission from the disk should be more extended. We find that by using a Briggs parameter of 0, the 6.8 mm emission is concentrated within a central point source, and is thus likely dominated by processes other than dust thermal emission. On the contrary, with a Briggs parameter of 0, the 9.2 mm image is extended, meaning that it still contains substantial dust emission. Thus, we include the 6.8 mm but not the 9.2 mm observations in the non-dust model.

From the combined 66 mm, 40 mm, 14 mm, and 6.8 mm image, we find that the wind/coronal emission is unresolved, as expected. We obtain a model consisting of only one point, located at the center of the disk, with a flux density of 150 µJy at 26 GHz and a spectral index of +0.46, consistent with moderately optically thick free-free emission (Rodmann et al. 2006). This allows us to predict the contributions from processes other than dust thermal emission at 9.2 mm and 6.8 mm, assuming that such emission did not vary significantly between 2021 October and December.

To remove coronal/wind contamination from the 6.8 mm and 9.2 mm images, we then insert the non-dust model images (the combined 66 mm, 40 mm, 14 mm, and 6.8 mm image) into the "model" column of the 6.8 mm and 9.2 mm visibilities, using the `ft` task, and subtract the model column from the data column, using `uvsub`. Finally, we produce a model image using the parameters presented in Section 2.1. We present the 6.8 mm and 9.2 mm images and the major axis cut before and after correction in Figure 10. These images clearly show that the correction significantly reduces the centrally peaked feature within one beam of the center, while keeping the outer regions unaffected. When modeling the disk (Section 4), we use the corrected 9.2 mm image, but consider only the regions beyond one beam size, to limit the effect of the potential variability of the wind/coronal emission.
<table><caption>**Table 5**
Small Grain Layer and Envelope Parameters</caption><thead><tr><th></th><th></th><th>Envelope and Cavity</th></tr></thead><tbody><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>0.005-0.5</td></tr><tr><td>R<sub>in</sub> - R<sub>out</sub></td><td>(au)</td><td>0.1-1000</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>1 · 10<sup>-5</sup></td></tr><tr><td>γ</td><td>(-)</td><td>-2</td></tr><tr><td>z<sub>0</sub>, r<sub>0</sub>, β<sub>env</sub></td><td>(au, au, -)</td><td>5, 1, 1</td></tr></tbody><thead><tr><th></th><th></th><th>Small Grain Disk Layer</th></tr></thead><tbody><tr><td>a<sub>min</sub> - a<sub>max</sub></td><td>(μm)</td><td>0.005-10</td></tr><tr><td>Rin - Rout</td><td>(au)</td><td>0.1-300</td></tr><tr><td>M<sub>dust</sub></td><td>(M<sub>☉</sub>)</td><td>5 · 10<sup>-7</sup></td></tr><tr><td>H<sub>sd,100au</sub>, β, p</td><td>(au, -, -)</td><td>6.7, 1.14, -1</td></tr></tbody></table>

**Note.** The parameters for the small grain and envelope regions allow the production of a more realistic temperature structure for the disk than in Section 4, but we emphasize that they are not expected to reproduce the previous scattered light observations and should not be considered as constraints on the envelope mass and shape.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Comparison of the SED of the best model from Section 4.2.1 (orange), that model plus an envelope and small grain disk layer (Table 5; red), and the data (blue). The 6.8 mm and 9.2 mm VLA fluxes are corrected from processes other than dust emission.</figcaption></figure>

## Appendix B
Simultaneous Match of the Thermal Part of the SED and ALMA/VLA Maps

Previous observations and modeling (Lucas & Roche 1997; Padgett et al. 1999; Wolf et al. 2003) have found that the scattered light is dominated by an envelope. In Section 4, however, we have omitted this part, in order to save some computational power in the generation of our grid. In this Appendix, we show that including a disk layer of small grains and an envelope does not significantly affect the shapes of the

0.89 mm, 2.1 mm, and 9.2 mm observations. Our goal here is to produce a more realistic temperature structure of the disk, which allows the matching of the millimeter fluxes and the thermal part of the SED ($\lambda \gtrsim 20 \mu m$).

In order to obtain a better temperature structure of our disk, we thus add a disk layer of small grains and an envelope to the layer of large grains obtained as best model from our fitting methodology. We first add a layer of small grains with the parameters indicated in Table 5. We then test a small number of envelope parameters. For simplicity, we assume that the envelope is spherically symmetric and that the density follows

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Comparison of the major axis profiles (top) and minor axis profiles (bottom) of the best model from Section 4.2.1 (orange), that model plus an envelope and small grain layer (Table 5; red), and the data (blue).</figcaption></figure>
a power law in radius: $\rho \propto r^{\gamma}$. Following Wolf et al. (2003), we also include a cavity in our modeling. We assume that no dust is present above a surface represented by $z = z_0(r/r_0)^{\beta_{\text{env}}}$. We set the inner radius to 0.1 au, the outer radius to 1000 au, and vary the following parameters: $a_{\text{max,env}} \in [0.5, 1]$, $\gamma \in [-1, -2]$, $z_0 \in [1, 3, 5]$, and $M_{\text{dust,env}} \in [10^{-5}, 10^{-4}, 10^{-3}]$. Finally, we determine the best set of parameters based on a $\chi^2_{\text{SED}}$ estimate, using the SED match between 20 and 10,000 µm (seven points) and $\phi^2$ of the 0.9 mm observations only.

In Table 5, we show the resulting parameters for the small grain and envelope layers. We also include comparisons with the best big-grain model ($H_d = 3$au, $i = 89^{\circ}$, $a_{\text{max}} = 1000$ µm, $\beta = 1.14$, $p = -1.5$, and $M_d = 0.0003 M_{\odot}$), with or without including small grains and an envelope, in Figure 11 for the major and minor axis cut, and in Figure 12 for the SED. As expected, we find that the thermal emission part of the model SED reproduces the data significantly better when small grain and envelope layers are included. On the other hand, no significant differences are found in the shapes of the millimeter/centimeter emission between the two models. Most of the differences lie on the outer edges of the major axis profile at 0.9 mm. This shows that while the small grain and envelope layers are mostly invisible in the shapes of the millimeter/centimeter observations, they contribute to the obtaining of a more realistic disk structure and a larger total millimeter/centimeter flux, by allowing the midplane to be warmer. We emphasize that the parameters for the small grain and envelope regions are not expected to reproduce the previous scattered light observations and should not be considered as constraints on the envelope mass and shape.

We find that adding the small grain and envelope layers allows a significant increase in the midplane temperature. At 100 au from the central star, it goes from 11 K, with large grains only, to 21 K, with the inclusion of the other layers. Using this midplane temperature value, it is possible to infer the gas scale height, assuming that it is at the hydrostatic equilibrium, following

$$H(r) = \sqrt{\frac{k_B T(r) r^3}{G M_* \mu m_p}}, \quad (B1)$$

where $k_B$ is the Boltzmann constant, $\mu = 2.3$ is the reduced mass, $m_p$ is the proton mass, and $G$ is the gravitational constant. Using a stellar mass of $M_* = 1.6 M_{\odot}$ (Z.-Y. D. Lin et al. 2023, in preparation), we find that the gas scale height of this disk is $H_{g,100au} \sim 7$ au.

## Appendix C
### Additional Model Comparisons

#### C.1. Residual Maps of the Best Model

We present the residual maps of the best model in Figure 13. The parameters for the best model are: $H_{ld,100au} = 3$ au, $i = 89^{\circ}$, $a_{\text{max}} = 1$ mm, $\beta = 1.14$, $p = 1.5$, and $M_{\text{dust}} = 3 \times 10^{-4} M_{\odot}$. We also note that the maps shown here do not have the small grain disk and envelope layers discussed in Appendix B.

#### C.2. Effects of p and i on $a_{\text{max}}$

While the maximum grain size has a clear impact on the shapes of the major axis profiles (Section 4.2.2; Figure 7), it is possibly degenerate with three other parameters: $M_{\text{dust}}$, i, and p. In our first fitting step, we determine the best mass that allows

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Left: data images. Middle: images of the best model. Right: residual maps, with deviations by more than 5σ indicated by the black contours. The central mask on the 9.2 mm maps indicates the region that is potentially affected by coronal/wind emission but not considered in our modeling.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Median $\bar{\chi}_{\text{maj}}^2$ values for each combination of ($a_{\text{max}}$, p) and ($a_{\text{max}}$, i). The darker colors indicate smaller $\bar{\chi}_{\text{maj}}^2$, and thus a better combination of parameters.</figcaption></figure>

to reproduce the major axis profile, thus the possible degeneracies between $M_{\text{dust}}$ and $a_{\text{max}}$ are already taken into account. In Figure 14, we show the detailed $\bar{\chi}_{\text{maj}}^2$ median matrix for each combination of ($a_{\text{max}}$, i) and ($a_{\text{max}}$, p). Both matrices show that, independently of the values of i and p, $a_{\text{max}} = 1$ cm (the lowest line) leads to significantly higher median $\bar{\chi}_{\text{maj}}^2$ values than lower maximum grain sizes. As indicated previously, this implies that the maximum grain size in the system cannot be 1 cm in the outer regions of the disk, and must be smaller.

For all $a_{\text{max}}$, it is also clear that higher inclinations and higher surface density exponents are preferred ($i > 88^{\circ}$ and
<figure><img src="image_18.png" /><figcaption>Figure 15. Effect of *p* on the major axis profiles, for *a*<sub>max</sub> = 100µm (top) and *a*<sub>max</sub> = 1 cm (bottom).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 16. Median $\bar{\phi}^2$ for each combination of ($H_d, i$) and ($H_d, \beta$). The darker colors indicate smaller $\bar{\phi}^2$, and thus a better combination of parameters.</figcaption></figure>

$p > 0.5$). This is because inclinations of 88° lead to models that are too centrally peaked. On the contrary, a surface density exponent of $p = 0.5$ implies a shallower radial profile, leading the models to be too bright and radially extended at 9.2 mm compared to the data. In Figure 15, we show the effect of the surface density exponent on the shapes of the major axis profiles, for grain sizes of $a_{max} = 100$ µm and $a_{max} = 1$ cm. For the lower maximum grain size ($a_{max} = 100$ µm), low values of surface density exponents are more strongly excluded than when $a_{max} = 1$ cm.

### C.3. Effects of $\beta$ and $i$ on $H_{ld,100\,au}$

The models presented in Figure 8 were computed for an inclination of $i = 90°$ and flaring of $\beta = 1.14$. However, both

the inclination and the flaring exponent $i$ and $\beta$ can also impact the apparent minor axis size of the disk. In Figure 16, we present the median $\bar{\phi}^2$ matrix for each combination of ($H_{ld,100\,au}$, $i$) and ($H_{ld,100\,au}$, $\beta$). These maps ($\bar{\phi}^2$) show that larger inclinations and flaring are favored. Independent of the inclination and flaring exponent, we also find that a scale height of $H_{ld,100\,au} \sim 3$ au is favored.

## ORCID iDs

M. Villenave https://orcid.org/0000-0002-8962-448X

L. Podio https://orcid.org/0000-0003-2733-5372

G. Duchêne https://orcid.org/0000-0002-5092-6464

K. R. Stapelfeldt https://orcid.org/0000-0002-2805-7338

C. Melis https://orcid.org/0000-0001-9834-7579

C. Carrasco-Gonzalez https://orcid.org/0000-0003-2862-5363

V. J. M. Le Gouellec https://orcid.org/0000-0002-5714-799X

F. Ménard https://orcid.org/0000-0002-1637-7393

M. de Simone https://orcid.org/0000-0001-5659-0140

C. Chandler https://orcid.org/0000-0002-7570-5596

A. Garufi https://orcid.org/0000-0002-4266-0643

C. Pinte https://orcid.org/0000-0001-5907-5179

E. Bianchi https://orcid.org/0000-0001-9249-7082

C. Codella https://orcid.org/0000-0003-1514-3074

## References

Andre, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Tucson, AZ: Univ. Arizona Press), 59

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41
Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

Armitage, P. J. 2015, arXiv:1509.06382

Barrière-Fouchet, L., Gonzalez, J. F., Murray, J. R., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185

Beck, T. L., Bary, J. S., & McGregor, P. J. 2010, ApJ, 722, 1360

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

CASA team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 1041

Codella, C., Cabrit, S., Gueth, F., et al. 2014, A&A, 568, L5

Doi, K., & Kataoka, A. 2021, ApJ, 912, 164

Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2022, arXiv:2203.09759

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

Dunham, M. M., Allen, L. E., Evans, N. J. I., et al. 2015, ApJS, 220, 11

Flores, C., Duchêne, G., Wolff, S., et al. 2021, AJ, 161, 239

Galli, P. A. B., Loinard, L., Bouy, H., et al. 2019, A&A, 630, A137

Gräfe, C., Wolf, S., Guilloteau, S., et al. 2013, A&A, 553, A69

Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323

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

Harris, R. J., Cox, E. G., Looney, L. W., et al. 2018, ApJ, 861, 91

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

Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353

Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117

Kristensen, L. E., & Dunham, M. M. 2018, A&A, 618, A158

Lee, C.-F. 2020, A&ARv, 28, 1

Lee, C.-F., Codella, C., Ceccarelli, C., & López-Sepulcre, A. 2022, ApJ, 937, 10

Lee, C.-F., Hirano, N., Zhang, Q., et al. 2014, ApJ, 786, 114

Lee, C.-F., Li, Z.-Y., Ho, P. T. P., et al. 2017, SciA, 3, e1602935

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

Liu, Y., Bertrang, G. H. M., Flock, M., et al. 2022, SCPMA, 65, 129511

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

Lucas, P. W., & Roche, P. F. 1997, MNRAS, 286, 895

Maret, S., Maury, A. J., Belloche, A., et al. 2020, A&A, 635, A15

Melis, C., Duchêne, G., Chomiuk, L., et al. 2011, ApJL, 739, L7

Mendigutía, I., Brittain, S., Eiroa, C., et al. 2013, ApJ, 776, 44

Michel, A., Sadavoy, S. I., Sheehan, P. D., Looney, L. W., & Cox, E. G. 2022, ApJ, 937, 104

Murillo, N. M., Harsono, D., McClure, M., Lai, S. P., & Hogerheijde, M. R. 2018, A&A, 615, L14

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

Nakatani, R., Liu, H. B., Ohashi, S., et al. 2020, ApJL, 895, L2

Ohashi, S., Nakatani, R., Liu, H. B., et al. 2022, ApJ, 934, 163

Padgett, D. L., Brandner, W., Stapelfeldt, K. R., et al. 1999, AJ, 117, 1490

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Podio, L., Garufi, A., Codella, C., et al. 2020, A&A, 642, L7

Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328

Rodmann, J., Henning, T., Chandler, C. J., Mundy, L. G., & Wilner, D. J. 2006, A&A, 446, 211

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

Rota, A. A., Manara, C. F., Miotello, A., et al. 2022, A&A, 662, A121

Sakai, N., Oya, Y., Higuchi, A. E., et al. 2017, MNRAS, 467, L76

Sheehan, P. D., Tobin, J. J., Li, Z.-Y., et al. 2022, ApJ, 934, 95

Silsbee, K., Akimkin, V., Ivlev, A. V., et al. 2022, ApJ, 940, 188

Teague, R., Bae, J., & Bergin, E. A. 2019, Natur, 574, 378

Tobin, J. J., Hartmann, L., Calvet, N., & D'Alessio, P. 2008, ApJ, 679, 1364

van't Hoff, M. L. R., Harsono, D., Tobin, J. J., et al. 2020, ApJ, 901, 166

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

Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11

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

Wolf, S., Padgett, D. L., & Stapelfeldt, K. R. 2003, ApJ, 588, 373

Wolf, S., Schegerer, A., Beuther, H., Padgett, D. L., & Stapelfeldt, K. R. 2008, ApJL, 674, L101

Wolf, S. G., Duchêne, G., Stapelfeldt, K. R., et al. 2021, AJ, 161, 238

Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588

Zhang, C.-P., Launhardt, R., Liu, Y., Tobin, J. J., & Henning, T. 2021, A&A, 646, A18