One Solution to the Mass Budget Problem for Planet Formation: Optically Thick Disks with Dust Scattering

Consider a flat disk region with a uniform temperature of T. The intensity emitted by this region has been calculated by Miyake & Nakagawa (1993). Here, we follow the derivation given by Birnstiel et al. (2018) and extend it further to very optically thick cases.

The general radiative transfer equation is

$$\begin{eqnarray}&&\displaystyle \frac{1}{c}\displaystyle \frac{\partial {I}_{\nu }}{\partial t}+{\boldsymbol{n}}\cdot {\rm{\nabla }}{I}_{\nu }=-({\sigma }_{\nu ,a}+{\sigma }_{\nu ,s}^{\mathrm{eff}}){I}_{\nu }+{j}_{\nu }+{\sigma }_{\nu ,s}^{\mathrm{eff}}{J}_{\nu }\end{eqnarray} \tag{ 3 }$$

where I_ν(${\boldsymbol{x}}$, t, ${\boldsymbol{n}}$) is the intensity at the position ${\boldsymbol{x}}$, time t, and along the direction of ${\boldsymbol{n}}$. ${J}_{\nu }={\left(4\pi \right)}^{-1}\int {I}_{\nu }d{\rm{\Omega }}$ and j_ν/σ_ν,a = B_ν, while σ_ν,a and σ_ν,s^eff are the absorption and effective scattering opacity at the frequency of ν. This equation implicitly assumes that the scattering is isotropic. Because the scattering is not isotropic for the dust with sizes (s) 2πs ≫ λ, we use the effective scattering coefficient to approximate the anisotropic scattering effect with ${\sigma }_{\nu ,s}^{\mathrm{eff}}=(1-{g}_{\nu }){\sigma }_{\nu ,s}$ where g_ν is the usual forward-scattering parameter. This approximation is valid for the optically thick disk (Ishimaru 1978) that is the focus of this work.

Assuming that the disk surface follows the 1D plane atmosphere geometry and has a time-independent radiation field, the radiative transfer equation throughout the disk is simplified to

$$\begin{eqnarray}&&\mu \displaystyle \frac{{{dI}}_{\nu }}{{dz}}=-({\sigma }_{\nu ,a}+{\sigma }_{\nu ,s}^{\mathrm{eff}}){I}_{\nu }+{j}_{\nu }+{\sigma }_{\nu ,s}^{\mathrm{eff}}{J}_{\nu },\end{eqnarray} \tag{ 4 }$$

where μ = cos (θ) and θ is the angle between ${\boldsymbol{n}}$ and the vertical direction (the z-direction). The 1D plane atmosphere geometry can be justified considering that the radio emission comes from a thin disk midplane (Pinte et al. 2016).

If we adopt $d{\tau }_{\nu }=-({\sigma }_{\nu ,a}+{\sigma }_{\nu ,s}^{\mathrm{eff}}){dz}$, we have

$$\begin{eqnarray}&&\mu \displaystyle \frac{{{dI}}_{\nu }}{d{\tau }_{\nu }}={I}_{\nu }-{S}_{\nu }\end{eqnarray} \tag{ 5 }$$

with

$$\begin{eqnarray}&&{S}_{\nu }=(1-{\omega }_{\nu }){B}_{\nu }(T)+{\omega }_{\nu }{J}_{\nu }({\tau }_{\nu }),\end{eqnarray} \tag{ 6 }$$

where the single scattering albedo ${\omega }_{\nu }={\sigma }_{\nu ,s}^{\mathrm{eff}}/({\sigma }_{\nu ,a}+{\sigma }_{\nu ,s}^{\mathrm{eff}})$ and $1-{\omega }_{\nu }={\sigma }_{\nu ,a}/({\sigma }_{\nu ,a}+{\sigma }_{\nu ,s}^{\mathrm{eff}})$.

With the Eddington approximation, the second moment of the radiative transfer equation becomes

$$\begin{eqnarray}&&\displaystyle \frac{1}{3}\displaystyle \frac{{\partial }^{2}J}{\partial {\tau }_{\nu }^{2}}=(1-{\omega }_{\nu })({J}_{\nu }-{B}_{\nu }(T)).\end{eqnarray} \tag{ 7 }$$

If the temperature of the plane slab is a constant and there is no incoming radiation field at the upper and lower disk surface, the solution of the equation can be derived using the two-stream approximation (Miyake & Nakagawa 1993) as

