Issue |
A&A
Volume 688, August 2024
|
|
---|---|---|
Article Number | A204 | |
Number of page(s) | 10 | |
Section | Interstellar and circumstellar matter | |
DOI | https://doi.org/10.1051/0004-6361/202449253 | |
Published online | 21 August 2024 |
Uncertainties of the dust grain size in protoplanetary disks retrieved from millimeter continuum observations
1
Purple Mountain Observatory & Key Laboratory of Radio Astronomy, Chinese Academy of Sciences,
10 Yuanhua Road, Qixia District,
Nanjing
210023,
PR
China
e-mail: yliu@pmo.ac.cn
2
School of Astronomy and Space Science, University of Science and Technology of China,
96 Jinzhai Road,
Hefei
230026,
PR
China
Received:
17
January
2024
Accepted:
13
June
2024
Context. Investigating the dust grain size and its dependence on substructures in protoplanetary disks is a crucial step in understanding the initial process of planet formation. Spectral indices derived from millimeter observations are used as a common probe for grain size. Converting observed spectral indices into grain sizes is a complex task that involves solving the radiative transfer equation, taking into account the disk structure and dust properties.
Aims. Under the assumption of vertically isothermal disks, the solution to the radiative transfer equation can be approximated with an analytic expression, with which the fitting procedure can be done very fast. Our work aims to investigate the applicability of this method to grain size retrieval.
Methods. We ran reference radiative transfer models with known disk properties, and generated four synthetic images at wavelengths of 0.8, 1.3, 3, and 7.8 mm, representing high-resolution continuum observations. Rings and gaps were considered in the setup. We fit the synthetic images using the analytic solution to investigate the circumstances under which the input grain sizes can be recovered.
Results. Fitting images at only two wavelengths is not sufficient to retrieve the grain size. Fitting three images improves the retrieval of grain size, but the dust surface density is still not well recovered. When taking all of the four images into account, degeneracies between different parameters are highly reduced, and consequently the best-fit grain sizes are consistent with the reference setup at almost all radii. We find that the inclination angle has a significant impact on the fitting results. For disks with low inclinations, the analytic approach works quite well. However, when the disk is tilted above ~60°, neither the grain size nor the dust surface density can be constrained, as the inclination effect will smooth out all substructures in the radial intensity profile of the disk.
Key words: radiative transfer / protoplanetary disks / circumstellar matter
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Dust grains in protoplanetary disks are the building blocks of planets. During the planet formation process, the size of dust particles increases by more than ten orders of magnitude (e.g., Armitage 2010). It is generally believed that micron-sized dust grains collide and stick together by van der Waals forces, forming large aggregates up to the scale of millimeters, and that the gravitational force between kilometer-sized planetesimals is enough to accrete surrounding dust, making the bodies grow to larger sizes (e.g., Chokshi et al. 1993; Dominik & Tielens 1997; Birnstiel et al. 2016; Drążkowska et al. 2023; Birnstiel 2023).
In existing theories, it is difficult for dust grains to grow from the scale of millimeters to that of meters; this is known as the meter-size barrier (Weidenschilling 1977). On the millimeter scale, the relative velocities between dust particles become significant, causing them to collide with higher energies. The increased energy in collisions often leads to bouncing instead of sticking, preventing the grains from further growth (e.g., Testi et al. 2014). Furthermore, radial drift is another obstacle to grain growth. Dust grains lose angular momentum due to the sub-Keplerian motion of gas in protoplanetary disks, and so fast migrate toward the central star (Whipple 1972; Adachi et al. 1976; Nakagawa et al. 1986; Youdin & Shu 2002; Brauer et al. 2007). Consequently, dust grains, especially millimeter pebbles, will drift to the inner region before they can grow into larger sizes (Takeuchi & Lin 2002; Brauer et al. 2008; Birnstiel et al. 2009). However, mechanisms that facilitate dust grains to accumulate and grow likely work at specific locations in the disk.
High-resolution Atacama Large Millimeter/submillimeter Array (ALMA) observations reveal that substructures, for instance rings and gaps, are common in protoplanetary disks (e.g., Long et al. 2018; Andrews et al. 2018; Andrews 2020; Bae et al. 2023). It has been proposed that the origin of these finely scaled features is related to various processes, such as (magneto-)hydrodynamic instabilities (e.g., Lesur et al. 2022), snow lines (e.g., Zhang et al. 2015), and planets (e.g., Kley & Nelson 2012; Paardekooper et al. 2022). From a theoretical perspective, disk substructures act as pressure bumps that can trap dust particles, which increases the density and facilitates grain growth (Whipple 1972; Pinilla et al. 2012b; Teague et al. 2018a,b; Rosotti et al. 2020). Investigating the relationship between grain size and substructures is crucial to understanding dust evolution and developing models of planet formation.
A common probe for grain size is the millimeter spectral index , where
and
are the measured flux densities at frequencies ν1 and ν2, respectively. In the (sub)millimeter regime, the mass absorption coefficient (κν) of dust grains can be approximated as
. In the optically thin case, the spectral index is linked to the opacity slope via βmm = αmm − log(Bν1/Bν2)/ log(ν1 / ν2), where Bν is the Planck function. With a given dust composition and shape, the opacity slope, βmm, is dependent on the maximum grain size, amax (see Fig. 1). For micron-sized solids, βmm does not change with amax, and is close to the average value of the interstellar medium dust; that is, ~1.6 ± 0.1 (Li & Draine 2001; Planck Collaboration Int. XLVIII 2016). The vertical dashed line in Fig. 1 denotes the boundary of amax ~0.4 mm, which divides the amax values into two regions according to the β(1.3–3 mm) curve. In the region on the left side, βmm monotonically increases with amax, whereas βmm monotonically decreases with amax on the right side of the dashed line (Draine 2006; Ricci et al. 2010b; Birnstiel et al. 2018). The amax boundary between the two different trends of βmm depends on the wavelengths used to calculate the opacity slope. It can be seen that a β value lower than ~1.7 is a strong indication of dust grains of millimeter sizes or even larger (e.g., Beckwith & Sargent 1991; Andrews & Williams 2005, 2007; Ricci et al. 2010a; Pérez et al. 2012; Tazzari et al. 2016; Ansdell et al. 2018).
Constraining dust grain size from observed spectral indices is a challenging task. This is because for a βmm value there may be two solutions of amax. The situation is even worse when the grain size in the disk is smaller than ~0.1 mm, since then βmm is no longer sensitive to amax (see Fig. 1). Moreover, translating spectral indices into grain sizes requires knowledge about the dust temperature distribution in the disk, which is poorly known. In order to self-consistently calculate the dust temperature and constrain the grain size, Monte Carlo radiative transfer models are built to fit the spectral energy distribution and millimeter observations (e.g., Pinte et al. 2016; Liu et al. 2017; Li et al. 2023). These sophisticated models, usually in two dimensions, can properly treat the effect of optical depth. However, this kind of methodology requires a lot of computational resources, and is therefore preferable for case studies. Sierra et al. (2019) derived an approximate solution to the radiative transfer equation for vertically isothermal disks. The emergent intensity from the disk is expressed as an analytic formula that involves the dust temperature, surface density, and grain properties, including amax. Parameter estimation based on the analytic solution is efficient and flexible, and has therefore been widely adopted in recent works (e.g., Carrasco-González et al. 2019; Sierra et al. 2021; Macías et al. 2021; Zhang et al. 2023; Ohashi et al. 2023).
Nevertheless, the analytic solution needs an input for the dust temperature to obtain βmm from αmm, and the vertical disk structure is ignored. How, and to which degree these simplifications affect the results needs to be investigated. In this work, we aim to examine the applicability of using the analytic models to constrain the grain size, and explore under which circumstances the grain size can be properly retrieved. We first build self-consistent reference disk models with known properties, and generate synthetic millimeter continuum images. Then, we fit the synthetic images to search for the best-fit parameter set. Finally, we compare the grain sizes of the best-fit model with the reference setup to assess the quality of retrieval. The paper is organized as follows. The establishment of the reference disk models is introduced in Sec. 2. Section 3 describes the fitting process. We present and discuss the results in Sec. 4. The paper ends up with a summary in Sec. 5.
![]() |
Fig. 1 Dust opacity slope, βmm, as a function of the maximum grain size, amax. The slopes between two wavelengths given in the legend are derived by calculating the mass absorption coefficients (κabs). To calculate κabs, we adopted the MIE scattering theory, assuming the DSHARP dust composition (see Sec. 2.3). The grain size distribution follows a power law, N(a) ∝ a−3.5, with the minimum dust size fixed to amin = 0.01 μm. The vertical dashed line marks the 0.4 mm boundary that divides the amax values into two regimes (a small size region and a large size region) according to the β(1.3–3 mm) curve. |
2 Reference models
We created two disk models that differ only in the radial profile of dust grain size. The model parameters are listed in Table 1. We used the RADMC-3D1 code (Dullemond et al. 2012) to calculate the dust temperature and simulate the millimeter continuum images. This section describes the model setup, including stellar parameters, density distribution, dust properties, and observational parameters.
2.1 Stellar parameters
We assume that the disk is passively heated by stellar irradiation. The stellar properties were set according to DS Tau, which is a T Tauri star. It has a spectral type of M0.4 and an effective temperature of 3800 K (Herczeg & Hillenbrand 2014). The stellar luminosity was fixed to 0.7 L⊙ (Li et al. 2023). The incident stellar spectrum was taken from the Kurucz atmosphere grid, assuming solar metallicity and a log g = 3.5 (Kurucz 1994). We solved the radiative transfer problem at 160 wavelengths that were logarithmically spaced between 0.1 μm and 10 000 μm.
2.2 Dust density distribution
The disk extends from an inner radius of Rin = 0.1 AU to an outer radius of Rout = 100 AU, and features a flared geometry. The volume density of the dust grains is given as
(1)
The slope, βH, describes the extent of disk flaring, and H100 represents the scale height at the radial distance of R = 100 AU from the central star. We prescribed the dust surface density as
(3)
where the power-law slope was set to −0.5, which has been found to be a reasonable value in previous works (Davis 2005; Andrews et al. 2009; Williams & Cieza 2011; Tazzari et al. 2017; Carrasco-González et al. 2019). The proportionality constant, ∑0, was derived by normalizing the total dust mass, Mdust, which we fixed to a typical value (1 × 10−4 M⊙) for protoplanetary disks (e.g., Williams & Cieza 2011; Miotello et al. 2023). To simulate a pair consisting of a gap and a ring, as is commonly observed in protoplanetary disks, a Gaussian depletion term and a Gaussian magnification term were multiplied. The position and width of the gap are denoted as RGap and σGap, respectively. Similarly, RRing and σRing refer to the position and width of the ring. The parameter, ηGap, represents the depletion factor of in the gap, while ηRing represents the magnification factor of in the ring. The black line in the upper panel of Fig. 2 shows the dust surface density of the reference model. We note that the dust surface densities of model 1 and model 2 are the same.
Parameter values of the reference models.
![]() |
Fig. 2 Properties of the reference models. Top panel: dust surface density. The dust surface densities of model 1 and model 2 are the same. Middle panel: radial profile of the maximum grain size of model 1 (red line) and model 2 (blue line). The horizontal dashed line marks the 0.4 mm boundary defined in Fig. 1. Bottom panel: optical depth at the wavelength of 1.3 mm. The color coding is identical to that in the middle panel. The τ1.3mm = 1 condition is indicated with the horizontal dashed line. |
2.3 Dust properties
We assume that dust grains are compact homogeneous spheres, composed of 20% water ice, 33% astronomical silicates, 7% troilite, and 40% refractory organics, with an average density of ρdust = 1.675 g cm−3. The percentages are given for the mass fraction of each component. This dust prescription was introduced to interpret the ALMA data from the Disk Substructures at High Angular Resolution Project (DSHARP, Andrews et al. 2018). The complex refractive indices were provided by Birnstiel et al. (2018), and the dielectric function and dust absorption or scattering coefficients were derived using the Bruggeman effective medium theory (Bruggeman 1935) and the Mie theory, respectively.
The grain size distribution is assumed to be a power law N(a) ∝ a−3.5, where a is the grain radius, and N(a) denotes the number of dust grains within the size interval [a, a + da]. We fixed the minimum dust size to amin = 0.01 μm and allowed the maximum grain size, amax, to vary as a function of R. The radial profile of amax is similar to that of ∑dust, expressed as:
(4)
where a0 is the amax at R = 1 AU. We set the power law index to −1, which has been found to be typical in protoplanetary disks (Trotta et al. 2013; Tazzari et al. 2016; Guidi et al. 2016; Carrasco-González et al. 2019; Sierra et al. 2021). The red and blue lines in the middle panel of Fig. 2 show the radial profiles of amax for models 1 and 2, respectively. The maximum grain size of model 1 at each radius is ten times larger than that of model 2. More importantly, all the amax(R) in model 1 are larger than the 0.4 mm boundary defined in Fig. 1, which implies that the model is in a parameter space with a monotonic dependence between amax and βmm if 1.3 mm and 3.0 mm images are used in the fitting procedure. On the contrary, amax(R) in model 2 spans the 0.4 mm boundary, and hence the model is located in a parameter space with a non-monotonic dependence between amax and βmm. This setup allows us to examine the performance of the retrieval framework with the analytic solution when a βmm value corresponds to two amax values (see Fig. 1). Once amax and ∑dust are known, we can calculate the optical depth, τν. As is shown in the bottom panel of Fig. 2, model 1 is optically thin at 1.3 mm, while model 2 is optically thin as well, except for the very inner disk and the ring center.
Numeric simulations of dust trapping have demonstrated that larger dust grains are more concentrated toward the ring center than the smaller ones (e.g., Pinilla et al. 2012a, 2015a,b). Therefore, we prescribed the radial profile of the grain size to be the same as the dust surface density, where the grain size increases in the ring and decreases in the gap. In order to have more experiments, we also ran the fitting procedure for three smooth disks without substructures. The dust surface density was fixed to ∑dust = ∑0R−0.5, and the radial profile of amax follows a power law, amax(R) = a0R−γ, with the power law index being γ = − 0.5, −1, and −1.5, respectively. The results of these additional tests are presented and discussed in Sec. 4.5.
To set up amax(R) in the radiative transfer simulation, we divided the disk into 40 radial bins that are linearly spaced with a width of 2.5 AU. The value of amax for each radial bin was determined based on the amax(R) profile; for example, Eq. (4). The radial bin size was chosen to be half of the beam size (i.e., 5 AU, see Sec. 2.4) of the synthetic images, ensuring that our model is sufficient to resolve the disk substructure. Using more radial bins needs more computational resources; however, the results are not significantly affected.
2.4 Observational parameters
The synthetic observations were generated assuming a distance of 100 pc. We simulated continuum images at 0.8, 1.3, 3.0, and 7.8 mm, and convolved the raw images from RADMC-3D with a circular beam of 0.05″, mimicking high-resolution ALMA and Very Large Array (VLA) observations. The beam size corresponds to a physical scale of 5 AU at the distance of 100 pc. The synthetic images of model 1 are shown in Fig. 3, in which the gap and ring are clearly seen. The disk inclination was set to i = 0°; a face-on orientation. We explore the effects of disk inclination on the fitting results in Sec. 4.4.
We extracted the radial intensity profile from the synthetic images by applying an azimuthal average on a series of concentric ellipses with shapes determined by the disk inclination. The radial intensity profiles are the observables to be fitted.
3 Fitting approach
To fit the synthetic images, we adopted the analytic solution to the radiative transfer equation, which was introduced by Sierra et al. (2019). It is now a consensus that dust scattering should be considered in interpreting millimeter observations of protoplanetary disks (e.g., Zhu et al. 2019; Liu 2019). When taking dust scattering into account, the dust continuum emission from the disk can be written as
(5)
where T is the dust temperature, and μ = cos(i) represents the cosine of the disk inclination (Miyake & Nakagawa 1993; Sierra et al. 2019). The dust surface density (∑dust), absorption coefficient (κν) and albedo (ων) work together to determine the optical depth, . In addition, the scattering correction term (F) is denoted as
(6)
The quantity ϵν equals . In Eq. (5), we have four free parameters: i, T, τν, and ων. The disk inclination (i) can be well constrained by high-resolution observations (e.g., Long et al. 2018; Huang et al. 2018). Hence, we do not explore this parameter throughout our work. In previous studies, the dust temperature (T) is either considered as a free parameter (Sierra et al. 2019; Macías et al. 2021; Zhang et al. 2023; Ohashi et al. 2023) or assumed to follow the expected radial profile for a passively irradiated disk (Kenyon & Hartmann 1987; Dullemond et al. 2001; Sierra et al. 2021; Doi & Kataoka 2023). In order to reduce the model degeneracy, we fixed T at each radius to the mass-averaged mean dust temperature that was obtained from the RADMC-3D simulation. The optical depth (τμ) was determined by ∑dust, κν, and ων. The absorption coefficient (κν) and albedo (ων) depend only on the maximum grain size (amax) once the dust composition and minimum grain size (amin) are given. Therefore, we are left with a parameter space that consists of ∑dust and amax.
For each parameter set, {∑dust, amax}, we calculated the likelihood via
(7)
is the chi-square statistic. The intensities extracted from the synthetic images at different wavelengths are denoted as , while
refers to the model intensity prescribed by Eq. (5). The parameter
represents the uncertainties of the observed intensity profiles. We assume that the uncertainties are dominated by the absolute calibration error, and took 10% for all of the four considered wavelengths. This assumption is reasonable according to ALMA and VLA observations of protoplanetary disks (e.g., Long et al. 2018; Andrews et al. 2018; Francis et al. 2020; Tobin et al. 2020).
The fitting procedure was performed at each integer value of the radius of the disk in units of AU. At each radius, we calculated a large grid of model intensities by sampling ∑dust and amax. For ∑dust, the range for exploration is from 10−4 to 10 g cm−2. In previous works, amax was commonly sampled at two different ranges: one is located in the small size regime, [10 μm, 0.4 mm], and the other is in the large size domain, [0.4 mm, 10 cm] (see Fig. 1). Then the fitting was conducted twice, separately within the small and large size boundaries. The final results were obtained by combining the two fitting solutions (Sierra et al. 2021; Macías et al. 2021; Zhang et al. 2023). The purpose of this two-round practice is to reduce the degeneracy between the small and large size regimes. However, combining the results of fitting in the small and large size regimes is a challenge, and is usually guided by empirical knowledge from theories of dust evolution. For instance, it is generally believed that dust grains grow faster in regions with higher densities; for example, in the inner disk. Moreover, millimeter- and centimeter-sized dust grains would rapidly drift toward the center of the disk (Weidenschilling 1977). As a result, the grain population in the inner disk is dominated by large dust, while small dust grains are more concentrated within the outer disk (e.g., Birnstiel et al. 2016; Birnstiel 2023). In addition, disk substructures, frequently observed with ALMA, have been interpreted as the outcome of large dust grains being trapped toward local gas pressure maxima (e.g., Pinilla et al. 2012a, 2015a,b). Therefore, dust grains in the rings are expected to be larger than those in the low-surface-density gaps. Nevertheless, we do not know the true amax(R) in real protoplanetary disks. Hence, in our work, we have treated the dust size as a black box, and sampled the parameter from 10 μm all the way to 10 cm. Both ∑dust and amax were sampled in a logarithmic manner with a stepwidth of 0.01 dex.
In the fitting procedure, the joint probability distribution, P(amax, ∑dust), was marginalized to obtain the probabilities of ∑dust and amax. As an example, the shaded area in panel a of Fig. 4 illustrates the regions in the parameter space that can reproduce the observations at R = 5 AU in model 1. The marginalized probabilities of amax and ∑dust are shown in panels b and c, respectively. The best-fit parameter set, indicated with the red star, is identified as the intersection where both ∑dust and amax reach their peak probability. Moreover, we used a Monte Carlo method to derive the final best-fit parameter set and its error in order to minimize the influence of the uncertainty of observations () on the results. We generated a set of 5000 intensity profiles at each wavelength by randomly adding Gaussian noise, with the standard deviation being 10% of the absolute value of the intensities. The final best model was determined as the 50th percentile among the results from fitting the 5000 random intensities, and the confidence interval was associated with the 25th and 75th percentiles. The effect of the Monte Carlo method on the parameter uncertainty is described in Sec. 4.3. Once the fitting procedure for each radius was done, we obtained the best-fit model profiles of ∑dust and amax, and could compare with the reference setup.
![]() |
Fig. 3 Synthetic images of model 1 at wavelengths of 0.8 mm, 1.3 mm, 3.0 mm, and 7.8 mm. The images are observed with a face-on orientation. The beam with a size of 0.05″ is indicated with the white dot in the lower left corner of each panel. |
![]() |
Fig. 4 Fitting process for model 1 at R = 5 AU. The normalized joint probability distribution is shown as the color map in panel a. The marginalized probabilities for amax and ∑dust are indicated with the black lines in panels b and c, respectively. The red star denotes the best-fit parameter set identified from the fitting procedure. |
4 Results and discussion
In this section, we present the results by fitting images at only two wavelengths, three wavelengths, and all four wavelengths. Moreover, we investigate the influence of the disk inclination on the fitting results.
4.1 Fitting to the 1.3 and 3.0 mm images
The spectral index, measured between 1.3 mm and 3.0 mm, is suitable for grain size studies (e.g., Ricci et al. 2010b,a; Andrews 2020; Tazzari et al. 2021). At these two wavelengths, disks have relatively strong thermal emission and low optical depth, which are favorable conditions for disk characterization given limited observation time. Optical depth is higher toward shorter wavelengths, whereas long wavelength observations take a much longer integration, and contaminations by other radiation processes, such as free-free emission, should be carefully treated (e.g., Macías et al. 2021; Rota et al. 2024). Therefore, we first fit the 1.3 mm and 3 mm images. The results for model 1 are presented in panels a–d of Fig. 5. The dashed red lines in panels a and b are the reference profiles for amax and ∑dust, respectively, while the cyan lines refer to the best-fit results. It should be noted that the reference model profiles were convolved with a same beam size (0.05″) as the observations. The gray color scale represents the probability distribution of parameter values obtained from fitting the 5000 Monte Carlo models at each radius. The vertical bar indicates the confidence interval of the best-fit result. For a better visualization, we only show the typical uncertainty that is the median of the uncertainties at all radii. A comparison of the intensities and spectral index between the reference setup and the best-fit results is given in panels c and d.
As can be seen from Fig. 5, the best-fit amax(R) based on the 1.3 mm and 3 mm images match with the reference values only for the inner 15 AU and the ring region. In other parts of the disk, both amax and ∑dust were not satisfactorily retrieved. The errors are pretty large in general. For model 2, the results are similar, and are shown in panels e–h of Fig. 5. Nevertheless, except for the inner 10 AU, the best-fit surface densities are close to the reference setup.
The mismatch between the model and synthetic observation is mainly due to the following reasons. Images of two wavelengths can derive only one spectral index; however, there are two solutions of amax for one particular value of spectral index even in the optically thin regime (see Fig. 1). One is located in the small size region, and the other is in the large size domain. This fact leads to an ambiguity in identifying the optimum amax, which will in turn affect the retrieval for ∑dust. Therefore, if amax cannot be well retrieved, ∑dust often deviates from the reference values at a significant level. Even at locations where amax is well retrieved, ∑dust may not be correctly obtained, directly illustrating the large uncertainty in the results based on images at only two wavelengths.
We also examined the results by fitting other pairs of images, such as a combination of 1.3 and 7.8 mm or 3.0 and 7.8 mm. These attempts arrive at the similar conclusion.
![]() |
Fig. 5 Results of fitting to the 1.3 mm and 3.0 mm images for model 1 (left column) and model 2 (right column). For model 1, the best-fit amax and ∑dust are indicated with the solid cyan lines in panels a and b, respectively, whereas the dashed red lines refer to the reference profiles convolved with a same beam (i.e., 0.05″) to that of the synthetic images. The gray color scale represents the probability distribution of parameter values returned from fitting the 5000 randomly generated intensity profiles (see Sec. 3). The vertical bar marks the typical confidence interval of the best-fit parameter set. In panels c and d, the reference profiles of intensity and spectral index are shown with dashed lines, while the best-fit results are indicated with solid lines. Panels e–h present the results for model 2. |
4.2 Fitting to the 0.8, 1.3 and 3.0 images simultaneously
We tested the accuracy of the fitting by considering three images, with the results shown in Fig. 6. We selected the 0.8, 1.3, and 3 mm images as the observable, because they correspond to Band 7, Band 6, and Band 3 ALMA observations that have commonly been carried out for protoplanetary disk (e.g., Pascucci et al. 2016; Andrews et al. 2018; Öberg et al. 2021). For model 1, compared to the fitting with only two images, the retrieval for amax in the gap and outer disk is significantly improved, and uncertainties of amax also decrease. The dust surface densities in the ring and outer disk are better recovered. The improvements in the parameter retrieval for model 2 are less than those for model 1.
We also considered a combination of 1.3, 3.0 and 7.8 mm images. The results are similar to those obtained by fitting the images of the three short wavelengths, but the constraints on the dust surface density are overall better since the longest wavelength of 7.8 mm yields the lowest optical depth.
4.3 Fitting to the 0.8, 1.3, 3.0 and 7.8 mm images simultaneously
As a further experiment, we simultaneously fit the four images. The results are presented in Fig. 7. For model 1, the best-fit amax(R) closely match the reference values at all radii. The quality of retrieval for ∑dust is also high. For model 2, the results are also satisfactory, though the best-fit model slightly overpredicts amax(R). It is evident that regardless of whether amax crosses the boundary between the small and large size regions, fitting images at four wavelengths constrain the dust grain size well. This can be understood because the number of data points is larger than that of the free parameters. One can derive only one spectral index with two wavelengths. Three spectral indices can be obtained by using three images. When four images are considered, there will be six different combinations. Figure 1 shows three spectral indices as an illustration, β(0.8–1.3 mm), β(1.3–3 mm), and β(3–7.8 mm), which vary differently with amax. This highly reduces the model degeneracy.
In the inner disk (i.e., R ~ 10 AU) and the ring, the best-fit amax and ∑dust of the 5000 Monte Carlo models still span a large range of values, similar to the outcome of model 1 when fitting with two images (see Fig. 5). Our work shows that the Monte Carlo approach described in Sec. 3 helps to obtain more accurate grain sizes and dust surface densities, because fitting to the observable only once probably returns amax and ∑dust profiles that fluctuate significantly with radius (see Fig. 12 in Sierra et al. 2021 and Fig. 4 in Zhang et al. 2023, for instance). The strong fluctuation is induced by the model degeneracy. Especially, when the observed spectral index, αmm, corresponds to an opacity slope, βmm, close to ~1.8 (see Fig. 1), there exists two solutions; that is, the low-surface-density large size solution and high-surface-density small size solution (see panel a of Fig. 4). With the Monte Carlo approach, we fit the intensities that were generated by randomly adding Gaussian noise to the reference intensity profiles, and obtained a result that favors the reference setup.
![]() |
Fig. 7 Same as in Fig. 5, but for simultaneously fitting to the 0.8 mm, 1.3 mm, 3.0 mm, and 7.8 mm images. |
![]() |
Fig. 8 Same as in Fig. 5, but for investigating the effect of disk inclination on the parameter retrieval. In the fitting procedure, we adopted the setup of model 1, and took all four images into account. |
4.4 Effects of the disk inclination
The above experiments were all conducted for face-on disks. In this section, we investigate the effect of disk inclination on the fitting results by tilting model 1 with different angles. Figure 8 displays the results for the inclination of 20°, 40°, and 60°, respectively. When the disk inclination is low – for example, ≲ 20° – the fitting quality is minimally affected, and both amax and ∑dust can be well retrieved. As the inclination increases to 40°, the best-fit profile for both parameters starts to deviate from the reference profile, with the gap-ring contrast being shallower than the reference. When the disk is highly inclined – for example, i = 60° – the retrieved results are not reliable for most parts of the disk, though radial intensity profiles can still be well recovered at all wavelengths.
Constraining amax and ∑dust for disks with high inclinations is challenging. The observed intensity is an integration of dust emission along the line of sight. Therefore, substructures will be smoothed out when the disk is tilted, leading to smooth intensity profiles rather than showing distinct gap and ring features. Accordingly, the retrieved profiles for amax and ∑dust are also flat. For face-on disks, the constrained amax and ∑dust trace the mean disk properties over a radius range determined by the beam size; 5 AU in our case. However, for inclined disks, the radius range for averaging is larger due to the inclination effect, making the parameter retrieval more complex and difficult. Therefore, special attention should be paid in the investigation of dust grain size in disks with high inclinations.
Parameter of reference models for smooth disks.
4.5 Fitting to reference models of smooth disks
The reference models used for the above tests have a gap and a ring (see Eqs. (3) and (4)). We also examined the retrieval process using three smooth disks without substructures. For this purpose, we prescribed the dust surface density as ∑(R) = ∑0 R−0.5. The three disks differ in terms of the radial grain size profile given by amax(R) = a0 R−γ. Table 2 tabulates the model parameters. For each of the three reference disks, we fit the 0.8, 1.3, 3.0, and 7.8 mm images simultaneously. The results are presented in Fig. 9.
These additional experiments demonstrate again that the analytic approach recovers the grain size well when four images are considered in the fitting process. As the disk inclination increases, the best-fit results appear inconsistent with the reference setup, particularly for disks with a large γ. Since a0 is fixed, a steeper grain size profile means that more disk regions are populated with amax < 0.1 mm. Figure 9 shows that it is difficult to recover grain sizes smaller than 0.1 mm. The reason is that all of the three βmm are close to ~1.8, and no longer sensitive to the grain size in this regime (see Fig. 1).
![]() |
Fig. 9 Same as in panels a and b of Fig. 5; however, the reference setup does not include substructures, mimicking smooth disks. The three reference models, denoted as ms1, ms2, and ms3, differ in terms of the radial grain size profile (see Sec. 4.5). |
5 Summary
Investigating the grain size and its correlation with disk substructures is crucial to understand the initial step of planet formation. Constraining the grain size from millimeter observations needs to self-consistently solve the radiative transfer equation, which is time-consuming and not efficient for parameter studies with large sample sizes. As an alternative, analytic models under the assumption of vertically isothermal disks have been proposed and widely used in the literature.
In this work, we created reference disk models with known properties, and simulated millimeter continuum images as synthetic observations. We fit the synthetic images using analytic models to investigate the applicability of the approach, and examined the conditions under which the input grain sizes can be recovered. Our conclusions are given as follows:
- (1)
Fitting images at only two wavelengths is not sufficient to retrieve the reference setup, mainly due to the degeneracy between the small and large size solutions (see Fig. 1). Fitting three images improves the retrieval of grain size; however, the obtained dust surface density still differs a lot from the reference model. The analytic models work well when all of the four images are considered in the fitting process. By combining each pair of images from the four, one can derive six spectral indices that behave differently as a function of amax (see three of the six spectral indices in Fig. 1). This highly reduces the model degeneracy.
- (2)
The inclination angle has a significant impact on the fitting results. When the inclination is low – for example, ≲ 40° – the reference grain sizes can be recovered well using the analytic models. However, when the disk is highly tilted, both the grain size and the dust surface density cannot be constrained. This is because the inclination effect will smooth out any substructure in the radial intensity profile of the disk, and the optical depth effect needs to be treated more properly.
- (3)
The uncertainties of the best-fit parameter set are relatively large in some regions; for example, in the inner disk and the ring. The Monte Carlo method introduced in Sec. 3 helps to obtain more accurate results.
Acknowledgements
We thank the anonymous referee for the constructive comments that improve our manuscript. Y.L. acknowledges financial supports by the National Natural Science Foundation of China (Grant Number 11973090), and by the International Partnership Program of Chinese Academy of Sciences (Grant Number 019GJHZ2023016FN). HW acknowledges the financial support by the National Natural Science Foundation of China (Grant Number 11973091).
References
- Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progr. Theor. Phys., 56, 1756 [Google Scholar]
- Andrews, S. M. 2020, ARA&A, 58, 483 [Google Scholar]
- Andrews, S. M., & Williams, J. P. 2005, ApJ, 631, 1134 [Google Scholar]
- Andrews, S. M., & Williams, J. P. 2007, ApJ, 671, 1800 [CrossRef] [Google Scholar]
- Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502 [CrossRef] [Google Scholar]
- Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
- Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [NASA ADS] [CrossRef] [Google Scholar]
- Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge University Press) [Google Scholar]
- Bae, J., Isella, A., Zhu, Z., et al. 2023, in ASP Conf. Ser., 534, 423, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [NASA ADS] [Google Scholar]
- Beckwith, S. V. W., & Sargent, A. I. 1991, ApJ, 381, 250 [NASA ADS] [CrossRef] [Google Scholar]
- Birnstiel, T. 2023, ARA&A, submitted [arXiv:2312.13287] [Google Scholar]
- Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41 [Google Scholar]
- Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 [CrossRef] [Google Scholar]
- Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Bruggeman, D. A. G. 1935, Ann. Phys., 416, 636 [Google Scholar]
- Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71 [Google Scholar]
- Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806 [Google Scholar]
- Davis, S. S. 2005, ApJ, 627, L153 [NASA ADS] [CrossRef] [Google Scholar]
- Doi, K., & Kataoka, A. 2023, ApJ, 957, 11 [CrossRef] [Google Scholar]
- Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [Google Scholar]
- Draine, B. T. 2006, ApJ, 636, 1114 [Google Scholar]
- Drążkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, in ASP Conf. Ser., 534, 717, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [NASA ADS] [Google Scholar]
- Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957 [NASA ADS] [CrossRef] [Google Scholar]
- Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool [Google Scholar]
- Francis, L., Johnstone, D., Herczeg, G., Hunter, T. R., & Harsono, D. 2020, AJ, 160, 270 [NASA ADS] [CrossRef] [Google Scholar]
- Guidi, G., Tazzari, M., Testi, L., et al. 2016, A&A, 588, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97 [Google Scholar]
- Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42 [NASA ADS] [CrossRef] [Google Scholar]
- Kenyon, S. J., & Hartmann, L. 1987, ApJ, 323, 714 [Google Scholar]
- Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211 [Google Scholar]
- Kurucz, R. 1994, Solar abundance model atmospheres for 0, 19 [Google Scholar]
- Lesur, G., Ercolano, B., Flock, M., et al. 2022, arXiv e-prints [arXiv:2203.09821] [Google Scholar]
- Li, A., & Draine, B. T. 2001, ApJ, 554, 778 [Google Scholar]
- Li, D., Liu, Y., Wang, H., Wang, Y., & Ma, Y. 2023, MNRAS, 518, 6092 [Google Scholar]
- Liu, H. B. 2019, ApJ, 877, L22 [NASA ADS] [CrossRef] [Google Scholar]
- Liu, Y., Henning, T., Carrasco-González, C., et al. 2017, A&A, 607, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17 [Google Scholar]
- Macías, E., Guerra-Alvarado, O., Carrasco-González, C., et al. 2021, A&A, 648, A33 [EDP Sciences] [Google Scholar]
- Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in ASP Conf. Ser., 534, 501, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [NASA ADS] [Google Scholar]
- Miyake, K., & Nakagawa, Y. 1993, Icarus, 106, 20 [NASA ADS] [CrossRef] [Google Scholar]
- Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [Google Scholar]
- Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1 [CrossRef] [Google Scholar]
- Ohashi, S., Momose, M., Kataoka, A., et al. 2023, ApJ, 954, 110 [NASA ADS] [CrossRef] [Google Scholar]
- Paardekooper, S.-J., Dong, R., Duffell, P., et al. 2022, arXiv e-prints [arXiv:2203.09595] [Google Scholar]
- Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125 [Google Scholar]
- Pérez, L. M., Carpenter, J. M., Chandler, C. J., et al. 2012, ApJ, 760, L17 [CrossRef] [Google Scholar]
- Pinilla, P., Benisty, M., & Birnstiel, T. 2012a, A&A, 545, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012b, A&A, 538, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pinilla, P., de Juan Ovelar, M., Ataiee, S., et al. 2015a, A&A, 573, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pinilla, P., van der Marel, N., Pérez, L. M., et al. 2015b, A&A, 584, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [Google Scholar]
- Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ricci, L., Testi, L., Natta, A., & Brooks, K. J. 2010a, A&A, 521, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ricci, L., Testi, L., Natta, A., et al. 2010b, A&A, 512, A15 [CrossRef] [EDP Sciences] [Google Scholar]
- Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173 [Google Scholar]
- Rota, A. A., Meijerhof, J. D., van der Marel, N., et al. 2024, A&A, 684, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Sierra, A., Lizano, S., Macías, E., et al. 2019, ApJ, 876, 7 [Google Scholar]
- Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14 [NASA ADS] [CrossRef] [Google Scholar]
- Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344 [NASA ADS] [CrossRef] [Google Scholar]
- Tazzari, M., Testi, L., Ercolano, B., et al. 2016, A&A, 588, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Tazzari, M., Testi, L., Natta, A., et al. 2017, A&A, 606, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Tazzari, M., Testi, L., Natta, A., et al. 2021, MNRAS, 506, 5117 [NASA ADS] [CrossRef] [Google Scholar]
- Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018a, ApJ, 860, L12 [NASA ADS] [CrossRef] [Google Scholar]
- Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018b, ApJ, 868, 113 [Google Scholar]
- Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339 [Google Scholar]
- Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130 [CrossRef] [Google Scholar]
- Trotta, F., Testi, L., Natta, A., Isella, A., & Ricci, L. 2013, A&A, 558, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
- Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211 [Google Scholar]
- Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [Google Scholar]
- Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 [NASA ADS] [CrossRef] [Google Scholar]
- Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7 [Google Scholar]
- Zhang, S., Zhu, Z., Ueda, T., et al. 2023, ApJ, 953, 96 [NASA ADS] [CrossRef] [Google Scholar]
- Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
![]() |
Fig. 1 Dust opacity slope, βmm, as a function of the maximum grain size, amax. The slopes between two wavelengths given in the legend are derived by calculating the mass absorption coefficients (κabs). To calculate κabs, we adopted the MIE scattering theory, assuming the DSHARP dust composition (see Sec. 2.3). The grain size distribution follows a power law, N(a) ∝ a−3.5, with the minimum dust size fixed to amin = 0.01 μm. The vertical dashed line marks the 0.4 mm boundary that divides the amax values into two regimes (a small size region and a large size region) according to the β(1.3–3 mm) curve. |
In the text |
![]() |
Fig. 2 Properties of the reference models. Top panel: dust surface density. The dust surface densities of model 1 and model 2 are the same. Middle panel: radial profile of the maximum grain size of model 1 (red line) and model 2 (blue line). The horizontal dashed line marks the 0.4 mm boundary defined in Fig. 1. Bottom panel: optical depth at the wavelength of 1.3 mm. The color coding is identical to that in the middle panel. The τ1.3mm = 1 condition is indicated with the horizontal dashed line. |
In the text |
![]() |
Fig. 3 Synthetic images of model 1 at wavelengths of 0.8 mm, 1.3 mm, 3.0 mm, and 7.8 mm. The images are observed with a face-on orientation. The beam with a size of 0.05″ is indicated with the white dot in the lower left corner of each panel. |
In the text |
![]() |
Fig. 4 Fitting process for model 1 at R = 5 AU. The normalized joint probability distribution is shown as the color map in panel a. The marginalized probabilities for amax and ∑dust are indicated with the black lines in panels b and c, respectively. The red star denotes the best-fit parameter set identified from the fitting procedure. |
In the text |
![]() |
Fig. 5 Results of fitting to the 1.3 mm and 3.0 mm images for model 1 (left column) and model 2 (right column). For model 1, the best-fit amax and ∑dust are indicated with the solid cyan lines in panels a and b, respectively, whereas the dashed red lines refer to the reference profiles convolved with a same beam (i.e., 0.05″) to that of the synthetic images. The gray color scale represents the probability distribution of parameter values returned from fitting the 5000 randomly generated intensity profiles (see Sec. 3). The vertical bar marks the typical confidence interval of the best-fit parameter set. In panels c and d, the reference profiles of intensity and spectral index are shown with dashed lines, while the best-fit results are indicated with solid lines. Panels e–h present the results for model 2. |
In the text |
![]() |
Fig. 6 Same as in Fig. 5, but for simultaneously fitting to the 0.8 mm, 1.3 mm, and 3.0 mm images. |
In the text |
![]() |
Fig. 7 Same as in Fig. 5, but for simultaneously fitting to the 0.8 mm, 1.3 mm, 3.0 mm, and 7.8 mm images. |
In the text |
![]() |
Fig. 8 Same as in Fig. 5, but for investigating the effect of disk inclination on the parameter retrieval. In the fitting procedure, we adopted the setup of model 1, and took all four images into account. |
In the text |
![]() |
Fig. 9 Same as in panels a and b of Fig. 5; however, the reference setup does not include substructures, mimicking smooth disks. The three reference models, denoted as ms1, ms2, and ms3, differ in terms of the radial grain size profile (see Sec. 4.5). |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.