© The Authors 2025

1 Introduction

The dust content of protoplanetary disks sets the most basic initial conditions for planet formation. The total dust mass in disks represents the budget available to form planetesimals, rocky planets, and the cores of giant planets (Drazkowska et al. 2023). Its evolution with time is a proxy for the timescales of disk evolution and planet formation. Moreover, the disk mass that can be estimated from the dust mass (by assuming a gas-to-dust ratio) offers insights into which dynamical processes are ongoing during planet-forming stages (e.g., gravitational instabilities; Longarini et al. 2022).

Unfortunately, this key ingredient in our planet formation theories is still highly uncertain. The classical mass measurement technique, based on integrated fluxes at (sub)millimeter wavelengths (Hildebrand 1983), leads to dust masses of Class II disks that do not reach the solid mass found in exoplanetary systems (Manara et al. 2018). While this could suggest that early- stage formation (Class I disks; Tychoniec et al. 2020) and/or latestage infall are at play (Gupta et al. 2024), new high-resolution observations hint that these values are severely underestimated (Liu et al. 2022). High-resolution observations from facilities such as the Atacama Large Millimeter/submillimeter Array (ALMA) or the Spectro-Polarimetric High-Contrast Exoplanet Research (SPHERE) instrument at the Very Large Telescope (VLT) (e.g., Andrews et al. 2018, Garufi et al. 2024) reveal that protoplanetary disks are characterized by various dust substructures (Andrews 2020), most of which concentrate millimeter- to centimeter-sized particles. The optical thickness of these dense regions at (sub)millimeter wavelengths could be hiding substantially large amounts of dust, leading to low measured dust masses. Although this flux-based mass measurement technique is easy to apply to a large sample of disks, it is unfit for accurately assessing the dust content of protoplanetary disks.

Another key missing piece in planetary formation theory is how millimeter- to centimeter-sized pebbles grow into kilometersized planetesimals. Planetesimal formation via direct growth is severely limited by fragmentation and radial drift (Drazkowska et al. 2023). The radial dust redistribution in disks is driven by aerodynamic drag from the surrounding gas, acting as a headwind that slows down particles, causing them to drift rapidly toward the disk center (Whipple 1972). Immersed in a turbulent gas medium, the relative velocities of grains increase, eventually leading to fragmentation rather than accretion (Birnstiel et al. 2012).

The ubiquitous detection of substructures in high-resolution ALMA observations suggests a way to ease this rapid depletion of a disk’s dust reservoir. If these substructures result from pressure maxima, they can act as dust traps, slowing drifting particles (Pinilla et al. 2012). In these dense dust regions, the dust’s back- reaction on the gas slows the local gas velocity, reducing the headwind experienced by the dust. For sufficiently high dust- to-gas ratios (Li & Youdin 2021; Lim et al. 2024), this positive feedback leads to the formation of self-gravitating clumps of particles in the mid-plane that rapidly grow into kilometer-sized planetesimals. This process, known as the streaming instability, is currently the leading theory for planetesimal formation (Youdin & Goodman 2005). By measuring the level of grain growth and dust density in substructures, we can test if the physical criteria required to trigger streaming instabilities are met (Scardoni et al. 2021).

Polarization observations at (sub)millimeter wavelengths have been proposed as a method for measuring grain sizes in disks. Recent measurements from ALMA suggest submillimeter maximum dust sizes (Kataoka et al. 2017), in contrast to the (sub)centimeter maximum particle sizes predicted by simulations of coagulation or fragmentation dust evolution with dust trapping (Birnstiel et al. 2012). However, different interpretations of the polarization pattern and fraction yield differing outcomes that may reconcile these discrepancies (e.g., a double dust population, Kataoka et al. 2017; particle porosity, Zhang et al. 2023). At present, the complex interpretation of these polarization observations hinder a clear determination of the maximum grain size.

Alternatively, multi-wavelength continuum observations of the dust thermal emission at (sub)millimeter wavelengths remain one of the most efficient tools to study both the particle size distribution and the dust surface density in protoplanetary disks. The spectral energy distribution (SED) of the thermal dust emission of disks depends on the dust temperature and, if optically thin detections are included, on the surface density and on the maximum particle size through the frequency dependence of the dust opacity (see Sect. 2.2). The few bright sources that have been observed at high resolution at various wavelengths all show evidence of accumulations of pebbles in substructures, in agreement with our current view on planet formation (e.g., Carrasco-González et al. 2019; Sierra et al. 2021; Macías et al. 2021). To confirm these preliminary findings, more detailed multi-wavelength studies of a larger unbiased sample of disks are necessary. However, differences in sensitivity and resolution achievable across various frequencies by advanced facilities, such as ALMA and the Jansky Very Large Array (VLA), currently pose significant limitations to this goal.

High-resolution disk observations outside the 211 GHz– 373 GHz range (ALMA Bands 6 and 7) are uncommon due to the challenges of achieving both high resolution and sensitivity at lower frequencies, where the emission is fainter (ALMA Bands 3 and 4), and at higher frequencies, where atmospheric conditions are more demanding (ALMA Bands 8, 9, and 10). Moreover, both high-resolution multi-wavelength analyses (e.g., Carrasco-González et al. 2019; Macías et al. 2021) and population synthesis models (Delussu et al. 2024)) have shown that dense disk regions, such as the inner disk or rings, can remain marginally optically thick even at wavelengths up to 3 mm (ALMA Band 3). Unfortunately, optically thin longer- wavelength observations currently offer coarser resolution than required to resolve substructures (<0.03″–0.05″). ALMA Band 1 (6–8.6 mm) data can achieve a maximum resolution of 0.1″ . Although the VLA can provide somewhat higher resolutions (0.05″–0.06″ at 7–9 mm), its observations are of lower quality compared to ALMA, due to the VLA’s limited UV coverage and poorer phase stability. Achieving sufficient sensitivity requires extensive observation times (e.g., ∼32 hours for the 9 mm VLA image of HL Tau; Carrasco-González et al. 2019). Additionally, joint ALMA and VLA observations are limited to a small sample of sources in a narrow region of the sky. As a result, only a limited number of SED analyses of disks incorporate high-resolution optically thick and thin data, which are crucial for accurate measurements of dust properties.

All of the previously mentioned limitations make multifrequency continuum observations among the most timeconsuming methods for studying protoplanetary disks. Thus, it is fundamental to investigate the most economic way of performing this dust characterization analysis.

For this work, we addressed this question by benchmarking the state-of-the-art dust characterization analysis on simulated multi-wavelength observations of protoplanetary disks. First, our aim is to identify the optimal combination of multi-wavelength observations that achieves accurate results, while minimizing observational demands on telescope resources. Since the best multi-band setup depends on the optical depth of the observations, we varied the total dust content of our model disks while performing this analysis. Additionally, we thoroughly tested how the differences in quality and spatial resolution and/or the absence of optically thin observations, which characterize the majority of the available spatially resolved SED analyses, affect the determination of the dust content of pro- toplanetary disks. Finally, we propose an innovative approach that enables the inference of dust properties at high angular resolution using economical medium-resolution ALMA Band 1 observations. This method not only exploits ALMA’s higher quality to capture even fainter objects compared to VLA, but also provides a more cohesive dataset for multi-band analyses by ensuring a consistent field of view across observations. Together, these advantages allow us to extend dust characterization to a larger sample of disks, enhancing the statistical robustness of our findings.

In Sect. 2, we describe the simulated multi-frequency dust thermal emission of disks. In Sect. 3 we introduce the state- of-the-art multi-wavelength SED analysis. The results of benchmarking this technique on simulated observations are shown in Sect. 4. In particular, in Sect. 4.1 we discuss how the dust thermal emission at various frequencies aid the reconstruction of the dust properties from an SED analysis. Instead, in Sect. 4.2 we discuss how the limited resolution of real observations affects the measurement of the dust properties. Finally, in Sect. 5, we propose new approaches to the dust characterization technique that allows us to extend this method to a larger and unbiased sample of disks.

2 Simulated multi-wavelength observations

To quantify the effectiveness of the current SED analysis procedure in constraining dust properties from real observations, we applied this technique to simulated multi-wavelength observations of disks.

We designed realistic observations of the dust thermal emission of a protoplanetary disk by first evolving an interstellar medium (ISM) dust population over 1 Myr. This evolution accounts for the coagulation and fragmentation processes expected within the disk, as well as dust transport (see Sect. 2.1). Then, we analytically computed the multi-band thermal emission of this evolved dust population, assuming that it can be described as a 1D vertically isothermal slab (see Sect. 2.2). We included in this analysis the emission at 0.45 mm (ALMA Band 9), 0.87 mm (ALMA Band 7), 1.29 mm (ALMA Band 6), 3.07 mm (ALMA Band 3), and 7.46 mm (ALMA Band 1 or VLA Q Band). We used Band 9 observations for optically thick detections of dust thermal emission, as Band 10 requires very strict weather conditions for observations. While including both Band 7 and 6 might not be necessary since they have similar expected optical depths, the combination of both these highly requested observations mitigates the uncertainties coming from the noise and the flux calibration errors. The emission at both Band 4 and Band 3 is expected to be marginally thin but, because of the slightly lower frequency, Band 3 observations are preferred. Both Band 8 and Band 5 detections are less common than Band 9/7 or Band 6/4. For the long-wavelength observations, we considered a frequency of 40 GHz, which can be observed by both ALMA and the VLA, enabling comparative analysis. Similar results can be obtained with the 9 mm emission observed by the VLA in Ka Band.

We designed two disk models: a compact disk (with a radius after 1 Myr of disk evolution of R_68% ~ 0.15″ = 22 au and R_99% ~ 0.39″ = 55 au at 140 pc, in Band 7 continuum emission), with a low contrast and small scale ringed morphology that resembles TW Hya (Macías et al. 2021); and an extended (with a radius after 1 Myr of disk evolution of R_68% ~ 0.35″ = 50 au and R_99% ~ 0.82″ = 115 au at 140 pc, in Band 7 continuum emission) double-ringed disk, similar to HD 163296 (Guidi et al. 2022). We simulated multi-wavelength radio-interferometric observations of these observations at three different resolutions: 0.05″, which represents the maximum resolution currently achievable at 7 mm by the VLA and at 3 mm by ALMA; 0.1″, the limit of ALMA Band 1 with ALMA’s most extended configuration; and 0.2″ , a practical resolution for most disks (see Sect. 2.3).

2.1 Simulated dust distributions

We employed the open-source code DustPy (Stammler & Birnstiel 2022) to simulate the evolution of the dust mass distribution in the two disks with substructures we are modeling. This code solves the Smoluchowski coagulation equation (Smoluchowski 1916), while accounting for dust fragmentation and transport. Following the results of Jiang et al. (2024), we assumed a population of fragile pebbles with a fragmentation velocity of 1 m/s. For the initial grain size distribution, we adopted the default MRN distribution observed by Mathis et al. (1977) in the interstellar medium (ISM): n(a) da ∝ a^−q da, with q = 3.5, a minimum grain size of 0.005 µm, and a maximum grain size of 1 µm.

The gas surface density evolves with the master equation of viscous accretion (Lynden-Bell & Pringle 1974), while the dust surface density evolution is solved by the advection-diffusion equation (Birnstiel et al. 2011). The relative velocities of dust particles prior to collisional growth accounts for thermal Brownian motion, vertical stirring and settling, turbulent mixing, as well as azimuthal and radial drift.

