The Disk Substructures at High Angular Resolution Project (DSHARP). V. Interpreting ALMA Maps of Protoplanetary Disks in Terms of a Dust Model

Dust thermal (sub-)millimeter emission from the outer regions ($r\gtrsim 10\,\mathrm{au}$) of protoplanetary disks has traditionally been considered to be optically thin, because it would require implausible amounts of dust mass to make the disk optically thick at these wavelengths out to many tens of au. Even if it were optically thick, it would produce much higher disk-integrated flux values than are observed (Ricci et al. 2012). If the assumption of low optical depth were true, it would aid the interpretation of submillimeter continuum maps in terms of the properties and dynamics of the dust grains, because the observed intensity I_ν would be directly proportional to the underlying dust surface density Σ_d.

Observational results from the past decade have shown that protoplanetary disks do not have simple monotonically decreasing surface density profiles, but consist of multiple narrow rings (e.g., ALMA Partnership et al. 2015; Andrews et al. 2016; Isella et al. 2016; Fedele et al. 2017, 2018; Huang et al. 2018), apparently single massive rings (e.g., Brown et al. 2009; Casassus et al. 2013), but also often non-axisymmetric features such as lopsided rings (e.g., van der Marel et al. 2013) and spirals (e.g., Pérez et al. 2016). These narrow or compact features present concentrations that may be optically thick or moderately optically thick. In between these features, the material is often optically thin.

This is both a curse and a blessing. Such optical depth effects make the interpretation of the data more difficult. In particular, the spectral slope variations at (sub-)millimeter wavelengths (${I}_{\nu }\propto {\nu }^{{\alpha }_{\mathrm{mm}}}$) across rings and gaps are strongly affected, perhaps even dominated, by these effects. But optical depth effects also provide new opportunities to measure the properties of the dust. An example of this is the scattering of its own thermal emission, and the induced polarized millimeter emission (Kataoka et al. 2015). Another example is when dust rings extinct part of the CO line emission from the back side of a disk (Isella et al. 2018).

But there is, of course, a major uncertainty in the dust opacity law. This is a long-standing problem (Beckwith & Sargent 1991) that still has not been resolved. This is in part because dust in protoplanetary disks, in particular the large grains seen as settled grains in a thin midplane layer, is very different from the dust in the interstellar medium. However, it is also due in part to uncertainties in the numerical and conceptual challenges in computing opacities. When comparing computed opacities with laboratory-measured opacities in the millimeter range, one often sees discrepancies of up to a factor of 10 or more (e.g., Demyk et al. 2017b).

By inspecting the (sub-)millimeter spectral slope ${\alpha }_{\mathrm{mm}}$ we can learn something about the grain size distribution and the opacity law (e.g., Beckwith & Sargent 1991; Testi et al. 2003; Wilner et al. 2005). With high angular resolution observations this can now be done as a function of radial coordinates in the disk (e.g., Isella et al. 2010; Guilloteau et al. 2011; Pérez et al. 2012), which shows that bright rings have shallower slopes than the dark annuli between them (e.g., ALMA Partnership et al. 2015; Tsukagoshi et al. 2016; Huang et al. 2018). However, changes in the spectral slope can also be caused by optical depth effects. The results of the Disk Substructures at High Angular Resolution Project (DSHARP; Andrews et al. 2018) show that these dust rings often have an optical depth close to unity (Dullemond et al. 2018; Huang et al. 2018; Isella et al. 2018).

With the detailed spatial information of the substructures in millimeter continuum and line maps of protoplanetary disks that the Atacama Large Millimeter/submillimeter Array (ALMA) is providing, it becomes increasingly important to discuss the details of the opacities used and the methods applied to translate the observations into information about the underlying dust grains.

This Letter is meant to give an overview of the methods used and choices made by the DSHARP collaboration to make this translation. We do not claim in any way that our choices and methods are better than those used by others, nor that we can resolve any of the uncertainties of the opacities. Instead, this Letter describes our methods and choices, and we provide an easy-to-use python module and a set of example calculations for the reproduction of these opacities and variants of them, as well as for handling some of the optical depth effects discussed in this Letter.

The structure of this Letter is as follows. In Section 2 we discuss our choices for computing the dust opacities, and present our Python tools that are publicly available. In Section 3 we apply these opacities to simple size distribution models, starting with a standard power-law model, ending in Section 4 with the analytic steady-state dust coagulation/fragmentation size distribution of Birnstiel et al. (2011, henceforth B11), again including the corresponding Python scripts. Finally, in Section 5 we present a very simple model to link the observed thermal emission to the observed extinction of back-side CO line emission (Isella et al. 2018).

As discussed above, the goal of this Letter is not to provide a "better" dust opacity model, but instead a transparent model based on open-source software that is easy to reproduce or to modify. To this end, we follow seminal works regarding protoplanetary disk composition and grain structure (Pollack et al. 1994, henceforth P94) that are widely used throughout the literature, and we use updated optical constants where available. To stay comparable to previously used opacities (and thus the resulting mass or surface density estimates), we chose to assume particles without porosity. This is a pragmatic choice instead of a realistic one, for protoplanetary disks because at least the initial growth phase involves larger porosities (Kempf et al. 1999; Ormel et al. 2007; Zsom et al. 2010; Okuzumi et al. 2012; Krijt et al. 2015).

