© ESO 2018

1 Introduction

The growth of planetary cores is a conceptually simple process in the planetesimal-driven core accretion scenario. An increasingly massive planetary embryo is impacted by km-sized impactors that add to its mass (Pollack et al. 1996; Hubickyj et al. 2005). Soon, the small protoplanet starts to gravitationally bind an envelope of gas. This envelope is initially very poor in high-Z material due to the limited interaction between the envelope and any impacting planetesimals (Podolak et al. 1988). Core growth eventually comes to a halt when the planet has no more solid material in its feeding zone left to accrete. Unfortunately, the planetesimal-driven core accretion scenario faces a time-constraint problem at distances beyond several AU. Simulations indicate that planet formation typically requires more time than the expected lifetimes of the disks from which their matter originates (Kobayashi et al. 2010; Levison et al. 2010; Bitsch et al. 2015). The underlying problem is the strong dynamical excitation a growing planet induces to the planetesimal disk, i.e., a negative feedback effect. Although planetesimal-driven planet formation scenarios may still be viable with optimal choices regarding planetesimal sizes and strengths (e.g., Kobayashi et al. 2016), an alternative approach is to investigate formation models that are not affected by this problem.

One such scenario is pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012). Instead of km-sized planetesimals, the protoplanet accretes smaller objects (cm- to m-sized) called pebbles. There are two main differences that originate from this size difference. Firstly, pebble accretion is expected to accrete mass approximately an order of magnitude more rapidly than planetesimals (Lambrechts et al. 2014), and operates particularly well when the pebbles are settled into the mid-plane regions (Ormel 2017). The reason for this is the more efficient capturing of pebbles due to their increased susceptibility to gas drag. Secondly, pebbles feature a stronger interaction with a planet’s atmosphere. Their higher area-to-mass ratio makes it easier for them to be slowed down by aerodynamic gas drag and reduces their ablation timescales (e.g., Love & Brownlee 1991; McAuliffe & Christou 2006). As a result, these pebbles can deposit high-Z vapor into even a small growing planet. This in turn can have a large effect on the planet’s further growth and evolution (Venturini et al. 2016). The presence of high-Z material locally increases the mean molecular weight, which drives up densities and temperatures, increasing the envelope mass.

So far, research on pebble accretion has mainly focused on the first difference by calculating accretion rates though disk evolution and drag considerations (e.g., Morbidelli & Nesvorny 2012; Chambers 2014; Ida & Guillot 2016). The evaporation effects have only recently been considered by Alibert (2017), based on impact simulations by Benz et al. (2006). In his simulations, direct core impacts terminate either when the planet becomes too hot (T_core > 1600 K) or the envelope becomes too massive (M_env > 10⁻⁴ M⊕) for pebbles to reach the core. He found these conditions to be met before the core grew to 1 M⊕ . However, heneglected any enrichment effects and did not compute the planet’s evolution beyond the point of full ablation. In reality, the core can continue to grow if ablated material oversaturates the gaseous atmosphere, and the leftover material rains out to the core (Iaroslavitz & Podolak 2007).

In this work, we investigate how cores grow by pebble accretion. We expand previous calculations of core growth by including enrichment effects from when they first occur and by considering further growth through the rainout of ablated material. Our code simulates planet formation along with pebble impacts in a way that quantitatively incorporates how the pebbles interact with the growing atmosphere. We calculate the effects of ablation and gas drag on impacting pebbles and relate these effects to the planet’s subsequent evolution by modifying the luminosity profile and mass deposition curve. Our results show that impacting pebbles will be fully ablated before the planet even reaches 0.5 M⊕ . Subsequent rainout of ablated material in a super-saturated atmosphere can add ~0.1 M⊕ to the core mass, but this direct core growth terminates well before the core mass has reached 1 M⊕.

This paper is organized as follows. Sections 2.1 and 2.2 introduce our impact and planet simulation models. Section 2.3 elaborates on our treatment of high-Z enrichment and on the rainout of ablated material to the core. Our pebble simulation results for SiO2 and H2 O are presented in Sects. 3.1–3.3. We compare planetesimal mass deposition in Sect. 4.1 and examine the core-surface conditions in Sect. 4.2. Section 5 covers the discussion and conclusions.

2 Model description

The goal of this work is to provide an estimate of the core sizes that planets can form directly by pebble accretion. To clarify our approach, we distinguish three distinct mechanisms by which planetary core growth can occur. Direct impact of solids reaching the core, the rainout of ablated material if the envelope is super-saturated, and a phase change of gas in the inner envelope under sufficiently high pressure. We consider only the first two mechanisms in this work, and show that the third is not yet relevant for our early model. This is what we mean by direct core growth. We limit ourselves to the solids accretion phase when the planet’s accretion is still dominated by solids instead of disk gas and Kelvin–Helmholtz contraction is not yet significant.

In our formation scenario, a growing protoplanet is constantly being impacted by small, cm- to m-sized objects. These impactors partially or fully evaporate during entry and contribute high-Z material (in our case SiO2, quartz) to the core and envelope. To simulate the planet’s growth, we used a new evolution code consisting of two components. The first is a time-dependent calculation of the planet’s structure and gas accretion rate. We modeled its evolution, as is typically done, by a series of quasi-hydrostatic models with a separate core and envelope (e.g., Mordasini et al. 2012; Piso & Youdin 2014). The second component of our code is a calculation of the interaction between impacting pebbles and the planet’s atmosphere. The accretion of solids is modeled by a single-body simulation code that uses the layers generated by the structure code. Impactors hit the planet head-on and lose mass by thermal ablation and friction. The resultant core accretion and mass deposition rates affect the planet’s further evolution. In principle, the two codes should be iterated to make sure they work self-consistently. In tests, one iteration was found to be sufficient to achieve this. Therefore, we ran them twice during every time step.

