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

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, ¹² 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 gencal 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 (∼02) 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.

Zoom In Zoom Out Reset image size

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 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.

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 08 × 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 λ ≤ 13μ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 (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 ≲01. 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 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.

The spectral index of millimeter–centimeter emission, defined as α in F_ν ∝ ν^α , 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 a spectral index of α_{14−66 mm} = 0.97 ± 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 α_{0.9–9.2 mm} = 2.54 ± 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; 030 × 024; 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.

Zoom In Zoom Out Reset image size

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 05 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.1.1. Model Description

4.1.2. Fixed and Varied Parameters

We use a stellar effective temperature of T_eff = 4500 K and a stellar radius of 3.7R_⊙ (Gräfe et al. 2013). We fix the inner disk radius to R_in = 0.1 au and the outer radius to R_out = 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}_{\min }=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_{ld,100 au}, the inclination i, the maximum grain size a_max, the surface density exponent p, the flaring exponent β, and the dust mass M_dust. 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.

Table 3. Overview of the Parameter Ranges of the Grid

Hld,100 au	(au)	[1, 3, 4, 5, 6]
i	(°)	[88, 89, 90]
amax	(μm)	[100, 1000, 10000]
p	(−)	[0.5, 1, 1.5]
β	(−)	[1.06, 1.14, 1.22, 1.3]
${\mathrm{log}}_{10}({M}_{\mathrm{dust}}/{M}_{\odot })$	(−)	[−5.0, −4.75, −4.5, −4.25, −4.0
		−3.75, −3.5, −3.25, −3.0, −2.75
		−2.5, −2.25, −2.0]

Note. H_{ld,100 au}, i, a_max, β, and M_dust 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.

Download table as:  ASCIITypeset image

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}_{\min }=10$ μm, with varied a_max; 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_{ld,100 au}, i, a_max, p, and β), we first determine the dust mass M_dust 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:

$$\begin{eqnarray}&&{\chi }^{2}=\sum \displaystyle \frac{1}{n}{\left(\displaystyle \frac{{F}_{d}-{F}_{m}}{\sigma }\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where F_d and F_m 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 ${\bar{\chi }}_{\mathrm{maj}}^{2}$, taken as the mean ${\chi }_{\mathrm{maj}}^{2}$ between 0.89 mm, 2.1 mm, and 9.2 mm. For each set of parameters (H_{ld,100 au}, i, a_max, p, and β), the first step determines the best-mass model (best M_dust), which corresponds to that with the lowest ${\bar{\chi }}_{\mathrm{maj}}^{2}$.

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 ${\bar{\chi }}_{\mathrm{maj}}^{2}$, we evaluate ${\bar{\chi }}_{\min }^{2}$ 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 ${\chi }_{\min }^{2}$ using Equation (1) along the minor axis profiles, and obtain ${\bar{\chi }}_{\min }^{2}$ 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}=\tfrac{{\bar{\chi }}_{\mathrm{maj}}^{2}+{\bar{\chi }}_{\min }^{2}}{2}$.

4.2.1. Grid Best Models

We show the distributions of the dust masses for each combination of (H_{ld,100 au}, i, a_max, p, and β), as obtained during the first step, in Figure 5. The median of the distributions is M_dust = 3.2 × 10⁻⁴ M_⊙, and the 25% and 75% quartiles reach values of 1.8 × 10⁻⁴ M_⊙ and 1.8 × 10⁻³ M_⊙, respectively. For comparison, we also indicate the dust masses of the top 10% of best models (based on ϕ²) in this figure. They fall within the range of the best dust masses of the distribution from the first step.

Zoom In Zoom Out Reset image size

In Figure 5, we can also see that 24 models saturate to the high mass limit of our grid, namely at M_dust = 10⁻² M_⊙. These are exclusively models with H_{ld,100 au} = 1 au and i = 88°, 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_d = 10⁻² M_⊙, 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_⋆ ∼ 1.6 M_⊙; Z.‐Y. D. Lin et al., in preparation). Thus, this first modeling step indicates that IRAS04302 likely does not have the combination of H_{ld,100 au} = 1 au and i = 88°.