We initialized the gas surface density as the similarity solution of Lynden-Bell & Pringle (1974), $${\text{Σ}}_{\text{g}}={\text{Σ}}_{\text{g},0}{\left(\frac{r}{r_{\text{c}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{r_{\text{c}}}\right)}^{2-\gamma }\right],$$(1)

with γ = 1 and a characteristic radius r_c = 100 au. For both models, we assumed an initial gas mass of 0.05 M_⊙ and an initial gas-to-dust ratio of 100. The disks were assumed to be vertically isothermal. This is a good approximations since we aim to analyze the (sub)millimeter emission of disks, which is tracing the geometrically thin midplane layer of the disk where large dust grains have settled (Pinte et al. 2016). We defined the temperature profile as $$T_{\text{dust}}=T_0{\left(\frac{r}{r_{\text{c}}}\right)}^{-0.45},$$(2)

with the typically assumed T₀ = 20 K at r_c = 100 au. To simulate the dust evolution in disks with substructures, we introduced pressure traps by radially modifying the gas viscosity as (Dullemond et al. 2018) $${\alpha }_{v,r}=\frac{{\alpha }_v}{F(r)},F(r)=\exp \left[-A\exp \left(-\frac{{\left(r-r_0\right)}^2}{2w^2}\right)\right],$$(3)

with α_ν = 10⁻⁴. The radial mixing (δ_r), turbulent mixing (δ_t) and vertical mixing (δ_υ) parameters are fixed to a value of 10⁻⁴. Our extended disk model is characterized by two deep gaps (A = 2.0) at positions r₀ = 50 au and r₀ = 100 au with a width of w = 1 au. On the contrary, to obtain the compact disk model, we introduced two shallow gaps (A = 0.1 au) at half the distance from the star (r₀ = 25 au and 50 au) and with half the width (w = 0.5 au).

The results of the DustPy simulations after 1 Myr are shown in Figure 1. The top panel shows the grain size distribution as a function of the radial distance from the star for the compact disk (left panel) and the extended disk (right panel). At each radius, what we identified as the maximum grain size a_max(r), which is the maximum grain size that has nonnegligible dust density, is highlighted in lime green. Here we defined it as the grain size corresponding to a dust mass density in the 10⁻⁶–10⁻⁸ g/cm³ range. This range was selected through visual inspection of Figure 1, as lower densities do not reflect the expected conditions in the midplane of a protoplanetary disk. We note that this choice is somewhat arbitrary but a different choice of threshold would simply move the maximum grain size profile slightly up or down, which does not affect the aim or results of this paper. In the bottom row, the corresponding gas and dust surface densities distributions are shown. The vertical gray lines represent the position of the bright (B) and dark (D) rings in the dust surface density profile.

The maximum grain size in the compact disk is ∼3 mm in the inner disk. It reaches a value of 1.5 mm in the inner ring before decreasing to submillimeter sizes in the outer disk. Instead, the maximum grain sizes at the two rings of the extended disk are respectively: ∼2.5 mm and 0.5 mm. This is in agreement with the results of SED analyses of the handful of T Tauri stars for which we have high-resolution multi-wavelength observations (Sierra et al. 2021; Carrasco-González et al. 2019; Tazzari et al. 2016).

Fig. 1 DustPy outputs at 1 Myr for the compact (left panel) and the extended disk models (right panel). Top panel: grain size distribution of the dust particles as a function of the radial distance from the star. At each radius, the maximum grain size, is highlighted in lime green. Bottom panel: radial profiles of gas (in blue) and dust (in orange) total column densities. The vertical gray lines represent the position of the bright (B) and dark (D) gaps in the dust surface density profile.

2.2 Analytical radiative transfer

To compute the continuum emission at (sub)millimeter wavelengths of the simulated grains population of our model disks, we assumed the disks to be azimuthally symmetric, vertically isothermal, and razor-thin. With these assumptions, the dust thermal emission from the disk mid-plane can be computed as a 1D vertically isothermal slab (Miyake & Nakagawa 1993; Sierra et al. 2019), $$I_v=B_v\left(T_{\text{dust }}\right)\left[\left(1-{\text{e}}^{-{\tau }_v/\mu }\right)+{\omega }_{v}F\left({\tau }_v,{\omega }_v\right)\right],$$(4)

where B_ν(T_dust) is the blackbody emission at the dust temperature T_dust and frequency v. τ_ν = Σ_dust χ_ν is the optical depth, defined as the product of the dust surface density Σ_dust and the total dust opacity χ_ν .

At millimeter wavelengths, emission due to dust scattering is not negligible. On the contrary, dust scattering can considerably reduce the intensity of an optically thick disk (Zhu et al. 2019). Thus, here we considered the total dust absorption χ_ν as the sum of the absorption opacity $k_v^{\text{a}}$ and the effective scattering opacity $k_v^{\text{s},\text{eff}}$. The effective scattering opacity is an approximation for the anisotropic scattering (Henyey & Greenstein 1941; Birnstiel et al. 2018) and is related to the scattering opacity $\left(k_v^s\right)$ as: $k_v^{\text{s,eff }}=(1-{𝑔}_v)k_v^{\text{s}}$ , where ɡ_ν is the forward scattering parameter. The albedo is defined as ${\omega }_v=k_v^{\text{s}}/\left(k_v^{\text{a}}+k_v^{\text{s}}\right)$. Lastly, µ = cos(i) is the cosine of the disk inclination angle i. Here, we assumed face-on (µ = 1) disks.

Since scattering is taken into account, the solution for the radiative transfer also depends on $$\begin{array}{l}F\left({\tau }_v,{\omega }_v\right)=\frac{1}{\exp \left(-\sqrt{3}{ϵ}_v{\tau }_v\right)\left({ϵ}_v-1\right)-\left({ϵ}_v+1\right)}\\ \times \left[\frac{1-\exp \left(-\left(\sqrt{3}{ϵ}_v+1/\mu \right){\tau }_v\right)}{\sqrt{3}{ϵ}_v\mu +1}+\frac{\exp \left(-{\tau }_v/\mu \right)-\exp \left(-\sqrt{3}{ϵ}_v{\tau }_v\right)}{\sqrt{3}{ϵ}_v\mu -1}\right],\end{array}$$(5)

where ${ϵ}_v=\sqrt{1-{\omega }_v}.$

Given a dust temperature radial profile T_dust(r) and the absorption and scattering opacities $k_v^{\text{a}}$ and $k_v^{\text{s}}$ , the continuum intensity in Equation (4) can be computed for the two dust surface density models (Σ_dust ) at all wavelengths. We defined T_dust(r) as in Equation (2). While this equation does not account for the thermal behavior of gaps and rings, performing a full radiative transfer calculation would be computationally expensive and would require detailed information on the vertical structure of disks, which is not simulated in DustPy. Additionally, a more physical description would not significantly affect the scope of our work.

Since the dust population in disks consists of a continuous size distribution instead of a single size, as highlighted by both real observations (Avenhaus et al. 2018) as well as models of grain growth with fragmentation (Weidenschilling 1984; Dullemond & Dominik 2005; Brauer et al. 2008; Birnstiel et al. 2010), we averaged at each radius the size-dependent opacities $k_v^{\text{a}}\left(a\right)$ and $k_v^{\text{s}}\left(a\right)$ over the particle mass distribution: $$k_v^{\text{a}/\text{s},\mathrm{tot}}(r)=\frac{{\displaystyle {\int }_{a_{\min }(r)}^{a_{\max }(r)}n}(a)m(a)k_v^{\text{a}/\text{s}}(a)\text{d}a}{{\displaystyle {\int }_{a_{\min }(r)}^{a_{\max }(r)}n}(a)m(a)\text{d}a}.$$(6)

In the previous equation m(a) is the mass of a grain of size a, amin/max (r) are the minimum/maximum grain sizes at each radius and n(a) da is the number of grains with sizes between a and a + da. To reduce computational costs and better isolate the effects of the other variables in the analysis, we again adopted the commonly used power-law approximation for the grain size distribution: n(a) da ∝ a ^−q da with q = 3.5, which is consistent with observations in protoplanetary disks (Doi & Kataoka 2023). We note that the resulting size distribution in the DustPy simulations is, in any case, similar to a power-law with q = 3.5, so the effects of this assumption are small. With this choice of q, the value of a_min(r) has little effect on the opacity at millimeter wavelengths and we fixed it to 10⁻⁵ cm at all radii. At each radius, a_max (r) is defined as the maximum grain size that has nonnegligible dust density (see Sect. 2.1).

To evaluate the dust opacities $k_v^{\text{a}}\left(a\right)$ and $k_v^{\text{s}}\left(a\right)$, as well as the particle masses, we adopted the “DSHARP” dust composition: a compact mixture of water ice, troilite, refractory organics, and astronomical silicates (Birnstiel et al. 2018). We computed the opacities with the Python package dsharp_opac. We do not discuss different dust opacities in this work but we emphasize how different dust materials, such as amorphous carbons/silicates results in very different opacities and dust properties, first of all dust masses (Zagaria et al., in prep.).

We note that, despite the approximations required to derive them, Zhu et al. (2019) showed that both the Equations (4) and (5) are in excellent agreement with more complex and time-consuming radiative transfers simulations. Thus, we follow this analytical approach without testing more computational expensive radiative transfer models.

2.3 Simulated observations

Through Equation (4), we evaluated the intensity profiles of the dust continuum emission of both the compact and extended disk model at the following wavelengths: 0.45 mm, 0.87 mm, 1.29 mm, 3.07 mm and 7.46 mm. We also produce their 2D emission map, after assuming axisimmetry. To simulate the effect of the finite resolution of real observations on our models, we then applied to the multi-frequency 2D emission maps the simobserve task of CASA (CASA Team 2022).

Given an input sky model, simobserve accounts both for discrete uv-sampling of the ALMA/VLA configurations and the thermal noise. Given the different requirements on sensitivity of our multi-wavelengths and multi-facility observations, we neglected the thermal noise in our images. By doing this, we assumed that the emission at all wavelengths is detected with high quality and a high signal-to-noise ratio, and the flux calibration error is the dominant source of uncertainties. Further comments about the limitations introduced in the SED analysis by limited sensitivity can be found in Sect. 5.2.

To test how the limited resolution and the beam smearing affect the results of the SED characterization, we produced mock observations with three different resolutions: 0.05″, the maximum resolution achievable by the VLA in Band Q; 0.1″ , the maximum resolution of ALMA Band 1; and 0.2″ , a reasonable and time-efficient resolution for ALMA Band 1 observations of disks.

We chose to position our mock sources at the distance (140 pc) and at the coordinates of HL Tau. At this declination (around +18:13:57), sources can be observed with a contained beam elongation by both ALMA and VLA.

The ALMA and VLA antennas configurations used to simulate the observations at the different frequencies and with different resolution are listed in Table 1. At each wavelength, we observed with both an extended and compact configuration. We follow the ALMA technical handbook for the time multipliers between configurations. We used a time multiplier of 0.25 between VLA configuration B and A. Each extended configuration has been observed for 3 hours in transit.

We imaged the simulated observations with the CASA task tclean. We emploied the mtmfs deconvolver (Rau & Cornwell 2011), assuming that the frequency dependency of the emission follows a Taylor expansion to first terms (i.e., nterms = 2). We used briggs weighting and multiple scales at 0 (point-source), 1, 3, 5 times the beam size. We defined the pixel size to be ∼1/7 of the beam size. Different values of the robust parameter, from −1.0 to 0.5, were explored to produce a clean beam as close as possible to the expected resolution. To minimize any uncertainty due to convolution in the image plane, we also made use of an uv-taper during the cleaning procedure to circularize the beam and approach the expected resolution. Finally, we smoothed the clean image to the chosen resolution (0.05″ or 0.1″ or 0.2″) using the imsmooth task of CASA, assuring that the emissions at different frequencies are all convolved to the same beam.

We obtained azimuthally averaged radial intensity profiles of the mock observations by averaging the emission in concentric annuli (here rings since we assume face-on disks). To minimize correlated information in adjacent annuli, we defined the radial width of each annulus as 1/3 of the beam size. The uncertainty of the radial profiles in each annulus was computed as the standard error of the mean: the standard deviation within the annulus (σ_Annulus (r)) divided by the square root of the number of beams N_Beams . In other words, $$\overline{{\sigma }_{\text{Annulus }}}(r)=\frac{{\sigma }_{\text{Annulus }}(r)}{\sqrt{N_{\text{Beams }}}}=\frac{{\sigma }_{\text{Annulus }}(r)}{\sqrt{A_{\text{Annulus }}/A_{\text{Beams }}}},$$(7)

where A_Annulus and A_Beams are the area of the annulus and the beam respectively.

3 SED analysis

We now aim to treat our synthetic multi-wavelength observations as real data. Specifically, we extracted information about the dust content from our simulations by analyzing their spatially resolved SEDs. We then compared these measured values with the initial prescriptions used in our models, in order to assess the reliability of this procedure usually applied to real data.

Through a Bayesian approach, we compared, at each radius, the mock intensity profiles at various wavelengths to the physical model in Equation (4), which is a function of the dust temperature T_dust(r), surface density Σ_dust(r) and maximum grain size a_max(r). We made the same assumptions on dust composition and optical properties as we did in Sect. 2.2.

We emploied the code emcee (Foreman-Mackey et al. 2019), that implements the affine-invariant Markov-Chain Monte Carlo (MCMC) Ensemble Sampler of Goodman & Weare (2010), to estimate the posterior distribution of the model parameters at each radius. We used a standard log-normal likelihood function $$\ln p\left(I_v(r)\mid \theta \right)=-\frac{1}{2}{\displaystyle \sum _v\left({\left(\frac{I_v(r)-I_v^{\mathrm{model}}(r,\theta )}{{\sigma }_{\text{tot},v}(r)}\right)}^2+\ln \left(2\pi {\sigma }_{\text{tot},v}^2(r)\right)\right)},$$(8)

where I_ν(r) is the azimuthally averaged intensity profile at radius r and $I_v^{\text{model }}(r,\theta )$ is the model intensity in Equation (4) evaluated for the vector of model parameters θ = [T_dust, Σ_dust, a_max].

In this work, we fixed the power-law index of the grain size distribution q, introduced in Equation (6), to 3.5. This is consistent with how we simulated the multi-wavelength dust thermal emission (see Sect. 2.2). Since this parameter is strongly correlated with a_max , we still performed a simple test of the ability of the SED analysis to retrieve both parameters (see Appendix A).

The total uncertainty σ_tot,ν (r)) is defined as $${\sigma }_{\text{tot},v}^2(r))={\sigma }_v^2(r)+{\left({\delta }_v{I}_v(r)\right)}^2,$$(9)

