© ESO 2018

1 Introduction

At early stages of planet formation, coagulation of dust and ice triggers solidification of bodies within the circumstellar gas disk that is present in the first few Myr. Even though the process behind the growth of such bodies is not well understood (Morbidelli et al. 2009; Johansen et al. 2014), if one accepts the core accretion scenario (Perri & Cameron 1974; Mizuno 1980), protoplanets with several Earth masses, M_⊕, have to exist at early disk stages. Once a planetary core reaches a sufficient mass, its gravitational potential dominates the local enthalpy of the surrounding disk gas and the gravitational accumulation of gas into an atmosphere begins.

In order to carry out time-dependent simulations for such atmospheres, it is a practical approach to assume hydrostatic and thermal equilibrium (Stökl et al. 2016). Hydrostatic equilibrium is fulfilled during most phases of atmospheric evolution, but thermal equilibrium is much more difficult to maintain. The absence of energy sources and sinks corresponds to a radially constant energy flow, and therefore a stationary atmospheric model requires a specification of the luminosity of a planet. The common approach in this instance is to associate the planetary luminosity with gravitational energy released by a flow of accreted planetesimals onto the surface of the protoplanet (Hayashi et al. 1979). Despite the plausibility of such a luminosity source, this requires a quantification of the accretion rate and the distributions of planetesimal size and material strength.

The maximum lifetime of a circumstellar gas disk is typically about 10 Myr (Williams & Cieza 2011), which implies a long-term mean mass accretion rate depending on the formation timescale ${\mathcal{M}}_{\text{acc}}$ of the system. Assuming ${\mathcal{M}}_{\text{acc}}$ to be the upper limit of protoplanetary mass accretion, the accretion rate for an individual protoplanet then only depends on the core mass M_C.

To assess the effects of the accretion of solid bodies onto the planet, both the total accretion rate and the size distribution of the incoming bodies must be estimated. Whereas the mass accretion rate is constrained from the formation timescale ${\mathcal{M}}_{\text{acc}}$, a planetesimal size distribution is difficult to obtain. Planetesimals form in the disk by collisional coagulation as dust aggregates drift towards the mid-plane of the disk and form a thin, dense layer of larger bodies (Weidenschilling 2000). Collisional growth continues and produces planetesimals with sizes up to 10 km within a timescale of a few thousand years in the habitable zone. The simulations of Weidenschilling (2000) suggest a size distribution that is dominated by a relatively small number of 100 km bodies after about 1 Myr at a distance of 1 AU. Adopting their results, we finds an power-law particle distribution that has a maximum number density for submillimetre-sized objects and some rare 100-km-sized planetesimals. Knowledge about the size distribution of planetesimals is important because small objects experience more deceleration as they travel through the atmosphere and hence are not able to penetrate into deep atmospheric layers. The amount of deceleration depends on the amount of atmosphere captured by the protoplanet. As a consequence, small objects dissipate their gravitational energy already in the optically thinner parts of the atmosphere and only sufficiently massive planetesimals contribute significantly to the luminosity of the planet in the deep stationary atmosphere. However, large objects are less common and therefore contribute less to impact statistics.

Planetesimal material properties are an important factor for how an object travels through the atmosphere given that the object can in some cases fragment and eventually dissolve. High temperatures close to the protostar in the centre of the early disk dissociate matter into its elements, especially into H, C, N, O, Mg, Fe, Si, and S. As the disk gas is cooling down, those atoms recombine to form molecules that eventually form into dust grains and subsequently into planetesimals (Gautier & Hersant 2005). Observational data, for example from the Infrared Space Observatory (ISO), show that interplanetary dust particles (IDP) mainly consist of silicates (Bouwman et al. 2001; Keller et al. 2002; Molster & Waters 2003; Rietmeijer 2004) and hence an appropriate conjecture for the composition of planetesimals is a stony, silicate mineral.

However, if we take the effects of planetesimal size and strength into account, we might find that the assumption of impact driven, constant luminosity is invalid. In this study, we attempt to obtain further information about how much energy from impacting planetesimals is dissipated in the atmosphere and therefore cannot be released at the planetary surface. We show that in the case of heavy protoplanets, even cooling of the atmosphere can occur. Additionally, we demonstrate that mass accretion of planetesimals is a highly efficient process and hence might require further consideration in future models of protoplanetary primordial atmosphere simulations.

In Sect. 2, we describe our model for the transport of solid bodies through a primordial, planetary atmosphere.In Sect. 3, we verify our model by comparing observational data and other verified models of actual meteorite infall events into the Earth’s atmosphere with our model data. In Sect. 4, we present the results for infalling planetesimals into primordial atmospheres of protoplanets of various sizes. Finally, in Sect. 5 we discuss the outcome of our simulations and give an overview about the scientific relevance of our results.

2 Model description

From the time-dependent models developed by Stökl et al. (2016), we see that primordial planetary atmospheres quickly develop into dense, hot, extended structures. It is interesting to examine the process of objects falling through such dense atmospheres in more detail. Our main aim is to investigate the energy dissipation for a specific type of bolide and to explore how much of the initial mass is lost while travelling through the atmosphere.

Assessing the planetesimal size distribution in the disk is difficult, especially at different stages of atmospheric evolution during the embedded phase. Hence, we have used time-dependent planetesimal simulations (Weidenschilling 2000) and picked a particle diameter spectrum from 1 ≤ D₀ ≤ 10⁷ cm that is plausible at all reasonable disk stages. The spectrum is approximately a power law and is described later in the text. Furthermore, we assume this distribution remains constant during the accretion phase.

Following Perri & Cameron (1974) and Mizuno (1980), we assume that protoplanets with masses between 0.1 and 5 M_⊕ exist within the inner disk early enough in order to maintain an atmosphere up the Hill radius. Moreover, as noted, planetesimal material properties are an important factor for interactions of planetesimals with the atmosphere. Hence, a strength value σ is required to describe the fragmentation of bolides properly. For our study, we set σ = 10⁸ dyn cm⁻², which represents an averaged value for a stony, silicate mineral (Bouwman et al. 2001; Keller et al. 2002; Molster & Waters 2003; Rietmeijer 2004).

2.1 Atmospheric model

When planetesimals cross the Hill radius R_H, they start to interact with the primordial planetary atmosphere of the protoplanet. Therefore, our model adopts the density and temperature stratification of a corresponding atmospheric model. Generally, a protoplanet can accumulate gas from a circumstellar disk into a planetary envelope. As the structure and properties of nebula-embedded protoplanetary atmospheres are an inherently time-dependent problem, Stökl et al. (2016) developed a 1D spherically symmetric hydrodynamics code in order to simulate the accretion process of disk gas onto planetary cores and the subsequent evolution of embedded atmospheres. These authors considered core masses between 0.1 and 5 M_⊕, situated in the habitable zone around a Sun-like star. It should be noted that the structure and composition of primordial atmospheres differ significantly from evolved planetary atmospheres that are no longer attached to the surrounding disk medium. Such early atmospheres are much more massive and can reach extremely high surface temperatures (up to 10 000 K). Figure 1 exhibits the density and temperature profiles for core masses with 0.1, 1, and 3 M_⊕. Such atmospheres contribute several percent to the mass of the planet, whereas the Earth’s atmosphere only adds about 10⁻⁶ %.

