Free Access
Issue
A&A
Volume 572, December 2014
Article Number A118
Number of page(s) 34
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201423702
Published online 08 December 2014

© ESO, 2014

1. Introduction

Thanks to the precise determination of the mass and radius of many extrasolar planets, it has recently become possible to study the bulk composition of a rapidly growing number of exoplanets. These planets have masses ranging from super-Earth to Jovian masses. This allows us to investigate statistically how the bulk composition depends on fundamental planetary parameters like the mass or the semi-major axis, and how this compares to the solar system with its basic types of planets (terrestrial, ice giants, and gas giants). A specific question that has recently attracted a lot of attention is the transition from solid to gas dominated planets, or in other words, how the hydrogen/helium mass fraction depends on a planet’s properties (e.g., Miller & Fortney 2011; Weiss et al. 2013; Lopez & Fortney 2014; Marcy et al. 2014).

Planet formation theory based on the core accretion paradigm (Perri & Cameron 1974; Mizuno 1980; Stevenson 1982; Bodenheimer & Pollack 1986) predicts that the H/He mass fraction is an increasing function of the planet’s mass (see Rogers et al. 2011; and Mordasini et al. 2012b). This is because the Kelvin-Helmholtz timescale for the contraction of the envelope which controls the gas accretion rate in the subcritical regime is a decreasing function of the core mass (Ikoma et al. 2000). Population synthesis simulations based on global core accretion models have therefore found a synthetic mass-radius relationship (or in other words, a bulk composition) that resembles the observed relationship (Mordasini et al. 2012b). These calculations also predict that the population of planets with radii larger than ~2 R can be characterized by planets possessing a H/He envelope, whereas other planetary types dominate at smaller radii. Similar conclusions were reached from observational data (Howard et al. 2012; Gaidos et al. 2012; Wolfgang & Laughlin 2012; Marcy et al. 2014).