2.1 Impact model

2.1.1 Gas drag

The pebbles are modeled as perfect spheres of uniform size and a constant density of 2.65 × 10³ kg m−3. We further assume that they impact the planet head-on. Pebbles that lack spherical symmetry or enter the planet’s atmosphere at an angle will experience additional ablation due to their increased travel time. In this sense, we study the limit in which pebbles can most easily reach the core. Their initial velocities are set to the lower of the terminal and escape velocities. The first is a drag-induced limitation, calculated by the velocity at which the local drag force is exactly equal in magnitude to the gravitational force. The second is the limitation of the total gravitational energy that an impactor can acquire as it moves towards the planet. As impactors approach the core, gravity accelerates them and gas drag slows them down. The drag force is given by (1) where ρ_g is the gas density, A_i the impactor’s frontal area, v_i its relative velocity to the gas, and Cd the drag constant. Most of the physics is contained in this drag constant. For our purposes, we want to be able to describe drag over a wide range of conditions for impactors varying in size from 1 cm to several kilometers. There has been extensive research on empirical expressions for different drag regimes. Weadopt an extended version of the general expression by Melosh & Goldin (2008) based on a review of this literature. Theirs is a continuous expression of Cd that is applicable to a wide range of conditions. It depends on the local mach number Ma = v_i∕c_{s, g}, where c_{s, g} is the local sound velocity; on the Reynolds number Re =v_iR_i∕η_g, where η_g is the dynamical viscosity of the gas and R_i is the impactor radius; and on Knudsen’s number (Kn), which is defined as the ratio between the Mach and Reynolds numbers, Kn = Ma∕Re. The full expression for C_d reads as (2)

where C₁, C₂ , and H are the auxiliary functions

and T_s and T_g are the impactor’s surface temperature and the local gas temperature, respectively; γ_g is the adiabatic index of the gas, explained in detail in Sect. 2.2; and C₁ is an extended version of the expression given by Melosh & Goldin (2008) as suggested by D’Angelo & Podolak (2015), based on an empirical review from Brown & Lawler (2003). Their expression only breaks down at Reynolds numbers exceeding approximately 3 × 10⁵ when the type of flow around the sphere becomes turbulent. We have modeled this turbulent part of the flow as a step-function decrease of Cd to 0.2 (e.g., Michaelides 2006). To illustrate the importance of including a varying drag constant, Fig. 1 shows the drag constant evolution of a 0.1 m SiO2 pebble impacting a 0.3 M⊕ planet at escape velocity. During this impact, Kn is near unity at the edge of the planet and drops close to zero as the pebble approaches the core. In this regime, Kn < 1 and Ma > 1 for most of the impact, resulting in values of C_d that are considerably different from the low-v drag expressions typically used in the planet formation literature (Weidenschilling 1977; Whipple 1972).

Fig. 1 Drag constant evolution of a 0.1 m SiO2 pebble impacting a 0.3 M⊕ planet at escape velocity. The blue line represents the drag constant, the green line indicates the impactor’s velocity relative to the surrounding gas, and the black dashed curve represents the core radius.

2.1.2 Ablation

When an impactor enters the planet’s envelope, it becomes surrounded by increasingly hot gas. We include two heating mechanisms, radiation and friction. Thermal radiation is approximated by perfect Planck curves. The net total of irradiated energy is then given by the Stefan–Boltzmann expression with the area of a sphere, (6)

where α_i is the pebble’s average absorption constant, assumed to be unity. We use ambient gas temperatures even when incoming objects are supersonic. Gas temperatures and densities arehighest near the core. We find that thermal ablation in combination with impactors being slowed down typically leads to very localized ablation in this inner region.

For large or fast impactors, frictional heating becomes important. This is generally the dominant heating mechanism in the cold outer region of a planet’s atmosphere where thermal ablation is minimal but the impactors travel at high velocities. These objects lose significant amounts of kinetic energy, a fraction of which (f_h ) gets converted into heat: (7)

Estimating the numerical value of f_h is a complex problem. Most works on impacts assume it to be equal to some constant of varying size (e.g., Podolak et al. 1988; Pinhas et al. 2016). We allow f_h to vary during the impacts and evaluate f_h based on gas characteristics with the expressions from Melosh & Goldin (2008), (8)

where r^′ is the recovery factor, evaluated as (Mills 1999). In addition to Re, f_h also depends on the Nusselt (Nu) and Prandtl (Pr) numbers where C_{p, g} is the specific heat of the gas, Nu_c and M^′ are auxiliary functions

and η_g is the dynamic viscosity of the gas, the components k of which we evaluate with the common expression from kinematic theory (e.g., Dean 1985) (13) where d_k and m_k are the molecular diameter and weight of the gas constituent and k_b is the Boltzmann constant. The values of η_k are summed by volume fraction to yield η_g . In the simulated impacts of pebbles, the friction factor is found to vary between approximately 0.5 and 0.05.

As an impactor heats up, its inner region is initially shielded from the heat. This leads to the formation of a temperaturegradient. The importance of this gradient can be inferred from the Biot number (14)

where h^′ is the characteristic heat transfer coefficient; L is a characteristic length equal to R∕3 for a sphere; and k_i is the thermal conductivity, which we take equal to 2 W K m⁻¹. This is a commonly used value, based on data from Powell et al. (1966) (e.g., Adolfsson et al. 1996). In the context of planetary entry, h^′ can be approximated as (Love & Brownlee 1991) (15) where ΔT is the temperature difference between the local gas and the impactor’s surface. Whenever the Biot number is less than 0.1, an object can roughly be considered as isothermal. A typical first order approach to modeling a temperature gradient is setting the Biot number equal to 0.1 and interpreting the resultant radius as an isothermal layer of temperature T_s with depth and mass (e.g., Love & Brownlee 1991; McAuliffe & Christou 2006)