$$\begin{eqnarray}&&\begin{array}{l}{J}_{\nu }({\tau }_{\nu })={B}_{\nu }(T)\\ \times \left(1-\displaystyle \frac{{e}^{-\sqrt{3(1-{\omega }_{\nu })}{\tau }_{\nu }}+{e}^{\sqrt{3(1-{\omega }_{\nu })}({\tau }_{\nu }-{\tau }_{\nu ,d})}}{{e}^{-\sqrt{3(1-{\omega }_{\nu })}{\tau }_{\nu ,d}}(1-\sqrt{1-{\omega }_{\nu }})+(\sqrt{1-{\omega }_{\nu }}+1)}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where τ_ν,d and τ_ν are the total and variable optical depth in the vertical direction.

With J_ν known, we can integrate Equation (5) throughout the disk to derive the emergent intensity (${I}_{\nu }^{\mathrm{out}}$). Based on the Eddington–Barbier approximation, the solution (Birnstiel et al. 2018) is

$$\begin{eqnarray}&&{I}_{\nu }^{\mathrm{out}}=(1-{e}^{-{\tau }_{\nu ,d}/\mu }){S}_{\nu }({\tau }_{\nu }=2\mu /3)\end{eqnarray} \tag{ 9 }$$

with S_ν(τ_ν) given in Equation (6) and J_ν(τ_ν) in Equation (8). If τ_ν,d < 4μ/3, τ_ν in S_ν is chosen as τ_ν,d/2.

With I_ν^out, we can define its deviation from the blackbody radiation using

$$\begin{eqnarray}&&\chi \equiv \displaystyle \frac{{I}_{\nu }^{\mathrm{out}}}{{B}_{\nu }}.\end{eqnarray} \tag{ 10 }$$

As alluded to in the introduction, χ < 1 can be due to either emission from the optically thin region or dust scattering in the optically thick region.

If we choose ${\tau }_{\nu }=2\mu {\tau }_{\nu ,d}/(3{\tau }_{\nu ,d}+1)$ in Equation (9) to approximate both the optically thick and thin cases, Equations (9) and (10) can be written out explicitly as

$$\begin{eqnarray}&&\begin{array}{l}\chi \equiv \displaystyle \frac{{I}_{\nu }^{\mathrm{out}}}{{B}_{\nu }}=(1-{e}^{-{\tau }_{\nu ,d}/\mu })\\ \times \,\left(1-{\omega }_{\nu }\tfrac{{e}^{-\sqrt{3(1-{\omega }_{\nu })}{\tau }_{\nu }}\,+\,{e}^{\sqrt{3(1-{\omega }_{\nu })}({\tau }_{\nu }-{\tau }_{\nu ,d})}}{{e}^{-\sqrt{3(1-{\omega }_{\nu })}{\tau }_{\nu ,d}}(1-\sqrt{1-{\omega }_{\nu }})+(\sqrt{1-{\omega }_{\nu }}+1)}\right)\\ \qquad \mathrm{with}\,{\tau }_{\nu }=2\mu {\tau }_{\nu ,d}/(3{\tau }_{\nu ,d}+1).\end{array}\end{eqnarray} \tag{ 11 }$$

For the optically thin region, Equations (9) or (11) reduces to ${I}_{\nu }^{\mathrm{out}}\to (1-{\omega }_{\nu }){\tau }_{\nu ,d}{B}_{\nu }/\mu$. The quantity of $(1-{\omega }_{\nu }){\tau }_{\nu ,d}$ is basically the disk optical depth calculated with the absorption coefficient or ${\tau }_{\nu ,d}^{\mathrm{abs}}$. Thus, when the disk is optically thin (τ_ν,d < 1), the emergent intensity reflects only the absorption opacity.

On the other hand, when the disk is optically thick with (1 − ω_ν)τ_ν,d ≫ 1, the emission for a disk with scattering is smaller than a disk without scattering (the blackbody radiation) by a factor of

$$\begin{eqnarray}&&\chi \equiv \displaystyle \frac{{I}_{\nu }^{\mathrm{out}}}{{B}_{\nu }}=1-\displaystyle \frac{{\omega }_{\nu }}{(\sqrt{1-{\omega }_{\nu }}+1){e}^{\sqrt{3(1-{\omega }_{\nu })}\cdot 2\mu /3}}.\end{eqnarray} \tag{ 12 }$$

Because $\chi \lt 1$ if ω_ν > 0, the optically thick scattering disk looks fainter than the blackbody radiation calculated using the same disk temperature.

If we expand the exponent in Equation (12) with the Taylor series, we can simplify the equation further to

$$\begin{eqnarray}&&\chi =\displaystyle \frac{\sqrt{3}+2\mu }{2\mu +\tfrac{\sqrt{3}}{\sqrt{1-{\omega }_{\nu }}}}.\end{eqnarray} \tag{ 13 }$$

In the extreme case with ${\omega }_{\nu }\to$ 1, Equation (13) becomes $\chi =(1+1.15\mu )\sqrt{1-{\omega }_{\nu }}\sim \sqrt{1-{\omega }_{\nu }}$, which is similar to the result based on the mean free path argument in the introduction. We have verified that Equation (13) only deviates from Equation (12) by less than 5% over the whole parameter space. This enables us to solve ω_ν analytically using χ, as

$$\begin{eqnarray}&&{\omega }_{\nu }=1-{\left(\displaystyle \frac{1}{\chi }+\left(\displaystyle \frac{1}{\chi }-1\right)\displaystyle \frac{2\mu }{\sqrt{3}}\right)}^{-2}.\end{eqnarray} \tag{ 14 }$$

If we just apply Equation (1) to calculating the disk optical depth using the emergent intensity from an optically thick disk (Equation (12)), we will derive an optical depth of

$$\begin{eqnarray}&&{\tau }_{\mathrm{obs}}=-\mathrm{ln}(1-\chi ),\end{eqnarray} \tag{ 15 }$$

even if the disk is very optically thick. Thus, another explanation for τ_obs = 0.6 in the DSHARP disks is that these disks are actually very optically thick (e.g., τ = 10⁴) but with χ = 0.45 due to dust scattering (by plugging τ_obs = 0.6 into Equation (15)). With μ = 1 or 0.5, χ = 0.45 corresponds to ω_ν = 0.93 or 0.89, respectively (Equation (14)), suggesting that the dust in these disks is highly reflective.

Figure 1 shows how χ changes with ω using Equation (9). As expected, the intensity drops when the disk becomes more optically thin. On the other hand, even if the disk is optically thick, stronger scattering can also lead to a smaller intensity. In the optically thick limit (τ = 100), Equation (12) (yellow curves) agrees with the full solution (Equation (9)) very well. For the marginally optically thick or optically thin disks (the green and red curves), the intensity increases when the disk is more inclined (dashed curves). This is simply because our line of sight passes through more column (1/μ factor) when the disk is inclined. On the other hand, for very optically thick disks (τ ≳ 5), the intensity decreases when the disk is more inclined due to dust scattering. Thus, inclined optically thick disks look even fainter than face-on disks. The change in χ with respect to the disk inclination in the optically thick limit is shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To study the effect of dust scattering, we can also compare the intensity from disks having strong scattering with the intensity from disks having zero albedo. We thus define

$$\begin{eqnarray}&&{\chi }_{t}=\displaystyle \frac{{I}_{\nu }^{\mathrm{out}}}{{B}_{\nu }(1-{e}^{-{\tau }_{\nu ,d}})}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\chi }_{a}=\displaystyle \frac{{I}_{\nu }^{\mathrm{out}}}{{B}_{\nu }(1-{e}^{-{\tau }_{\nu ,d}^{\mathrm{abs}}})}.\end{eqnarray} \tag{ 17 }$$