The maximum lifespan of a disk is estimated to be ~ 10 Myr or between 1 and 3 Myr for inner disks, respectively (Williams & Cieza 2011). As we consider an embedded protoplanet, we utilized a snapshot at 1 Myr for various protoplanets ranging from 0.1 ≤ M_C ≤ 5M_⊕ from the atmospheric models of Stökl et al. (2016). Choosing snapshots at 1 Myr is, on the one hand, consistent with the existence of a disk (in the range of the common hypothesis of disk lifetime) and, on the other hand, such atmospheres appear to be in a dynamically stable state for a reasonable amount of time; there is no runaway accretion of disk gas over outer boundary and no significant turbulent flows in atmosphere (compare Fig. 1). The parameters that we use from the models of Stökl et al. (2016) are the atmospheric temperature, T_atm, and the atmospheric density, ρ_atm as a function of the radius up to R_H.

In the embedded phase, that is the phase during which the planet is embedded within the disk gas, the relative velocity of planetesimals with respect to the protoplanet is low because of the Keplerian motion of early disks. Hence, we only consider head-on infall events with a small initial planetesimal velocity (relative velocity) v₀ at Hill radius. In order to agree with the literature (Kokubo & Ida 2012), we set v₀ to the individual escape velocity at the surface for each bolide class. Different initial velocities affect the results but our simulations only show negligible variations in the outcome when setting the planetesimal infall speed to their surface escape velocity. We also note that for later phases, when the disk has vanished, the motion of planetesimals within R_H is dominated by the effects of the planet’s gravity, and hence such infall scenarios do not apply correctly. Additionally, post-disk phase infall events require more detailed information about infall velocities and infall trajectories for appropriate simulation results.

Fig. 1 Density (upper panel) and temperature (lower panel) profiles of primordial atmospheres as a function of altitude, for core masses of 0.1 M⊕ (black), 1 M⊕ (green), and 3 M⊕ (red). The dashed lines represent the individual Hill radii of the cores. These data represent a snapshot at 1 Myr from the time-dependent atmospheric evolution simulations and a star with 1 M⊙ (see Stökl et al. 2016).

2.2 Time-dependent model

Consider a bolide with mass m_ptm, at R_H with its initial velocity v₀ for a protoplanetary mass m_pl. As the body accelerates towards the protoplanet’s surface, it collides with atmospheric gas, and the resulting atmosphericdrag has two consequences for the bolide: first, the body heats up due to friction with the atmospheric gas, and second, the body breaks up if the ram pressure exceeds the bolide strength. The exact fraction f of heat transfer from atmospheric gas to the bolide is uncertain, but the upper limit of f is given by f = c_D∕2, where c_D denotes the drag coefficient of the planetesimal (Allen 1962). Podolak et al. (1988) argued that this value may in fact be considerably smaller (f = c_D∕20), and it certainly has to fulfil the condition 0 < f_ene < 1. By analysing observed meteorite falls (see Sec. 3), we found that a value of f_ene = 0.1 best corresponds with the observational results; if f_ene > 0.1 the meteors were burned up too early, and if f_ene < 0.1, they did notshow comparable mass loss when hitting the Earth’s surface.

Friction heating not only causes the body to heat up but also results in a mass loss of surface material via sublimation and via melting and evaporation. This process is called thermal ablation and is discussed in detail by Podolak et al. (1988). Thermal ablation only affects the surfaces of infalling bodies because heat conductivity is slow compared to the ablation rate and thus leaves their core regions almost unaffected. For simplicity, and because the melting and evaporation temperatures of silicates are within the same order of magnitude, we do not distinguish between sublimation and melting followed by evaporation.

Furthermore, we consider all bolides to be spherical and the material properties that we assume are listed in Table 1, where H_f denotes the amount of energy necessary to change a mass of 1 g from the solid to the liquid state, without changing its temperature.

Deriving a set of equations that fully describes the energetic state of the bolide is straight forward. The position of the body depending on the time spent travelling through the atmosphere is given by the velocity of the bolide $$\frac{\text{d}r}{\text{d}t}=v_{\text{ptm}},$$(1)

where r represents the distance from the centre of the protoplanet.

In order to derive the equation of motion, we consider two forces acting on the body. These are the gravitational force F_G and the resistance due to interaction with the atmosphere F_D. Additionally, we assume the bolide velocity vector to point in the same direction as the gravitational force, which leads to $$F_{\text{G}}+F_{\text{a}}-F_{\text{D}}=0,$$

where $$\begin{array}{l}{F}_{\text{G}}=-G\frac{m_{\text{ptm}}M}{r^2}=-GM\frac{\pi D_{\text{ptm}}^3{\rho }_{\text{ptm}}}{6r^2}\\ F_{\text{a}}=-m_{\text{ptm}}a=-\frac{1}{6}\pi D_{\text{ptm}}^3{\rho }_{\text{ptm}}\frac{\text{d}{v}_{\text{ptm}}}{\text{d}t}\\ F_{\text{D}}=\frac{1}{2}{\rho }_{\text{atm}}{v}^2{c}_{\text{D}}{A}_{\text{proj}}=\frac{1}{8}\pi D_{\text{ptm}}^2{c}_{\text{D}}{\rho }_{\text{atm}}{v}_{\text{ptm}}^2,\end{array}$$

where $A_{\text{proj}}=\frac{1}{4}\pi D_{\text{ptm}}^2$ denotes the projected surface of the bolide and $P_{\text{ram}}=\frac{1}{2}{\rho }_{\text{atm}}{v}^2$ is the ram-pressure. Combining these, we can obtain the equation of motion $$\frac{\text{d}{v}_{\text{ptm}}}{\text{d}t}=-\frac{GM}{r^2}+\frac{3}{4}\frac{{\rho }_{\text{atm}}}{{\rho }_{\text{ptm}}}\frac{c_{\text{D}}}{D_{\text{ptm}}}{v}_{\text{ptm}}^2.$$(2)

Concerning the energy balance of bolides and their exchange with the atmosphere, we consider the following processes:

• friction heating,

• radiation of bolide,

• radiation of atmosphere,

• thermal ablation,

whereas friction heating and the radiation of the bolide contribute to the heating of the atmosphere (energy input, E_in) and the radiation of the atmosphere itself plus thermal ablation are cooling and therefore take energy out (E_out) of the protoplanetary atmosphere.

2.2.1 Energy input into the atmosphere

We consider two mechanisms for the heating of the atmosphere by the infalling body. These are radiations from the body and frictional heating. Since bolides radiate approximately as black bodies, the heating rate due to radiation is given by the Stefan–Boltzmann law $$e_{T_{\text{ptm}}}={\sigma }_{\text{S}}{T}_{\text{ptm}}^4\pi D_{\text{ptm}}^2.$$(3)

The frictional heating is given by the drag force multiplied by the speed of the body, i.e. $$e_{\text{D},\text{atm}}=F_{\text{D}}{v}_{\text{ptm}}(1-f_{\text{ene}})=\frac{1}{8}\pi D_{\text{ptm}}^2{c}_{\text{D}}{\rho }_{\text{atm}}{v}_{\text{ptm}}^3(1-f_{\text{ene}}).$$(4)

The kinetic energy that the body loses from the drag force is given both to the atmosphere and to the body as thermal energy. The factor f_ene is the fraction of this energy that is given to the body.

2.2.2 Planetesimal heating

We consider two mechanisms for planetesimal heating. These are from the radiation of the atmosphere and the frictional drag. We estimate the radiative heating by assuming that the atmosphere radiates as a black body, meaning $$e_{T_{\text{atm}}}={\sigma }_{\text{S}}{T}_{\text{atm}}^4\pi D_{\text{ptm}}^2,$$(5)