P94 and subsequent work by D'Alessio et al. (2001) chose a mixture of water ice, astronomical silicates, troilite, and refractory organic material. The water fraction that was used in those works (around 60% by volume) was, however in disagreement with typically observed disk spectral energy distributions (SEDs) as pointed out by D'Alessio et al. (2006) and Espaillat et al. (2010), who reduced the water fraction to 10% of the value used in P94. As comets are thought to be a relatively pristine sample of the planet-forming material, we chose a water fraction of 20% by mass, in agreement with measurements of comet 67P/Churyumov–Gerasimenko (Pätzold et al. 2016).

To calculate mass absorption or scattering coefficients ${\kappa }_{\nu }^{\mathrm{abs}}$ or ${\kappa }_{\nu }^{\mathrm{sca}}$, we assume a vacuum as the embedding medium and furthermore need the complex refractive indices $m(\lambda )\,=n(\lambda )+i\,k(\lambda )$, which are functions of wavelength λ. These refractive index data will be called optical constants for simplicity.

The optical constants of water used in D'Alessio et al. (2001) were from Warren (1984), who gave tables for various temperatures. Henning & Stognienko (1996) used optical constants from Hudgins et al. (1993; amorphous ice at at 100 K) between 2.5 and 200 μm and the constants from P94 for the remaining wavelengths ranges. In this Letter, we use the more recent data from Warren & Brandt (2008). However, the differences with previous works are small (see Figure 1).

Zoom In Zoom Out Reset image size

The astronomical silicates used in D'Alessio et al. (2001) took a constant k value for λ > 800 μm. Henning & Stognienko (1996) argued (their Section 5.1) that $k\propto {\lambda }^{-1}$ is usually assumed, but that Campbell & Ulrichs (1969) indeed measured a high value of k = 0.05 at 2.7 mm (see Figure 1). Nevertheless, we use the opacities from Draine (2003) for astronomical silicates without increasing k.

For troilite and refractory organics, we use the optical constants from Henning & Stognienko (1996). For troilite, the constants are partly based on Begemann et al. (1994; in the range of 10–500 μm), with longer and shorter ranges taken from P94. Henning & Stognienko (1996) also performed Kramers–Kronig analysis on the organics optical constants of P94, which yielded little differences. The Henning & Stognienko (1996) data sets for troilite and refractory organics are available online, and are included in our opacity module with kind permission from Thomas Henning. The optical constants used in this Letter are shown in Figure 1.

Deriving optical constants for a mix of materials is a challenging task as it depends on the detailed structure of the composite particle. No general solution can be given. For more complex setups, computationally expensive numerical models need to be used. For some limiting cases analytical expressions can be given. These are typically called the Maxwell–Garnett rule (valid for inclusions in a background "matrix") and the Bruggemann rule, for a homogeneous mix without a dominant matrix. Details can be found in Bohren & Huffman (1998), whose notation we will follow.

If f_i denotes the volume fractions of the N inclusions ($i=1\ldots N$) and ${\epsilon }_{i}={m}_{i}^{2}$ the dielectric functions of the inclusions¹¹ (with refractive indices m_i), while f_m and ${\epsilon }_{m}$ are the corresponding values for the matrix, then the Maxwell–Garnett rule for spherical inclusions yields the mixed dielectric function