where σ_ν (r) is the standard error of the mean obtained while azimuthally averaging the intensity profiles (see Sect. 2.3). Instead δ_ν is the flux calibration error. Following the ALMA Technical Handbook¹ and the Guide to Observing with the VLA², we set δ_ν = 10% for ALMA Band 9, 7, 6, 1 and VLA Q Band; instead, we adopted δ_ν = 5% for ALMA Band 3 data.

While sampling, we allowed the parameters to vary within uniform priors of ranges: 0 ≤ Tdust /K ≤ 100, 0 ≤ (Σ_dust/gcm⁻²) ≤ 30 and 0 ≤ (a_max/cm) ≤ 100. We sampled each parameter in linear – instead of logarithmic – scale. It should be noted that this might not be the most appropriate choice for the maximum grain size. Grains larger than 10 cm do not emit efficiently in the millimeter fluxes we are analyzing. Therefore, a logarithmic sampling biased toward lower values or a narrower uniform prior are a better representation of our physical expectations. Here, we still chose a larger upper limit for the maximum grain size and a normal sampling to highlight when the parameters are unconstrained in our results.

As an example, Appendix B shows that, if optically thin observations are included in an SED analysis, both a linear and a logarithmic sampling produce consistent results. Instead, if only marginally thick observations are considered, the complexity of the parameter space can cause the samplers to find a local convergence within a more confined prior or with a logarithmic sampling, leading to a solution that is reasonable but not necessarily accurate. To test the stability of the solution of an SED analysis, we suggest to double check the results with both a linear and logarithmic sampling or to explore the parameter space with other sampling techniques, such as a No-U-Turn Sampler (NUTS).

Following Macías et al. (2021), we also imposed a temperature prior based on the expected temperature profile of a passively irradiated flared disk in radiative equilibrium: $$T_{\text{dust}}(r)={\left(\frac{\psi L_{\star }}{8\pi r^2{\sigma }_{\text{SB}}}\right)}^{0.25}.$$(10)

In the previous equation σSB is the Stefan-Boltzmann constant and ψ is the flaring angle that we vary uniformly between 0.01– 0.06 (Huang et al. 2018; Dullemond et al. 2001). We assume a normal distribution for the stellar luminosity L_⋆ centered at 5 L_⊙ and with a conservative standard deviation of 50% L_⋆ .

To assess convergence we followed Goodman & Weare (2010) by introducing in our code criteria based on the integrated autocorrelation time. After a burn-in of 1000 steps, we performed a maximum of 50 sampling, each with a total number of steps equal to 1000. At each sampling cycle, we evaluated the autocorrelation time. We defined a chain converged if: (i) they are longer than 100 times the estimated autocorrelation time (i.e., the walkers are independent); (ii) for each parameter, the autocorrelation time at different cycles changes by less than 1% (i.e., the convergence is stable). We chose to initialize 12 walkers (equal to 4 times the dimension of the problem), a larger number of walkers in this parameter space results in increasing autocorrelation times. We tested different moves of the ensamble and we concluded that the default stretch move of Goodman & Weare (2010) leads to faster convergence.

The result of an emcee analysis are the marginalized probability distributions of the parameters. In the following, we define the results of an SED analysis as the value corresponding to the 50th percentiles of the distribution and, as errors, the difference between the 50th and the 16th and 84th percentiles: $50{\text{th}}_{-(50\text{th}-16\text{th})}^{+(84\text{th}-50\text{th})}$.

4 Benchmark of the SED analysis

4.1 Different band combinations

In the wide range of frequencies spanned by ALMA and VLA, the dust thermal continuum emission of a protoplanetary disk is characterized by very different optical depths that aid the reconstruction of the dust properties from a disk SED differently. In a simplified case where the dust opacity is dominated only by absorption and scattering is negligible (at millimeter /cm wavelengths this assumption is valid only for small particles a_max ≤ 10⁻³ cm and might not be appropriate for protoplan- etary disks; Birnstiel 2011), Equation (4) simplifies to $I_v=B_v\left(T_{\text{dust }}\right)\left(1-{\text{e}}^{\left(-{\tau }_v^a/\mu \right)}\right)$, where ${\tau }_v^a$ is the optical depth derived by the absorption opacity alone.

In this simplified scenario, optically thick emissions help get an accurate constraint on T_dust . If ${\tau }_v^a\gg 1$, the previous equation reduces to $I_v^{\text{thick }}=B_v\left(T_{\text{dust }}\right)$ and we are able to measure the dust temperature free from degeneracies with the other dust properties. Instead, optically thin observations are sensitive also to the dust mass and particle sizes. If ${\tau }_v^a\ll 1,I_v^{\text{thin }}\simeq B_v\left(T_{\text{dust }}\right){\tau }_v^a\simeq B_v\left(T_{\text{dust}}\right){\text{Σ}}_{\text{dust}}{k}_v^a$.

At millimeter wavelengths, the spectral behavior of the absorption coefficient is usually parameterized as k_ν ∝ ν^β . Therefore, in the Rayleigh–Jeans approximation, valid for disks at low optically thin frequencies, the optically thin intensity scales as $I_v^{\text{thin }}\propto v^{2+\beta }$, and the spectral index of the emission between two millimeter wavelengths, α = log (I_ν1 /I_ν2)/ log (ν₁/ν₂) spans from 2 (optically thick) to 2+β (optically thin). Thus, in the thin case, the maximum grain size can be easily inferred from the measure of β as β = α − 2 (Testi et al. 2014; Ricci et al. 2010).

How the nonnegligible dust scattering affects the previous picture is not trivial, since the intensity of an optically thick disk includes a dependence from the albedo if scattering is considered: $I_v^{\text{thick }}~B_v\left(T_{\text{dust }}\right)\cdot \sqrt{1-{\omega }_v}$ (Zhu et al. 2019). Therefore, a complete SED analysis including absorption and scattering effects is important.

Even in a single spectral window the spectral properties of a disk dust emission are not uniform. The most recent observations of protoplanetary disks with ALMA at high angular resolution have revealed that dust substructures are ubiquitous, and the emission in the central parts of the disks, as well as in some of the bright rings, could be very optically thick even at long wavelengths (e.g., Jin et al. 2016; Pinte et al. 2016; Liu et al. 2017). In other words, for certain disk morphologies and dust masses, the 3 mm dust continuum emission could still be marginally thick and insufficient to the measurement of the dust properties from the SED. In these cases, longer-wavelength, but unfortunately, lower-resolution observations are to be included in the analysis (such as the 7 mm emission from ALMA Band 1 or VLA Band Q observations).

Here, we assessed, for the two different disk morphologies we are simulating (compact, with shallow rings; and extended, double-ringed disk models), which ranges of dust masses result in marginally thick emission at 3 mm and 7 mm and how these optical depths affect the measurement of dust properties from the disk’s SED. We note that we are assuming a reasonable distance of 140 pc (the distance of nearby star-forming regions such as Taurus and Ophiuchus), and a typical particle size distribution (millimeter-sized in the inner disk and rings and submillimeter in the outer disk).

To vary the optical depths of our simulated observations, we kept the maximum grain size profiles a_max(r), measured from our DustPy simulations (see Sect. 2.1), fixed while manually rescaling the dust surface densities Σ_dust(r). We reproduced three scenarios: i) the case where the disks have enough dust mass that their emission in the inner disk (<15 au)/brighter ring is marginally thick at both 3 and 7 mm (optical depths τ_3,7mm > 1–2); ii) the more expected case of less massive but still detectable disks where the inner disk/brighter rings are marginally thick at 3 mm but thin at longer wavelengths (τ_3mm > 1–2, τ_7mm < 1); and iii) a last case where the disks emission is thin at both 3 mm and 7 mm even in the inner disk and in the most dense substructures (τ_3,7mm < 1).

For both the compact disk (top figure) and the extended one (bottom figure), the optical depths of the multi-wavelength observations in the three scenarios are shown in the top rows of Figure 2. The right panel show the dust-rich disk case (i), obtained by increasing the dust mass until the inner disk (for the compact disk) or the brighter ring (for the extended disk) have marginally thick 7 mm emission. Instead, the central and left panels, represent, respectively, the average dust mass case (ii) and the dust-poor case (iii).

For both disk morphologies, the dust mass ranges that result in the three scenarios are almost the same. This is expected since the simulated disks have the same age, distance, initial mass and have a similar grain size profile. Since we chose typical values for all these properties, similar ranges are to be expected for the majority of protoplanetary disks. The ranges of dust masses that lead to marginally thick 3 mm and thin 7 mm emissions are: 0.05 < M_dust/M_Jup < 1.07 for the compact disk and 0.03 < M_dust/M_Jup < 1.52 for the extended one. With the current ALMA and VLA sensitivities a disk that falls in the dust-poor picture (with masses smaller than 0.03/0.05 M_Jup) is typically too faint to have detected 7 mm continuum emission at a distance of ∼140 pc. On the other hand, recent dynamical mass measurements of the total mass of protoplanetary disks (Martire et al. 2024), suggest that – for a standard gas-to-dust ratio of 100 – even the brightest disks do not typically exceed a dust mass of 1–2 M_Jup . Consequently, we expect the majority of the detected protoplanetary disks to have dense marginally thick 3 mm emission regions and, instead, optically thin 7 mm emission.