The heating rate due to the drag force is $$e_{\text{ptm},\text{D}}=F_{\text{D}}{v}_{\text{ptm}}{f}_{\text{ene}}=\frac{1}{8}\pi D_{\text{ptm}}^2{c}_{\text{D}}{\rho }_{\text{atm}}{v}_{\text{ptm}}^3{f}_{\text{ene}}.$$(6)

We only take thermal ablation due to melting the body’s surface regions into account. This gives the energy input rate that corresponds to the mass loss of the bolide as follows: $$e_{\text{m}}=-H_{\text{f}}\frac{\text{dm}}{\text{d}t}=-\frac{1}{2}{H}_{\text{f}}{\rho }_{\text{ptm}}\pi D_{\text{ptm}}^2\frac{\text{d}{D}_{\text{ptm}}}{\text{d}t}.$$(7)

Combining Eqs. (3), (5)–(7) and considering the fact that energy input and energy output have to be in equilibrium, we can derive an equation that represents the diameter change due to mass loss of the bolide as follows: $$\frac{\text{d}{D}_{\text{ptm}}}{\text{d}t}=-\frac{1}{{\rho }_{\text{ptm}}{H}_{\text{f}}}\underset{(\star )}{\underbrace{\left[\frac{1}{4}{\rho }_{\text{atm}}{c}_{\text{D}}{f}_{\text{ene}}{v}_{\text{ptm}}^3+2{\sigma }_{\text{S}}\left(T_{\text{atm}}^4-T_{\text{ptm}}^4\right)\right]}}.$$

As the surface of the body can only heat up until it reaches melting temperature before mass loss occurs, the surface temperature of the bolide has to be limited by its melting temperature T_m. Therefore, we set (⋆) = 0 (boundary condition) and write for T_ptm $$T_{\text{ptm}}=\min \left({\left[T_{\text{atm}}^4+\frac{1}{8{\sigma }_{\text{S}}}{\rho }_{\text{atm}}{c}_{\text{D}}{f}_{\text{ene}}{v}_{\text{ptm}}^3\right]}^{\frac{1}{4}},T_{\text{m}}\right).$$

Furthermore, we assume that the infalling material is composed of silicate, hence we set the melting temperature in our simulations to T_m = 1800 K (Podolak et al. 1988).

In addition to thermal ablation, the body can also break up if the ram pressure exceeds the internal bolide strength (strength criterion). We assume that break-up is a discrete event where the body separates into a certain specified number of fragments n_sep whenever the strength criterion is fulfilled. In our simulations, we set this fragmentation number to the lower limit of n_sep = 2. After a break-up, the fragments do not act as individual bodies right away but only after a period of time in which the fragments drift apart (distance criterion). We consider that a body has broken apart when $${\sigma }_{\text{ptm}}\le \frac{1}{2}{\rho }_{\text{atm}}{v}_{\text{ptm}}^2$$

and we consider the bodies to be separate bodies when $$x_{\text{H}}\ge f_{\text{dist}}{D}_{\text{ptm},}$$

where x_H is the horizontal distance between the fragments and f_dist is an arbitrary factor that describes how far the fragments need to drift apart (in bolide diameters) until they are seen as individual objects. Obviously, an accurate description of how the fragments drift apart is not trivial (see e.g. Shuvalov et al. 2014). As our model does not need detailed information about the drift process, but only needs information about the drift velocity of the fragments, we followed Shuvalov et al. (2014) who suggested that a good assumption for the horizontal drift velocity is $$v_{\text{H}}=\left|v_{\text{ptm}}\right|\sqrt{\frac{{\rho }_{\text{atm}}}{{\rho }_{\text{ptm}}}}$$

and a distance parameter of f_dist = 2 is a reasonable value to treat fragments as separate objects.

After complete break-up, we have n_sep bolides with diameters of $D_{\text{ptm}}=D_{\text{ptm},\text{0}}/\sqrt[3]{n_{\text{sep}}}$, where D_ptm,0 denotes the bolide diameter before break-up. Furthermore, since the horizontal drift velocity is evaluated by the radial bolide velocity, it can be easily seen that the change of v_H is very small compared to the time it takes for the fragments to achieve the distance criterion. Thus, we can assume that it is constant during the drift event, which implies that the diameter change can be assumed to be linear with separation distance x_H. Hence, the actual bolide diameter is $$D_{\text{ptm}}(t)=D_{\text{ptm},\text{0}}-\frac{1-n_{\text{sep}}^{-1/3}}{f_{\text{dist}}}{x}_{\text{H}},$$

with the corresponding time derivative $$\frac{\text{d}{D}_{\text{ptm}}}{\text{d}t}=-\frac{1-n_{\text{sep}}^{-1/3}}{f_{\text{dist}}}{v}_{\text{H}}.$$

Adding this finding to the equation of bolide mass loss we get $$\begin{array}{ll}\frac{\text{d}{D}_{\text{ptm}}}{\text{d}t}=\hfill & -\frac{1}{{\rho }_{\text{ptm}}{H}_{\text{f}}}\max \left[\frac{1}{4}{\rho }_{\text{atm}}{c}_{\text{D}}{f}_{\text{ene}}{v}_{\text{ptm}}^3+2{\sigma }_{\text{S}}\left(T_{\text{atm}}^4-T_{\text{m}}^4\right),0\right]\hfill \\ \hfill & -\frac{1-n_{\text{sep}}^{-1/3}}{f_{\text{dist}}}{v}_{\text{H}}.\hfill \end{array}$$(8)

2.2.3 Final set of equations

Summarising all the results from this section the set of equations that have to be solved takes the following form $$\begin{array}{ll}\frac{\text{d}r}{\text{d}t}=\hfill & v_{\text{ptm}},\hfill \end{array}$$(9) $$\begin{array}{l}\frac{\text{d}v}{\text{d}t}=-\frac{GM}{r^2}+\frac{3}{4}\frac{{\rho }_{\text{atm}}}{{\rho }_{\text{ptm}}}\frac{c_{\text{D}}}{D_{\text{ptm}}}{v}_{\text{ptm}}^2,\\ \frac{\text{d}{D}_{\text{ptm}}}{\text{d}t}=-\frac{1}{{\rho }_{\text{ptm}}{H}_{\text{f}}}\max \left[\frac{1}{4}{\rho }_{\text{atm}}{c}_{\text{D}}{f}_{\text{ene}}{v}_{\text{ptm}}^3+2{\sigma }_{\text{S}}\left(T_{\text{atm}}^4-T_{\text{m}}^4\right),0\right]\\ -\frac{1-n_{\text{sep}}^{-1/3}}{f_{\text{dist}}}{v}_{\text{H}},\\ \frac{\text{d}{x}_{\text{H}}}{\text{d}t}=\{\begin{array}{ll}{v}_{\text{H}}\hfill & \text{at break-up,}\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}\end{array}$$

which can be solved by numerical integration. We adopted a 4th order Runge–Kutta scheme with adaptive step-size control as described in Press et al. (2007, page 554).

2.3 Energy dissipation into a protoplanetary atmosphere