These comparisons assume that μ = 1. Without dust scattering, we have χ_t = χ_a = 1. With dust scattering, we have χ_t = (1 − ω_ν) and χ_a = 1 in the optically thin limit. Figure 3 shows how both χ_t and χ_a change with τ_ν,d and ${\tau }_{\nu ,d}^{\mathrm{abs}}$ if ω = 0.9. When τ_ν,d < 1, χ_t = (1 − ω_ν) is a good approximation. When ${\tau }_{\nu ,d}^{\mathrm{abs}}\gt 1$, Equation (12) is a good approximation. Figure 3 also shows an interesting phenomenon: dust scattering does not always reduce intensity. When ${\tau }_{\nu ,d}^{\mathrm{abs}}\sim 1$, the disk with scattering is actually brighter than the disk without scattering, as long as these two disks have the same absorption optical depth. Due to this complex phenomenon at ${\tau }_{\nu ,d}^{\mathrm{abs}}\sim 1$, in the rest of this Letter when we call a disk "optically thin" we refer to τ_ν,d < 1, and when we call a disk "optically thick" we refer to ${\tau }_{\nu ,d}^{\mathrm{abs}}\gt 1$, unless otherwise stated.

Zoom In Zoom Out Reset image size

We also want to study how "wrong" the derived optical depth can be if we use the traditional method (Equation (1)) to measure the optical depth of a disk with scattering. We first calculate the intensity emitted by the disk with the absorption optical depth of ${\tau }_{\nu ,d}^{\mathrm{abs}}$ and the scattering albedo (ω_ν) using Equation (9). Then, we use Equation (1) to derive the observed optical depth (τ_obs), assuming that we know the actual disk temperature. The relationship between the observed disk optical depth and the true disk optical depth is shown in Figure 4 for different disk albedos. Clearly, even if ω is only 0.1, an extremely optically thick disk can be misidentified as a disk with the optical depth of order unity.

Zoom In Zoom Out Reset image size

To validate the approximations used in the previous section (e.g., the Eddington, two-stream, and Eddington–Barbier approximations), we have carried out direct radiative transfer calculations using the radiation module in Athena++ (Jiang et al. 2014). It solves the radiative transfer equation explicitly using the short characteristics method. Here we only solve the radiative transfer equation without evolving the hydrodynamics.

We set up a plane-parallel atmosphere with a density profile of

$$\begin{eqnarray}&&\rho ={\rho }_{0}{e}^{-({z}^{2}-{z}_{\min }^{2})/2{H}^{2}},\end{eqnarray} \tag{ 18 }$$

to represent the disk vertical density structure, where ρ₀ = 1 and H = 0.05 in the code unit. The simulation domain extends from the midplane at z_min = 0 to z_max = 0.35 with 256 uniform grid cells. For the radiation field, the reflecting boundary condition (which flips the z-direction of the intensity rays) has been adopted at the disk midplane z_min, considering that the disk is symmetric with respect to the midplane. The vacuum boundary condition has been adopted at z_max to simulate the outflowing radiation field. We vary the opacity to control the optical depth of the disk but keep ω = 0.9 for all the simulations. We solve the radiative transfer equation along 40 different angles.

Figure 5 shows the density structure and the radiation field from our fiducial simulation. The rightmost panel shows the intensity in the direction that is perpendicular to the disk surface. The red curves are from disks with ω = 0.9 and τ_ν,d = 40,000, while the black curves are from disks without scattering and τ_ν,d = 4000. Both disks have ${\tau }_{\nu ,d}^{\mathrm{abs}}=4,000$. The dotted curves in the middle panel are calculated using the analytical theory (Equation (8)). As clearly shown, the analytical theory reproduces the radiation field in the simulation very well, and scattering decreases the mean field J at the disk surface. The J panel also demonstrates why scattering can reduce the emergent intensity in the optically thick limit. Scattering couples the emergent intensity with J (Equations (6) and 9), which deceases at the disk surface. For the intensity coming out of the disk at other angles and the intensity from disks with different optical depths, the simulation data are plotted against the analytical theory in Figures 2 and 3. It seems that the analytical theory can explain the simulation results reasonably well, although it can underpredict χ by up to 15% in the optically thick limit. This is probably due to the approximations used in the analytical calculations.

