UTILITARIAN OPACITY MODEL FOR AGGREGATE PARTICLES IN PROTOPLANETARY NEBULAE AND EXOPLANET ATMOSPHERES

Jeffrey N. Cuzzi; Paul R. Estrada; Sanford S. Davis

doi:10.1088/0067-0049/210/2/21

1. INTRODUCTION

The total extinction opacity of a medium κ_e (cm² g⁻¹) is the ratio of its volume extinction coefficient (cm⁻¹ of path) to its volume mass density. Equivalently, it is the cross-section per unit mass along a path. In protoplanetary nebulae, the relatively small fraction (by mass) of solid particles provides the bulk of the opacity (D'Alessio et al. 1999, 2001) unless the temperature is so large (above 1500 K) that common geological solids evaporate and molecular species dominate (Nakamoto & Nakagawa 1994; Ferguson et al. 2005). In exoplanet atmospheres, which are generally much denser, the particle contributions are more situation-dependent (Marley et al. 1999; Sudarsky et al. 2000, 2003; Ackerman & Marley 2001; Tsuji 2002; Currie et al. 2011; de Kok et al. 2011; Madhusudhan et al. 2011; Morley et al. 2012, 2013; Vasquez et al. 2013).

While typical interstellar grain size distributions (Mathis et al. 1977; Draine & Lee 1984) have been assumed in many nebula and planetary atmospheric opacity models (Pollack et al. 1985, 1994; D'Alessio et al. 1999, 2001; Bodenheimer et al. 2000; Dullemond et al. 2002), coagulation models in the nebula context, going back to Weidenschilling (1988, 1997), and most recently Ormel & Okuzumi (2013), suggest that such small grains grow to 100 μm size or larger on short timescales. This is because small and/or fluffy grains collide at very low relative velocity and stick readily. Other recent studies of grain growth in the nebula context include Brauer et al. (2008), Zsom et al. (2010), Schraepler et al. (2012), and Birnstiel et al. (2010, 2012) to give only a few examples. In cloud formation models for gas giant planets (Movshovitz et al. 2010) and for brown dwarfs and exoplanets (Helling et al. 2001; Kaltenegger et al. 2007; Kitzmann et al. 2010; Zsom et al. 2012). A review by Helling et al. (2008) shows how exoplanet and brown dwarf clouds appear at different altitudes and with different particle sizes, formed from a variety of constituents from volatile ices to silicate or even iron metal. In a recent review, Marley et al. (2013) illustrate how some materials condense as liquids (forming droplets that are appropriate for the usual Mie theory models), while other materials condense as solids, forming irregular, porous aggregates as they grow. To our knowledge, no current exoplanet radiative transfer models have explored the implications of porous aggregate cloud particles. In this paper we present a simple way of accounting for radiative transfer involving aggregates of arbitrary size and porosity. Some giant planet atmosphere models (Podolak 2003; Movshovitz & Podolak 2008; Movshovitz et al. 2010; Helled & Bodenheimer 2011) and one other application (Amit & Podolak 2009) have already incorporated the opacity model we describe here. Using these opacity models, Movshovitz et al. (2010) showed that realistic grain coagulation leads to smaller opacities, and shorter formation times, for gas giant planets than previously appreciated in the context of the core-instability scenario (see Section 5.4).

Miyake & Nakagawa (1993), Pollack et al. (1994), and D'Alessio et al. (2001) give a few examples of the dramatic effect of particle growth on opacity across a broad range of sizes of interest in the nebula case. Because the aggregating constituents are not always liquid, the growing particles can have a fairly open, highly porous, nature (Meakin & Donn 1988; Donn 1990; Dominik & Tielens 1997; Beckwith et al. 2000; Dominik et al. 2007; Blum 2010). Growth beyond mm–cm sizes in protoplanetary nebulae remains an active area of study (see reviews by Cuzzi & Weidenschilling 2006; Dominik et al. 2007, and models by Brauer et al. 2008; Zsom et al. 2010; Hughes & Armitage 2012; Birnstiel et al. 2010, 2012; and others), yet none of these latter models actually incorporates equally realistic and self-consistent treatments of opacity either in their radiative transfer or in predictions of spectral energy distributions (SEDs), and only Pollack et al. (1985, 1994) and Miyake & Nakagawa (1993) have addressed the issue of large particle porosity in the implications for thermal radiative transfer. In both nebula and atmospheric cases, the vertical temperature structure is determined by the heating source (external and/or internal) and the opacity—usually expressed as a temperature-weighted mean such as the Rosseland mean opacity κ_R (Pollack et al. 1985, 1994; Henning & Stognienko 1996; Semenov et al. 2003); the Planck mean can also be used for optically thin applications and is even simpler to calculate (Nakamoto & Nakagawa 1994).

Interpretation of SEDs, or intensity as a function of wavelength, requires an understanding of the wavelength-dependent opacity as well as the detailed vertical temperature distribution of the nebula or planet's atmosphere. Traditionally, regions sampled are thought to be optically thin, and particles poor scatterers being much smaller than the wavelength, so only the absorption part of the extinction opacity is needed (Miyake & Nakagawa 1993). Nebula gas masses have often been inferred by combining the mass of particulates (inferred from the observed emission and theoretical opacity) with the canonical mass ratio of roughly one percent for solids to total mass (e.g., Andrews & Williams 2007; Williams & Cieza 2011). Theoretical analyses of observed mm–cm wavelength opacities at least since D'Alessio et al. (2001) (e.g., review by Natta et al. 2007; also Isella et al. 2009; Birnstiel et al. 2010; Ricci et al. 2010, and others) use full Mie theoretical analyses of particle opacities, which we show in Section 5.3 is, in fact, essential when the particle size and wavelength are comparable. de Kok et al. (2011) modeled emission spectra of exoplanets, where the primary spectral variation is that of the gas, but where the more gradually varying extinction properties of cloud particles still provide important variations across the near-IR passband because the thermal wavelengths are not extremely large compared to the particle size.

Below we present a very simple, but general and flexible, model which captures the essential physics of opacity in ensembles of aggregate particles of arbitrary size and porosity, and can be easily made part of nebula or planetary atmosphere thermal structure and/or evolution models. We give some examples of how growth and porosity affect the opacity of realistic aggregate particles. The opacity model is simple to combine with any grain coagulation/aggregation model, and can be used to study the temporal/spatial dependence of the SED and Rosseland mean opacity κ_R on varying particle size, density, porosity, and composition, or to calculate thermal equilibrium profiles. The model is not ideally suited to interpreting mm–cm wavelength spectral slopes in protoplanetary disks, where the dominant size is often found to be comparable to the wavelength.

2. PRIOR WORK AND OVERVIEW OF PAPER

Draine & Lee (1984) studied small, independent particles (r ≪ λ, where r is particle radius and λ is wavelength) in the dipole approximation; Pollack et al. (1985) modeled various size distributions of mineralogically realistic, spherical, solid grains using exhaustively compiled refractive indices and Mie scattering; Wright (1987) modeled grains of various fractal dimension (basically, porosity) using the discrete dipole approximation or DDA (Purcell & Pennypacker 1973; Draine & Goodman 1993). Wright (1987) showed that nonsphericity could cause a significant enhancement of opacity, if grains had fairly large refractive indices. More recent work has elaborated upon this (Fabian et al. 2001; Min et al. 2008; Lindsay et al. 2013), showing how this caveat applies to detailed spectral analysis of highly nonspherical particles (elongated, flattened, or even having sharp edges) through strong absorption bands. Bohren & Huffman (1983, p. 140) also give constraints on how particle size and refractive index must both be considered, in cases where the refractive index is large. Our model is most suitable for roughly equidimensional particles with moderate refractive indices (high porosity and realistic compositional mixtures usually preclude high average indices).

Ossenkopf (1991) studied the effect of metallic inclusions on equidimensional porous particles, emphasizing small particles and the role of randomly connected iron grains in leading to effective dipole-like structures. For the low volume densities of iron grains in protoplanetary nebula aggregates, at least for the abundances assumed here, the effect is small. Pollack et al. (1985, 1994) assumed materials which are cosmically abundant (with an eye towards molecular cloud and protoplanetary nebula applications); of these, water ice, silicates, and refractory organics have refractive indices that are too small for nonsphericity effects to be especially important in the regime where particle size is smaller than a wavelength. Other common materials, specifically iron metal and iron sulfide (Troilite, FeS) do have refractive indices that are high enough to show an effect, but nonsphericity was not treated by Pollack et al. (1985, 1994).

Regarding distributions of larger particles, Pollack et al. (1985, 1994) and Miyake & Nakagawa (1993) calculated both monochromatic opacities and Rosseland Mean opacities, and D'Alessio et al. (2001) calculated monochromatic opacities and employed them in a wide range of nebula structure models. Pollack et al. (1985) showed selected results in the separate limits of large and small particles, but did not connect them. Mizuno et al. (1988) used a combination of Rayleigh scattering for small particles and diffraction-augmented geometric optics (Q_e = 2) for large ones. Miyake & Nakagawa (1993) used a combination of full Mie calculations and geometrical optics, along with an ad hoc powerlaw particle size distribution, to demonstrate the effects of particle growth. To match confidently with the geometric optics limit, they carried out Mie calculations to impressively large values of r/λ ∼ 10⁵, with implications for array size and computational time that would be hard for the typical user to achieve even today, in multidimensional and/or evolutionary nebula or exoplanet atmosphere models. D'Alessio et al. (2001) also seem to have used full Mie calculations for all their particles, up to 10 cm radius, treating each material as a separate species. Pollack et al. (1985), Mizuno et al. (1988), Miyake & Nakagawa (1993), and D'Alessio et al. (2001) all assumed solid objects; Pollack et al. (1994) looked briefly at particles of higher porosity, using an Effective Medium Theory (EMT; see Section 4) and Mie calculations. Rannou et al. (1999) used a semi-empirical model of their own design to capture the properties of aggregates in certain size and compositional regimes.

In this paper, we build on and generalize this prior work by showing how arbitrary size, porosity, and compositional distributions can be handled easily by users wanting to explore their own choices for grain properties, and simply enough to incorporate within evolutionary models or in iterative analysis of observed SEDs for applications where grain growth into sizeable aggregates has occurred. We show how porosity trades off with size in determining emergent SEDs and determining Rosseland mean opacities. We use accurate material properties for a realistic, temperature-dependent, compositional suite, and simplified but realistic scattering and absorption efficiencies which avoid the need for numerical Mie calculations. In this sense the approach is similar in spirit to some previous studies, but emphasizes a utilitarian approach and ease of general applicability. The physics and simplifications are described in Section 3 below; in Section 4 we show some validation tests of the model, and in Section 5 we describe the behavior of porous aggregates. Some basic derivations are presented in Appendices A–C. The code itself is quite simple, taking negligible CPU time compared to a Mie code, and a version is available from the authors.

3. OPACITY MODEL

In this section we describe the theoretical basis for our calculations of particulate opacity. Our goal is to capture all of the significant physics in the simplest fashion possible, so the model can be incorporated into evolutionary models at little computational cost. Realistic material refractive indices are included for a cosmic abundance suite of likely nebula solids: water ice, silicates, refractory organics, iron sulfide, and metallic iron. These refractive indices, and relative abundances as a function of temperature are taken from Pollack et al. (1994; see Table 1). We chose these specific values to allow better validation and comparison with previous work, but they are widely used; alternate tabulations are easily incorporated, such as found in Draine & Lee (1984) or Henning et al. (1999).

Table 1. Our Assumed Compositional Mixture, Adapted from Pollack et al. (1994)

Material	Density	α_j	T_evap	β_j < 160 K	β_j < 425 K	β_j < 680 K	β_j < 1500 K
Material	Density	α_j	(K)	β_j < 160 K	β_j < 425 K	β_j < 680 K	β_j < 1500 K
Water Ice	0.9	5.55e−3	160	6.11e−1	0	0	0
Organics	1.5	4.13e−3	425	2.73e−1	7.00e−1	0	0
Troilite	4.8	7.68e−4	680	1.58e−2	4.07e−2	1.36e−1	0
Silicates	3.4	3.35e−3	1500	9.93e−2	2.55e−1	8.51e−1	9.84e−1
Iron (<680 K)	7.8	1.26e−4	1500	1.60e−3	4.11e−3	1.37e−2	⋅⋅⋅
Iron (>680 K)	7.8	6.15e−4	1500	⋅⋅⋅	⋅⋅⋅	0	1.58e−2

Notes. For simplicity we have merged their two kinds of silicates and their two kinds of organics. When Troilite (FeS) decomposes at 680 K, we follow Pollack et al. in assuming the liberated iron adds to the existing iron metal and the S remains in the gas phase. The parameters α_j and β_j define the fractional mass of a compositional species, and the fractional number of its particles if particles are segregated compositionally. For their definitions see Section 4 and Equations (C17) (or (19)) respectively.

Download table as: ASCII Typeset image

Our emphasis is on how size and porosity effects can be dealt with simply and seamlessly as grains aggregate and grow into a size regime where details of composition and refractive indices are perhaps secondary. Our basic approach will be to separate particles having an arbitrary size distribution into two regimes which both have closed-form solutions in practical cases. We use the fact that particle interactions with radiation of any wavelength λ can be systematically addressed in terms of the optical size of the particle x = 2πr/λ, where r is the particle radius. The model presented here treats "small" particles (x ≪ 1) as volume absorbers/scatterers and "large" particles (x ≫ 1) as geometrical optics absorbers/scatterers. We will occasionally refer to Van de Hulst (1957, henceforth VDH), Hansen & Travis (1974, henceforth HT), and Bohren & Huffman (1983).

3.1. Extinction, Absorption, and Scattering Efficiencies

Radiation is removed from any beam at a rate:

$\begin{equation} I = I_0 e^{- \kappa _e \rho l} \end{equation} \tag{ 1 }$

where ρ is the volume mass density of the gas and dust mixture, κ_e (cm² g⁻¹) is the total opacity, and l is the path length. The monochromatic opacity, for a gas-particle mixture containing particles of radius r with number density n(r), is formally defined as:

$\begin{equation} \kappa _{e,\lambda }= \frac{1}{\rho _g } \int \pi r^2 n(r) Q_e(r,\lambda) dr, \end{equation} \tag{ 2 }$

where Q_e is the extinction efficiency for a particle of radius r, at wavelength λ. The Rosseland mean opacity, which controls the flow of energy through an optically thick medium, is derived by an appropriate weighting of κ_{e, λ} over wavelength (see Appendix A). Extinction is due to a combination of pure absorption (reradiated as heat) and scattering (redirection of the incident beam). Rigorously, these components are additive; that is Q_e = Q_a + Q_s, where Q_a is the absorption efficiency and Q_s the scattering efficiency. The single-scattering albedo of the particle is ϖ = Q_s/Q_e. To the extent that the particles do scatter radiation without absorbing it (Q_s ≠ 0), the angular distribution of the scattered component (the phase function P(Θ)) is relevant. The first moment of the phase function g describes the degree of forward scattering:

$\begin{equation} g = \left\langle {\rm cos} \Theta \right\rangle ={\int P(\Theta) {\rm cos}\Theta {\rm sin} \Theta d \Theta \over \int P(\Theta) {\rm sin} \Theta d \Theta }. \end{equation} \tag{ 3 }$

Isotropic scattering results in g = 0; very small particles (Rayleigh scatterers) have a phase function proportional to cos²Θ, which also leads to g = 0. This is why scattering by tiny particles is often approximated as isotropic. Large particles have a strong forward scattering lobe due to diffraction; for such particles g may approach unity (cf. HT) and can even be approximated as unscattered (Irvine 1975). In many cases of thermal emission when r ≪ λ (Section 5.3), only the absorption/emission component is important:

$\begin{equation} \kappa _{a,\lambda }= \frac{1}{\rho _g } \int \pi r^2 n(r) Q_a(r,\lambda) dr. \end{equation} \tag{ 4 }$

Analytical expressions for Q_a and Q_s have been known for over a century in certain asymptotic regimes; the well-known Rayleigh scattering regime for particles much smaller than the wavelength is one example. Similarly, scattering properties of very large objects have long been well understood in terms of Lambertian and related scattering laws. Considerable amounts of computation have been devoted in recent years to determination of efficiencies and phase functions in the intermediate Mie scattering regime, where particle size is comparable to the wavelength.

One important point for our purposes, that is demonstrated by HT, is that many of the exotic fluctuations in scattering and absorption properties which characterize the Mie regime vanish if one averages over a broad distribution of particle sizes. Further smoothing occurs given a combination of realistic particle nonsphericity and random particle orientations. Furthermore, a large number of experimental studies (see also Bohren & Huffman 1983; Pollack & Cuzzi 1980 for references) have shown that integrated parameters such as Q and g are much less affected by shape effects than, for example, the phase function itself, and the parameters of interest for thermal emission and absorption (Q_a, Q_s, and g) vary smoothly between the "Rayleigh" and "geometric optics" regimes. This observed behavior provides a certain comfort level for the simple assumptions we make here. Our model does not treat P(Θ) at all, but constrains g directly based on a fairly well-behaved dependence on the size and composition of the scatterer (Section 3.2). We use g to correct Q_a and Q_s for energy which is primarily scattered forward, and thus does not participate significantly in directional redistribution of energy. Thus, our model would not be appropriate for applications where scattering dominates absorption and the directional distribution of scattered energy is of primary interest. Fortunately, primitive aggregate particles made of ice, silicates, carbon-rich organics and tiny metal grains are good absorbers and poor scatterers across most size ranges, but Q_s does contribute to κ_λ in some regimes.

3.2. The Model: Calculation of Efficiencies

A good exposition of the form of the particle absorption and scattering coefficients is presented by Draine & Lee (1984, henceforth DL); our expressions may be derived directly from those published by DL, and we will not repeat much of their presentation. DL describe the radiative properties in terms of absorption and scattering cross sections C_{a, s} = Q_{a, s}πr², which are directly related to particle volume in the limit r ≪ λ. The efficiencies may then be expressed in terms of the electric polarizability (VDH p. 73), which is a function of the (complex) particle dielectric constant epsilon = epsilon ₁ + i epsilon ₂ (DL Equation (3.11)). The particle refractive index m = n_r + in_i is the square root of the dielectric constant. In general, Q_a and Q_s depend on the particle shape and orientation as well as the material refractive indices. The "small dipole" limit, in which the particles are both highly nonspherical and have extremely large refractive indices, has been discussed by DL, Wright (1987), Fabian et al. (2001), and Min et al. (2008). Lindsay et al. (2013) show that even the shape of tiny solid grains can be diagnostic in high spectral resolution observations of strong absorption bands (Section 2). After a small amount of algebra, DL equations (3.3), (3.6), and (3.11) lead directly to values of Q_ak = C_ak/πr² and Q_sk = C_sk/πr², where subscript k refers to alternate orientations of a nonspherical particle relative to the wave vector:

$\begin{eqnarray} Q_{ak} &=& {2 \over r^2 \lambda }{V \over L_k^2} { \epsilon _2 \over \left(L_k^{-1} + \epsilon _1 -1\right)^2 + \epsilon _2^2 } \nonumber\\ &=& \frac{4}{3L_k^2} \left(\frac{2 \pi r}{\lambda }\right){ \epsilon _2 \over \left(L_k^{-1} + \epsilon _1 -1\right)^2 + \epsilon _2^2 } \end{eqnarray} \tag{ 5 }$

and

$\begin{eqnarray} Q_{sk} &=& {8 \over 3r^2}\left({2 \pi \over \lambda }\right)^4\left({V \over 4 \pi L_k}\right)^2 { (\epsilon _1 - 1)^2 + \epsilon _2^2 \over \left(L_k^{-1} + \epsilon _1 -1\right)^2 + \epsilon _2^2 } \nonumber\\ &=& \frac{8}{27L_k^2}\left(\frac{2 \pi r}{\lambda }\right)^4 { (\epsilon _1 - 1)^2 + \epsilon _2^2 \over \left(L_k^{-1} + \epsilon _1 -1\right)^2 + \epsilon _2^2 }, \end{eqnarray} \tag{ 6 }$

where V is the particle volume and

$\begin{equation} \epsilon _1 = n_r^2 - n_i^2; \hspace{36.135pt} \epsilon _2 = 2 n_r n_i. \end{equation} \tag{ 7 }$

Only the expression for Q_a is explicitly given in DL (compare our Equation (5) with their Equation (3.12)). For the case at hand, however, once any significant amount of particle accumulation occurs, aggregates are likely to be roughly equidimensional and, especially if the particles are porous, the average refractive indices are not large (see below), so we will treat Q_a and Q_s as independent of particle orientation and, setting L_k = 1/3, obtain

$\begin{equation} Q_{a} = 12\left(\frac{2 \pi r}{\lambda }\right){ \epsilon _2 \over (\epsilon _1 +2)^2 + \epsilon _2^2 } \end{equation} \tag{ 8 }$

and

$\begin{equation} Q_{s} = \frac{8}{3} \left(\frac{2 \pi r}{\lambda }\right)^4 { (\epsilon _1 - 1)^2 + \epsilon _2^2 \over (\epsilon _1 +2)^2 + \epsilon _2^2 } \end{equation} \tag{ 9 }$

(VDH p. 71, DL84; see also Equation (14) of Miyake & Nakagawa 1993). Neither Miyake & Nakagawa (1993) or Pollack et al. (1994) discuss Q_s, but it becomes important for porous, weakly absorbing particles in the transition regimes we wish to bridge with our model, where wavelength and particle size are not extremely different. In the geometrical optics regime, nonspherical shapes can increase the surface area per unit mass, and thus the extinction efficiency, but for moderate nonsphericity the effect is only some tens of percent (Pollack & Cuzzi 1980). This effect is also neglected for the present. For plausible aggregates, magnetic effects are not important (see however Appendix B).

In the limit where n_i ≪ n_r − 1, which is reasonable for ices, silicates, organics, and their ensemble aggregates, Equations (8) and (9) reduce to the simpler and more familiar forms:

$\begin{equation} Q_a = { 24 x n_r n_i \over \left(n_r^2 + 2\right)^2} \end{equation} \tag{ 10 }$

and

$\begin{equation} Q_s = {8 x^4 \over 3} {\left(n_r^2 - 1\right)^2 \over \left(n_r^2 + 2\right)^2}, \end{equation} \tag{ 11 }$

where x = 2πr/λ. These are essentially the classical expressions given by VDH (p. 70).

In our code we actually use the full equations (8) and (9) above for Q_a and Q_s for the radiative behavior of particles much smaller than the geometric optics regime; these equations are valid for arbitrary m (DL). Equation (8) for Q_a (linear in particle radius) is used for all sizes. However, in the Mie transition region (in which the particle size is comparable to a wavelength) we modeled Q_s using a function valid for n_r of order unity (the Rayleigh–Gans regime; VDH pp. 132–133, p. 182). That is, at a transition value x₀ we change from the DL expression (Equation (9)) to the Rayleigh–Gans regime expression (VDH p. 182 and final equation of Section 11.23):

$\begin{equation} Q_s = \frac{1}{2} (2x)^2(n_r - 1)^2 \left(1 + \left({n_i \over n_r -1}\right)^2\right). \end{equation} \tag{ 12 }$

The transition value x₀ = 1.3 comes from setting the two expressions equal to each other. This transition from x⁴ to x² dependence bridges the region between the Rayleigh and geometric optics regimes. It is shown in VDH (p. 177) that a coupled dependence on the so-called "phase shift" 2x(n_r − 1) (for n_i ≪ n_r − 1), which includes both n and r, captures the peak and shape of Q_s. These expressions are valid as long as n_r − 1 is not much larger than 1/x (Bohren & Huffman 1983, p. 140); this is the expected regime for aggregates and mixtures of plausible materials. Our expressions do not capture the Fresnel-like oscillatory behavior of Q shown by the more complete series expansions in VDH and the full Mie theory, but do capture the correct small- and large-x behavior, and the proper transition size/wavelength. This is an increasingly good approach as the particle size distribution gets broader and the particles become more absorbing (HT), in which case the resonance behavior near 2x(n_r − 1) ∼ several diminishes and indeed vanishes (see, e.g., Cuzzi & Pollack 1978). Aggregate particles considered here are very good absorbers in general, having very broad size distributions, but we show below that the technique works surprisingly well in even more challenging regimes.

A key simplifying assumption of our model is that, while Q_s and Q_e grow from small values as x increases, following Equations (8) and ((9) or (12)), they are truncated at values representing the geometrical optics limit, as discussed below. See for instance Figure 1 (left), where in the geometric optics limit 2πr/λ → ∞, Q_e = Q_s → 2. The figure represents lossless particles; if absorption is present, Q_s (which includes diffraction) trends downwards and asymptotes to a value closer to unity (HT Figure 9). Note that the exotic ripples are primarily seen for monodispersions and narrow size distributions, and vanish as broader size distributions are used.

**Figure 1.** Left: scattering efficiency *Q_s* of a particle with refractive indices *n_r* = 1.33, *n_i* = 0 as a function of its optical size x = 2πa/λ using Mie theory (figure reproduced directly from Hansen & Travis (1974), by permission, where particle radius = a). The parameter b varies the width of the size distribution; as the size distribution gets broad, exotic ripples dues to interference average out and vanish. Right: dependence on the optical size and material properties of the scattering asymmetry parameter g = 〈cos Θ〉, for narrow size distributions of different optical size. Note there are essentially two well-defined regimes in *n_i*. Figure 2 gives more detail on how g varies with x.
Download figure:
Standard image High-resolution image

Although the full radiative transfer problem must be solved when the angular variation of the intensity is of primary concern, so that the phase function P(Θ) is important, certain valuable simplifications are common when only angle-integrated quantities such as energy or flux are of concern. In such cases, one often merely corrects the extinction efficiency for the overall degree of forward- or back-scattering. We will adopt a common scaling (e.g., Van de Hulst 1980) sometimes referred to as the radiation pressure scaling, and adjust our efficiencies to take account of nonisotropic scattering as follows:

$\begin{equation} Q^{\prime }_e = Q_e - g Q_s = Q_e(1 - \varpi _0 g) = Q_a + Q_s(1 - g). \end{equation} \tag{ 13 }$

That is, to the degree that scattered radiation is concentrated into the forward direction, it may be regarded as unremoved from the beam. For instance, even a large scattering efficiency contributes nothing to the extinction if the scattering is purely forward directed (g = 1). Negative values of g (preferential backscattering) increase the extinction efficiency; in this case, it is even harder for radiation to escape than for thermalized or isotropically scattered radiation. This same correction is applied by Pollack et al. (1985, 1994); see Section 5.4 for more discussion. Modeling efforts relying to any degree on the Eddington approximation, in which scattered radiation is nearly isotropic, are perhaps better conducted using extinction efficiencies that are scaled in the way shown above, instead of using the formal Mie calculations (see Irvine 1975). For instance, de Kok et al. (2011) show how forward scattering by plausible exoplanet cloud particles can influence emergent near infrared fluxes; the above treatment of forward scattering partially accounts for this behavior by truncating away strong forward scattering and allowing the non-truncated scattered component to be treated as isotropic. In the approach presented here, the normal optical depth is defined using a value of κ_λ that has been adjusted for forward scattering effects following Equations (2) and (13) (and (17) if appropriate):

$\begin{equation} \tau _{\lambda }(z) = \int _z^{\infty } \kappa _{\lambda }(z) \rho _g(z) dz. \end{equation} \tag{ 14 }$

Instead of performing full Mie calculations to obtain the scattering asymmetry parameter g = 〈cos Θ〉, we make use of the crudely partitioned behavior illustrated by HT (reproduced in Figure 1 (right)). We have explored several simple ways of doing this, and have chosen the following determination of g, consistent with Figures 1 and 2 (notice that Figure 1 (right) is for a fairly low n_r = 1.33):

$\begin{equation} {\rm for}\,\, n_i \,{<}\, 1: g=0.7(x/3)^2 \hspace{7.22743pt} {\rm if} \hspace{7.22743pt} x \,{<}\, 3, \,\,{\rm and}\,\, g \approx 0.7 \hspace{7.22743pt} {\rm if} \hspace{7.22743pt} x \,{>}\, 3, \end{equation} \tag{ 15 }$

$\begin{equation} {\rm for}\,\, n_i \,{>}\, 1: g \approx -0.2 \hspace{7.22743pt} {\rm if} \hspace{7.22743pt} x \,{<}\, 3, \,\,{\rm and}\,\, g \approx 0.5 \hspace{7.22743pt} {\rm if} \hspace{7.22743pt} x \,{>}\, 3. \end{equation} \tag{ 16 }$

The value of the g-asymptote at large x, and the transition value of x, vary slightly with n_r (Figure 2); the values we selected apply to the EMT values of n_r ∼ 1.7–2.1 and n_i ≪ n_r, most relevant to our mixed-aggregate applications. We have not attempted to fine-tune either the asymptote or the transition value, and leave the large-n_i definition of g in its crude bimodal form (the large-n_i recipe has almost no effect on the results). Further fine-tuning is somewhat beyond what we expect of a simple utilitarian model, but could be done easily.

The final piece of the model involves matching the growing values of Q_s and Q_a in the Rayleigh and/or Rayleigh–Gans regime to constant values in the geometric optics regime, simply by limiting their magnitudes:

$\begin{equation} Q_a < 1 \hspace{14.45377pt} {\rm and} \hspace{14.45377pt} Q_s(1-g) < 1. \end{equation} \tag{ 17 }$

This can be thought of as limiting the absorption cross section, and the diffraction-corrected scattering cross section, to the physical cross section. D'Alessio et al. (2001) did not adjust their Q, albedo, or κ values for g as do we and Pollack et al. (1994), instead taking them straight from their Mie code. If one then carries on to solve the radiative transfer equations in an Eddington-like approximation (that is, as in their Equation (3), assuming isotropic scattering, or otherwise without detailed treatment of the angular dependence of often strongly forward-scattered radiation), indeed it would be better to first perform the truncation following Equation (13). A complete multidimensional treatment can capture all this behavior of course. D'Alessio et al. (2006) compared forward scattering with isotropic scattering, showing very little difference (their Figure 7),³ but at their longer wavelengths r ≪ λ so scattering is negligible in the first place, and at shorter wavelengths the emission arises from the photosphere of an opaque disk. A detailed discussion of this issue is beyond the scope of this paper.

One final detail must be mentioned, relevant for applications where particles of pure metal are important. Equations (8) and ((9) or (12)) include only the electric dipole interaction terms that are appropriate for everything but metals, and thus for all of our aggregate particle models where refractive indices are not extremely large. However, when we compare our model with the heterogeneous-composition grain mixture of Pollack et al. (1994), in which pure iron metal grains play a role, we have used a crude correction to Q_a for magnetic dipole terms as well, following DL84, Pollack et al. (1994), and others (see Appendix B); these calculations are all compared with the full Mie theory below.

3.3. Material Properties

For ease of comparison with previously published results, we adopt refractive indices, material abundances, and stability regimes for the condensible constituents as published by Pollack et al. (1994). We have merged the two different silicates (olivine and pyroxene) and the two different organics (higher and lower volatility) of Pollack et al. (1994) into a single population of each for simplicity. The exact composition of protoplanetary nebulae is, after all, poorly known, but these are consistent with known properties of grains in primitive meteorites and comets and cover an appropriate range of evaporation temperature. The components, their assumed relative abundances, and their evaporation temperatures are shown in Table 1. We will show that particle growth leads to at best a weak dependence on detailed stipulations of composition. It would be simple for anyone to select an alternate mixture.

Our refractive indices are shown in Figure 3. The materials can be classified into two groups—the water ice, silicate, and organic material have real refractive index n_r on the order of unity, and imaginary refractive index n_i much less than unity. On the other hand, metallic iron and iron sulfide particles have both refractive indices larger (often much larger) than unity.

4. MIXTURES OF GRAINS WITH PURE COMPOSITION

In our primary intended application to protoplanetary nebulae, particles are porous composites or granular mixtures of much smaller grains of all compositions that are solid at local temperature T, and we will calculate their ensemble refractive indices (n_r, n_i) using the Garnett EMT as described in Appendix C and Section 5.1 below. However, in some applications, such as exoplanet atmospheres, grains of a particular composition—even perhaps grains of iron metal—might be isolated within their own cloud layer.

To assess our approach in this regime, we first calculate opacities for a compositionally heterogeneous ensemble of micron-size grains adopted by Pollack et al. (1985, 1994) in their exact Mie scattering calculations. That is, distinct populations of non-porous grains with distinct compositions and refractive indices are assumed, and the opacities so determined are averaged as described below. This is a more challenging regime than a situation in which aggregates are of mixed composition and thus do not have extreme values of refractive index (Section 5.1).

For particles of species j, the monochromatic opacity at wavelength λ (relative to the gas mass density ρ_g) is:

$\begin{equation} \kappa _{e,j,\lambda } = {\beta _j \over \rho _g} \int n(r) \pi r^2 Q^{\prime }_{e,j}(r, \lambda) dr, \end{equation} \tag{ 18 }$

where n(r) = ∑n_j(r) and n_j(r) = β_jn(r); β_j is the fraction by number of particles of species j; that is, as assumed by Pollack et al. (1994) for this particular case, all species are assumed to have the same functional form for their size distribution. The Pollack size distribution used is a differential powerlaw for the particle number volume density in some radius bin dn(r, r + dr) = n(r)dr where n(r) = n_or^−p, p = −3.5 for 0.05 μm <r < 1 μm, p = −5.5 for 1 μm <r < 5 μm, and zero otherwise. Note that this size distribution is actually quite narrow: it covers roughly one decade in radius nominally (the piece from 1–5 μm contains very little mass or area), and is sufficiently steep for most of the surface area to be provided by the smallest members. In Appendix C we show that

$\begin{equation} \beta _j = { \alpha _j / \rho _j \over \sum (\alpha _j / \rho _j)}, \end{equation} \tag{ 19 }$

where α_j is the fractional mass of constituent j relative to the gas, and ρ_j is the density of solid constituent j. Note that β_j is the fractional number of particles of species j, and is subtly different in general from the fractional volume weighting factor f_j that is used when compositions are mixed together in the same particle as in our EMT approach (for compact or non-porous particles they are the same). The above definition of β_j differs slightly from that given in Pollack et al. (1985); their Equations (2)–(4) capture the right essence of the solution, but are not rigorously correct as written. By stepping back and incorporating the Q factors (including the correction for scattering asymmetry) inside the integrals of their Equation (4), the same result can be obtained:

$\begin{eqnarray} \kappa _{e,\lambda } &=& \frac{1}{\rho _g} \sum _j \beta _j \int n(r) \pi r^2 Q^{\prime }_{e,j}(r, \lambda) dr \nonumber\\ &=& {1 \over \rho _g} \int n(r) \pi r^2 \left(\sum \beta _j Q^{\prime }_{e,j}(r, \lambda)\right) dr, \end{eqnarray} \tag{ 20 }$

where $Q^{\prime }_{e,j}$ includes the forward scattering correction (Equation (13)). This is probably the most challenging test we could pose—an ensemble of heterogeneous particles, including some with very large refractive indices. Some are good absorbers and some are very good and rather anisotropic scatterers.

4.1. Wavelength-dependent Opacities

The fundamental calculation is of the effective, forward-scattering-corrected (or "radiation") extinction efficiency for each compositional type $Q^{\prime }_{ej}(r,\lambda)$ (Equation (13), using Equations (8) and (9) or (12)). If the situation calls for a heterogeneous mixture of grains with distinct compositions, $Q^{\prime }_{ej}(r,\lambda)$ must be calculated separately for each composition.

We have compared our model calculations directly with Mie calculations for five separate compositions, in each case averaged over the Pollack et al. (1994) particle size distribution (which is fairly narrow) as weighted by particle number density and geometric area. That is:

$\begin{equation} Q^{\prime }_{ej}(\lambda) = { \int n(r)\pi r^2 Q^{\prime }_{ej}(r,\lambda) dr \over \int n(r)\pi r^2 dr}. \end{equation} \tag{ 21 }$

These results are shown in Figure 4. The model does surprisingly well overall, except for metallic particles (FeS has a near-metallic behavior at short wavelengths). Note the different model curves in the iron metal panel; the solid curve attempts to correct for the effects of magnetic dipole interactions with the DL84 approximation also used by Pollack et al. (1994; Section 3.2). This term is only the leading term in an expansion, here used outside its formal realm of validity (see Appendix B). We would not advocate use of our simple model in cases where pure metal particle clouds might be encountered, such as in some exoplanet atmospheres (Lodders & Fegley 2002; Marley et al. 2013). However, it does improve agreement (especially at long wavelengths) with the full Mie calculations of similarly size-and-abundance averaged values of both $Q^{\prime }_e(\lambda)$ (bottom right panel of Figure 4) and κ_R (Figure 5 below), for the Pollack et al. (1994) heterogeneous grain mixtures. In fact, all panels in Figure 4 include calculations done with and without magnetic dipole terms, but the differences are insignificant for all materials except iron metal (there is a small effect for FeS), and in the averaged values (lower right panel). These individual compositional efficiencies and/or cross-sections are then combined into an overall opacity (Equation (20)), or simply an average efficiency, as weighted by their abundances:

$\begin{equation} Q^{\prime }_{e}(\lambda) = \sum \beta _j Q^{\prime }_{ej}(\lambda). \end{equation} \tag{ 22 }$

A plot of $Q^{\prime }_{e}(\lambda)$ is shown in the bottom right panel of Figure 4, also compared with the similarly summed Mie-derived values.

**Figure 4.** Comparison with Mie theory of our mean extinction efficiency $Q^{\prime }_{ej}(\lambda)$ for five pure materials, averaged and weighted over the Pollack et al. (1994) size distribution (see Section 4.1). Our model results are shown in solid curves and Mie calculations are shown in dashed curves. The lower right panel is the abundance-weighted average value $Q^{\prime }_e(\lambda)$ . The solid curve in the iron panel includes an (imperfect) correction for magnetic dipole effects (Appendix B), and the dotted (*Q_a*) and dot-dash (*Q_s*) curves do not. Actually, all the solid curves incorporate this correction term, but it makes a difference only for iron.
Download figure:
Standard image High-resolution image

**Figure 5.** Comparison between the Rosseland mean opacities as determined here, with much more elaborate full Mie scattering calculations (see text for discussion). The Pollack model has silicates, organics, and water ice evaporating at slightly different temperatures than assumed here (in fact they are a function of pressure); these are easily adjusted if desired.
Download figure:
Standard image High-resolution image

The more plausible case for protoplanetary nebulae, where all the particles are single (mixed) composition aggregates, with refractive index given by the Garnett EMT (see Appendix C and Section 5 below), is even simpler. A single efficiency $Q^{\prime }_{e}(r,\lambda)$ , associated with the EMT refractive indices, is simply integrated over the size distribution to give some wavelength-dependent opacity Q_e(λ). Whether determined from a compositionally heterogeneous or homogeneous mixture, $Q^{\prime }_e(\lambda)$ is the basis of extinction calculations such as Rosseland opacities (Equations (1) and (2)).

4.2. Rosseland Mean Opacities

We initially apply our models to calculate Rosseland mean opacities κ_R (Appendix A) for the standard Pollack et al. (1994) heterogeneous grain size distribution. Figure 5 shows that our results agree quite well with the exact Mie scattering calculations of Pollack et al. (1994), especially considering the simplicity of our approach. The Pollack et al. size distribution is an interstellar medium (ISM) distribution, with a steep powerlaw n(r) between radii of .005–1.0 μm and a steeper powerlaw from 1.0–5.0 μm (Section 4 above). Pollack et al. (1994) introduce two types of moderately refractory organics (with the same refractive indices but different abundances and evaporation temperatures) to provide greatly increased opacity between 160 K and 425 K, which we combine into a single component; thus we lack one evaporation boundary at about 260 K; a second difference is that we have only used a single silicate component. In view of these differences in detail, and for a situation like this (five distinct species including iron metal, with grain sizes in the resonance regime at the shorter wavelengths at high temperatures) which is far more challenging than our primary intended application (well-mixed aggregates of moderate refractive index, many in the geometrical optics limit), the agreement is surprisingly good.

5. AGGREGATE PARTICLES, EACH OF MIXED COMPOSITION

In the remainder of the paper, we explore the radiative properties of more mature and (at least for protoplanetary nebula applications) probably more realistic particle porosity and size distributions between a micron and arbitrarily large sizes. Evidence from meteorite and interplanetary dust particle (IDP) samples indicates that accumulated particles of these larger sizes are heterogeneous aggregates of all candidate materials, with individual submicron-to-micron size "elements" being composed of one mineral or other. Even mm-size chondrules are each a mixture of 1–10 μm size mineral grains of silicate, iron, and iron sulfide, and the ubiquitous meteorite matrix is a mix of smaller grains of all these materials, perhaps originally aggregates. So, for our aggregates, we assume only one (mixed) grain composition; the assumption may be shaky at the very smallest sizes which may be monomineralic, but coagulation models show that the tiniest particles are quickly consumed by growing aggregates (Ormel & Okuzumi 2013). Then, we are back to Equation (2):

$\begin{equation} \kappa _{e,\lambda } = \frac{1}{\rho _g} \int n_0(r)\pi r^2 Q^{\prime }_e(r,\lambda)dr. \end{equation} \tag{ 23 }$

As the gas density (and grain number density) increases, the collision rate increases accordingly. Grains of micron size are well coupled to the gas, and have fairly low collision velocities relative to each other (most recently Ormel & Cuzzi 2007). The low relative velocities imply that sticking will be fairly efficient at least until mm–cm sizes are reached (Dominik et al. 2007; Güttler et al. 2010; Zsom et al. 2010; Birnstiel et al. 2010, 2012). Several laboratory simulations of this process have shown that the resulting aggregates are of fairly low density (Donn 1990; Beckwith et al. 2000). It has been noted that aggregates of this sort are fractals in the sense that their volume depends on their mass to an arbitrary power. Normally this relationship is expressed as m∝r^s, where m is an individual particle mass and s is the fractal dimension (Wright 1987; Beckwith et al. 2000; Dominik et al. 2007). If the dimension is less than 3, their internal density decreases as the mass increases. Experimental results suggest that the internal densities and volume fractions of such particles are likely to be quite small, with porosities approaching 70% even after they begin compacting each other at increasing relative velocities (Weidenschilling 1997; Ormel et al. 2008 and references therein). Growth of this sort is certain to be robust in giant planet atmospheres as well (Ackerman & Marley 2001; Helling et al. 2008; Marley et al. 2013 and references therein).

5.1. Refractive Indices of Aggregate Particles

We envision each particle as compositionally heterogeneous—made up of much smaller constituents of specific composition and/or mineralogy, as seen in meteorites, IDP's, and cometary dust. That is, our particles are aggregates of subelements of species j, each smaller (and usually much smaller) than any relevant wavelength. The average particle refractive indices can then be calculated using an EMT approach. The two best known EMTs are due to Garnett and Bruggeman (Bohren & Huffman 1983). In the Garnett model, there is presumed to be one pervasive "matrix" in which distinct grains of other materials are embedded. In the Bruggeman theory, there is no structural distinction between domains of different refractive index. Bohren & Huffman (1983) believe that the Garnett rule is fundamentally to be preferred for aggregate particles where there is a well-defined matrix and it is vacuum (even, in principle, as the porosity gets very small). As this is our application, we use the Garnett EMT. The expressions and some additional discussion are given in Appendix C, and a set of effective refractive indices for each temperature range (each ensemble of condensed solids) is shown in Figure 6.

**Figure 6.** Left: refractive indices of aggregate particles, obtained using the Garnett EMT, as functions of wavelength for four different temperature regimes in which all materials are present (black), water ice has evaporated (red), organics have evaporated (green) and troilite has evaporated (blue). The temperatures at which these transitions occur are shown in Table 1. Right: ratio of solid/porous particle EMT refractive indices in the four temperature regions, for a porosity of 90%. For porosities in this range, these refractive indices could simply be used in our basic Equations (8) and ((9) or (12)), along with Equation (13) as constrained by Equations (15)–(17). These results would be used in our model to calculate particle opacities at temperatures where one or more of the five basic constituents has evaporated or decomposed (this is, of course, implicit in calculation of Rosseland mean opacities as a function of temperature, but the monochromatic opacities shown here are for a temperature where all candidate solids are condensed). Tabulated values available from the authors and will be posted online.(A color version and supplemental data of this figure are available in the online journal.)
Download figure:
Standard image High-resolution image

In the limit where all of the refractive indices are of order unity (this holds away from absorption bands for all species but the iron and troilite), an intuitively simple linear volume average captures the sense of the effect and might serve acceptably well in some regimes:

$\begin{equation} n_i = \frac{1}{f}\sum _j f_j n_{ij} \end{equation} \tag{ 24 }$

and

$\begin{equation} n_r = 1 + \frac{1}{f}\sum _j f_j (n_{rj} - 1). \end{equation} \tag{ 25 }$

Here, f is the volume fraction of all solids in a given particle, and f_j is the volume fraction for each species j. For low mass density particles such as some of those seen in our model distributions, the average imaginary index can get quite small and the real index can approach unity. However, in Appendix C and Section 5.2 below, we demonstrate that such a simple volume average greatly overestimates the contribution of even a small volume fraction of high-refractive-index grains (i.e., iron and troilite) in the mixture; for the purpose of this paper only the full Garnett theory (described in detail in Appendix C) is used.

5.2. Wavelength-dependent Opacities

In Figure 7 we compare weighted averages of $Q^{\prime }_e(\lambda)$ over heterogeneous mixtures, as described in Section 4, with values obtained for the same size distribution of aggregates in which the same material is mixed within each particle, and refractive indices are determined using the Garnett EMT (Section 5.1 and Appendix C). As expected for particles much smaller than the wavelength (Section 3.2), and as found by other investigators in the past, there is very little difference in the wavelength-dependent efficiency between a calculation in which each material is treated separately, and one in which they are merged together and treated with the EMT. The Pollack size distribution, while formally extending to 5 μm radius, is very steep and in reality has nearly all the cross section and mass at radii smaller than, if not much smaller than, 1 μm (Section 4).

The extremely high degree of agreement between the Garnett EMT and the heterogeneous mixture shown in Figure 7 was initially a surprise to us, given the presence of materials of high refractive index (see Appendix C). However, we believe this agreement actually validates not only the Garnett EMT itself (which hardly needs more validation) but also our numerical implementation of it. Consider the fact that for tiny particles with r ≪ λ, as described in Section 3.2, Q_e becomes proportional to the total volume or mass of material, regardless of the specific grain size distribution. This suggests that Q_e also becomes independent of the configuration of the grains; they can be dispersed or combined into clumps, as long as the clumps themselves remain much smaller than the wavelength, without affecting the result. Indeed this is the configuration implicit in the Garnett-averaged results for the Pollack size distribution; tiny particles are rearranged into aggregates of still-tiny particles. The good agreement testifies to the validity of the Garnett EMT in calculating a single set of effective refractive indices (at each wavelength) from which to calculate Q_e. That this is not trivial is demonstrated by the poor agreement obtained from the intuitively simpler linear volume mixing approximation to the mean refractive index (Appendix C; blue curve in Figure 7). Figure 12 in Appendix C shows in another way how badly volume mixing performs when even small amounts of high-refractive-index material are involved, as they are here. Of course, as particles grow larger, this equivalence will fail (a hint of divergence is seen at short wavelength in Figure 7). In this regime we will continue to rely on the Garnett EMT.

5.3. Monochromatic Opacities and Disk Masses

Perhaps the hardest challenge for our model is calculating monochromatic opacities for particle size distributions dominated by sizes comparable to the wavelength. Applications to wavelengths which are either much longer than the particle sizes in question (e.g., most of Figures 4 or 7), or much shorter (geometrical optics), are very reliable and straightforward. To better determine the limitations of our model, we explore monochromatic opacity for broad size distributions where many of the particles are comparable to the wavelengths of greatest interest to millimeter–centimeter observations of disks. The results (Figure 8) are comparable to those shown by Miyake & Nakagawa (1993, MN93) and D'Alessio et al. (2001, D01). Because MN93 and D01 each adopt somewhat different choices for refractive indices and compositional regimes, we do not compare our results directly to theirs, but instead conduct our own Mie calculations using the Pollack et al. (1994) refractive indices. D01 assumed a heterogeneous particle distribution (for instance, particles of pure ice, pure silicate, or pure troilite occur up to the mm–cm maximum sizes), whereas MN93 more typically assumed heterogeneous mixtures of aggregate particles with variable porosity, which we feel is more plausible under realistic conditions and also assume below.

In Figure 8 we compare monochromatic profiles of the true absorption opacity κ_{a, λ} from our model (Equation (4)) with full Mie calculations. MN93 have argued that scattering (which is important for extinction and thermal equilibrium modeling) is negligible in observations of this type and that pure absorption dominates thermal emission, and D01 agree.

Results are shown for differential size distributions n(r) = n_or^−s, with s = 2.0 and 3.1, with smallest radius of 1 μm and largest radius r_max. The value s = 2 (heavy curves) allows the larger particles to dominate the area and mass, while s = 3.1 (light curves) gives a more equitable distribution of area across the range of sizes. For each powerlaw slope, r_max is varied from 10 μm to 1 cm (the Mie calculations bog down for larger sizes). In the left panel, we show results for solid particles (internal density = 1.38 g cm⁻³). In the right panel, we show particles with 90% porosity. The solid curves are from our model and the dashed curves are full Mie calculations.

The agreement for porous particles is extremely good for the full range of sizes and size distributions. For solid particles (Figure 8 left), the Mie calculations exhibit an enhanced Q_a in the resonance region, covering perhaps a decade of wavelength around 2πr/λ ∼ 1 (see, e.g., HT). This effect is difficult for a model lacking complete physical optics to reproduce. Such "bumps" in κ are visible in Figure 5 of MN93, for the solid particle case, but lacking in their Figure 8 (for 90% porous particles), in good agreement with our results where all combinations of r_max and powerlaw slope reach the same long-wavelength asymptote rather quickly. The same "bump" is seen, for instance, in Figure 3 of Ricci et al. (2010), where the particles are actually not very porous (∼40%?). It is this contribution which carries through at shorter wavelengths, leaving our opacities 20% low or so relative to Mie values; this small difference, in a wavelength-independent regime, is probably insignificant for most purposes when the maximum particle size is probably close to the specific observing wavelength and the porosity is not large, as discussed more below.

The more shallow 2.0 powerlaw, dominated in area and mass by particles at or near the upper size cutoff, abruptly transitions to wavelength-independent opacity at the wavelength where, roughly, Q_a(r_max, λ) = 1. At longer wavelengths, a universal curve is followed, characterized by the total mass in the system. Notice that, for the 2.0 powerlaws, the inflection point moves to shorter wavelengths and the short-wavelength opacity increases as the upper radius cutoff decreases or as the porosity increases. This is a direct implication of Equations (8) and (7), which together imply that Q_a∝rn_i, and in most cases, in spite of the imperfection of the linear volume mixing model, n_i is roughly proportional to (1 − ϕ) (Figure 12). Meanwhile, for wavelengths shorter than the Q_a(r_max, λ) = 1 turnover, each particle of radius r is less massive than a solid particle by the factor (1 − ϕ), so the number of particles is larger (for a given total mass) by a factor 1/(1 − ϕ) (see also Section 5.4 below). That is, fluffy cm-size particles have ten times the opacity of solid particles of the same radius (in the short-wavelength limit) because their 10 times lower mass per particle allows there to be ten times as many of them than the solid particles of the same size. Their mass per particle is still 100 times larger than that of a mm-size solid particle, but their cross-sectional area per particle is 100 times larger, so the curve for porous, cm-radius particles lies on top of the curve for solid, mm-radius particles even at short wavelengths. This behavior can be seen tracing the differences between the heavy red and black curves, through solid, dashed, and dotted manifestations.

The lightweight curves, for the steeper 3.1 powerlaw size distribution having the same radius limits, contain more small particles and thus show more spectral signatures at mid-infrared wavelengths. They also show broader slope transitions in $\kappa ^a_{\lambda }$ around Q_a(r_max, λ) = 1. Note the envelope of slope for all combinations of powerlaw, upper size limit, and porosity, for which Q_a(r_max, λ) = 1 has not been reached. Opacities (fluxes) at these wavelengths are independent of particle size, and thus capture all the mass even if particles have grown to cm size. Different powerlaw slopes or other details of the size distributions could lead to slightly different functional forms in the transition region, and can thus be constrained by high-quality observations. These comparisons of our model monochromatic opacities with full Mie theory provide further support for the Rosseland mean opacities which are based on them, certainly for porous particles, and illustrate the extent of the limitations for solid particles. In applications such as mm–cm monochromatic SED slope analysis, with the goal of determining largest particle sizes, full Mie theory should be used, and are not burdensome here because the wavelengths of interest are not much smaller than the particle sizes of interest.

These scalings with size and porosity are further illustrated in Figure 9, comparing our runs for porous particles, which agree well with Mie calculations but are easier to extend to large sizes, with actual Mie calculations for solid particles. In Figure 9 we show a widely used opacity rule (Beckwith et al. 1990; Williams & Cieza 2011): κ_ν = 0.1(ν/10¹² Hz)^β, with β = 1. Note that the canonical Beckwith et al. (1990) opacity could be fit by a suite of porous particles extending to 100 cm radius, especially with a slope of 3.1 (or perhaps slightly steeper). As noted by previous authors, centimeter-size solid particles also provide a fairly good match to the slope. All models that match the slope have opacity quantitatively smaller than the canonical B90 opacity, however, by a factor of at least several (also found by Birnstiel et al. 2010) possibly suggesting larger inferred disk masses.

In outer disks where gas densities are very low, particles of even mm size and low-to-moderate porosity are dynamic radial migrators under the influence of gas drag (Takeuchi & Lin 2005; Brauer et al. 2008; Hughes & Armitage 2010, 2012; Birnstiel et al. 2010, 2012) and observations apparently show radial segregation of gas and particles in more than one outer disk (Andrews et al. 2012; Pérez et al. 2012). Moreover, such particles also couple to the large, high-velocity eddies and would be expected to have fairly high collision velocities. In this regime, large particles, even with 90% porosity, would collide at significant velocities, inconsistent with retaining their postulated high porosities. This is because the aerodynamic coupling of a particle to turbulence is determined by the product of its radius and density (see, e.g., Völk et al. l980; Cuzzi & Hogan 2003; Dominik et al. 2007; Ormel & Cuzzi 2007). Thus it seems more plausible that the particles matching the B90 opacity in shape (or something like it) are indeed moderately compact, cm-radius, objects, and not meter-size, high-porosity puffballs; thus full Mie theory is needed for analyses of these observations.

5.4. Rosseland Mean Opacity: Effects of Size and Porosity

In this section we continue to assume well-mixed aggregates of material in each particle and extend our calculations of Rosseland mean opacities to particles with a wider range of size and porosity. First we show how growth affects the Rosseland mean opacities introduced in Section 3.1 and Appendix A. Clearly, for particle size much larger than the wavelength, the opacity varies as the ratio of cross section to mass, or as 1/r (see Appendix of Pollack et al. 1994 for a longer exposition). As shown by Pollack et al. (1985) and Miyake & Nakagawa (1993), growth from microns to centimeters implies a decrease in κ_R by four orders of magnitude, dwarfing any uncertainties regarding the actual composition of the particles and rendering the small differences seen in Figures 4–7 somewhat moot. We illustrate the effect using Figures 10 and 11. The size distributions in Figure 10 are monodispersions, which we have cautioned about previously, but in these calculations the large range of wavelengths over which efficiencies are integrated mimics the effects of a broad size distribution.

**Figure 10.** Rosseland mean opacities for various particle sizes and porosities, all calculated using our approach. Green: heterogeneous ISM distribution of Pollack et al. (1994); black curves: solid aggregate particles of no porosity with the same mass; solid: 10 μm radius monodispersion; dashed: 100 μm radius monodispersion; dotted: 1000 μm radius monodispersion (multiplied by 10). Also shown are red: 100 μm radius monodispersion with the same mass but porosity of 70%; blue: 100 μm radius monodispersion with the same mass but porosity of 90% (discussed in Section 5.2).
Download figure:
Standard image High-resolution image

**Figure 11.** Rosseland mean opacities for the same powerlaw size and porosity distributions as used in Figure 8. Both linear and logarithmic forms are presented for the same results. Green: heterogeneous ISM distribution of Pollack et al. (1994); black curves: solid aggregate particles of no porosity with the same mass and powerlaw size distributions of differential slope 3.1; solid: 1–10 μm radius; dashed: 1–100 μm radius; dotted: 1 μm–1 mm radius; dot-dash: 1 μm–1 cm radius (shown are red: 1 μm–1 cm radius with porosity of 90%.
Download figure:
Standard image High-resolution image

Particle growth from the ISM distribution of Pollack et al. (containing many submicron grains which are effective at blocking shortwave radiation characteristic of the higher temperatures shown), rapidly decreases the Rosseland opacity, even when growth is only to radius of 10 μm (solid black curve). The decrease continues as particles grow to 100 μm (dashed black) and 1 mm (dotted black, multiplied by 10) radii. The opacity curves for the larger particles have a characteristic stepped appearance, with the steps representing changes in solid mass fraction at different evaporation temperatures. The opacity is nearly constant within a step, because particles much larger than a wavelength are in the constant-Q_e regime independent of wavelength. Comparing the black dashed and dotted curves, where the dotted curve is 10 times the value of κ_R for 1 mm radius particles, shows that the 1/r scaling is nearly exact for these sizes. The scaling is only slightly less easily explained between the 10 μm and 100 μm radius solid particles. At the lower temperatures (longer wavelengths) the 10 μm particles have not yet reached the geometrical optics limit, which occurs roughly for temperatures higher than 300 K where the blackbody peak wavelength is around 10 μm. Thus, the scaling between these two sizes is less than a factor of ten by an amount that depends on temperature (wavelength).

Figure 10 also shows the expected effects of porosity (Section 5.3); the red and green curves, for 100 μm radius monodispersions of particles having porosities of 70% and 90% respectively, show good agreement with the expected scaling by 1/(1 − ϕ); that is, to conserve mass if the particle density decreases to $(1-\phi)\overline{\rho }$ , we must increase their number density by a factor of 1/(1 − ϕ). In this regime where the particles are already much larger than the most heavily weighted wavelengths in question (except at the far left of the plot), the Q_e per particle does not change, so the net opacity of the ensemble increases by a factor of 1/(1 − ϕ). So, the particle growth and porosity effects are robust, easily explained using simple physical arguments, and captured by the model.

The Rosseland mean opacities for several wide powerlaw size distributions are shown in Figure 11. The behavior is similar overall to behavior seen previously and explained above. These powerlaws have a smallest particle radius (1 μm) which is slightly larger than the typical size on the Pollack et al. size distribution, so κ_R is somewhat larger at low temperatures (long effective wavelengths) but slightly lower and flatter at high temperature (short effective wavelengths). Powerlaws extending to larger size reduce the overall opacity at higher temperature, at least, by tying up more mass in larger particles. Porosity increases the overall opacity everywhere for the powerlaws extending to the larger sizes, because the dominant wavelengths involved are generally tens of microns or less. However, some wrinkles in the details, and differences between Figures 10 and 11, arise from the scattering terms due to Q_s(1 − g) in κ_R, which become important for particles with low n_i when the particle size and wavelength are comparable.

It is these drastic decreases in opacity with grain growth (Movshovitz & Podolak 2008) that led Movshovitz et al. (2010) to conclude that gas giant formation could have happened much earlier than previous estimates which had already assumed an arbitrary 50× cut in opacity relative to the Pollack et al. (1994) baseline (Hubickyj et al. 2005; Lissauer et al. 2009). The reason is that growth in the most important part of the radiative zone proceeds to centimeter-size, implying opacities even 10× smaller there than the smallest (for solid mm-size particles) shown in Figures 10 and 11. If the particles are porous, their opacity could be increased, however, as shown in those figures. This would seem to be an example of how the destiny of the great can be determined by the behavior of the small.

6. CONCLUSIONS

We outline a very simple, closed-form radiative transfer model which incorporates and connects well-understood asymptotic behavior for particles smaller than and larger than the wavelength. The approach is simple enough to provide good physical insight and to include in evolutionary models requiring radiative transfer. The model is easily adapted to arbitrary combinations of particle size, composition, and porosity across the range of plausible protoplanetary nebula and exoplanet cloud particle properties (excepting highly elongated particles), and yields values of Rosseland mean opacity which are in good agreement with more sophisticated but more time consuming Mie or DDA calculations. Planck opacities are even simpler to calculate (they are straight Blackbody-weighted means over wavelength). We illustrate the significant roles of particle growth and porosity in determining opacity. The model is not recommended even for Rosseland opacities in cases where ensembles of large, pure metal particles are expected. The method even gives very good approximations to monochromatic opacities unless the particles are solid and the specific wavelength of interest is comparable to the dominant particle size, where the model is unable to track "resonance" behavior of the absorption efficiency. Thus for detailed spectral index analysis of millimeter–centimeter wavelength protoplanetary nebula emission spectra, in cases where particles in the millimeter–centimeter radius range might be solid and contribute significant mass and area, full Mie theory should be used. It appears that canonical millimeter-wavelength spectral slopes are more plausibly explained by solid, centimeter-size particles than larger, but more porous, aggregates.

To summarize, the model consists of Equations (8) and ((9) or (12)), along with Equation (13) as constrained by Equations (15)–(17). These equations lead to a monochromatic opacity κ_{e, λ}, which is then used to calculate a Rosseland mean κ_R using Equation (A5). For calculation of mm–cm SEDs, only $\kappa ^a_R$ should be used (Equation (4)). Our treatment can assume either a heterogeneous "salt and pepper" particle size and compositional mix, as in Pollack et al. (1994) or D'Alessio et al. (2001) (with Equations (18)–(20)), or a distribution of internally mixed, homogeneous aggregate particles having arbitrary porosity (using Equations (C6)–(C11) in Equations (8)–(17), and integrating with Equation (23)), as for instance modeled by Miyake & Nakagawa (1993).

J.N.C. benefitted greatly from many discussions of radiative transfer with Jim Pollack over the years. We thank Ted Roush for providing tabular values of the refractive indices used by Pollack et al. (1994). We thank Pat Cassen, Ke Chang, Tom Greene, Lee Hartmann, Stu Weidenschilling, and Diane Wooden for helpful conversations. We thank Kees Dullemond and Naor Movshovitz for encouragement to make this work more widely available. We thank our reviewer for helpful comments. J.N.C. and P.R.E. were supported under a grant from the Origins of Solar Systems Program. S.S.D. was partially supported by Grants from the LASER and Astrobiology programs.

Note added in proof. A recent submission takes an approach similar to ours (Kataoka et al. 2014).

APPENDIX A: ROSSELAND MEAN OPACITY

While derivations of the Rosseland mean opacity can be found in various sources, we provide here a quick derivation of κ_R because it is surprisingly difficult to find good ones, and we make reference to some of its specific aspects. Essentially one derives an expression for the monochromatic flux in a high-opacity thermal radiation field, then integrates over frequency, and from this identifies the associated mean opacity. Start with the radiative transfer equation for monochromatic intensity I_ν in a plane layered medium: μdI_ν/dτ = S_ν − I_ν where S_ν is the source function and τ_ν is monochromatic optical depth measured normal to the layer, and μ = cos θ where θ here is the angle from the layer normal and dτ_ν = κ_νdz where dz is an increment of thickness in the layer. In a medium of high optical depth where thermal radiation dominates everything else, I_ν ≈ S_ν ≈ B_ν where B_ν is the Planck function. Then we can rewrite the above equation to first order as I_ν = B_ν − dB_ν/dτ. We then determine the local energy flux through the layer $F_{\nu } = 2\pi \int _{-1}^1 I_{\nu } \mu d\mu$ . In the presence of the gradient derived above, substituting for I_ν:

$\begin{equation} F_{\nu } = 2 \pi \left[ \int _{-1}^1 B_{\nu }\mu d\mu - \int _{-1}^1{ \mu dB_{\nu } \over \kappa _{\nu } dz} \mu d\mu \right]. \end{equation} \tag{ A1 }$

The first integral vanishes because B_ν = constant; the flux becomes

$\begin{equation} F_{\nu } = -\frac{2 \pi }{\kappa _{\nu }} \frac{dB_{\nu }}{dz}\int _{-1}^1\mu ^2d\mu = -\frac{4 \pi }{3 \kappa _{\nu }}\frac{dB_{\nu }}{dz}. \end{equation} \tag{ A2 }$

The standard trick is to set dB_ν/dz = (dB_ν/dT)(dT/dz). The monochromatic flux is then integrated over frequency, after rearranging terms:

$\begin{equation} F = \int F_{\nu } d\nu = -\frac{4 \pi }{3}\frac{dT}{dz}\int \frac{1}{\kappa _{\nu }}\frac{dB_{\nu }}{dT}d\nu. \end{equation} \tag{ A3 }$

One then merely asserts that the frequency integrated flux can be written in the same form, except with a mean opacity κ_R:

$\begin{equation} F = -\frac{4 \pi }{3}\frac{dT}{dz}\frac{1}{\kappa _R}\int \frac{dB_{\nu }}{dT}d\nu; \end{equation} \tag{ A4 }$

and after setting the two expressions equal, we obtain the definition of κ_R:

$\begin{equation} \frac{1}{\kappa _R} = \frac{\int \frac{1}{\kappa _{\nu }}\frac{dB_{\nu }}{dT}d\nu }{\int \frac{dB_{\nu }}{dT}d\nu }; \end{equation} \tag{ A5 }$

essentially, we are obtaining the weighted average of 1/κ_ν, thereby emphasizing spectral regions where energy "leaks through," and where the integral of the weighting function dB_ν/dT in the denominator may, if we like, be further simplified as d/dT(∫B_νdν) = 4σ_SBT³, where σ_SB is the Stefan–Bolzmann constant.

APPENDIX B: METAL PARTICLES

The small particle expansions of DL84, leading to our primary Equations (8) and (9), incorporate only electric dipole terms. Similar expressions can be derived for magnetic dipole terms, which are more cumbersome in their full glory (see, e.g., Tanner 1984 or Ossenkopf 1991) and are usually approximated. For instance, DL84 (Equation (3.27)), citing Landau & Lifshitz (1960, Sections 45, 72, and 73), give a handy first-order correction factor for Q_a with the caveat that it is valid for small x, but without specific guidance as to what is "small." The correction is simply $Q^{\prime }_a(r,\lambda) = Q_a(r,\lambda)(1 + F)$ , where

$\begin{equation} F = \left(\frac{2 \pi r}{\lambda }\right)^2{(\epsilon _1+2)^2 +\epsilon _2^2 \over 90}, \end{equation} \tag{ B1 }$

and $\epsilon _1 = n_r^2 -n_i^2$ and epsilon ₂ = 2n_rn_i (Equation (7)). Similar expressions appear in Ossenkopf (1991) and Tanner (1984), and surely elsewhere. The factor F gets very large when the refractive indices are large. In the case of iron metal, where n_r and n_i are nearly proportional to λ, the spectral behavior of Equation (8) is flattened from x⁻³ nearly to x⁻¹, as can be seen by expanding terms in the large-refractive-index limit.

This magnetic dipole correction is insignificant in the mixed-material, aggregate particle case, as it only enters when the refractive indices of a particle are ≫λ/r, so readers interested only in aggregate grains need not be concerned further with this term. Yet, for those who might be interested in clouds of metallic particles, the correction might be of interest. Unfortunately, it is formally only valid for particles that are small compared to the wavelength inside the particle, or equivalently the skin depth of the wave, which is given by $x \sqrt{\epsilon } \ll 1$ , and for the large epsilon of metals (see Figure 3) the allowed x is much smaller than the usual condition x = 2πr/λ ≪ 1. The correction is mentioned by Pollack et al. (1994) but it is unclear from that paper just how it was used; their Figure 2(a) shows results for metal particles having their standard size distribution in a wavelength range where it is not formally valid, and they give no comparison with Mie theory. We have found that if Equation (B1) is blithely applied for 2πr/λ < 1, well out of its formal domain of validity, and the ensuing Q_a is subject to our overall constraint Q_a < 1, the agreement with Mie calculations for pure metal is improved from nonexistent to marginal (Figure 4). Tanner (1984) has found similar behavior, in that use of the "approximation" outside of its formal domain of validity gives surprisingly better agreement with observed behavior than application of the more complete theory. This is not an argument for general acceptance of the approximation, and in any application where metal particles dominate the situation the complete Mie theory is probably required. Nevertheless, we have employed it in our model.

APPENDIX C: GARNETT EFFECTIVE MEDIUM THEORY

In this appendix, we give an overview of the Garnett theory for the average complex dielectric constant ( epsilon ) of an inhomogeneous medium. The reader is referred to Bohren & Huffman (1983) for a more complete exposition with background.

A number of exhaustive and sophisticated studies have compared several different kinds of EMT models to rigorous, brute-force numerical DDA models (Perrin & Lamy 1990; Ossenkopf 1991; Stognienko et al. 1995; Voshchinnikov et al. 2005, 2006; see also Semenov et al. 2003 and references therein); the differences are generally small and composition-dependent for small particles which can be treated using DDA (and much smaller than the huge differences due to growth which our model is primarily intended to capture). For instance, Figure 16 of Voshchinnikov et al. (2005) shows a nearly insignificant difference between Garnett and Bruggeman models relative to DDA calculations, in the regime where all scattering elements in the aggregates are truly small compared to the wavelength. Interestingly, they show that a model of their own device does a better job matching certain very porous aggregates containing monomers with a distribution of sizes (at least, at long wavelengths). We believe that both traditional EMT theories fail to match the "size distribution of inclusions" DDA results of Voshchinnikov et al. (2005) Figure 16, because the aggregates contain numerous embedded wavelength-sized monomers, which violate the assumptions of both EMT models (some of the monomer inclusions have diameter as large as the wavelength). For this situation to occur at wavelengths, particles, and temperatures of interest for Rosseland mean opacities of our paper would require monomers in the few-to-tens of micron size, that are in turn embedded by assumption in much larger particles—which our model would treat in the geometric optics limit in any case. In Sections 5.3 and 6, we note deviations from our model when a significant fraction of the particles of interest are wavelength-sized (we anticipate this to be a problem mostly for millimeter–centimeter wavelength observations). The behavior is similar to that seen in Figure 16 of Voshchinnikov et al. (2005).

We derive the components of epsilon = epsilon ' + i epsilon '' and convert the results to the perhaps more familiar complex refractive index n' + in''. We then assess the realm of validity for a linear approximation of refractive index as a function of material density by comparing the two expressions for a range of component materials.

Interactions among different constituents of an inhomogeneous medium make the determination of an average dielectric constant problematic. The problem is generally insoluble, save by brute force (DDA) or approximation methods, and different approximations inevitably lead to different expressions. As an example, one choice that can be found in the literature as far back as 1850 is the Rayleigh expression that relates the density and dielectric constant of a powder to that of the corresponding solid. This result follows from the Clausius–Mosotti law that the quantity epsilon − 1/ epsilon + 2 is proportional to the density of a material:

$\begin{equation} \frac{1}{\rho }\frac{\epsilon -1}{\epsilon +2} = \frac{1}{\rho _{o}}\frac{\epsilon _{o}-1}{\epsilon _{o}+2}, \end{equation} \tag{ C1 }$

where epsilon _o, ρ_o are the complex dielectric constant and density of the solid material, and epsilon , ρ are those of the powder. This expression has been shown to be fairly accurate for powders made from various geological materials (e.g., Campbell & Ulrichs 1969). A more complete expression, which reduces to the Rayleigh formula in the case of one component, was derived by Maxwell Garnett (see, e.g., Bohren & Huffman 1983). Garnett's model is that of inclusions of dielectric constant epsilon _o embedded in a homogeneous medium of dielectric constant epsilon _m. In the simplest version, the inclusions are assumed to be identical in composition, but may vary in shape, size, and orientation. Assuming that all the inclusions are spherical, the expression for the average dielectric constant epsilon is given by

$\begin{equation} \epsilon = \epsilon _m \left[ 1 + 3f \left(\frac{\epsilon _o - \epsilon _m}{\epsilon _o + 2\epsilon _m} \right) \left(1 - f \left(\frac{\epsilon _o - \epsilon _m}{\epsilon _o + 2\epsilon _m} \right) \right) ^{-1} \right], \end{equation} \tag{ C2 }$

where f is a mass fraction which is defined below. For epsilon _m = 1, this reduces to the Rayleigh formula with f = ρ/ρ_o (for a more detailed derivation, including nonsphericity effects, see Bohren & Huffman 1983). The advantage of the Garnett equation is that it can be generalized to a multiple component medium. The general result can be cast into the same form as Equation (C2) except now f and epsilon _o become f_j and epsilon _j, while the numerator and denominator of Equation (C2) are summed over j species (Bohren & Huffman 1983, Section 8.5).

Noting that epsilon = epsilon ' + i epsilon '', we expand Equation (C2) into its complex parts, convert to fractional form, and separate the real and imaginary components:

$\begin{equation} \epsilon = \epsilon ^{\prime }+i\epsilon ^{\prime \prime } = \frac{1+ 2\sum _{j} f_{j}\sigma _{j} + i6\sum _{j} f_{j}\gamma _{j}}{1 - \sum _{j} f_{j}\sigma _{j} - i3\sum _{j} f_{j}\gamma _{j}}, \end{equation} \tag{ C3 }$

where we have defined σ_j and γ_j to be

$\begin{equation} \sigma _{j}=\frac{(\epsilon ^{\prime }_{j}-1)(\epsilon ^{\prime }_{j}+2) + \epsilon ^{\prime \prime 2}_{j}}{(\epsilon ^{\prime }_{j}+2)^{2} + \epsilon ^{\prime \prime 2}_{j}}, \end{equation} \tag{ C4 }$

and

$\begin{equation} \gamma _{j}=\frac{\epsilon ^{\prime \prime }_{j}}{(\epsilon ^{\prime }_{j}+2)^{2} + \epsilon ^{\prime \prime 2}_{j}}, \end{equation} \tag{ C5 }$

respectively. We apply the complex conjugate to Equation (C3), determine the real and imaginary parts of epsilon , and combine like terms to get

$\begin{equation} \epsilon ^{\prime }=\frac{1+\sum _{j} f_{j}\sigma _{j} - 2\sum _{i} \sum _{j} f_{i}f_{j}(\sigma _{i}\sigma _{j}+9\gamma _{i}\gamma _{j})}{1 - 2\sum _{j} f_{j}\sigma _{j}+\sum _{i} \sum _{j} f_{i}f_{j}(\sigma _{i}\sigma _{j}+ 9\gamma _{i}\gamma _{j})}=\frac{N_{R}}{D}, \end{equation} \tag{ C6 }$

for the real part, and

$\begin{equation} \epsilon ^{\prime \prime }=\frac{9\sum _{j} f_{j}\gamma _{j}}{1 - 2\sum _{j} f_{j}\sigma _{j}+\sum _{i} \sum _{j} f_{i}f_{j}(\sigma _{i}\sigma _{j}+ 9\gamma _{i}\gamma _{j})}=\frac{N_{I}}{D} \end{equation} \tag{ C7 }$

for the imaginary part. A simple check will show that this generalized Garnett formula reduces to Equation (C1), the Rayleigh formula, for i = j = 1.

Conversion of Equations (C6) and (C7) to an expression in terms of refractive indices poses no great problem. Making the usual assumption that magnetic permeability is of order unity (see, however, Appendix C) we define

$\begin{equation} \epsilon ^{\prime }+i\epsilon ^{\prime \prime }=(n^{\prime }+in^{\prime \prime })^{2}, \end{equation} \tag{ C8 }$

where now, Equations (C4) and (C5) become

$\begin{equation} \sigma _{j}=\frac{\left(n^{\prime 2}_{j}-n^{\prime \prime 2}_{j}-1\right) \left(n^{\prime 2}_{j}-n^{\prime \prime 2}_{j}+2\right) + 4n^{\prime 2}_{j}n^{\prime \prime 2}_{j}}{\left(n^{\prime 2}_{j} - n^{\prime \prime 2}_{j}+2\right)^{2} + 4n^{\prime 2}_{j}n^{\prime \prime 2}_{j}}, \end{equation} \tag{ C9 }$

and

$\begin{equation} \gamma _{j}=\frac{2n^{\prime }_{j}n^{\prime \prime }_{j}}{\left(n^{\prime 2}_{j} - n^{\prime \prime 2}_{j}+2\right)^{2} + 4n^{\prime 2}_{j}n^{\prime \prime 2}_{j}} \end{equation} \tag{ C10 }$

respectively. Finally, we express the average refractive index of the inhomogeneous medium in terms of epsilon as

$\begin{equation} n^{\prime }+in^{\prime \prime }=\left[\frac{\sqrt{\epsilon ^{\prime 2} + \epsilon ^{\prime \prime 2}}+\epsilon ^{\prime }}{2}\right]^{1/2} + i\left[\frac{\sqrt{\epsilon ^{\prime 2} + \epsilon ^{\prime \prime 2}}-\epsilon ^{\prime }}{2}\right]^{1/2}. \end{equation} \tag{ C11 }$

The Garnett formula defines the f_j as the volume fraction of inclusions of species j within a particle of some overall volume V and total mass M; then the total solid volume fraction in the particle is f = ∑_jf_j = 1 − ϕ where ϕ is the porosity of the particle. We can determine f_j in terms of mass fractions as follows: v_kj is the kth volume element of species j so that the volume fraction of inclusions of species j is f_j = ∑_kv_kj/V. We define α_j as the mass of all component j per unit nebula gas mass, and α = ∑_jα_j is the mass fraction of all solids per unit nebula gas mass. With these definitions in mind,

$\begin{equation} f_{j}=\frac{v_{j}}{V}=\frac{m_{j}}{\rho _{j}V}= \frac{m_{j}\rho }{\rho _{j}M}=\frac{\rho }{\rho _{j}}\frac{\alpha _{j}}{\alpha } = \frac{(1-\phi) \overline{\rho } \alpha _j}{\rho _j \alpha } \end{equation} \tag{ C12 }$

where ρ is the average mass density of a composite particle; ρ can vary from nearly zero to the "solid" average value $\overline{\rho }$ at ϕ = 0. Using the last expression of Equation (C12) in ∑_jf_j = f = 1 − ϕ gives

$\begin{equation} \overline{\rho }=\frac{\alpha }{\sum _{j} (\alpha _{j}/\rho _{j})}. \end{equation} \tag{ C13 }$

Under certain conditions, a simple linear approximation for the refractive index at a particular wavelength λ can greatly simplify calculations (see Equations (22) and (23)). Generally, such an approximation works well when n is of order unity (such as visible wavelengths).

Equations (22) and (23) follow from a simple, but approximate, extension of the so-called "Wiener rule" (Voshchinnikov et al. 2005): 〈 epsilon 〉 = Σf_j epsilon _j, where epsilon _j are the dielectric constants of the different materials j and the relation assumes lossless materials (n_i = 0). Equations (22) and (23) are simply the analogous linear volume averages applied to "moderate" refractive indices not far from unity, as is true for most silicates and ices. Assuming that the Wiener rule can be applied to a case where n_i ≠ 0, Equations (22) and (23) can be derived in an approximate way, in the limit that n_i ≪ 1 and n_r = 1 + δ with δ² ≪ 2δ. Given 〈 epsilon 〉 = Σf_j epsilon _j, separate the real and imaginary parts of the sum recalling that epsilon = (n_r + in_i)², giving $\left\langle \epsilon _i\right\rangle = \Sigma f_j (2 n_i n_r) \sim 2\overline{n_r}\Sigma f_j n_{ji}$ . The LHS of the equation can also be approximated as $2\overline{n_r}\left\langle n_i\right\rangle$ , so 〈n_i〉 ∼ Σf_jn_ji as in Equation (23). Similarly $\left\langle \epsilon _r\right\rangle = \left\langle (1+\overline{\delta })^2\right\rangle \sim 1+2\overline{\delta }$ , and on the RHS Σf_j(1 + δ_j)² ∼ Σf_j(1 + 2δ_j). Separating the sum on the RHS and noting that Σf_j = 1, we get 〈δ〉 ∼ Σf_jδ_j as in Equation (22). Of course, these are only approximations and valid only under the conditions stated.

Figure 12 shows a direct comparison of both the linear approximation and the generalized Garnett function, using the real and imaginary refractive indices of several combinations of likely materials at λ = 100 μm as a function of the density ratio $\rho /\overline{\rho }$ . Curve (a) corresponds to a single component of ice. Curves (b) and (c) are two-component mixtures of ice/rock and ice/iron. The fraction of iron is considerably less than the ice (see Table 1). Curve (d) is a three component mixture with the addition of troilite to ice and rock. Troilite also has a fairly large refractive index at this wavelength, which affects the linear approximation significantly. Curve (e) is a five component mixture with the contents of (d) augmented by iron and organics. The Garnett dielectric constant is unaffected by the minute amount of iron, as has been found by others (Ossenkopf 1991), but the linear approximation diverges considerably because of iron's very high dielectric constant. This shows that metals, in general, are well beyond the realm of validity for the linear approximation.

**Figure 12.** Plots of the real (left) and imaginary (right) components of average refractive index of a material of density ρ as a function of the density ratio $\rho /\overline{\rho }$ . Curves marked (a) correspond to pure ice. Curves (b and c) correspond to two component mixtures of ice and rock (b) and ice and iron (c). Curve (d) corresponds to a three component mixtures of ice, rock, and troilite, while curve (e) contains these three as well as organics and iron. Note that the linear and Garnett correlate very poorly for mixtures that contain even small mixing ratios of elements with very high refractive indices. This message is also conveyed by Figure 7 and the associated discussion.
Download figure:
Standard image High-resolution image

The fractional mass and volume factors α_j and f_j discussed above are simply related to the particle number density weighting factors β_j of Section 4 (e.g., Equation (16)ff). Recall that n(r) is the total number density of particles of all species, and if all the species have the same size distribution, the number density of particles of species j is n_j(r) = β_jn(r), where ∑β_j = 1. In the limit where the scattering and absorption by each species needs to be separately treated (physically, for particles small enough to be monomineralic), ρ_j(r) = ρ_j = constant. Then the total nebular mass density in species j is

$\begin{equation} \rho _{pj} = \int \frac{4}{3} \pi r^3 \rho _j \beta _j n(r) dr = \rho _j \beta _j \int \frac{4}{3} \pi r^3 n(r) dr = \rho _j \beta _j \xi, \end{equation} \tag{ C14 }$

where ξ is the nebula volume fraction of all solids. Recall that ρ_p = ∑ρ_pj = αρ_g, where α is the total mass fraction solid material. Therefore,

$\begin{equation} \beta _j = {\rho _{pj} \over \xi \rho _j }= {\alpha _j \rho _g \over \xi \rho _j} = {\alpha _j \rho _g \overline{\rho } \over \rho _p \rho _j}, \end{equation} \tag{ C15 }$

where the last equality follows from $\rho _p \equiv \alpha \rho _g = \overline{\rho } \xi$ , since for the heterogeneous particle case we assume the particle porosity ϕ = 0. We previously defined the volume-averaged mass density of a compact particle as Equation (C13):

$\begin{equation} \overline{\rho }=\frac{\alpha }{\sum _{j} (\alpha _{j}/\rho _{j})} = \frac{\rho _p}{\xi }; \end{equation} \tag{ C16 }$

then

$\begin{eqnarray} \beta _j &=& \frac{\alpha _j}{\rho _j} \frac{\rho _g}{\xi } = \frac{\alpha _j}{\rho _j} \frac{\rho _p}{\xi \alpha } = \frac{\alpha _j}{\rho _j} \frac{\overline{\rho }}{\alpha } = \frac{\alpha _j}{\alpha } \frac{\overline{\rho }}{\rho _j}\nonumber\\ & = & {\alpha _j \over \rho _j \sum _j (\alpha _j / \rho _j) } = \frac{f_j}{1-\phi }, \end{eqnarray} \tag{ C17 }$

so it is simply verified that ∑_jβ_j = 1. Moreover, Equations (C13) and (C17) can be combined to give $\alpha _j \overline{\rho } = \alpha \beta _j \rho _j$ , and summing both sides over j leads to the intuitively obvious alternative expression $\overline{\rho } = \sum _j \beta _j \rho _j$ . All of the above quantities are obtained from the presumed compositional makeup of the protoplanetary nebula.

UTILITARIAN OPACITY MODEL FOR AGGREGATE PARTICLES IN PROTOPLANETARY NEBULAE AND EXOPLANET ATMOSPHERES

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. PRIOR WORK AND OVERVIEW OF PAPER

3. OPACITY MODEL

3.1. Extinction, Absorption, and Scattering Efficiencies