Since we are interested in the total energy dissipation into the atmosphere, further considerations are necessary. The total energy dissipation is given by the sum of the energies dissipated by each of the bodies that travel through the atmosphere. Since the size and mass of an object determines how much energy it brings, it is necessary to consider the distribution of object sizes that enter the atmosphere. This distribution evolves over time and has typically many millimetre-sized and a few kilometre-sized objects. In order to describe the particle size distribution, we followed Weidenschilling (2000), whose simulations show that because of the growth of solids, the distribution is approximately following a power law. Using this result, we extract a number density per bolide diameter, per volume of surrounding disk material N(D_ptm) $$\log (N(D_{\text{ptm}}))=k\log (D_{\text{ptm}})+d,$$(13)

where k is the slope and d is the offset in the logarithmic size distribution. Referring N(D_ptm) on an initial bolide diameter D_ptm,0 we get $$N(D_{\text{ptm}})=N_{\text{ptm},\text{0}}{\left(\frac{D_{\text{ptm}}}{D_{\text{ptm},\text{0}}}\right)}^k.$$

Adopting discrete mass bins (given by the bolide-diameter grid, where the index i is used to indicate the quantity is for the ith diameter bin), the total mass per bolide class per volume of surrounding disk material is given by $${\rho }_{\text{ptm},i}=N_{\text{ptm},i}{m}_{\text{ptm},i}=\frac{1}{6}\frac{N_{\text{ptm},\text{0}}}{D_{\text{ptm},\text{0}}^k}\pi {\rho }_{\text{B}}{D}_{\text{ptm},i}^{k+3}.$$

Mass accretion onto the protoplanet depends on the number of planetesimals within a given distance of the orbit of the planet. We estimate that the total mass of planetesimals that can effectively be accreted by the protoplanet is distributed in an annulus around the protostar. According to the work of Greenzweig & Lissauer (1990), the maximum radial extend of the accretion zone due to gravitational focussing is approximately 4 R_H (radial component) and the scale height H_p of the planetesimal disk can be estimated by H_p ~ 0.01 H (Lambrechts & Johansen 2012), where H ~ 0.05 R_P is the scale height of the gaseous disk (Armitage 2011). This leads to a volume of the annulus of $$V=0.004\pi R_{\text{P}}^2{R}_{\text{H}},$$

where R_H is the Hill radius and R_P = 1 AU is the distance of the protoplanet to the protostar. Multiplying this by the individual densities gives the total masses within each diameter bin within the annulus. Summing the masses of all mass bins gives the total mass that is available for accretion onto the protoplanet as follows: $$m_{\text{tot}}={\sum }_i{M}_{\text{ptm},i}={\sum }_i{\rho }_{\text{ptm},i}V=\alpha \frac{N_{\text{ptm},\text{0}}}{D_{\text{ptm},\text{0}}^k}{\sum }_i{D}_{\text{ptm},i}^{k+3},$$

where $\alpha =2{\pi }^2{\rho }_{\text{B}}{R}_{\text{P}}^2{R}_{\text{H}}/3000$.

Mass accretion does not happen instantaneously but instead happens over a certain period of time. Thus, we can define a (mass-bin specific) frequency factor ν_i by dividing the planet’s mass accretion rate ${\mathcal{M}}_{\text{acc}}$ by m_tot and multiplying with the weighted factor β_i = m_ptm,i∕m_tot, which gives $${\nu }_i=\frac{{\mathcal{M}}_{\text{acc}}}{m_{\text{tot}}}\frac{m_{\text{ptm},i}}{m_{\text{tot}}}=\frac{{\mathcal{M}}_{\text{acc}}}{m_{\text{tot}}}{\beta }_i,$$(14)

where ν_i describes the total infall rate per mass bin. The mass accretion rate for a specific planetary core can be determined by multiplying the core mass M_C (in Earth masses) with the mean mass accretion rate ${\mathcal{M}}_{\text{acc},\oplus }={10}^{-7}{M}_{\oplus }{\text{yr}}^{\text{-1}}$.

The total amount of energy dissipated in the protoplanetary atmosphere is given by the frictional heating plus the energy exchange rate between the bolide and the atmosphere (Sect. 2.2). As the object might break up while travelling through the atmosphere, the number of fragments must be taken into account. This value can be obtained by simply multiplying the number of fragments by the individual energy input rates per bolide mass bin, $$\begin{array}{l}{e}_{\text{fric},i}=N_{\text{ptm},i}\left(1-f_{\text{ene}}\right)\frac{1}{8}{\rho }_{\text{atm}}\pi c_{\text{D}}{v}_{\text{ptm},i}^3{D}_{\text{ptm},i}^2\\ e_{\text{rad},i}=N_{\text{ptm},i}{\sigma }_{\text{S}}\left(T_{\text{ptm},i}^4-T_{\text{atm}}^4\right)\pi D_{\text{ptm},i}^2\\ e_{\text{diss},i}=e_{\text{fric},i}+e_{\text{rad},i},\end{array}$$

where the number of fragments N_ptm,i is given by $$N_{\text{ptm},i}=\frac{6m_{\text{ptm},i}}{D_{\text{ptm},i}^3\pi {\rho }_{\text{ptm}}},$$

which can be obtained from the individual bolide mass using $$\begin{array}{lll}\frac{\text{d}{m}_{\text{ptm},i}}{\text{d}t}\hfill & =\hfill & -\frac{\pi }{H_{\text{f}}}{D}_{\text{ptm},i}^2\hfill \\ \hfill & \hfill & \times \max \left[\frac{1}{8}{\rho }_{\text{atm}}{c}_{\text{D}}{f}_{\text{ene}}{v}_{\text{ptm},i}^3+{\sigma }_{\text{S}}\left(T_{\text{atm}}^4-T_{\text{m}}^4\right),0\right].\hfill \end{array}$$(18)

It should be noted that to derive Eq. (18), melting in terms of thermal ablation is the only mass loss process we consider (no mass loss on break-up), so the above equation represents the conservation of mass for a particular bolide class.

Equations (15)–(17) describe an energy input rate, hence, to compute the corresponding total amount of dissipated energy per bolide mass bin, these equations have to be integrated over time. For integration, we transform Eq. (17) to spatial coordinates according to d ∕d r = v⁻¹d∕dt. Doing so, results in the total amount of dissipated energy by a single bolide mass bin being given by $$\begin{array}{l}{E}_{\text{tot},i}={\displaystyle {\int }_R{N}_{\text{ptm},i}}\\ \underset{(\star \star )}{\underbrace{\left[\left(1-f_{\text{ene}}\right)\frac{1}{8}{\rho }_{\text{atm}}{c}_{\text{D}}{v}_{\text{ptm},\text{i}}^2{D}_{\text{ptm},\text{i}}^2+\frac{{\sigma }_{\text{S}}}{v_{\text{ptm},\text{i}}}\left(T_{\text{ptm}}^4-T_{\text{atm}}^4\right)\right]}}\text{d}r,\end{array}$$(19)

where the term in braces is the total force per bolide diameter. We integrate this equation numerically using a Simpson integration scheme.

In order to obtain the total amount of energy input into the atmosphere, the individual total energies (Eq. (19)) have to be summed. Multiplying by Eq. (14), leads to the energy input rate evaluated with the distribution function and mass accretion rate, given by $$e_{\nu }=\nu E_{\text{tot}}={\sum }_i{\nu }_i{E}_{\text{tot},i}.$$

We note that when we use too few mass bins, e_ν depends on the number of bins, which is undesirable. As we increase the number of bins, our calculated e_ν changes until N_MB > 700, at which point we find that e_ν is independent of N_MB. We therefore assume N_MB = 701, which give us an error less than 10%, which can be considered as an appropriate result in the limit of our model. Therefore, the computations in the next sections are performed with N_MB = 701.

3 Model verification