$$\begin{eqnarray}&&\bar{\epsilon }=\displaystyle \frac{(1-f)\,{\epsilon }_{m}+{\displaystyle \sum }_{i=1}^{N}{f}_{i}\,{\beta }_{i}\,{\epsilon }_{i}}{1-f+{\displaystyle \sum }_{i=1}^{N}{f}_{i}\,{\beta }_{i}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{\beta }_{i}=\displaystyle \frac{3\,{\epsilon }_{m}}{{\epsilon }_{i}+2\,{\epsilon }_{m}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&f=\displaystyle \sum _{i=1}^{N}{f}_{i}.\end{eqnarray} \tag{ 3 }$$

For the case of the Bruggeman rule, the mixed material itself acts as matrix and the mixed dielectric function can be calculated by solving for $\bar{\epsilon }$ in the relation

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{N}{f}_{i}\,\displaystyle \frac{{\epsilon }_{i}-\bar{\epsilon }}{{\epsilon }_{i}+2\,\bar{\epsilon }}=0.\end{eqnarray} \tag{ 4 }$$

Both the Bruggeman and the Maxwell–Garnett rules are implemented in our opacity module. However, for the compact mixture of materials specified above and in Table 1, the Bruggeman rule is the appropriate choice. For this dust model, the resulting effective medium optical constants are shown in Figure 2.

Zoom In Zoom Out Reset image size

Our opacity module includes a subroutine to perform Mie opacity calculations using a Fortran90 subroutine for performance. This Fortran90 code is based on a Fortran70 version by Bruce Draine,¹² which itself is derived from the original Mie code published by Bohren & Huffman (1998). To avoid strong and artificial Mie interferences, we do not use single-grain-size opacities, but instead calculate the opacity for 40 linearly spaced bins within each grain size bin and average over those opacity values to calculate an averaged opacity value for every size bin (each bin is 0.035 dex in size). The resulting absorption and scattering opacities are shown in Figure 3 and are available from the module repository.¹³

Zoom In Zoom Out Reset image size

It is known that dust grains in protoplanetary disks are not "mono-disperse;" i.e., at a given radius in the disk the dust does not consist of only a single size, or a narrow size distribution. Perhaps the most spectacular evidence of this is found in the source IM Lup. When observed at near-infrared wavelengths, this disk shows a strongly flaring geometry (Avenhaus et al. 2018). Clearly, a substantial amount of fine-grained dust is suspended several scale heights above the midplane, and is continuously replenished by turbulent stirring. When observed at millimeter wavelengths, however, we see a disk with a small-scale ring and spiral substructure that can only be explained if the geometry of this dust layer is vertically geometrically thin (e.g., Pinte et al. 2008; Huang et al. 2018). This must be of a grain population that is vastly larger (and/or more compact) than the grains seen in the near-infrared. Posed more precisely: IM Lup features at least two dust grain populations, one with a very small Stokes number, and therefore vertically extended, and one with a much larger Stokes number, and therefore vertically flat due to settling.

It is reasonable to expect that the dust population in fact consists of a continuous size distribution instead of just two distinct sizes. This is what is expected from models of dust coagulation that include fragmentation (Weidenschilling 1984; Dullemond & Dominik 2005; Brauer et al. 2008; Birnstiel et al. 2010). These populations change with time as the grains drift and grow at different rates. The complexity of this process makes it hard to define a simple "one-size-fits-all" dust opacity model to be used for interpreting millimeter continuum maps of protoplanetary disks. On the other hand, detailed coagulation/fragmentation modeling is numerically expensive, and it is not feasible to analyze all data with such full-fledged models.

Many authors, therefore, compromise by assuming that the dust grain size distribution follows a simple power law with a cutoff at small and large grain sizes,

$$\begin{eqnarray}n(a)\propto \,\left\{\begin{array}{ll}{a}^{-q} & \mathrm{for}\,{a}_{\min }\leqslant a\leqslant {a}_{\max }\\ 0 & \mathrm{else},\end{array}\right.\end{eqnarray} \tag{ 5 }$$

where the total dust density is defined as ${\rho }_{{\rm{d}}}={\int }_{0}^{\infty }\,n(a)m(a){da}$ with m(a) being the mass of a dust particle of radius a. The resulting opacity at (sub-)millimeter wavelengths is found to be less sensitive to the minimum grain radius ${a}_{\min }$, but much more so to the maximum particle size ${a}_{\max }$ as well as the power-law index q (Draine 2006; Ricci et al. 2010).

The index q was found to be around 3.5 for interstellar extinction measurements (e.g., Mathis et al. 1977, henceforth MRN), which is consistent with collisional cascades (Dohnanyi 1969; Tanaka et al. 1996) and also found to be consistent with submillimeter observations of debris disks (Ricci et al. 2015). However, the physics of debris disks is very different from gaseous protoplanetary disks. One way out is to use simplified dust coagulation/fragmentation models. For instance, B11 presented an analytic multi-power-law fit to the results of the full-fledged numerical dust coagulation/fragmentation models. A summary of the size distributions and the acronyms used throughout this Letter can be found in Table 2.

In the following, we will use the DSHARP opacity model of Section 2 and apply it to a simple power-law size distribution. We compare those results to the ones obtained if the analytic coagulation model fits of B11 are used. The Python script for creating the resulting size-averaged opacities is publicly available in the module repository.

The total absorption opacity ${\kappa }_{\nu }^{\mathrm{abs},\mathrm{tot}}$ of a particle size distribution n(a) at frequency ν is calculated from the size-dependent opacity ${\kappa }_{\nu }^{\mathrm{abs}}(a)$ via

$$\begin{eqnarray}&&{\kappa }_{\nu }^{\mathrm{abs},\mathrm{tot}}=\displaystyle \frac{{\int }_{{a}_{\min }}^{{a}_{\max }}n(a)\,m(a)\,{\kappa }_{\nu }^{\mathrm{abs}}(a)\,{da}}{{\int }_{{a}_{\min }}^{{a}_{\max }}n(a)\,m(a)\,{da}}.\end{eqnarray} \tag{ 6 }$$

The top panel in Figure 4 shows the total absorption and scattering opacities at a wavelength of 1 mm for a particle size distribution with ${a}_{\min }={10}^{-5}\,\mathrm{cm}$ as function of ${a}_{\max }$. The bottom panel shows the spectral index $\beta =\partial \mathrm{ln}{\kappa }_{\nu }^{\mathrm{abs},\mathrm{tot}}/\partial \mathrm{ln}\nu$. Similar trends as in Ricci et al. (2010) are observed: changes in the size distribution index q mainly affect the asymptotic behavior at long wavelengths and the strength of the Mie interference at ${a}_{\max }\sim \tfrac{\lambda }{2\pi }$. Figure 4 also shows that for size distributions that extend up to ${a}_{\max }\gtrsim 100\,\mu {\rm{m}}$, the scattering opacity ${\kappa }_{\nu }^{\mathrm{sca},\mathrm{tot}}$ exceeds the absorption opacity ${\kappa }_{\nu }^{\mathrm{abs},\mathrm{tot}}$.

Zoom In Zoom Out Reset image size

Three different particle size distributions that have the same maximum particle size ${a}_{\max }$ of 1 mm were chosen; one follows the MRN-like size-exponent of q = 3.5, while the other two distribution are steady-state size distributions where continuous particle growth and fragmentation lead to a stationary size distribution. The first of these (orange line in Figure 5) is from detailed analytical fits to numerical simulations from B11. The second steady-state distribution, shown in green in Figure 5, is a simplified version of these fits. This implements only the piecewise power laws from B11 and ignores finer details. This avoids calculating collision velocities for all particle sizes and thus makes the calculation easier and faster (see Appendix A). This simplified fit still captures the important aspects of the simulated distributions much better than the two-power-law fits used in Birnstiel et al. (2015) and Ormel & Okuzumi (2013). Especially for large particle sizes, the two-power-law fits can underpredict the number of small particles available.

Zoom In Zoom Out Reset image size

Wavelength-dependent opacities for all three distributions have been calculated and are shown in Figure 6 in comparison with opacities used in the literature. It can be seen that the overall behavior is—by construction—similar to the opacities in D'Alessio et al. (2001) or Andrews et al. (2009), with slightly different behavior at long wavelengths. Differences in the μm wavelength range are mainly due to the different amounts of small grains present in the distributions due to the knee in the steady-state distributions. This comes from the fact that smaller particles that move at higher Brownian motion velocities are more efficiently incorporated into larger particles. It can be seen that the two different fitting methods (labeled B11 and B11S) yield virtually identical opacities. The small differences to the simple MRN power law stems from the fact that parameters were chosen to yield the same ${a}_{\max }$. As seen from Figure 4, ${a}_{\max }$ is the most important parameter influencing the size-averaged opacity.

Zoom In Zoom Out Reset image size

The simplified fits labeled B11S are compared to the more detailed fits from B11 in Figure 7. Given the uncertainty in the details of collision models (see, for example, Güttler et al. 2010; Windmark et al. 2012), and for ease of reproduction, we will be using the B11S fits in the following. They are explained in Appendix A. Throughout this Letter, we assume a dust-to-gas mass ratio of 0.01 and consider size distributions integrated over height; settling will cause the size distribution to depend on the vertical position above the midplane.

Zoom In Zoom Out Reset image size

Figure 7 shows how the particle size distribution in steady state varies with the gas temperature T, the gas surface density ${{\rm{\Sigma }}}_{{\rm{g}}}$, the fragmentation threshold velocity v_frag, and the turbulence parameter α (Shakura & Sunyaev 1973). It can be seen that the position in the knee of the distribution at size a_BT (compare Equation (37) in B11) has only a weak dependence on those parameters, but that the maximum particle size ${a}_{\max }$ is a strong function of these parameters (quadratic in v_f, linear in all others). As such, it also affects the Planck and Rosseland mean opacities

$$\begin{eqnarray}&&{\bar{\kappa }}_{{\rm{P}}}(T)=\displaystyle \frac{{\int }_{0}^{\infty }{\kappa }_{\nu }^{\mathrm{abs},\mathrm{tot}}{B}_{\nu }(T)d\nu }{{\int }_{0}^{\infty }{B}_{\nu }(T)d\nu },\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\bar{\kappa }}_{{\rm{R}}}(T)={\left(\displaystyle \frac{{\int }_{0}^{\infty }\tfrac{1}{{\kappa }_{\nu }^{\mathrm{ext},\mathrm{tot}}}\tfrac{{{dB}}_{\nu }(T)}{{dT}}d\nu }{{\int }_{0}^{\infty }\tfrac{{{dB}}_{\nu }(T)}{{dT}}d\nu }\right)}^{-1},\end{eqnarray} \tag{ 8 }$$

because now not only the Planck spectrum ${B}_{\nu }(T)$ is temperature dependent, but also the size-averaged opacities ${\kappa }_{\nu }^{\mathrm{abs},\mathrm{tot}}$, and ${\kappa }_{\nu }^{\mathrm{ext},\mathrm{tot}}$. This means that the mean opacities are additionally dependent on other physical parameters, ${{\rm{\Sigma }}}_{{\rm{g}}},{v}_{\mathrm{frag}}$, and α.

As an example, Figure 8 shows the Rosseland and Planck mean opacities for two cases: a MRN-size distribution (as in Figure 5) and the steady-state distributions B11 and B11S as a function of temperature. It can also be seen that differences between the two steady-state distributions are small, allowing the simpler model to be used without caveats. For the fiducial values of ${M}_{\star }={M}_{\odot }$, $r=1\,\mathrm{au}$, ${v}_{\mathrm{frag}}=100\,\mathrm{cm}\,{{\rm{s}}}^{-1}$, and α = 10⁻³, high-temperature mean opacities are generally higher for steady-state distributions owing to the fact that the knee at μm sizes produces more small grains for the same ${a}_{\max }$ than a single power-law size distribution. Furthermore, high temperatures tend to produce smaller ${a}_{\max }$ than 1 mm, which additionally increases the amount of small particles that contribute most to the mean opacities. The shaded areas in Figure 8 cover the ranges $1\,{\rm{K}}\leqslant T\leqslant 1500\,{\rm{K}}$, $50\,\mathrm{cm}\,{{\rm{s}}}^{-1}\,\leqslant {v}_{\mathrm{frag}}\leqslant 3\times {10}^{3}\mathrm{cm}\,{{\rm{s}}}^{-1}$, 10⁻⁵ ≤ α ≤ 10⁻¹, 1 g cm⁻² ≤ Σ_g ≤ 10⁻⁴ g cm⁻². It should be noted that the distributions discussed here only apply to those parts of the disk where particles reach the fragmentation barrier a_frag—this is only possible if (1) collision velocities are high enough (i.e., the root in Equation (28) is real), and (2) fragmentation is more important in limiting particle growth than drift, which is the case if (see Birnstiel et al. 2015)

$$\begin{eqnarray}&&\alpha \gt \displaystyle \frac{\left|\gamma \right|\,{{\rm{\Sigma }}}_{{\rm{g}}}}{3\,{{\rm{\Sigma }}}_{{\rm{d}}}}\left(\displaystyle \frac{{v}_{\mathrm{frag}}}{{V}_{{\rm{K}}}}\right),\end{eqnarray} \tag{ 9 }$$

where γ is the radial logarithmic pressure gradient $\partial \mathrm{ln}P/\partial \mathrm{ln}r$ and V_K the Keplerian orbital velocity. If the drift limit applies, the size distribution will contain more mass at the largest sizes and is more strongly dependent on global redistribution of particles. These non-local processes can be approximated as in Birnstiel et al. (2015); however, the accuracy of these approximations are likely only good enough for applying them to the long-wavelength opacity. Short wavelengths are too sensitive to small changes in the amount of small grains, which in turn depend sensitively on radial mixing and details of the collisional model.

Zoom In Zoom Out Reset image size

For protoplanetary disks, it is mostly assumed that only the absorption opacity ${\kappa }_{\nu }^{\mathrm{abs}}$ matters, and not the scattering opacity. For optically thin dust layers this is indeed appropriate. Recently, the importance of scattering and its effects on (sub-)millimeter polarization of disks was pointed out by Kataoka et al. (2015). In the DSHARP campaign we have seen that the optical depths are not that low ($0.1\lesssim \tau \lesssim 0.6$; see Dullemond et al. 2018). Furthermore, the CO line extinction data of HD 163296 discussed by Isella et al. (2018) suggest that the dust layer has an extinction optical depth close to unity (see also Guzman et al. 2018). Even if the absorption optical depth is substantially below 1, the total extinction (absorption + scattering) can easily exceed unity if the grains are of similar size to the wavelength. In fact, for a ≃ λ/2π = 0.13 cm/2π = 0.02 cm the albedo of the grain can be as high as 0.9 (see Figure 4 and Appendix B).

The inclusion of scattering complicates the radiative transfer equation enormously. Strictly speaking, a full radiative transfer calculation, for instance with a Monte Carlo code, is necessary. However, in the spirit of this Letter we wish to find a simple approximation to handle this without resorting to complex numerical simulations.

There are two issues to be solved. What thermal emission will a non-optically thin dust layer produce if scattering is taken into account? Also, how do we compute the extinction coefficient to be used in the CO line extinction analysis?

For the first issue, we will outline here a simple two-stream radiative transfer approach to the problem. We will assume that the dust seen in the ALMA observations is located in a geometrically thin layer at the midplane, so we can use the 1D slab geometry approach. We will assume that the scattering is isotropic. This may be a bad approximation, especially for $2\pi a\gg \lambda$. To reduce the impact of this approximation we replace the scattering opacity ${\kappa }_{\nu }^{\mathrm{sca}}$ with

$$\begin{eqnarray}&&{\kappa }_{\nu }^{\mathrm{sca},\mathrm{eff}}=(1-{g}_{\nu }){\kappa }_{\nu }^{\mathrm{sca}}\end{eqnarray} \tag{ 10 }$$

where g_ν is the usual forward-scattering parameter (the expectation value of $\cos \theta$, where θ is the scattering angle). According to Ishimaru (1978), this approximation works well in optically thick media.

We will now follow the two-stream/moment method approach from Rybicki & Lightman (1991) to derive the solution to the emission/absorption/scattering problem in this slab. The slab is put between $z=-\tfrac{1}{2}{\rm{\Delta }}z$ and $z=+\tfrac{1}{2}{\rm{\Delta }}z$ and we assume a constant density of dust between these two boundaries. The mean intensity ${J}_{\nu }(z)$ of the radiation field then obeys the equation

$$\begin{eqnarray}&&\displaystyle \frac{1}{3}\displaystyle \frac{{d}^{2}{J}_{\nu }}{d{\tau }_{\nu }^{2}}={\epsilon }_{\nu }\left({J}_{\nu }-{B}_{\nu }({T}_{d})\right)\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{\tau }_{\nu }={\rho }_{d}({\kappa }_{\nu }^{\mathrm{abs}}+{\kappa }_{\nu }^{\mathrm{sca},\mathrm{eff}})z\equiv {\rho }_{d}\,{\kappa }_{\nu }^{\mathrm{tot}}z\end{eqnarray} \tag{ 12 }$$

with ${\rho }_{d}$ being the dust density, and

$$\begin{eqnarray}&&{\epsilon }_{\nu }^{\mathrm{eff}}=\displaystyle \frac{{\kappa }_{\nu }^{\mathrm{abs}}}{{\kappa }_{\nu }^{\mathrm{abs}}+{\kappa }_{\nu }^{\mathrm{sca},\mathrm{eff}}}.\end{eqnarray} \tag{ 13 }$$

The boundary conditions at $z=\pm \tfrac{1}{2}{\rm{\Delta }}z$ are

$$\begin{eqnarray}&&\displaystyle \frac{{{dJ}}_{\nu }}{d{\tau }_{\nu }}=\mp \sqrt{3}{J}_{\nu }.\end{eqnarray} \tag{ 14 }$$

This leads to the following solution:

$$\begin{eqnarray}&&\displaystyle \frac{{J}_{\nu }({\tau }_{\nu })}{{B}_{\nu }({T}_{d})}=1-b\left({e}^{-\sqrt{3{\epsilon }_{\nu }^{\mathrm{eff}}}\left(\tfrac{1}{2}{\rm{\Delta }}\tau -{\tau }_{\nu }\right)}+{e}^{-\sqrt{3{\epsilon }_{\nu }^{\mathrm{eff}}}\left(\tfrac{1}{2}{\rm{\Delta }}\tau +{\tau }_{\nu }\right)}\right)\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Delta }}\tau ={\rho }_{d}\,{\kappa }_{\nu }^{\mathrm{tot}}{\rm{\Delta }}z$, and b is