Zoom In Zoom Out Reset image size

Our calculations above assume that the scattering is isotropic. To test anisotropic scattering, we have carried out MCRT calculations using RADMC-3D.¹¹ We set up the disk with Toomre Q = 1, which is the most massive disk possible. The disk temperature is prescribed as in Equation (2), using L_* = L_⊙. The disk is locally isothermal along the z-direction at a given r. We also assume M_* = M_⊙. The disk scale height H(r)/r is calculated from H(r)/r = c_s/v_ϕ where ${c}_{s}^{2}={RT}/\mu$ and μ = 2.35. To keep Q = 1 throughout the disk, the gas surface density

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{g}(r)=330{\left(\displaystyle \frac{r}{20\mathrm{au}}\right)}^{-1.75}{\rm{g}}\,{\mathrm{cm}}^{-2}.\end{eqnarray} \tag{ 19 }$$

To be in hydrostatic equilibrium, the disk has a Gaussian density profile along the z-direction,

$$\begin{eqnarray}&&{\rho }_{g}(r,z)=\displaystyle \frac{{{\rm{\Sigma }}}_{g}(r)}{\sqrt{2\pi }H(r)}\exp \left(-\displaystyle \frac{{z}^{2}}{2H{\left(r\right)}^{2}}\right).\end{eqnarray} \tag{ 20 }$$

The disk is truncated at r_in = 3 au and r_out = 150 au in the radial direction. At 20 au, H(r)/r = 0.047 and T = 28 K. The dust density is set to be 1/100 of the gas density and the dust scale height is one-fifth of the gas scale height. The DSHARP opacity (with water ice) is adopted, with s_min = 0.1 μm, s_max = 1 mm, and n(s) ∼ s^−3.5. The composition and optical constants are the same as in Table 1 and Figure 2 of Birnstiel et al. (2018), and can be obtained by dsharp_opac.¹² For full anisotropic scattering, the Müeller matrices are calculated using Mie theory, specifically the Bohren–Huffman program (Bohren & Huffman 1983). For isotropic scattering calculations, the opacity that is normalized to the dust density is ${\kappa }_{\mathrm{abs}}{| }_{d}$ = 2.1 cm² g⁻¹ and ${\kappa }_{\mathrm{sca}}{| }_{d}$ = 19.5 cm² g⁻¹ at 1.25 mm (so ω = 0.9). To compare with the rest of the Letter, where the opacity is normalized to the gas density, we can derive the gas density normalized opacity κ_abs = 0.021 cm² g⁻¹ and κ_sca = 0.195 cm² g⁻¹. Here, ρ_dκ_abs,d or 100ρ_dκ_abs are basically σ_ν,a in Section 2.1. For the face-on disk, τ ∼ 7 at 20 au. For every disk inclination, we have run three simulations: one without scattering, one with isotropic scattering, and one with full anisotropic scattering treatment. All of these three simulations have the same absorption opacity. We have used 5 × 10⁸ photon packages. The resolution in the radial, poloidal, and azimuthal directions are 512, 2048, and 32 cells, respectively. The cell size in the radial direction is uniform in logarithmic space, while the cell size in the poloidal and azimuthal directions are uniform in linear space from 0° to 50°, and from 0 to 2π, respectively. The reflecting boundary condition is used at the disk midplane. Such high resolution in the poloidal direction is crucial for treating the scattering process properly.

The results are shown in Figure 6. The top panels show the 2D intensity maps at 1.25 mm for the disks with anisotropic scattering. Intensity maps from MCRT calculations without scattering and with isotropic scattering are also generated. We cut through the horizontal major axis in the images to derive the 1D profiles, which are shown in the middle and bottom panels. The middle panels show the 1D profiles of the brightness temperature that is converted from the measured intensity. When the disk is face on, $\chi \sim 0.6$; this is consistent with our analytical estimate using ω = 0.9. The bottom panels show the derived optical depths using Equation (1). Different colored curves show disks with different scattering treatments. We do not show the measured optical depth for no-scattering cases, because when the optical depth becomes very large (e.g., >10), Equation (1) cannot provide an accurate estimate of the optical depth. Compared with isotropic scattering, full anisotropic scattering treatment does not change the results qualitatively. Clearly, if the disks are highly optically thick but have scattering, the measured optical depths are ≲1 using Equation (1). If the scattering disk is more inclined (to the right panels), the measured brightness temperature is smaller and thus the derived optical depth becomes smaller.

Zoom In Zoom Out Reset image size

Figure 6 also suggests that even a Q = 1 disk will become optically thin beyond 50 au, where the brightness temperature decreases much faster than the midplane temperature. Thus, if the measured τ at the outer disk is small, the disk is probably truly optically thin instead of optically thick with strong scattering (the possible tests using the spectral index are presented in Section 3.3).

Our proposed optically thick disk with scattering scenario implies that protoplanetary disks can potentially hide a large amount of solids in disks (within 50 au for the Q ∼ 1 disk) from ALMA. Previous submm surveys may underestimate the solids in disks significantly. Class II disks may still have a significant amount of dust to form planetesimals or planets later. Such disks with more solids are more consistent with the amount of solid mass in exoplanets. Similar to our simulations, Evans et al. 2017 have carried out MCRT simulations for gravitationally unstable disks and concluded that such disks can hide a factor of 3–30 dust mass to ALMA observations, even if dust scattering has not been included in these calculations. Furthermore, the massive disks are more consistent with the fact that Class II disks are still accreting at moderate rates (Hartmann et al. 2006), implying at least ${10}^{-8}{M}_{\odot }\,{\mathrm{yr}}^{-1}\times {10}^{6}\mathrm{yr}\sim 0.01\,{M}_{\odot }$ gas in disks.