Even though bolide dynamics in the Earth’s atmosphere is not the primary objective of our model, observed meteorite falls are one valuableway to validate the reliability of our approach. Therefore, a number of Earth-related meteorite events with good observational records and/or simulation data were chosen and are listed in Table 2. Some of these data are based on estimates as neither remains of the objects were found nor any video material was available (e.g. Tunguska). Furthermore, as stated by some authors (Chyba et al. 1993; Borovicka & Kalenda 2003), even if video material is available, it is difficult to determine exact bolide parameters due to image saturation and/or fuzziness of the image.

Table 2 shows the estimated initial parameters of the test meteorites, where R₀ is the initial height of the object and mainly corresponds to the starting point of observation, M₀ the initial mass, v₀ the initial velocity, ρ the density, and σ is the strength of the meteorite. The other parameters are h_AB, which represents the observed or modelled airburst height, N_frag is the number of recovered fragments, and M_rec is the total recovered mass of the bolide. These are the main indicators for our model verification.

Fig. 2 Airburst heights of the different meteorites; plot shows the number of fragments Ni as a functionof altitude. A number of fragments, Ni = 2, indicate that both the strength criterion and distance criterion have been fulfilled for the first time. A number of fragments between 1 and 2 indicate that the strength criterion has been saisfied, but the distance criterion has not, and therefore weconsider that the body has broken apart, but the fragments have not separated far enough to be considered separate bodies.

3.1 Verification results

There is no observational record (apart from seismic wave data and physical destruction of the landscape in the impact area) of the Tunguska event and therefore the initial parameters have to be chosen seperately. Because of this, we followed Chyba et al. (1993) where different impact scenarios of the meteorite were studied. As our model only supports the perpendicular motion of bolides through the atmosphere, we only looked at those simulations with comparable set-ups. As there is no specific starting point mentioned for the Benešov bolide, we assume the initial height to be at Hill radius R_H in our computations.

For the model of the Earth’s atmosphere, we used the International Standard Atmosphere (ISA). Since ISA is only defined up to a height of 86 km and as there are some bolides with an initial height greater than that value, we isothermally extrapolated the atmospheric ISA values up to R_H. In our case, the choice of the atmospheric model (ISA, NRLMSISE, etc.) is not important because the density profiles in the lower atmosphere (up to 50 km) are comparable and densities above that altitude are low enough not to affect the meteorites aerodynamically (see Fig. 2).

Figure 2 shows airburst heights for the individual test candidates. The value N_i = 2 indicates the first complete separation of two fragments (see Sect. 2.2.2). One could assume that the actual fragmentation has already started as soon as the strength criterion is fulfilled (as soon as N_i > 1). But as the drift velocity of the fragments scale with the radial velocity of the bolide, the separation criterion is fulfilled soon after the strength criterion, compared to the distance the meteorite has travelled through the atmosphere. Hence, we assume a fragmentation number of 2 to be the indicator forthe airburst event.

Comparing numerical results computed by our model with airburst heights stated by the various authors listed in Table 2, we can see that they are in good agreement, indicating that our model is a valid description of the actual physical processes. Almost all candidates show a break-up within the observed height or the values are close to the comparison models. Only the Morávka main body and the Benešov meteorite are slightly off. In the case of the Benešov event, Artemieva & Shuvalov (2001) referred to a reference particle with a slightly higher density of ρ = 4.2 × 10⁷ dyn cm⁻², which mightaffect the airburst height of our model. The computational results are listed in Table 3 and the highest deviation from literature data for an airburst is registered for the Morávka main body with 3.1%, which is still a good result considering a totally different model approach.

Comparing observational data other than airburst height must be treated with caution. This is because the number of found fragments N_frag or the total recovered mass M_rec are not reliable numbers as it is unlikely that all bolide fragments were found and recovered. In our case, it can be a good indicator formodel verification in terms of whether a bolide impacts the surface or completely dissolves in the atmosphere, as we are not interested in detailed impact scenarios for the various meteorites.

Figure 3 shows the bolide diameter, velocity, and mass while travelling through the atmosphere. All values are normalised to their initial values to obtain a more assessable plot. Steep slopes in the diameter plots represent runaway break-up, when multiple break-up events rapidly take place. This means that the break-up condition (strength criterion) is again fulfilled right after the fragments reach separation distance (distance criterion) and consequently trigger another fragmentation process. If the individual fragments cannot slow down quickly enough this process continues until the whole body completely dissolves in the atmosphere.In Fig. 3, the Benešov event clearly shows these characteristics. In the simulation, the bolide diameter drops quite rapidly down to submillimetre-sized dust particles that remain suspended in the atmosphere and eventually rain out after a certain amount of time. This corresponds to a fall in which no fragments of the meteorite were found.

The Tunguska event clearly shows an impact of the bolide with a rather high velocity. Since the initial size of the bolide was large and the initial velocity rather low, there was no significant mass loss at the time of impact. Our model computes an impact energy (kinetic energy) of 4.85 × 10²³ erg that is in good agreement with the literature where the event is liberating between 4 × 10²³ erg (Hunt et al. 1960)and 4 × 10²⁵ erg (Turco et al. 1982).

Our model suggests that in both of the Morávka events the remaining fragments would be around millimetre size. This does not agree with the observational data, but, as stated before, we can only simulate perpendicular motion, which does not compare with the literature as the trajectory is supposed to be around 20° (Borovicka & Kalenda 2003). Furthermore, Borovicka & Kalenda (2003) describe the body as rather flat and this causes a higher deceleration in the upper atmosphere. In this case, our method does not allow for a detailed description of bolide behaviour after the first fragmentation event is not possible. Nevertheless, our model predicts a total recoverable mass of 760 kg, distributed over a large amount of fragments for the Morávka event. Simulation results for the Morávka main body (Borovicka & Kalenda 2003) predict a total mass of 17.1 kg at an altitude of 21.5 km. Comparing this result with our model data, we see that our model predicts a total mass of about 25 kg at the same altitude, which seems to be an acceptable result considering a completely different model approach.

Looking at the Neuschwanstein meteorite our model clearly predicts an impact that is in accordance with the literature. Our model suggests a total mass, possibly hitting the Earth’s surface, of about 1.03 kg. This is a bit less than the mass of the recovered object with 1.75 kg (Spurný et al. 2003). Again the authors describe the atmospheric trajectory of the bolide to be 49.5°, which will have effects on the mass loss of the body.

The Bunburra Rockhole meteorite fall is a well observed event by the Australian Desert Fireball Network (DFN; Spurný et al. 2012). Our model suggests a catastrophic break-up for the meteorite down to submillimetre-sized particles. Considering the weak composition of the body, this again seems to be a reasonable result. Moreover, the heavy fragmentation triggers a large deceleration of the individual fragments, leading to a negligible amount of friction heating and thus no mass loss from thermal ablation. This seems to be in good agreement with the literature because only three rocks with a total mass of 339 g were found. Furthermore, the fragmentation events of the meteorite are described as “explosive,” which substantiates very small fragments predicted by our model. Our model suggests a total mass, possibly hitting the surface, of ~20 kg distributed over a large amount of submillimetre particles.

Summarising the above results, we emphasise that our model lacks a detailed reproduction of fragmentation after the initial break-up although it produces acceptable results with airburst heights for all cases. It does not allow us to analyse different atmospheric trajectories and does not take a shape factor into account, which has an impact on the aerodynamics of the object. It is not our aim in this paper to analyse detailed impact scenarios of bolides, but our focus lies on the energy dissipation in dense, protoplanetary atmospheres. With that in mind, our results are robust with respect to numerical and physical parameters and even deliberately choosing extreme values for planetesimal parameters only moderately affects the results.