The total heat contents of the impactor are therefore approximated by , where the temperature of the interior region is assumed to be equal to that of the disk, T_disk. When the temperature of the surface layer has increased sufficiently, atoms start to evaporate into the envelope. This process can be characterized by the Langmuir formula (18) where M_i is the impactor mass, m_g is the mean molecular weight of the gas, and P_v is the vapor pressure of the impactor (see Sect. 2.3). Evaporation prevents the impactor from heating indefinitely. The surface temperature starts to decrease again when the energy loss through evaporation exceeds the incoming heat. Using the conservation of energy, the impactor’s surface temperature can be obtained from (19)

In summary, the terms in Eq. (18) denote contributions respectively from radiative cooling, frictional heating, evaporative cooling, and internal energy. This temperature balance is invoked to yield the equilibrium temperature of the isothermal layer T_s with a resolution of 10⁻³ K at any time step. Such a precise temperature determination is necessary to ensure numerical stability despite the strong temperature dependence of the vapor pressure.

2.2 Planetary structure

We integrate the planet’s atmosphere from the outside in, starting at either the Hill or the Bondi radius, depending on which is smaller. At this radius the envelope’s composition and conditions are taken to be the same as the disk’s. Our code solves for two boundary conditions at the core surface: the mass condition and the luminosity condition (20)

where χ_{i, core} and v_{i, core} are the surviving mass fraction and velocity of impactors that reach the core and Ṁ_acc is the solids accretion rate. We consider accretion to be the only energy source at the core and neglect any additional contributions, such as radioactive heating or core contraction. The total luminosity and envelope mass are adjusted such that these conditions are matched at every time step. Disk gas flows into the envelope to fit the increased interior mass.

The structure of the planet is divided into two parts: a solid core and a gaseous envelope. We use a typical constant density for the core of 3.2 × 10³ kg m−3 (e.g., Pollack et al. 1996; Hubickyj et al. 2005). The core grows only when solid material reaches its boundary. In our model, this can occur in two ways. Either by direct impacts of accreting material or by the rainout of ablated silicates when the vapor saturation of the envelope is exceeded. The envelope’s structure is calculated with the one-dimensional spherically symmetric stellar structure equations where G is the gravitational constant, P is pressure, T the temperature, M the total mass interior to radius r, ρ the envelope density, and r the distance to the planet’s center. We use the Schwarzschild criterion to determine the thermal gradient ∇_th as , where ∇_conv is approximated by the adiabatic gradient ∇_ad , which we take as a constant based on the local gas composition (24)

where γ_g is the average adiabatic index of the gas, defined as the ratio of the specific heat under constant pressure and volume C_{p, g}∕C_{v, g}. We use the fitting formula from Chase (1998) to evaluate the components C_{p, k} of C_{p, g} as a function of temperature

where A_k − E_k are constants corresponding to the component k (here k = (H₂ , He, SiO₂ , H₂ O)). To calculate the components C_{v, k} = C_{p, k}∕γ_k, we assume constant adiabatic indices γ_k equal to , , and to find (27)

where f_{g, k} are the local gas mass fractions. For the radiative outer zone of the envelope, we use (28)

where κ is the opacity and L the luminosity. Following Ormel (2014), the opacity is assumed to originate only from the gas and is approximated by the analytical expression from Bell & Lin (1994), (29)

valid for a solar-type composition. In the inner, highly enriched regions Eq. (29) would no longer apply, but these regions are adiabatic and thus κ does not enter the structure equations. We do not assume a grain opacity, as these are likely to be very small as a result of coagulation and sedimentation processes (Movshovitz et al. 2010; Mordasini 2014; Ormel 2014).

The luminosity has four components, kinetic energy deposition by the impactors, (negative) latent heat release upon ablation, high-Z rainout, and contraction, (30)

where K_i is the kinetic energy of an impactor, Ṁ _rain is the local rainout rate, and V is the specific gas volume 1∕ρ. The core luminosity is determined by the mass fraction (χ_{i, core}) of the pebbles that reaches the core and their kinetic energy (see Eq. (20)). The structure equations are supplemented by the ideal gas equation of state for a mixed composition: (31)

The range of conditions we experience in our simulations is still well approximated by the ideal gas equation (see Sect. 3.2). Temperatures inside the planet’s inner region can exceed those required for the dissociation of hydrogen (Lee et al. 2014), but these regions are modeled as silicate-dominated and contain little hydrogen. The effective difference with using a more sophisticated equation of state is therefore minor. An explicit test of this is included in Sect. 3.2.

The structure equations are integrated using a fourth-order Runge–Kutta scheme. We have set the time step Δt = ΔM∕Ṁ_acc such that the planet accretes 0.005 M⊕ at every interval. The grid consists of set logarithmic distances to provide enough resolution near the core and consists of approximately 20 000 layers per integration. The reason for this high number is to ensure a smooth temperature profile of the inner region of the planet. This reduces the number of iterations needed in the impact code and speeds up run-times.

2.3 High-Z enrichment and rainout

We consider the planet to be embedded in a gaseous disk with mass fractions of 75% molecular hydrogen and 25% helium. The initial envelope and the accreted gas do not contain any high-Z materials. However, as the protoplanet is impacted by pebbles, some of their material ablates into the atmosphere before it reaches the core. We consider the effects this has on the planet’s evolution and subsequent impacts by locally increasing the mean molecular weight where high-Z material is present in vapor form. The maximum fraction (by mass) of high-Z vapor that the gas can contain depends on the local silicate vapor and total pressure (32)

where is the local mass fraction of silicates in vapor phase. Equation (32) can be interpreted as the high-Z saturation curve of the gas, depending on pressure and temperature. Vapor pressures can be expressed by the Clausius–Clapeyron equation (33)