$$\begin{eqnarray}&&b={\left[(1-\sqrt{{\epsilon }_{\nu }^{\mathrm{eff}}}){e}^{-\sqrt{3{\epsilon }_{\nu }^{\mathrm{eff}}}{\rm{\Delta }}\tau }+1+\sqrt{{\epsilon }_{\nu }^{\mathrm{eff}}}\right]}^{-1}.\end{eqnarray} \tag{ 16 }$$

Given this solution for the mean intensity ${J}_{\nu }({\tau }_{\nu })$ we can now numerically integrate the formal transfer equation along a single line of sight passing through the slab at an angle θ:

$$\begin{eqnarray}&&\mu \displaystyle \frac{{{dI}}_{\nu }({\tau }_{\nu })}{d{\tau }_{\nu }}={\epsilon }_{\nu }^{\mathrm{eff}}{B}_{\nu }({T}_{d})+(1-{\epsilon }_{\nu }^{\mathrm{eff}}){J}_{\nu }({\tau }_{\nu })-{I}_{\nu }({\tau }_{\nu })\end{eqnarray} \tag{ 17 }$$

where $\mu =\cos \theta$. We start at ${\tau }_{\nu }=-\tfrac{1}{2}{\rm{\Delta }}\tau$ with ${I}_{\nu }=0$, and integrate to ${\tau }_{\nu }=+\tfrac{1}{2}{\rm{\Delta }}\tau$. The resulting ${I}_{\nu }^{\mathrm{out}}={I}_{\nu }(\tfrac{1}{2}{\rm{\Delta }}\tau )$ is the intensity that is observed by the telescope. An approximation for ${I}_{\nu }^{\mathrm{out}}$, which works well to within a few percent, is the following modified version of the Eddington–Barbier approximation:

$$\begin{eqnarray}&&{I}_{\nu }^{\mathrm{out}}\simeq \left(1-{e}^{-{\rm{\Delta }}\tau /\mu }\right){S}_{\nu }\left(\left(\tfrac{1}{2}{\rm{\Delta }}\tau -{\tau }_{\nu }\right)/\mu =2/3\right)\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&{S}_{\nu }({\tau }_{\nu })={\epsilon }_{\nu }^{\mathrm{eff}}{B}_{\nu }({T}_{d})+(1-{\epsilon }_{\nu }^{\mathrm{eff}}){J}_{\nu }({\tau }_{\nu })\end{eqnarray} \tag{ 19 }$$

is the source function. In the optically thin case, when ${\rm{\Delta }}\tau /\mu \lt 2/3$, the value of S_ν is taken at the edge of the slab. The results are shown in Figure 9.

Zoom In Zoom Out Reset image size

For small optical depth (${\rm{\Delta }}\tau \ll 1$) the role of scattering vanishes, and the solution approaches: ${I}_{\nu }^{\mathrm{out}}\to {\epsilon }_{\nu }^{\mathrm{eff}}{\rm{\Delta }}\tau {B}_{\nu }({T}_{d})/\mu$. This is the same limiting solution as when ${\kappa }_{\nu }^{\mathrm{sca}}$ is set to zero but ${\kappa }_{\nu }^{\mathrm{abs}}$ is kept the same. For high optical depth the outcoming intensity does not saturate to the Planck function, but a bit below, if the albedo is non-zero. This is the well-known effect that scattering makes objects appear cooler than they really are.