Figure 2 also shows the spectral indices of the multiwavelength intensity profiles computed in the three scenarios (bottom row). To disentangle optical depth effects from the ones due to the limited resolution of the observations, here we computed the radial intensity profiles without applying simobserve. Thus, these spectral indices have infinite resolution. Again, we note that the spectral indices obtained in the intermediate mass case (ii) are in agreement with the few real measurement available. A spectral index between Band 7 and 6 (0.87 mm and 1.29 mm) of ∼2 in the dense disk regions have been commonly observed in T Tauri stars (see, e.g., HL Tau, Carrasco-González et al. 2019; CI Tau, Zagaria et al., in prep.; and TW Hya, Macías et al. 2021). A similar spectral index is expected between Band 7 and Band 9 (0.45 mm) since, at these frequencies, the emission is optically thick. Sierra et al. (2021) shows spectral indices between 100 GHz and 226 GHz (Band 3 to Band 6) between 2.5 and 3.5 within 50 au for the MAPS disks. Spectral indices between Band 4 (2.1 mm) and 8 mm (VLA Q+Ka Band) of 3.3 (at 50 au) and 3.5 (at 80 au) have been reported for HL Tau by Carrasco-González et al. (2019).

To evaluate how the emissions at different optical depths aid the reconstruction of the dust properties, we performed, for each dust mass scenario, three SED analyses including, each time, a different combination of the intensity profiles at various wavelengths. In particular, we analyzed the following Band combinations: (a) from 0.45 mm to 7.46 mm; (b) from 0.87 mm to 7.46 mm; and (c) from 0.87 mm to 3.07 mm. We note that we always included at least three wavelengths in the SED analysis to assure a good retrieval of the dust parameters. Degeneracies in the results of SED analyses have been observed by (Li et al. 2024) when including only two wavelengths.

Case (a) is a well-studied source for which we have very thick emission (Band 9 here but the same apply for ALMA Band 10/8 observations), intermediate emission (ALMA Band 7 and 6), marginally thin emission (ALMA Band 3 or Band 4) but also optically thin emission (ALMA Band 1 or VLA Q/Ka Band). Compared to this previous picture, in case (b) we miss the very thick emission. Case c) is the more common case of a disk with available high-resolution Band 7, 6, and 3 observations. It should be noted that while this band combination is more realistic (since it requires the smallest wavelength range), it is still uncommon to have these high-resolution multi-wavelength observations: there are only 11 disks with available high-resolution (<0.05″) Band 7, 6, and 3 observations.

The results of the three SED analyses for the compact (left) and the extended (right) disks are shown in Figure 3. The three rows represents, in order: the dust-rich case (i), the average dust mass case (ii), and the dust-poor case (iii). Each subplot displays the dust temperature (top), dust surface density (middle) and maximum grain size (bottom) measured from an SED analysis of: (a) all the bands (in blue, 0.45 mm → 7.46 mm); (b) Bands 7, 6, 3 and 1 (in orange, 0.87 mm → 7.46 mm); and (c) without the long wavelengths (in green, 0.87 mm → 3.07 mm). The black line indicates the dust properties of the models.

The conclusions of Figure 3 are clear: in a dust-poor disk in which the 3 mm emission is already optically thin, longer- wavelength observations are not needed to accurately reconstruct the dust properties from the disk SED. Both the broad (blue) and narrow-band (green) analysis in the last row of Figure 3 recover the expected dust temperature, surface density and maximum grain size almost everywhere. Only in the inner disk (<15 au) our simulations have marginally thick 3 mm emission and we can distinguish the results obtained with and without including the 7 mm emission. Instead, for a protoplanetary disk with an average dust mass (second row), the difference between including or not the 7 mm emission is nonnegligible. The dust content in the dense regions of the disk, such as the inner disk and the rings, is robustly probed only including this optically thin emission. Lastly, in a dust-rich scenario (first row), even lower frequencies than 40.0 GHz (7 mm) needs to be included in the SED analysis and none of the three band combinations we are testing is sufficient to accurately measure dust properties in the whole disk. This discussion is independent from the disk morphology for similar optical properties.

Figure 3 also highlights that by including a conservative but reliable prior (see Equation (10)) the dust temperature is always well constrained in an SED analysis independently from the optical depths of the observations. Contrary to the expectations of an absorption-only scenario (see the beginning of this section), we found that completely optically thick observations do not aid in measuring the dust temperature. For a disk with an average dust mass, this suggests that Band 9 observations are not required. At millimeter wavelengths, the main contribution to opacity in B9 is scattering, making it less effective as a temperature tracer despite its high optical depth.

We note that in all three dust mass scenarios the maximum grain size in the first and second gaps of the extended disk morphology are not well constrained. This has a simple explanations: the dust population in these gaps has sizes of ∼10–100 µm, one order of magnitude lower than the submillimeter wavelengths. Grains that are much smaller than the mean free path of the photons (a < λ/2π) act as Rayleigh scatterers. Consequently, the opacities are insensitive to the grain sizes and we are not able to constrain them from an SED analysis (see Figure 4 of Birnstiel et al. 2018).

On top of including flux calibration errors while analyzing the multi-band SED, in Appendix C, we also investigated the potential impact of these errors by performing 50 separate SED analyses. In each analysis, the intensity profiles at the various wavelengths were rescaled within their respective flux calibration error margins. We observe that the spread in dust temperature measurements is less than 17%. Additionally, the dust surface density and maximum grain size are nicely constrained. Except near the very inner disk, the order of magnitude for both of these parameters is consistently matched.

Fig. 2 Optical depths (upper row) and spectral indices (lower row) computed for the compact (top figure) and the extended disk (bottom figure) for the three scenarios: (i) dust-poor disks with optically thin emission at both 3 mm and 7 mm (left); (ii) dust-average disks with optically thick emission at 3 mm and thin at 7 mm (middle); and (iii) dust-rich disks with optically thick emission at both 3 and 7 mm (right). The dust masses that give rise to these optical properties (for fixed maximum grain size and distance) are (from left to right) 0.06 MJup, 0.14 MJup, and 1.10 MJup for the compact disk; and 0.03 MJup, 0.13 MJup, and 1.53 MJup for the extended disk. The vertical gray lines represent the position of the bright (B) and dark (D) rings in the dust surface density profile. The horizontal lines in the spectral indices plots indicate α = 3.7 (typical value in the ISM) and α = 2.0 (optically thick emission in the Rayleigh-Jeans limit). The dark gray shaded region in the optical depth plots indicate τ > 2, the light gray region τ > 1.

4.2 Different spatial resolutions

After testing how the different optical depths of the multiwavelength observations contribute to the measurement of dust properties from a disk SED, now we focus on the effect of the finite resolution of disks observations. As discussed in the previous subsection, we expect the majority of protoplanetary disks to be characterized by marginally thick 3 mm emissions regions. Therefore, if we want to accurately measure the dust content even of the denser substructures, we need to include longer- wavelength observations, such as the 7 mm emission, detected by ALMA in Band 1 or by VLA in Q Band. This comes with a limit: currently, ALMA Band 1 has a maximum resolution of 0.1″, insufficient to resolve substructures in most disks (<0.03″– 0.05″). High resolution is partially reached at high frequencies by VLA (0.05″ for Q Band, 0.06″ for Ka) but with poorer image quality.

In this section, we compare the results of SED analyses of observations with different resolution. In particular, we focus on the following three resolutions: 0.05″, the maximum resolution currently achievable by ALMA in Band 3 and by VLA in Q Band; 0.1″, the maximum resolution achievable in ALMA Band 1; and 0.2″, a more realistic and time-efficient resolution for the current capability of ALMA. For a bright disk with a relatively low elevation from the ALMA site, such as HL Tau, the maximum resolution achievable in ALMA Band 1 is 0.144″. In order to detect the outer disk of HL Tau with a signal-to-noise ratio greater than 5, a total observing time of 5.46 hours is required. By decreasing the resolution to 0.2″, the total observing time is reduced by a factor of ∼4. Less bright disks, such as Elias 2-24, Tw Hya, or PDS 70, can all be detected with a resolution of 0.2″ and an rms of 5 µJy in ALMA Band 1 in less than 6 hours. With this resolution, longer-wavelength observations are in principle achievable for a sample of average brightness disks.

Following the conclusions of the previous subsection, here we analyze, for both the compact and the extended disks, at all three resolutions, the dust thermal emission from 0.87 mm to 7.46 mm of an average dust-massive disk (case ii). To simulate the effects of the observations on the multi-wavelength intensity profiles that we computed from Equation (4), we applied the simobserve task of CASA as explained in Sect. 2.3. The 7.46 mm emission at an angular resolution of 0.05″ was simulated using a VLA antennas configuration, at lower resolution, we emploied an ALMA Band 1 antennas configuration (see Table 1).

The dust temperature, surface density and maximum particle size measured from the SED analyses for the compact (top plot) and extended (bottom plot) disks are reported in Figure 4. Each plot displays the three solutions obtained from observations with resolutions: 0.05″ (in orange); 0.1″ (in green); and 0.2″ (in purple). No solution is shown near the location of the second gap of the extended disk because the simulated observations of this deep substructure present negative intensity profiles. Thus, we excluded this region from the analysis.

Figure 4 clearly highlights the need for high-resolution observations: for both the compact and extended disk, a mediumresolution SED analysis is insufficient to get an accurate measurement of the dust properties. The top plot of Figure 4 shows that the dust properties of the compact disk measured from 0.2″ multi-wavelength observations generally underestimate those from the model in the inner disk and, at the same time, overestimate them in the outer parts. This is expected: at a distance of 140 pc, a resolution of 0.2″ corresponds to 28 au, which means that the radial extent of the compact disk is resolved in only two beams. As a result of beam smearing, the flux of the inner disk is smeared toward the outer disk. Therefore, we have a reduction of the emission in the central parts while the emission of the outer disk is contaminated by the smeared flux. This leads to an underestimated dust temperature in the inner ∼20 au, that is otherwise consistent with expectation. The dust surface density is also only underestimated in the inner 20 au and robustly constrained elsewhere (with limited resolution); however, the maximum grain size has the opposite behavior: it is adequately constrained in the inner parts but severely overestimated outside. This results in underestimated optical depths in the inner disk and overestimated in the outer parts.

The beam smearing effect is reduced by increasing the resolution of the SED analysis. However, it is to be noted that even for the highest resolution (0.05″, orange line), we still see an underestimation of the dust temperature and surface density in the inner 20 au. An even higher resolution would be required to resolve both the inner disk and the substructures in this compact disk with shallow gaps. The gap widths are ∼5 au for the first ring and ∼7 au for the outer one and only marginally resolved even in the higher-resolution case (0.05″ = 7 au at 140 pc).

The conclusions for the beam analysis comparison for the extended disk are different (see bottom plot in Figure 4). This morphology presents a fainter inner disk and brighter dust rings than the compact disk and is almost three times more extended. Here, the beam smearing mainly affects the rings: because of the limited resolution the peak emission of the rings contaminates the flux coming from the gaps. This results, as evident from the purple line (0.2″) in the Figure, in uniform optical depths in the gaps and in the rings. In particular, we are not able to distinguish the dust surface density or the maximum grain size across the rings and gaps. A resolution of 0.1″ is a great improvement to the SED characterization for this disk morphology (see green line in Figure 4). While we still notice a nonnegligible amount of flux contamination and a consequent overestimation of the dust surface density and maximum grain size in the gaps, here we witness different optical depths coming from the dark and bright regions of the disk. As previously discussed, this more ideal scenario is now achievable only for a few selected disks for which we have high-resolution long-wavelength observations. With high-sensitivity (≳50 hours) VLA Band Q observations, we can even push to a resolution of 0.05″ (see orange line in Figure 4), which for a disk with a morphology and extent similar to the one we are simulating, provides accurate and robust measurements of the disk dust content.

It is to be noted that, independently from the resolution, we are always able to estimate the maximum grain size at the peak of the substructures and the maximum level of grain growth in a disk. The same is true for the dust mass (see the labels in Figure 4): even though we miss spatial accuracy in the dust surface density, the total dust mass obtained by integrating the measured density profiles over the analyzed area remains accurate. This is true only if the optically thin observations are included in the analysis (see for a comparison Figure 3). A discussion on the uncertainties of dust mass measurements from optically thick observations is carried out in Sect. 6.1.