where a₀ = 29.5, a₁ = 46071.4 K, and a₂ = 58.9 K are constants determined by the fit published on NIST, based on original data from Stull (1947). We iterate on μ_g at every time step to ensure an accurate estimation. When the gas temperature increases such that , we limit to unity. In this case, the atmospheric layer consists entirely of silicate gas. The temperature dependency of the vapor pressure limits the presence of silicate vapor to the envelope’s warmer interior layers. Finding the partial pressure also requires identifying how much total silicate mass is present in a given layer. This is affected by the amount of material that has been ablated up to that point, but also depends on what happens to the silicates after ablation (i.e., outward mixing or settling). Tosimplify matters, we look at two contrasting scenarios:

• 1.

  In the mixing case, we assume that all the ablated material mixes uniformly through the entire envelope. This includes the material in solid and in vapor form. We track the total high-Z mass that has ablated up to that point in time, M_Z,abl. Meanwhile, the planet has also attracted hydrogen and helium gas masses and M_He from the disk. Complete mixing then leads to one global silicate mass fraction (vapor + solids) that is calculated as (34)

• 2.

  In the outer layers, and most of the silicate mass is in solid form. This trend continues inward up to the point that the temperature of a layer is sufficiently high that . Here, mixing limits the silicate vapor mass fraction to of Eq. (34). Direct core growth in this scenario only continues as long as some of the impactors can directly reach the core.

• In the rainout case, any material in excess of what can be contained in vapor (see Eq. (32)) falls to the layer below. The resulting rainout of high-Z material increases the planet’s luminosity as descending silicate mass loses gravitational energy on its way to the core. If the envelope cannot absorb all the ablated material, the rest of it adds to the core mass. This allows the core to grow beyond the point that impactors fail to reach the core directly. Direct core growth ceases when all the accreting solids can be fully absorbed by the envelope as vapor.

Vapor pressures rise quickly as the local gas temperature increases. Close to the core, conditions can be such that the vapor pressure exceeds the total pressure, and the high-Z mass fraction approaches one. The high mean molecular weight of this high-Zlayer surrounding the core leads to very steep pressure and density curves, as can be seen in Fig. 4.

3 Results

3.1 Standard model

The simulations shown here use the parameters from Table 1 unless specifically stated otherwise. We simulate the planet at a distance of 5.2 AU fromthe central star, embedded in a disk with a typical local temperature and density of 150 K and 5 × 10⁻⁸ kg m−3 (Hubickyj et al. 2005). We use a standard pebble accretion rate of 10⁻⁵ M_⊕ yr⁻¹ (Lambrechts et al. 2014). The goal is to estimate to what mass a core can grow before direct core accretion stops. This is either when impacting pebbles fully evaporate in our mixing scenario, or when the atmosphere can absorb all ablated material in the rainout case. Figure 2 shows the result of our SiO2 pebble simulations for both of these assumptions.

In general, we can distinguish three phases of core growth:

• 1.

  When the core mass is still below ~0.23–0.39 M⊕, depending on impactor size, all pebbles can reach the core without experiencing significant thermal ablation.

• 2.

  After this point, ablation becomes increasingly severe and a decreasing mass fraction impacts the core directly. The planet’s temperature and pressure are still too low for the envelope to be able to retain a significant amount of silicate vapor, causing the rest to rain out. If the ablated material mixes outwards instead, core growth slows down.

• 3.

  At around 0.50 M⊕, the planet’s envelope mass and temperature have increased sufficiently for absorption of high-Z vapor in the envelope’s inner region to become significant. Impactors can be fully ablated in the atmosphere and direct core impacts terminate. The planet becomes heavily enriched during this third phase, leading to the formation of a high-Z layer around the core with a very high density and temperature.

To further illustrate the core growth process, we have plotted the envelope-to-core ratio of a growing planet along with its absorption rates in Fig. 3. The absorption rate is defined as the fraction of ablated accreting material that stays in the atmosphere in vapor form. The planet’s envelope-to-core ratio starts to increase rapidly after its core mass reaches ~0.5 M⊕. This coincides with the formation of a high-Z layer near the core, as can be seen by the simultaneous increase in the absorption rate. The envelope mass is at first mostly due to hydrogen and helium gasses. When absorption becomes significant, the envelope’s metal rich inner region becomes very hot and dense and its constituting silicate vapor starts to dominate the envelope mass.

The temperature, density, pressure, and high-Z enrichment curves of the resulting planet are plotted in Fig. 4. The gas of the enriched inner region has a much highermean molecular weight in the rainout case, causing it to become very hot and dense.

Fig. 2 Core growth by pebble accretion. Left panel: mass fractions of solids impacting the core directly (green curves) or reaching the core by rainout of ablated material (blue curve). The differently spaced green curves correspond to three sizes of pebbles: 0.1 m (solid), 1 m (dashed), 0.01 m (dash-dotted). Right panel: growth of the core over time. In this figure, the rainout scenario is indicated by the blue curve and the mixing case is shown by the green curve. Both curves in panel b correspond to 0.1 m impactors.

Fig. 3 Absorption fraction (blue) and envelope-to-core ratio (ecr; green) of a planet impacted by 0.1 m pebbles.The green curves indicate the ecr, defined as Menv∕Mcore. The dash-dotted curve is its silicate mass fraction and the dashed curve is its hydrogen and helium mass fraction. They sum to the total ecr, shown by the solid green line. The blue absorption curve indicates the mass fraction of the accreting solids that can be contained by the envelope as vapor.

Fig. 4 Temperature (a), density (b), pressure (c), and high-Z enrichment (d) interior curves of a planet formed by pebble accretion. They are plotted at the point in the planet’s evolution when core growth is equal to 10% of the pebble accretion rate. This corresponds to core masses of 0.58 M⊕ in the rainout scenario (blue curves) and 0.41 M⊕ in the mixing case (green curves). The dotted vertical lines indicate the position of the core.