The Kelvin-Helmholtz timescale is, however, not only a function of the core mass, but also of the opacity κ in the protoplanetary envelope (e.g., Ikoma et al. 2000)1. The latter is thought to be mainly due to tiny dust grains suspended in the protoplanet’s outer radiative zone. The lower the opacity, the shorter the cooling timescale, so that for a fixed core mass a more massive envelope can be accreted during the lifetime of the protoplanetary gas nebula. This eventually translates into different planetary mass-radius relationships, because at a given total mass, the H/He mass fraction will be higher if κ was low during formation which causes a larger planetary radius during the evolutionary phase after the dispersion of the nebula (e.g., Fortney et al. 2007; Valencia et al. 2013). Interestingly, it means that it is, at least in principle, possible to relate an observable quantity like the mass-radius relationship to the microphysics of grain growth during formation (for other effects, see Vazan et al. 2013). This link was studied in Mordasini et al. (2014, hereafter Paper I, with the conclusion that in order to reproduce the observed planetary bulk compositions, the grain opacity κgr in protoplanetary atmospheres must be much smaller than in the interstellar medium (ISM).

The magnitude of the grain opacity is essential not only for determining the H/He content of low-mass planets, it is also a crucial factor in giant planet formation because the critical core mass beyond which rapid runaway gas accretion occurs, depends on the opacity2. For example, as already found in the early work of Mizuno (1980), the critical core mass can vary from 1.5 M at grain-free opacity to 12 M at full ISM opacity. This is a very significant difference.

The potentially very important consequences of different grain opacities outlined above motivated Podolak (2003) to develop a detailed but numerically expensive model for the microphysics of the grains in protoplanetary atmospheres that yields a physically motivated κgr. It numerically solves the Smoluchowski equation in each atmospheric layer, taking into account the effects of grain growth, settling, and vaporization. This model was further elaborated in Movshovitz & Podolak (2008, hereafter MP08, and was applied to an atmospheric structure of a protoplanet calculated by Hubickyj et al. (2005). In Movshovitz et al. (2010, hereafter MBPL10, and later Rogers et al. (2011) this grain evolution model was finally self-consistently coupled to the giant planet formation model of Bodenheimer and collaborators. In particular, in MBPL10 simulations of the in situ formation of Jupiter as in Pollack et al. (1996) were presented, but using a physically motivated grain opacity instead of an arbitrarily scaled ISM opacity.

These calculations showed that first, the conditions in protoplanetary atmospheres are such that the aforementioned processes may substantially alter the properties of the grains. Second, it was found that the grain opacity and resulting optical depth of the atmosphere can be substantially reduced relative to the ISM case. This results in a short formation timescale of giant planets, and allows low-mass cores of only about 4 M to trigger gas runaway accretion during the typical lifetime of a protoplanetary disk. For comparison, at full ISM opacity, this would take several hundred Myr (Paper I). Finally, the radial structure of the opacity in the protoplanet’s atmosphere is characterized by a high, ISM-like opacity in the outer layers (κgr ~ 1 cm2/g) that falls to very low values (κgr ~ 10-3 cm2/g) in the deep layers close to the radiative-convective boundary. This substantially differs from the radial opacity structure that is obtained when simply scaling the ISM opacity (see Sect. 3.1.5).

Despite this, with the exception of the aforementioned works, it has been customary in the literature to model the effect of a lower grain opacity by simply reducing the ISM opacity by some arbitrary uniform reduction factor fopa. This is the case for very different kinds of studies ranging from analytical work to 3D hydrodynamic simulations (e.g., Mizuno et al. 1978; Stevenson 1982; Pollack et al. 1996; Papaloizou & Nelson 2005; Hubickyj et al. 2005; Tanigawa & Ohtsuki 2010; Levison et al. 2010; Hori & Ikoma 2011; Ayliffe & Bate 2012). Because of the lack of knowledge of even just the magnitude of the grain opacity κgr, very large ranges for fopa were often considered, typically from 10-3 to 1. The same basic approach was also taken in Paper I, with the improvement that fopa was calibrated by determining which reduction factor leads to the same formation timescales for Jupiter as in MBPL10. A best-fitting value of fopa ≈ 0.003 was found in this way, and used as the nominal value in population synthesis calculations. It is clear, however, that one global reduction factor that was only calibrated for certain conditions like a semi-major axis, core mass, or pressure and temperature in the protoplanetary disk is, in principle, not applicable to an entire population of planets that cover a very wide range of planetary properties.

This shows that it is desirable to develop a better analytical understanding of the mechanisms governing the grain evolution in protoplanetary atmospheres, and to be able to estimate the expected magnitude of the grain opacity. The goal of this work is therefore to derive a first (and simple) analytical model for the grain opacity based on microphysical processes. Despite the simplicity, this model should be able to predict κgr in such a way that it also quantitatively agrees with the results of the numerical model of MP08 and MBPL10. On one hand, this fosters the physical understanding of how important parameters like the planet’s mass or the dust-to-gas ratio in the disk (which could be related to the stellar [Fe/H]) influence the grain opacity and therefore the growth history of a protoplanet. On the other hand, from a practical point of view the model makes it possible to efficiently calculate a physically motivated κgr because the expressions for κgr derived below can be easily integrated in the calculation of the standard internal structure equations of a protoplanetary envelope without requiring special, time consuming measures like substepping or iterations. This in turn means that this model can be used in numerical simulations of the accretion of gas by protoplanets in general, and in population synthesis calculations in particular. In a contemporaneous manuscript, Ormel (2014) has presented a model of intermediate complexity. Similar to the analytical model presented here, it assumes that per layer there are only one or two typical grain sizes. Similar to the numerical model, the grain size is found numerically by integrating a fifth equation in addition to the normal planetary structure equations. This in particular makes it possible to take into account the effect of a bimodal grain size distribution and of a radially varying grain input from planetesimal ablation.

In future work, we will therefore run population syntheses using the analytical grain opacity model and the work of Ormel (2014), and compare these with the results of Paper I. We will also take into account effects like the concurrent formation of several embryos (Alibert et al. 2013), the effect of envelope enrichment by heavy elements (Fortney et al. 2013), or the evaporation of the primordial H/He envelope of close-in planets (Jin et al. 2014). This will help to better understand the role of the grain opacity in shaping the observed mass-radius relationship and thus bulk composition of the extrasolar planets.

The contents of this paper are as follows. In Sect. 2 we first write down the general equations describing the dynamics of the grains due to dust settling, growth (by coagulation due to Brownian motion or differential settling), evaporation, and due to the contraction of the gaseous envelope itself. In Sect. 2.2.2 we discuss the two mechanisms that bring new grains into the protoplanetary atmosphere, which are the accretion together with newly accreted nebular gas, and second the breakup of planetesimals flying through the atmosphere. We then derive analytical expressions for κgr in Sect. 2.5, treating the five relevant regimes (Epstein/Stokes drag, Brownian coagulation/differential settling, and grain advection) separately. In Sect. 3, the analytical model is applied to calculate the grain opacity in the same atmospheric structure of Hubickyj et al. (2005) that was already considered in MP08. A detailed analysis of the processes regulating the growth of the grains is made, as well as different comparisons of the results of the analytical and numerical model. In Sect. 4, we show the results of coupling the analytical grain opacity model with our core accretion model. This section corresponds to the work of MBPL10, and extensive comparisons are made between the analytical and numerical result. Finally, in Sect. 5 we summarize our results and present our conclusions.

2. Analytical model

The basic approach of the analytical model is to calculate for each layer the typical size of the grains by comparing the timescales of the governing microphysical processes. Based on the typical size of the grains, and the (supposed) conservation of the radial flux of the dust grains across the atmosphere (due to settling), it is then possible to calculate the grain opacity in each layer. This fundamental approach is very similar to the one used by Rossow (1978) who calculated the microphysics of cloud formation in the atmospheres of different planets in the solar system. This approach was later used by Cooper et al. (2003) to study cloud formation and the resulting opacity in atmospheres of brown dwarfs. In the context of grain opacity in protoplanetary atmospheres, estimations based on timescale arguments were also made in Podolak (2004); Movshovitz et al. (2010); Nayakshin (2010); and Helled & Bodenheimer (2011). In the context of grain growth in protoplanetary disks, timescale arguments were used in Birnstiel et al. (2012), for example. Here we develop a simple, but complete model to dynamically calculate the grain opacity in protoplanetary atmospheres as a function of local atmospheric properties. It can be used in numerical simulations of (giant) planet formation (Sect. 4).

2.1. Relevant timescales

For the calculation of the grain size, the most important timescales are the settling and growth timescale. Additional timescales that must be considered are the advection timescale of the grains (because the gas envelope is itself not static) and the evaporation timescale. While the specific form of the settling and growth timescale depends on the drag regime (Epstein and Stokes regime) and the growth mechanism (Brownian motion coagulation or differential settling), they share the property that the growth timescale increases with increasing grain size, while the settling timescale decreases with increasing grain size. This means that in a given layer, very small grains will only be found in small quantities, because they quickly grow to larger sizes. Very large grains will also be present only in small quantities because they quickly fall out of the layer. The typical size is therefore found by equating the two timescales.

In this section, we first define the general equations necessary for the analytical model, including the different aerodynamic regimes, the expression for the grain accretion rate, and the general expressions for the opacity. We then calculate the typical grain size for four regimes: Brownian coagulation in the Epstein regime, Brownian coagulation in the Stokes regime, differential settling in the Epstein regime, differential settling in the Stokes regime. And finally, we study the fifth regime, the advection regime, where the radial motion due to the contraction of the gas envelope is dominant over the settling of the grains. Once the grain size (and abundance) is known, the opacity can be calculated. We then address the effect of grain evaporation, give an approximation for the extinction coefficient, and show how the different regimes are coupled to obtain the final expression for κ.

2.2. General equations

The settling timescale is the time it takes the grains to cross a typical length scale in the atmosphere if they are settling at a velocity vset relative to the gas. The natural length scale in an atmosphere is the scale height H, so that the settling timescale τset can be estimated in the limit of a static gas envelope as τset = H/vset. The scale height in the atmosphere is (1)where kB is the Boltzmann constant, T the temperature of the gas, μ its mean molecular weight (2.4 for solar composition and molecular hydrogen), mH the mass of a hydrogen atom, and g is the gravitational acceleration that is given as g = GM(R) /R2. In this equation, G is the gravitational constant, R the radial distance measured from the center of the planet (hereafter called the “height”), and M(R) the mass within R (core plus envelope gas). The structure of a protoplanetary envelope usually consists of a deep convective zone and an upper radiative zone, even though more complicated structures with several convective zones occur (see, e.g., Bodenheimer & Pollack 1986,or Fig. 11). Our model for the grain opacity applies to the upper radiative zone that we call (as in MBPL10) the atmosphere.

Typically, the mass in the radiative zone is much lower than the mass of the core and the envelope mass in the convective zone, so that M(R) is nearly constant. However, a constant M(R) is not required for the analytical model for κgr.

2.2.1. Aerodynamic regimes

The settling and growth of the grains depend on the aerodynamic regime they are in. Two dimensionless numbers characterize the aerodynamic regimes. First, the Knudsen number Kn describes whether the grains are in the free molecular flow regime where the grains interact with single molecules (large Kn), or if the gas acts as a continuous, hydrodynamic fluid (small Kn). The Knudsen number is given as Kn = /a, where a is the radius of a grain and is the mean free path of a gas molecule that is calculated as (2)In this expression, ρ is the gas density and σm the collisional cross section of a molecule, . For simplicity, we set the molecular diameter dm equal to that of molecular hydrogen H2 at standard conditions 2.7 × 10-8 cm (Waldmann 1958). It is clear that, in reality, dm would depend on both the composition of the gas as well as its pressure and temperature. We set the critical Knudsen number for the transition from free molecular to continuum flow to Kn = 4/9 (e.g., Stepinski & Valageas 1996).

The second dimensionless number is the Reynolds number. It describes whether the flow is laminar (small Re) or turbulent (large Re) and is given as Re = 2aρv/η, where v is the settling velocity of the grains and η is the dynamic viscosity of the gas. It is given in the approximation of hard-sphere molecules with a Maxwellian velocity distribution as ρvth/ 3, where vth is the mean thermal velocity of the gas, (3)More accurate expressions for η based on, e.g., the Lennard-Jones potential instead of hard spheres exist (see, e.g., Podolak 2003), but given the other simplifications in the model, we stick to this basic expression.

We find below that the grains are typically in the free molecular flow regime in the outer parts of the atmosphere, and in the continuum flow in the inner parts. The Reynolds number on the other hand is always small (typically much less than unity) throughout the atmosphere, so that the grains are always in the laminar flow regime (see Fig. 5). This justifies the choice of drag laws used below.

In the derivation below, we will also assume that the grains always settle at the terminal velocity where the drag and gravitational force balance each other. This holds for particles with a short stopping timescale and can be quantified by a third dimensionless number which is the ratio of the stopping time τstop of the particles to the characteristic flow time3. In the current situation, the background motion is the grain settling, therefore we define RSS = τstop/τset where a small RSS means that the terminal velocity is rapidly reached. In this expression, the stopping time is τstop = mgrv/FD, where FD is the drag force and mgr is the mass of a dust grain. As in MP08 we assume for simplicity spherical dust grains (which might not be justified, see Sect. 2.9) so that mgr = 4/3πρgra3, where ρgr is the material density of the grains. Like MP08 we also only consider silicate grains in this work, and follow them in setting the material density of the grains to ρgr = 2.8 g/cm3.

We find below (Sect. 3.1.2) that the grains in protoplanetary atmospheres have very small RSS ~10-7, and very short stopping times of a few seconds to minutes. The assumption of settling at the terminal velocity in the settling regime or of a motion at the velocity of the gas in the advection regime (Sect. 2.5.5) is therefore justified.

2.2.2. Planetesimal mass deposition, grain accretion rate, and dust-to-gas ratio

In this section we specify the rate at which grains are accreted into the protoplanetary atmosphere, gr. Two different sources exist: first, grains that are brought into the atmosphere together with the newly accreted nebular gas and second, grains that are due to the accretion of planetesimals. During their flight through the envelope, planetesimals undergo mass loss via thermal ablation. Large impactors are additionally aerodynamically disrupted by large dynamic pressures (see Podolak et al. 1988 and Mordasini et al. 2006). The aerodynamic disruption breaks big impactors into small fragments that are then easily digested by thermal ablation since fragmentation greatly increases the ablating surface.

Grains that are accreted by the gas are injected into the planet’s envelope at its outer radius. Grains originating from planetesimal ablation can, in contrast, be injected into the envelope at different heights depending on the aforementioned mechanisms that determine the radial mass deposition profile of the planetesimals. This mass deposition profile can be complex, and depends on the planetesimal size, the impact velocity, the tensile strength of the body, the impact geometry, and so on. In this work, we consider two idealized, limiting cases.

First, in the “deep deposition” case, it is assumed that the planetesimals fly through the atmosphere (outer radiative zone) without losing any mass. The planetesimals deposit their mass only in the deep parts of the envelope. If the mass deposition occurs at a radius where the grain-free gas opacity dominates over the grain opacity in all cases, and/or where the temperature in the background envelope is higher than the evaporation temperature of the grains, then the grains brought in by the planetesimals do not influence the opacity and therefore the atmospheric structure. The second depends on the opacity; therefore, the mass deposition profile and atmospheric profile are interdependent in both directions. The only source of grains in the deep deposition case is then the accretion of gas that occurs at a rate XY. Grains are mixed into this gas at a dust-to-gas mass ratio in the protoplanetary disk, fD/G,disk. Typically fD/G,disk ≈ 0.01, but this ratio could also be much lower if most solids have already been incorporated into large bodies. The dust-to-gas ratio in the disk could also scale with the stellar [Fe/H], establishing an interesting link between opacity (and in the end the planetary H/He content, i.e., bulk composition) and stellar properties. We discuss the impact of [Fe/H] on κgr in Sect. 3.1.7.

Second, in the “shallow deposition” case, it is instead assumed that the planetesimals break up into grains at the very top of the atmosphere, so that the full accretion rate of planetesimals Z contributes to the grain accretion rate. In equations, the two cases are (4)The larger the planetesimals, the deeper they usually deposit their mass, even though the actual behavior is complex because of different mechanisms leading to mass loss (pure thermal ablation versus aerodynamic disruption; see Mordasini et al. 2006). Still, as a first approximation, the deep (shallow) deposition case can be associated with the accretion of big (small) planetesimals. This is illustrated in Fig. 1. The figure shows the radial mass deposition profile (ablated mass per unit length) of planetesimals of different sizes as a function of height in the protoplanetary atmosphere in the MBPL10 comparison calculation (see Sect. 4.3.2 below for a description). The profiles were obtained by simulating the impact of a planetesimal with an impact model similar to Podolak et al. (1988), but using an updated model for the heat transfer coefficient, and a multi-staged aerodynamic fragmentation model (see Mordasini et al. 2006 and Alibert et al. 2005 for a short description; a full model description will be given in a future work).

thumbnail Fig. 1

Mass deposition profile (ablated mass per unit length) as a function of height for 100 km planetesimals (solid blue line). The brown dash-dotted line is the mass deposition for 1 km planetesimals, multiplied by a factor 106, while the green dashed line shows the result for 10 m planetesimals, multiplied by a factor 1012, so that the total mass deposited is identical in the three cases. The gray vertical lines show the radius where the envelope temperature becomes so high that silicate grains evaporate (“evap”) and where the molecular opacities alone become larger than the opacity due to dust in the deep deposition case (“κ-DG”). The steps in the right part of the lines are a numerical artifact without physical meaning.

Open with DEXTER

Figure 1 shows the mass deposition by a planetesimal with an initial size of 100 km, 1 km, and 10 m (see also Appendices A and B). The curves for 1 km and 10 m were scaled so that the total deposited mass is the same in the three cases (on the order of 1021 g). The planetesimals hit head-on with an initial velocity approximately equal to the planet’s escape velocity. Material parameters appropriate for water ice were used which is, strictly speaking, inconsistent with the assumption of silicate grains made above. The consequences of the planetesimal composition on the mass deposition profile is studied in Appendix B. It is found that rocky planetesimals deposit their mass in deeper layers, as expected. However, the extent of the difference between icy and rocky impactors depends significantly on the impactor size, with a larger difference for small planetesimals.

One notes that the bigger the planetesimal, the deeper the peak mass deposition, as expected. For the 100 and 1 km sized planetesimals, the mass deposition is very low in the outer parts since pure thermal ablation is not efficient for large bodies (Svetsov et al. 1995), but then shoots up by many orders of magnitudes as soon as the impactor is aerodynamically fragmented. This mass (and energy) deposition concentrated over a comparatively small radial domain (0.05 R in the 100 km case) is reminiscent of the terminal explosion that was seen in the simulations of the collision of comet SL-9 with Jupiter (e.g., Mac Low & Zahnle 1994) even if the absolute scales are much longer here. The mass deposition profile of the 10 m planetesimal is, in contrast, much smoother because no violent fragmentation occurs. Much more mass is deposited in the outer layers. Figure 1 also contains two vertical lines. The right line shows where the grain-free molecular opacities for a solar composition gas becomes dominant over the grain opacity under the assumption of deep deposition (Sects. 3.1.1 and 4.3.1). The left line corresponds to the radius where the temperature in the undisturbed envelope is sufficiently high to evaporate silicate grains on a short timescale (see Sect. 2.6). The radiative-convective boundary is found at nearly the same radius, meaning that most of the mass deposited by 100 and 1 km planetesimals goes in vaporized form into the deep convective zone without contributing much to the grain opacity in the outer atmosphere. The 10 m planetesimals deposit a significant fraction of their mass in the atmosphere, contributing in an important way to the grain influx into the atmosphere. It is clear that in reality, most impacts will not be head-on as assumed here. The consequences of different impact geometries are studied in Appendix A. As expected, off-center collisions typically increase the mass deposition in the upper layers. However, for 100 km planetesimals, the vast majority of the mass is still deposited in the terminal explosion the height of which varies by ~1 R depending on the geometry.

The plot also indicates that the two extreme cases for gr given by Eq. (4) are useful limiting cases bracketing the actual behavior, even if they are highly simplified representations of the actual mass deposition profile. Since we assume that the planetesimals have a size of 100 km in the nominal model, we will consider deep deposition as the nominal case, but we will also investigate the “shallow” case, which should be more appropriate for the accretion of small bodies (pebbles, see Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012).

We will see below (Sects. 2.5.3 and 2.5.4) that somewhat surprisingly, in the most important grain growth regime, the assumption of a deep or shallow deposition is not important. We will even find that the grain opacity is generally independent of gr for the dominant range of conditions. This is clearly the result of this study that has the most important implications for general planet formation theory (Sect. 3.1.7).

Since the grain opacity is important in the outer radiative zone that typically contains only a very small fraction of the total envelope mass (e.g., Bodenheimer & Pollack 1986), we assume that grains do not accumulate locally in the atmosphere, but that all grains that enter the atmosphere at its outer edge will also settle out of it at the radiative-convective boundary. As demonstrated by MP08, we can also assume that the processes governing grain evolution (growth, settling, evaporation) are so fast that the grains almost instantaneously assume a steady state. This is justified by the fact that these processes occur on timescales of just 1–100 years (see Fig. 4), which is typically much shorter than the timescales on which the properties of the protoplanetary envelope itself change.

We can then assume that the instantaneous grain accretion rate is radially constant in the atmosphere, i.e., dgr/ dR = 0. Mass conservation then directly yields the dust-to-gas ratio fD/G (by mass) in an atmospheric layer as (5)where vset is the vertical settling velocity of the grains. Because both ρ and vset vary with R, fD/G will, in general, also be a function of the radius inside the atmosphere and will differ from the value at the top fD/G,disk, except for the advection regime (Sect. 2.5.5).

2.2.3. Basic expression for opacity

Equating the growth and settling time yields the typical grain size in each layer. Once this grain size a is found, the settling velocity can also be calculated so that, together with Eq. (5) and the assumption of spherical grains, one can calculate the number density of the grains ngr = fD/Gρ/mgr. The contribution of grains of size a to the wavelength dependent grain opacity is then (MP08) (6)where Q(x) is the scattering efficiency. We note that in this work, with grain opacity we always mean the opacity per unit mass of envelope material (gas and grains), and not per unit mass of grains only. The scattering efficiency depends on grain material properties and on the size parameter x = 2πa/λ where λ is the wavelength of the impinging photons. The calculation of Q is described below in Sect. 2.7. In the limit of large grains relative to the typical wavelength (x ≫ 1), Q(x) is simply equal to 2. This regime is found below to be the most important one.

In our simple analytical model, there is only one grain size per layer. In addition, if the wavelength dependency of Q over the relevant wavelength domain is weak, then κ given by Eq. (6) can be directly identified with the Rosseland mean opacity which is the quantity needed to solve the internal structure equations. With (7)we can also write (8)This shows that for a nearly wavelength independent Q only the ratio of the dust-to-gas ratio to the typical particle size enters the opacity.

2.3. Drag regimes

We now give the equations for the dust-to-gas ratio and the settling velocity in the Epstein and Stokes drag regimes (see, e.g., Weidenschilling 1977; Stepinski & Valageas 1996), assuming that the Reynolds number is always small (Sect. 2.2.1). We treat the two regimes separately, but we note that it is in principle possible to write the drag law in a combined form that approaches the two separate expressions in the limit of large or small Knudsen numbers (Nayakshin 2010). The reason for treating the regimes separately is that it is then possible to solve analytically for a. The combined case is treated in Appendix C.

We find that the impact of treating the two regimes separately has no negative consequences, because the transition behaves in a benign way even if, a priori, it is not known which regime applies to a given atmospheric level. Operationally, one therefore derives the grain size for both drag regimes, and then checks a posteriori which regime is self-consistent (i.e., if the a priori assumed Kn regime is the one that is actually found a posteriori for the resulting grain size). One finds that not only is the resulting grain radius continuous across the transition between the two regimes at Kn = 4/9 (as it must be), but also that the overlapping radial domain where both the Epstein and Stokes regime lead to self-consistent results (so that it is unclear which regime actually applies) is small.

2.3.1. Epstein: Kn 4/9 (small grains)

In the Epstein regime which applies to small grains in the upper part of the atmosphere, the drag force is given as (9)This leads to a stopping time of (10)The terminal settling velocity found by equating the drag and gravitational force in the Epstein regime and the corresponding settling timescale through an atmospheric scale height H is (11)Together with Eq. (5), the dust-to-gas ratio is (12)

2.3.2. Stokes: Kn < 4/9 (large grains)

The Stokes regime applies to grains that are large in comparison to the free mean path of the gas molecules. The drag force in this regime is (13)leading to a stopping time of (14)The terminal settling velocity and the timescale to cross one scale height is (15)Finally, the dust-to-gas ratio is given as (16)

2.4. Brownian coagulation

As do MP08, we consider two processes that lead to grain growth in protoplanetary atmospheres: first, coagulation due to the Brownian motion of the grains, and second, growth due to the differential settling speed of grains of different sizes. The effects of convection that were included by Podolak (2003) and MBPL10 (and of turbulence in the atmosphere in general) are neglected since MBPL10 find that including or neglecting convection in the grain calculations has a negligible effect.

The growth timescale due to Brownian motion coagulation4 can be derived (e.g., Helled & Bodenheimer 2011) by writing the coagulation timescale τcoag as τcoag = gr/Vth,gr, where gr is the mean free path of a grain between collisions with another grain, while Vth,gr is its mean thermal velocity. The mean free path is calculated as (17)where σgr = 4πa2 is the collisional cross section. Assuming that the temperature of the grains and gas is identical, we also have (18)Combining these equations gives the coagulation timescale for Brownian motion, (19)which is identical to the expression in Helled & Bodenheimer (2011). The coagulation timescale thus decreases with an increasing dust-to-gas ratio as 1 /fD/G and increases with grain size as a5/2. We also see that fluffy grains with a low material density ρgr coagulate faster, as expected.

2.5. Calculation of the typical grain size and opacity

With the expressions for the settling and growth timescale at hand, we can proceed to the central part of the analytical model which is the calculation of the typical grain size a. We consider five different regimes separately where we always equate the appropriate growth and settling timescales: Brownian coagulation in the Epstein regime, Brownian coagulation in the Stokes regime, differential settling in the Epstein regime, differential settling in the Stokes regime, and, as a special case, Brownian coagulation and growth by differential settling in the grain advection regime which is due to the contraction of the gaseous envelope.

In the advection regime it is sometimes found that the grain size resulting from equating the growth and advection timescale is so small that it is smaller than the radius thought to be typical for monomers, amono. In such a case (and also if this should occur in any other regime), we replace the grain radius found from the timescale argument by the assumed monomer size. As MP08, we consider monomer sizes of 1, 10, and 100 μm (see Sect. 3.1.8). It is, however, found that this condition rarely happens, and that amono therefore only has a small effect on the global evolution of the planets (Sect. 4.2.3), provided that the monomer size is not very large (more than ~100 μm). Very large grains could in contrast lead to a significant reduction of the formation timescale.

2.5.1. Brownian coagulation: Epstein regime

In the first regime, we insert Eq. (12) into Eq. (19) to find a coagulation timescale of (20)Setting this timescale equal to the settling timescale in the Epstein regime (Eq. (11)) τset,E = τcoag,E, and solving for the grain radius aC,E (coagulation, Epstein) yields (21)This expression is identical to the one derived in MBPL10 if their typical length scale is associated with the atmospheric scale height. Inserting into Eq. (12) gives the dust-to-gas ratio (22)The opacity is finally found by inserting aC,E and fD/G,C,E into Eq. (8), yielding (23)or, in terms of the fundamental local properties of the atmospheric layer (and gr), (24)We see that κC,E increases as . The dependency on the position inside the atmosphere R and on M, and therefore the planet’s core or total mass, is not straightforward to see because T and ρ also depend on these quantities.

2.5.2. Brownian coagulation: Stokes regime

In the second regime, we insert Eq. (16) into Eq. (19) to find a coagulation timescale of (25)Setting this equal to the settling timescale in the Stokes regime (Eq. (15)), τset,S = τcoag,S, and solving for the grain radius, we find (26)Inserting back into Eq. (16) yields the dust-to-gas ratio (27)and the opacity (28)which we can also write as (29)In this regime we have κC,E increasing as . The dependency on R and M is again not obvious because of the interdependency with T and ρ.

2.5.3. Differential settling (coalescence): Epstein regime

Larger grains can sweep up smaller ones because of their higher sedimentation speed, leading to growth by differential settling. This process is called coalescence in Rossow (1978). In protoplanetary disks, differential settling is the dominant growth mode for larger grains, and means that the grains grow more quickly than with Brownian coagulation alone (e.g., Brauer et al. 2008). We will find below that in protoplanetary atmospheres, differential settling is also the dominant growth process. To model this growth mode, we use (as did Cooper et al. 2003) the growth timescales determined by Rossow (1978) who generalizes the results for grains of different sizes to determine the growth timescale when (formally) only one grain size is present per layer.

For small grains in the Epstein drag regime (Kn > 4/9), the timescale for growth by differential settling is (30)When using the expression for vth, ngr, fD/G, and the expression for the settling velocity in the Epstein regime, this corresponds to(31)One thus finds that τcoal,E ∝ 1 /fD/Ga/gr. Equating this with the settling timescale in the Epstein regime τcoal,E = τset,E gives a typical grain size of (32)The dust-to-gas ratio is then (33)Finally, the opacity is (34)The opacity in this regime has thus a simple functional form and the remarkable property that it is independent of gr. The result for κgr is therefore, in particular, also the same in the deep or shallow deposition regime for planetesimal impacts (Sect. 2.2.2). The lack of a dependence on gr is a consequence of the fact that κgrfD/G/a, and that in this regime both fD/G and a scale in the same way with . Physically it means that if we increase gr, we increase the dust-to-gas ratio which would in principle increase the opacity, but the increased dust-to-gas ratio also leads to the formation of larger grains, and this decreases the opacity. In the current regime, the two effects compensate for each other. It is found that this regime (in the upper layers) and differential settling in the Stokes regime (in the lower layers) are the dominant ones for the opacity in protoplanetary atmospheres (see Fig. 4). At the same time, the grains are typically large in comparison to the wavelength, so that Q ≈ 2. Therefore, in the outer parts of a protoplanetary atmosphere, the opacity is simply 27/(4).

The finding that κgr is independent of gr can of course not hold ad infinitum. For example, if gr approaches zero, then the grain opacity must also approach zero, at least in our limit that the steady state is assumed instantaneously. What actually happens is the following: if gr becomes very small, then τcoal,E becomes long, to an extent that it becomes longer than the advection timescale of the gas. Therefore, the growth regimes changes from the differential settling regime into the advection regime, where there is a dependency on gr. We will study the impact of gr (which can be studied by varying fD/G,disk, because of Eq. (4)) in Sect. 3.1.6.

With the definition of H and g, and the ideal gas law5, we can also write (35)where P is the atmospheric pressure. This form is reminiscent of Eddington’s photospheric boundary condition, where κ = 2g/ (3P).

2.5.4. Differential settling (coalescence): Stokes regime

For large particles with Kn < 4/9 at low velocities (Re 70), the timescale for growth by differential settling is given as (Rossow 1978) (36)With the definitions of ngr and fD/G, this can be written as (37)which has the same dependency on fD/G, a, and gr as in the last regime. Equating this timescale with the settling timescale for Stokes drag, τset,S = τcoal,S, gives a typical grain radius (38)The dust-to-gas ratio is (39)which finally leads (with Eq. (8)) to an opacity (40)This is, except for a factor on the order of order unity (1.69), the same result as for differential settling in the Epstein regime. It is therefore a very simple expression that is in particular independent of the grain accretion rate, . The reason is the same as before, except that in this regime, fD/G and a both scale as instead of . Using again the fundamental variables describing the atmosphere, we can also write the opacity as (41)This equation and the corresponding one for the Epstein regime (Eq. (35)) are the most important results of this study.

2.5.5. Advection regime: effect of envelope contraction

Up to this point, we have completely neglected that the atmosphere is itself not static. In general, this is a good approximation, because the timescale on which the atmosphere changes is much longer than the timescale on which the grains evolve.

However, under some circumstances, namely at or after crossover, it is necessary to take into account that in reality, the gaseous envelope of the planet is continuously contracting. This means that the gas itself will also have a (small) downward motion through the atmosphere because the gas that is newly accreted at the outer radius settles down into the deep parts of the envelope where the mass accumulates (in the convective zone). This leads to a velocity of the gas vgas that must be taken into account to describe the net velocity of the grains relative to the Eulerian atmospheric p-T structure at a fixed radius (MP08).

Typically, only a very small gas mass is contained in the radiative zone, so we can assume that the gas accretion rate, XY, is radially constant in the atmosphere; of course, this does not hold in the inner parts of the envelope where the gas accumulates. The vertical velocity of the gas can then be estimated as(42)and the associated timescale is simply τadv = H/vgas. We have called this timescale the advection timescale for the following reason: In the limit that vgas becomes large, the gas starts to entrain the grains, so that the grains get externally advected through a scale height by the gas flow. This holds because the grains are well coupled to the gas as is made clear by the small stopping times (see, e.g., Fig. 4).

In principle it is possible to take the effect of the gas flow into account exactly by writing the drag force proportional to the relative velocity. At equilibrium, when the gravitational force FG and drag force FD balance each other, we again have FG + FD = 0, but FD now has the form Ccouple(vvgas), where the coupling constant Ccouple depends on the drag regime (Eqs. (9) and (13)). Solving for the velocity of the grains, we have v = (FG + Ccouplevgas) /Ccouple. The first regime occurs when the second term is negligibly small (this is the previously treated static case), while the second regime occurs when the first term is negligible, so that v = vgas, i.e., when the grains are advected by the gas flow. In this first analytical model, we take the effect of envelope contraction into account only in the case of complete advection where the grains settle exactly at vgas (but we show in Appendix C how the relative velocity can be included self-consistently. This has, however, the consequence that the grain size and thus opacity can no longer be found analytically. Instead, a root must be determined numerically in each layer.)

When the grains settle at vgas, the gas accretion rate and scale height alone define the advection timescale τadv. To determine the typical grain size, we can ask how big the grains can grow while they are advected through one scale height. This can again be estimated by setting the advection timescale equal to the growth timescale. Therefore, the basic approach to determine the grain size remains the same as in the settling regime. As in the previous regimes, we consider grain growth by Brownian motion coagulation and differential settling separately.

2.5.6. Advection regime: Brownian motion coagulation

In contrast to the settling regime, here we do not have to consider the Epstein and Stokes case separately. This is a result of the important difference that (in the case of complete advection we consider) τadv is independent of grain size. The equation that gives the grain size in this regime is thus τadv = τcoag,A, where(43)from Eqs. (19) and (5). This yields a grain size of (44)while the dust-to-gas ratio is simply fD/G,A = gr/XY which is, in the case of deep deposition (Eq. (4)), simply equal to fD/G,disk as it must be from mass conservation. The opacity is finally given as (45)which is for deep deposition (46)This shows that the opacity increases with both the gas accretion rate and the dust-to-gas ratio in the disk, in contrast to growth by differential settling in the settling regime.

2.5.7. Advection regime: differential settling

For differential settling we have to consider the Epstein and Stokes regimes separately as the expressions for the growth timescales differ. The latter regime is, however, not important, as advection only occurs in the upper atmospheric layers, and there the grains are in the Epstein regime.

When the grains are moving at vgas, the growth timescale (Eq. (30)) becomes (47)This growth timescale is independent of the grain size a. As the advection timescale is also independent of the grain size, this means that our normal approach of determining the typical grain size by equating τadv = τcoal,A,E and solving for a is not applicable here. This is a sign that in this regime, the grain size cannot be determined locally in one layer. Instead one needs to follow the growth as grains sink from the top of the atmosphere down to the layer in question.

To model this growth, we can use the definition of the growth timescale a/ȧ and the fact that the grains travel in this regime at vgas. This yields the equation that links the change of the grain size with the change of the radius R:(48)As shown in Ormel (2014), a related approach can be used to derive a general framework where the grain size is found by integrating numerically the expression for da/dR. Here we only search for an approximate solution to Eq. (48). For this, the following finding is useful. As a general result, the advection regime is found to occur first at or after the crossover point (Sect. 4.3.2) when the gas accretion rate starts to increase. It then first occurs in the very top layers at the outer boundary of the atmosphere. After crossover, the thickness of the advection layer then grows inwards as the gas accretion rate further increases.

We can therefore estimate the growth of a grain as it traverses the layer where the advection regime occurs by integrating Eq. (48) from the outer radius of the atmosphere at R0 down to the current radius R. At R0, the grains have an initial size a0 (this will typically be the same as the assumed monomer size amono). Inserting the equation for the growth timescale and vgas and integrating then yields the grain size(49)where we have assumed that the (contained) mass M is constant in the atmosphere. An analytical approximation of the integral is discussed in Appendix D. From the grain size, the opacity can then be calculated in an analogous way as in the other regimes. In Appendix D we compare Brownian motion and differential settling in the advection regime. It is found that the two mechanisms occur on comparable timescales leading therefore to similar grain sizes and opacities, in contrast to the normal case where growth by differential settling is much more efficient. In Sect. 4 we therefore only consider Brownian motion in the advection regime. It is clear, however, that this regime should be further investigated with a numerical approach that is more suitable than the local analytical description used here.

For completeness we finally give the grain size for advection and growth by differential settling in the Stokes regime. As mentioned, this regime does not usually have a practical meaning. It would occur when the gas accretion rate is very high, so that advection also occurs in the deeper atmospheric layers. This should only occur well after the crossover point, when the planet approaches the disk limited gas accretion rate. This phase is extremely short, and it is questionable whether the 1D radially symmetric approximation still holds. One finds a growth timescale (50)This timescale depends on the grain size, therefore it is possible to estimate the typical grain size by equating τcoal,A,S with τadv. This yields a grain size (51)

2.5.8. Relevance of the advection regime

The timescale arguments indicate that it is, in principle, plausible that at high gas accretion rates, the advection of grains becomes important, leading to high opacities in the outer layers. Indications of this are also seen in the numerical simulation of MP08 (see discussion in Sect. 3.1.6). In the current analytical model, only the limiting cases are considered, i.e., that growth and settling occur as if the atmosphere is completely static, or (the other extreme) where all grains move at vgas.

We find below that the model therefore predicts a relatively sudden increase of κ in the advection regime (Fig. 12). This is understood from the fact that the growth timescale in the differential settling and advection regime can be quite different (Fig. 4). In the first regime, the grains quickly reach sizes of ~0.01 cm, while in the advection regime, they can stay very small at a = amono (Fig. D.1). While the decrease of Q to values much smaller than 1 partially compensates for the increase of the absorbing surface for such small grains (Eq. (8)), we still find that significant jumps in κgr occur of up to one order of magnitude at the transition between the settling and the advection regime. Solving numerically for the grain size to include intermediate states for the relative velocity of grains and gas (Appendix C), as well as including growth by differential settling in the advection regime (Appendix D) partially changes this, but the basic prediction of a high opacity is found to remain6. On the other hand, the analytical model is built on a number of idealizations which in reality probably do not apply in the top layers. First, the outer layers partially participate in the general gas flow in the protoplanetary disk (Lissauer et al. 2009; Ormel 2013). Therefore, the layers could be turbulent and/or show patterns of atmospheric circulation, in contrast to the simplification made here of a purely radial laminar gas motion. Second, the grains that are accreted from the nebula into the planetary atmosphere are probably not monodisperse. Potentially, this would both make a faster growth due to differential settling possible, in particular by gas motions driving small grains against the large ones (Weidenschilling 1984). Such effects related to a bimodal size distribution are difficult to study in the context of the analytical model here, but should be addressed in future work with numerical models (Podolak 2003; Ormel 2014).

On the other hand, from a point of view that is mostly interested in the final outcome of the formation process (i.e., in knowing the final mass of a planet) rather than the atmospheric structure during all phases of the formation, the details of the advection regime are of secondary importance. This is because of the following (Sect. 4): the advection regime occurs (at least in the cases we studied) only at and (with increasing importance) after crossover, where it influences the time between crossover and the moment the disk-limited gas accretion rate is reached. Including or completely neglecting advection is found to lead to variations of this time of about 105 years (Fig. 10). This is short in comparison with typical disk lifetimes, which in turn means that the impact on the final properties of a planet will, in general, only be small.

2.6. Grain evaporation

A last effect we need to take into account is that grains evaporate at sufficiently high temperatures, so that the grain opacity disappears (Lenzuni et al. 1995). The mass loss rate of a grain due to evaporation can be calculated with the Knudsen-Langmuir formula (e.g., Bronsthen 1983) as (52)In this equation, μgr is the mean molecular weight of an evaporated grain molecule which was set to 33 (Pollack et al. 1986). We have assumed that the grains evaporate from the entire 4πa2 surface, that the partial pressure of the rock vapor in the atmosphere is small, and that the surface of the grains has the same temperature as the gas. The vapor pressure pvap is approximately given as pvap = pvap,0eCvap/T, where pvap,0 and Cvap are material properties. We take the values suggested by Podolak et al. (1988) for rocky material, so that pvap,0 = 1.5 × 1013 dyn/cm2 and Cvap = 56 655 K. These parameters mean that grains evaporate at a temperature of roughly 1300 K. It is worth noting that instead of the pure stoichiometric evaporation considered here, minerals can also decompose by reacting with the surrounding H2 gas. This can increase the mass loss rate (Boss et al. 2012).

The timescale of grain evaporation is (53)In principle, one could imagine that in certain parts of the atmosphere, grain growth (by Brownian coagulation or differential settling) is balanced by evaporation, so that the two processes must be combined to obtain the actual timescale on which the grain size changes. However, due to pvap there is a very strong exponential dependency of τevap on the temperature. In practice it means that either the evaporation timescale is many orders of magnitudes longer than the growth timescale at low temperatures, or many orders of magnitudes shorter as soon as the temperature is sufficiently high. This means that the radial domain in the envelope where evaporation changes from completely negligible to completely the dominating process is very small. This in turn means that there is a well-defined evaporation temperature and evaporation radius Revap.

Therefore, it is possible to simplify the treatment of evaporation. If τevap is longer than the relevant growth timescale τgrowth (given by one of the regimes of Sect. 2.5, typically differential settling regime with Stokes drag), it is possible to completely neglect evaporation. We need to take it into account as soon as τevap<τgrowth, in such a way that the grain opacity decreases rapidly towards the interior. Operationally, we therefore multiply in the latter case κgr calculated without evaporation by a factor of τgrowth/τevap. This factor very quickly approaches zero with decreasing radius within the evaporation radius. Clearly, this particular setting is used for convenience rather than physical reasons. However, the particular setting for the reduction of κgr because of evaporation remains in any case without consequences because of the following interesting result.

Because of grain growth, it is found (e.g., Fig. 11) that the grain opacity decreases with decreasing height in the atmosphere by about three orders of magnitude (as already found by MP08 and MBPL10). For this reason the opacity of the grain-free gas (molecular and atomic opacities) becomes larger than the grain opacity at a certain radius. This radius (Rκ − DG) is found to be larger than the radius where the grains evaporate Revap (Fig. 1). This means that grain evaporation remains without consequences for the actual total opacity. In other words, the grains are so large in the deeper layers that before they evaporate, their contribution to the total opacity becomes smaller than that of the grain-free gas. Our simplified treatment of evaporation is therefore justified.

It is clear that besides the destruction by evaporation, grains can also break up in mutual collisions. We follow in this paper earlier works on grain growth in protoplanetary atmospheres (Podolak 2003, MP08, and MBPL10) and currently neglect this effect. Instead, we simply assume perfect sticking. For the cases studied in detail here (see Sects. 3.1.4 and 4.3.4), it is found that the relative velocities of the grains are likely too small for fragmentation. Bouncing could, in contrast, be important. Future work should explore the importance of imperfect sticking.

thumbnail Fig. 2

Extinction efficiency Q as a function of xmax = 2πa/λmax. The coefficient is shown for grain sizes of a = 1, 10, and 100 μm, computed with Mie theory (solid lines) and with the fit of Eq. (54) (dashed lines). In all three cases, λmax corresponds to the wavelength of the maximum in the blackbody spectrum.

Open with DEXTER

2.7. Extinction coefficient Q

Once the grain size and dust-to-gas ratio is calculated, in order to get the opacity, we need to calculate the extinction coefficient Q that gives the ratio of the extinction cross section of a grain for radiation of wavelength λ to its geometric cross section πa2.

The extinction coefficient is the sum of the absorption and scattering efficiency. For spherical grains, it can be calculated with Mie theory. It is then a function of the real nr and imaginary ni refractive indices of the grain material (e.g., Dorschner et al. 1995) and of the size parameter x = 2πa/λ. As discussed in Podolak (2003) and MP08, instead of full Mie theory, it is possible to use approximative expressions for Q as a function of x, nr, and ni, derived by Dr. J. Cuzzi7.

In the context of protoplanetary atmospheres, it is first of interest to see which values of x occur. The typical wavelength of the radiation in the atmosphere at a temperature T can be estimated with Wien’s displacement law as CWien/T. The temperature in the envelope can vary from ~100 K in the outer parts of the atmosphere to about 1500 K, the temperature where the grains evaporate. The grain size will vary from about 1 μm, the monomer size, to about 0.1 cm in the deep layers (MP08, Fig. 8). This means that x roughly runs from 0.2 to 3000.

MP08 investigated the impact of different grain materials (tholine, olivine, and iron) that have different refractive indices. They found that the resulting opacity in the atmosphere is not very sensitive to the specific material type, i.e., to nr and ni. We therefore restrict ourselves in this study to tholins, the nominal case in MP08. We note that this is, strictly speaking, inconsistent with the assumption of silicate grains for the evaporation. In Fig. 2, the solid lines show Q(xmax) for grains with a = 1, 10, and 100 μm. The size parameter was calculated as xmax = 2πa/λmax(T), and T ran in each case from 10 K to 2000 K.

In Fig. 2, the extinction efficiency was calculated using a table of wavelength dependent refractive indices for tholins given by Khare et al. (1984) and the Mie code of Bohren & Huffman (1983). We then derived a simple empirical fit by comparison with the actual Q in Fig. 2, finding (54)This fit is shown with dashed lines in Fig. 2. For small x the fit approaches the absorption efficiency for small particles that is linear in x as found for isotropic Rayleigh spheres. It dominates over scattering (that would scale as x4) because of the complex refractory index (e.g., Hansen & Travis 1974). Expressed with nr and ni, the absorption efficiency can be approximated as (MP08)(55)For nr = 1.5 and ni = 0.1, the typical values mentioned by MBPL10, one finds Qabs = 0.20x, while using nr = 1.6 and ni = 0.16 (which we find to provide a somewhat better fit to the tholine data of Khare et al. 1984) yields Qabs = 0.28x. This is close to the value given by the fitting function that was derived by comparing visually a fitting function and the Mie results. On the other hand, for very large x, the fit approaches the geometrical limit of Q = 2.

Equation (8) gives the wavelength-dependent opacity of a grain of size a for radiation of wavelength λ. For the calculation of the atmospheric structure, one needs in contrast the Rosseland mean opacity. In the light of an absence of a grain size distribution in the simple analytical model, and, more importantly, the fact that under most circumstances aλ, so that Q ≈ 2 (i.e., without a significant wavelength dependency for the contributing part of the Planck function), here we use the following simplification: instead of calculating the actual Rosseland mean, we use the opacity as found with Q(x) evaluated at xmax, i.e., at the wavelength λmax(T) of peak flux in the Planck distribution, where T is the temperature of the gas in the atmospheric layer.

The comparison with the results of MP08 and MBPL10 (where the actual Rosseland mean opacity is calculated) indicates that this simplification does not have a significant impact on the predicted grain opacity, at least for the cases studied, because grain growth is so efficient that we are in the geometrical limit. This is confirmed by Fig. 10 in Cuzzi et al. (2014): from the figure it becomes clear that the wavelength-independent geometrical limit for Q (and therefore our approximation) is applicable for grains of 10μm if the temperature is larger than 300 K, and for grains of 100 μm (or larger) if T ≳ 50 K. These conditions are usually met in protoplanetary atmospheres except for the advection regime where the grains remain small in the outer layers of the atmosphere. Fortunately, this regime occurs only when the planet is already at or close to the critical core mass, so that the impact on the total formation history of a planet remains small (Sect. 4). Nevertheless, this is a point that will be further addressed in future work.

2.8. Combination of the different regimes, gas opacity, and final expression for κ

In the previous sections, we derived the timescales of several processes that govern the evolution of grains and the associated opacity. To get the final expression of κ, we combine them, following the general approach that the process with the shortest timescale is taken as the dominant one (Cooper et al. 2003).

Operationally, we first calculate all four growth (and by definition equal settling) timescales for Brownian coagulation and for differential settling in both drag regimes. The appropriate drag regime is then chosen by checking which regime leads to a self-consistent Knudsen number (see Sect. 2.2.1). Next, the timescale of growth by Brownian motion and differential settling are compared, and the one with the shorter timescale is retained. This timescale is then compared with the advection timescale, and again the shorter one is retained. Finally, this timescale is compared with the evaporation timescale and if it is shorter, the opacity is reduced as described in Sect. 2.6.

Up to this point, we have only considered the opacity due to grains κgr. To get the total opacity, we also need to take into account the molecular and atomic opacity of the grain-free gas, κgas. This gives the opacity in the deep convective interior of the envelope, but also in the deeper parts of the atmosphere because the grain opacity becomes so small there due to grain growth (see Fig. 7). The gas opacity is obtained from the combination of the analytical expressions of Bell & Lin (1994) for T> 4000 K and the tables of Freedman et al. (2008) at lower temperatures, assuming a solar-composition gas. As already noted by MBPL10, the details of the opacities at high temperatures are not important, because the energy transport is due to convection.

For simplicity, we set the resulting total opacity equal to κ = MAX(κgr,κgas). It is found that the specific way the gas and grain opacity are combined is not very important, because at the radius where the two become equal, κgr rapidly decreases towards the interior, whereas κgas rapidly increases (Fig. 7). This means that the radial domain where the two are of the same magnitude is small.

2.9. Limitations of the model

This first simple analytical model has a number of limitations. An obvious one is the absence of a grain size distribution in each atmospheric layer. As explained in MP08, it is not assured that the opacity in a layer is always dominated by the grains with the most frequent size. For example, the analytical model predicts that in deep layers, there are only large grains because it assumes that small grains are brought into the atmosphere only at its outer boundary. If in reality, planetesimal impacts deposit many small grains in the lower layers, these could dominate the opacity, even if they do not dominate the size distribution. On the other hand, in our comparisons with the results of MP08 and MBPL10, the consequences of a grain size distribution instead of one single size becomes apparent only in one rather minor feature, which turns out not to have consequences on the overall opacity (Sect. 3.1.1). Nevertheless, more work is needed to further investigate and improve this point. It could be possible to extend the model by including size distributions centered on the typical size, with a shape and width typical for either Brownian coagulation or growth by differential settling.

A second, related important limitation is that the analytical model assumes that the total flux of grains remains radially constant. This makes it impossible to capture the consequences for κgr if planetesimal impacts locally increase the dust input. Our results indicate that this should typically not have very strong consequences for the overall opacity (and therefore optical depth) for two reasons. First, the opacity in the most common regime (Eqs. (35) and (41)) is independent of gr. This picture is, roughly speaking, confirmed by the results of MP08 (their Fig. 11): a very strong increase in the planetesimal input (by a factor of 100) in a certain layer leads locally to a higher opacity, but the total optical depth only changes by a factor of a few (except in the special case that the increased deposition occurs just at the bottom of the atmosphere). Second, at least for km sized planetesimals, most mass deposition should occur deep enough not to influence the grain opacity much (see Fig. 1).

A third limitation is that we assume that the opacity of the grains for radiation at λmax can be used as the Rosseland mean opacity. In the most important growth regime (differential settling) and in most parts of the envelope, aλmax, so that the extinction efficiency is largely independent of the wavelength for the relevant part of the Planck function. In this regime, the difference between κ(λmax) and the Rosseland mean is therefore small. This is confirmed in tests where we calculated the Rosseland mean. On the other hand, in the advection regime, there are small grains, and the wavelength dependency becomes more important. Tests show that the Rosseland mean and κ(λmax) can differ in this regime by up to one order of magnitude. On the other hand, the advection regime only occurs at or after crossover, so that it is less important (Sect. 4).

A fourth important limitation is the assumption of perfect sticking despite the fact that bouncing may be important (Sects. 3.1.4 and 4.3.4). This is a limitation that is shared with the numerical models of MP08 and MBPL10. Additional limitations that are shared with these numerical models are the assumptions that the grains are spherical, that no turbulence occurs in the envelope, and that only silicate grains are considered. The effects of convection are also ignored in the analytical model because MBPL10 showed that convection hardly influences their results.

3. Comparison with numerical models

We now compare the analytical expression for the opacity derived in the previous sections with the results of MP08, and then of MBPL10 (Sect. 4). Building on the initial work of Podolak (2003), MP08 and MBPL10 use a complex (and computationally very heavy) numerical model to simulate the temporal evolution of grains under the influence of vertical settling and growth by Brownian coagulation and differential settling. In contrast to the analytical model, at each layer and moment in time they have a full distribution of grain sizes, expressed in the number of grains of a given size per unit volume. The temporal evolution of this quantity is obtained by solving numerically the Smoluchowski equation. The Smoluchowski equation contains the terms representing coagulation and settling, plus a source term. The source term has two contributions: first, grains that are deposited by impacting planetesimals, and second, grains that are accreted together with the gas. The radial distribution of the grain input due to planetesimals is given directly by the radial mass deposition profiles of the impacting planetesimals obtained with the code of Podolak et al. (1988). This is in contrast to the analytical model where gr has to be radially constant. The grains accreted with the gas are injected into the envelope at the top layer (which is the same as in the analytical model). In their nominal model, MP08 assume that the monomer size both of grains accreting with the gas and resulting from planetesimal deposition is approximately 1 μm.

The effects of envelope contraction which advects small grains (Sect. 2.5.5) and of convection are also included in some simulations. Once the grain size distribution is calculated, the wavelength-dependent opacity is computed using an approximation to Mie theory (Cuzzi et al. 2014). Finally, the Rosseland mean is calculated using the Planck function at the local temperature of the gas. It is clear that the numerical model is itself simplified in a number of points, of which the most important were mentioned in the last section. This means that the results of the numerical model must be regarded as an approximation, too.

Nevertheless, the results of the numerical model are crucial in this study, as they are the best currently available information on how grains evolve in protoplanetary atmospheres. Therefore, they serve as the benchmark cases. Ideally, the analytical model would perfectly reproduce the results of the numerical model. This is of course impossible because of the numerous simplifications. Still, we find below that the numerical model and the analytical approximation agree surprisingly well and share many common important features, also in a quantitative sense. A litmus test for the analytical model is its prediction for the opacity if the dust-to-gas ratio in the accreted gas fD/G,disk is varied. We will see that the analytical model passes this test very well, and that this is a very telling result (Sects. 3.1.6 and 3.1.7).

The analytical model yields in particular a radial opacity structure that is much closer to the numerical result than the one obtained by reducing the ISM grain opacity by some constant factor (Fig. 7), as was often done in many previous works of different groups as mentioned in the introduction. This is one of the goals of this work. The most important input parameters for the comparison cases are given in Table 1.

Table 1

Settings for the comparison with MP08.

3.1. Comparison with MP08

The first numerical simulation to which we compare our analytical model is the calculation of MP08 for a giant planet forming at 5.2 AU. This simulation is very important for the validation of the analytical model because it is the only published numerical simulation with fairly complete information on the grain dynamics. This means that not only is the radial profile of the opacity known (as in MBPL10), but also the shape of the grain size distribution at different heights in the atmosphere, the mass density and settling velocity of grains as a function of height, the impact of different monomer sizes, and, very importantly, the effect of different dust-to-gas ratios in the newly accreted gas. The last group of information is important as it helps us to understand which growth mechanism (Brownian motion or differential settling) is dominant.

thumbnail Fig. 3

MP08 comparison case. Left panel: atmospheric gas density (solid line) and temperature (dashed) as a function of height (distance from the planet’s center) in the protoplanet’s atmosphere (outer radiative zone). Right panel: opacity as a function of height in the numerical model of MP08 (dashed line: grain opacity only) and in the analytical model of this work (thick solid line: combined grain and gas opacity). The thin blue solid (dashed) line is the gas (grain) opacity alone.

Open with DEXTER

The atmospheric temperature and density structures used by MP08 were taken from Hubickyj et al. (2005) and correspond to a moment in time when the planet is in the middle of phase II (see Pollack et al. 1996 for the different phases occurring during classical in-situ giant planet formation simulations). At this moment in time, the mass of the core is 12.61 M, and the envelope mass is 2.73 M. The outer radiative zone/atmosphere of the planet contains very little mass, therefore we assume that the mass interior of the atmosphere driving the settling is equal to the sum of the two. In MP08 the gas accretion rate at this moment (which we need in order to calculate the grain accretion rate, Eq. (4)) is not explicitly stated, but the simulations of Hubickyj et al. (2005) indicate that it should be approximately 8 × 10-6M/yr. We assume that the planetesimals deposit their mass deep into the convective zone so that the advection of grains with the newly accreted gas is the only source of grains. Our results are neither sensitive to the exact value of XY nor to the assumption of deep/shallow deposition (Sect. 3.1.6). In MP08, grains are brought into the atmosphere both by gas accretion in the top layer as in the analytical model, but also at other heights by planetesimal ablation (which is not included in the analytical model). It is therefore not clear a priori if a useful comparison can still be made. However, as discussed in MP08, the planetesimal source is important only in the deepest layers because these grains are released mostly in the deep atmosphere (below ~1011 cm, see their Fig. 9). Therefore, it is possible to compare the analytical and numerical result.

We note that the calculation of κgr in MP08 is not self-consistent in the sense that the temperature and density structures were obtained by Hubickyj et al. (2005) with a pre-specified (and different) opacity, and not the one obtained in the grain calculations. In contrast, the calculations of MBPL10 discussed in Sect. 4 are self-consistent. For our purposes of comparison of analytical and numerical models, however, this has no negative impact.

3.1.1. Atmospheric structure and opacity

The left panel of Fig. 3 shows the atmospheric gas density (solid blue line, left y-axis) and the temperature (brown dashed line, right y-axis) as a function of height (distance from the planetary center). These lines are simply taken from Fig. 1 in MP08. As in MP08, only the outer radiative zone of the gaseous envelope is considered. The density increases from about 10-10 g/cm3 at the outer radius (which is comparable to typical nebular densities) to about 10-6 g/cm3 at the inner boundary. The temperature structure consists of an approximately isothermal outer part at 150 K (again the nebular temperature), followed by an increase to about 670 K at the inner boundary. This temperature is so low that silicate grains do not evaporate in the atmosphere.

The right panel of Fig. 3 directly shows the main result of the calculations which is the opacity as a function of height. The results of MP08 and of the analytical model are both shown. One sees a radial opacity structure that is characteristic for all calculations, at least during phase II: At the outer boundary, the opacity is high, on the order of 1 cm2/g, and therefore comparable to ISM opacities. Towards the interior, it decreases strongly, going down to about 10-3 cm2/g at the inner boundary. It is the consequence of the increase of the typical grain size a and of the decrease of fD/G with decreasing height (see Fig. 8 below). Clearly, this is a consequence of grain growth.

This general structure is found both in the numerical and analytical model. We note that for the former, only the opacity due to grains κgr is shown, while for the latter, the thick line shows the total opacity, including the contributions of gas and grains. The thin dashed line shows the grain opacity alone in the analytical model. We see that in the analytical model, κgr is a strictly monotonically increasing function of the height. In the numerical model, κgr has, in contrast, a local minimum at about 1011 cm. Inside this point, it increases again. The reason for this is explained in MP08. At the minimum, the drag law for the large grains changes from Epstein into the Stokes regime, leading to an increase in the settling velocity (we find the same behavior in the analytical model, see Fig. 6). The small grains (which are brought deep into the atmosphere by local planetesimal mass deposition, see MP08), in contrast, still remain in the Epstein regime. This causes a shift in the grain size distribution, resulting in an increase of κgr. Because we do not consider planetesimal ablation as a local source of small grains deep in the envelope (as mentioned, gr needs to be radially constant in the analytical model), and because we only have one grain size (which corresponds to the big grains in MP08), we cannot have such a behavior in the analytical model. This is the only point where a difference between a grain size distribution, and only one representative size becomes apparent, at least in the comparisons made here.

thumbnail Fig. 4

MP08 comparison case. Left panel: timescale of growth by differential settling/coalescence (blue solid line) and Brownian coagulation (green dash-dotted). These timescales are by construction equal to the associated settling timescale. The timescale of gas advection is shown by the brown dashed line. In this atmosphere, the evaporation timescale of the grains is always many orders of magnitude longer because of the low temperatures. The thin gray dashed and dash-dotted lines show the result if the velocity of the gas is directly included to calculate the growth timescales of differential settling and Brownian coagulation, respectively (Appendix C). Right panel: ratio of the stopping to the settling timescale RSS (solid) and stopping timescale (dashed line) of the grains.

Open with DEXTER

Interestingly, the difference between the numerical and analytical model decreases for the following reason. If one considers the total opacity as predicted by the analytical model including the opacity due to the grain-free gas which is the physically relevant quantity, then we see that the total opacity in the analytical model also has a local minimum at a radius 6 × 1010 cm. Within this distance, the total opacity also increases rapidly, but now due to the molecular opacities which were not considered in MP08. Clearly, this is a priori a chance result, but we find in the other comparisons (Sect. 4.3.1), that this is actually a general behavior: in the deeper parts of the atmosphere, the grain opacity decreases so much with decreasing altitude that the gas opacities become dominant at some height. Below this height, the opacity increases again, but now due to the molecules. This means that there is a local minimum in the opacity. It is found that the analytical model predicts both the location of the local minimum and its value (typically a few 10-3 cm2/g) in a way that agrees well with the results from the numerical calculations (see Figs. 11 and 13). This makes it possible to use the analytical model also for quantitative calculations.

We also see that the radial structure of the opacity in the two models can differ by up to one order of magnitude at some heights. Interestingly, the shape of κ agrees better in the simulations where the atmospheric structure is calculated self-consistently (see Figs. 11 and 13). In the simulation here, at an intermediate height of about 2.5 × 1011 cm, the opacity is about a factor of 10 higher in the analytical model. The total optical depth of the atmosphere on the other hand agrees to a factor of 2–3 in the two models (Fig. 7). The reason for the difference could be an overestimation of the efficiency of grain growth in the outer and intermediate layers in the analytical model, leading to grains of 10 μm in size which cause a relatively high opacity (see MP08 and the discussion below). In the nominal numerical model of MP08 that assumes a monomer size of 1 μm, the grains in the outer layers are smaller, causing a lower κgr because of their low Q. On the other hand, if a monomer size of 10 μm is used in the numerical model, then the results for κgr are in better agreement with the analytical result (see Fig. 9).

3.1.2. Timescales

We now discuss the microphysical grain dynamics that lead to the opacity. The analytical model is based on timescales; therefore, we first show in the left panels of Fig. 4 the timescales of growth by differential settling, Brownian motion coagulation, and gas advection. By construction, in all three regimes, the growth and settling timescales are the same.

Movshovitz & Podolak (2008) found that the grain evolution approaches a steady state in only 100 yrs, and that after 1000 yrs, the steady state is fully established. Similarly short timescales are also found in the analytical model. The timescale for Brownian coagulation varies between about 500 years and 2000 years. The timescale for gas advection is also about 1000 years in the outer layers and increases to about 6000 years at in the inner layers which is, however, irrelevant. The process with the shortest timescale is growth by differential settling. It occurs on a timescale of about 300 years in the outer layers, and of less than 20 years at the radiative-convective boundary. This is indeed a very short timescale.

thumbnail Fig. 5

MP08 comparison case. Left panel: grain radius (solid) and dust-to-gas ratio (dashed line). Right panel: Knudsen (solid) and Reynolds number (dashed line) of the grains as a function of height. The horizontal line indicates the critical Knudsen number of 4/9 where the drag regime changes.

Open with DEXTER

As did Cooper et al. (2003), we assume that the process with the shortest timescale is the dominant one. The figure shows that, in the entire atmosphere, differential settling is the dominant process, meaning that the opacity is given by either Eq. (35) or (41), depending on the Knudsen number. The dominance of differential settling is a general result, which applies to most atmospheres (see Sect. 4.3.2). Coagulation is about one order of magnitude slower in the outer layers, and a factor of 30 slower in the inner layers. Furthermore, gas advection is clearly slower, and therefore not important. This can be seen from the negligible change of the timescale for differential settling if the effect of the gas velocity is directly included (gray dashed line), as described in Appendix C. This will be partially different in the simulation discussed in Sect. 4.3.1. There, the core is more massive at the crossover point, so that the gas accretion rate is higher, and therefore the gas advection timescale shorter.

Regarding the importance of Brownian coagulation versus growth by differential settling, one could argue that for differential settling, grains of different sizes must already exist. In the deeper atmosphere, this should always be the case because grains of different sizes rain into it from above. At the outer radius, in the limit that the initial size of the grains is identical (which is probably not realistic as grains already grow in the protoplanetary disk), this is not the case. Therefore, at the outer boundary, it is likely that first the slower process, i.e., Brownian coagulation, sets as the bottleneck the timescale up to a moment when differential settling kicks in. This would mean that Brownian motion needs first to act as the trigger of the coagulation process as is the case in disks (Weidenschilling 1984; Brauer et al. 2008). This is an effect that we do not capture in the analytical model. Therefore, we likely overestimate the efficiency of grain growth in the top layers, and therefore underestimate the opacity. We will indeed see indications of such a behavior when studying the effect of different dust-to-gas ratios in the disk (Sect. 3.1.6).

We have assumed that grains always settle at the terminal velocity. The right panel of Fig. 4 shows that this is justified. It shows that the stopping time of the grains is very short, 10–1000 s. As the settling timescales are on the order of 10–100 years, this means that the ratio of stopping to settling timescale RSS is tiny. As shown by the figure, it has values on the order of 10-8 to 10-7, meaning that the grains very quickly reach the equilibrium between drag and gravitational force.

3.1.3. Grain size, dust-to-gas ratio, and aerodynamic regime

Now that differential settling has been identified as the process determining the typical grain size, local dust-to-gas ratio, and opacity, it is interesting to study these quantities. The left panel of Fig. 5 shows the typical size of the grains a as a function of height. At the outer boundary, the grains are relatively small (3 × 10-3 cm). This size is a factor of 30 larger than the assumed monomer size. In the analytical model, the size found by equating the timescales is almost always larger than the (typical) monomer size (~1 μm). Exceptions are the rarely occurring advection regime or when large monomer sizes are assumed (100 μm). The results of the analytical model are therefore usually independent of the assumed monomer size. This only holds because we do not also inject grains at the monomer size into the deeper layers.

Farther down into the atmosphere the grains then grow to a size of about 3 × 10-2 cm at a height of slightly less than 1011 cm. Within this height, the size remains approximately constant. This change at 1011 cm comes from a transition in the drag regime, as we will see below in this section.

It is interesting to compare the typical size predicted by the analytical model with the grain size distribution found by MP08 where the grain size distribution is shown for four layers in the atmosphere (their Fig. 4). At the outer boundary, their size distribution has only one global maximum which is at the assumed monomer size (1 μm). Inside the atmosphere, the size distribution is bimodal. One larger maximum is always found at the monomer size (1 μm). This peak is likely due to locally injected monomers originating from planetesimal ablation. The second maximum in contrast shifts to larger grain sizes as we move inwards, as in the analytical model. It is likely that this peak corresponds to grains which have already settled through a significant part of the atmosphere while growing. They therefore correspond to the grains studied here. In Table 2 we compare the grain size of these maxima with the typical grain size found in the analytical model.

Table 2

Grain size a in cm in MP08 and the analytical model as a function of height in the atmosphere.

The table shows that except for the top layer, the maximum of the grain size distribution in the numerical model and the typical size found here agree relatively well, namely to a factor of two or less. This indicates that despite its simplicity, the analytical model is able to describe the essence of the grain evolution quite well, also quantitatively.

thumbnail Fig. 6

MP08 comparison case. Left panel: settling and by construction identical growth timescale in years (solid line) and settling velocity (dashed line) as a function of height. The thin dash-dotted line is the settling velocity found by MP08. As they only give the relative velocity vset(R) /vset(Rout), we have normalized their curve to the value predicted by the analytical model at Rout. Right panel: dust mass density (fD/G × ρ) found by the analytical model (blue solid line) and by MP08 (thin dash-dotted line). The two brown lines are the extinction efficiency, calculated with the fit of Eq. (54) (dashed) and directly with Mie theory (dash-dotted line).

Open with DEXTER

The left panel of Fig. 5 also shows the local dust-to-gas ratio fD/G as a function of height. It scales as 1/(R2ρvset) (Eq. (5)). The plot shows that fD/G decreases by two orders of magnitude across the atmosphere, showing that the effect of an increasing gas density and settling velocity of the grains (owing to grain growth) dominates over the decrease of R. This decrease of fD/G and the increase of the grain size a both conspire to the strong decrease of κgrfD/G/a in the deeper parts of the atmosphere.

When calculating gr, it was assumed that fD/G,disk = 0.01 as in MP08. The plot shows that the model predicts a sudden jump down to about 0.003 in the outermost layer. Such a non-continuous behavior is in principle not physical. It is a consequence of the simplification that the growth by differential settling sets in instantaneously already in the outermost shell, which overestimates the growth as discussed in Sect. 3.1.2 and therefore artificially reduces fD/G.

The right panel of Fig. 5 shows the Knudsen and Reynolds number as a function of height. In the outer parts, the Knudsen number is much larger than 4/9, therefore the Epstein regime applies. As we move deeper into the atmosphere, the mean free path of a gas molecule decreases while the grain size increases, so that at a height of about 1011 cm, the drag changes into the Stokes regime. This affects the settling velocity as we will see in the next section, and therefore also affects the grain size.

The Reynolds number is very small in the outer atmospheric layers (~10-7) and increases to about ~10-2 at the boundary of the convective zone. This shows that the grains always settle in a laminar flow, meaning that the expressions for the drag laws we use are appropriate.

3.1.4. Settling velocity, dust density, and extinction efficiency

The left panel of Fig. 6 shows the settling velocity. We see that the settling velocity is about 25 cm/s in the top layers. It then decreases to about 5 cm/s at a height of 1011 cm. After this local minimum, it increases again, reaching about 18 cm/s at the boundary of the convective zone. This increase is related to the change of the drag regime, as mentioned by MP08. In the Epstein regime, the settling velocity scales as , while in the Stokes regime it scales as without a dependency on the gas density. This means for the Stokes regime that the settling velocity increases faster with increasing particle size a (or in other words, that it becomes more difficult to grow larger, as visible in Fig. 5 ), and that the settling velocity does not decrease with increasing density as we go deeper into the atmosphere.

The settling velocity found in the analytical model can be compared with Fig. 8 in MP08. A direct quantitative comparison is difficult because MP08 plot only the normalized velocity vset(R) /vset(Rout) relative to the value in the top layer at Rout for two specific grain sizes (0.5 and 0.005 cm). In Fig. 6, we also include their result for 0.5 cm grains, normalized to the settling velocity given by the analytical model at Rout. Nevertheless, we can see that the shape of the curves is similar in both models with the maximum value at the top, the local minimum at about 1011 cm for the large grains, and the re-increase below this height.

The settling velocity is also important from another point of view. In this simple model presented here, we implicitly assume that all particle collisions lead to perfect sticking. The very high number of studies regarding this subject in the context of grain as well as planetesimal growth in protoplanetary disks (e.g., Blum & Münch 1993; Güttler et al. 2010; Okuzumi et al. 2012; Windmark et al. 2012 to name just a few), shows that this is an optimistic assumption, which is justified only for certain collisional regimes. In other regimes, fragmentation or bouncing occurs. The outcome of a collision depends on several quantities (see, e.g., Blum & Wurm 2008), namely the particle sizes, the collisional velocity, the particle inner structure (fractal dimension, compact versus porous), and the material type (ice versus silicates).

Laboratory experiments indicate (see Güttler et al. 2010) that the critical velocity for collisional fragmentation of similar sized particles made of SiO2 monomers is approximately 100 cm/s. This is more than the settling velocity found in Fig. 6 (the collision velocities should be a fraction of settling velocity, depending on the relative size), so that fragmentation should not occur.

thumbnail Fig. 7

MP08 comparison case. Left panel: opacity as a function of height. The brown dashed is the result of MP08. Other lines show the total opacity found in this work (blue solid line), the gas opacity only (green dash-dotted line), the full ISM opacity (thin black solid line) and the ISM opacity × 0.003 (thin black dashed line). Right panel: corresponding optical depths.

Open with DEXTER

However, when comparing the grain mass-settling velocity relation predicted by the analytical model for this atmosphere with Fig. 11 in Güttler et al. (2010; see also Sect. 4.3.4), one sees that the collisions could fall into the bouncing regime, except if the collision velocities are only on the order of 1% of the settling velocity, or if mostly particles of significantly different sizes collide. In this case, the collisions would again lead to growth. We note that there is currently a disagreement between numerical simulations (where little bouncing is observed) and laboratory experiments (Windmark et al. 2012). The uncertainty associated with the impact of bouncing should be kept in mind when interpreting the result of the present work. To better understand this mechanism, future models of grain growth in protoplanetary atmospheres should build on the apparatus that has been developed in the past to study grain growth in protoplanetary disks (Ormel & Mordasini, in prep.).

The left panel of Fig. 6 shows again the timescale of growth by differential settling (coalescence) which is the relevant regime. It is clear that a sticking efficiency of less than unity would also affect the timescales. Movshovitz & Podolak (2008) investigated the effect of uniformly reduced sticking efficiencies. At a lower efficiency of 0.4, the opacity increased mildly. Lower efficiencies also mean that the time until a steady state is reached increases. At an efficiency of 0.1, the steady state is established after ~3500 years instead of less than 1000 years.

The right panel of Fig. 6 shows the mass density of dust (ρ × fD/G) as a function of height. Movshovitz & Podolak (2008) show the same quantity in their Fig. 5, and their curve is also included here. We see that the two lines have a similar shape, but that there is a certain offset. Even relatively fine structures like the small bump close to 1011 cm are seen in both models. Interestingly, the bump only appears clearly in the product of ρ and fD/G, but not in the two quantities alone.

There are also differences: in MP08, the dust density increases outside of about 4.5 × 1011 cm to reach at the outer boundary the value that is simply given by the nebular density times the dust-to-gas ratio in the disk, corresponding to about 10-12 g/cm2. In the analytical model, the density only decreases towards the exterior. This is probably a consequence of the mentioned shortcoming of the analytical model that predicts that differential settling dominates everywhere, while in the numerical model, Brownian motion dominates in the outer layers (MP08). The comparison allows us to deduce that the thickness of the layer where Brownian motion dominates is about 0.5 × 1011 cm. A similar length scale is also seen in the numerical model when the dust-to-gas ratio in the disk is varied (Sect. 3.1.6).

A second difference is that the analytical model predicts in the middle and lower parts a mass density of dust that is about a factor of 3 larger than in MP08. This could be the reason why the opacity predicted by the analytical model is also (for most of the atmosphere) higher than in the numerical model (Fig. 3).

In the right panel of Fig. 6 the extinction efficiency Q also is given. For the dashed line Q was calculated using the fit from Sect. 2.7, while for the dash-dotted line it was calculated directly with Mie theory. We see that in both cases, Q is approximately equal to 2 in the entire atmosphere. This is the consequence of aλmax. This is again a general result. We find that in most atmospheres, grain growth is so efficient that Q is to good approximation equal to 2. This means that the expressions derived for κgr in Eqs. (35) and (41) are useful ones as, typically, they can be directly used for the Rosseland mean opacity (except for the outermost layers). Finally, we have tested if the results for Q change if we calculate the Rosseland mean opacity instead of evaluating Q only at the wavelength of maximum flux as described in Sect. 2.7. We found that the results are again very similar because the grains are so large and the temperature sufficiently high that we are in the wavelength-independent regime (Cuzzi et al. 2014).

3.1.5. Comparison with scaled ISM opacity and optical depth of the atmosphere

In the past, numerous studies have used the ISM grain opacity multiplied by a constant reduction factor to mimic the effects of grain evolution (see Paper I). It is therefore interesting to compare the radial structure of the ISM opacity and the opacity found with the grain evolution calculations. Figure 7, left panel, shows the opacity as found by MP08 and the analytical model together with the full ISM opacity (from Bell & Lin 1994) and the ISM opacity multiplied by 0.003. This reduction factor was determined in Paper I as the fitting factor that leads to a duration of phase II that is similar as found in Movshovitz et al. (2010). The figure also contains the grain-free gas opacity for a solar composition gas from Freedman et al. (2008). The right panel shows the corresponding optical depths of the atmosphere, τatmo.

thumbnail Fig. 8

MP08 comparison case. Left panel: opacity as a function of height for eight different dust-to-gas ratios in the accreted gas fD/G,disk indicated in the plot. Right panel: corresponding size of the grains (solid lines) and local dust-to-gas ratio (dashed lines) for the four dust-to-gas ratios in the accreted gas indicated in the plot.

Open with DEXTER

The plot shows that the full ISM opacity is more or less radially constant at a high value (~1.3–4 g/cm2), in clear contrast to the dynamically calculated opacity. The ISM opacity is therefore up to a factor of 2000 higher (in the deeper layers) than the opacity predicted both by the numerical and analytical model. This is a very large difference. The corresponding total optical depth is also about a factor of ~400–800 larger than predicted by the grain calculations. The optical depth for full ISM opacity is approximately 4 × 104, while the numerical and analytical models predict τatmo of 90 and 50, respectively. They therefore agree to a factor of ~2.

The ISM opacity multiplied by 0.003 is lower than the result from the microphysical calculations in the outer layers, but higher in the inner layers, meaning that it has a clearly different (roughly speaking constant) radial shape. The differences are up to one and a half orders of magnitude. The total optical depth (110) on the other hand agrees well with the numerical result. This is probably not surprising, since the reduction factor was derived by fitting similar simulations (Paper I).

The grain-free opacity is also different. It only decreases towards the outer layers, which is also in contrast to the grain calculations. This means that it predicts a vanishing optical depth for the atmosphere in the range 1–5 × 1011 cm, whereas the calculations including the grain dynamics predict an optical depth of at least 10 for this part, which is not negligible. In the deeper layers, the grain free opacity dominates over the grains as described earlier. This means that the total optical depth is about 25 even without grains. While this is significantly lower than the result from the grain calculations, it still yields an interesting insight: using grain-free opacities (see Hori & Ikoma 2010) provides an optical depth that is closer to the result from the grain calculations than the full ISM value. On the other hand, the substantially different radial structure and the too low total optical depth mean that a completely grain-free opacity also differs from the result obtained when the grain dynamics are considered.

In retrospect it means that classical studies like Pollack et al. (1996) should have considered the grain-free case as a meaningful alternative to the case that grains exist in the atmosphere. Potentially, this would have meant that the classical timescale problem of the core accretion theory would not have been regarded as an important drawback of this theory for several decennia.

3.1.6. Impact of the dust-to-gas ratio in the disk

The impact of different dust-to-gas ratios in the disk fD/G,disk and therefore different gr = fD/GXY on the opacity is the litmus test showing which grain growth mechanism (Brownian coagulation versus differential settling) dominates. For Brownian coagulation, κgr increases for a fixed gas accretion rate XY as and in the Epstein and Stokes regime, respectively (Eqs. (23), (28)). While these are not very strong dependencies, they are nevertheless in clear contrast to the growth by differential settling where κgr is independent of fD/G,disk both for the Epstein and Stokes regimes (Eqs. (35), (41)).

Movshovitz & Podolak (2008) investigated the impact of varying fD/G,disk between 0.001 and 0.02, i.e., by a factor of 20. For Brownian coagulation, this would lead to a variation of κgr by a factor of approximately five according to the analytical estimates. But MP08 actually found that in the dominant part of the atmosphere (between a height of about 1–4.5 × 1011 cm), the opacity was nearly independent of fD/G,disk. The analytical model now gives an explanation for this behavior. It is also a very clear sign that differential settling is indeed the dominant growth mode and not Brownian coagulation, as indicated by the analytical timescale criteria. This is a result with important wider implications which we discuss in Sect. 3.1.7. Within 1× 1011 cm, MP08 found that fD/G,disk is slightly anti-correlated with κgr (variation by a factor of 3). This is likely a consequence of a change in the grain size distribution, an effect that we cannot capture in the analytical model (see Sect. 3.1.1). In contrast, outside 4.5× 1011 cm, the opacity increases with fD/G,disk (variation by a factor of ~20). This is again a clear sign that in the outermost layers, Brownian coagulation is important (Movshovitz & Podolak 2008).

The left panel in Fig. 8 shows the opacity for fD/G,disk between 0.0001 and 0.5, a very wide (and not necessarily realistic) range. One sees that for all fD/G,disk except 10-4, the model predicts a very similar opacity. The variation that still exists is a consequence of the decrease of Q(xmax) that approaches 2 (Fig. 2) as the grains become larger (see next). Only for fD/G,disk = 10-4 does the opacity decrease in the outer parts of the atmosphere. The reason is that at such low grain concentrations, the timescale of growth by differential settling becomes longer than the gas advection timescale. The growth regime therefore changes in the outer layers into the advection regime, resulting in low opacities because of the very low fD/G. The analytical model thus agrees with the numerical result that the opacity is nearly independent of fD/G,disk in the dominant part of the atmosphere. It differs by predicting that this independence holds for the entire atmosphere, including the top layers. This is attributable to the overestimation of the efficiency of differential settling also in the outermost layers that was mentioned earlier (Sect. 3.1.2).

The right panel illustrates the reason for the independence of κgr on fD/G,disk. The figure shows the corresponding grain size a and the local dust-to-gas-ratio fD/G as a function of height. For an approximately constant Q (because of the large grains relative to the wavelength), κgr scales simply fD/G/a (Eq. (8)). One sees that high fD/G,disk indeed lead to a higher fD/G which would increase the opacity, but at the same time these higher grain concentrations also lead to larger grains, reducing the opacity. As fD/G and a depend in the settling regime in the same way on fD/G,disk, these two effects cancel. We have checked that even for the largest grain, the Reynolds number stays below unity, meaning that the applied drag laws remain valid. We note that for the very large grains expected for high fD/G, bouncing and fragmentation could become important (Sect. 4.3.4). This could limit the efficiency of grain growth in very metal-rich envelopes.

Up to now, we have discussed the effect of varying the grain input in the top layer that is caused by a different fD/G,disk in the accreted gas. Movshovitz & Podolak (2008) also studied the effect of varying the grain input due to the ablation of planetesimals. For an increase of the grain input in a certain layer, the general result is that around this layer, the opacity is also increased, while in the remaining parts of the atmosphere, κ stays approximately unaffected. This means that around the layer of increased grain input, Brownian motion dominates, while far away, it is differential settling. This is the same general behavior as for an increase because of a higher fD/G,disk. The simulations of MP08 also show that the thickness of the layer with an increased opacity becomes smaller if the grains are injected in deeper layers. This is expected, because the timescale of growth decrease with depth (Fig. 4).

3.1.7. Implications for the metallicity effect and pebble accretion

The independence of the opacity on the grain accretion rate is the result of this study that has the most important implications for planet formation theory. It means that even if more grains are added at the top of the atmosphere, the κgr does not significantly increase, as one may naively expect. This is important in two contexts.

First, in the context of the metallicity effect, which is the observational finding that giant planets are clearly more frequent around host stars with a high [Fe/H] (Gonzalez 1997; Santos et al. 2001; Fischer & Valenti 2005). Within the core accretion theory, this can be explained if a higher stellar [Fe/H] means a higher associated [Fe/H] in the protoplanetary disk. This likely leads to a higher surface density of planetesimals (Kornet et al. 2005), which in turn leads to a faster growth of more massive cores that can trigger gas runaway accretion (e.g., Pollack et al. 1996). This mechanism means that planet population synthesis models built on the core accretion paradigm reproduce the metallicity effect (Ida & Lin 2004; Mordasini et al. 2012a), but one could also argue that at a higher stellar [Fe/H], not only the planetesimal surface density, but also the fD/G,disk should be higher, which could (via a higher opacity) be detrimental to giant planet formation. The analytical model (and the numerical simulations of MP08) show that this is not the case. A higher stellar [Fe/H] should not strongly increase the grain opacity in the protoplanetary atmospheres. This is an important insight for the core accretion theory. Furthermore, in this light, the a priori questionable assumption of past population synthesis calculations (e.g., Ida & Lin 2004; Mordasini et al. 2012a) that the opacity is independent of the disk [Fe/H] can be physically justified, at least as a first-order approximation.

Second, the independence is important in the context of the recently proposed mechanism that the accretion of small objects (pebbles instead of 100 km planetesimals) is an essential element for the formation of giant planet cores (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012). Such small objects make a rapid formation of the core possible thanks to their low random velocities and the large, drag-enhanced capture radius of the protoplanet. This mechanism could, however, suffer from a disadvantage. In contrast to 100 km planetesimals (which deposit most of their mass in the convective zone or even directly on the core), small objects will substantially increase the grain input in the outer layers of the atmosphere (as illustrated by Fig. 1). This could potentially increase the opacity, making giant planet formation harder. The study here shows that this is, at least to the first order, not the case. An accretion of the core via small bodies should not substantially increase the opacity. We have further tested this by re-calculating the opacity for the case of shallow deposition (Eq. (4)), adding grains in the top layer at a rate of 10-6M/yr (MP08). We found that the opacity again remains virtually identical. This means that the small bodies would not substantially increase the opacity, but they would increase the mean molecular weight of the envelope gas in the layers that are sufficiently hot for grain evaporation (see Fortney et al. 2013 for estimates of the mean molecular weight in protoplanetary atmosphere). A high mean molecular weight is known to reduce the critical core mass (Stevenson 1982; Hori & Ikoma 2011). The low random velocity, the large capture radius, the high mean molecular mass, and the low opacity thus conspire to facilitate giant planet formation.

We have found an absence of a strong dependency of κgr on the dust-to-gas ratio in the dominant differential settling regime. However, the quantity that enters the equations is not fD/G,disk alone, but fD/G,diskXY. This means that for a fixed dust-to-gas ratio, the opacity is also not explicitly dependent on the gas accretion rate. It is clear on the other hand that the quantities on which κgr does depend (local atmospheric pressure P and gravity g) and XY are not independent. The absence of an explicit dependency on XY still means that if P and g remain similar when XY changes, then also κgr remains similar. This could be part of the explanation why the global structure of κ remains similar for different core masses (Sect. 4.3.2).

3.1.8. Impact of the monomer size

thumbnail Fig. 9

MP08 comparison case. Impact of the monomer size on the opacity. Blue: 1 μm. Brown: 10 μm. Green: 100 μm. Dashed lines are from MP08, solid lines from the analytical model where the 1 and 10 μm cases lead to an identical opacity so that the lines lie on top of each other.

Open with DEXTER

As did MP08, we finally also investigate how the grain opacity changes if the assumed size of the monomers is varied. As mentioned (Sect. 2.5), in atmospheres where advection is not important as is typically the case, the size of the monomers enters the analytical model only if the size of the grains calculated from the timescales is smaller than amono. Figure 9 shows the opacity as a function of height for monomer sizes of 1, 10, and 100 μm. The results of both the analytical model and that of MP08 are shown.

For the analytical model, the opacity is identical for monomer sizes of 1 and 10 μm. This is simply because the size of the grains found from the timescale arguments is always larger than 10 μm (Fig. 5). Only for amono = 100μm, is there an impact, and over a large part of the outer atmosphere the opacity becomes radially nearly constant at a low value of about 0.05 cm2/g. This is quite similar to the result in the numerical model, even if the numerical value of κgr differs by roughly a factor of five. In the numerical model, there is also a difference between the 1 μm and 10 μm cases. The 10 μm case has a higher opacity and is clearly closer to the result of the analytical model than the amono = 1 μm case. This is not surprising; it is again a manifestation of the result that the analytical model overestimates the efficiency of grain growth in the outermost layers. Therefore, in the analytical model, the grains in the outer regions also have a size of 30 μm. This means that the purely local and monodisperse description of the grain dynamics by the comparison of timescale alone is not sufficient to capture all effects seen in the numerical model, where in particular the relative contributions from different sized grains matter. As explained by MP08, grains with a size of about 10 μm are more efficient at causing a high opacity than both smaller and larger grains. If amono is set to 10 μm instead of 1 μm in the numerical model, the opacity becomes higher, and the radial shape is now more similar to the analytical model.

The results illustrate that it is important to know the typical size of the grains when they are accreted into the protoplanet. If the grains in the protoplanetary disk are typically already larger than ~100 μm, this could reduce the opacity in the protoplanetary atmosphere. A coupling of the grain evolution in the disk and in the protoplanetary atmospheres should therefore be included in future planet formation models.

thumbnail Fig. 10

Core (solid line) and envelope mass (dash-dotted line) as a function of time for Jupiter in situ formation. The panels show the results for an initial planetesimal surface density of 10 (left), 6 (middle), and 4 g/cm2 (right). The thick blue lines are calculated with the analytical model for the grain opacity, while the thinner green lines are from the numerical model of MBPL10. The thin red lines are also from the analytical model, but without taking the advection regime into account.

Open with DEXTER

4. Simulations of Jupiter’s formation coupled with the analytical grain growth model

We now present the results of coupling the analytical grain evolution model with (giant) planet formation simulations. This is analogous to the work of MBPL10 where the numerical grain evolution model of Podolak (2003) and MP08 was coupled with the giant planet formation model of Pollack et al. (1996). Our primary interest is to compare the analytical and numerical results for the giant planet formation timescale and the radial structure of the envelope.

4.1. Core accretion model and initial conditions

The core accretion model used here was first introduced in Alibert, Mordasini & Benz (2004), while a detailed description can be found in Alibert et al. (2005) and Mordasini et al. (2012c). Here we only give a short summary. In Mordasini et al. (2012c), simulations can be found that resemble closely the ones presented here with the difference that in the earlier work, the grain opacity was found by scaling the ISM opacity. Our code is based on the core accretion paradigm (see, e.g., Mordasini et al. 2010, for an overview) and contains three computational modules that are very similar to the ones used in the Pollack et al. (1996) model. Whenever possible, we have chosen parameters and settings that are identical to the ones in MBPL10.

The first computational module calculates the core accretion rate Z. It uses the same prescriptions for the planetesimal random velocities as Pollack et al. (1996). The gravitational capture radius that includes the effect of the star is taken from Greenzweig & Lissauer (1992). The size of the planetesimals is 100 km. The solid core is assumed to have a fixed density of 3.2 g/cm3, in contrast to Mordasini et al. (2012b) where the core density is variable. All these setting are very similar or identical to those in MBPL10.

The second computational module calculates the trajectories of the impacting planetesimals as they fly through the protoplanetary envelope (Mordasini et al. 2006). This model is equivalent to the one of Podolak et al. (1988). It yields the drag enhanced capture radius of the protoplanet which can be significantly larger than the core radius (Inaba & Ikoma 2003). The second output from the impact model, namely the radial profiles of mass and energy deposition, are currently not used. Instead, all mass is directly added to the core using the sinking approximation for the calculation of the core luminosity. The enrichment of the gaseous envelope is thus not considered, while it could be very important for a reduction of the critical core mass (Hori & Ikoma 2011), because high enrichments of Z ≈ 0.8 are expected for planets with masses between 1 and 10 M (Fortney et al. 2013).

In our simulations, the only source of grains in the gr term is the accretion of grains with the gas in the top layer of the atmosphere, i.e., we assume a deep deposition of the planetesimals (Eq. (4)) for the reasons shown in Sect. 2.2.2. In MBPL10, dust grains are, in contrast, injected both by gas accretion and also planetesimal ablation, but MBPL10 we note that the planetesimals deposit their mass typically fairly deep in the gaseous envelope. In any case, based on the analytical results (Eqs. (35) and (41)), a weak (or no) dependency of κgr on gr is expected if the dominant grain growth mechanism is differential settling. This is indeed found to be the case, as will be discussed below. The grains are added assuming a dust-to-gas ratio in the disk fD/G,disk = 0.01 for all planetesimal surface densities as in MBPL10. We thus ignore potential correlations of fD/G,disk and the planetesimal surface density. We have also run non-nominal simulations that use a fD/G,disk of 0.001 and 0.1. We find that this only has a 8% effect on the formation timescale (see Sect. 4.2.3). The monomer size is 1 μm, but this only matters in the advection regime, which is itself not important (see Fig. 10); MBPL10 use a amono = 1.26μm.

The third module solves the internal structure equations for the gaseous envelope in the quasi-static 1D approximation (see Bodenheimer & Pollack 1986) using the SCvH EOS (Saumon et al. 1995). We simplify the equations by assuming that the luminosity L is constant with radius r inside the envelope and we solve the temporal evolution based on an energy conservation approach (Mordasini et al. 2012c). This does not significantly influence the results, since ∂L/∂r is very small in the radiative zones (see Mordasini et al. 2012c for details). The outer boundary conditions (radius, pressure, temperature) are the same as in MBPL10. The opacity is calculated with the analytical model described above, it thus contains the contributions from the grains, plus the molecular and atomic opacities from Freedman et al. (2008) and Bell & Lin (1994). The grain opacity can be calculated dynamically during a radial integration of the envelope structure based on the local properties in a layer and on gr during one timestep. No special iterations or recursion on the opacity in the last timestep are necessary. This means that the computational time is only insignificantly increased compared to the case where tabulated or fitted (ISM) opacities are used.

As in MBPL10, the planet forms in situ at 5.2 AU in a disk with a gas surface density of 700 g/cm2. Disk evolution, orbital migration (Alibert et al. 2005), and the interaction with other embryos (Alibert et al. 2013) are all neglected. The planetesimal accretion rate is calculated as in Pollack et al. (1996), therefore it is probably too high for 100 km planetesimals (e.g., Fortier et al. 2013). This means that the simulations are not intended to show a necessarily realistic formation scenario of Jupiter. In any case, our interest here is to compare the analytical results with the numerical model of MBPL10, and therefore to run simulations that are as similar as possible.

Finally, we consider three initial surface densities of planetesimals at the planet’s position: 10, 6, and 4 g/cm2, again as in MBPL10, calling the three simulations Σ10, Σ6, and Σ4, respectively.

4.2. Time evolution of the core and envelope mass

The basic evolution of the core and envelope mass as a function of time for in situ formation has been described in many previous works (e.g., Pollack et al. 1996). The main result is that there are three different phases which we here only summarize very briefly. During phase I, the solid core forms rapidly (because of the assumption of low planetesimal random velocities). The luminosity of the planet is high and is caused by the accretion of the planetesimals. The envelope mass is low. This first phase ends once the isolation mass (Lissauer 1993) is reached, so that phase II starts.

Phase II is a plateau phase that is characterized by a relatively slow increase in envelope and core mass, with the envelope growing somewhat faster than the core. The duration of the phase II depends on the grain opacity (see, e.g., Paper I), because in order to expand its feeding zone of planetesimals, the envelope mass needs to grow, and this is only possible if the gas already contained within the planet’s Hill sphere can cool and contract. Phase II is the longest phase, hence its duration also determines the duration of the overall formation process. Phase II ends at the so-called crossover point when the core and envelope mass are equal and phase III begins.

In this third phase, the gas accretion rate increases rapidly because the planet’s cooling timescale decreases with increasing mass, leading to runaway gas accretion. Shortly after the crossover point (typically a few 104 yrs for a 10 M core), the gas accretion rate found by solving the internal structure equations becomes larger than the limiting gas accretion rate at which the disk can supply gas to the protoplanet (typically on the order of 10-2M/yr). At this point, the planet detaches from the disk and contracts rapidly (Lissauer et al. 2009; Mordasini et al. 2012c; Ayliffe & Bate 2012).

This fundamental behavior can be seen in Fig. 10. It shows the core and envelope mass as a function of time for the three different surface densities. The results obtained with the model presented here and the results of MBPL10 are included. In our model, the protoplanet starts with a tiny initial mass of 0.01 M. We have set the embryo’s starting time in a way that the evolution of the core mass during phase I overlaps in both models. This first phase is very similar in both simulations, as exactly the same equations are integrated, and only the evolution of the core (but not of the envelope and thus the opacity) matter.

Table 3

Characteristics of the Jupiter in situ formation simulations at the crossover point.

In the ideal case, the predicted temporal evolution found with the analytical and numerical model would be identical. A crucial figure that can be used to quantify the degree of agreement is the time of crossover, because the time to reach the crossover point depends directly on the opacity. It is also an important quantity in determining the fundamental outcome of the planet formation process: if the time to reach crossover is shorter then the disk lifetime, then a gas dominated Jovian planet formation is likely to form, because the gas accretion after the crossover point occurs on a short timescale of 104 − 105 years which is one to two orders of magnitude shorter than a typical disk lifetime. In the other case, only a solid dominated Neptunian planet will form. It is therefore obvious that the crossover time as predicted with the analytical grain growth model is an important quantity if the model is used in planetary population syntheses. We therefore list in Table 3 the time of crossover, the crossover mass, and the luminosity at this point for the different simulations. Since phase I is nearly identical in both models, the variation in the time of crossover is given by a variation in the duration of phase II. This is the quantity that was already studied in Paper I.

4.2.1. The nominal Σ10, Σ6, and Σ4 simulations

The left panel of Fig. 10 shows the Σ10 simulation. For a full ISM grain opacity, a crossover time of about 6–7.5 Myr has been found in previous studies (Pollack et al. 1996; Hubickyj et al. 2005; Mordasini et al. 2012c). Taking into account the grain physics, the analytical model predicts a crossover time of only 0.81 Myr, while the numerical model of MBPL10 predicts 0.97 Myr. This corresponds to a difference of only about 16%. Given the simplifications in the analytical approach, this is a good agreement. It indicates that the analytical and numerical models predict grain opacities (or at least total atmospheric optical depths) that agree relatively well during all of phase II. In the following section, we will compare the radial structure of the envelope including the opacity. Figure 10 and Table 3 show that the crossover mass (16 M) is also very similar in the two models (difference of only 2%). This is expected (Paper I) because the crossover mass is independent of the opacity for calculations similar to those performed by Pollack et al. (1996). The luminosity is higher in the analytical model, which is expected given the shorter accretion timescale.

The middle panel of Fig. 10 shows the Σ6 case. At a full ISM grain opacity, a crossover time of more than 60 Myr is expected (see Paper I) which clearly exceeds protoplanetary disk lifetimes. With grain evolution, the crossover time is found to be 2.35 and 1.63 Myr in the analytical and numerical model, respectively. In contrast to the Σ10 case, here the analytical model leads to a crossover time that is longer than that found by MBPL10 by about 0.7 Myr or about 40%. While this is a larger difference, it is a consequence of the stronger sensitivity of the formation time on κ at a lower core mass. This can be seen by the fact that in the Σ10 simulation, the reduction factor relative to the full ISM opacity is only a factor of about 6, whereas here the reduction factor is about 30. This means that the details of the calculation of κ are more important. The crossover masses of about 7.6 M agree again very well, while the luminosity reflects the difference in the formation time.

The left panel finally shows the Σ4 simulation. The core mass is now very low, only about 3 M at isolation, so that we are studying gas accretion onto super-Earth-type planets. Such planets could be the progenitors of the recently observed low-mass, low-density planets (e.g., Rogers et al. 2011; Lissauer et al. 2013). At full ISM opacity, the crossover point would only be reached after several hundred Myr (Paper I), which means that the formation of giant planets would be excluded for such disk conditions. Comparing the analytical and numerical model, a similar tendency to the Σ6 case is seen, but the difference in the crossover time is now larger (7.1 instead if 3.6 Myr). This corresponds to a disagreement by a factor of 2, and since it is of the same order as typical disk lifetimes, it can mean that in some cases the analytical and numerical model would disagree on whether giant planet formation is possible or not. On the other hand, except for the works of Movshovitz et al. (2010) and Rogers et al. (2011), in basically all other planet formation calculations (like Mizuno et al. 1978; Stevenson 1982; Pollack et al. 1996; Papaloizou & Nelson 2005; Hubickyj et al. 2005; Tanigawa & Ohtsuki 2010; Levison et al. 2010; Hori & Ikoma 2011; Ayliffe & Bate 2012) the grain opacity was essentially an unconstrained parameter which was varied over two or even three orders of magnitude. This would translate into a variation of the crossover time of the same magnitude in the linear regime (see Paper I). In this sense, a difference of a factor 2 found from a first simple analytical model is still a relatively small difference, but one should keep in mind that the analytical model has a tendency to overestimate the grain opacity in the atmosphere of low-mass protoplanets (see Fig. 13). The reason for this is currently being investigated (Ormel & Mordasini, in prep.). One possible reason is that the analytical model overestimates the impact of the planet mass on the grain size. In the analytical model presented here, at crossover, the grains are in the Σ10 case a factor of about 3 to 6 (depending on height) bigger than in the Σ4 case. If the actual dependency is weaker, this could explain why the formation timescale is underestimated in the Σ10 case, but overestimated in the two other cases with a smaller core mass.

Table 4

Impact of fD/G,disk on the crossover time in Myr.

4.2.2. Impact of the advection regime

Figure 10 also shows simulations where the advection regime (Sect. 2.5.5) is not taken into account. One sees that this only appreciably affects the accretion after the crossover point. Thus, the crossover time is nearly identical in runs with or without the advection regime (see also Table 4). After crossover, the gas accretion rate increases rapidly, so that the advection timescale becomes shorter (Sect. 2.5.5) to an extent that in the outer atmospheric layers, it is less than the timescale of growth by differential settling. Therefore, the growth regime changes. This has the consequence that an outer layer of high opacity forms (see Fig. 11 below). This lengthens the time between the crossover point and the moment when the limiting accretion rate is reached. For the Σ6 case, for example, the time between crossover and the moment when the limiting gas accretion rate of 10-2M/yr is reached is about 0.38 and 0.29 Myr for the simulation with and without advection, respectively. In MBPL10, the time between crossover and reaching a similar limiting gas accretion rate is 0.24 Myr. This indicates that the effect of advection is indeed overestimated in the analytical model as discussed in Sect. 2.5.5. In any case, differences of about 105 yrs are usually not important for the global growth history and final mass of a planet, because it is one order of magnitude less than typical disk lifetimes.

thumbnail Fig. 11

Radial structure of the envelope in the Σ10 case shortly before crossover. The plot shows the opacity (blue), temperature (brown), and gas density (green) as a function of height. Solid lines are the results from this work with convective parts indicated by a thicker line. Dashed lines show the corresponding results from MBPL10.

Open with DEXTER

4.2.3. Impact of the disk dust-to-gas ratio and monomer size

In Sects. 2.5.3 and 2.5.4 we have found that opacity in the dominating differential settling regime does not depend on the dust accretion rate or dust-to-gas ratio in the disk. To quantify the impact of fD/G,disk on the overall formation history, we have recalculated the Σ10 and Σ4 cases also for other fD/G,disk besides the nominal 0.01. We find that the impact on the formation history is small. Table 4 shows the crossover times for different fD/G,disk for the Σ10 and Σ4 cases, and for simulations with/without the advection regime. For the Σ10 case with advection, the time of crossover is 0.88, 0.81, and 0.80 Myr for fD/G,disk = 0.001, 0.01, and 0.1. This means that a lower dust-to-gas ratio even leads to a slightly longer formation timescale. Similar effects were already seen by MP08. In all cases, the difference corresponds to just 9% of the formation timescale for a change of fD/G,disk by a factor of 100. This means that the dependency is very weak. We note that the analytical model tends to overestimate the efficiency of differential settling in the outermost layers (Sect. 3.1.2). If these layers are in reality dominated by Brownian motion, then at least some positive correlation of fD/G,disk and the formation timescale is expected. The same is true for the impact of shallow deposition. In the current analytical model, it remarkably only has a negligible effect in the Σ10 simulation, and even leads to a mild reduction of the crossover time in the Σ4 case as grains grow more massive. These large grains cause a Q ≈ 2 in the entire atmosphere, while otherwise Q ≈ 3 in the outer layers. The occurrence of Brownian motion (or potentially collisional grain fragmentation, see Sect. 4.3.4) could, however, significantly alter this finding.

The table also shows the consequences of fD/G,disk for the Σ4 case. We see that the impact is clearly larger, with a reduction of the crossover time in simulations without advection from 7.09 in the nominal case to 5.49 Myr for fD/G,disk = 10-4 (25% difference). This is due to the following: in the Σ4 case, at low fD/G,disk, the grains in the outer layers stay so small (~1 μm) that their extinction efficiency Q decreases to 0.1. The leads to a radially roughly constant low κ ≈ 0.07 cm2/g in the outer parts of the atmosphere (the specific value might, however, be affected by our simplistic way of calculating Q). For higher fD/G,disk, Q is always between 2 and 3. In the Σ10 case, the grains also become smaller at lower fD/G,disk as seen in Sect. 3.1.6. But in contrast to the Σ4 case, for this more massive planet they still remain large enough for Q to be in the geometrical limit. As the grain growth mode is differential settling, this means that fD/G,disk does not really affect κ (Eqs. (34), (40)), so that the crossover time is almost independent of it. Including advection slightly reduces the crossover time to 5.44 Myrs in the Σ4 simulation. This is due to a modification of the radial opacity structure by the same mechanism as observed for fD/G,disk = 10-4 in the MP08 comparison case (Sect. 3.1.6). The impact on the overall formation time is still not extremely large, as the reduction of κ only occurs in the outer layers. The density there is low, so that the associated reduction of the total optical depth is not very high.

We have also tested the implications of different monomer sizes amono. For the MP08 comparison case, it was found that increasing the monomer size from 1 μm to 10 μm has no impact, while setting it to 100 μm leads to a reduction of the opacity in the outer layers (Fig. 9). Simulating the Σ10 case with amono = 10 and 100 μm instead of the nominal 1 μm has similar consequences: for amono = 10μm, the formation is virtually the same as in the nominal case. For amono = 100μm, the time of crossover is shorter by about 10%. In particular, a higher gas mass for a given core mass is seen in phase I, because in phase I the grain size found from the analytical estimates can be relatively small in the outer layers (less than 100 μm)

4.3. Radial temperature, density, and opacity structure

The ISM opacity scalings presented in Paper I lead with a calibrated grain opacity reduction factor of 0.003 to formation timescales that agree with those obtained with the numerical model of MBPL10. However, a simple scaling of the ISM opacity suffers from the limitation that it does not lead to the same (or similar) radial structure of the opacity. Only the total optical depth is similar (see Fig. 7). This means that it is unclear if the reduction factors that were calibrated for certain initial conditions are also applicable for other situations.

In this respect, the analytical growth model has a significant advantage as it calculates the opacity dynamically based on microphysical considerations. It is therefore very interesting to compare what the radial structures of the envelopes look like, and how they compare with the results of MBPL10.

4.3.1. The Σ10 case: structure at crossover

Figure 11 shows the structure of the gaseous envelope in the Σ10 case at a moment slightly before the crossover point. The core mass is 15.98 M while the envelope mass is 14.47 M. The temperature, density, and opacity are shown. We note that in contrast to the situation in Sect. 3 where the density and temperature were externally given and only the opacity was then calculated, here the entire structure is the result of solving the coupled internal structure equations. This means that in contrast to the results in Sect. 3, T, ρ, and κ are self-consistent.

thumbnail Fig. 12

Temporal evolution of the radial profile of κ as a function of height (distance from the planet’s center) in the Σ10 simulation. The left end of the lines corresponds to the core-envelope boundary while the right end is the planet’s outer radius. Radiative layers are red, convective parts blue. The two numbers at the left give the corresponding core and envelope mass. The left panel shows the evolution during phase I. A sequence of 12 structures is shown, covering a time interval of 0.27 Myr. During this time the core mass increases from 0.6 M to 11.97 M while the envelope mass grows from 6 × 10-6M to 0.9 M. The right panel shows the evolution during phase II. A sequence of 15 structures is shown, covering a time interval of 0.32 Myr. During this time the core mass increases from 12.5 M to 16 M while the envelope mass grows from 1.9 M to 14.5 M. This last structure is the same as shown in Fig. 11 and thus occurs just before crossover. It is the only structure where the advection regime occurs in the top layers, leading to an increased opacity.

Open with DEXTER

The plot also shows the corresponding result of MBPL10 (their Fig. 4). In their simulation, the planet is at crossover and has a core mass of 16.09 M. One sees that the two models predict an internal structure that agrees well. In MBPL10 the temperature at the core-envelope boundary is 18 600 K, while we find a temperature of 18 040 K, corresponding to a difference of 3%. The temperature structure consists of an outer zone that is radiative and roughly speaking isothermal (temperature increase by a factor of ~2), and a large inner convective zone where the temperature increases strongly (there is a tiny inner radiative zone in the range of about 4–4.6 × 1010 cm, but it remains without visible consequences for ρ and T). In the density structure, the radiative and convective parts can also be distinguished.

For the opacity, the following layers can be seen, starting at the top: in the outermost parts of the atmosphere (R ≳ 3.6×1011 cm), the growth occurs in the advection regime. This is a regime that only occurs at/after crossover (Fig. 10). The grains remain relatively small, causing an opacity that is rather high, at least in the top layer where κ ≈ 3.6 cm2/g. This is comparable to typical ISM opacities. The opacity structure of MBPL10 shows a feature that might be related, but this is not entirely clear (see also Fig. 11 in D’Angelo et al. 2014). For radii between 7.5 × 1010 cm and 3.6 × 1011 cm, the growth mode is differential settling. This is the typical regime in the atmospheres. A slight change in the slope of κgr can be seen at R ≈ 1.2 × 1011 cm. This is where the drag regime changes from Epstein to Stokes drag. The grain opacity falls strongly with decreasing height (but is still larger than κgas) which is very important for the ability of the planet to cool as noted by MBPL10. This is clearly a consequence of the growth of the grains. The opacity reaches a global minimum of only about 7 × 10-4 cm2/g at R = 7.5 × 1010 cm. Within this radius, the opacity is given by the molecular opacities rather than the grains, because now κgr<κgas. This is the same behavior as seen in Fig. 7, and was already found by MBPL10. The structure of the opacity below R = 7.5 × 1010 cm is thus given by the Freedman et al. (2008) data. This leads in particular to a second local minimum in κ at about 4 × 1010 cm, but in general, the opacity now increases with depth. The results of MBPL10 are very similar in this part of the envelope. This is not surprising as they use the same molecular opacity tables.

The radius where the molecular opacities become larger than the grain opacity is farther out than the radius where the grains evaporate. This grain evaporation radius is found to be at about 4 × 1010 cm as estimated from the evaporation timescale (Sect. 2.6). This means that grain evaporation is not important for the resulting opacity, which is an interesting result. The radius where the molecular opacities become dominant over the grain opacity is also larger than the radius where convection first sets in (at about 6 × 1010 cm). This means that no convective regions exist where the opacity is given by the grains. This is important, since such regions would make a special treatment necessary (MBPL10). This result is probably related to the finding of MBPL10 that including or neglecting the effects of convection on the grain growth process leads to negligible differences for the planet’s evolution.

4.3.2. The Σ10 case: temporal evolution of κ(R)

Figure 11 shows the radial structure of κ only at one specific moment in time. In order to understand how general this structure is, we plot in Fig. 12 the opacity as a function of height (covering the entire gaseous envelope) at 27 different moments in time. The left panel shows the temporal evolution during phase I, while the right panels shows phase II. A large range of core and envelope masses (and of luminosities) is covered, as indicated in the figure.

thumbnail Fig. 13

Radial profile of the temperature, density, and opacity in the Σ4 simulation. Solid lines are from the analytical model with convective parts indicated by thicker lines. Dashed lines show the corresponding results of MBPL10. Left panel: when the core and envelope mass is 3.27 and 0.74 M, respectively (3.21 and 0.74 M in MBPL10). Right panel: at a moment just before crossover. The core and envelope mass are 4.18 and 4.17 M, respectively (4.09 M in MBPL10).

Open with DEXTER

We see that nevertheless, the general structure of κ is, roughly speaking, similar in all states. It is on the order of 1 cm2/g at the top, and then decreases to ~10-3 cm2/g with increasing depth. Within this minimum, the molecular opacities start to dominate, so that κ now increases as we move farther in. Somewhat within the κ minimum, the energy transport mechanism changes from radiative diffusion to convection. There are some variations around this general picture. Very low-mass planets, for example, have an outer convective zone. It is therefore possible that the grain opacity calculated in these parts is affected by the fact that we do not take convection into account in the analytical model. These low-mass planets are also found to have a relatively large outer part where the opacity is on the order of 1 cm2/g, followed by a rapid decrease in κ.

The last radial structure shown in the right panel of Fig. 12 is the same as in Fig. 11, i.e., it is at crossover. As mentioned, the outer part is radiative down to R ≈ 6 × 1010 cm. A second tiny inner radiative zone occurs in the range 4.0–4.6 × 1010 cm. It is clearly related to the local minimum in the molecular opacities in this region. Within 4.0 × 1010 cm, the envelope is again convective. This sequence of convective and radiative zones is very similar to Fig. 3 of MBPL10 (which is, however, for a smaller envelope and core mass). At crossover, MBPL10 do not find the detached convective zone, and the radiative-convective boundary is at about 3.7 × 1010 cm, corresponding to our inner radiative-convective boundary at 4.0 × 1010 cm. A comparison of the radiative and convective gradients in this region shows, however, that they are very similar, meaning that a small difference in κ can easily change the classification of a layer as being radiative or convective, but the temperature gradient on the other hand will remain nearly the same, independent of the energy transport mechanism.

In their Σ10 case (but not in the Σ4 simulation), MBPL10 find a convective zone in the outermost layers. In the layers where the advection regime applies, we also see a steepening of the radiative gradient because of the increased opacity, but it remains a factor of ~3 shallower than the convective gradient so that we do not have such a top convective zone. On the other hand, farther down where the differential settling determines the low opacity, the radiative gradient is shallower than the convective one by a factor of 30, so that there is no doubt that these deeper layers are radiative. This means that in the analytical model also, the formation of a top convective layer could, in principle, occur for relatively small changes of κ (or the EOS).

Furthermore, the following general properties are found. First, during phase II, the outer radiative zone is always, roughly speaking, isothermal (see Rafikov 2006). Second, the grain growth mechanism is always differential settling and not Brownian motion. This also holds for the Σ6 and Σ4 case, but we caution that the analytical model tends to overestimate the efficiency of differential settling in the outer layers, as mentioned in Sect. 3.1.2. Third, the grain size generally increases as the planet grows. During the early (later) stages of phase I, the grains have sizes on the order of 1 μm (10 μm) in the outermost layers, and of 10 μm (100 μm) in the inner atmospheric layers. During phase II, the grain size is on the order of 100 μm at the top, and 1000 μm (0.1 cm) at the place where the grains evaporate.

4.3.3. The Σ4 case

This simulation is particularly interesting because first, it corresponds to a surface density that is only a factor of 1.6 higher than in the MMSN of Hayashi (1981), and second, the low core mass of only 3–4 M means that this case is interesting also in the context of the formation of low-mass low-density planets which were recently discovered in large numbers (e.g., in the Kepler-11 system, Lissauer et al. 2013; see also Rogers et al. 2011).

Figure 13 shows the radial structure of the envelope at two moments in time as found by the analytical model and MBPL10. The global structure of the opacity is similar to the Σ10 case. The opacity is on the order of unity in the top layers and then decreases substantially to a minimum where κ is between a few times 10-4 cm2/g and 10-3 cm2/g. Within this minimum, the gas opacities dominate over the grains, and the opacity now increases rapidly towards the interior. Somewhat within the minimum, the energy transport changes from radiative diffusion to convection. The small inner radiative zone that was seen in the Σ10 case (and also the top convective layer in MBPL10) is not found here. The envelope therefore simply consists of a deep convective part where almost all envelope mass resides (98% at crossover) below a radiative atmosphere where the temperature only changes by a factor of ~2. This is in agreement with MBPL10.

The figure in the left panel shows the envelope when the core and envelope mass is 3.27 and 0.74 M, respectively. This is the same envelope mass as in MBPL10, while the core mass is slightly lower (3.21 M). Movshovitz et al. (2010) find a temperature at the core-envelope boundary of 6230 K, while we find 6290 K which is very similar. However, the plot also shows that the opacity in the atmosphere calculated by the analytical model is always larger than in the numerical model (by a varying factor of ~2–5). This is very likely the reason for the longer formation timescale found in Sect. 4.2.1. In the structure at crossover shown in the left panel, the agreement is better, but there is still a certain tendency of the analytical model to overestimate κ. For this structure, MBPL10 find a core/envelope interface temperature of 7680 K, while we have 7816 K corresponding to a difference of about 2%. The reason for the higher opacity is difficult to pinpoint without further information on the state of the grains in the numerical model (like the grain size or dust-to-gas ratio as a function of height), but will be further investigated in upcoming work (Ormel & Mordasini, in prep.).

4.3.4. Perfect sticking versus bouncing and fragmentation

An important simplification in the model is the one of perfect sticking. For grain growth in protoplanetary disks it is well known that collisions can in reality also lead to bouncing or even fragmentation (see the discussion in Sect. 3.1.4). It is therefore important to understand to what extent the assumption of perfect sticking is justified in protoplanetary atmospheres. Figure 14 shows an overlay of the relation of settling velocity versus grain mass predicted by the analytical model and the outcome of laboratory experiments on dust collisions (a simplified version of Fig. 11 in Güttler et al. 2010). Most of these experiments were conducted with dust aggregates made of spherical monodisperse SiO2 monomers of about 1 μm in size. The phase II of the Σ10 and Σ4 simulation are shown. The actual collision velocity will be a fraction of the settling velocity. This overlay allows to study the fundamental outcome of the collisions of the grains (growth, bouncing, fragmentation).

thumbnail Fig. 14

Overlay of the settling velocity-mass relation and the outcome of laboratory experiments on dust collisions. The red and blue lines show the settling velocity versus the grain mass in the radiative part of the atmosphere during phase II for the Σ10 and Σ4 case, respectively. For the Σ10 case, these are the same 12 structures shown in the right panel of Fig. 12. The black lines show in a simplified way the outcome of collision experiments (Güttler et al. 2010) in the collision velocity-grain mass plane. The actual outcome also depends on the grain porosity and the mass ratio of the projectile and target.

Open with DEXTER

The collision outcome depends, however, not only on the grain mass and collision velocity. The mass ratio of the projectile and target as well as the porosity of the two collision partners is equally important, which means that there are, in principle, eight different panels necessary to show the collisional outcomes for one fixed material type (see Güttler et al. 2010, for details). Figure 14 therefore only gives a rough overview of the general behavior. A general result is first that low velocity collisions of small bodies lead to growth (triangle in the lower-left corner labeled “Growth”). This regime is important for grain growth in phase I (not shown in the figure). For grain growth to fall into this regime also in phase II, the typical collisional velocities need to be approximately 1–10% of the settling velocity for the Σ4 case, and ~0.1 % for the Σ10 case. Second, if the collisional velocity is increased, bouncing becomes important (inverted triangle in the middle labeled “Bouncing, Growth (pP)”. However, growth can still occur if the collisions occur mainly between bodies of a different size that are both porous (represented as “pP”) and/or if small porous objects hit large compact ones (represented as “pC”). If the collision velocities are of the same order as the settling velocities, this regime is dominant for the Σ4 case. The larger grains in the Σ10 simulation partially fall in this regime, but also in the one labeled “Bouncing, Transition (pP)”. In this regime, the collisional outcome transitions from growth to bouncing also in the pP case. This is an indication that bouncing is more important in more massive protoplanets. Third, if the impact velocity is further increased (beyond ~100 cm/s in most regimes, ~350 cm/s in the pP case), fragmentation occurs, except in the pC regime. In the Σ10 case this regime is marginally reached at crossover in the outer layers as shown by the rightmost structure which is the one at crossover. After crossover, the settling velocities continue to increase, reaching ~600 cm/s at the onset of the disk limited gas accretion. For the Σ4 case such high velocities are only reached very shortly before the disk limited gas accretion. Together with grain advection, this could lead to an increased opacity in phase III, lengthening the duration of this phase.

In summary, we conclude that the collision velocities in protoplanetary atmospheres are likely too low for grain fragmentation to be important, at least before crossover. The suppression of further grain growth owing to bouncing could in contrast be important especially for more massive cores. The actual outcome will depend on the typical mass ratio of the colliding grains, their porosity, and composition. Further studies describing these quantities (Podolak 2003; Ormel et al. 2007; Zsom et al. 2010) are thus needed to clarify to what extent the assumption of perfect sticking is applicable.

5. Summary and conclusions

In this work, we have derived a first simple but complete analytical model for the opacity due to grains suspended in the atmosphere (outer radiative zone) of forming protoplanets. The grain opacity κgr is a poorly known but central quantity for planet formation as it controls the protoplanet’s gas accretion rate (e.g., Ikoma et al. 2000). This determines the final bulk composition of low-mass planets (their H/He mass fraction) as well as the time until runaway gas accretion occurs and a giant gaseous planet forms. Therefore, the magnitude of the grain opacity leads to potentially observable imprints in the planetary mass-radius relationship (Paper I).

Our general approach is similar to the one of Rossow (1978) and consists in the comparison of the timescales of the microphysical processes that govern the grain dynamics. We consider the following processes (Sect. 2.1): dust settling in the Epstein and Stokes drag regime, growth by Brownian motion coagulation and differential settling, grain evaporation, and the advection of grains due to the contraction of the gaseous envelope. With these timescales it is possible to determine the typical size of the grains in each atmospheric layer as well as the dust concentration. This allows us to calculate the opacity. We derive analytical expressions for the grain opacity separately in five different regimes (Sect. 2.5).

We find that the dominating growth regime is differential settling. In this regime, the grain opacity takes a simple form and is given as 27Q/ 8 in the Epstein drag regime and 2Q/ in the Stokes regime where H is the atmosphere scale height and ρ the gas density (Sects. 2.5.3 and 2.5.4). Typically, the grains grow much larger than the relevant wavelengths of radiation, therefore the extinction coefficient Q is simply equal to 2 under most circumstances (Sect. 2.7). These two equations are the most important direct results of this study. These expressions are in particular independent of the protoplanet’s grain accretion rate gr.

This first analytical model is simplified in many aspects (Sect. 2.9): important simplifications are a radially constant mass flux of grains through the atmosphere, the replacement of the grain size distribution by only one typical grain size per layer, the assumption of perfect sticking, and the approximation of the Rosseland mean opacity by the opacity evaluated at the wavelength associated with the maximum of the Planck function.

Despite this, we find in general fair and sometimes even surprisingly good agreement with the numerical results of Movshovitz & Podolak (2008) and Movshovitz et al. (2010), also on a quantitative level. Building on the model of Podolak (2003), these works determine the grain opacity with a complex but computationally expensive numerical model. Numerous comparisons of the analytical and numerical results can be found in Sect. 3. They also allow us to understand the limitations of the analytical model, like an overestimation of the efficiency of growth in the top layers at least when advection is not important (Sects. 3.1.2, 3.1.8), or the absence of a high opacity zone around a layer of high planetesimal ablation (3.1.6).

But in general, we confirm their most important finding that grain growth is an efficient process in protoplanetary atmospheres. Starting with μm sizes in the outermost layers, the grains can grow up to ~0.1 cm in the deeper layers (Sect. 3.1.8). Grain growth therefore leads to a characteristic radial opacity structure (Sects. 3.1.1, 4.3.1): in the outermost layers, κgr takes high ISM-like values (~1 cm2/g), but, as we move deeper into the atmosphere, it decreases sharply to very low values (~10-3 cm2/g). At this point, the grain opacity becomes so low that the grain-free molecular opacity of the gas becomes dominant. This general structure is found for many different core and envelope masses (Sect. 4.3.2). This radial structure differs from the one obtained from scaling the ISM opacity by one uniform (and a priori arbitrary) reduction factor (Sect. 3.1.5). This was, nevertheless, the approach taken in many previous works (e.g., Mizuno et al. 1978; Stevenson 1982; Pollack et al. 1996; Papaloizou & Nelson 2005; Hubickyj et al. 2005; Tanigawa & Ohtsuki 2010; Levison et al. 2010; Hori & Ikoma 2011; Ayliffe & Bate 2012; and Paper I).

A consequence of the grain growth is a total optical depth of the atmosphere that is much smaller than for a full ISM opacity, allowing rapid gas accretion. For example, in the Movshovitz & Podolak (2008) comparison case (Mcore = 12.6 M, Menv = 2.7 M), the optical depth is found to be about 90 and 50 in the numerical and analytical model, respectively (Sect. 3.1.5). For a full ISM grain opacity, it is, in contrast, about 4 × 104. For comparison, in a grain-free atmosphere, an optical depth of about 25 is found. This indicates that grain-free opacities provide an optical depth that is closer to the grain calculations than the full ISM value. In retrospect, this suggests that classical studies like Pollack et al. (1996) should have considered the grain-free case to be as equally meaningful as the full ISM opacity case. Potentially, this would have meant that the classical timescale problem of the core accretion theory would not have been seen as an important drawback of this paradigm for several decennia.

Analogous to the calculations presented in Movshovitz et al. (2010), we then coupled the analytical grain opacity model to our giant planet formation code. We studied the in situ formation of Jupiter for different planetesimal surface densities (Sect. 4.2.1) and compared the predicted formation timescales with the results of Movshovitz et al. (2010). For core masses of 16, 8, and 4 Earth masses, the analytical (numerical) model predicts a crossover time of 0.81 (0.97), 2.35 (1.63), and 7.09 (3.62) Myr. For comparison, at a full ISM grain opacity, the formation times would be about 6, 60, and 100 Myr, respectively. This means that the formation timescale predicted by the analytical model agrees very well with the numerical model for larger core masses. It is overestimated for low-mass cores, but, compared to the previous situation where not even the order of magnitude of κgr was analytically known, this is still an improvement.

We have studied the impact of different dust-to-gas ratios in the accreted gas on the grain opacity (Sects. 3.1.6, 4.2.3). As indicated by the aforementioned equations for κgr, we only found a very weak (or no) dependency, as already observed by Movshovitz & Podolak (2008). Our analytical results yield the physical explanation which is that a higher dust input leads to a higher dust-to-gas mass ratio (which increases the opacity), but also larger grains (which decreases the opacity). In contrast to Brownian coagulation, these effects cancel each other out in the dominating growth regime of differential settling. The dependency of κgr on different dust-to-gas ratios is therefore the litmus test to understand which growth process is dominant.

This independence of κgr on the grain accretion rate is the result of this study with the most important implications for planet formation (Sect. 3.1.7): first, it means that a higher stellar [Fe/H] (which presumably leads to more dust in the protoplanetary disk) should not strongly increase the grain opacity in protoplanetary atmospheres. This is an important insight for the core accretion theory. It means that a higher stellar [Fe/H] only favors giant planet formation as it presumably also leads to a higher surface density of planetesimals without being detrimental to it because of an increased opacity. This corroborates the result that the core accretion paradigm can explain the observed increase of the giant planet frequency with stellar [Fe/H].

Second, it is important for the recently suggested formation of giant planet cores by the accretion of pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012). Such small objects make high solid accretion rates possible. However, in contrast to classical 100 km planetesimals that deposit most of their mass in the convective zone or even directly on the core, such small objects strongly increase the grain input in the outer layers of the atmosphere because of ablation (Sect. 2.2.2, Appendices A, B). This could potentially increase the opacity, making giant planet formation harder. Our study shows that this is –at least to the first order– not the case. An accretion of the core via small bodies should not strongly increase the opacity. Instead, the low random velocity, the large capture radius, the resulting high mean molecular weight of the envelope, and the low opacity should conspire to facilitate giant planet formation with small bodies. We note, however, that higher grain opacities are expected if the protoplanet accretes bodies with a wide range of sizes: this leads to an input of small grains by ablation in all layers, which increases the opacity (Movshovitz & Podolak 2008; D’Angelo et al. 2014).

The analytical grain opacity can be incorporated into numerical models calculating the gas accretion at a computational cost that is only insignificantly higher than tabulated opacities. The contemporaneous work of Ormel (2014) also allows to estimate the grain opacity in an efficient fashion by adding a fifth equation to the normal planetary structure equations. This equation is numerically integrated in concert with the other four structure equations. In contrast to the analytical model here, it allows to take into account a radially varying grain input, and the effects of a bimodal grain size distribution. The basic conclusions obtained through this approach agree well with those found in the analytical model.

Therefore, in upcoming work, we will repeat our population synthesis calculations presented in Paper I, but instead of scaled ISM opacities, we will use physically motivated values of κgr. This will make it possible to get a better understanding of the role of grain opacity in controlling the planetary bulk composition, and the associated consequences for the planetary mass-radius relationship. This is very important in a period where it becomes for the first time possible to study observationally the bulk composition of many exoplanets (e.g., Lopez & Fortney 2014), and in particular the transition from solid to gas-dominated planets (e.g., Marcy et al. 2014).


1

In this work, “opacity” always indicates the Rosseland mean, except when otherwise stated.

2

This holds for strictly static calculations like in Mizuno (1980) or the corresponding analytical solution Stevenson (1982). In time dependent calculations as in Pollack et al. (1996), the crossover mass (which is the equivalent of the critical mass) is independent of κ, but the time until the crossover mass is reached depends on it (Paper I). This means that qualitatively, the implications are the same.

3

A similar dimensionless number is the Stokes number. It is usually defined in the context of turbulent flows (e.g., Birnstiel et al. 2012), while here we are dealing with a situation where in principle there is no background flow of gas (at least if we neglect the effect of envelope contraction). Therefore, we do not call this quantity directly the Stokes number despite its similar nature.

4

We find below that growth by differential settling is the dominant process. We include Brownian motion mainly for completeness.

5

One can also use H-1 = (dP/ dR)(1 /P) and dP/ dr = −ρg from the hydrostatic equilibrium.

6

During the preparation of this work we became aware of the results of D’Angelo et al. (2014). In their Fig. 11, such an effect can be seen for the most massive core they consider.

7

During the preparation of this paper, we became aware of the work of Cuzzi et al. (2014). In future work it will be possible to use the results for Q from this model instead of the less general approximations presented here.

Acknowledgments

I thank Th. Henning, G.-D. Marleau, P. Mollière, K.-M. Dittkrist, H. Klahr, Y. Alibert, W. Benz and C. W. Ormel for enlightening discussions and valuable input. I thank the Max-Planck-Gesellschaft for the Reimar-Lüst Fellowship. I also thank the referee Dr. Morris Podolak for an interesting and constructive report that helped to improve the manuscript.

References

Online material

Appendix A: Planetesimal mass deposition: impact geometry

thumbnail Fig. A.1

Left panel: trajectories of 100 km icy planetesimals in the envelope of the protoplanet in the Σ10 case at crossover for two different impact parameters. The critical case (red) and the mean case (blue) are shown. The three concentric black circles are the Hill sphere, capture, and core radius (from large to small). Right panel: zoom-in onto the central regions of the two trajectories. The color code gives the instantaneous mass loss rate along the trajectory. Almost all of the mass is deposited in the terminal explosion (red parts of the trajectory). The black circles are again the capture and core radius.

Open with DEXTER

The radial mass deposition profiles due to impacting planetesimals shown in Fig. 1 were calculated for an exact head-on impact geometry. However, it is clear that off-center impacts are much more likely than the (near) head-on case because of the larger collisional cross section. In this Appendix we investigate the consequences of non-central impact geometries on the mass deposition profile.

The dissipative effect of gas drag increases the collisional cross section of a protoplanet that has a gaseous envelope beyond the cross section that is caused by gravitational focussing only (Podolak et al. 1988; Inaba & Ikoma 2003). The largest impact parameter that leads to accretion is defined as the critical impact parameter pcrit. It results in a planetesimal trajectory with an apocenter after the first encounter that is equal to the protoplanet’s Hill sphere radius (Pollack et al. 1996). Such an orbit must inevitably lead to the accretion of the planetesimal (at least in the two-body approximation of protoplanet and planetesimal that is used here) while planetesimals with a larger impact parameter leave again the planet’s Hill sphere and return under the dominant gravitational control of the star. Assuming that planetesimals enter the protoplanet’s Hill sphere isotropically then leads to a mean impact parameter of pmean = pcrit/.

Figure A.1, left panel, shows the trajectories of 100 km planetesimals in the envelope of the protoplanet in the Σ10 simulation at the crossover point when the core and envelope mass is approximately 16 M. The core radius (about 0.27 R), the capture radius (the pericenter of the critical orbit on the first encounter, about 6.9 R), and the Hill sphere radius (about 348 R) are also shown.

In the critical case, the planetesimal is destroyed after about 17 orbits around the protoplanet. In the mean case, the planetesimal is destroyed already during the first encounter. This indicates that the annulus of impact parameters that lead to spiraling trajectories is only relatively narrow. The right panel of the figure shows a zoom-in onto the central regions of the two trajectories. The color of the lines shows the instantaneous mass loss rate.

One sees that the mass loss rate is on the order of 1016 g/s during the first few pericenter passages. However, as soon as the aerodynamical disruption starts (at the very end of the trajectories), the mass loss rate increases very rapidly by approximately four orders of magnitudes. This is caused by the rapid increase in the planetesimals’s cross section as the planetesimal is deformed by the aerodynamic load into a flat “pancake”, leading to a terminal explosion (Zahnle 1992; Chyba et al. 1993). This final destruction phase where more than 90% of the mass is deposited has a duration of only ~100 s. It is a violent process where an energy deposition rate of ~1031 erg/s is reached.

Figure A.2 shows the corresponding radial mass deposition profiles. The curve for the head-on case is the same as shown in Fig. 1. All three curves share the property that the mass loss rate is relatively small in the outer layers where only thermal ablation occurs, followed by a sharp upturn when the aerodynamic disruption starts. The details differ, however, and can be understood from the shape of the trajectories. One first notes that the aerodynamical disruption begins in the head-on and mean cases at approximately the same altitude (about 4 R). This is due to the fact that both planetesimals approach the planet’s core at a velocity v approximately equal to the escape velocity. This means that they are exposed to approximately the same stagnation pressure (1/2ρv2 where ρ is the atmospheric gas density). However, the head-on case penetrates about 0.5 R deeper, since its velocity vector only has a radial component towards the core, while the mean case has a substantial tangential (along-orbit) component. The critical case differs in the sense that its velocity is given by the Keplerian velocity around the planet and not the escape velocity, i.e., it is about a factor smaller. This means that the dynamic pressure becomes equal to the planetesimal’s tensile strength only in a deeper layer. Because of this, the terminal explosion occurs at a somewhat lower altitude (about 3.5 R), but since the radial component of the velocity vector is only very low, the radial range over which the mass is deposited is very thin, as can be seen in Fig. A.2. In this mass deposition profile, one also clearly sees the consequences of the orbits around the planet before the final destruction. It leads to an increase in the mass deposition of about one, and, at special radii, two orders of magnitude. There is for example a local maximum at about 7 R. Clearly, this corresponds to the region around the capture radius, where the planetesimal makes many orbital passes.

It is interesting to compare the cumulative mass deposition profile with the altitude Revap where the ambient temperature of the atmosphere is sufficiently high to evaporate silicate grains. This altitude occurs at about 5 R in the deep deposition case (see the “evap” line in Fig. 1). We neglect here that this distance is itself dependent on the opacity and thus the mass deposition profile. Planetesimal material that is deposited above this height can contribute to grain formation. We find that about 0.2% and 8% of the planetesimal’s initial mass is deposited above this distance in the mean and critical impact geometry, respectively. At crossover, the gas accretion rate is about a factor of 5 higher than the solid accretion rate (Pollack et al. 1996). Assuming an fD/G,disk = 0.01, this means that the ratio of the effective grain accretion rate due to planetesimals relative to the grain accretion rate together with the gas is about 0.04 and 1.6, again for the mean and critical geometry. We conclude that in the typical impact geometry, grain accretion together with the gas is clearly dominant so that the deep deposition assumption is indeed justified, while for the critical geometry both sources have approximately the same importance in the lower atmospheric layers. It is clear that we have here only analyzed an individual case. A systematic study is certainly warranted for future work.

thumbnail Fig. A.2

Consequence of the impact geometry on the radial mass deposition profile of 100 km icy planetesimals. The mass deposition per length as a function of height above the planet’s center is shown for a head-on collision (same line as in Fig. 1) and for the critical and mean impact parameter (corresponding to the trajectories shown in Fig. A.1).

Open with DEXTER

Appendix B: Planetesimal mass deposition: planetesimal composition

In this section we investigate the consequences of the planetesimal composition (ices versus silicates) for the radial mass deposition profile. One expects that icy planetesimals deposit much more mass in the upper layers compared to rocky impactors (Podolak et al. 1988). In this section we confirm this general result, but we also find that the extent and cause of the difference can vary substantially depending on the impactor size.

As in the previous section, we study the radial mass deposition profile in the Σ10 case at the crossover point. We note that the atmospheric structures at this moment are not identical for different planetesimal sizes and compositions as the planet’s accretion history depends on these quantities (e.g., Pollack et al. 1996). Figure B.1 shows the radial profiles for planetesimals made of silicates for a head-on and critical impact geometry and for initial planetesimal sizes of 100 km, 1 km, and 10 m. The curves for the last two were scaled such that the total mass deposited is the same as in the 100 km case. The profiles of head-on icy impactors, the same as in Fig. 1, are also shown for comparison.

thumbnail Fig. B.1

Radial mass deposition profiles of 100 km (left), 1 km (center), and 10 m planetesimals (right panel). The 1 km and 10 m curves are scaled so that the total mass deposition is the same as in the 100 km case. The planetesimal material and impact geometry is indicated in the plot. The head-on impact of icy planetesimals is the same as shown in Fig. 1. The curves for 1 km and 10 m rocky planetesimals lay partially on top of each other for the two impact geometries.

Open with DEXTER

Three material properties are mainly responsible for the differences: the density (set to 1 and 3.2 g/cm2 for ice and silicate, respectively), the tensile strength (4×106 and 3.5 ×108 dyn/cm2, Svetsov et al. 1995), and the surface temperature where thermal ablation becomes important. This temperature depends on the vapor pressure law and is in the range 1300–2000 K for silicates (Opik 1958) and 180–220 K for ices (Iseli et al. 2002).

The left panel shows the mass deposition of 100 km impactors. For both compositions, the impact ends in a terminal explosion. We see that it occurs for rocky planetesimals at about 2 R instead of about 3 R as for the icy impactors. This is a direct consequence of the different tensile strengths. The difference between the head-on and critical geometry is quantitatively equivalent as found in Fig. A.2 for icy planetesimals. One also notes that in contrast to icy bodies, no mass at all is lost outside of about 10 R. This is a direct consequence of the higher temperature necessary for thermal ablation. For such large bodies, the main energy input driving ablation is shock wave radiation since before the terminal explosion, such big bodies always travel at a hypersonic velocity. The terminal explosion happens far below the altitude where the background atmosphere is hot enough to vaporize silicate grain (Revap ≈ 5 R), and pure thermal ablation is inefficient for such large bodies. This is due to the low heat transfer coefficient and the low surface-to-mass ratio. Therefore, only about 0.02% and 2% of the initial mass of the planetesimal could potentially contribute to grain formation. Grain accretion together with the gas is therefore clearly the dominant source in this scenario and the deep deposition scenario applies.

The case of 1 km planetesimals (middle panel) is different. Here, the change of the material properties (namely the higher tensile strength) leads to a fundamentally different character of the impact. For icy composition, the 1 km planetesimal gets aerodynamically disrupted as is the case for the 100 km impactors, so that it does not penetrate deeper than 3–4 R. Initially, the rocky impactor also travels at a hypersonic velocity. Therefore the thermal mass loss starts at the same height as in the 100 km case (10 R), but due to the lower surface-to-mass ratio, the 1 km planetesimal undergoes a stronger deceleration at higher altitudes. This means that it never flies at a very high velocity in the dense lower layers. Therefore, the dynamic pressure never overcomes the tensile strength of the silicate impactor (but it does so for ice). The silicate impactor thus does not undergo aerodynamic disruption. Instead, it is slowed down to terminal velocity where the drag and gravitational force are in equilibrium. It then sinks down at subsonic velocity, undergoing thermal ablation only due to the ambient radiation. Pure thermal ablation is still relatively inefficient for a 1 km sized body. Therefore, this impactor penetrates surprisingly deep into the envelope, namely to about 0.7 R, much deeper than larger (and smaller) bodies. Such 1 km silicate planetesimals thus belong to the class of intermediate-sized objects that are too small to undergo violent aerodynamical fragmentation, but too large to be efficiently thermally ablated. They therefore have a surprisingly large penetration depth and form a “core-hit finger” (see Mordasini et al. 2006).

The critical impact makes about nine orbits around the protoplanet before destruction. This is, however, not visible in the mass deposition profile as the temperature is too low in this part of the orbit. Ablation only sets in once the planetesimal is slowed down to terminal velocity and then radially sinks towards the core. Therefore, the mass deposition profile is nearly identical to the head-on case.

The 10 m impactor is shown in the right panel of Fig. B.1. One sees that the mass deposition profile of silicate impactors significantly differs from that of icy bodies. Mass loss only occurs in a small annulus between 4–5 R. Such small bodies are slowed down to terminal velocity relatively high in the atmosphere (both for a head-on and critical impact geometry) and then radially sink down undergoing thermal ablation due to ambient radiation. This is the same fundamental type of impact without aerodynamic disruption as for 1 km silicate planetesimals. The difference is that the 1 km body is initially heated by shock wave radiation and thus also loses some mass farther up in the atmosphere (but still only within 10 R). The 10 m object in contrast starts to lose mass only at about 5 R where the undisturbed atmosphere become hot enough. It is not a coincidence that this is the same altitude as where silicate dust grains start to evaporate (Revap). This means that this small silicate body does not contribute any material that can form grains (at least in our approximation that no mechanical mass loss occurs at aerodynamic loads smaller than the tensile strength). This is an interesting result. It indicates that for protoplanets inside the iceline that accrete pebbles, accretion of grains together with the gas is the only source of material for grain opacity. Thus, again the deep deposition approximation would apply here.

The large diversity of planetesimal impact scenarios in terms of destruction mechanisms, penetration depths, or the importance of the impact geometry and the implications for the mass deposition profile will be studied in a systematic way in future work.

Appendix C: Determination of the grain size including the gas-grain relative velocity

In the analytical model, we consider two limiting cases, namely that the gas envelope is completely static (settling regime), or that the grains move at the full velocity of the gas (advection regime). This has the advantage that the grain size and thus opacity can be determined completely analytically. In this paragraph we demonstrate how the relative velocity of the gas to the grains can be included directly in the determination of the grain size. The disadvantage is that the grain size can then only be determined by a numerical root search in each layer. On the other hand, it is then also possible to smoothly combine the two drag regimes into one single equation.

The comparison of the equations for the drag force in the Epstein and Stokes regimes shows that the two can be combined into one expression that reduces to the two separate expressions in the limits of very high or low Knudsen numbers (see Nayakshin 2010). Writing this drag law in terms of the relative velocity vvgas and equating with the gravitational force then yields terminal grain velocity vtrav. It is given as (C.1)From this we get the timescale on which grains traverse one scale height, τtrav = H/vtrav. This expression is valid both in the Epstein and Stokes regimes, and for settling and advection that had to be treated separately above. One thus has (C.2)For vgas = 0, this reduces to Eqs. (11) and (15) in the appropriate limits of the Knudsen number, while for a sufficiently high gas velocity it approaches τadv = H/vgas.

thumbnail Fig. C.1

Left panel: grain size as a function of height for different growth mechanisms in the MP08 comparison case. The plot shows the size caused by growth by Brownian motion coagulation (dotted) and differential settling (solid) as predicted by the analytical expressions and by solving numerally for the grain size. Right panel: corresponding grain opacity.

Open with DEXTER

Combining the expression for the growth timescale in the Brownian motion regime (Eq. (19)) with the expression for vtrav and assuming again a radially constant grain flux yields the growth timescale (C.3)For vgas = 0, this equation reduces to τcoag,S (Eq. (20)) for small grains and τcoag,E (Eq. (25)) for large ones, while for high gas velocities it approaches Eq. (43).

For growth by differential settling, we can combine the growth timescale in the Epstein and Stokes regime into the equation (C.4)with Eqs. (5), (7), and (C.1) this becomes (C.5)It is straightforward to show that this expression reduces to Eq. (31), (37), (47), and (50) in the appropriate limits.

To determine the grain size, one again sets Eq. (C.2) equal to Eq. (C.3) and (C.5) and numerically determines the grain size a in the two regimes. Once it is determined, the process with the shorter timescale is chosen and the other quantities like in particular the opacity can be calculated in an analogous way as in the analytical model.

The left panel of Fig. C.1 shows the grain size in the atmosphere of the MP08 comparison case (Fig. 3), comparing the size as found by the analytical expressions and the numerical solution for a. In the entire atmosphere, growth by differential settling is found to occur on a shorter timescale than by Brownian motion (Fig. 4), leading to grains that are about one order of magnitude larger. The grain size found analytically for differential settling is therefore the same as shown in Fig. 8. In this atmosphere, the advection regime as defined by a τadv that is smaller than the growth timescale for vgas = 0 does not occur. Therefore, advection is completely neglected in the analytical model. The figure shows that this leads to grain sizes that are slightly overestimated relative to the result when numerically solving for a. Here the small, but non-zero vgas = 0 is properly taken into account with the consequence that grains cross one scale height on a somewhat shorter timescale. This reduces their grain size. We also see that the transition from the Epstein to the Stokes regime (at about 7 × 1010 cm) occurs smoothly in the numerical solution, while at least for Brownian coagulation, a small kink is visible in the analytical model. Both these differences between the analytical and numerical models are expected.

The right panel of the figure shows that the effect on the resulting opacity (that given by growth due to differential settling) is, however, very small. Only a small increase occurs in the outer layers which does not contribute much to the total optical depth. This justifies our neglecting of advection in the analytical model for this atmosphere. The grain opacity in the upper layers as predicted from Brownian motion coagulation is lower in the numerical model. This comes from the decrease of Q for smaller grains. In the inner parts, the opacity is almost identical. This is a consequence of the fact that the non-zero vgas not only decreases the grain size, but also reduces the grain concentration.

To understand the global effects, we have recalculated the Σ10 and Σ4 cases, determining a numerically. We found that the difference is very small. For the Σ10 case, the crossover time is somewhat longer relative to the analytical model as expected, but only by about 10%. For the Σ4 simulation, the difference virtually vanishes (less than 1%). This shows again that advection is more important for more massive cores. It also means that it is not necessary to numerically solve for a, at least before crossover and for cores that are not much more massive than in the Σ10 case.

Appendix D: Differential settling and Brownian motion in the advection regime

Figure D.1 shows the grain sizes a function of height at crossover in the Σ10 case. The density and temperature of MBPL10 shown in Fig. 11 is used to calculate the grain sizes together with a gas accretion rate of XY = 1.8 × 10-4M/yr. In this atmosphere, the analytical model finds that the growth timescale by differential settling τcoal,E becomes longer than the advection timescale τadv outside of about 3.8 × 1011 cm. This is indicated by a vertical gray line. Grain growth by differential settling neglecting the effect of envelope contraction (blue dashed line) leads to quite large grains of about 0.01 cm. In contrast, grain growth by Brownian motion coagulation in the advection regime (blue dash-dotted line) is inefficient, so that grains remain at their smallest allowed size of 1 μm. The black line shows the combined interpolated grain size. The brown solid line shows the grain size if no lower limit of 1 μm is imposed. Then the grains would have an even slightly smaller size of 0.7 μm.

The green solid line shows the grain size according to Eq. (49) for growth by differential settling in the advection regime. It shows that at least according to the description in Sect. 2.5.6, it does not seem to be an extremely efficient process. While it does predict larger grain sizes, the difference to Brownian motion is only about a factor of 2 or 3. In the normal regime (no advection), the grains formed by differential settling (coalescence) are, in contrast, typically one order of magnitude larger than for Brownian motion. This means that even if coalescence is included in the advection regime, there is no large effect (reduction) of the opacity predicted by the analytical model in the outermost layers, at least with the description used here. While the effect of advection before runaway is small, we will nevertheless investigate this regime further (Ormel & Mordasini, in prep.).

thumbnail Fig. D.1

Grain size as a function of height for different growth mechanisms at crossover in the MBPL10 case (initial planetesimal surface density of 10 g/cm2).

Open with DEXTER

Equation (49) gives an estimate of the grain size because of differential settling in the advection regime involving an integral. For completeness, we also give an analytical expression in this regime using that the temperature is roughly constant in the outermost layers (see Figs. 3 and 11) while the density decreases rapidly towards the exterior roughly as r-4. One can then calculate the integral (which gives the column density) to find (D.1)The grain size obtained with this equation is also shown in Fig. D.1 with the green dashed line where it provides a good approximation to the numerical result. It is, however, clear that this equation is not generally applicable, as it requires an a priori knowledge of the temperature and density structure.

All Tables

Table 1

Settings for the comparison with MP08.

Table 2

Grain size a in cm in MP08 and the analytical model as a function of height in the atmosphere.

Table 3

Characteristics of the Jupiter in situ formation simulations at the crossover point.

Table 4

Impact of fD/G,disk on the crossover time in Myr.

All Figures

thumbnail Fig. 1

Mass deposition profile (ablated mass per unit length) as a function of height for 100 km planetesimals (solid blue line). The brown dash-dotted line is the mass deposition for 1 km planetesimals, multiplied by a factor 106, while the green dashed line shows the result for 10 m planetesimals, multiplied by a factor 1012, so that the total mass deposited is identical in the three cases. The gray vertical lines show the radius where the envelope temperature becomes so high that silicate grains evaporate (“evap”) and where the molecular opacities alone become larger than the opacity due to dust in the deep deposition case (“κ-DG”). The steps in the right part of the lines are a numerical artifact without physical meaning.

Open with DEXTER
In the text
thumbnail Fig. 2

Extinction efficiency Q as a function of xmax = 2πa/λmax. The coefficient is shown for grain sizes of a = 1, 10, and 100 μm, computed with Mie theory (solid lines) and with the fit of Eq. (54) (dashed lines). In all three cases, λmax corresponds to the wavelength of the maximum in the blackbody spectrum.

Open with DEXTER
In the text
thumbnail Fig. 3

MP08 comparison case. Left panel: atmospheric gas density (solid line) and temperature (dashed) as a function of height (distance from the planet’s center) in the protoplanet’s atmosphere (outer radiative zone). Right panel: opacity as a function of height in the numerical model of MP08 (dashed line: grain opacity only) and in the analytical model of this work (thick solid line: combined grain and gas opacity). The thin blue solid (dashed) line is the gas (grain) opacity alone.

Open with DEXTER
In the text
thumbnail Fig. 4

MP08 comparison case. Left panel: timescale of growth by differential settling/coalescence (blue solid line) and Brownian coagulation (green dash-dotted). These timescales are by construction equal to the associated settling timescale. The timescale of gas advection is shown by the brown dashed line. In this atmosphere, the evaporation timescale of the grains is always many orders of magnitude longer because of the low temperatures. The thin gray dashed and dash-dotted lines show the result if the velocity of the gas is directly included to calculate the growth timescales of differential settling and Brownian coagulation, respectively (Appendix C). Right panel: ratio of the stopping to the settling timescale RSS (solid) and stopping timescale (dashed line) of the grains.

Open with DEXTER
In the text
thumbnail Fig. 5

MP08 comparison case. Left panel: grain radius (solid) and dust-to-gas ratio (dashed line). Right panel: Knudsen (solid) and Reynolds number (dashed line) of the grains as a function of height. The horizontal line indicates the critical Knudsen number of 4/9 where the drag regime changes.

Open with DEXTER
In the text
thumbnail Fig. 6

MP08 comparison case. Left panel: settling and by construction identical growth timescale in years (solid line) and settling velocity (dashed line) as a function of height. The thin dash-dotted line is the settling velocity found by MP08. As they only give the relative velocity vset(R) /vset(Rout), we have normalized their curve to the value predicted by the analytical model at Rout. Right panel: dust mass density (fD/G × ρ) found by the analytical model (blue solid line) and by MP08 (thin dash-dotted line). The two brown lines are the extinction efficiency, calculated with the fit of Eq. (54) (dashed) and directly with Mie theory (dash-dotted line).

Open with DEXTER
In the text
thumbnail Fig. 7

MP08 comparison case. Left panel: opacity as a function of height. The brown dashed is the result of MP08. Other lines show the total opacity found in this work (blue solid line), the gas opacity only (green dash-dotted line), the full ISM opacity (thin black solid line) and the ISM opacity × 0.003 (thin black dashed line). Right panel: corresponding optical depths.

Open with DEXTER
In the text
thumbnail Fig. 8

MP08 comparison case. Left panel: opacity as a function of height for eight different dust-to-gas ratios in the accreted gas fD/G,disk indicated in the plot. Right panel: corresponding size of the grains (solid lines) and local dust-to-gas ratio (dashed lines) for the four dust-to-gas ratios in the accreted gas indicated in the plot.

Open with DEXTER
In the text
thumbnail Fig. 9

MP08 comparison case. Impact of the monomer size on the opacity. Blue: 1 μm. Brown: 10 μm. Green: 100 μm. Dashed lines are from MP08, solid lines from the analytical model where the 1 and 10 μm cases lead to an identical opacity so that the lines lie on top of each other.

Open with DEXTER
In the text
thumbnail Fig. 10

Core (solid line) and envelope mass (dash-dotted line) as a function of time for Jupiter in situ formation. The panels show the results for an initial planetesimal surface density of 10 (left), 6 (middle), and 4 g/cm2 (right). The thick blue lines are calculated with the analytical model for the grain opacity, while the thinner green lines are from the numerical model of MBPL10. The thin red lines are also from the analytical model, but without taking the advection regime into account.

Open with DEXTER
In the text
thumbnail Fig. 11

Radial structure of the envelope in the Σ10 case shortly before crossover. The plot shows the opacity (blue), temperature (brown), and gas density (green) as a function of height. Solid lines are the results from this work with convective parts indicated by a thicker line. Dashed lines show the corresponding results from MBPL10.

Open with DEXTER
In the text
thumbnail Fig. 12

Temporal evolution of the radial profile of κ as a function of height (distance from the planet’s center) in the Σ10 simulation. The left end of the lines corresponds to the core-envelope boundary while the right end is the planet’s outer radius. Radiative layers are red, convective parts blue. The two numbers at the left give the corresponding core and envelope mass. The left panel shows the evolution during phase I. A sequence of 12 structures is shown, covering a time interval of 0.27 Myr. During this time the core mass increases from 0.6 M to 11.97 M while the envelope mass grows from 6 × 10-6M to 0.9 M. The right panel shows the evolution during phase II. A sequence of 15 structures is shown, covering a time interval of 0.32 Myr. During this time the core mass increases from 12.5 M to 16 M while the envelope mass grows from 1.9 M to 14.5 M. This last structure is the same as shown in Fig. 11 and thus occurs just before crossover. It is the only structure where the advection regime occurs in the top layers, leading to an increased opacity.

Open with DEXTER
In the text
thumbnail Fig. 13

Radial profile of the temperature, density, and opacity in the Σ4 simulation. Solid lines are from the analytical model with convective parts indicated by thicker lines. Dashed lines show the corresponding results of MBPL10. Left panel: when the core and envelope mass is 3.27 and 0.74 M, respectively (3.21 and 0.74 M in MBPL10). Right panel: at a moment just before crossover. The core and envelope mass are 4.18 and 4.17 M, respectively (4.09 M in MBPL10).

Open with DEXTER
In the text
thumbnail Fig. 14

Overlay of the settling velocity-mass relation and the outcome of laboratory experiments on dust collisions. The red and blue lines show the settling velocity versus the grain mass in the radiative part of the atmosphere during phase II for the Σ10 and Σ4 case, respectively. For the Σ10 case, these are the same 12 structures shown in the right panel of Fig. 12. The black lines show in a simplified way the outcome of collision experiments (Güttler et al. 2010) in the collision velocity-grain mass plane. The actual outcome also depends on the grain porosity and the mass ratio of the projectile and target.

Open with DEXTER
In the text
thumbnail Fig. A.1

Left panel: trajectories of 100 km icy planetesimals in the envelope of the protoplanet in the Σ10 case at crossover for two different impact parameters. The critical case (red) and the mean case (blue) are shown. The three concentric black circles are the Hill sphere, capture, and core radius (from large to small). Right panel: zoom-in onto the central regions of the two trajectories. The color code gives the instantaneous mass loss rate along the trajectory. Almost all of the mass is deposited in the terminal explosion (red parts of the trajectory). The black circles are again the capture and core radius.

Open with DEXTER
In the text
thumbnail Fig. A.2

Consequence of the impact geometry on the radial mass deposition profile of 100 km icy planetesimals. The mass deposition per length as a function of height above the planet’s center is shown for a head-on collision (same line as in Fig. 1) and for the critical and mean impact parameter (corresponding to the trajectories shown in Fig. A.1).

Open with DEXTER
In the text
thumbnail Fig. B.1

Radial mass deposition profiles of 100 km (left), 1 km (center), and 10 m planetesimals (right panel). The 1 km and 10 m curves are scaled so that the total mass deposition is the same as in the 100 km case. The planetesimal material and impact geometry is indicated in the plot. The head-on impact of icy planetesimals is the same as shown in Fig. 1. The curves for 1 km and 10 m rocky planetesimals lay partially on top of each other for the two impact geometries.

Open with DEXTER
In the text
thumbnail Fig. C.1

Left panel: grain size as a function of height for different growth mechanisms in the MP08 comparison case. The plot shows the size caused by growth by Brownian motion coagulation (dotted) and differential settling (solid) as predicted by the analytical expressions and by solving numerally for the grain size. Right panel: corresponding grain opacity.

Open with DEXTER
In the text
thumbnail Fig. D.1

Grain size as a function of height for different growth mechanisms at crossover in the MBPL10 case (initial planetesimal surface density of 10 g/cm2).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.