Fig. 3 Results of the three SED analyses for the compact (left) and the extended (right) disks. The three rows of panels represent from top to bottom: the dust-rich case, the average dust mass case, and the dust-poor case (see Sect. 4.1). Each subpanel displays the dust temperature (top), dust surface density (middle), and maximum grain size (bottom) measured from an SED analysis that includes noiseless infinite resolution emissions from: 0.45 mm → 7.46 mm (in blue); 0.87 mm → 7.46 mm (in orange); and 0.87 mm → 3.07 mm (in green). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models.

Fig. 4 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from SED analyses of the 0.87 mm → 7.46 mm emissions for the average dust mass compact disks (top plot) and extended disks (bottom plot). Each plot displays three solutions obtained from observations with different resolution: 0.05″ (in orange); 0.1″ (in green); and 0.2″ (in purple). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters; the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

5 New approaches to the SED analysis

In the previous section we commented on the necessity of including optically thin high-resolution observations in probing the dust content of protoplanetary disks. This limits the number of sources to which we can apply this dust characterization technique. While higher-frequencies and higher-resolution (<0.05″) observations are becoming more and more common for an increasingly wide sample of disks (Andrews et al. 2018), high- quality ALMA observations still have a maximum resolution of 0.05″ at 3 mm and 0.1″ at 7 mm. With this resolution, we are unable to resolve the local dust accumulations in disks where the conditions for triggering streaming instabilities and forming planets are likely to occur.

In this section, we test on simulated observations techniques designed to handle a combination of lower- and higher-resolution intensity profiles and increase the spatial resolution of SED analyses.

5.1 The double-resolution analysis technique

Here, we propose a new approach to the dust characterization method tailored for when a fraction of the intensity profiles of our multi-wavelength observations have significantly lower resolutions than the others.

A requirement of SED analyses is to convolve all the observations to the lowest resolution. However, if some of the intensity profiles have medium to low resolution (>0.05″-0.1″), the results of the SED analysis are spatially unresolving substructures (see Figure 4 and the discussion in Sect. 4.2). To take advantage of the emission at all frequencies while resolving dust properties, we can use the unresolved results of this broadband medium-resolution analysis, in particular the maximum grain size, as a prior of a second narrower-band analysis that includes in the SED only the highest resolution profiles.

Following the technical capabilities of ALMA, we tested the realistic case of a disk with high-resolution (0.05″) Band 7, 6, and 3 and medium-resolution (0.2″) Band 1 observations. The double-resolution SED analysis technique for this scenario is sketched in Figure 5: we performed a medium-resolution broadband SED analysis (0.87 mm → 7.46 mm with a resolution of 0.2′′) to measure the maximum grain size in the disk; with this information, we defined a prior on a_max(r) and implemented it in a second analysis of the SED including only the 0.05″ resolution 0.87 mm → 3.07 mm observations. We defined the prior as an asymmetric Gaussian centered at the 50th percentiles of the marginalized probability distribution of a_max and with widths equals to two times the errors of the results (i.e., the differences between the 50th and 16th and the 84th and 50th percentiles). Defining the prior by interpolating the full density probability of the maximum grain size is equally valid but computationally expensive.

This double-resolution approach can be applied even if we have more than one lower-resolution observation, as long as the remaining high-resolution intensity profiles still have a high enough signal-to-noise ratio for the results to not be affected by the decreased amount of data in the SED.

The full potential of this double-resolution technique is shown in Figure 6. Figure 6 shows in purple the dust properties measured from the broadband lower-resolution SED (0.87 mm → 7.46 mm with 0.2″ of resolution) for both the compact (top plot) and the extended disk (bottom plot). As discussed in Sect. 4.2, brightness peaks are smeared because of the limited resolution. This results in almost uniform dust surface density and maximum grain size over the disk extent (for the compact disk) or across gaps and rings (for the extended disk). Conversely, we succeed in radially resolving the dust content of both models if we only analyze high-resolution Band 7, 6, and 3 observations (yellow solution). However, as we concluded in Sect. 4.1, the dense regions of protoplanetary disks (i.e., the inner disks and the rings) have marginally thick emission even at Band 3, thus without optically thin information we are not able to probe their dust content. This changes by incorporating in this high-resolution analysis the optically thin information of the unresolved broadband SED (in purple) through a prior. As depicted in orange in Figure 6, the optically thin information contained in the prior aids in measuring the dust content even in these denser regions without compromising resolution.

Fig. 5 Sketch of the double SED analysis technique designed to handle a combination of lower- and higher-resolution multi-wavelength observations of protoplanetary disks.

Fig. 6 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from SED analyses for the compact (top plot) and the extended (bottom plot) disks. Each plot displays three solutions. In purple we show the results of analyses from 0.87 mm → 7.46 mm with a resolution of 0.2″; in yellow from 0.87 mm → 3.07 mm with a resolution of 0.05″ ; in orange from 0.87 mm → 3.07 mm with a resolution of 0.05″, but including a prior on the maximum grain size that is based on the results of the medium-resolution analyses (in purple; see Sect. 5.1). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the position of the bright (B) and dark (D) gaps in the dust surface density profile.

5.2 Enhancing the resolution with Frank

Alternatively, we can increase the spatial resolution of an SED analysis by individually enhancing the resolution of each radial intensity profile of the multi-wavelength observations. Unlike the double-resolution method, this approach does not necessarily require to reduce the number of profiles used in the SED analysis, allowing, in principle, for a better noise management.

To extract higher angular resolution intensity profiles from the multi-wavelength dust thermal emissions, we performed a nonparametric visibility modeling of the continuum visibilities using Frankenstein (Frank) (Jennings et al. 2020). This tool, designed to recover axisymmetric disk structures at a sub-beam resolution, presents itself as an alternative to the traditional technique that we use to recover the brightness profile of our models from their interferometric observations: deprojecting and azimuthally averaging their CLEAN model images (see Sect. 2.3). During the cleaning procedure, a disk image is convolved with a beam that, as previously discussed and observed, smears and reduces the brightness of all the features of the source of size comparable to or smaller than the beam. In contrast, Frank infers the unconvolved brightness distribution by fitting the Real part of the visibilities as a function of the uv-distance. Because of the hypothesis of axisymmetry, the Imaginary part of the visibilities is not analyzed.

Again, we followed the current capabilities of ALMA by modeling the interferometric data of the 0.05″ Band 7, 6 and 3 and 0.2″ Band 1 observations.

So far, we simulated noiseless observations. This allowed us to focus on the intrinsic contribution of each frequency to the SED analysis without the limitations introduced by the different noise levels of the various bands. Since the performance of a Frank fit is deeply influenced by the noise level, here we employed the add_vis_noise task of Frank to simulate the effect of thermal noise on the visibilities. This function introduce a spread in the visibilities by adding Gaussian noise according to the data weights. After calibration, the weight of a couple of antennas (i, j) is defined as $1/{\sigma }_{i,j}^2$, where σ_i,j is the rms noise of a given visibility³. If the antennas systemic temperature can be approximated to an average temperature for the interferometer, this noise can be related to the expected rms noise of the disk image (σ_rms) as $${\sigma }_{i,j}={\sigma }_{\text{rms}}\sqrt{\frac{N(N-1)}{2}\frac{t_{\text{source}}}{t_{i,j}}{n}_{\text{pol}}}={\sigma }_{\text{rms}}\sqrt{N_{\text{vis}}{n}_{\text{pol}}},$$(11)

where N is the number of antennas, N_vis the total number of visibilities, n_pol the number of polarization (here =1), t_i,j the integration time, and t_source the time on source.

To manually introduce noise in the visibilities, we followed the weights definition of Equation (11). We employed the ALMA Sensitivity Calculator⁴ to evaluate the expected noise level of our observations for a source at HL Tau’s declination after 3 hours of observing time. We adopted the same observing time while producing our noiseless simulations (see Sect. 2.3).

We fitted at each frequency the noisy visibilities in logarithmic brightness space, to reduce oscillatory artifacts and prevents negative values in the reconstructed brightness profile. We set the disk radius out to which the fit is performed (R_max) equal to 1″ for the compact disk and 2″ for the extended disk, as recommended in Jennings et al. (2020), this is greater than 1.5 times the disks outer edge. To ensure a sub-beam fit, we defined the number of collocation points used in the fit (N) equal to 300. To minimize artifacts in the brightness profiles, we set α = 1.30 and w_smooth = 10⁻¹. These two hyper-parameters mimic, respectively, the signal-to-noise ratio (S/N) threshold below which the data are not used when fitting the visibilities, and the smoothness of the fit to the power-spectrum. Jennings et al. (2020) suggest varying α between 1.05 and 1.30 and w_smooth between 10⁻⁴ and 10⁻¹, thus α = 1.30 ensures we are not including noise-dominated features in the fit, and w_smooth = 10⁻¹ that we are not overfitting.

To check whether we reached a sub-beam resolution, we computed the intrinsic angular resolutions of the Frank radial profiles. Following Sierra et al. (2024), we generated synthetic visibility observations of a delta function centered at a radius of r₀ = R_max/2 with the same uv-coverage of each dataset. A nonzero center was chosen to avoid boundary effects. However, Sierra et al. (2024) report that negligible differences are found when using a different r0. We then fitted the input profile with the same method and choice of hyper-parameters. We defined the intrinsic resolution as the broadening (i.e., the full width half maximum, FWHM) of the brightness profile reconstructed by Frank.

The Frank fits of the multi-band observations of the compact (left) and extended (right) disk models with the expected noise level after 3 hours of observing time are shown in green in Figure 7. The light gray shadow in each panel represents the corresponding radially averaged noise level. In the same Figure, the model intensity profiles are plotted in black and the tclean profiles (i.e., the profiles obtained by averaging the clean simulated observations) (see Sect. 2.3) in blue.

Various considerations follow from Figure 7. The sensitivity requirements to reach a satisfying signal-to-noise ratio are extremely different in the various ALMA bands. With three hours of observing time, the outer ring of the compact disk is detected with a radially averaged signal-to-noise ratio of ~330 in Band 7, ~180 in Band 6, ≲5 in Band 3, and undetected in Band 1. It should be noted that even for nondetection Frank still produces a result that for an unresolved compact disk could be confused with a physical one (see bottom left panel of Figure 7). Similar S/Ns are observed for the extended disk morphology.

Following these nondetection, we visually inspected the minimum radial signal-to-noise ratio in the outer ring of both the compact and extended disk required to produce Frank profiles from the Band 1 observations that approximate the model. For both morphologies, signal-to-noise ratios lower than 15 were insufficient. With the current ALMA capabilities, these low noise levels in Band 1 observations requires times on source of respectively 35 h for the compact disk and 27 h for the extended one. Such expensive observations are still not enough for an accurate reconstruction of the 7.46 mm intensity profile (see the orange solutions in the bottom panels of Figure 7). The stricter requirement in the S/N with respect to the Band 3 profiles is likely caused by the lower uv-coverage of the Band 1 observations. We note that we are increasing the signal-to-noise ratio of the observations by increasing the weights, but the number of uv-points are fixed to the one produced from a 3 h observation. This is equivalent to binning the visibilities, an action that is typically adopted in faint observations to increase the S/N (see Sierra et al. 2024).

Artifacts in the Frank profiles are seen even for the higher- resolution Band 3, 6, and 7 observations, independently from the signal-to-noise ratio and despite the conservative values of the α hyper-parameter. The bright Band 7 and 6 observations are characterized by a bump in the intensity profile in the inner disk of both the extended and compact morphologies, suggesting that this is probably the high-resolution response of Frank to the abrupt inner cut-off of the disk model at ~10 au. Less expected bumps are found at the center of the compact disk in the Band 3 emission and in the first ring of the extended disk in Band 7 and 6. Additionally, we observed that different choices of hyperparameters result in different profiles.