Now let us turn to the second issue to be considered: how to calculate the actual extinction coefficient for the CO line extinction analysis. At first sight the answer is simple:

$$\begin{eqnarray}&&{\kappa }_{\nu }^{\mathrm{ext}}={\kappa }_{\nu }^{\mathrm{abs}}+{\kappa }_{\nu }^{\mathrm{sca}}\end{eqnarray} \tag{ 20 }$$

where we use ${\kappa }_{\nu }^{\mathrm{sca}}$, not ${\kappa }_{\nu }^{\mathrm{sca},\mathrm{eff}}$, for the scattering. This is, in fact, the correct answer for the case in which the emitting CO line emitting layer is far behind the extincting dust layer, where "far" is defined in comparison to the width of the extincting dust ring. However, if this condition is not met, then the scattering will not reduce the CO line intensity as much as naively expected. CO line photons that are heading elsewhere might then, in fact, get scattered into the line of sight. This effect is exacerbated for the case of small-angle scattering. One can therefore argue that this effect reduces the extinction of the CO line emission by the dust layer. For the case of HD 163296 (Isella et al. 2018) it appears that these effects are not too strong, so a first analysis without accounting for this is in order. A final answer, however, may require a full treatment of 3D radiative transfer.

In this Letter we present methods for translating observed dust emission and extinction used in the DSHARP campaign, which can also be used by other work. The DSHARP opacities presented here are merely a choice, based on reasonable assumptions. They provide the standard used within the DSHARP campaign.

We further explore how steady-state size distributions in a coagulation-fragmentation equilibrium affect dust opacities by introducing dependencies on temperature, surface density, turbulence, and material properties. We provide simplified fits to analytical functions and show the range of Rosseland and Planck mean opacities covered for a wide range of parameter choices.