To study the true amount of dust in protoplanetary disks, we may need to go to longer wavelengths. The Very Large Array (VLA), Next Generation Very Large Array (ngVLA; Murphy et al. 2018), and Square Kilometre Array (SKA; Testi et al. 2015) will be quite powerful in this regard, although these long wavelength observations will get significant contamination from free–free emission from the star or jet. Assuming that the absorption opacity changes with the frequency as ${\kappa }_{\mathrm{abs}}=0.021{\left(\nu /240\mathrm{GHz}\right)}^{\beta }$ cm² g⁻¹, κ_abs is 3.6 × 10⁻⁵ cm² g⁻¹ at 10 GHz with β = 2. Using Equation (19), we can derive that the Q = 1 disk will be optically thin (${\tau }_{\nu ,d}^{\mathrm{abs}}\lt 1$) beyond 2 au at 10 GHz. On the other hand, if the dust is big (Section 3.3) and β = 1, κ_abs is 8.8 × 10⁻⁴ cm² g⁻¹ at 10 GHz, and the Q = 1 disk will be optically thin beyond 10 au at 10 GHz.

To illustrate how much dust mass can be hidden from ALMA observations at short millimeter wavelengths in our Q = 1 disk with dust scattering, we plot the measured and real dust mass in Figure 7. The disk intensity is from the MCRT calculations in Section 2.3. In addition to the 1.25 mm band in Section 2.3, we also carry out MCRT calculations at 7 mm band. Because the input disk is as massive as a disk can get (Q = 1), Figure 7 shows the maximum amount of dust mass that can be hidden from observations. When the disk is large (e.g., 100 au), most of dust mass is at the outer disk, which is optically thin at 1.25 mm; therefore, 1.25 mm observations will only underestimate the disk mass by a factor of ∼2 for these extended disks. When the disk is compact (<30 au), ALMA 1.25 mm observations can easily underestimate the real dust mass by a factor of 10. ALMA protoplanetary disk surveys suggest that most disks are actually compact disks (Figure 2 in Ansdell et al. 2016), which can be due to dust radial drift. Thus, these surveys may underestimate the dust mass significantly for the whole sample.

Zoom In Zoom Out Reset image size

On the other hand, for VLA 7 mm observations, the disk is optically thin beyond several au. Thus, VLA observations provide a much more accurate mass estimate. We notice that the estimated dust mass is slightly larger than the real dust mass when the disk is optically thin, which is due to the fact that the brightness temperature is slightly higher with scattering included (Figure 6). We suspect that this is due to the intensity-enhancing effect in Figure 3. Recent VLA observations by Tychoniec et al. (2018) suggest that Class 0/I objects are much more massive than Class II objects. However, we caution that this large difference may be due to the fact that observations for Class 0/I and II disks are carried out at different bands. Similar VLA surveys for Class II disks are required to probe the real dust mass in these disks.

Previously, it has been proposed that we can measure the disk temperature by observing the dust continuum at higher frequencies (e.g., ALMA Bands 8, 9), because the disk is optically thick at those bands (e.g., Kim et al. 2019). We point out that if the disk has scattering, then we cannot use this method to measure the disk temperature as the intensity from thermal radiation is affected by the scattering. Furthermore, we cannot use ALMA to measure the disk mass accurately if the disk is optically thick.

On the other hand, the simple relationships between the reduced emission (χ) and albedo (Equations (12)–(14)) provide an unique opportunity to study dust properties at the τ ∼ 1 surface. These relationships are independent of any particular disk model (e.g., whether turbulent or not), which is why they are powerful enough to constrain the dust properties. We can measure χ if we know the disk temperature and use χ to constrain the dust albedo directly. If we have multi-band observations, we can also use the spectral index to constrain the change in albedo (Section 3.3).

Although optically thick disks with ω ∼ 0.9 can explain τ_obs ∼ 0.6 in DSHARP observations, it is crucial to understand if dust in protoplanetary disks can have such a high albedo. Assuming that the dust follows the n(s) ∝ s^−3.5 size distribution with s_max = 1 mm, we use the DSHARP opacity (Birnstiel et al. 2018) to calculate the dust opacity and albedo at different wavelengths, as shown in Figure 8. Clearly, the albedo can be as high as 0.9 for radio observations at mm to cm. Dust scattering is most efficient when 2πs ∼ λ. With s_max = 1 mm, the strongest scattering (ω ∼ 0.96) occurs at 4 mm. With these albedos, we can calculate χ and τ_obs at different wavelengths, assuming that the disk is optically thick (with Equations (12) and (15)), as shown in the bottom panels of Figure 8. Clearly, for mm-cm observations (shaded region), χ is less than 0.7 and τ_obs is less than 1 with our assumed dust population. Thus, this assumed dust population can naturally explain the τ_obs < 1 in DSHARP observations.

Zoom In Zoom Out Reset image size