In each panel of Figure 7, we are also reporting the intrinsic resolution of the Frank profiles. We tested the robustness of these values against changes of the N and R_max parameter. Since here the disks geometry is well known, the uncertainties on these values are low. The improvement in resolution compared to the tclean profile of the Band 7 and 6 observations is significant: ≳3 for the compact disk and ≳2.5 for the extended one. This is not the case for the fainter Band 3 observations (a factor ~1.4 improvement for the compact disk, and ~1.2 for the extended one). With a high signal-to-noise ratio but limited uv-coverage, the Frank Band 1 profiles have an intrinsic resolution of 0.0782″ (~2.5) and 0.1341″ (~1.5) for the compact and extended disks.

Since the increase in resolution of the Frank fits of the optically thin observations is still insufficient to resolve substructures, we applied the double-resolution technique on the Frank intensity profiles with Band 1 again as the lowest resolution profile (see Sect. 5.1). In this scenario, we are still reducing the amount of data in the SED. The results are shown in green in Figure 8. For a comparison, we also displayed the results of the double-resolution technique applied on the tclean profiles (blue).

For the compact disk, the difference between analyzing the Frank and the tclean profiles is negligible. The improvement in resolution of the higher-frequency profiles (0.037″ instead of 0.05″, the intrinsic resolution of the Band 3 Frank fit) is not significant in resolving the substructures of this small disk. Unexpectedly, the important artifacts introduced in the Band 1 Frank fit are not affecting the dust properties determination. This might suggest that the Band 1 flux – instead of the spatially resolved observation – is sufficient in helping reconstructing the dust properties in the higher-resolution SED analysis. This is also supported by the conclusions of the previous section (see Figure 6): lower-resolution Band 1 observations can still aid the dust content reconstruction when used as a prior of a narrow band higher-resolution analysis.

Similar results are found for the more complex morphology of the extended disk model. Despite the artifacts of the Frank profiles, the Band 1 information included in the prior are aiding the reconstruction of the dust properties in the non-noise dominated rings. In the noise-dominated regions of this disk, the optical depth is not constrained and no information on the maximum grain size or on the dust density can be gathered.

Fig. 7 Frank fits of the multi-band observations of the compact (left) and extended (right) disk models. The fits to the visibilities with the expected noise level after 3 hours of observing time are shown in green. The light gray shadow in each panel is the corresponding radially averaged noise level. The model intensity profiles are plotted in black and the tclean profiles in blue. For Band 1 (bottom panel), the fits to the visibilities with an azimuthally averaged signal-to-noise ratio of 15 in the outer ring are also shown in orange and the corresponding noise level is represented as a dark gray shadow. The vertical gray lines in each panel represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile. In each panel the intrinsic resolution of the Frank fits are reported.

Fig. 8 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from SED analyses for the compact (top plot) and the extended (bottom plot) disks. Each plot displays the results of the double-resolution technique presented in Sect. 5.1 applied on the Frank intensity profiles (in green) and the tclean intensity profiles (in blue). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters; the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the position of the bright (B) and dark (D) gaps in the dust surface density profile.

6 Discussion

6.1 Measuring dust masses from millimeter emissions

More than 1000 protoplanetary disks in several nearby starforming regions have been surveyed with ALMA in the 211 GHz–373 GHz frequency range, where the dust emission and the CO isotopologues are bright (e.g., Pascucci et al. 2016; Ansdell et al. 2016; Barenfeld et al. 2016). These measured millimeter fluxes have been converted into dust masses via (Hildebrand 1983), $$M_{\text{dust}}=\frac{1}{\mu }\frac{F_v{D}^2}{\overline{k_v}{B}_v\left({\overline{T}}_{\text{dust}}\right)},$$(12)

where F_ν is the dust flux at frequency ν, B_ν (${\overline{T}}_{\text{dust }}$) is the blackbody emission at the dust temperature ${\overline{T}}_{\text{dust }}$ and frequency ν, µ is the cosine of the disk inclination which has been assumed equal to 1 for most of the published surveys, D is the distance of the source and ${\overline{k}}_v$ is the dust absorption opacity. The overline indicates that the dust properties are averaged over the disk extent.

Equation (12) can be easily derived from the thermal emission of a 1D razor-thin vertically isothermal slab (Equation (4)) under the assumption of optically thin emission, but spectral indices between Band 7 and 6 of ~2, consistent with optically thick emission, have been commonly observed in T Tauri disks (e.g., Carrasco-González et al. 2019; Macías et al. 2021).

Through self-consistent radiative transfer models at 0.88 mm and 1.33 mm, Liu et al. (2022) observed that this traditional dust mass measurement method can underestimate the mass by a factor of a few to several hundreds, depending on the optical depth along the line of sight set mainly by the true dust mass, disk size, and inclination. Here, we complement Liu et al. (2022) results with an analysis of dust mass estimates obtained via the multiwavelength analysis of simulated disks with different optical depths and within different frequency ranges.

In Figure 9, we show the correction factor (i.e., the ratio of the measured dust mass to the model dust mass) that we obtain from our SED analyses. To investigate this technique ability to accurately measure dust masses for different optical depths, we again refer to the three scenarios we introduced in Sect. 4.1: (i) a dust-rich (M_dust ≳ 1.6M_Jup), (ii) a dust-average (0.06M_Jup ≲ M_dust ≲ 1.6 M_Jup), and iii) a dust-poor (M_dust ≲ 0.06 M_Jup) disk model. The results obtained from SED analyses of the emission from 0.87 mm to 3.07 mm are displayed with a circle in blue for the compact disk model and in pink for the extended disk morphology. With a black-bordered square we are showing the results of SED analyses that also include the 7.46 mm emission. We note that the dust masses have been measured from SED analyses of multi-band observations with infinite resolution. As discussed in Sect. 4.2, the total dust mass is always recovered in SED analyses even if the intensity profiles are unresolved. The gray shadow highlights the range of dust masses that we can typically expect for ALMA observations of T Tauri disks (scenario ii; see Sect. 4.1). The horizontal dashed line marks an accurately estimated dust mass.

For a comparison, we are also showing in Figure 9 the correction factors that we obtain by measuring the dust masses from the 0.87 mm (green cross) and 7.46 mm (orange cross) fluxes through Equation (12). We follow Andrews et al. (2013) in defining a stellar-luminosity based average dust temperature: ${\overline{T}}_{\text{dust }}=25{\left(L_{\star }/L_{\odot }\right)}^{0.25}$ K. In the literature, the dust average opacities are typically computed as the absorption opacity of a grain of 1 mm size (Andrews et al. 2013; Beckwith et al. 1990). For the “DSHARP” composition and assuming a power-law index of the grain size distribution q = 3.5 (see Sect. 2.2), this equals 3.5 cm²/g at 0.87 mm and 0.1 cm²/g for 7.46 mm. We vary the flux within typical flux calibration errors (10% for both 0.87 mm and 7.46 mm) to estimate the uncertainties of these masses.

From Figure 9, it is evident that, independently from the optical depth of the observations and the disk morphology, classical flux-based mass measurements are not accurate. This is expected for the 0.87 mm flux (green cross). The 0.87 mm emission remains marginally optically thick even for a dust-poor disk, concealing a significant amount of dust and leading to a severely underestimated dust mass. On the other hand, we expect the optically thin 7.46 mm flux to provide a better constraint of the dust content. This is not observed in Figure 9. The masses derived from the 7.46 mm flux (orange cross) are even more underestimated than the ones derived from Band 7. This highlights another important simplification of this flux-based mass measurement technique: the assumption of a single arbitrary dust opacity instead of a more consistent value modeled at each radius. By assuming a grain size of 1 mm (close to the model’s maximum grain size of ~3 mm), we are severely overestimating the dust opacity in Band 1, which in turn results in an underestimated dust mass. According to Mie theory of dust opacity, the 7.46 mm flux presents a peak in the dust opacity at a grain size of the order of λ/2π ~ 1 mm (Birnstiel et al. 2018).

Very different dust masses are measured from SED analyses. If optically thin observations are not included in the SEDs, the dust mass is overestimated – instead of being underestimated as in the case of flux-based masses, and more importantly, not robustly constrained. This is the case for the dust masses derived for the dust rich disk models with any combination of wavelengths, as well the 0.87 mm → 3.07 mm analysis for the dust-average scenario. This behavior can be understood by looking at Figure 3, which shows the probability distribution of dust parameters obtained from SED analysis for both compact and extended disks across all three dust mass scenarios. When only optically thick observations are considered, the SED becomes insensitive to dust opacities, leading to different combinations of Σ_dust and a_max having similar probabilities. This implies a broad marginalized distribution for Σ_dust, and thus an unconstrained dust mass. The overestimation of the masses is a consequence of linear sampling and the limited prior assumptions.

On the contrary, optically thin observations in the SED analysis results in dust masses that are both accurate and robustly constrained. This is the case for the masses derived from the 0.87 mm → 7.46 mm (square) SED analyses of the dust-average and dust-poor disk models as well as the 0.87 mm → 7.46 mm (circle) SED analysis of the extended dust-poor disk (pink). For the compact disk (blue) in the dust-poor scenario, the dust mass obtained from the 0.87 mm → 3.07 mm SED analysis is slightly overestimated. As shown in Figure 3, the dense inner disk of this compact morphology remains marginally optically thick even at low dust masses.

To conclude, the results presented in Figure 9 demonstrate that the SED analysis is a significantly more effective method for measuring dust masses in protoplanetary disks, provided that optically thin observations are included. While this analysis is more expensive than simply measuring dust masses from millimeter fluxes, the fact that the spatial resolution of these optically thin observations does not critically affect the derived dust mass estimates (see Sect. 4.2) makes the method more feasible for samples of disks, as high-resolution observations are not strictly required. At the same time, it is also evident that an analysis that relies only on optically thin fluxes at a single wavelength is insufficient (see orange crosses in Figure 9) since this method relies on assumptions about an average dust temperature and, in particular, an average dust opacity.

Fig. 9 Correction factor, i.e., the ratio of the measured to the model dust mass, for various measurement techniques. The circles refer to the masses measured from SED analyses including the emission from 0.87 mm to 3.07 mm, in blue for the compact disk model and in pink for the extended disk morphology (see Sect. 2). With a black-bordered square we show the masses measured from SED analyses that also include the 7.46 mm emission. The correction factors that we obtain by measuring the dust mass from the 0.87 mm (green cross) and 7.46 mm (orange cross) fluxes through Equation (12) are also displayed. The gray shadow highlights the range of dust masses that we can typically expect from ALMA detections of T Tauri disks (scenario ii; see Sect. 4.1). The horizontal dashed line marks a accurately estimated dust mass.

6.2 Extending SED analyses with the ALMA WSU and the ngVLA

In most of the published SED analyses of protoplanetary disks, the longest wavelengths incorporated are high- to mediumresolution (0.05″–0.1″) marginally thick observations at 2–3.6 mm (ALMA Band 3/4). Alongside the more frequently observed Band 6/7, fully optically thick data from ALMA Bands 9/10 are occasionally included.

The benchmarking of the dust characterization method on simulated observations, as shown in Sect. 4, indicates that fully optically thick observations do not significantly enhance the measurement of dust properties if an accurate prior can be included. This is in contrast with what was observed by (Kim et al. 2019) while performing a synthetic ALMA multi-band analysis to find the best band configuration for constraining the dust content of TW Hya. The authors found that the dust properties of the disk were better constrained when including ALMA Band 10/9 observations. However, it is to be noted that their radiative transfer model does not account for scattering. As discussed in Sect. 4.1, at millimeter wavelengths the main contribution to opacity in B9 is scattering. Consequently, in the presence of scattering, the emission in Band 9 could be not as good as a temperature tracer as suggested by its high optical depth. Nevertheless, we note that including Band 9/10 observations in a dust characterization analysis might aid noise management by providing an additional point in the SED (see Sect. 4.1).