Given the large albedo at (sub-)millimeter wavelength ranges, we derive solutions to the radiative transfer equation for a homogeneous medium with scattering and absorption. Together with the DSHARP dust model, and based on the measurements of Isella et al. (2018), we find that the particles in the rings of HD 163296 should be at least 0.2 cm in size.

Along with this Letter, we present publicly available Python scripts that contain the optical data of many literature materials. In addition to that, functions are available for mixing optical constants with effective medium theory, for calculating opacities using Mie theory, and for averaging opacities over particle size distributions. Implementations of the steady-state particle size distributions discussed in this Letter are also included. The online material includes these python modules, scripts for generating the results and figures of this Letter, and the opacity tables. These materials will likely be extended in the future but the version used in this Letter is available at Birnstiel (2018). Additional material will be described in appendices of future papers, as they become available.

We like to thank Ryo Tazaki, Akimasa Kataoka, and Satoshi Okuzumi for helpful discussions, and Thomas Henning for providing the optical constants data used in this Letter. T.B. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 714769. C.P.D. acknowledges support by the German Science Foundation (DFG) Research Unit FOR 2634, grants DU 414/22-1 and DU 414/23-1. Z.Z. and S.Z. acknowledge support from the National Aeronautics and Space Administration through the Astrophysics Theory Program with grant No. NNX17AK40G and Sloan Research Fellowship. Simulations are carried out with the support from the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG- AST130002. S.A. and J.H. acknowledge funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 $\_$2-0012 issued through the Exoplanets Research Program. J.H. acknowledges support from the National Science Foundation Graduate Research Fellowship under grant No. DGE-1144152. A.I. acknowledges support from the National Aeronautics and Space Administration under grant No. NNX15AB06G issued through the Origins of Solar Systems program, and from the National Science Foundation under grant No. AST-1715719. L.P. acknowledges support from CONICYT project Basal AFB-170002 and from FCFM/U. de Chile Fondo de Instalación Académica.

Software: Numpy (Van Der Walt et al. 2011), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration et al. 2013).

The simplified version of the B11 steady-state distributions used in this Letter are defined as a broken but continuous power law as a function of particle size

$$\begin{eqnarray}\sigma (a)=\left\{\begin{array}{ll}{a}^{p-1} & \mathrm{for}\quad a\leqslant {a}_{\mathrm{frag}}\\ 0 & \mathrm{else}\end{array}\right.\end{eqnarray} \tag{ 21 }$$

where the exponents p are changing at specific sizes. They are chosen according to the algorithm

$$\begin{eqnarray}\begin{array}{rcl} & & \mathrm{if}(a\lt {a}_{\mathrm{BT}}):\\ p & = & \displaystyle \frac{3}{2}\quad \mathrm{or}\quad \displaystyle \frac{5}{4}\\ & & \mathrm{elseif}(a\leqslant {a}_{12}):\\ p & = & \displaystyle \frac{1}{4}\quad \mathrm{or}\quad 0\\ & & \mathrm{elseif}(a\leqslant {a}_{\mathrm{frag}}):\\ p & = & \displaystyle \frac{1}{2}\quad \mathrm{or}\quad \displaystyle \frac{1}{4}.\end{array}\end{eqnarray} \tag{ 22 }$$

Here the first value of p corresponds to the case $a\leqslant {a}_{\mathrm{set}}$, and the second (after the "or") applies to sizes above a_set. The sizes ${a}_{\mathrm{BT}},{a}_{12},{a}_{\mathrm{set}}$ are calculated according to Equations (37), (40), and (27) in B11,

$$\begin{eqnarray}&&{a}_{\mathrm{set}}=\displaystyle \frac{2\,\alpha \,{{\rm{\Sigma }}}_{{\rm{g}}}}{\pi \,{\rho }_{{\rm{s}}}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{a}_{\mathrm{BT}}={\left[\displaystyle \frac{8{{\rm{\Sigma }}}_{{\rm{g}}}}{\pi {\rho }_{{\rm{s}}}}\cdot {\mathrm{Re}}^{-\tfrac{1}{4}}\cdot \sqrt{\displaystyle \frac{\mu {m}_{{\rm{p}}}}{3\pi \alpha }}\cdot {\left(\displaystyle \frac{4\pi }{3}{\rho }_{{\rm{s}}}\right)}^{-\tfrac{1}{2}}\right]}^{\tfrac{2}{5}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{a}_{12}=\displaystyle \frac{1}{{y}_{a}}\,\displaystyle \frac{2{{\rm{\Sigma }}}_{{\rm{g}}}}{\pi \,{\rho }_{{\rm{s}}}}\cdot {\mathrm{Re}}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 25 }$$

and the particle Reynolds number $\mathrm{Re}$ is

$$\begin{eqnarray}&&\mathrm{Re}\approx \displaystyle \frac{\alpha \,{{\rm{\Sigma }}}_{{\rm{g}}}\,{\sigma }_{{H}_{2}}}{2\,\mu \,{m}_{{\rm{p}}}}.\end{eqnarray} \tag{ 26 }$$

Here, m_p is the proton mass μ = 2.3 the mean molecular mass in atomic units, ${y}_{a}\simeq 1.6$ (Ormel & Cuzzi 2007), ${\sigma }_{{{\rm{H}}}_{2}}\,\simeq 2\times {10}^{-15}\,{\mathrm{cm}}^{2}$ the atomic hydrogen cross section. The fragmentation limit is given by

$$\begin{eqnarray}&&{a}_{\mathrm{frag}}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{\pi \,{\rho }_{{\rm{s}}}\,b}\,\sqrt{1-4\,{b}^{2}},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&b=\displaystyle \frac{1}{3\,\alpha }{\left(\displaystyle \frac{{v}_{\mathrm{frag}}}{{c}_{{\rm{s}}}}\right)}^{2}.\end{eqnarray} \tag{ 28 }$$

If ${a}_{\mathrm{frag}}\lt {a}_{12}$, the fragmentation limit in the first turbulent regime needs to be calculated from

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{frag}}={\mathrm{Re}}^{-1/4}\,\displaystyle \frac{{v}_{\mathrm{frag}}}{{c}_{{\rm{s}}}}\sqrt{\displaystyle \frac{2}{3\,\alpha }}.\end{eqnarray} \tag{ 29 }$$

The distribution $\sigma (a)$ is normalized to the total dust surface density

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}}}={\int }_{-\infty }^{\infty }\sigma (a)\mathrm{ln}a.\end{eqnarray} \tag{ 30 }$$