On the other hand, we caution that the DSHARP opacity adopted here is for compact spheres without porosity, and that the mixture of different compositions are handled with the Bruggeman rule. Changing composition or porosity can change the opacity and albedo dramatically, as alluded to in Birnstiel et al. (2018). For example, the dust becomes more reflective with more water ice. Porosity reduces the resonant opacity features, making the opacity curve smoother. Carbonaceous materials have a substantial effect on the dust opacity. The maximum absorption feature can be shifted to much longer wavelengths than 2πs_max with some choices of carbonaceous materials. On the other hand, the peak of albedo is still at 2πs_max with different choices of carbonaceous materials (Birnstiel et al. 2018), making both χ and α measurements (Section 3.3) less sensitive to the choices of carbonaceous materials. More detailed calculations exploring different dust compositions, porosities, and size distributions are needed in future.

Just based on the DSHARP opacity, we want to explore how we can constrain the dust size distributions using the measured χ and the spectral index α. The spectral index α is

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{\mathrm{ln}\left(\tfrac{{I}_{{\nu }_{1}}}{{I}_{{\nu }_{2}}}\right)}{\mathrm{ln}\left(\tfrac{{\nu }_{1}}{{\nu }_{2}}\right)}=\left\{\begin{array}{l}\tfrac{\mathrm{ln}\left(\tfrac{{\chi }_{{\nu }_{1}}{B}_{{\nu }_{1}}}{{\chi }_{{\nu }_{2}}{B}_{{\nu }_{2}}}\right)}{\mathrm{ln}\left(\tfrac{{\nu }_{1}}{{\nu }_{2}}\right)}\sim \tfrac{\mathrm{ln}\left(\tfrac{{\chi }_{{\nu }_{1}}}{{\chi }_{{\nu }_{2}}}\right)}{\mathrm{ln}\left(\tfrac{{\nu }_{1}}{{\nu }_{2}}\right)}+2,\,\mathrm{if}\,\mathrm{thick}\\ \tfrac{\mathrm{ln}\left(\tfrac{{\kappa }_{\mathrm{abs},{\nu }_{1}}{B}_{{\nu }_{1}}}{{\kappa }_{\mathrm{abs},{\nu }_{2}}{B}_{{\nu }_{2}}}\right)}{\mathrm{ln}\left(\tfrac{{\nu }_{1}}{{\nu }_{2}}\right)}\sim \tfrac{\mathrm{ln}\left(\tfrac{{\kappa }_{\mathrm{abs},{\nu }_{1}}}{{\kappa }_{\mathrm{abs},{\nu }_{2}}}\right)}{\mathrm{ln}\left(\tfrac{{\nu }_{1}}{{\nu }_{2}}\right)}+2,\,\mathrm{if}\,\mathrm{thin},\end{array}\right.\end{eqnarray} \tag{ 21 }$$

where the approximation on the right is with the Rayleigh–Jeans approximation. Clearly, the measured α depends on totally different quantities in the optically thick and thin regimes. In the more familiar optically thin regime, it depends on how the absorption coefficient changes with the frequency. With κ_abs ∝ ν^β, α is β + 2. In the optically thick regime, α depends on χ. Assuming that χ ∝ ν^γ, α is γ + 2. As shown in Figure 1, ω changes monotonically with 1/χ that changes as λ^γ. Thus, if the albedo ω increases with wavelength (γ > 0), α = γ + 2 measured at this wavelength span will be larger than 2, and vice versa. To be more specific, if the disk is optically thick, we can use the measured α from observations to derive γ. Then we can constrain the relationship between ω_ν1 and ${\omega }_{{\nu }_{2}}$ using

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\nu }_{1}}{{\nu }_{2}}\right)}^{\gamma }=\displaystyle \frac{2\mu +\tfrac{\sqrt{3}}{\sqrt{1-{\omega }_{{\nu }_{2}}}}}{2\mu +\tfrac{\sqrt{3}}{\sqrt{1-{\omega }_{{\nu }_{1}}}}},\end{eqnarray} \tag{ 22 }$$

which is based on Equation (13).

Figure 9 shows how κ_abs, χ and the spectral index α change with different dust populations for different ALMA and VLA bands, using the DSHARP opacity. The panels in the second row show χ in the optically thick limit. The third row shows α in the optically thin limit (thus denoted as α_thin), while the bottom row is α in the optically thick limit (α_thick). Here, we simply use β + 2 or γ + 2 to calculate α, and thus do not consider the deviation of the blackbody radiation from the Rayleigh–Jeans limit (Huang et al. 2018a). χ and both αs have sharp transitions around s_max ∼ λ/2π. When s ≪ λ/2π, scattering is not important and χ = 1, α_thick = 2, α_thin ∼ 4. When s ≳ λ/2π, χ becomes less than 1, α_thick becomes larger than 2, and α_thin becomes less than 4. When the dust population has a lot of small dust (left panels), χ behaves more like a step function. When the disk is more populated with bigger particles (moving to the right panels), χ and both αs have sharper peaks around the transition. The peak of α_thick is due to the rapid change of ω with wavelength for particles at bigger sizes. If ω increases with wavelength, χ will decrease with wavelength and α_thick measured at this wavelength span will be larger than 2, and vice versa. For example, for the n(s) ∝ s^−3.5 and s_max = 1 mm case in Figure 8, ω increases with λ at λ < 4 mm and decreases with λ at λ > 4 mm, so the measured α_thick should be larger than 2 using observational bands <4 mm and smaller than 2 using observational bands >4 mm. This is shown in Figure 9. At s_max = 1 mm, all α_thick curves for the n(s) ∝ s^−3.5 case have values larger than 2, except the green curve for 10-7 mm.

Zoom In Zoom Out Reset image size