More importantly, our analysis in Subsections 4.1 and 4.2 leads to a clear conclusion: to accurately constrain the dust content of protoplanetary disks through SED analyses, high- resolution optically thin observations must be included. For the typical T Tauri disks we detect with ALMA, (0.07 M_Jup ≲ M_dust ≲ 1.4 M_Jup), and under the assumption of a “DSHARP” dust composition (see Sect. 2.2), this requires wavelengths longer than ~3 mm and resolutions ≲0.05″. The wavelength at which the emission transitions from optically thin to optically thick can be approximately identified from the spectral indices. As discussed at the beginning of Sect. 4.1, the spectral index of the emission at millimeter wavelengths spans from 2 for optically thick emission, to 2 + β for optically thin emission, with β being the power index of the absorption coefficient. This picture is complicated in a nontrivial way by dust scattering. However, spectral indices above a ~3 threshold, which are sensitive to substructures, are expected for completely optically thin observations (see Figure 2).

As previously discussed, the lack of longer wavelength observations of protoplanetary disks in the SED analyses can be traced back to the fact that dust thermal emission at these wavelengths can only be observed with high resolution by the VLA. The limited sensitivity and image quality of this facility constrains the analysis to bright disks visible at VLA declinations. SED analyses incorporating high-resolution VLA data are currently available for only two protoplanetary disks: HD 163296 (Guidi et al. 2022) and HL Tau (Carrasco-González et al. 2019).

The new capabilities of ALMA Band 1 to produce high- quality and time-efficient long-wavelength observations of average-brightness sources open a new window in dust characterization of protoplanetary disks. Although the maximum resolution of ALMA Band 1 (0.1″) remains insufficient for resolving substructures, the double-resolution analysis we are proposing in Sect. 5.1 effectively integrates optically thin information from these medium-resolution observations into a higher-resolution framework, allowing for high-resolution results.

With the upcoming Wide Sensitivity Upgrade (WSU), ALMA will further enhance our ability of measuring dust properties even of the faintest protoplanetary disks. The initially double, and eventually quadruple bandwidth upgrade will boost the continuum sensitivity of ALMA’s receivers by a factor of 3 up to 6 (Carpenter et al. 2020). This is crucial for easy detection of faint optically thin emission. Moreover, the WSU comes with a new access to long wavelengths: ALMA Band 2. Band 2 is designed to have the widest frequency coverage of any interferometric (sub)millimeter facility, from 67 GHz to 116 GHz, encompassing the Band 3 frequency range and extending to a maximum wavelength of 4.48 mm. This new Band is key for characterizing dust of faint disks, allowing for marginally thick (~3 mm) observations even of the faintest sources. While it does not extend as far as Band 1 in terms of wavelengths, a maximum wavelength of 4.48 mm can still provide sufficiently optically thin observations for faint sources with the advantage of increased sensitivity and resolution.

An ongoing study led by the European Southern Observatory (ESO) also addresses the feasibility of extending ALMA’s maximum baselines to as long as 32 km (Carpenter et al. 2020). The corresponding increase in resolution at higher frequencies (e.g., ALMA Bands 6 and 7) will complement the milliarcsecond resolution that the next generation VLA (ngVLA) will provide at longer wavelengths (~3–9 mm) (Wilner et al. 2024). This synergy between the two facilities will be critical, as it will allow us to resolve and characterize even the more compact and fainter protoplanetary disks, which constitute the majority of the population, including systems similar to the one that formed our solar system.

The increase in resolution provided by these facilities will also allow us to probe the dust content of even the most inner regions of protoplanetary disks, where rocky planets supposedly form. However, we know that at wavelengths longer than 6 mm, emission from ionized gas or other non-dust mechanisms can contaminate the dust emission originating from the inner parts of the disk (Carrasco-González et al. 2019). To remove this contribution, high-resolution observations over a wide range of long wavelengths are needed. The VLA and the next-generation Very Large Array (ngVLA) will be fundamental in providing these observations in the Northern Hemisphere, in the same way the Square Kilometre Array (SKA) will allow us to access the Southern Hemisphere.

With this boost in resolution across a wide wavelength range, we will finally have the capability to probe the fine structures and dust distributions in these smaller, less luminous disks, enabling a deeper understanding of the early stages of planet formation in a broader range of environments.

7 Summary and conclusions

Spatially resolved multi-wavelength studies of the dust thermal emission of protoplanetary disks are essential for understanding the initial stages of planet formation. By accurately studying both the particle size distribution and the dust surface density in these disks, we can identify where and when the conditions necessary to trigger streaming instabilities are met.

In this work we benchmarked the state-of-the-art dust characterization technique, based on SED analyses at different radial locations, on simulated multi-band data. We included completely optically thick observations (ALMA Band 9, 0.45 mm), commonly observed 0.87 mm and 1.29 mm data (ALMA Band 7 and 6); marginally thick 3 mm emissions (ALMA Band 3); and a longer-wavelength 7 mm intensity profile. We designed two disks models: a compact disk with shallow gaps (similar to TW Hya) and an extended double-ring morphology (a proxy for HD 163296).

We first tested which bands aid the reconstruction of the dust properties from the SED, varying the optical depths of the observations by manually rescaling the dust mass. We find that wavelengths longer that 3 mm are needed to accurately reconstruct the dust properties in dense disk regions such as the inner disks and the rings. Surprisingly, fully optically thick observations do not significantly enhance the measurement of dust properties, but they might aid with noise by providing an additional point in the SED. We verified this for disks simulated at a distance 140 pc, assuming a DSHARP dust composition, and for dust masses between ~0.06 M_Jup and 1.6 M_Jup. Higher dust masses are rare even for the brightest disks; masses lower than ~0.1 MJup are not detected with the current sensitivities of the more advanced radio facilities such as ALMA and VLA.

As a second test, we compared the dust properties derived from SED analysis with different resolutions: 0.05″, the maximum resolution achievable by ALMA in Band 3 and by VLA in Q Band (~7–8 mm); 0.1″ , the maximum resolution of ALMA Band 1 in Cycle 11; and 0.2″, a more realistic and time-efficient resolution for the current capabilities of ALMA. We observe in order to resolve substructures and characterize the radial variation of dust properties in the rings and gaps, that high-resolution (≲0.05″) observations are necessary. For coarser estimates of the total dust mass and maximum grain sizes at ring locations, lower- resolution (0.2″) SED analyses are sufficient. Thus, an analysis of multi-band fluxes can still be relevant to study the mass budget and the level of grain growth of samples of protoplanetary disks, but only the advent of milliarcsecond resolution facilities such as the ngVLA will allow us to actually resolve substructures and probe the planets’ birth environments, even in the more compact and fainter disks, which constitute the majority of the population.

Because of the limited sensitivity and image quality achievable at long wavelengths (7–8 mm) by the VLA, we propose a new approach to dust characterization that effectively integrates optically thin information from medium-resolution (0.2″) ALMA Band 1 observations into a higher-resolution framework. We demonstrate that the information obtained from a medium-resolution SED analysis that includes long wavelength observations can be used as a prior of a second narrower-band SED analysis of only higher-frequency higher-resolution intensity profiles. This technique opens a new window in the dust characterization of faint disks that will be further expanded with the upcoming ALMA’s Wide Sensitivity Upgrade. This upgrade will boost the sensitivity of ALMA’s receivers. Additionally, ALMA Band 2 will extend up to 4.85 mm, which can still provide sufficiently optically thin observations for faint sources with the advantage of increased sensitivity and resolution.

We also tested how the flux reconstruction tool Frank can help increase the resolution of an SED analysis by individually enhancing the resolution of each radial intensity. We note that, even for high signal-to-noise ratios, Frank introduces artifacts in the fitted profiles. Thus, these analyses should be carefully treated.

Finally, we discussed the total dust mass that we derive from the SED analyses and compared it with the traditional method of deriving dust masses from (sub)millimeter fluxes. Accurate dust mass measurements can be retrieved from the SED only if optically thin tracers are included. If not, the dust masses are not robustly constrained and are generally overestimated. On the contrary, dust masses measured from the 0.87 mm fluxes can be underestimated by one order of magnitude since the high optical depths at these wavelengths can hide substantially large amount of dust. We also observed that the accuracy of the results is independent from the spatial resolution of the SED within the explored range (0.05″–0.2″). We observed that mass derived from the 7.46 mm fluxes are also not accurate. A singlewavelength flux-based mass measurement method still relies on assumptions about average dust temperatures and opacities. Due to the varying morphologies, ages, and levels of grain growth in protoplanetary disks, dust opacities can differ significantly. A more complex and tailored approach to averaging opacities is needed, along with a more accurate method for estimating uncertainties.

Acknowledgements

We thank the anonymous referee for the helpful report. T.C.Y. is supported by JSPS KAKENHI Grant Number JP23KJ1008. P.C. acknowledges support by the ANID BASAL project FB210003.

Appendix A Fitting q: The power index of the grain size distribution

As discussed in Sect. 2.2, a commonly used power-law approximation for the grain size distribution in protoplanetary disks is n(a) da ∝ a^−q da.

Mathis et al. (1977) observed a power law index of the grain size distribution (q) equal to 3.5 in the interstellar medium. In protoplanetary disks, however, grain growth and accretion processes can modify this value, both globally and, more significantly, locally. Higher values of q are often found in disk gaps, which are populated by small dust particles, while in the rings q values closer to 3.5 are observed, consistent with fragmentation processes (e.g., Macías et al. 2021).

In this work, we assume a uniform value of q = 3.5, as this approximation does not affect the conclusions of our benchmarking analysis. Nonetheless, because of the strong correlation between a_max and q, here we assess the ability of the SED analysis method to retrieve both parameters. We do this by comparing the results of two analyses: one with q fixed at 3.5, and another with q as a free parameter.

Fig. A.1 Dust temperature (first panel), dust surface density (second panel), maximum grain size (third panel), and power index of the grain size distribution (fourth panel) measured from an SED analysis of the compact disk model (see Sect. 4.1). The SED includes emission from 0.87 mm to 7.46 mm, simulated with infinite resolution. The blue solution is obtained considering the power index of the grain size distribution (q) as one of the parameter of the SED analysis. Instead, the orange solution is the result of an SED analysis fixing q = 3.5, the expected value in our simulations. The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

Fig. A.2 Same as Figure A.1, but for the extended disk model (see Sect. 4.1).

Figure A.1 and A.2 show the dust temperature, dust surface density, maximum grain size and power law index of the grain size distribution (from top to bottom) measured from an SED analysis when fixing (orange solution), or fitting (blue solution) q. To not introduce limited resolution or optical depths effects in this comparison, we included in the SED intensity profiles from 0.87 mm to 7.46 mm with infinite resolution (see Sect. 4.1). Figure A.1 shows the comparison for the compact disk model, Figure A.2 for the extended disk model (see Sect. 2). The dust temperature, dust surface density and maximum grain size measured in the two SED analyses are consistent with each other. Moreover, the bottom panel of Figure A.1 and A.2 suggests that this dust characterization method efficiently deals with the correlation between a_max and q. The only regions where we are not able to reproduce the expected value of q are the marginally thick inner disk of the compact disk and the deep gaps of the extended disk. In these deep gaps, grains sizes are much smaller than the mean free path of the photons (a < λ/2π). In this Rayleigh regime, the opacities are insensitive to the grain sizes and we are not able to constrain them from an SED analysis (see Figure 4 of Birnstiel et al. 2018).

Appendix B Linear versus logarithmic sampling

As discussed in Sect. 3, in this work we estimate the dust properties of our disk models by comparing, through a Bayesian approach, their multi-band intensity profiles to a multiwavelength 1D slab model (equation 4). In particular, we employ the code emcee (Foreman-Mackey et al. 2019) to estimate the posterior distribution of the parameters of the model: the dust temperature, dust surface density, and maximum grain size.