Under the assumption of vertically well-mixed dust (which is not applicable in most parts of the disk), $\sigma (a)$ and n(a) are directly proportional to each other for all particle sizes. Vertical settling will reduce the vertical scale height for larger particles. The local densities of each particle size can be calculated under the assumption of a settling-mixing equilibrium, as for example in Fromang & Nelson (2009). An numerical implementation is included in the python module.

As explained in Section 2, the DSHARP opacities are based on several approximations or assumptions, based on practical choices. As such, they are meant to be used as a reference choice along the lines of previous literature values, and not to be seen as the last word on the subject. To demonstrate how some of these choices affect the resulting opacity values, we will explore the effects of mixing rules/porosity, water abundance, and the choice of carbonaceous material. However, these are by far not the only uncertainties. Far-infrared or (sub-)mm opacities were also found to be affected by temperature dependencies (Boudet et al. 2005; Coupeaud et al. 2011; Demyk et al. 2017b, 2017a). Instead of being compact and porous, particles could also be fractal instead (Tazaki et al. 2016; Tazaki & Tanaka 2018), and the composition and shape of the particles are largely unknown. Exploring all of these possible influences, however, is beyond the scope of this Letter, and we instead refer the reader to dedicated studies of this subject (for example, Draine 2006; Kataoka et al. 2014, 2015; Min et al. 2016; Tazaki et al. 2016; Woitke et al. 2016; Tazaki & Tanaka 2018).

In the following, we will start with the DSHARP opacities (labeled as default in Figures 10 and 11), as explained in Section 2 and then change some of those assumptions individually: using the same relative volume fractions, we include 80% porosity. In this case, the Maxwell–Garnett mixing rule (Equation (1)) is used and the resulting optical properties are shown in Figures 10 and 11, labeled as porous. It can be seen that in the porous case, the millimeter-range opacities are much lower, of the order of 0.3 cm² g⁻¹, the highest scattering opacity is shifted to larger ${a}_{\max }$, and the reduced Mie-interferences also result in a flat spectral index profile, as discussed in Kataoka et al. (2015). It should be noted that the absorption opacity at millimeter wavelengths can be enhanced by a factor of about 2 for silicate particles and a factor of 4 for amorphous carbon due to the interaction of monomers, which is ignored in the Maxwell–Garnett Mie theory, as shown in Tazaki & Tanaka (2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the next example, it is not the vacuum volume fraction (i.e., the porosity) that is increased, but rather the water volume fraction is raised to 60% (labeled high-water). This results in very small changes in the optical properties of the distribution, including slightly increasing the water feature around 3 μm.

Very significant changes are found, if the material termed "organics" is exchanged for other carbonaceous materials: as an example, we used the carbonaceous analog pyrolized at T = 800 K of Jäger et al. (1998; labeled Jäger), the graphite optical constants from Draine (2003; sample a = 0.01 μm, perpendicular alignment, labeled Draine), and the cosmic carbon analogs (sample ACH) from Zubko et al. (1996; labeled as Zubko), that are also widely used in the literature, for example in Ricci et al. (2010). Figures 10 and 11 show that those carbonaceous materials cause the strongest variations. They tend to give higher absorption opacities (see also Min et al. 2016), their maximum absorption can be significantly shifted away from ${a}_{\max }\sim \lambda /(2\,\pi )$ and the resulting spectral index is consequently affected significantly as well. Especially at millimeter wavelength these different compounds affect the absorption opacity the most (see Figure 11).

The bottom panel of Figure 10 shows the absorption probability ${\epsilon }_{1\ \mathrm{mm}}^{\mathrm{eff}}$ (which is $1-{\eta }_{\nu }$, where ${\eta }_{\nu }$ is the single scattering albedo) at a wavelength of λ = 1 mm averaged over the MRN-like size distribution with varying ${a}_{\max }$. It can be seen that, despite the strong changes in the opacities or the spectral index, this quantity has a very similar behavior in all cases: it is close to unity for particles smaller than λ/(2 π) and then drops quite sharply for sizes larger than that. Only the value reached for large ${a}_{\max }$ is very sensitive to the choices of the opacity model. Similar to the polarization fraction of scattered thermal dust emission discussed in Kataoka et al. (2015), ${\epsilon }_{\mathrm{eff}}$ can thus be used to constrain the maximum particle size. For the absorption and extinction optical depths measured in Isella et al. (2018), these considerations indicate that the particles present should be at least of a size of ∼0.2 mm.