Thus, in this optically thick disk with scattering scenario, α can be lower than 2, which is not likely to happen for optically thin disks (e.g., the panels in the third row of Figure 9). If observations show that α is less than 2, it could be a strong indication that the disk is optically thick and dust scattering plays an important role; albedo decreases with wavelength. Using the DSHARP opacity, α_thick panels in Figure 9 suggest that, if the observed α from ALMA bands at λ ≲ 3 mm is less than 2, then s_max is from 30 to 300 μm. When s_max is larger than λ/2π, the probed ω is on the left side of the ω peak in Figure 8 so that α becomes larger than 2. This transition is relatively quick. After an α peak around s_max ∼ 0.5λ, α plateaus around a value slightly larger than 2.

Consider a disk that is optically thick at the inner disk and optically thin at the outer disk: the spectral index will be around 2 at the inner disk and suddenly change to 3–4 when τ < 1. This is simply because α in the optically thick and thin regimes is determined by different physical mechanisms. To illustrate this point, we calculate I_ν/B_ν and α for a Q = 1 disk (Equation (19)) at both ALMA and VLA wavelength bands (Figure 10). At the inner disk where the disk is optically thick, χ and α flatten out. More disk region becomes optically thin with observations at longer wavelengths. At VLA bands of 7 mm and 1 cm, the disk is optically thin even down to 5 au as long as s_max ≲ 1 mm. With s_max = 1 cm, which has the maximum opacity at ∼1 cm, VLA observations can still probe the disk down to 20 au before the disk becomes optically thick. HL Tau observations by Carrasco-González et al. (2016) seem to indicate that α derived from ALMA bands changes from 2 to 3 relatively quickly from 40 to 60 au, which may indicate that HL Tau is optically thick within 40 au in ALMA observations. On the other hand, α derived from 3 mm ALMA and 7 mm VLA bands is larger than 2 throughout the disk, which indicates that the disk may be optically thin for VLA observations. Figure 10 also suggests that we may want to measure the spectral index using every combination of two different bands, because α can change dramatically at different wavelengths, especially at wavelengths close to the Mie resonances of the dust opacity.

Zoom In Zoom Out Reset image size

Note that each panel in Figure 10 assumes the same dust size distribution throughout the disk. In reality, dust size distributions vary both radially and vertically. When the disk is optically thick, the emission is determined by the τ_ν ∼ 1 surface. Consequently, χ and α measured at the inner optically thick disk only inform us the dust size distribution at the disk surface. The dust at the midplane could be a lot larger, which can only be probed by observations at longer wavelengths. When the disk is optically thin, the measured α informs us the dust size distribution at the midplane directly (see Section 3.5 and Figure 11).

Zoom In Zoom Out Reset image size

Previous spatially resolved α measurements have shown that α decreases toward the inner disk (Pérez et al. 2012, 2015; Carrasco-González et al. 2016; Tazzari et al. 2016). Yet α is always larger than 2 in these observations. TW Hya has some indications that α can be smaller than 2 at inner 20 au (Tsukagoshi et al. 2016; Huang et al. 2018a). Although this deviation from 2 can in part be explained by the fact that the Rayleigh–Jeans approximation deviates from the blackbody radiation there (Huang et al. 2018a), we discuss the possibility that the disk has strong scattering and is optically thick at the inner disk. Under the optically thick scattering disk scenario, α < 2 within 20 au implies that we have a large population of dust with ∼100 μm sizes there at the τ_ν ∼ 1 surface (based on the blue curves in the bottom row of Figure 9). At 25 au, observations suggest that there is an α peak reaching α = 2.5, which is similar to the α_thick peak in Figure 9. Thus, one explanation for this peak is that dust size increases with radius, or that we are probing deeper large dust layers in the disk as the disk starts to become optically thin; at 25 au there is a large population of ∼300 μm dust leading to the α_thick peak. Further out at 30 au, the dust becomes even bigger and α_thick decreases and plateaus, which is similar to the observations. The increase of α at the outer disk may be due to the change in dust size or the whole disk becoming optically thin.

Although this story may be too complicated, it indeed highlights that dust scattering in optically thick disks can also lead to gaps and rings. ALMA continuum observations have revealed many gaps and rings in protoplanetary disks. Although some of these features are very prominent, some of them are very weak, with only ∼20% fluctuations (e.g., Huang et al. 2018b). These shallow features may also be explained by optically thick disks with radially varying dust scattering properties. If the dust becomes more reflective at a particular distance from the star, the disk will look like it has a gap there. If the dust becomes less reflective there, the disk will look like it has a ring. Changing the scattering properties means changing the dust compositions or distributions at the τ ∼ 1 surface, so the spectral index should also change at these gaps/rings. The spectral index can be either higher or lower at the rings, depending on how the particle size changes (e.g., Figure 9). On the other hand, we caution that dust scattering is unlikely to explain deep gaps (more than a factor of 10 deep) observed in some systems (e.g., AS 209 Guzmán et al. 2018), as a factor of 10 intensity reduction requires the albedo of 0.998 (based on Equation (14)), which is extremely high.

