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/00046361/201423702  
Published online  08 December 2014 
Grain opacity and the bulk composition of extrasolar planets
II. An analytical model for grain opacity in protoplanetary atmospheres^{⋆}
MaxPlanckInstitut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
email: mordasini@mpia.de
Received: 24 February 2014
Accepted: 25 June 2014
Context. We investigate the grain opacity κ_{gr} in the atmosphere (outer radiative zone) of forming planets. This is important for the observed planetary massradius relationship since κ_{gr} affects the primordial H/He envelope mass of lowmass planets and the critical core mass of giant planets.
Aims. The goal of this study is to derive a simple analytical model for κ_{gr} and to explore its implications for the atmospheric structure and resulting gas accretion rate.
Methods. Our model is based on the comparison of the timescales of the most important microphysical processes. We consider grain settling in the Stokes and Epstein drag regime, growth by Brownian motion coagulation and differential settling, grain evaporation in hot layers, and grain advection due to the contraction of the envelope. With these timescales and the assumption of a radially constant grain flux, we derive the typical grain size, abundance, and opacity.
Results. We find that the dominating growth process is differential settling. In this regime, κ_{gr} has a simple functional form; it is given as 27Q/ 8Hρ in the Epstein regime in the outer atmosphere and as 2Q/Hρ for Stokes drag in the deeper layers. Grain growth leads to a typical radial structure of κ_{gr} with high ISMlike values in the outer layers but a strong decrease towards the deeper parts where κ_{gr} becomes so low that the grainfree molecular opacities take over.
Conclusions. In agreement with earlier results, we find that κ_{gr} is typically much lower than in the ISM. In retrospect, this suggests that classical giant planet formation models should have considered the grainfree case to be as equally meaningful as the full ISM opacity case. The equations also show that a higher dust input in the top layers does not strongly increase κ_{gr}. This has two important implications. First, for the formation of giant planet cores via pebbles, there could be the adverse effect that pebbles tend to increase the grain input high in the atmosphere because of ablation. This could in principle increase the opacity, making giant planet formation difficult. Our study indicates that this potentially adverse effect is not important. Second, it means that a higher stellar [Fe/H] which presumably leads to a higher surface density of planetesimals only favors giant planet formation without being detrimental to it because of an increased κ_{gr}. This corroborates the result that core accretion can explain the observed increase of the giant planet frequency with stellar [Fe/H].
Key words: opacity / planets and satellites: formation / planets and satellites: atmospheres / planets and satellites: interiors / planets and satellites: individual: Jupiter / methods: analytical
Appendices are available in electronic form at http://www.aanda.org
© 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 superEarth to Jovian masses. This allows us to investigate statistically how the bulk composition depends on fundamental planetary parameters like the mass or the semimajor 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 KelvinHelmholtz 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 massradius 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 KelvinHelmholtz 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 massradius 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 massradius 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 lowmass 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 opacity^{2}. For example, as already found in the early work of Mizuno (1980), the critical core mass can vary from 1.5 M_{⊕} at grainfree 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 selfconsistently 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 lowmass 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, ISMlike opacity in the outer layers (κ_{gr} ~ 1 cm^{2}/g) that falls to very low values (κ_{gr} ~ 10^{3} cm^{2}/g) in the deep layers close to the radiativeconvective 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 f_{opa}. 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 f_{opa} were often considered, typically from 10^{3} to 1. The same basic approach was also taken in Paper I, with the improvement that f_{opa} was calibrated by determining which reduction factor leads to the same formation timescales for Jupiter as in MBPL10. A bestfitting value of f_{opa} ≈ 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 semimajor 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 dusttogas 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 closein planets (Jin et al. 2014). This will help to better understand the role of the grain opacity in shaping the observed massradius 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 v_{set} 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/v_{set}. The scale height in the atmosphere is (1)where k_{B} is the Boltzmann constant, T the temperature of the gas, μ its mean molecular weight (≈2.4 for solar composition and molecular hydrogen), m_{H} the mass of a hydrogen atom, and g is the gravitational acceleration that is given as g = GM(R) /R^{2}. 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 d_{m} equal to that of molecular hydrogen H_{2} at standard conditions ≈2.7 × 10^{8} cm (Waldmann 1958). It is clear that, in reality, d_{m} 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 hardsphere molecules with a Maxwellian velocity distribution as ρv_{th}ℓ/ 3, where v_{th} is the mean thermal velocity of the gas, (3)More accurate expressions for η based on, e.g., the LennardJones 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 time^{3}. 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} = m_{gr}v/F_{D}, where F_{D} is the drag force and m_{gr} 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 m_{gr} = 4/3πρ_{gr}a^{3}, 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/cm^{3}.
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 dusttogas 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 grainfree 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 dusttogas mass ratio in the protoplanetary disk, f_{D/G,disk}. Typically f_{D/G,disk} ≈ 0.01, but this ratio could also be much lower if most solids have already been incorporated into large bodies. The dusttogas 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 multistaged 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).
Fig. 1 Mass deposition profile (ablated mass per unit length) as a function of height for 100 km planetesimals (solid blue line). The brown dashdotted line is the mass deposition for 1 km planetesimals, multiplied by a factor 10^{6}, while the green dashed line shows the result for 10 m planetesimals, multiplied by a factor 10^{12}, 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 10^{21} g). The planetesimals hit headon 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 SL9 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 grainfree 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 radiativeconvective 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 headon as assumed here. The consequences of different impact geometries are studied in Appendix A. As expected, offcenter 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 radiativeconvective 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., dṀ_{gr}/ dR = 0. Mass conservation then directly yields the dusttogas ratio f_{D/G} (by mass) in an atmospheric layer as (5)where v_{set} is the vertical settling velocity of the grains. Because both ρ and v_{set} vary with R, f_{D/G} will, in general, also be a function of the radius inside the atmosphere and will differ from the value at the top f_{D/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 n_{gr} = f_{D/G}ρ/m_{gr}. 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 dusttogas ratio to the typical particle size enters the opacity.
2.3. Drag regimes
We now give the equations for the dusttogas 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 selfconsistent (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 selfconsistent 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 dusttogas 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 dusttogas 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 coagulation^{4} can be derived (e.g., Helled & Bodenheimer 2011) by writing the coagulation timescale τ_{coag} as τ_{coag} = ℓ_{gr}/V_{th,gr}, where ℓ_{gr} is the mean free path of a grain between collisions with another grain, while V_{th,gr} is its mean thermal velocity. The mean free path is calculated as (17)where σ_{gr} = 4πa^{2} 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 dusttogas ratio as 1 /f_{D/G} and increases with grain size as a^{5/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, a_{mono}. 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 a_{mono} 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 a_{C,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 dusttogas ratio (22)The opacity is finally found by inserting a_{C,E} and f_{D/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 dusttogas 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 v_{th}, n_{gr}, f_{D/G}, and the expression for the settling velocity in the Epstein regime, this corresponds to(31)One thus finds that τ_{coal,E} ∝ 1 /f_{D/G} ∝ a/Ṁ_{gr}. Equating this with the settling timescale in the Epstein regime τ_{coal,E} = τ_{set,E} gives a typical grain size of (32)The dusttogas 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 κ_{gr} ∝ f_{D/G}/a, and that in this regime both f_{D/G} and a scale in the same way with . Physically it means that if we increase Ṁ_{gr}, we increase the dusttogas ratio which would in principle increase the opacity, but the increased dusttogas 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/(4Hρ).
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 f_{D/G,disk}, because of Eq. (4)) in Sect. 3.1.6.
With the definition of H and g, and the ideal gas law^{5}, 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 n_{gr} and f_{D/G}, this can be written as (37)which has the same dependency on f_{D/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 dusttogas 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, f_{D/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 v_{gas} that must be taken into account to describe the net velocity of the grains relative to the Eulerian atmospheric pT 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/v_{gas}. We have called this timescale the advection timescale for the following reason: In the limit that v_{gas} 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 F_{G} and drag force F_{D} balance each other, we again have F_{G} + F_{D} = 0, but F_{D} now has the form C_{couple}(v − v_{gas}), where the coupling constant C_{couple} depends on the drag regime (Eqs. (9) and (13)). Solving for the velocity of the grains, we have v = (F_{G} + C_{couple}v_{gas}) /C_{couple}. 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 = v_{gas}, 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 v_{gas} (but we show in Appendix C how the relative velocity can be included selfconsistently. 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 v_{gas}, 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 dusttogas ratio is simply f_{D/G,A} = Ṁ_{gr}/Ṁ_{XY} which is, in the case of deep deposition (Eq. (4)), simply equal to f_{D/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 dusttogas 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 v_{gas}, 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 v_{gas}. 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 R_{0} down to the current radius R. At R_{0}, the grains have an initial size a_{0} (this will typically be the same as the assumed monomer size a_{mono}). Inserting the equation for the growth timescale and v_{gas} 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 v_{gas}.
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 = a_{mono} (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 remain^{6}. 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 disklimited gas accretion rate is reached. Including or completely neglecting advection is found to lead to variations of this time of about 10^{5} 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 KnudsenLangmuir 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πa^{2} 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 p_{vap} is approximately given as p_{vap} = p_{vap,0}e^{− Cvap/T}, where p_{vap,0} and C_{vap} are material properties. We take the values suggested by Podolak et al. (1988) for rocky material, so that p_{vap,0} = 1.5 × 10^{13} dyn/cm^{2} and C_{vap} = 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 H_{2} 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 p_{vap} 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 welldefined evaporation temperature and evaporation radius R_{evap}.
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 grainfree 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 R_{evap} (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 grainfree 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.
Fig. 2 Extinction efficiency Q as a function of x_{max} = 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 dusttogas 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 πa^{2}.
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 n_{r} and imaginary n_{i} 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, n_{r}, and n_{i}, derived by Dr. J. Cuzzi^{7}.
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 C_{Wien}/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 n_{r} and n_{i}. 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(x_{max}) for grains with a = 1, 10, and 100 μm. The size parameter was calculated as x_{max} = 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 x^{4}) because of the complex refractory index (e.g., Hansen & Travis 1974). Expressed with n_{r} and n_{i}, the absorption efficiency can be approximated as (MP08)(55)For n_{r} = 1.5 and n_{i} = 0.1, the typical values mentioned by MBPL10, one finds Q_{abs} = 0.20x, while using n_{r} = 1.6 and n_{i} = 0.16 (which we find to provide a somewhat better fit to the tholine data of Khare et al. 1984) yields Q_{abs} = 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 wavelengthdependent 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 x_{max}, 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 wavelengthindependent 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 selfconsistent 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 grainfree 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 solarcomposition 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 wavelengthdependent 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 dusttogas ratio in the accreted gas f_{D/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.
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 dusttogas 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.
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 insitu 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^{6}M_{⊕}/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 ~10^{11} 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 selfconsistent in the sense that the temperature and density structures were obtained by Hubickyj et al. (2005) with a prespecified (and different) opacity, and not the one obtained in the grain calculations. In contrast, the calculations of MBPL10 discussed in Sect. 4 are selfconsistent. 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 yaxis) and the temperature (brown dashed line, right yaxis) 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/cm^{3} at the outer radius (which is comparable to typical nebular densities) to about 10^{6} g/cm^{3} 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 cm^{2}/g, and therefore comparable to ISM opacities. Towards the interior, it decreases strongly, going down to about 10^{3} cm^{2}/g at the inner boundary. It is the consequence of the increase of the typical grain size a and of the decrease of f_{D/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 10^{11} 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.
Fig. 4 MP08 comparison case. Left panel: timescale of growth by differential settling/coalescence (blue solid line) and Brownian coagulation (green dashdotted). 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 dashdotted 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 grainfree 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 × 10^{10} 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} cm^{2}/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 selfconsistently (see Figs. 11 and 13). In the simulation here, at an intermediate height of about 2.5 × 10^{11} 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 radiativeconvective boundary. This is indeed a very short timescale.
Fig. 5 MP08 comparison case. Left panel: grain radius (solid) and dusttogas 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 dusttogas 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, dusttogas ratio, and aerodynamic regime
Now that differential settling has been identified as the process determining the typical grain size, local dusttogas 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 10^{11} cm. Within this height, the size remains approximately constant. This change at 10^{11} 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.
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.
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 dashdotted line is the settling velocity found by MP08. As they only give the relative velocity v_{set}(R) /v_{set}(R_{out}), we have normalized their curve to the value predicted by the analytical model at R_{out}. Right panel: dust mass density (f_{D/G} × ρ) found by the analytical model (blue solid line) and by MP08 (thin dashdotted line). The two brown lines are the extinction efficiency, calculated with the fit of Eq. (54) (dashed) and directly with Mie theory (dashdotted line). 

Open with DEXTER 
The left panel of Fig. 5 also shows the local dusttogas ratio f_{D/G} as a function of height. It scales as 1/(R^{2}ρv_{set}) (Eq. (5)). The plot shows that f_{D/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 f_{D/G} and the increase of the grain size a both conspire to the strong decrease of κ_{gr} ∝ f_{D/G}/a in the deeper parts of the atmosphere.
When calculating Ṁ_{gr}, it was assumed that f_{D/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 noncontinuous 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 f_{D/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 10^{11} 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 10^{11} 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 v_{set}(R) /v_{set}(R_{out}) relative to the value in the top layer at R_{out} 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 R_{out}. 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 10^{11} cm for the large grains, and the reincrease 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 SiO_{2} 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.
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 dashdotted 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 masssettling 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 (ρ × f_{D/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 10^{11} cm are seen in both models. Interestingly, the bump only appears clearly in the product of ρ and f_{D/G}, but not in the two quantities alone.
There are also differences: in MP08, the dust density increases outside of about 4.5 × 10^{11} cm to reach at the outer boundary the value that is simply given by the nebular density times the dusttogas ratio in the disk, corresponding to about 10^{12} g/cm^{2}. 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 × 10^{11} cm. A similar length scale is also seen in the numerical model when the dusttogas 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 dashdotted 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 wavelengthindependent 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 grainfree gas opacity for a solar composition gas from Freedman et al. (2008). The right panel shows the corresponding optical depths of the atmosphere, τ_{atmo}.
Fig. 8 MP08 comparison case. Left panel: opacity as a function of height for eight different dusttogas ratios in the accreted gas f_{D/G,disk} indicated in the plot. Right panel: corresponding size of the grains (solid lines) and local dusttogas ratio (dashed lines) for the four dusttogas 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/cm^{2}), 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 × 10^{4}, 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 grainfree 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 × 10^{11} 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 grainfree 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 grainfree 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 grainfree 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 dusttogas ratio in the disk
The impact of different dusttogas ratios in the disk f_{D/G,disk} and therefore different Ṁ_{gr} = f_{D/G}Ṁ_{XY} 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 f_{D/G,disk} both for the Epstein and Stokes regimes (Eqs. (35), (41)).
Movshovitz & Podolak (2008) investigated the impact of varying f_{D/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 × 10^{11} cm), the opacity was nearly independent of f_{D/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× 10^{11} cm, MP08 found that f_{D/G,disk} is slightly anticorrelated 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× 10^{11} cm, the opacity increases with f_{D/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 f_{D/G,disk} between 0.0001 and 0.5, a very wide (and not necessarily realistic) range. One sees that for all f_{D/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(x_{max}) that approaches 2 (Fig. 2) as the grains become larger (see next). Only for f_{D/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 f_{D/G}. The analytical model thus agrees with the numerical result that the opacity is nearly independent of f_{D/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 f_{D/G,disk}. The figure shows the corresponding grain size a and the local dusttogasratio f_{D/G} as a function of height. For an approximately constant Q (because of the large grains relative to the wavelength), κ_{gr} scales simply ∝ f_{D/G}/a (Eq. (8)). One sees that high f_{D/G,disk} indeed lead to a higher f_{D/G} which would increase the opacity, but at the same time these higher grain concentrations also lead to larger grains, reducing the opacity. As f_{D/G} and a depend in the settling regime in the same way on f_{D/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 f_{D/G}, bouncing and fragmentation could become important (Sect. 4.3.4). This could limit the efficiency of grain growth in very metalrich envelopes.
Up to now, we have discussed the effect of varying the grain input in the top layer that is caused by a different f_{D/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 f_{D/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 f_{D/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 firstorder 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, dragenhanced 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 recalculating the opacity for the case of shallow deposition (Eq. (4)), adding grains in the top layer at a rate of 10^{6}M_{⊕}/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 dusttogas ratio in the dominant differential settling regime. However, the quantity that enters the equations is not f_{D/G,disk} alone, but f_{D/G,disk}Ṁ_{XY}. This means that for a fixed dusttogas 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
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 a_{mono}. 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 a_{mono} = 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 cm^{2}/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 a_{mono} = 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 a_{mono} 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.
Fig. 10 Core (solid line) and envelope mass (dashdotted 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/cm^{2} (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/cm^{3}, 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 dusttogas ratio in the disk f_{D/G,disk} = 0.01 for all planetesimal surface densities as in MBPL10. We thus ignore potential correlations of f_{D/G,disk} and the planetesimal surface density. We have also run nonnominal simulations that use a f_{D/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 a_{mono} = 1.26μm.
The third module solves the internal structure equations for the gaseous envelope in the quasistatic 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/cm^{2}. 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/cm^{2}, 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 socalled 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 10^{4} 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^{2}M_{⊕}/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.
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 10^{4} − 10^{5} 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 superEarthtype planets. Such planets could be the progenitors of the recently observed lowmass, lowdensity 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 lowmass 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.
Impact of f_{D/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^{2}M_{⊕}/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 10^{5} 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.
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 dusttogas 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 dusttogas ratio in the disk. To quantify the impact of f_{D/G,disk} on the overall formation history, we have recalculated the Σ10 and Σ4 cases also for other f_{D/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 f_{D/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 f_{D/G,disk} = 0.001, 0.01, and 0.1. This means that a lower dusttogas 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 f_{D/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 f_{D/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 f_{D/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 f_{D/G,disk} = 10^{4} (≈25% difference). This is due to the following: in the Σ4 case, at low f_{D/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 cm^{2}/g in the outer parts of the atmosphere (the specific value might, however, be affected by our simplistic way of calculating Q). For higher f_{D/G,disk}, Q is always between 2 and 3. In the Σ10 case, the grains also become smaller at lower f_{D/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 f_{D/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 f_{D/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 a_{mono}. 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 a_{mono} = 10 and 100 μm instead of the nominal 1 μm has similar consequences: for a_{mono} = 10μm, the formation is virtually the same as in the nominal case. For a_{mono} = 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 selfconsistent.
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 coreenvelope 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^{6}M_{⊕} 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 coreenvelope 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 × 10^{10} 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×10^{11} 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 cm^{2}/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 × 10^{10} cm and 3.6 × 10^{11} 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 × 10^{11} 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} cm^{2}/g at R = 7.5 × 10^{10} 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 × 10^{10} cm is thus given by the Freedman et al. (2008) data. This leads in particular to a second local minimum in κ at about 4 × 10^{10} 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 × 10^{10} 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 × 10^{10} 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.
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 cm^{2}/g at the top, and then decreases to ~10^{3} cm^{2}/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 lowmass 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 lowmass planets are also found to have a relatively large outer part where the opacity is on the order of 1 cm^{2}/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 × 10^{10} cm. A second tiny inner radiative zone occurs in the range 4.0–4.6 × 10^{10} cm. It is clearly related to the local minimum in the molecular opacities in this region. Within 4.0 × 10^{10} 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 radiativeconvective boundary is at about 3.7 × 10^{10} cm, corresponding to our inner radiativeconvective boundary at 4.0 × 10^{10} 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 lowmass lowdensity planets which were recently discovered in large numbers (e.g., in the Kepler11 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} cm^{2}/g and 10^{3} cm^{2}/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 coreenvelope 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 dusttogas 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 SiO_{2} 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).
Fig. 14 Overlay of the settling velocitymass 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 velocitygrain 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 lowerleft 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 (p →P)”. However, growth can still occur if the collisions occur mainly between bodies of a different size that are both porous (represented as “p →P”) and/or if small porous objects hit large compact ones (represented as “p →C”). 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 (p →P)”. In this regime, the collisional outcome transitions from growth to bouncing also in the p →P 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 p →P case), fragmentation occurs, except in the p →C 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 lowmass 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 massradius 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/ 8Hρ in the Epstein drag regime and 2Q/Hρ 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 ISMlike values (~1 cm^{2}/g), but, as we move deeper into the atmosphere, it decreases sharply to very low values (~10^{3} cm^{2}/g). At this point, the grain opacity becomes so low that the grainfree 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 (M_{core} = 12.6 M_{⊕}, M_{env} = 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 × 10^{4}. For comparison, in a grainfree atmosphere, an optical depth of about 25 is found. This indicates that grainfree 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 grainfree 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 lowmass 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 dusttogas 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 dusttogas 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 dusttogas 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 massradius 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 gasdominated planets (e.g., Marcy et al. 2014).
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.
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.
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.
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 MaxPlanckGesellschaft for the ReimarLüst Fellowship. I also thank the referee Dr. Morris Podolak for an interesting and constructive report that helped to improve the manuscript.
References
 Alibert, Y., Mordasini, C., & Benz, W. 2004, A&A, 417, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ayliffe, B. A., & Bate, M. R. 2012, MNRAS, 427, 2597 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Münch, M. 1993, Icarus, 106, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bodenheimer, P. H., & Pollack, J. B. 1986, Icarus, 67, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles (New York: Wiley) [Google Scholar]
 Boss, A. P., Alexander, C. M. O., & Podolak, M. 2012, Earth Planet. Sci. Lett., 345, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bronsthen, V. A. 1983, Physics of Meteoric Phenomena (Dordrecht: D. Reidel Publishing Company) [Google Scholar]
 Chyba, C. F., Thomas, P. J., & Zahnle, K. J. 1993, Nature, 361, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Cooper, C. S., Sudarsky, D., Milsom, J. A., Lunine, J. I., & Burrows, A. 2003, ApJ, 586, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Cuzzi, J. N., Estrada, P. R., & Davis, S. S. 2014, ApJS, 210, 21 [NASA ADS] [CrossRef] [Google Scholar]
 D’Angelo, G., Weidenschilling, S. J., Lissauer, J. J., & Bodenheimer, P. 2014, Icarus, 241, 298 [NASA ADS] [CrossRef] [Google Scholar]
 Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503 [NASA ADS] [Google Scholar]
 Fischer, D. A., & Valenti, J. A. 2005, ApJ, 622, 1102 [NASA ADS] [CrossRef] [Google Scholar]
 Fortier, A., Alibert, Y., Carron, F., Benz, W., & Dittkrist, K.M. 2013, A&A, 549, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Gaidos, E., Fischer, D. A., Mann, A. W., & Lépine, S. 2012, ApJ, 746, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, G. 1997, MNRAS, 285, 403 [NASA ADS] [CrossRef] [Google Scholar]
 Greenzweig, Y., & Lissauer, J. L. 1992, Icarus, 100, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hansen, J. E., & Travis, L. D. 1974, Space Sci. Rev., 16, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Helled, R., & Bodenheimer, P. 2011, Icarus, 211, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Hori, Y., & Ikoma, M. 2010, ApJ, 714, 1343 [NASA ADS] [CrossRef] [Google Scholar]
 Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419 [NASA ADS] [CrossRef] [Google Scholar]
 Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Hubickyj, O., Bodenheimer, P. H., & Lissauer, J. J. 2005, Icarus, 179, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Ida, S., & Lin, D. N. C. 2004, ApJ, 616, 567 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 Inaba, S., & Ikoma, M. 2003, A&A, 410, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Iseli, M., Küppers, M., Benz, W., & Bochsler, P. 2002, Icarus, 155, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Khare, B. N., Sagan, C., Arakawa, E. T., et al. 1984, Icarus, 60, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Kornet, K., Bodenheimer, P. H., Rozyczka, M., & Stepinski, T. F. 2005, A&A, 430, 1133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lenzuni, P., Gail, H.P., & Henning, T. 1995, ApJ, 447, 848 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Thommes, E. W., & Duncan, M. J. 2010, ApJ, 139, 1297 [Google Scholar]
 Lissauer, J. J. 1993, ARA&A, 31, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Hubickyj, O., D’Angelo, G., & Bodenheimer, P. H. 2009, Icarus, 199, 338 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., JontofHutter, D., Rowe, J. F., et al. 2013, ApJ, 770, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Mac Low, M.M., & Zahnle, K. J. 1994, ApJ, 434, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, N. K., & Fortney, J. J. 2011, ApJ, 736, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, H. 1980, Prog. Theor. Phys., 64, 544 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Mizuno, H., Nakazawa, K., & Hayashi, C. 1978, Prog. Theor. Phys., 60, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Alibert, Y., & Benz, W. 2006, in Tenth Anniversary of 51 Pegb: Status of and prospects for hot Jupiter studies, eds. L. Arnold, F. Bouchy, & C. Moutou (Paris: Frontier Group), 84 [Google Scholar]
 Mordasini, C., Klahr, H., Alibert, Y., Benz, W., & Dittkrist, K.M. 2010 [arXiv:1012.5281] [Google Scholar]
 Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012a, A&A, 541, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Alibert, Y., Georgy, C., et al. 2012b, A&A, 547, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012c, A&A, 547, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Klahr, H., Alibert, Y., Miller, N., & Henning, T. 2014, A&A, 566, A141 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Movshovitz, N., & Podolak, M. 2008, Icarus, 194, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Movshovitz, N., Bodenheimer, P. H., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616 [NASA ADS] [CrossRef] [Google Scholar]
 Nayakshin, S. 2010, MNRAS, 408, 2381 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Opik, E. J. 1958, Physics of meteor flight in the atmosphere (New York: Interscience Publishers) [Google Scholar]
 Ormel, C. W. 2013, MNRAS, 428, 3526 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W. 2014, ApJ, 789, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Papaloizou, J. C. B., & Nelson, R. P. 2005, A&A, 433, 247 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perri, F., & Cameron, A. G. W. 1974, Icarus, 22, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Podolak, M. 2003, Icarus, 165, 428 [NASA ADS] [CrossRef] [Google Scholar]
 Podolak, M. 2004, in Rev. Mex. Conf. Ser., 22, eds. G. GarciaSegura, G. TenorioTagle, J. Franco, & H. W. Yorke, 104 [Google Scholar]
 Podolak, M., Pollack, J. B., & Reynolds, R. T. 1988, Icarus, 73, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J. B., Podolak, M., Bodenheimer, P., & Christofferson, B. 1986, Icarus, 67, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P. H., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R. 2006, ApJ, 648, 666 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, L. A., Bodenheimer, P. H., Lissauer, J. J., & Seager, S. 2011, ApJ, 738, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Rossow, W. B. 1978, Icarus, 36, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Santos, N. C., Israelian, G., & Mayor, M. 2001, A&A, 373, 1019 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Stepinski, T. F., & Valageas, P. 1996, A&A, 309, 301 [NASA ADS] [Google Scholar]
 Stevenson, D. J. 1982, Planet. Space Sci., 30, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Svetsov, V. V., Nemtchinov, E. V., & Teterev, A. V. 1995, Icarus, 116, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Tanigawa, T., & Ohtsuki, K. 2010, Icarus, 205, 658 [NASA ADS] [CrossRef] [Google Scholar]
 Valencia, D., Guillot, T., Parmentier, V., & Freedman, R. S. 2013, ApJ, 775, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283 [NASA ADS] [CrossRef] [Google Scholar]
 Waldmann, L. 1958, Transporterscheinungen in Gasen von mittlerem Druck. In Handbuch der Physik, Band XII, Thermodynamik der Gase, ed. S. Flügge (Berlin, Göttingen, Heidelberg: SpringerVerlag) [Google Scholar]
 Weidenschilling, S. 1977, Astrophys. Space Sci., 51, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1984, Icarus, 60, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wolfgang, A., & Laughlin, G. 2012, ApJ, 750, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Zahnle, K. J. 1992, JGR, 97, 10243 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Online material
Appendix A: Planetesimal mass deposition: impact geometry
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: zoomin 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 headon impact geometry. However, it is clear that offcenter impacts are much more likely than the (near) headon case because of the larger collisional cross section. In this Appendix we investigate the consequences of noncentral 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 p_{crit}. 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 twobody 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 p_{mean} = p_{crit}/.
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 zoomin 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 10^{16} 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 ~10^{31} erg/s is reached.
Figure A.2 shows the corresponding radial mass deposition profiles. The curve for the headon 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 headon 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ρv^{2} where ρ is the atmospheric gas density). However, the headon 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 (alongorbit) 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 R_{evap} 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 f_{D/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.
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 headon 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 headon 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 headon icy impactors, the same as in Fig. 1, are also shown for comparison.
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 headon 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/cm^{2} for ice and silicate, respectively), the tensile strength (4×10^{6} and 3.5 ×10^{8} dyn/cm^{2}, 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 headon 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 (R_{evap} ≈ 5 R_{♃}), and pure thermal ablation is inefficient for such large bodies. This is due to the low heat transfer coefficient and the low surfacetomass 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 surfacetomass 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 intermediatesized 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 “corehit 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 headon 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 headon 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 (R_{evap}). 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 gasgrain 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 v − v_{gas} and equating with the gravitational force then yields terminal grain velocity v_{trav}. It is given as (C.1)From this we get the timescale on which grains traverse one scale height, τ_{trav} = H/v_{trav}. 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 v_{gas} = 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/v_{gas}.
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 v_{trav} and assuming again a radially constant grain flux yields the growth timescale (C.3)For v_{gas} = 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 v_{gas} = 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 nonzero v_{gas} = 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 × 10^{10} 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 nonzero v_{gas} 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^{4}M_{⊕}/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 × 10^{11} 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 dashdotted 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.).
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/cm^{2}). 

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
Characteristics of the Jupiter in situ formation simulations at the crossover point.
All Figures
Fig. 1 Mass deposition profile (ablated mass per unit length) as a function of height for 100 km planetesimals (solid blue line). The brown dashdotted line is the mass deposition for 1 km planetesimals, multiplied by a factor 10^{6}, while the green dashed line shows the result for 10 m planetesimals, multiplied by a factor 10^{12}, 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 
Fig. 2 Extinction efficiency Q as a function of x_{max} = 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 
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 
Fig. 4 MP08 comparison case. Left panel: timescale of growth by differential settling/coalescence (blue solid line) and Brownian coagulation (green dashdotted). 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 dashdotted 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 
Fig. 5 MP08 comparison case. Left panel: grain radius (solid) and dusttogas 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 
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 dashdotted line is the settling velocity found by MP08. As they only give the relative velocity v_{set}(R) /v_{set}(R_{out}), we have normalized their curve to the value predicted by the analytical model at R_{out}. Right panel: dust mass density (f_{D/G} × ρ) found by the analytical model (blue solid line) and by MP08 (thin dashdotted line). The two brown lines are the extinction efficiency, calculated with the fit of Eq. (54) (dashed) and directly with Mie theory (dashdotted line). 

Open with DEXTER  
In the text 
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 dashdotted 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 
Fig. 8 MP08 comparison case. Left panel: opacity as a function of height for eight different dusttogas ratios in the accreted gas f_{D/G,disk} indicated in the plot. Right panel: corresponding size of the grains (solid lines) and local dusttogas ratio (dashed lines) for the four dusttogas ratios in the accreted gas indicated in the plot. 

Open with DEXTER  
In the text 
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 
Fig. 10 Core (solid line) and envelope mass (dashdotted 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/cm^{2} (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 
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 
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 coreenvelope 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^{6}M_{⊕} 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 
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 
Fig. 14 Overlay of the settling velocitymass 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 velocitygrain 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 
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: zoomin 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 
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 headon 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 
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 headon 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 
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 
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/cm^{2}). 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.