3.2 Model sensitivities and limitations

To investigate the sensitivity of our results to the variation of several parameters, we performed the same simulations with different pebble sizes, disk orbital radii, solid accretion rates, or SiO ₂ characteristics (see Table 2). The disk conditions for different planet positions correspond to an extrapolation of the conditions at 5.2 AU. We used MMSN scaling, so T_disk α r^−0.5 and ρ_disk α r^−2.75 (e.g., Weidenschilling 1977; Hayashi 1981). The last columns of Table 2 show the final core masses for various sets of parameters. These represent the points at which pebbles are either fully evaporated (mixing) or all ablated material can be absorbed as vapor (rainout).

To zeroth order, the computed final core masses shown in Table 2 can be predicted by the point at which the core surface temperature first exceeds about 1600 K, and SiO2 starts evaporating. Our general picture of limited core growth doesnot change if we vary the pebble size, disk conditions, or the solids accretion rate. The inner atmosphere always reaches 1600 K at a core mass of ~0.3 M⊕. As expected, smaller pebbles evaporate sooner than larger ones, thus leading to a lower core mass in the mixing case. Final core masses after rainout terminates, are less significantly influenced by pebble size because the envelope’s absorption rates stay the same.

These absorption rates do change when the accretion rate is varied, but then final core masses in the mixing and rainout cases do not differ much from 0.4 and 0.6 M⊕ , respectively.This can be explained by noting that increasing the accretion rate of pebbles causes two competing effects on pebble impacts and absorption rates. Firstly, higher accretion rates lead to increased luminosity and thus to higher temperatures, causing more ablation. At the same time however, these higher temperatures reduce gas densities and thus inhibit gas drag. The result is that pebbles fully evaporate at a similar core mass, but at a significantly lower envelope mass (not shown in table). Similarly for absorption rates, higher accretion rates mean that the temperatureat which silicates vaporize starts further from the core, but at a lower density. The same effects also influence the final core mass as a function of position, where a planet’s envelope closer to the central star is hotter and less dense. In our model, final core masses tend to decrease as the planet’s location moves closer to the star.

We find that the core formation process is most sensitive to changes in the mean molecular weight of the high-Z constituent SiO2, which determines the density scaling of the high-Z layer near the core and therefore has a great influence on how much vapor this layer can contain. Variation in this parameter is possible, as the pressure and temperature conditions near the core may allow for SiO2 dissociation (e.g., Medvedev 2016), especially in the presence of hydrogen (e.g., El-Sayed et al. 2015; Soubiran et al. 2017). We vary the mean molecular weight of the high-Z component all the way down to . Our simulations show that in this case, the final core mass of the rainout case increases to 1.49 M⊕ . Other equation-of-state effects are not found to be as significant. Table 2 shows that final core masses do not change when we switch from an ideal gas to a more sophisticated equation of state for a mixture of hydrogen, helium, and SiO2 of Vazan et al. (2013) that does not include SiO2 dissociation. We explain this with the observation that, while the temperatures and densities near the core can exceed those required for non-ideal gas effects such as hydrogen dissociation, this hot inner region consists mainly of silicates (see Fig. 4). We have also included a variation in the adiabatic index because this value is less well known. It is found not to have a large impact on the resultant core masses.

Finally, we performed an additional simulation without considering any enrichment effects, where all ablated material falls to the core. In this case we found that the core grew to 5.1 M⊕, at which point gas accretion became dynamically unstable and our simulation was stopped. This value for the critical core mass is consistent with previous calculations that assumed low (grain-free) opacities (e.g., Mizuno 1980; Hori & Ikoma 2010).

3.3 Icy pebbles

In reality, pebbles outside the iceline are expected to be supplemented by H2O ices (e.g., Schoonenberg et al. 2017; Ida & Guillot 2016). The impactor’s composition is a very important parameter in our core formation model. It determines how rapidly the impactors ablate, how much of the ablated material can be absorbed by the envelope, and what the envelope’s interior region looks like. As a comparison, we repeated the previously described simulations with entirely icy H2O pebbles. We find that these icy pebbles can be fully ablated very early in the planet’s evolution, before the planet’s core has even grown to 0.05 M⊕. This is the case even if we simulate a planet that forms well beyond the H2O iceline at 10 AU at a local disk temperature of 108 K. The subsequent absorption of the water vapor is also more efficient than the absorption of SiO2 vapor. The main determinant of these processes is the vapor pressure, which for water already becomes significant at temperatures of around 200 K, compared to 1600 K for SiO2. This leads to the envelope becoming dominated by water vapor at a low mass (Venturini et al. 2015). If we assume that the inner region of the envelope stays in vapor form, we find that direct core growth terminates before the core reaches 0.1 M⊕. However, we observe that conditions are such in the inner region that the water vapor near the core will liquify, leading to a more complex evolution of the planet, similar to that described by Chambers (2017).

4 Comparison with planetesimal impacts

4.1 Massdeposition

In order to show the main differences in core growth between pebbles and planetesimals, we perform a further analysis on planetesimal mass deposition rates. We consider 1 km planetesimals for two reasons. Firstly, they serve as a limiting case of relatively severe thermal ablation compared to more massive impactors. Secondly, when planetesimals are sufficiently massive, they become dynamicallystable due to their self-gravity. For rocky impactors this can happen when their radius exceeds about 70–100 km (e.g., Podolak et al. 1988; Mordasini et al. 2015). We use a typical planetesimal accretion rate of 10⁻⁶ M⊕yr⁻¹ (e.g., Pollack et al. 1996). Fig. 5a shows the mass fraction that 1 km planetesimals have lost upon core impact. Ablation rates remain limited to one percent even up to a core growth of 2 M⊕ , which is consistent with previous calculations where planetesimals are assumed to reach the core intact.