Recent surveys (Tripathi et al. 2017; Andrews et al. 2018b) confirm the previous hinted linear relationship (Andrews et al. 2010) between the submm continuum luminosity and its emitting surface area. Andrews et al. (2018b) also confirmed the disk luminosity and stellar luminosity relationship found in Andrews et al. (2013), Ansdell et al. (2016), and Pascucci et al. (2016). Andrews et al. (2018b) explored the optically thick disk scenario and found that the shapes of both relationships can be reproduced under this scenario. However, the optically thick disks generate too much emission. To reduce the luminosities of optically thick disks so that they are consistent with observations, a filling factor of 0.3 in the disk is needed. Tripathi et al. (2017) and Andrews et al. (2018b) suggested that substructures (e.g., rings, gaps) can lead to this filling factor. Here, we suggest that, in addition to substructures, dust scattering can also decrease the luminosity for optically thick disks. Instead of the filling factor, a high albedo at the τ_ν ∼ 1 surface may also explain the observations.

Reducing blackbody intensity due to scattering is known in various astronomical communities (e.g., electron scattering reduces the radiation from accretion disks around compact objects). For the protoplanetary disk study, Miyake & Nakagawa (1993) solved the radiative transfer equation for an isothermal disk (as summarized in Section 2). These results have been mentioned in many subsequent works (e.g., D'Alessio et al. 2001; Birnstiel et al. 2018). Sierra et al. (2017) applied this to dusty vortices and showed that the optically thick vortex center becomes fainter if dust scattering is considered.

On the other hand, dust scattering is largely ignored in radio observations because the protoplanetary disk is thought to be optically thin so that dust scattering is not important. Recent works by Kataoka et al. (2015) and Yang et al. (2016a) have suggested that dust scattering may be crucial for explaining submm polarization measurements, although other mechanisms may still be needed to explain the observations (Yang et al. 2016b; Kataoka et al. 2017). However, there is a strong tension between the dust size constrained by polarization measurements and submm-cm continuum spectral index measurements (Kataoka et al. 2016). Here, we suggest that such tension may be due to the optically thin assumption in both polarization and submm continuum studies. If the disk is optically thick, the spectral index is normally smaller than that from an optically thin disk. Assuming that the disk is optically thin, submm continuum observations can overestimate the particle size significantly. In reality (the schematic diagram from Figure 11), the small α by submm observations could simply reflect that the disk is optically thick for these observations and dust at the τ_mm ∼ 1 surface has a typical size of 0.1–1 mm with strong scattering. For longer wavelength observations by VLA (Pérez et al. 2015), the disk is likely to be optically thin beyond 10 au. The small α measured from these observations could indeed imply the presence of large particles at deeper layers (likely the midplane) in the disk where VLA is probing, which is a natural outcome from dust settling. Note that the absorption opacity at 7 mm is 10–50 times smaller than the opacity at 1.25 mm. Thus, VLA probes a much deeper layer in the disk. MCRT calculations for disks with such vertically varied dust size distribution will be presented in S. Zhang et al. (2019, in preparation).

Another line of evidence showing that most protoplanetary disks are optically thick at ALMA bands is that inclined disks are systematically less massive than face-on disks based on ∼1 mm observations, assuming that the disks are optically thin (Figure 4 in Garufi et al. 2018). If the disks are optically thin, the measured dust mass should not depend on the disk inclination. Thus, this inclination dependence suggests that the disks are optically thick at ALMA bands. Furthermore, with scattering included, we suggest that inclined optically thick disks will look even fainter (Figure 2). This effect may explain the very low temperature derived from the edge-on disk Flying Saucer (Guilloteau et al. 2016).

Although the optically thick scattering disk scenario seems to be promising for explaining several observations, there is evidence that the disks are not optically thick everywhere. The first is that the rings beyond 40 au in the DSHARP sample can be well fitted with Gaussian profiles along the radial direction, instead of flat-topped profiles (Dullemond et al. 2018). Because the distribution of dust trapped in turbulent disks with rings should follow a Gaussian profile, the intensity profile will have a flat top instead of a Gaussian profile if these rings are optically thick at the ring center. The second line of evidence against the optically thick rings is that CO emission coming from the back side of HD 163296 is dimmer at the location of the bright dust rings at 67 and 100 au (Isella et al. 2018). This dimming is due to dust extinction by the rings. Because the CO emission does not disappear completely, the dust rings cannot be optically thick. However, these two lines of evidence only apply to the rings at the outer disk beyond 40 au. This is actually consistent with our Figure 6, where even a Q = 1 disk will be optically thin beyond 50 au. The optically thick scattering disk scenario only applies to the inner disk within 50 au. We need to carry out similar tests or use other methods (e.g., Harsono et al. 2018; Powell et al. 2017, 2019; Casassus et al. 2019) to study the optical depth and mass of the inner disk in future.

One assumption in this work is that the disk is isothermal in the vertical direction. This is not quite correct, as the dusty disk intercepts the stellar irradiation causing the temperature inversion at the disk atmosphere (Calvet et al. 1991). On the other hand, large dust in protoplanetary disks settles to the disk midplane (probed by radio observations), while small dust is still suspended at the disk atmosphere (probed by near-infrared observations). Small dust intercepts the stellar irradiation, which determines the disk temperature structure. Small dust radiates energy vertically toward the disk midplane to warm up the large dust. Because large dust sits in the thermal bath generated by small dust, we expect that it should be approximately isothermal. The MCRT simulations for such disk configuration will be presented in S. Zhang et al. (2019, in preparation). Furthermore, the emission reduction argument in the abstract and our preliminary simulations show that, if the disk is optically thick but not isothermal, the emission reduction still applies and is mainly determined by the dust scattering properties at τ_ν ∼ 1. Thus, if the disk that has large dust is optically thick but not isothermal, the emission reduction probes the albedo of the τ_ν ∼ 1 surface.