Fig. 3 Combined plot for various meteorites. Panel a: normalised diameter reduction while travelling through the Earth’s atmosphere,panel b: normalised velocity, and panel c: normalised mass loss as a function of altitude.

4 Results for protoplanetary atmospheres

Our results for infalling planetesimals into primordial protoplanetary atmospheres reveal several interesting outcomes, some of which are seemingly counterintuitive. This underlines the demand of a reconsideration of the classical idea that all the kinetic energy of planetesimals is deposited on the surface of the planet. In this section, we apply our model to the case of bodies travelling through the thick protoatmospheres of planets embedded in the disk. We run a series of infall models for various body sizes and planetary masses.

4.1 Effects on planetesimals

Intuitively one might think that the more massive an infalling body, the more likely it will impact on the planetary surface and the higher its impact velocity, but inspecting Fig. 4, the reality appears to be more complicated.

Figure 4 shows a simulation for bolides with initial diameters ranging from 1 ≤ D₀ ≤ 10⁷ cm for a planetary core mass with 1 M_⊕. All panels show eight logarithmically distributed test bodies starting from the lower limit (1 cm; black lines) to the upper limit (100 km; redlines). Even though our simulations were carried out all the way to the planetary surface, we cut all plots in Fig. 4 at an altitude of 10 km since below this altitude, there are no significant changes in the variables.

In Fig. 4a, the 10 km test bolides eventually cross lines with bolides with lower masses and diameters, which implies that initially larger bodies do not always result in a higher impact velocity and vice versa. We see that small objects up to a size of 10 cm experience adeceleration already in the outer, rather tenuous atmosphere. Once the particle reaches the hot and dense layers (at around h = 10⁹ cm, cf. Fig. 1), gravitational acceleration and deceleration by friction, balance out and cause a flattening in the velocity curve. The relatively low speeds (~ 10 m s⁻¹) cause the particle to remain in the atmosphere’s high-temperature regions longer, which leads to a heating of the surface of the bolide and eventually, once the surface reached melting temperature (Fig. 4d), to a complete dissolution by thermal ablation (Figs. 4b and c).

Bodies with initial diameters from about 1 m to approximately 10 m are less affected by the outer atmosphere, thus reaching higher velocities. The higher velocity results in more friction heating leading to a higher surface temperatureof the bolide. When reaching the inner, dense layers of the atmosphere, the planetesimal experiences a large deceleration. At this point, the temperature of the body is higher than the surrounding gas temperature and the decreasing efficiency of friction heating causes a radiation cooling of the bolide (negative gradients in Fig. 4d). As the velocity decreases further and the atmospheric temperature rises, the surface of planetesimal reaches melting temperature and, similar to smaller objects, the bolide disintegrates owing to thermal ablation.

Planetesimals with sizes around 100 m ≤ D₀ ≤ 1 km are less affected by a deceleration in the outer atmosphere, but they experience more friction heating because of their larger effective surfaces. The velocity of such particles (in the denser atmosphere) is not enough to cause fragmentation but also high enough to cause significant friction heating and therefore extensive mass loss (Fig. 4c).

Large planetesimals (10 km ≤ D₀ ≤ 100 km) do not experiance much deceleration in the outer atmosphere, but again more friction heating owing to a larger cross section of the body. As one can see in Fig. 4d, such bolide surfaces reach melting temperature rather early but owing to their high velocities, they do not remain in the hot regions long enoughfor a significant mass depletion (Fig. 4c). Comparing the 10-km-sized with the 100-m-sized test bolide, we note that the 10 km object experiences a heavy deceleration to almost the same final velocity as the 100 m object. This behaviour is again attributable to the size of the object as the larger effective surface causes the planetesimal to experience a runaway break-up until the individual pieces are slow enough to dissatisfy the strength criterion (Sect. 2).

If we investigate corresponding impact velocities (Fig. 5), we see again that bolide size is not necessarily an indicator for impacts. Some planetesimals are affected by a runaway break-up and might eventually dissolve in the atmosphere, whereas others are small enough to never exceed the strength criterion and, at the same time, are large enough so that mass loss due to thermal ablation is negligibly small.

Certainly, these processes depend very much on the structure and size of the protoplanetary atmosphere. Heavy planetary cores are surrounded by denser atmospheres that cause an earlier break-up for two equally sized bodies. Figure 5 shows impacting bodies for three protoplanets 0.1, 1, and 3 M_⊕. Shaded areas represent the range of impactors at the protoplanetary surface with a speed higher than drift velocity v_drift. The drift velocity is used as a criterion to stop the simulation when a particle reaches such low speeds that the time it remains in a quasi-isothermal atmospheric layer is comparable to the timescale of heat conduction. Protoplanetary atmospheres near the surface reach temperatures much higher than melting temperatures of silicates (Stökl et al. 2016), which cause a stony particle travelling at drift velocity to completely melt at one point (Fig. 4c; bolide sizes up to some 10 m). Additionally, it appears to be reasonable that convective processes in such atmospheres are able to capture rather lightweight particles and trap them long enough to completely disintegrate the body. According to that, we set v_drift = 1 cm s⁻¹, which seems to be an acceptable value taking the above statements into account. For the results represented in Fig. 5 we ran simulations with the same, previously defined bolide grid but using a much finer grid of planetesimal sizes (ΔD = 10⁻² cm).

A protoplanet with 0.1 M_⊕ holds a much thinner atmosphere, at the same evolutionary phase, than a core with 1 M_⊕. The smaller core mass accumulates less atmosphere during its evolution and has a smaller Hill radius. Therefore, not only is the initial gravitational energy for equally sized bolides less than for a larger protoplanet, but atmospheric drag is also weaker as a result of thethinner atmosphere. We see that in this case, the full range of planetesimals reaches the surface of the planet. Whereas large bodies (10–100 km) are almost unaffected, small bodies experience a deceleration to about v_ptm ~ 10³ cm s⁻¹ (Fig. 5a).

For a protoplanets with 1 M_⊕, small bolides with sizes up to D₀ ~ 45.7 m, completely disperse in the dense, inner atmosphere (Fig. 5b). Larger planetesimals, up to the upper limit of the planetesimal grid, reach the surface either as a whole or as fragments with a smaller impact velocity (compare with Fig. 4a).

Figure 5c shows possible impact candidates for a protoplanet with 3 M_⊕. As a result of a more extended, denser atmosphere, we are only left with a small gap of planetesimals ranging from D₀ = 13.2 km to D₀ = 28.2 km. As stated before, we count bodies as impactors if their surface velocity is higher than v_drift. In the case of this scenario, all bodies (and their fragments) show relatively low impact speeds and hence cannot efficiently contribute to a core luminosity as the majority of their gravitational energy is already dissipated in the atmosphere.

Fig. 4 Figure comparing bolide velocity (panel a), diameter (panel b), mass loss (panel c), and temperature (panel d) as a function of altitude in the atmosphere for a planetary core mass of 1 M⊕ for various initial diameters ranging from D0 = 100 cm to D0 = 107 cm for one particular arbitrary chosen protoplanetary atmosphere.

4.2 Effects on protoplanetary atmospheres

