© The Authors 2026

1 Introduction

In recent years, the advent of extreme adaptive optics (AO) has made it possible to observe disks at high resolution in visible and near infrared (NIR) light (Benisty et al. 2023). In this wavelength range, we see the stellar light that scatters off the small submicron to micron (sub-µm to µm) sized grains, whose scattering opacity dominates (Hashimoto et al. 2011; Garufi et al. 2014) over the large grains. These small grains are fundamental in shaping the disk evolution, as they control how stellar radiation is absorbed and re-emitted, shaping the thermal structure of the disk (Chiang & Goldreich 1997; Dullemond et al. 2001; D’Alessio et al. 1999) and its chemistry (Öberg et al. 2023). Their abundance and size distribution therefore offer key insights into dust growth, fragmentation, and the early stages of planet formation (Drążkowska et al. 2023).

Large surveys have been conducted and are still ongoing, such as SEEDS (Tamura 2016), DARTTS-S (Avenhaus et al. 2018), Gemini-LIGHTS (Rich et al. 2022), and SPHERE-DESTINYS (Ginski et al. 2021), meaning that a large sample of observations is becoming available (more than 100 disks observed (Garufi et al. 2026). In these observations, disks exhibit a flared surface. The height of the surface can be measured with geometrical methods (Avenhaus et al. 2018; de Boer et al. 2016; Ginski et al. 2016) or with detailed radiative transfer modeling (Villenave et al. 2019). When rings are present and sufficiently symmetric their projected shapes on the sky can be fitted with ellipses to geometrically estimate the scattering surface height.

While these observations are becoming routine, understanding the physics that sets the scattering surface becomes important for interpreting the structure and dust content of disks. A natural approach is to use numerical radiative transfer models, as done in the literature (Pohl et al. 2017; Villenave et al. 2019; Rich et al. 2022; Muro-Arena et al. 2018; Franceschi et al. 2023). However, given the broad parameter space, full numerical radiative transfer models are computationally expensive and feasible only for individual disks. On the other hand, semi-analytic models would allow us to capture the key physical processes while remaining computationally efficient.

Protoplanetary disks dust masses are usually estimated through the mm-emission, in the optically thin approximation (e.g., Andrews et al. 2018). While this approach is useful to constrain the bulk dust mass, it is insensitive to the population of small µm-sized grains that dominate the disk surface layers and as a result, the mass in small grains is poorly constrained, despite the important role they play in setting the disk thermal structure, in controlling the reaction rates and tracing the early stages of grain growth. In this work, we link the scattering height with the physical properties of the disk, in particular the small dust mass, with a simple semi-analytic model. Starting from radiative transfer considerations, we develop and validate a model that successfully reproduces the emission surface of mid-inclination protoplanetary disks observed in scattered light. We then invert the model and apply it on scattering height measurements from the literature, deriving the mass in small dust for a sample of 10 disks. Section 2 develops the theoretical model, Sect. 3 validates our model against the numerical radiative transfer code MCFOST. Section 4 describes the inversion and mass estimates. Section 5 discusses the implications for disk chemistry and grain evolution, while Sect. 6 summarizes our findings.

2 Methods

In this section, we investigate the basic physics governing the location of the scattering surface, namely, the layer from which most of the emission in scattered light comes. The general expectation that the scattering surface lies close to the τ∗ ≈ 1 surface with respect to stellar photons, where τ∗ is the optical depth from the star to the scattering point (e.g., Dullemond et al. 2001) is well established. However, we aim to identify the key disk parameters that influence the location of the surface and to determine relations between these parameters and the surface height. In Sect. 2.1 we build a semi-analytical model under extremely simplifying assumptions, and in Sect. 2.2 we test the possible approximations and refine the model.

2.1 Semi-analytical model

We considered a thermally stratified disk in hydrostatic equilibrium with inclination i with respect to the observer’s line of sight. The vertical density profile is given by (e.g., Rosotti et al. 2020; Martire et al. 2024) $$\rho (z)={\rho }_0\frac{c_{s,\text{mid-plane}}^2}{c_s^2(z)}\exp (-{\int }_0^z\frac{{\mathrm{\Omega }}_k^2(z^{\prime })z^{\prime }d{z}^{\prime }}{c_s^2(z)}),$$(1)

where ρ₀ is the midplane density, Ω_k = GM_*/(r² + z²)^3/2 is the Keplerian frequency, and the sound speed cs(z) is $$c_s(z)=\sqrt{\frac{k_{B}T(z)}{\mu m_H}}.$$(2)

The disk is a passive irradiated disk, for which the temperature structure is set by stellar irradiation rather than viscous heating. In this regime, which is appropriate for Class II protoplanetary disks at radii ≥1 AU (D’Alessio et al. 1998), the upper layers intercept more stellar radiation than the midplane, resulting in T increasing with height and a cold midplane (Chiang & Goldreich 1997). In the radial direction, T decreases with distance from the star. This vertical and radial temperature structure is directly observed through optically thick CO line emission (Pinte et al. 2018; Law et al. 2021, 2022, 2023; Galloway-Sprietsma et al. 2025). Disks typically show a smooth radial dependence of CO temperature, and following Galloway-Sprietsma et al. (2025) and Law et al. (2021), the temperature is parametrized as $$T(z)=\{\begin{array}{l}{T}_{atm}+(T_{mid}-T_{atm}){\cos }^2(\frac{\pi }{2}\frac{z}{z_q})\text{for\hspace{0.17em}}z<z_q\\ T_{atm}\text{for\hspace{0.17em}}z<z_q,\end{array}$$(3)

where $$T_{atm}(r)=T_{atm,0}{\left(\frac{r}{r_0}\right)}^{q_{atm}},$$(4) $$T_{mid}(r)=T_{mid,0}{\left(\frac{r}{r_0}\right)}^{q_{mid}},$$(5)

and $$z_q(r)=z_0(\frac{r}{r_0}{)}^{\beta }.$$(6)

Here, we are dealing with typical temperatures on the order of 30 to 60 K in the atmosphere, and 10–30 K in the midplane. While our model is general, for the examples of this paper we adopt as a reference the vertical temperature profile of the LkCa15 disk measured by Galloway-Sprietsma et al. (2025). This choice serves as a representative example and all results hold for disks with different thermal structures. For the surface density profile we adopt the Lynden-Bell & Pringle (1974) profile, $$\mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0\left(\frac{r_c}{r}\right)\exp (-\frac{r}{r_c}),$$(7)

where Σ₀ is the normalization and r_c is the cutoff radius. The consequences of this assumption are discussed in Appendix E, where we show that assuming a different slope of the surface density has only a marginal effect on the results. We chose a range of total disk masses 0.001 M_⊙ ≤ M ≤ 0.1 M_⊙.

Scattered light observations are conducted in J (1.2 µm), H (1.6 µm) and K (2.2 µm) bands: at these wavelengths the opacity is dominated by grains with size up to ≈1 µm. To confirm this we consider a grain size distribution n(a) ∝ a^−q with q = 3.5 and a_max = 1 mm; and we compute the grain size a₉₀ that contributes for the 90% of the opacity by solving $$\frac{{\int }_{a_{\text{min}}}^{a_{\text{90}}}\kappa (a)a^{3}n(a)da}{{\int }_{a_{\text{min}}}^{a_{\text{max}}}\kappa (a)a^{3}n(a)da}=0.90,$$(8)

where n(a) is the grain size distribution, κ(a) is the opacity as a function of grain size, and the factor a³ weights each size bin by its mass contribution (since the opacity κ is defined per unit mass of dust). We compute κ(a) with Optool (Dominik et al. 2021), using the composition from Ossenkopf & Henning (1994), and find a₉₀ ≃ 1 μm at a wavelength λ = 1.2 μm. Henceforth, we refer to grains with sizes up to a_max = 1 μm and a grain size distribution exponent of q = 3.5 as small dust. We take the opacity per gram of small dust to be κ = 10⁴ cm²g⁻¹ at λ = 1.2 μm and we assume the dust to gas ratio, f_{gas to dust ≤1 µm}, to be 10⁻³–10⁻⁵. Since we are considering micron-sized grains, we neglected dust settling; this approximation is tested in Sect. 2.2.

To model the stellar radiation propagation through the disk, we adopted the approach of D’Alessio et al. (1999), with the single isotropic scattering approximation (we test this approximation in Sect. 3). In this framework, the intensity along a given ray is given integrating the radiative transfer equation, $$\frac{d{I}_{\nu }}{ds}={\sigma }_{\nu }{B}_{\nu }(T_{\ast })W(r)e^{-{\tau }_{\nu 1}}{e}^{-{\tau }_{\nu 2}(s)},$$(9)

where s is the coordinate along the ray, σ_ν [cm⁻¹] is the scattering coefficient, Bν is the Planck function evaluated at the temperature of the star, and T_*, W(r) = (R_*/r)², is the geometric dilution factor of stellar radiation in the point-source approximation. Then, τ_ν 1 is the optical depth from the star to the scattering point and τ_ν 2(s) is the optical depth from the scattering point to the observer. As shown in Fig. 1, we can divide the path into three contributions to the differential intensity:

• the path from the star to the scattering point (r₀, z₀), characterized by the optical depth ${\tau }_1={\int }_0^{s_0}{\kappa }_{\nu }\rho (s)ds$ (τ_* in Dullemond et al. (2001);

• the scattering at (r₀, z₀), characterized by the cross-section, σ_ν(r₀, z₀) = κ_vρ(r₀, z₀);

• the path from the scattering point to the observer at infinity, characterized by the optical depth ${\tau }_2={\int }_{s_0}^{\mathrm{\infty }}{\kappa }_{\nu }\rho (s)ds$.

We are not concerned with the absolute magnitude of the intensity and, instead, we want to identify whether dI_ν/ds as a function of height in the disk has a maximum, which identifies the scattering surface. Specifically, computing dIν/ds numerically as a function of the height in the disk z, we find a clear maximum (see Fig. 2) and we identify the height z0 of the scattering surface above the midplane as a function of radius.

Next, we address what sets the scattering surface. At a fixed radius, it is possible to analyze how different contributions to the optical path vary with height, as illustrated in Fig. 3, which shows how these three contributions vary with height in the disk at a fixed radius of 150 AU. The left panel displays the optical depths, τ1 and τ2 , the central panel shows the corresponding attenuation factors, e^−τ₁ and e^−τ₂, and the right panel presents the scattering cross-section, σν . When plotting the scattering height, z₀, it becomes evident that this height coincides with the location where τ₁ ≈ 1 (see left panel). To understand this result, we highlight that there is a competition in setting the scattering surface between the cross-section, σ_ν (right panel) and the transmission factor, e^−τ₁ (middle panel). At heights lower than the scattering surface, τ₁ > 1 and the intensity of the stellar radiation is highly attenuated by the exponential decay, implying that there is little radiation to scatter towards the observer. At heights higher than the scattering surface, the transmission factor flattens out, and the reduced cross-section implies that although these locations are highly irradiated by the star, most of this radiation simply passes through, with the same net result that little radiation reaches the observer. Consequently, the scattering height is where τ1 ≈ 1, maximizing the cross-section, while maintaining sufficient transmission. We note that in general τ₂ ≪ τ₁ and at the scattering surface τ₂ ≪ 1; therefore, its effect on the total intensity from the scattering surface is negligible.

To further confirm this result, we computed the scattering height varying the disk mass (ranging from 0.001 M_⊙ to 0.1 M_⊙) and the inclination. We find that this result holds for all the disk masses and for inclinations, i ≤ 60°. For higher inclinations, τ₂ becomes closer to one and obscures the scattering surface.

The scattering height can therefore be computed solving the implicit equation, τ₁ = 1 : $$\begin{gathered} {\tau }_1={\kappa }_{\nu }\int \rho (r,z)ds={\kappa }_{\nu }{\int }_0^{s_0}\frac{\mathrm{\Sigma }(r)}{\sqrt{2\pi }H(r)}\cdot \\ \frac{c_{s,\text{mid-plane}}^2}{c_s^2(z)}\exp (-{\int }_0^z\frac{{\mathrm{\Omega }}_k^2(z^{\prime })z^{\prime }d{z}^{\prime }}{c_s^2(z)})ds=1, \end{gathered}$$(10)

where s is the coordinate along the path of the ray, such that r = s cos θ and z = s sin θ, with θ = arctan(z₀/r₀). The midplane scale height H(r) is computed taking into account the midplane temperature via $$H(r)=\sqrt{\frac{k_B{T}_{\text{mid}}(r)}{\mu m_{H}G{M}_{\ast }}{r}^3}.$$(11)

The scattering height is therefore determined solving the following equation for s₀, namely, $$\begin{gathered} {\mathrm{\Sigma }}_0{r}_c\sqrt{\frac{\mu m_{H}G{M}_{\ast }}{k_B}}{\int }_0^{s_0}\frac{\sqrt{T_{\text{mid}}(r)}}{T(r,z)r^{5/2}}\exp (-\frac{r}{r_c})\cdot \\ \cdot \exp [-\frac{\mu m_{H}G{M}_{\ast }}{k_b}{\int }_0^z\frac{z^{\prime }d{z}^{\prime }}{(r^2+z^{\prime 2}{)}^{3/2}T(r,z^{\prime })}]ds=1. \end{gathered}$$(12)

Finally, it is important to note that the determination of the scattering surface is a nonlocal problem. The height of the surface at a given radius depends not only on the local column density but also on the integrated density along the entire optical path from the star to the scattering point. This highlights the role of the global disk structure in shaping the observed scattered light distribution.

Fig. 1 Model of the path of a ray through the disk.

Fig. 2 Differential intensity as a function of z/H for r = 150 AU.

2.2 Testing approximations

In Sect. 2.1, we introduce several approximations that need to be tested. In this section, we evaluate the impact of dust settling. In Sect. 3, we test the single scattering approximation. A test of the isotropic scattering is given in Appendix A and it is found to be irrelevant for the analysis conducted here. Finally, we show that the isothermal approximation is not appropriate in Appendix B.

Here, we focus on testing whether accounting for dust settling would affect the scattering height. The vertical distribution of the dust in protoplanetary disks is set by the balance between gravity, which pulls the grains towards the mid-plane, and turbulent diffusion, that stirs them up. Small grains that are more coupled to the gas are diffused more efficiently, while bigger grains, which are less coupled to the gas, are settled towards the mid-plane (Birnstiel 2024; Villenave et al. 2020). We treated settling following Dullemond & Dominik (2004), who state that the gas to dust ratio g(m, z) for a particles of mass, m, as a function of height can be computed integrating the following equation, namely, $${\rho }_{\text{gas}}\frac{\mathrm{\partial }g(m,z)}{\mathrm{\partial }t}-D(m,z)\frac{\mathrm{\partial }g(m,z)}{\mathrm{\partial }z}+g(m,z)v_{\text{settle}}(m,z)=0,$$(13)

where, in the stationary case, this can be reduced to $$g(m,z)=g_0(m,z)\exp \left[{\int }_0^z\frac{v_{\text{settle}}(m,z^{\prime })}{D(m,z^{\prime })}d{z}^{\prime }\right],$$(14)

with the settling velocity (Birnstiel 2024) via $$v_{\text{settle}}(m,z)=-\frac{3}{4}\sqrt{{\mathrm{\Omega }}_k}\frac{z}{\rho <v_{th}>}\frac{m}{\sigma }.$$(15)

D(m, z) is the diffusion coefficient (Youdin & Lithwick 2007) via $$D(m,z)=\frac{\alpha <v_{th}>H}{1+{\text{St}}^2},$$(16)

where < v_th > is the mean thermal velocity, α is the (Shakura & Sunyaev 1973) parameter describing the magnitude of turbulence, and St is the Stokes number, which in the Epstein regime (Epstein 1924) is $$\text{St}=\frac{3}{4}\frac{m}{\sigma }\frac{{\mathrm{\Omega }}_k}{<v_{th}>}\frac{1}{\rho }.$$(17)

All the quantities involved depend on the radius, so that the dust to gas ratio, g(m, z), is a function of radius. We computed the scattering surface taking settling into account (i.e. weighting the dust densities with g(m, z)) in the isothermal and nonisothermal case, for grains with a = 1 μm, using m = 4π/3 ρ_sa³ and ρ_s = 1.67 g/cm³. The results are shown in Fig. 4, where the scattering height at r = 100 AU as a function of α (which changes the settling efficiency). We chose to test also the vertically isothermal model because Dullemond & Dominik (2004) found that in the isothermal case settling lowers the scattering surface and we find the same for very low values of the viscosity such as α = 10⁻⁶, even μm grains are settled. In the nonisothermal case; however, we find that the scattering surface is independent of α and therefore settling does not influence the scattering surface. This is because the height at which settling becomes important occurs higher up than the scattering surface. To understand this, we look to Eq. (14). The change in the dust to gas ratio g(z) is due to the ratio between the settling velocity and the diffusion coefficient; since v_settle ∝ ρ⁻¹ and D ∝ ρ, the whole factor is proportional to ρ⁻². The density profile in the nonisothermal case is more vertically extended (a consequence of the higher temperature in the upper layers), as shown in the top panel of Fig. B.1, so that settling is less important. We tested the effect of settling for the other disks with an observationally measured temperature structure for the sample of Galloway-Sprietsma et al. (2025) and we conclude that 1 µm grains are not settled in scattered light observations, as is typically assumed in observational studies (Villenave et al. 2019). Therefore, we find that ignoring settling when computing the scattering surface is a reasonable approximation.

Fig. 3 Contributions to the differential intensity as a function of z in the disk at fixed radius r = 150 AU. Left : optical depth τ1 and τ2 as a function of z/H. Middle : attenuation factors e−τ1 and e−τ2 as a function of z/H. Right : scattering cross-section, σν, as a function of z/H.

Fig. 4 Scattering height at r = 100 AU for different α parameters in the isothermal and non-isothermal case, taking dust settling into account.

3 Comparison with the radiative transfer code MCFOST

To assess the validity of the simplifying assumptions underlying our semi-analytic model, in particular, single scattering, we compared our results against the radiative transfer code MCFOST (Pinte et al. 2006, 2009). This comparison serves as an independent validation of our approach, allowing us to identify the regime of inclinations where our model remains reliable. For simplicity, we compared the numerical radiative transfer results with an isothermal model, as the important physics is the same. We set up a disk model with a tapered profile, as in Eq. (7), and a scale height of $$H(r)=0.1{\left(\frac{r}{100}\right)}^{0.25}.$$(18)

We used astronomical silicates, with the opacities given Fig. 3 in Draine & Lee (1984), where that of the J band (λ ≈ 1.2 μm) does not diverge from the one adopted in our semi-analytic model (see Appendix C). The grain size distribution as the same as in Mathis (1997) n(a) ∝ a^−3.5 and with a maximum grain size of 1 mm. We ignore dust settling and consider isotropic scattering. We ray-traced the image at λ = 1.2 μm.

First we ran the model and used the command -optical_depth_to_cell to extract the optical depth from the star to every point in the disk. We find that the τ₁ ≈ 1 surface in our model is consistent with the one computed with MCFOST. To obtain the scattering height from the computed images, we developed a simple algorithm similar to what is seen in observations (see, e.g., Ginski et al. 2016; de Boer et al. 2021; Avenhaus et al. 2018), but also applicable to smooth disks that do not show rings. In this case, we expect the isophotes to have the same role as ring structures; in fact, since we are considering isotropic scattering, the brightness of a point depends only on its distance from the star, and therefore each isophote traces a single radius. We fit ellipses to the isophotes and in this way we can measure the offset u from the ellipses’ centre to the star and then their height, $$z=\frac{u}{\sin i},$$(19)

with i being the disk’s inclination.

We find that the scattering surface computed with MCFOST and ellipses fitting is comparable with the one of our model and with MCFOST’s τ₁ ≈ 1 surface, as shown in Fig. 5, serving as an empirical validation of our method.

Furthermore, we tested the single scattering hypothesis, which is a critical one and one that makes the problem treatable semi-analytically. To include multiple scattering, we used the radiative transfer numerical model MCFOST turning on and off multiple scattering. We computed the scattering surface fitting the isophotes with ellipses as before. We find that for low inclinations, i ≤ 60°, multiple scattering does not affect the scattering height. For higher inclinations, however, the scattering surface is shifted to higher z. This behavior can be understood in terms of the optical depth structure of the disk. The location of the scattering surface is defined by τ₁ = 1, which marks where stellar photons undergo their first interaction with the disk surface. Since regions with τ₁ > 1 are optically thick to the incoming stellar radiation, scattering events occurring deeper in the disk do not alter the position of the surface. When multiple scattering is accounted for, however, the optical depth to a second scattering, τ₂, becomes relevant. Once the condition τ₂ ≪ 1 no longer holds, additional scattering events occur, producing a net upward shift of the effective scattering surface. Consequently, the inclination at which multiple scattering becomes nonnegli-gible coincides with the breakdown of the τ₂ ≪ 1 assumption, which in our models occurs at approximately i ≃ 60°, Fig. 6 shows an example: we can see that as the inclination increases τ₂ becomes higher (bottom panel), and consequently the difference between the scattering height computed with the single or multiple scattering approximations increases as well (top panel). We point out that our analysis does not account for the case in which observations are performed in polarized light. A more detailed study should consider radiative transfer for scattered polarized light. However, we expect that the dominant contribution to the polarized signal comes from single scattering, as multiple scattering strongly reduces the polarization fraction. This is consistent with our approach, since our model is valid precisely in the single scattering regime, and in fact, observations in polarized light naturally select this regime, as the multiple scattering contribution is suppressed.

Fig. 5 Comparison between the scattering surface from our semi-analytic model, the one computed with MCFOST, and the τ1 = 1 surface from MCFOST, at an inclination of i = 40°.

4 Results

Section 2 presents the direct problem: we determined the scattering height based on the relevant disk parameters. When dealing with observations, however, we are typically more interested in the inverse problem, namely to infer the disk parameters from the observed scattering surface. This is the problem we consider in this section. In particular, from the observed scattering height, we aim retrieve the mass in small (µm sized) dust. Inverting Eq. (12), we obtain $$\begin{array}{r}\kappa M_{sd}(R)=(2\pi {)}^{3/2}{r}_c(1-e^{-R/r_c})\sqrt{k_b}(\mu m_{H}G{M}_{\ast }{)}^{-1/2}\cdot \\ \cdot {\left({\int }_0^{s_0}\frac{\sqrt{T_{mid}(r)}}{T(z,r)r^{5/2}}{e}^{-r/r_c}{e}^{-\frac{\mu m_{H}G{M}_{\ast }}{k_b}{\int }_0^z\frac{z^{\prime }d{z}^{\prime }}{(r^2+z^{\prime 2}{)}^{3/2}T(z^{\prime },r)}}ds\right)}^{-1}.\end{array}$$(20)

Here, the cumulative mass is contained in small dust up to the radius, R; we note that there is a degeneracy between the small dust mass, M_sd, and the opacity, κ. Some considerations on opacity models and uncertainties are in Appendix C. Additionally, we point out that this estimate assume a radially smooth distribution of small dust grains, without accounting for gaps or rings.

Fig. 6 Top: difference in height between the scattering surface computed in the single scattering approximation or accounting for multiple scattering, at 200 AU. Bottom: optical depth from the scattering point to the observer as a function of disk’s inclination.

4.1 Sample

We now wish to build an observational sample of disks we can apply our model to. To apply Eq. (20), we need the following observables:

• coordinates of the scattering surface (r₀, z₀), inferred from the morphology of the rings in scattered light;

• thermal structure of the disk, which is derived from spatially resolved observations of CO rotational lines;

• stellar mass;

• characteristic radius, r_c, that marks the onset of the exponential taper in the surface density profile. In our sources this parameter is constrained by fitting rotation curves.

Therefore, constraints on the properties of small dust grains in protoplanetary disks require multi-wavelength observational data that include scattered light imaging and molecular line emission. From the literature, we were able to find ten disks that have all the required observations. The main limitation consists in having the first two. Table 1 reports the references for the literature measurements of scattering heights, temperature structures and stellar masses, for our disk sample. The cut-off radius r_c has been taken from Martire et al. (2024) and Longarini et al. (2025), however, this parameter does not have a big impact on our results. The top panel of Fig. 7 reports the measured scattering surface from the literature for our sample. The typical uncertainties in the scattering height measurements are on the order of 10–30%, with some outliers (Byrne et al. 2026); we discuss the temperature data in Sect. 5.4.

4.2 Mass in small grains

We derived the small dust mass for the sample described in Sect. 4.1. We computed the small dust masses applying Eq. (20) with κ = 10⁴ cm²/g (see Sect. 2.1 and Appendix C). Since the expression for the mass involves integrals that must be evaluated numerically, a standard error propagation is not straightforward. To estimate the uncertainty on the derived mass, we randomly perturbed the input parameters within their respective observational uncertainties and recompute the resulting M_sd values for each realization. The uncertainty on the mass is then taken as the standard deviation of the resulting distribution in logarithmic space, namely, log(M_sd). Itis important to point out that the scattering surface data points correspond to ring locations; therefore, with our model we estimate the cumulative mass for grains with a ≤ 1 µm up to the ring where the scattering height has been measured. When multiple rings are present, we treat each ring independently, that is, we computed the small dust mass enclosed within each ring without subtracting the contribution from inner rings. This is particularly evident in the case of HD 97048, where the scattering surface is not a monotonic function of radius and consequently neither is the inferred cumulative mass.

The results are presented in Fig. 7: in the top panel, there are the measured scattering surfaces from the literature and in the middle panel, our measurements are given as a function of radius. The small dust masses are also reported in Table 2.

From these results, the following considerations can be made: most of the masses are in the range 10⁻² M_⊕ ≤ M ≤ 1 M_⊕ and the disks with more than one ring in general exhibit an increasing trend of the mass with the radius (with the previously mentioned exception of HD 97048); this is expected because the measured mass is the cumulative mass up to the ring radius. SY Cha has a particularly large uncertainty in its small dust mass estimate. This is primarily due to the uncertainty on the location of the scattering surface, which dominates the overall error. Finally, even if HD 97048 has a high and flared scattering surface, its mass is comparable to the other disks’ mass, and the geometry of the surface can be explained by his unusually high temperature (Pezzotta et al., in prep.).

Fig. 7 Top: measured scattering surfaces from the literature. Middle: small dust masses as a function of radius inferred from the inverse problem in our work. Bottom: small dust mass fraction for all the disks in the sample.

4.3 Small dust ratio

We now try to give additional constraints on the dust distribution combining our results with data from continuum submillimeter (submm) observations. Large grains dominate the mass budget of the grain size distribution; therefore, the dust mass estimated from the continuum sub-mm data can be considered as the total dust mass, M_d. We computed the small dust ratio as $${ϵ}_{sd}=\frac{M_{sd}}{M_d}.$$(21)

It is important to highlight a key distinction: when measuring the small dust mass, we consider only the material within the radius of a given ring, while the large dust mass measurements are integrated over the entire disk. The ideal scenario would be to have very outer rings to capture the total small dust mass, ϵ_sd, or alternatively, to have radius-dependent large dust mass measurements. The latter is only available for the MAPS sample (IM Lup, GM Aur, HD 163296, MWC 480), that has multiwavelength data, and for which Sierra et al. (2021) were able to fit simultaneously for the surface density Σ(r) and the maximum grain size a_max(r). For these disks, we integrated the surface density profile obtaining the mass as a function of radius, $$M(r)={\int }_0^r\mathrm{\Sigma }(r)2\pi rdr.$$(22)

For the rest of the sample, we computed ϵ_sd based on the available total dust masses from Curone et al. (2025); Walsh et al. (2016). These masses are computed under the optically thin assumption, using the equation, $$M_d=\frac{d^2{F}_{\nu }}{{\kappa }_{\nu }{B}_{\nu }}$$(23)

with κ_ν = 3.5 cm² g⁻¹·870 μm/λ (Beckwith et al. 1990). In cases where the scattered light measurements are limited to small-radius rings, the results may be underestimated and those cases are reported in Fig. 7 with pointing up arrows on the error bar. The small dust ratio (i.e. the mass fraction of micron-sized grains) computed for all the disks in the sample is reported in Table 3 and shown in the bottom panel of Fig. 7, where it can be seen that the majority of the ratios are in the range −3.5 ≤ log(ϵ_sd) ≤ −2.5, and all the disks fall in the range −4 ≤ log(ϵ _sd) ≤ −2. We present a discussion of the information that can be extracted from these results in the following subsections.

5 Discussion

5.1 Comparison with exoALMA XV

Rosotti et al. (2025) recently developed a model for the CO emitting height, within the exoALMA sample (Teague et al. 2025), which is very similar to what we present here. They adopted a semi-analytical approach to link the CO emitting height to the disk temperature structure and gas column density, assuming a fixed CO abundance. In the model of Rosotti et al. (2025), the important quantity is the optical depth between the observer and the CO emitting height. This requires integrating vertically along the z direction and makes the problem local, allowing them to infer the local gas surface density. In contrast, reflecting the different physical mechanisms underlying CO emission and scattered light, in our model the relevant optical depth, τ1 , is integrated from the star to the scattering point (i.e. radially). As a consequence, it depends on the full radial and vertical structure of the disk. The problem is therefore nonlocal and our method does not allow us to derive a local surface density, but an integrated quantity corresponding to an enclosed mass instead.

5.2 Impact on chemistry

The abundance of small grains plays an important role in the chemical modeling of protoplanetary disks because it controls both the available surface area for grain-surface reactions and the penetration of UV radiation into the deeper layers of the disk. Since small grains dominate the total surface area, they enhance the adsorption of gas-phase molecules, the formation of ices and surface chemistry reactions, such as hydrogenation (Aikawa & Nomura 2006). At the same time, they dominate the opacity at at shorter wavelengths, so that a higher fraction of small grains increases UV shielding, modifying the photo-dissociation rates and the gas thermal balance (Woitke et al. 2016). In general, a smaller amount of small grains in the disk atmosphere can alter the chemical abundances through the disk and lead to a chemically more active disk interior (Vasyunin et al. 2011). The values we derived for the small dust mass ratio are valuable constraints for chemical models (e.g., DALI; Bruderer et al. 2012; Bruderer 2013) that typically take the small dust mass ratio as input.

5.3 Grain size distribution

We go on to interpret the small dust ratios we have obtained in light of typical grain size distributions used in the protoplanetary disk literature. In the interstellar medium, in debris disks or asteroids, dust follows an MRN distribution via n(a) ∝ a^−q, with q = 3.5 (Mathis 1997; Tanaka et al. 1996), whereas dust in protoplan-etary disks evolves under both coagulation and fragmentation, leading to more complex distributions. Some typical limit cases include fragmentation-limited distributions, with an exponent of 3.5 (Williams & Wetherill 1994), drift-limited ones, with q = 2.5 (Birnstiel et al. 2012), or even steeper distributions, produced when fragmentation is induced by radial drift rather than turbulence, with q = 4 (Birnstiel et al. 2015). Understanding these distributions is essential for interpreting observations, modeling dust evolution, and constraining planet formation processes.

To this point, we computed the small dust mass ratio using reference values of the opacity (i.e., a value of 10⁴ cm²/g for the small dust and a value of 3.5 cm²/g for the large dust). However, in reality, we would expect the two opacities to depend on the properties of the grain size distribution (as well as on the dust composition and porosity). In practice, this means that the value we have found of the small dust ratio is not necessarily consistent with the opacity value we have assumed. To solve this problem, we attempted to use both the submm and scattered light constraint at the same time in a self-consistent manner. As we discuss in the following shortly, solutions consistent with the observational constraints are possible (for a given choice of dust composition and porosity) only for specific grain size distributions. Our results can therefore be interpreted as constraints on the grain size distribution. We detail our procedure to find a self-consistent solution below.

Given an assumed power-law grain size distribution n(a) = c a^−q, we defined the mass in small dust as the mass in grains with a_min ≤ a ≤ a_thresh and the one in big grains as the mass in grains with a_thresh ≤ a ≤ a_max. The small dust ratio therefore becomes $${ϵ}_{sd}=\frac{a_{\text{thresh}}^{4-q}-a_{\text{min}}^{4-q}}{a_{\text{max}}^{4-q}-a_{\text{thresh}}^{4-q}}.$$(24)

We chose a_thresh = 1 µm to be the grain size that accounts for the 90% of the opacity in the J band, as explained in Sect. 2.1, and amin = 0.1 µm (noting that a change in this last parameter only has a minor effect on the results). We compute ϵsd for 10⁻⁴ cm ≤ a_max ≤ 10³ cm and 2 ≤ q ≤ 4, defining this as the true small dust ratio, ${ϵ}_{\text{sd}}^{\text{true}}(a_{\text{max}},q)$.

As previously explained, to compare ${ϵ}_{\text{sd}}^{\text{true}}$ with the value derived from our measurements, we must account for the opacity assumptions underlying both our model and literature submm measurements. Since the opacity depends on the grain size distribution, we have to recompute the opacity used in our measurements for each grain size distribution. Thus, rather than using reference values for the opacities, we computed a measured small dust ratio, ${ϵ}_{\text{sd}}^{\text{measured}}(a_{\text{max}},q)$, for each grain size distribution, consistent with the opacities used in each case.

For each disk in our sample, given a fixed dust composition, we recomputed the dust mass as follows. Using Optool, we calculated the opacity for a_max = 1 µm, while varying q between 2 and 4. We then derived M_sd via Eq. (20). For the large dust mass, we started from the literature values (Martire et al. 2024; Curone et al. 2025; Walsh et al. 2016), computed under the optically thin assumption, with Eq. (23). We then recomputed the opacity, κ′_ν, with Optool for each pair (a_max, q) and obtained a grid of large dust masses as $$M_d=M\frac{{\kappa }_{\nu }}{{\kappa }_{\nu }^{\prime }}.$$(25)

This allowed us to derive ${ϵ}_{\text{sd}}^{\text{measured}}$(a_max, q) self-consistently with the choice of the opacity.

Finally, we compared ${ϵ}_{\text{sd}}^{\text{true}}$ and ${ϵ}_{\text{sd}}^{\text{measured}}$. Only models where the two values are the same are acceptable and, thus, this allows us to evaluate $$f=\frac{{ϵ}_{\text{sd}}^{\text{measured}}}{{ϵ}_{\text{sd}}^{\text{true}}}$$(26)

and determine the regions in the (a_max, q) plane where f = 1, thereby constraining the grain size distribution.

Figure 8 shows the resulting constraints for LkCa15 under various opacity models (Draine & Lee 1984; Pollack et al. 1994; Ossenkopf & Henning 1994; Min et al. 2016; Kataoka et al. 2015). These models differ in terms of both the composition and porosity: the models by Draine & Lee (1984) and Pollack et al. (1994) represent compact grains made of astronomical silicates and carbon, while the model by Ossenkopf & Henning (1994) accounts also for thin icy mantles. The model from Min et al. (2016) represents grains with moderate (30%) porosity, while the model by Kataoka et al. (2015) represents highly porous grains with p = 90%. The corresponding plots for other disks are presented in Appendix D.

Small maximum grain sizes (a_max < 100 μm) are not compatible with typical exponents. Our results then imply that some level of grain growth is necessary to simultaneously interpret submm and scattered light observations. This is consistent with multiwavelength results at radio and mm wavelengths, but we remark that in this case the constraint is provided by a different wavelength and methodology. For some disks, values of q ≈ 3.5 are difficult to reproduce across all opacity models, implying a shallower grain size distribution, with the exception of the Kataoka et al. (2015) model, which assumes highly porous grains. While composition has little effect on the f = 1 contour, porosity plays an important role.

Figure 9 illustrates the impact of porosity on the f = 1 contour for the Pollack et al. (1994) composition. The variation in porosity, significantly alters the contour shape, introducing a major source of uncertainty. A similar effect is seen in Fig. 10, where we compare a grain population composed of compact small grains and highly porous large aggregates with the corresponding zero-porosity case. In both setups, porosity significantly impacts the inferred constraints. However, porosity can be constrained observationally through the dust scattering phase function, which encodes the dust grain properties, especially the porosity and the fractal dimension of the aggregates (Tazaki et al. 2019). The phase function, derived from the intensity of scattered light as a function of angle, can be inferred from the measured scattering height and brightness distribution. Such analyses have been performed for individual disks (Ginski et al. 2016; Tazaki et al. 2023; Columba et al. 2026) and for larger samples (Ginski et al. 2023). A similar analysis is beyond the scope of this work, but we point out that the uncertainty due to the porosity can be reduced in this way.

Finally, we highlight a few caveats of this analysis: we considered a single opacity for either the big or the small grain populations, without considering variations in the vertical or radial opacity. Additionally, these constraints refer to the disk-integrated small-dust ratio. They do not capture radial variations, and substructures within the disk may host grain-size distributions that differ substantially from the global average. Finally, we assume that the small and large grain populations follow the same underlying size distribution, which allows us to directly compare their properties and draw conclusions on the global grain size distribution. While in reality the two populations may evolve independently, relaxing this assumption would make the comparison between the two populations ill-defined, preventing any constraint on the overall grain size distribution.

Fig. 8 Constraints on the maximum grain size and power-law exponent for LkCa15 obtained using different opacity models.

Fig. 9 Constraints on the maximum grain size and power-law exponent for LkCa15 using the Pollack et al. (1994) composition and varying porosity.

Fig. 10 Constraints on the maximum grain size and power-law exponent for LkCa15 using Draine & Lee (1984), Pollack et al. (1994), and Ossenkopf & Henning (1994) opacities. The dotted lines represent the same composition, but considering a porous aggregates scenario, as big grains have a high porosity (p ~ 0.8).

5.4 The importance of the vertical temperature gradient

It is common when modeling proto-planetary disks to assume that they are vertically isothermal. Our results, detailed in Appendix B and Sect. 2.2 show that, at least as far as the scattered light height is concerned, such an assumption has a significant impact on the results. On the one hand, including the vertical temperature gradient pushes the scattering surface significantly higher up than in an isothermal model; on the other hand, it also makes the effect of dust settling unimportant. As a consequence, in this paper we were only able to apply our method to a rather limited disk sample, employing the strong constraint that each individual disk needs to have a measured temperature profile from molecular submm emission, using, for example, the method from Law et al. (2021). There is indeed a considerable difference from source to source in the parameters describing how the temperature varies with the vertical coordinate and one cannot simply extrapolate from the midplane temperature. If more observations of this kind become available in the future, it will enable us to extend our analysis to a larger sample.

It should also be kept in mind that even though we have taken the temperature profiles measured by CO ALMA observations as input for the purposes of this work, there is considerable uncertainty associated to them. Fundamentally, this comes from the fact that the emission of a given transition is associated with a single emitting layer and, therefore, at each radial distance, it only provides a temperature measurement at a single height. The more transitions are available-following the approach of Pezzotta et al. 2025, who employed seven different transitions to extract multiple rotation curves in the disk around HD 163296, and extending it to also measuring the temperature-the better one can allow us to measure how the temperature varies with the vertical coordinate. However, given that most existing measurements are based only on two transitions (most often ¹²CO and ¹³CO), there is a significant amount of extrapolation in the profiles reported in the literature. In addition, it is likely that the temperature structure of real disks is more complex than the simple radial power-law behavior often assumed, as recently shown by Fehr & Andrews (2025) through a nonparametric method and as can be inferred directly using tomography on edge-on disks (Dutrey et al. (2017); Flores et al. (2021) as well as the upcoming results from the DiskStrat ALMA large program).

5.5 Comparison with literature models

Radiative transfer models of protoplanetary disks in scattered light are already available in the literature. For instance, Villenave et al. (2019) used MCFOST to model SPHERE observations of J1852 and J1608, analyzing the different radial extents of small and large dust grains. Similarly, Muro-Arena et al. (2018) reproduced the ALMA and SPHERE observations of HD 163296 with the radiative transfer code MCMax3D (Min et al. 2009) to investigate the distinct radial morphologies of the disk in the millimeter and NIR regimes. Rich et al. (2021) focused instead on the relative heights of the scattering surface and the CO emission surface, using them to constrain the dust-to-gas ratio in HD 163296, HD 97048, and IM Lup. These studies differ from ours because the dust properties were fixed; in particular, the grain size distribution was assumed a priori, while other physical parameters such as the spatial distribution of grains were varied to reproduce the observations. However, from the assumed grain size distribution, the small dust ratio is 0.08 for HD 163296 and 0.02 for J1852, higher and not compatible with what we found. By contrast, Franceschi et al. (2023) employed RADMC-3D¹ to constrain the turbulence and the vertical and radial distribution of large and small grains in IM Lup, fitting for the grain size distribution. They found that reproducing the scattered-light emission required a mass fraction of micronsized grains of ∼10⁻⁴, slightly lower than our median value of 10⁻³.5±0.4⁵ between the rings, but still comparable to our result. Recently, Swastik et al. (2026) modeled LkCa15 simultaneously in the millimeter and scattered-light regimes, constraining the scattering phase function as part of the fit. They found that the model that best matches the observations is characterized by highly porous aggregates, with a porosity p= 87% and a local grain-size distribution exponent of 2.3 between 12 µm and 2 mm. We can compare this with what we found for LkCa15, considering Fig. 10 we find that for q = 2.3 and p = 90% this corresponds to a_max ≈ 5 mm, which is comparable with the result of Swastik et al. (2026).

It is worth highlighting that our approach differs from those presented in the cited works because we have focused solely on the height of the scattering surface, using a semi-analytical modeling framework, rather than on reproducing in detail the intensity of the emission. Accepting this downside, our method does not require detailed modeling of individual disks and, therefore, it can be efficiently extended to larger samples.

6 Conclusions

We used scattered-light surface height measurements to directly constrain the mass of small dust grains in protoplanetary disks. To do so, we developed a semi analytical model for the scattering surface. Starting from radiative transfer considerations, we modeled the path of the stellar radiation and its interaction with the disk to understand what is the basic physics that sets the height of the surface. The results of our theoretical investigation are as follows:

• The scattering surface coincides with the τ_* = 1 surface, and the scattering height can be expressed as a function of the disk’s parameters, z₀(r₀,κ, M_{small dust}, T(r, z)). This provides a direct link between observable surface heights and physical disk parameters;

• The determination of the scattering surface is a nonlocal problem. The height of the surface at a given radius depends not only on the local column density, but also on the integrated density along the entire optical path from the star to the scattering point.

• We tested this prediction against the radiative transfer model MCFOST and consistency was found;

• When considering disks in scattered light, the isothermal approximation is not correct. Ignoring the temperature gradient (a common simplification) leads to significant underestimation of the surface height (>30%). This necessitates spatially resolved temperature measurements for accurate mass estimates;

• In contrast, dust settling and anisotropic scattering have negligible effects on the scattering surface location and can be safely ignored. Multiple scattering becomes important only for high inclinations (i ≥ 60°).

We then approached the inverse problem: using measured scattering surface heights to constrain the properties of small dust grains (a ≤ 1 μm). The inversion is degenerate in terms of dust opacity and mass, so by assuming an opacity model, we were able to estimate the small dust mass for a sample of ten pro-toplanetary disks. From these measurements, we were able to obtain the following results:

• We used the combined scattered-light data (probing small grains) with submm data (probing large grains) to derive the small dust ratio: the fraction of small dust mass- finding typical values around 10⁻³;

• By requiring self-consistency between scattered-light (small dust) and submm (total dust) observations, we constrained the grain size distribution. For moderate grain growth a_max ≥ 0.1 mm, we retrieved common power-law indices of q ~ 33.5;

• Dust porosity significantly affects these constraints and represents a major uncertainty. This can be reduced through scattering phase function analyses, which independently constrain porosity and grain structure, as in Tazaki et al. (2023); Ginski et al. (2023).

Taken together, these results show that scattered-light surface heights, combined with thermal structure measurements, can provide constraints on the small dust grain population in pro-toplanetary disks, in particular on its mass. By combining this diagnostic with submm data, one can obtain constraints on the grain size distribution. While porosity remains one of the major sources of uncertainty, this method, combined with scattering phase function analyses, can provide a statistical characterization of the small dust population in disks.

Acknowledgements.

We thank Anibal Sierra for kindly providing us with the MAPS fits used in this work. We acknowledge support from the European Union (ERC Starting Grant DiscEvol, project number 101039651) and from Fondazione Cariplo, grant No. 2022-1217. MB has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement No. 101002188). S.F. acknowledges financial contributions from the European Union (ERC, UNVEIL, 101076613) and from PRIN-MUR 2022YP5ACE. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

Appendix A Nonisotropic scattering

Fig. A.1 Scattering height considering forward scattering and isotropic scattering.

Here, we test the isotropic scattering approximation by verifying whether considering nonisotropic scattering, with various degree of anisotropy, significantly affects the location of the scattering surface. Scattering anisotropy is described by the scattering phase function, i.e. the probability of scattering into the direction μ = cos ψ, where ψ is the deflection with respect to the photon incoming direction. If the photon is scattered backward ψ = π, if forward ψ = 0.

We use the Henyen-Greenstein parametrization of the phase function (Henyey & Greenstein 1941) via $$\beta (g,\psi )=\frac{1}{4\pi }\frac{1-g^2}{(1+g^2-2g\cos \psi {)}^{3/2}},$$(A.1)

where g ∈ [-1,1] parametrizes the degree of anisotropy. Considering the scattering of the stellar light with the protoplanetary disk, the maximum scattering angle corresponds to the near side of the disk, and the maximum angle to the far side. Therefore, we expect to observe the maximum difference in the height between the far side and the near side of the disk. Consequently, assuming an anisotropy factor g intrinsic of the dust, we compare the scattering height between the far side and the near side and see if there is a significant difference with respect to the isotropic scattering case.

The contribution of the phase function can be inserted into Eq. 9, which becomes $$\frac{d{I}_{\nu }}{ds}={\sigma }_{\nu }\beta (g,\psi )B_{\nu }(T_{\ast })W(r)e^{-{\tau }_{\nu ,\ast }}{e}^{-{\tau }_{\nu }(s)}.$$(A.2)

We compute the scattering surface for various degrees of anisotropy g, both in the case of forward and backward scattering. Figure A.1 shows an example in the case of forward scattering for i = 20°. The scattering surface is essentially unaffected by the phase function, even with high degrees of anisotropy. The reason behind this irrelevance is that τ₁ is such a steeply varying function of z (Fig. 3) that the phase function does not move the τ₁ ≈ 1 surface by much. In general, the scattering surface depends on azimuth, but the difference between the two sides, and the difference between the non isotropic and isotropic scattering is at most 2 AU, so it is negligible, compared to the typical errors in measuring the scattering height. In the case of higher inclinations this effect can become more important, however we already discussed that our model breaks down at high inclinations. It is therefore reasonable to ignore the scattering anisotropy.

Appendix B Isothermal approximation

Here we test whether it would be reasonable to assume that the disk is vertically isothermal. We model the disk as isothermal: the scale height can be computed from the midplane temperature as in Eq 11, and the density profile becomes $$\rho (z)={\rho }_0\exp [-\frac{z^2}{2H^2(r)}].$$(B.1)

Fig. B.1 Comparison between the isothermal and nonisothermal case. Top : Vertical density profile at r = 100 AU. Bottom : Scattering surface.

Figure B.1 shows the comparison between the stratified temperature case and the isothermal approximation for LkCa15: in the top panel we show the vertical density profile at 100 AU, and we point out that in the non isothermal case the disk is more flared. Consequently the scattering surface is higher, as shown in the second panel of Fig. B.1. The difference accounts for more than 20 AU, so the isothermal approximation while convenient is not adequate. Practically speaking, this means that when modeling the observations described in Sect. 4, we would require observational measurements of the temperature structure.

Appendix C Opacity

When inverting the model, there is a clear degeneracy between dust mass and opacity. In our work, we used the Ossenkopf & Henning (1994) opacity model. Figure C.1 shows the dust opacity in the J (λ ≈ 1.2, µm), H (λ ≈ 1.6, µm), and K (λ ≈ 2.2, µm) bands, computed with Optool (Dominik et al. 2021), considering grains up to 1 µm. The comparison includes five widely adopted opacity models with different compositions and porosities (Draine & Lee 1984; Ossenkopf & Henning 1994; Pollack et al. 1994; Kataoka et al. 2015; Min et al. 2016).

This comparison highlights the significant uncertainties associated with dust opacities, which can vary between models. As a result, it is generally not meaningful to directly compare the scattering heights derived from observations of different disks taken in different bands, since the variations in opacity can dominate over wavelength-dependent effects. However, if multiband scattered-light observations are available for the same disk, differences in the derived scattering height across bands can provide valuable constraints on the grain size distribution.

Fig. C.1 Opacity for different dust models.

Appendix D Grain size distribution plots

Here, we present the constraints on the grain size distribution for the rest of the disks.

Fig. D.1 f = 1 curve as computed in Sect. 5.3 for GM Aur.

Fig. D.2 f = 1 curve as computed in Sect. 5.3 for IM Lup.

Fig. D.3 f = 1 curve as computed in Sect. 5.3 for J1615.

Fig. D.4 f = 1 curve as computed in Sect. 5.3 for J1852.

Fig. D.5 f = 1 curve as computed in Sect. 5.3 for HD 163296.

Fig. D.6 f = 1 curve as computed in Sect. 5.3 for HD 97048.

Fig. D.7 f = 1 curve as computed in Sect. 5.3 for MWC 480.

Fig. D.8 f = 1 curve as computed in Sect. 5.3 for SY Cha.

Fig. D.9 f = 1 curve as computed in Sect. 5.3 for HD 34282.

Appendix E Surface density profile

In Eq. 20, we assume that the surface density profile scales as r⁻¹ before the exponential cutoff. Here, we quantify the impact of this assumption. Considering a more general Lynden-Bell & Pringle (1974) profile, $$\mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0{\left(\frac{r}{rc}\right)}^{-\gamma }\exp [-{\left(\frac{r}{r_c}\right)}^{2-\gamma }],$$(E.1)

we can compute the small dust mass varying the exponent γ between 0.5 and 1.5. The mass difference amounts to a factor of ≈2.5 and here we report the small dust mass measurements for γ = 0.5 and γ = 1.5.

All Tables

All Figures

	Fig. 1 Model of the path of a ray through the disk.
In the text

	Fig. 2 Differential intensity as a function of z/H for r = 150 AU.
In the text

	Fig. 3 Contributions to the differential intensity as a function of z in the disk at fixed radius r = 150 AU. Left : optical depth τ1 and τ2 as a function of z/H. Middle : attenuation factors e−τ1 and e−τ2 as a function of z/H. Right : scattering cross-section, σν, as a function of z/H.
In the text

	Fig. 4 Scattering height at r = 100 AU for different α parameters in the isothermal and non-isothermal case, taking dust settling into account.
In the text

	Fig. 5 Comparison between the scattering surface from our semi-analytic model, the one computed with MCFOST, and the τ1 = 1 surface from MCFOST, at an inclination of i = 40°.
In the text

	Fig. 6 Top: difference in height between the scattering surface computed in the single scattering approximation or accounting for multiple scattering, at 200 AU. Bottom: optical depth from the scattering point to the observer as a function of disk’s inclination.
In the text

	Fig. 7 Top: measured scattering surfaces from the literature. Middle: small dust masses as a function of radius inferred from the inverse problem in our work. Bottom: small dust mass fraction for all the disks in the sample.
In the text

	Fig. 8 Constraints on the maximum grain size and power-law exponent for LkCa15 obtained using different opacity models.
In the text

	Fig. 9 Constraints on the maximum grain size and power-law exponent for LkCa15 using the Pollack et al. (1994) composition and varying porosity.
In the text

	Fig. 10 Constraints on the maximum grain size and power-law exponent for LkCa15 using Draine & Lee (1984), Pollack et al. (1994), and Ossenkopf & Henning (1994) opacities. The dotted lines represent the same composition, but considering a porous aggregates scenario, as big grains have a high porosity (p ~ 0.8).
In the text

	Fig. A.1 Scattering height considering forward scattering and isotropic scattering.
In the text

	Fig. B.1 Comparison between the isothermal and nonisothermal case. Top : Vertical density profile at r = 100 AU. Bottom : Scattering surface.
In the text

	Fig. C.1 Opacity for different dust models.
In the text

	Fig. D.1 f = 1 curve as computed in Sect. 5.3 for GM Aur.
In the text

	Fig. D.2 f = 1 curve as computed in Sect. 5.3 for IM Lup.
In the text

	Fig. D.3 f = 1 curve as computed in Sect. 5.3 for J1615.
In the text

	Fig. D.4 f = 1 curve as computed in Sect. 5.3 for J1852.
In the text

	Fig. D.5 f = 1 curve as computed in Sect. 5.3 for HD 163296.
In the text

	Fig. D.6 f = 1 curve as computed in Sect. 5.3 for HD 97048.
In the text

	Fig. D.7 f = 1 curve as computed in Sect. 5.3 for MWC 480.
In the text

	Fig. D.8 f = 1 curve as computed in Sect. 5.3 for SY Cha.
In the text

	Fig. D.9 f = 1 curve as computed in Sect. 5.3 for HD 34282.
In the text