Fig. B.1 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from an SED analysis of the dust-average (top plots) and dust-rich (bottom plots) compact (left) and extended (right) disk models (see Sect. 4.1). The SED includes emission from 0.87 mm to 7.46 mm, simulated with infinite resolution. The orange solution is obtained by linearly sampling the parameter space, the blue solution with a logarithmic sampling. The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

Emcee samples the N-dimensional parameter space, where N is the number of parameters, to find the N-tuple that best reproduces the data, while ensuring the exploration remains within the boundaries set by a prior. In this work we choose to define uniform priors and sample the parameter space in linear scale.

For the maximum grain size, this might not be the most appropriate choice. Our choice of letting a_max vary between 0 cm and 100 cm does not actually reflect the grain size distribution we expect to detect at millimeter wavelengths. The peak of intensity at millimeter wavelengths is reached for grain sizes of the order of ~λ/2π, thus grains ≳ 10 cm emit at lower frequencies than the ones included in the SED. Furthermore, given the wide range of grain sizes expected in protoplanetary disks, spanning various orders of magnitude, a uniform prior in the logarithmic space is a better representation of our physical expectations.

Figure B.1 shows the dust temperature, dust surface density, and maximum grain size (from top to bottom) measured from SED analyses of the dust-rich (top plot) and dust-average (bottom plot) compact (left) and extended (right) disk models (see Sect. 4.1). In each panel, two solutions are shown: in orange, the results obtained if we sample a_max in linear space within 0 and 100 cm; in blue, the solution obtained with a uniform prior in logarithmic space. The SED analyses include emission from 0.87 mm to 7.46 mm, simulated with infinite resolution.

From the top plots of Figure B.1 it is clear that, if both optically thick and thin observations are included in the SED analyses, the dust properties of a disk are well reproduced independently from the sampling method. The two sampling methods yield different results only within the deep gaps of the extended disk model. These regions are optically thin across all wavelengths, and under the Mie theory of opacity, their dust opacity can be matched by two distinct dust populations (refer to Figure 4 of Birnstiel et al. 2018). Linear sampling tends to favor larger grain sizes, whereas logarithmic sampling identifies smaller grain sizes as a solution. While this serves as an interesting exercise, it is important to note that such regions, characterized by optically thin emission even at high frequencies (Bands 7 and 6), are not commonly targeted by ALMA and are challenging to detect.

On the other hand, when the observations included in the SED analysis are not optically thin, we see inconsistency in the dust properties measured with the two sampling method (see bottom plots of Figure B.1). The 7.46 mm emission of the dust-rich (M_dust ≳ 1.4 M_Jup) disk models are marginally thick in the inner disk and in the rings. In these dense regions, the dust properties measured by sampling the parameter space in linear (orange) of logarithmic (blue) scale differs. While neither of the SED analyses produce robustly constrained results, in a linear sampling the optical thickness of the observations is clearly highlighted by the walkers converging toward the maximum of the prior. Instead, a uniform prior in logarithmic space favors lower a_max values, resulting in more reasonable solutions that are closer to the expected values but not necessarily accurate: in the inner disk of the compact disk (bottom left), the measured maximum grain sizes can be ~2 order of magnitude larger than the true values.

To highlight difficult convergences and optical thickness, in our work we choose to linearly sample the parameter space. To test the stability of the solution of SED analyses of real observations, we suggest to double check the results with both a linear and logarithmic sampling or to explore the parameter space with other sampling techniques, such as a No-U-Turn Sampler (NUTS).

Appendix C Flux calibration errors

Something that we need to consider while performing dust characterization analyses is that, i addition to the uncertainties on the dust properties derived from an emcee analysis of a disk’s SED, uncertainties also arise from flux calibration errors that affect each multi-band intensity profile.

In the SED analyses of our simulated disk models, performed using the Ensemble Sampler emcee (Foreman-Mackey et al. 2019), we have already incorporated flux calibration errors. We define the error on the multi-wavelength intensity profiles, and thus the ’allowed’ distance between the model intensity and the data in the likelihood function, as the quadrature sum of two components: the standard error of the mean, obtained from azimuthally averaging the intensity profiles, and the flux calibration errors. We followed the guidelines from the ALMA Technical Handbook (ALMA Partnership 2023) and the Guide to Observing with the VLA (NRAO 2023), by setting flux calibration errors of 10% at 0.87 mm, 1.29 mm, and 7.46 mm, and 5% at 3 mm.

However, in our simulated multi-band disk observations, the fluxes are set to their expected values. Now, we directly estimate the spread introduced by flux calibration errors in the dust properties derived from the SEDs.

To do so, we perform 50 separate SED analyses. In each iteration, the intensity profiles at the different wavelengths are rescaled within their respective flux calibration error margins. The dust temperatures, dust surface densities and maximum grain sizes resulting from these 50 SED analyses including emission from 0.87 mm to 7.46 mm for the compact disk model are shown in Figure C.1.

By varying the multi-band fluxes within their calibration errors, we observe that the spread in dust temperature measurements is less than 17%. Interestingly, the dust surface density and maximum grain size are also surprisingly well constrained. Except near the very inner disk, the order of magnitude for both of these parameters is consistently matched.

Fig. C.1 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from a SED analyses of compact disk model (see Sect. 4.1). The SED includes emission from 0.87 mm to 7.46 mm, simulated with infinite resolution. In each panel, we are showing 50 different solutions obtained by rescaling each of the multi-band intensity profiles within typical flux calibration errors of radio facilities such as ALMA and VLA (10% at 0.87 mm, 1.29 mm, and 7.46 mm; 5% at 3.07 mm). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

References

All Tables

All Figures

Fig. 1 DustPy outputs at 1 Myr for the compact (left panel) and the extended disk models (right panel). Top panel: grain size distribution of the dust particles as a function of the radial distance from the star. At each radius, the maximum grain size, is highlighted in lime green. Bottom panel: radial profiles of gas (in blue) and dust (in orange) total column densities. The vertical gray lines represent the position of the bright (B) and dark (D) gaps in the dust surface density profile.

In the text

Fig. 2 Optical depths (upper row) and spectral indices (lower row) computed for the compact (top figure) and the extended disk (bottom figure) for the three scenarios: (i) dust-poor disks with optically thin emission at both 3 mm and 7 mm (left); (ii) dust-average disks with optically thick emission at 3 mm and thin at 7 mm (middle); and (iii) dust-rich disks with optically thick emission at both 3 and 7 mm (right). The dust masses that give rise to these optical properties (for fixed maximum grain size and distance) are (from left to right) 0.06 MJup, 0.14 MJup, and 1.10 MJup for the compact disk; and 0.03 MJup, 0.13 MJup, and 1.53 MJup for the extended disk. The vertical gray lines represent the position of the bright (B) and dark (D) rings in the dust surface density profile. The horizontal lines in the spectral indices plots indicate α = 3.7 (typical value in the ISM) and α = 2.0 (optically thick emission in the Rayleigh-Jeans limit). The dark gray shaded region in the optical depth plots indicate τ > 2, the light gray region τ > 1.

In the text

Fig. 3 Results of the three SED analyses for the compact (left) and the extended (right) disks. The three rows of panels represent from top to bottom: the dust-rich case, the average dust mass case, and the dust-poor case (see Sect. 4.1). Each subpanel displays the dust temperature (top), dust surface density (middle), and maximum grain size (bottom) measured from an SED analysis that includes noiseless infinite resolution emissions from: 0.45 mm → 7.46 mm (in blue); 0.87 mm → 7.46 mm (in orange); and 0.87 mm → 3.07 mm (in green). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models.

In the text

Fig. 4 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from SED analyses of the 0.87 mm → 7.46 mm emissions for the average dust mass compact disks (top plot) and extended disks (bottom plot). Each plot displays three solutions obtained from observations with different resolution: 0.05″ (in orange); 0.1″ (in green); and 0.2″ (in purple). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters; the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

In the text

	Fig. 5 Sketch of the double SED analysis technique designed to handle a combination of lower- and higher-resolution multi-wavelength observations of protoplanetary disks.
In the text

Fig. 6 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from SED analyses for the compact (top plot) and the extended (bottom plot) disks. Each plot displays three solutions. In purple we show the results of analyses from 0.87 mm → 7.46 mm with a resolution of 0.2″; in yellow from 0.87 mm → 3.07 mm with a resolution of 0.05″ ; in orange from 0.87 mm → 3.07 mm with a resolution of 0.05″, but including a prior on the maximum grain size that is based on the results of the medium-resolution analyses (in purple; see Sect. 5.1). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the position of the bright (B) and dark (D) gaps in the dust surface density profile.

In the text

Fig. 7 Frank fits of the multi-band observations of the compact (left) and extended (right) disk models. The fits to the visibilities with the expected noise level after 3 hours of observing time are shown in green. The light gray shadow in each panel is the corresponding radially averaged noise level. The model intensity profiles are plotted in black and the tclean profiles in blue. For Band 1 (bottom panel), the fits to the visibilities with an azimuthally averaged signal-to-noise ratio of 15 in the outer ring are also shown in orange and the corresponding noise level is represented as a dark gray shadow. The vertical gray lines in each panel represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile. In each panel the intrinsic resolution of the Frank fits are reported.

In the text

Fig. 8 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from SED analyses for the compact (top plot) and the extended (bottom plot) disks. Each plot displays the results of the double-resolution technique presented in Sect. 5.1 applied on the Frank intensity profiles (in green) and the tclean intensity profiles (in blue). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters; the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the position of the bright (B) and dark (D) gaps in the dust surface density profile.

In the text

Fig. 9 Correction factor, i.e., the ratio of the measured to the model dust mass, for various measurement techniques. The circles refer to the masses measured from SED analyses including the emission from 0.87 mm to 3.07 mm, in blue for the compact disk model and in pink for the extended disk morphology (see Sect. 2). With a black-bordered square we show the masses measured from SED analyses that also include the 7.46 mm emission. The correction factors that we obtain by measuring the dust mass from the 0.87 mm (green cross) and 7.46 mm (orange cross) fluxes through Equation (12) are also displayed. The gray shadow highlights the range of dust masses that we can typically expect from ALMA detections of T Tauri disks (scenario ii; see Sect. 4.1). The horizontal dashed line marks a accurately estimated dust mass.

In the text

Fig. A.1 Dust temperature (first panel), dust surface density (second panel), maximum grain size (third panel), and power index of the grain size distribution (fourth panel) measured from an SED analysis of the compact disk model (see Sect. 4.1). The SED includes emission from 0.87 mm to 7.46 mm, simulated with infinite resolution. The blue solution is obtained considering the power index of the grain size distribution (q) as one of the parameter of the SED analysis. Instead, the orange solution is the result of an SED analysis fixing q = 3.5, the expected value in our simulations. The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

In the text

	Fig. A.2 Same as Figure A.1, but for the extended disk model (see Sect. 4.1).
In the text

Fig. B.1 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from an SED analysis of the dust-average (top plots) and dust-rich (bottom plots) compact (left) and extended (right) disk models (see Sect. 4.1). The SED includes emission from 0.87 mm to 7.46 mm, simulated with infinite resolution. The orange solution is obtained by linearly sampling the parameter space, the blue solution with a logarithmic sampling. The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

In the text

Fig. C.1 Dust temperature (top panel), dust surface density (middle panel), and maximum grain size (bottom panel) measured from a SED analyses of compact disk model (see Sect. 4.1). The SED includes emission from 0.87 mm to 7.46 mm, simulated with infinite resolution. In each panel, we are showing 50 different solutions obtained by rescaling each of the multi-band intensity profiles within typical flux calibration errors of radio facilities such as ALMA and VLA (10% at 0.87 mm, 1.29 mm, and 7.46 mm; 5% at 3.07 mm). The continuum line represents the 50th percentiles of the marginalized probability distribution of the parameters, the shaded area extends from the 16th to the 84th percentiles. The black line indicates the dust properties of the models. The vertical gray lines represent the positions of the bright (B) and dark (D) gaps in the dust surface density profile.

In the text