After performing both fitting steps, we find that the parameters of our best model are H_{ld,100 au} = 3 au, i = 89°, ${a}_{\max }=1$ mm, β = 1.14, p = 1.5, and M_dust = 3 × 10⁻⁴ M_⊙. This model is associated with the following metrics: ϕ² = 8.1, ${\chi }_{\mathrm{maj}}^{2}=5.5$, and ${\chi }_{\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 ϕ², ${\chi }_{\min }^{2}$, and ${\chi }_{\mathrm{maj}}^{2}$ values for each H_{ld,100 au}, i, a_max, p, and β value of the grid, showing them in Figure 6. We find that H_{ld,100 au} is relatively well constrained, while there are trends for the values of a_max and i. Specifically, the models indicate that 1au < H_{ld,100 au} < 6au and suggest that ${a}_{\max }\lt 1$ cm and i > 88°. For the scale height, we note that ${\bar{\chi }}_{\min }^{2}$ and ${\bar{\chi }}_{\mathrm{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 β and p appear relatively flat, and these parameters are thus not constrained.

Zoom In Zoom Out Reset image size

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_{ld,100 au} = 3au, i = 90°, p = −1, and β = 1.14 (the best values from Figure 6). The corresponding dust masses of each model (as determined in our first modeling step) are M_dust = 10⁻⁴ M_⊙ for the models with ${a}_{\max }=100$ μm and ${a}_{\max }=1$ mm, and M_dust = 3 × 10⁻⁴ M_⊙ for the model with ${a}_{\max }=$ 1 cm.

Zoom In Zoom Out Reset image size

In Figure 7, we find that the models with ${a}_{\max }=100\mu {\rm{m}}$ and ${a}_{\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}_{\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_{ld,100 au}, i, p, and β) also exclude a maximum grain size of ${a}_{\max }=1$ cm, as can be seen in Figure 6. In Appendix C, we show a similar result, when considering the ${\bar{\chi }}_{\mathrm{maj}}^{2}$ values for different pairs of (a_max, p) and (a_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}_{\max }=100\mu {\rm{m}}$ and ${a}_{\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}_{\max }=100\mu {\rm{m}}$, p = −1, β = 1.14, i = 90°, and varying H_{ld,100 au}. The models with H_{ld,100 au} = 1, 3, and 6 au have dust masses of M_dust = 3 × 10⁻⁵, 10⁻⁴, and 3 × 10⁻⁴ M_⊙, respectively. Due to the high inclination of the models, these cuts are dominated by the vertical extent of the disk.

Zoom In Zoom Out Reset image size

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 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 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_{ld,100 au} or β. In the two figures, we fix i = 90°, ${a}_{\max }=100$ μm, and p = −1. The top figure shows the impact of the dust scale height on the minor axis size variation, for β = 1.14, and the bottom figure shows the effect of the flaring on the shape of the minor axis size variation, for H_{ld,100 au} = 6 au ¹³ .

Zoom In Zoom Out Reset image size

We find that while the H_{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 β 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.

In Section 4.2.3, we find that the large dust scale height is constrained within 1au < H_{ld,100 au} < 6au. 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, ¹⁴ which translates into a gas scale height of H_g,100au ∼ 7 au for a stellar mass of M_⋆ = 1.6 M_⊙ (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 α/St ratio, 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:

$$\begin{eqnarray}&&{H}_{d}={H}_{g}{\left(1+\displaystyle \frac{\mathrm{StSc}}{\alpha }\right)}^{-1/2},\end{eqnarray} \tag{ 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_{ld,100 au} > 1 au and H_g ∼ 7 au, we find that [α/St]_100au > 2 × 10⁻². This value is higher than previous estimates for Class II systems (Table 4), which we discuss in Section 5.2.

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 ¹²CO lines are used, while rarer species such as C¹⁷O, HCO⁺, 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 ${\dot{M}}_{\mathrm{acc}}=\tfrac{{L}_{\mathrm{bol}}{R}_{\star }}{{{GM}}_{\star }(1-{R}_{\star }/{R}_{\mathrm{in}})}$ (Gullbring et al. 1998), assuming that the bolometric luminosity is dominated by the accretion luminosity. We assume a stellar radius of R_⋆ = 2R_⊙ (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):

$$\begin{eqnarray}&&{t}_{\mathrm{settle}}=\displaystyle \frac{\rho }{{\rho }_{m}}\displaystyle \frac{{v}_{\mathrm{th}}}{a}\displaystyle \frac{1}{{{\rm{\Omega }}}_{K}^{2}}\propto \rho \displaystyle \frac{\sqrt{T}}{a}\displaystyle \frac{{R}^{3}}{{M}_{\star }},\end{eqnarray} \tag{ 3 }$$

where ρ_m is the dust material density, ρ is the gas density, a is the particle size, ${v}_{\mathrm{th}}=\sqrt{8{k}_{B}T/\pi \mu {m}_{{\rm{H}}}}$ is the thermal speed of the molecules, and ${{\rm{\Omega }}}_{K}=\sqrt{{{GM}}_{\star }/{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)—ρ = 10⁻¹² g cm⁻³, ρ_m = 3 g cm⁻³, and T = 50 K—we find that the settling timescale for grains of 1 mm is t_settle,100au ∼ 0.4–0.09 Myr at 100 au, and t_settle,50au ∼ 0.05–0.01 Myr at 50 au, for M_⋆ = 0.45–2M_⊙. 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 α/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.

We present the residual maps of the best model in Figure 13. The parameters for the best model are: H_{ld,100 au} = 3 au, i = 89°, ${a}_{\max }=1$ mm, β = 1.14, p = 1.5, and M_dust = 3 × 10⁻⁴ M_⊙. We also note that the maps shown here do not have the small grain disk and envelope layers discussed in Appendix B.

Zoom In Zoom Out Reset image size

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_dust, i, and p. In our first fitting step, we determine the best mass that allows to reproduce the major axis profile, thus the possible degeneracies between M_dust and a_max are already taken into account. In Figure 14, we show the detailed ${\bar{\chi }}_{\mathrm{maj}}^{2}$ median matrix for each combination of (${a}_{\max },i$) and (${a}_{\max },p$). Both matrices show that, independently of the values of i and p, ${a}_{\max }=1$ cm (the lowest line) leads to significantly higher median ${\bar{\chi }}_{\mathrm{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.

Zoom In Zoom Out Reset image size

For all a_max, it is also clear that higher inclinations and higher surface density exponents are preferred (i > 88° and 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.

Zoom In Zoom Out Reset image size

The models presented in Figure 8 were computed for an inclination of i = 90° and flaring of β = 1.14. However, both the inclination and the flaring exponent i and β can also impact the apparent minor axis size of the disk. In Figure 16, we present the median ϕ² matrix for each combination of (H_{ld,100 au}, i) and (H_{ld,100 au}, β). These maps (ϕ²) 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} ∼ 3 au is favored.

Zoom In Zoom Out Reset image size