Instead of ablation, the breakup of planetesimals during impact is typically expected to be the main cause of their mass deposition (Podolak et al. 1988; Mordasini et al. 2015). To calculate breakup points, we checked where in the envelope the dynamical pressure first exceedsthe compressive strength of the planetesimal S_c . Unfortunately, compressive strengths for planetesimals are largely unknown. Modeling of fragmentation events in the Earth’s atmosphere and testing of the impactors post-impact yields compressive strengths in the broad range of 1–500 MPa (e.g., Chyba et al. 1993; Svetsov et al. 1995; Petrovic 2002; Popova et al. 2011; D’Angelo & Podolak 2015). This large variation is due in part to differences in density, size, and the number of pre-existing faults in the rock. Generally, strength decreases with impactor size as the volume expands in which faults can occur (Benz & Asphaug 1999; Stewart & Leinhardt 2011). Around 100 m radius the opposite effect is true, and self-gravity effectively begins to strengthen the object. We use 1 km planetesimal strengths of 1 and 100 MPa for weak and strong rock, respectively. Compressive strength for ice is approximately an order of magnitude below this (Petrovic 2003; D’Angelo & Podolak 2015), and gravitational strengthening, though it still occurs, is lessened by the lower density. The results for weak and strong rock are plotted in Fig. 5b, along with a curve to signify the core radius for comparison. In the weak limit, breakup can happen in envelopes with core masses as small as 0.2 M⊕ . In contrast, strong 1 km planetesimals are found to be resilient to breakup until the core mass exceeds 1.2 M⊕ .

Even if breakup occurs, it does not necessarily imply an instantaneous deposition of the planetesimal’s mass in high-Z vapor. The fragmentation event does not alter the kinetic energy of the post-breakup material until the shape of its bow shock has changed. There are various ways to model this, i.e., by assuming that the impactor flattens out as in the pancake model (Mordasini et al. 2015) or that it expands radially, to subsequently split into smaller fragments when the radius of the post-breakup material becomes too large to support a common bow shock (Hills & Goda 1993). We aim only to get an impression of the importance of this effect. Therefore, we choose a very simple model, similar to Hills & Goda (1993), where it is assumed that the dynamical pressure works to spread the post-impact material radially. Neglecting self-gravity and equating the work done by the dynamical pressure to the time derivative of the kinetic energy of the radially expanding cloud, we find that (35)

where P_dyn is the dynamical pressure. This can be written as (36) (37)

which is very similar to typical expressions of post-fragmentational spreading, the only difference being our factor 5∕4 compared to the use of a constant of order unity as in Mordasini et al. (2015). To get an estimate of the velocity of radial spreading at the breakup point, we substitute P_dyn = S_c, to find (38)

The corresponding timescale to spread out the planetesimal is hence (39)

which may be regarded as the time for a planetesimal’s fragments to seperate. Once initiated, the spreading of the planetesimal will accelerate because of the linear dependence of R_i . We find that these timescales are small compared to the travel time of planetesimals between breakup and core impact. The crosses in Fig. 5b indicate the onset of this condition (t_break + t_spread < t_core). Therefore, (catastrophic) planetesimal fragmentation is likely, at least for small planetesimals. Sufficiently large planetesimals can be strengthened by self-gravity and may directly reach the core, in line with the standard assumption in the literature (e.g., Podolak et al. 1988; Mordasini et al. 2015).

Finally, even when planetesimals fragment and spread, it does not necessarily imply efficient deposition of high-Z (silicate) vapor. When the fragmentation takes place at distances close to the core, which is usually the case (see Fig. 5b), the fragmentation is very localized and our spherically symmetric structure model would greatly overestimate the absorption rates. Indeed, as planetesimals are impacting at velocities close to the escape velocity, the impact cone will be very collimated. Then, it is likely that any high-Z vapor in the impact plume will quickly become supersaturated. In conclusion, we envision that core growth by planetesimal accretion can proceed in a similar way to pebble accretion (with limited direct core growth, see Sect. 3.1), provided that the impacting plantesimals are

• 1.

  sufficiently weak (stronger plantesimals extend the duration of impacts, and thus direct core growth);

• 2.

  sufficiently small (self-gravity will prevent efficient fragmentation beyond radii of 70–100 km);

• 3.

  able to be absorbed efficiently by the envelope (localization of deposited mass in large fragmentation events can prevent this from happening).

A more precise model, which is beyond the scope of this work, should account for these effects to investigate the standard assumptions of planetesimals reaching the core (see Mordasini et al. 2015 for the present state of the art). On the other hand, if we assume that planetesimals reach the core intact (barring ablation), we can simulate the planet’s growth in exactly the same way as the pebble case. We then find that the envelope becomes dynamically unstable due to the accretion of disk gas when the core has grown to several M⊕, i.e., the critical core mass: 3.8 M⊕ for the rainout case and 4.4 M⊕ for mixing¹.

Fig. 5 Mass deposition by planetesimals. Left panel: ablation curves for 1 km planetesimals (red line) and 0.1 m pebbles (blue line) impacting a growing protoplanet, ignoring fragmentation. Planetesimal ablation rates are limited to approximately 1% up to a core mass of 2 M⊕ , whereas the pebbles are fully evaporated before 0.5 M⊕ . Right panel: breakup distances for strong (100 MPa, red solid curve) and weak rock (1 MPa, red dashed curve). The black dashed curve follows the core radius. The crosses represent the points at which the planetesimals have enough timeto radially spread to twice their initial radius. All plotted curves are almost identical for both enrichment assumptions (rainout plotted).

4.2 Core surface conditions