Our main interest is the amount of energy that is dissipated while such objects travel through a protoplanetary atmosphere. Figure 6 shows therelative energy dissipation for the same distribution of infalling bolides into a protoplanetary atmosphere with a planetary mass of 1 M_⊕. The grid is represented by N_B = 701 mass bins with 1 cm (purple) ≤ D₀ ≤ 10⁷ cm (yellow). We plot the relative total energy dissipation per mass bin, meaning that these results are not affected by the size distribution function. Small planetesimals (up to some 1 m; purple to blue lines) reach their maximum energy dissipation limit higher up in the atmosphere as they are more affected by atmospheric drag. Once they reach the inner dense layers, those objects start to disintegrate owing to thermal ablation. Since the transfer from the solid to liquid state of the bolide’s surface layers requires a certain amount of energy (H_f, Sect. 2), the energy input into the atmosphere reduces by around 20% for centimetre objects but there is no reduction for metre-sized objects. This finding represents the previously discussed, first gap (Fig. 5b) of non-impactors. Considering that such objects are small, the energy that is necessary for complete fusing is small and thus leading to a positive budget in energy input.

Larger bodies (D₀ > 4.57 × 10³ cm, Fig. 5b) appear with a negligible reduction of total energy input until we reach the lower limit of the second gap of non-impactors (D₀= 1.32 × 10⁶ cm). This time, the total mass of the bolide is much higher and therefore much more energy for complete fusing is necessary. In the case of this specific bolide range (1.32 × 10⁶ ≤ D₀ ≤ 2.82 × 10⁶ cm), the total energy input is negative, resulting in an effective cooling of the atmosphere.

Once planetesimal size passes that second gap, the energy budget remains positive again. Even though those planetesimals break up into individual bodies, fragmentation occurs lower in the atmosphere and the individual fragments are fast enough to reach the surface without getting completely disintegrated.

Despite the fact that some of the large planetesimals locally extract energy from the atmosphere because of melting, the total energy input (folded with the size distribution function and evaluated with the relative mass accretion rate; Sect. 2) for an M_C = 1 M_⊕ core is still positive (see Fig. 7a, blue line) as the energy input is dominated by the largest bolides. Figure 7a shows the relative total energy input rate e_ν,tot (see Sect. 2) divided by the maximum energy input rate e_ν,max as a function of altitude for four core masses: 0.1 M_⊕ (black), 1 M_⊕ (blue), 3 M_⊕ (green), and 5 M_⊕ (red). Whereas e_ν,tot is overall positive for a planetary core with M_C = 1 M_⊕, this is not true for larger planetary cores. Heavier cores hold more atmosphere and consequently narrow the range of possible impactors (Fig. 5c) because of fragmentation. Hence, more planetesimal fragments remain in the atmosphere and heat up. Once the surface of the body reaches melting temperature, energy is required for phase transition, thus leading to a local cooling of the atmosphere. Our simulations show that protoplanetary atmospheres with cores larger than M_C > 3 M_⊕, capture 100% of infalling planetesimals, which results in a significant cooling of lower atmospheric layers.

Finally, we discuss the results on how much energy is dissipated within the different layers of the atmosphere. Figure 7b shows the total heating rate e_ν,tot in values of e_ν,max with respect to the optical depth τ, for core masses with M_C= 0.1 M_⊕ (black), M_C= 1 M_⊕ (blue), M_C= 3 M_⊕ (green) and M_C= 5 M_⊕ (red). We see that, regardless of core size, the maximum energy is dissipated in the optical thick part of the atmosphere, which implies that the interaction of planetesimals with protoplanetary atmospheres is an important factor of atmosperic evolution; we note that we have cut the plot at τ = 10 as below that no significant changes in the values appear. Whereas with small planeteary cores most of the initial gravitational energy of the planetesimals is released as kinetic energy on impacts (compare Fig. 5), atmospheres around large cores cause heavy fragmentation even resulting in a local cooling of the atmosphere. As the optical thick layers of the atmosphere are fully convective, the atmosphere must react with expansion in case of heating and with contraction in case of cooling.

Fig. 5 Diameter range of impacting planetesimals as a function of altitude for protoplanet core masses of 0.1 M⊕ (panel a), 1 M⊕ (panel b), and 3 M⊕ (panel c) as a function of altitude. Grey areas denote ranges between the lower and upper limit of impacting mass bins (solid lines). The plot shows bolide velocities vi with respect to altitude h with an impact velocity vh=0 > vdrift. All simulations are carried out with a grid of initial bolide diameters ranging from 1 ≤ D0 ≤ 107 cm.

Fig. 6 Energy dissipation per mass bin in units of Emax into the atmosphere of a protoplanet with 1 M⊕. The colours represent a bolide grid ranging from D0 = 1 cm (purple) to D0 = 100 km (yellow) with NB = 701 mass bins. The total energy input is dominated by large bodies. Some bolide cool the atmosphere owing to phase transitions from solid to liquid.

5 Summary and discussion

In this paper, we explore interactions of stony planetesimals with disk-embedded primordial protoplanetary atmospheres. Modelling a stationary atmosphere requires a specification of the luminosity of the planet and a common approach is to associate luminosity with an impact energy of a stream of accreted planetesimals on the surface. Our main interest is to investigate how much gravitational energy is already depleted in the atmosphere and thus cannot contribute to a luminosity in the deep planetary atmosphere.Hence, our simulations do not provide detailed impact statistics but give an overview about the total energy dissipation of accreting bodies of different initial sizes. We not only investigate the individual behaviour of such infalling objects but also the total energy dissipation as a result of a designated planetesimal distribution in the surrounding disk.

In Sect. 3, we point out that our model lacks a detailed description of the fragmentation process of observed asteroids interacting with the Earth’s atmosphere, although it produces acceptable results with airburst heights for all tested cases. Our model does not take non-radial atmospheric trajectories nor a shape factor into account, both of which can influence the aerodynamics of the object. Keeping this in mind, the results in Sect. 3 are robust with respect to numerical and physical limitations and even investigating the extreme of the possible planetesimal parameters only moderately changes the results.

Protoplanets of different sizes accumulate different amounts of atmospheric gas during the disk phase. Primordial atmospheres are dense structures compare to secondary atmospheres (e.g. Earth’s present atmosphere) and extend up to the Hill radius. Time-dependent simulations show us that the structure of such a primordial atmosphere consists of a thin, cold outer region and a dense, hot inner layer around the protoplanet. If we consider a constant planetary luminosity that is generated by undisturbed infalling planetesimals, it is important to consider how much of the gravitational energy of the planetesimals is actually lost by friction heating as a consequence of interactions with the atmosphere and subsequently cannot contribute to core luminosity. Our results clearly show that, depending on the core mass, a significant amount of energy is released within the atmosphere. This can be as much as 100% for test bodies with diameters ranging from 1 ≤ D₀ ≤ 10⁷ cm, for a core mass with M_C > 3 M_⊕, with a primordial atmosphere at 1 Myr (as discussed in Sect. 1). One could comment that for an Earth-sized core, the maximum energy dissipation is reached at only around h = 50 km (see Fig. 7; blue line) but this certainly does not apply for larger protoplanets. Furthermore, even though energy dissipation is low within the upper atmosphere, the dynamical effects are complicated and need further investigation.