Another way to characterize the envelope of our core growth model is to study the core surface conditions. In Fig. 6 we present a SiO2 phase-diagram, using data from Vazan et al. (2013) and Poirier (2000), to explore these conditions and to verify that our model of a solid core surrounded by a gaseous envelope still holds. The figure shows comparison between 0.1 m pebbles accreting a core up to 0.59 M⊕, and indestructible planetesimals accreting up to 3.8 M⊕. Despite the difference in core mass between the pebble and planetesimal cases, we observe that the conditions at the core reach a pressure of 1 GPa and a temperature of 10⁴ K. In the pebble case, the high-Z inner region of the planet drives up the internal pressure and temperature, whereas in the planetesimal case the envelope contains a much more extended envelope with less enrichment: the envelope is never fully saturated with vapor, as only a small fraction of the impacting planetesimal mass ablates, and the rest reaches the core. Their initial difference is due to the higher accretion rate we use for pebbles, leading to increased luminosity (see Sect. 4).

While the core-envelope boundary conditions remain in the gaseous part of the phase-diagram in both the pebble and plantesimal simulations, the inner pressure and temperature increase rapidly as a function of core mass, especially in the pebble case. The core-envelope boundary conditions approach that of the liquidus towards the end of the simulation at 0.59 M⊕. Extrapolating the trend, this suggests that the high-Z vapor may rain out, resulting in a liquid outer core when the planet continues to accrete more mass.

Fig. 6 Core-surface conditions during core growth by pebble and (indestructible) planetesimal accretion. The core masses in the pebble and planetesimal cases range from 0.1 M⊕ to 0.59 M⊕ and 3.8 M⊕ , respectively. The pluses on the pebble line indicate 0.1 M⊕ core-mass intervals; the crosses on the planetesimal line indicate 1 M⊕ intervals. The locationof the 0.59 M⊕ point is indicated by the circles.

5 Discussion and conclusions

Planet formation by pebble accretion is an alternative to the planetesimal-driven core accretion scenario. The two main differences are a higher expected solids accretion rate and more rapid ablation in the case of pebbles. In this work, we focused on the second difference and analyzed its effect on core growth. To accomplish this, we simulated the early formation of planets for both thepebble and planetesimal-driven core accretion scenarios, using a code consisting of two components: an impact model and a planet evolution model. We find that the rapid ablation of pebbles changes the process of core growth, which means that pebbles can be prevented from directly reaching thecore of even a small growing protoplanet. Instead, their ablated material can either rain out to the core or be absorbed by the envelope as vapor. We have incorporated these enrichment effects from the moment that pebbles first start ablating by modifying the mean molecular weight and gas characteristics locally where high-Z material is present in vapor form.

We find that impacting SiO2 pebbles will fully ablate upon impact with planets of mass less than 0.5 M⊕ . The amount of high-Z material that can stay in vapor form is heavily temperature dependent. This means that only the envelope’s inner region is able to retain this material, leading to the formation of a high-Z layer surrounding the core where temperatures are highest. At first this region is small and unable to contain large amounts of vapor. Most ablated material then rains out to the core. Core growth slows down as the high-Z region expands outwards over time, increasing ablation for subsequent impactors and preventing their mass from reaching the core by absorbing it as vapor. This process limits direct core growth by pebbles to about 0.6 M⊕ . Our findings are relatively insensitive to the planet’s position in the accretion disk, pebble sizes up to 1 m, or the solids accretion rate. They depend most sensitively on the mean molecular weight of the high-Z material and the adopted expression for the vapor pressure (see Eq. (33)). The mean molecular weight determines the steepness of the temperature and density curves in the inner region, while vapor pressure determines its size and influences pebble ablation profiles.

By repeating our simulations with 1 km impacting planetesimals, we have shown that enrichment due to ablation is minimal if the impactors are sufficiently massive. Instead, breakup close to the core becomes an important mass deposition mechanism for impactors that are not sufficiently large to be held together by self-gravity. Plantesimals of 1 km are found to fragment when the planet’s core mass exceeds 0.2–1.2 M⊕, depending on their compressive strength. We estimate that breakup occurs sufficiently far away from the planetary core for the post-breakup material to separate, but the localization of this mass deposition likely still means that only a small portion of the mass can be absorbed by the envelope. This process should be studied in more detail in further work, as the ability of impactor mass to reach the core-envelope boundary is whatdetermines the size to which planetary cores can grow. If efficient absorption of post-fragmentation mass is not prevented by its localization, only planetesimals large enough to prevent fragmentation can directly grow cores of several M⊕.

In our standard model, we use impactors with a uniform rocky (SiO2) composition. In reality, impactors outside the iceline constitute of a mix of rocky and icy pebbles. Pebbles made of water ice will evaporate farther away from the core than rocky ones, after which they can be absorbed by the atmosphere as vapor. We find that this happens before the core has grown to 0.1 M⊕. The inner envelope conditions are such that liquid water can form around the core, as described by Chambers (2017). In the planet’s subsequent evolution, the H2O vapor can either remain in the envelope where it increases the mean molecular weight and speeds up the envelope collapse (Venturini et al. 2015) or it can recycle outwards into the disk (Ormel et al. 2015; Cimerman et al. 2017; Lambrechts & Lega 2017). In both cases, H2O does not contribute to any further core growth. Rocky materials require higher temperatures to vaporize and can therefore continue to grow the core. An implication of this finding is that pebble accretion is only able to form rocky cores beyond 0.1 M⊕, whereas icy cores must thus be formed by impacts of larger, planetesimal-sized impactors. It also offers an explanation for the assertion that close-in planets are predominantly rocky (Owen & Wu 2017; Jin & Mordasini 2018) based on an observed gap in their radius distribution (Fulton et al. 2017).