An interesting result is that parts of the atmosphere of large protoplanets can also be cooled, rather than heated, by planetesimal infall events. This clearly requires further analysis in future time-dependent models. Furthermore, partial cooling of the atmosphere could lead to contraction, leading to an accelerated inflow of surrounding disk gas and consequently to an accelerated mass increase. This scenario might be a contributing factor to why large gaseous planets already exist at early stages of disk formation. In order to get a better understanding of the dynamical effects in protoplanetary atmospheres caused by infalling planetesimals, implementing our model in the time-dependent simulations of Stökl et al. (2016) is necessary to investigate this hypothesis. Such an approach will give us the opportunity to study the time-dependent evolution of atmospheres including the back-reaction of planetesimals on the gas.

We want to emphasise that the atmospheric evolution models do not take planetary core growth into account (Stökl et al. 2016). Our model gives an overview about how much mass during various infall scenarios contributes to the atmosphere and how much reaches the surface. Additionally, we see that mass accretion for small bodies (up to D₀ = 1 m for M_C= 1 M_⊕) is a very efficient process (Fig. 4a) as those planetesimals experience drag already in the upper atmosphere and consequently slow down. This suggests that small-sized planetesimals remain trapped within the atmosphere regardless of infall speeds and infall trajectories. Considering the fact that the number density of such particles in the disk is significantly high, this mass gain should not be neglected. Additionally, we see a huge mass loss even for larger bolides (up to D₀ = 10⁵ cm) while travelling through the atmosphere for an Earth-sized protoplanet (Fig. 4c). This mass is indeed captured within the atmosphere and as we examine silicates, this mass cannot leave the atmosphere under normal circumstances. Hence, those captured silicates eventually rains out onto the core once the atmosphere is cooling down. Looking into such scenarios gets even more important when studying larger protoplanets with M_C> 1 M_⊕.

Recent core growth simulations of giant gaseous planets (Helled & Schubert 2008) suggest that planetesimal infall can drastically increase the results of the final mass of the planet due to adding a siginificant amount of refractory material. As our model not only considers the infall rate of plantesimals but also can determine at which layer the planetesimal mass is depleted, our findings can be used to improve such core growth models.

Since core growth is an important question for planetary formation, including our model in time-dependent atmospheric evolution models will eventually result in a more accurate description of primordial protoplanetary atmosphere evolution. Such a consistent time-dependent model is not only relevant during the disk phase but also afterwards as long as the planet retains a primordial atmosphere.

In this context, we also see similarities with the work of Movshovitz et al. (2010), which focusses on the growth rate of giant planets due to coagulation of dust grains within the atmosphere. These authors consider infall rates due to planetesimal accretion, which is thesource term for their coagulatiuon model. Nevertheless, they have no information about how much planetesimal mass is depleted in the various layers of the atmosphere. Our model predicts mass loss as a function of altitude, which could indeed enhance such coagulation models and give a better understanding of giant planet formation.

Fig. 7 Panel a: relative heating capacity eν in values of eν,max as a functionof altitude for 4 different core masses MC =  0.1 M⊕ (black), 1 M⊕ (blue), 3 M⊕ (green), and 5 M⊕ (red). Here, eν is the heating/cooling rate of the atmosphere evaluated with 701 mass bins, folded with the size distribution function ν and weighted with the planet’s mass accretion rate ${\mathcal{M}}_{\text{acc}}={10}^{-7}$ MC [M⊕ yr −1]. Whereas forsmaller planets planetesimals cause a general heating of the atmosphere, for larger cores local cooling can occur (compare MC = 3 and MC = 5 M⊕). Panel b: eν in values of eν,tot as a functionof optical depth τ. Again we compare 4 different core masses with MC = 0.1 M⊕ (black), 1 M⊕ (blue), 3 M⊕ (green), and 5 M⊕ (red). The dashed lines represent the individual protoplanetary surfaces. Planetesimal infall results in a maximum energy input/output within the optical thick layers of the atmosphere.

Acknowledgements

The authors acknowledge support by the FWF NFN project S116 “Pathways to Habitability: From Disks to Active Stars, Planets and Life”, and the related subprojects S11601-N16 “Hydrodynamics and Radiation in Young Star Disk Systems” and S11604-N16 “Radiation & Wind Evolution from T Tauri Phase to ZAMS and Beyond”. This publication is supported by the Austrian Science Fund (FWF). Furthermore, we want to thank the referee for clarifying the presentation of some of our results.

References

All Tables

All Figures

	Fig. 1 Density (upper panel) and temperature (lower panel) profiles of primordial atmospheres as a function of altitude, for core masses of 0.1 M⊕ (black), 1 M⊕ (green), and 3 M⊕ (red). The dashed lines represent the individual Hill radii of the cores. These data represent a snapshot at 1 Myr from the time-dependent atmospheric evolution simulations and a star with 1 M⊙ (see Stökl et al. 2016).
In the text

Fig. 2 Airburst heights of the different meteorites; plot shows the number of fragments Ni as a functionof altitude. A number of fragments, Ni = 2, indicate that both the strength criterion and distance criterion have been fulfilled for the first time. A number of fragments between 1 and 2 indicate that the strength criterion has been saisfied, but the distance criterion has not, and therefore weconsider that the body has broken apart, but the fragments have not separated far enough to be considered separate bodies.

In the text

	Fig. 3 Combined plot for various meteorites. Panel a: normalised diameter reduction while travelling through the Earth’s atmosphere,panel b: normalised velocity, and panel c: normalised mass loss as a function of altitude.
In the text

	Fig. 4 Figure comparing bolide velocity (panel a), diameter (panel b), mass loss (panel c), and temperature (panel d) as a function of altitude in the atmosphere for a planetary core mass of 1 M⊕ for various initial diameters ranging from D0 = 100 cm to D0 = 107 cm for one particular arbitrary chosen protoplanetary atmosphere.
In the text

Fig. 5 Diameter range of impacting planetesimals as a function of altitude for protoplanet core masses of 0.1 M⊕ (panel a), 1 M⊕ (panel b), and 3 M⊕ (panel c) as a function of altitude. Grey areas denote ranges between the lower and upper limit of impacting mass bins (solid lines). The plot shows bolide velocities vi with respect to altitude h with an impact velocity vh=0 > vdrift. All simulations are carried out with a grid of initial bolide diameters ranging from 1 ≤ D0 ≤ 107 cm.

In the text

	Fig. 6 Energy dissipation per mass bin in units of Emax into the atmosphere of a protoplanet with 1 M⊕. The colours represent a bolide grid ranging from D0 = 1 cm (purple) to D0 = 100 km (yellow) with NB = 701 mass bins. The total energy input is dominated by large bodies. Some bolide cool the atmosphere owing to phase transitions from solid to liquid.
In the text

Fig. 7 Panel a: relative heating capacity eν in values of eν,max as a functionof altitude for 4 different core masses MC =  0.1 M⊕ (black), 1 M⊕ (blue), 3 M⊕ (green), and 5 M⊕ (red). Here, eν is the heating/cooling rate of the atmosphere evaluated with 701 mass bins, folded with the size distribution function ν and weighted with the planet’s mass accretion rate ${\mathcal{M}}_{\text{acc}}={10}^{-7}$ MC [M⊕ yr −1]. Whereas forsmaller planets planetesimals cause a general heating of the atmosphere, for larger cores local cooling can occur (compare MC = 3 and MC = 5 M⊕). Panel b: eν in values of eν,tot as a functionof optical depth τ. Again we compare 4 different core masses with MC = 0.1 M⊕ (black), 1 M⊕ (blue), 3 M⊕ (green), and 5 M⊕ (red). The dashed lines represent the individual protoplanetary surfaces. Planetesimal infall results in a maximum energy input/output within the optical thick layers of the atmosphere.

In the text