Recently, Alibert (2017) has suggested that, if recycling operates vigorously and the envelope becomes fully mixed, any high-Z material present in the envelope will recycle back to the disk. In that case, direct core growth is the only growth mechanism in pebble accretion, core masses would be limited to 0.6 M⊕, and enrichment of the envelope would be very low. However, there are two arguments against his hypothesis of efficient mixing. First, recent hydrodynamical simulations involving radiation transport have shown that hydrodynamical recycling operates less efficiently in the inner-most regions where we expect the high-Z vapors to be concentrated (Cimerman et al. 2017; Lambrechts & Lega 2017). Second, the high-Z vapor will congeal when it is transported to cooler regions. These grains will subsequently rain out to the core, unless they can be very efficiently mixed with the outer envelope, i.e., by regular (eddy) convection. Conceivably, luminosity sources from radioactive heating or core cooling (Vazan et al., in prep.) could provide the required luminosity.

However, this mixing is only effective if the planet contains a large convective region. In our work, we find that the convective region is comparatively small and shielded by a large isothermal, radiative zone (see Fig. 4). Furthermore, in order to mix grains into the outer envelope, the (convective) mixing times must be shorter than the settling timescale (τ_mix < τ_settle). Given that the high-Z layer is situated directly on top of the core, settling timescales are very short even for submicron-sized grains. Therefore, we argue that the planet is likely to retain its high-Z vapor against recycling.

Under these conditions the formation of a high-Z layer, situated between the core and the H/He envelope, is an integral part of planet formation by pebble accretion. It can be considered as a dilute extension of the core, with densities approaching the core density and a metal-poor surrounding envelope. Upon further accretion of pebbles, the high-Z region will keep expanding outwards, increasing in mass along with the H/He-region. As it does in the planetesimal-driven scenario, this will end when hydrostatic balance can no longer be maintained. However, accretion of pebbles may alreadyterminate before this point. In that case the protoplanet will evolve through Kelvin–Helmholtz cooling of the envelope (e.g., Lee et al. 2014; Coleman et al. 2017). Because of the cooling, we then expect the high-Z layer to rainout to the core. This is the way we envision core growth to super-Earth and mini-Neptune sizes by pebble accretion.

To summarize, we find that pebble accretion can only directly form rocky cores up to 0.6 M⊕ , and is unable to form icy cores larger than 0.1 M⊕ . This contrasts with planetesimal-driven core accretion, which can directly produce more massive cores with various compositions if the planetesimals are sufficiently strong and large to deposit their mass at the core as is, hitherto, the standardapproach in the literature. The reason for limited direct core growth in pebble accretion is that pebbles ablate more rapidly upon impact. Their subsequent absorption by the surrounding gas prevents their mass from reaching the core when the envelope becomes sufficiently hot and massive. Core growth after this point may proceed through new indirect processes if the planet is able to retain its high-Z material. The expected localization of SiO2 vapors near the core can prevent it from being recycled into the disk and facilitate future core growth when the planet cools down. Further research is required to explore these later stages of core growth, and we will study it in a future work.

Acknowledgements

C.W.O. is supported by The Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422). We thank the referee, Michiel Lambrechts, for a very constructive report.

References

All Tables

All Figures

	Fig. 1 Drag constant evolution of a 0.1 m SiO2 pebble impacting a 0.3 M⊕ planet at escape velocity. The blue line represents the drag constant, the green line indicates the impactor’s velocity relative to the surrounding gas, and the black dashed curve represents the core radius.
In the text

Fig. 2 Core growth by pebble accretion. Left panel: mass fractions of solids impacting the core directly (green curves) or reaching the core by rainout of ablated material (blue curve). The differently spaced green curves correspond to three sizes of pebbles: 0.1 m (solid), 1 m (dashed), 0.01 m (dash-dotted). Right panel: growth of the core over time. In this figure, the rainout scenario is indicated by the blue curve and the mixing case is shown by the green curve. Both curves in panel b correspond to 0.1 m impactors.

In the text

Fig. 3 Absorption fraction (blue) and envelope-to-core ratio (ecr; green) of a planet impacted by 0.1 m pebbles.The green curves indicate the ecr, defined as Menv∕Mcore. The dash-dotted curve is its silicate mass fraction and the dashed curve is its hydrogen and helium mass fraction. They sum to the total ecr, shown by the solid green line. The blue absorption curve indicates the mass fraction of the accreting solids that can be contained by the envelope as vapor.

In the text

Fig. 4 Temperature (a), density (b), pressure (c), and high-Z enrichment (d) interior curves of a planet formed by pebble accretion. They are plotted at the point in the planet’s evolution when core growth is equal to 10% of the pebble accretion rate. This corresponds to core masses of 0.58 M⊕ in the rainout scenario (blue curves) and 0.41 M⊕ in the mixing case (green curves). The dotted vertical lines indicate the position of the core.

In the text

Fig. 5 Mass deposition by planetesimals. Left panel: ablation curves for 1 km planetesimals (red line) and 0.1 m pebbles (blue line) impacting a growing protoplanet, ignoring fragmentation. Planetesimal ablation rates are limited to approximately 1% up to a core mass of 2 M⊕ , whereas the pebbles are fully evaporated before 0.5 M⊕ . Right panel: breakup distances for strong (100 MPa, red solid curve) and weak rock (1 MPa, red dashed curve). The black dashed curve follows the core radius. The crosses represent the points at which the planetesimals have enough timeto radially spread to twice their initial radius. All plotted curves are almost identical for both enrichment assumptions (rainout plotted).

In the text

Fig. 6 Core-surface conditions during core growth by pebble and (indestructible) planetesimal accretion. The core masses in the pebble and planetesimal cases range from 0.1 M⊕ to 0.59 M⊕ and 3.8 M⊕ , respectively. The pluses on the pebble line indicate 0.1 M⊕ core-mass intervals; the crosses on the planetesimal line indicate 1 M⊕ intervals. The locationof the 0.59 M⊕ point is indicated by the circles.

In the text