© O. Voelkel et al. 2021

1 Introduction

The core-accretion scenario is currently the most widely used theory for planet formation. It states that at first, planetary cores form in protoplanetary disks, which then continue to grow by various forms of accretion (Pollack et al. 1996). The formation of these planetary cores clearly shapes the general picture of planet formation. To fully model the process of planet formation, we need to track the different growth processes involved, beginning from dust coagulation, pebble and dust dynamics, the formation of planetesimals, the formation of planetary embryos, and their subsequent growth until the circumstellar disk has vanished. A global model of planetesimal formation (Lenz et al. 2019) that is regulated by the local pebble flux (Birnstiel et al. 2012) was introduced into a global model of planet formation (Emsenhuber et al. 2020) in Voelkel et al. (2020). While this approach tracks the consistent formation and accretion of planetesimals on planetary embryos, the embryos themselves remain an ad hoc assumption. In this paper we investigate the formation of planetary embryos using N-body simulations (Levison et al. 2012), based on the evolution of the planetesimal surface density. Additionally, we construct an analytic one-dimensional parameterized prescription of planetary embryo formation that can be included in a global model of planet formation. In our companion paper, we add the effect of pebble accretion on the formation of planetary embryos.

Global models for planet formation that study planetary growth by solid accretion (Mordasini et al. 2012; Emsenhuber et al. 2020; Bitsch et al. 2015; and Ida & Lin 2004, to mention just a few) generally begin with the initial presence of massive objects in the circumstellar disk. These objects are mostly referred to as embryos. Once an embryo has formed, it can grow by the accretion of solids and eventually, by the accretion of gas. While the accretion of gas onto planets begins to be important at higher masses of around 10 M_⊕ (Pollack et al. 1996), these 10 M_⊕ objects in the disk are the consequence of a previous phase of solid accretion onto smaller embryos. For clarity, we define planetary embryos as objects of at least the mass of the Earth’s Moon (M = 0.0123 M_⊕). These objects are massive enough to accrete planetesimals and pebbles from their surrounding orbits, but they are far from massive enough to effectively accrete gas.

The growth of planetary embryos depends on the local disk environment, for instance, the availability of planetesimals and pebbles. These quantities change throughout the disk evolution and depend on the global evolution of the disk. Understanding where and when planetary embryos form based on the circumstellar disk evolution is of vital relevance because the evolution stage of the disk determines the subsequent growth of the embryos. While the size range from lunar-mass embryos to gas giant cores already spreads approximately 4 orders of magnitude in mass, we recall that these lunar-mass embryos themselves are the product of long-term planetesimal growth (Kokubo & Ida 1998; Kobayashi et al. 2011; Walsh & Levison 2019; Clement et al. 2020). The manner, the location, and the time at which these planetary embryos form out of much smaller planetesimals are the main subjects of this paper.

Despite all the uncertainties regarding the initial sizes at which planetesimals were formed (Schlichting et al. 2013; Schäfer et al. 2017; Walsh et al. 2017; Morbidelli et al. 2009) based on observational or theoretical arguments, here we assumed for simplicity that planetesimals all formed with a diameter of 100 km (Morbidelli et al. 2009). Even though this planetesimal size is much larger than the one inferred by other studies (Schlichting et al. 2013), we find 5 orders of magnitude in mass between a lunar-mass object and that of a 100 km planetesimal. Large planetesimals of 100 km are currently favored to explain the size distribution of asteroids and other minor bodies of the Solar System (Morbidelli et al. 2009). One hundred kilometers also seems to be the most likely size in simulations of planetesimal formation (Klahr & Schreiber 2020; Johansen et al. 2009; Abod et al. 2019), and as recent work suggests, this size is limited by diffusion (Klahr & Schreiber 2020). While 100 km planetesimals from gravitational collapse are larger than the small pebbles out of which they form, they are not massive enough to undergo pebble accretion. The formation of lunar-mass objects from 100 km planetesimals is therefore far from trivial and lays the foundation of subsequent planetary growth.

Formingmassive planetary cores of 10 M_⊕ at larger distances to the star within the lifetime of a gaseous disk is currently a challenge for planetesimal accretion models (Johansen & Bitsch 2019). A solution to this conundrum has appeared in the form of pebble accretion onto distant planetary embryos (Klahr & Bodenheimer 2006; Ormel & Klahr 2010; Lambrechts & Johansen 2012; Bitsch et al. 2015; Ndugu et al. 2017). This process describes the accretion of vastly smaller objects that radially drift toward the star, and it has been shown to bean effective planetary growth mechanism, even at larger distances up to the so-called pebble isolation mass (Lambrechts et al. 2014).

Similar to the case in which gas is accreted onto a 10 M_⊕ core, the accretion of pebbles also requires the presence of a massive body to effectively be accreted (Ormel & Klahr 2010). The previously discussed planetesimal sizes of up to 100 km are not thought to be large enough for significant pebble accretion. Assuming that pebble accretion can grow a lunar-mass object over 4 orders of magnitude to the mass of a gas giant core therefore requires the ad hoc assumption of an initial planetary embryo at a given location. This approach is commonly used in planet formation studies that form gas giants from pebble accretion, but it lacks any description of the initial solid evolution of a circumstellar disk that would form the necessary embryo. While pebble accretion requires an active radial pebble flux that is believed to decay faster than the gas disk due to radial drift (Birnstiel et al. 2012), we face a similar conundrum as before.

The question is under which circumstances a planetary embryo at a given radial distance can form within the lifetime of a radial pebble flux. To answer this question, we need a global study that models the formation of planetesimals from pebbles and tracks their following growth up to the size of lunar-mass objects. Because lunar-mass objects are not thought to form from the spontaneous collapse of a pebble cloud, numerous studies have investigated the growth from planetesimals to planetary embryos in a circumstellar disk. The timescales of planet formation from a disk of planetesimals have been estimated first by Safronov & Zvjagina (1969) and Lissauer (1987). It was later shown that the growth of planetary embryos can be split into different growth phases, such as runaway growth (Kokubo & Ida 1996) and eventually oligarchic growth (Kokubo & Ida 1998) after the embryo enhances the eccentricity of its surrounding planetesimals, effectively decreasing the accretion onto the planet. Not only is the accretion of planetesimals suppressed in the runaway regime, the embryos also arrange themselves around stable orbital separations when expressed in their mutual Hill radii (Kokubo & Ida 1998; Walsh & Levison 2019). While more has been done on the formation of planetary embryos, their growth-timescales and orbital separation are main subject of our study.

Planetary embryo formation depends on the spatial distribution of planetesimals within the circumstellar disk because these are the building blocks of planetary embryos. Models for the viscous evolution of the gas suggest a shallow planetesimal surface density (Σ_P ) profile of Σ_P ∝ r^−0.9 (Shakura & Sunyaev 1973). The minimum mass solar nebula hypothesis suggests a steeper density profile of Σ_P ∝ r^−1.5 (Weidenschilling 1977; Hayashi 1981). However, when we consider that planetesimal formation is proportional to the radial pebble flux, the surface density profile can be as steep as Σ_P ∝ r^−2.1 (Lenz et al. 2019).

The effect on planet formation of these different distributions under the assumption of initial embryo placement has recently been studied. It was suggested that the global planetesimal surface density distribution has major consequences for planet formation (Voelkel et al. 2020). Studying the formation of planetary embryos based on the planetesimal surface density slope is therefore the next logical step. Our goal is to determine the effect of the planetesimal surface density evolution on planetary embryo formation and derive an analytic recipe for planetary embryo formation. To do this, we conduct N-body simulations and model the dynamical evolution, growth, and fragmentation of planetesimals with an initial size of d = 100 km. Our study ranges from the initial gas and dust distribution to pebble and planetesimal formation up to the finally formed planetary embryos within 0.5 and 5 au of a protoplanetary disk around a solar-type star. In order to make this possible, we have connected a one-dimensional model for pebble-flux-regulated planetesimal formation (Lenz et al. 2019) with the N-body code LIPAD(Levison et al. 2012). This setup enables us to study the growth over multiple orders of magnitude in mass over 10⁶ yr with reasonable computational effort, allowing multiple simulations that cover a range of initial parameters. Based on analytic assumptions and numerical results, we present a one-dimensional model for the formation of planetary embryos as a function of the local planetesimal surface density evolution.

In the following section, we explain the physical models that we used in our study and our prescription on planetary embryo formation (Sect. 2). The connection between the one-dimensional planetesimal formation model and LIPAD can be found in Sect. 3, where we also explain the connection. Results and their discussion can be found in Sects. 4 and 5. Section 6 contains a brief summary of our study and an outlook on how we proceed with our results.

2 Planetesimal and embryo formation

Our goal is to consistently model the growth-timescales of planetary embryos from an initial disk of gas and dust. While this endeavor ranges over multiple orders of magnitude in mass, we have chosen to split it into two components. First we formed planetesimals of 100 km in diameter, using a one-dimensional parameterized description while considering pebble-flux-regulated planetesimal formation. Then we modeled the growth and fragmentation of the planetesimals in N-body simulations. Because both processes take place at the same time, it was necessary to connect our one-dimensional parameterized model with the N-body simulation, as described in Sect. 3. We focus on the description of the one-dimensional planetesimal formation and disk evolution model, as well as on the equations of the following planetesimal growth.

2.1 Disk evolution and planetesimal formation

We have chosen to use a one-dimensional viscous disk with an α prescription for turbulence (Shakura & Sunyaev 1973), to which we added a two-population model for solids (Birnstiel et al. 2012). Based on the radial drift of the solids, we formed planetesimals with a parameterized efficiency. An exact description of the two population model can be found in Birnstiel et al. (2012). We outline the basic equations and assumptions. The model uses a fixed mass relation between a smaller and a larger population of dust grains. These two populations are distinguished by whether particle growth is limited by radial drift or fragmentation, respectively. Each time step solves one advection diffusion equation of the combined solid density. $$\frac{\partial {\Sigma }_{\text{s}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left[r\left({\Sigma }_{\text{s}}\overline{u}-D_{}{\Sigma }_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{s}}}{{\Sigma }_{\text{g}}}\right)\right)\right]=0,$$(1)

with Σ_s the solid density, Σ_g the gas density, and D = D_gas∕(1 + St²) the diffusion coefficient. This contains the diffusion coefficient of the gas D_gas and the Stokes number St of the particles. ū is given as the weighted velocity of the two populations. The weighted velocity is given as $$\overline{u}=(1-f_{\text{m}}(r))\cdot u_0+f_{\text{m}}(r)\cdot u_1,$$(2)

where f_m(r) is given as themass relation described above that separates the two populations with their corresponding velocities u₀ and u₁. The individual populations are then given as $${\Sigma }_0(r)={\Sigma }_{\text{s}}(r)\cdot (1-f_{\text{m}}(r))$$(3) $${\Sigma }_1(r)={\Sigma }_{\text{s}}(r)\cdot f_{\text{m}}(r).$$(4)

The mass relation f_m was derived by fitting the two-population model to more sophisticated simulations of dust coagulation by Birnstiel et al. (2010). The values that showed the best results are given as $$f_{\text{m}}=\{\begin{array}{ll}0.97,\hfill & \text{drift-limited case}\hfill \\ 0.75,\hfill & \text{fragmentation-limited case}\text{.}\hfill \end{array}$$(5)

The decisionwhether a particle is within Σ₀ or Σ₁ was made according to its Stokes number (Birnstiel et al. 2012). Σ₀ contains particles with a small Stokes number (St ≪ 1). The motionof these particles is coupled to the motion of the gas. Σ₁ contains larger particles with St ≥ 1, which are no longer coupled to the gas. In the following we refer to Σ₀ as dust and to Σ₁ as pebbles. Planetesimals are formed based on the radial drift of the solid material in our disk. A detailed description of the planetesimal formation model can be found in Lenz et al. (2019). For our purpose, we assumed planetesimals to form with an initial size of 100 km in diameter. This choice is supported by numerical simulations of planetesimal formation by Klahr & Schreiber (2020) and observationsof asteroid and Kuiper belt objects (Morbidelli et al. 2009; Schäfer et al. 2017; Walsh et al. 2017).

The formation of planetesimals as described in Lenz et al. (2019) occurs in trapping zones in which disk instabilities can trigger planetesimal formation. These zones are distributed within the whole disk. The formation rate of planetesimals is then given proportional to the radial pebble flux and can be written as $${\dot{\Sigma }}_{\text{p}}(r)=\frac{ϵ}{d(r)}\frac{{\dot{M}}_{\text{peb}}}{2\pi r,}$$(6)

with ϵ the formation efficiency, d(r) the radial separation of pebble traps, and r the radial distance to the star. We chose d(r) to be five gas pressure scale heights and ϵ = 0.05 in our simulations. Ṁ_peb is the radial pebble flux, which in our model is defined as $${\dot{M}}_{\text{peb}}:=2\pi r{\sum }_{{\text{St}}_{\text{min}}\le \text{St}\le {\text{St}}_{\text{max}}}|v_{\text{drift}}(r,\text{St})|{\Sigma }_{\text{s}}(r,\text{St}),$$(7)

with v_drift as the radial drift velocity. The pebble-flux-regulated model for planetesimal formation results in a steeper radial planetesimal surface density profile (Σ_p ∝ r^−2.1, Lenz et al. 2019), as suggested by the minimum mass solar nebula hypothesis (Σ_P ∝ r^−1.5) or the gas surface density profile of a viscously evolving disk (Σ_P ∝ r^−0.9).

Becausewe did not specify the physical process that forms planetesimals (e.g., streaming instability or Kelvin-Helmholtz instability), this one-dimensional planetesimal formation description can be considered independent of a model. The formation of planetesimals is regulated by the local pebble flux, which in turn is regulated by dust coagulation and disk evolution (Birnstiel et al. 2012). This approach enables us to connect the timescales of the dynamical pebble evolution of the disk with the timescales of planetesimal formation. In our study, we chose to focus on three planetesimal surface density profiles while applying the formation rate from the pebble flux regulated model because it connects the viscous timescales of the disk with the formation of planetesimals.

The gas disk in our study used an initial profile as in Lynden-Bell & Pringle (1974), using α = 2 × 10⁻³, a cutoff radius of r_c = 46 au, and an initial gas surface density profile of Σ_g ∝ r^−0.9. The relative mass gain of the planetesimal disk between 0.5 and 5 au was then scaled to match a total mass of 6 M_⊕, 13 M_⊕ and 27 M_⊕ after 1 Myr, respectively. This was done in a post-processing fashion. The same post-processing was applied to the density slopes of the planetesimal surface density that is given to the N-body simulation. A more detailed description is found in Sect. 3.2.

2.2 Planetesimal growth and embryo formation

In the following we describe the one-dimensional analytical model that determines where and when lunar-mass planetary embryos are formed (based on the local evolution of the planetesimal surface density). The model connects analytic growth rates with the orbital separation of planetary embryos in the oligarchic regime. The mass of the largest object at an orbital distance r to the star at a time t is given as M_p (r, t). When planetesimals have formed at a time t₀ at a distance r, we introduce $$M_{\text{p}}(r,t_0)=M_{100\text{km}.}$$(8)

We set the initial mass to that of a planetesimal with a diameter of 100 km with a solid density of ρ_s = 1.0 g cm⁻³. During the evolution of the planetesimal disk, we integrated the mass growth rate of M_p within a swarm of planetesimals at every time step. The local mass growth rate is then given as (Lissauer 1993) $$\text{begin{align}}\frac{\text{d}{M}_p(r,t)}{\text{d}t}=\frac{\sqrt{3}}{2}{\Sigma }_{\text{P}}(r,t)\cdot \Omega (r)\pi r\text{{\_}b\text{}^}2\left(1+\frac{v_{\text{esc}}^2(M_{\text{p}},r_{\text{b}})}{v_{\infty }^2(r,t)}\right),\text{end{align}}$$(9)

with Ω as the orbital Kepler frequency, v_esc the escape velocity of an object with mass M_p and radius r_b, and v_∞ as the mean dispersion velocity within the swarm of planetesimals. We chose v_∞ (r) = e(r, t) ⋅ v_k(r) with e(r, t) as the local mean planetesimal eccentricity in our analytical model computation. v_k (r) is the Keplerian velocity at an orbital distance r. Equation (9) was integrated at every time step with the updated values for Σ_P, v_esc, and v_∞, therefore new planetesimals form over time. When M_p has reached the minimum mass of a planetary embryo M_emb (which in our study is given as a lunar mass) at a distance r, we determined this to be the location at which a lunar-mass planetary embryo can be formed. We did not track the subsequent evolution of the embryo. Our approach is solely designed to estimate the local timescales involved in forming an embryo-mass object within an evolving planetesimal disk. The eccentricity for the analytical model computation is given as e(r, t) = 5 × 10⁻⁴(1 + r^0.8), which results in a good fit to the numerical simulations. It is known that the size of planetesimals has a significant effect on the accretion rate (Fortier et al. 2007). The planetesimal size appears in v_∞, v_esc, and M_p (r, t = 0). Our model runs in Sect. 4 considered all planetesimals, including M_p (r, t = 0) to be 100 km in diameter. Equation (9) is still valid for different planetesimal sizes by adapting v_∞, v_esc, and M_p (r, t = 0), however.

Our one-dimensional embryo formation model can be described by two criteria. The first criterion refers to the necessary growth time at a distance r as a function of the planetesimal surface density evolution. The second criterion concerns the orbital separation to already present embryos. Criterion I for the embryo formation model can be written as $$M_{\text{p}}(r,t)\ge M_{\text{emb}\text{.}}$$(10)

The second criterion for the formation of a planetary embryo at r_i is the orbital separation to other planetary embryos at r_j. As suggested in numerical studies by Kokubo & Ida (1998), Kobayashi et al. (2011), Walsh & Levison (2019) and Clement et al. (2020), we find an orbital separation of planetary embryos in the oligarchic growth regime of Δr_orbit ~ 10–20R_Hill. We chose a randomized Gaussian distribution for the orbital separation of about 17 R_Hill with a standard deviation of σ_Δr = 2.5R_Hill in our analytic model runs. The mass for the computation of the Hill radius is always given as the mass of the embryos that have already been placed. Criterion II is then given as $$\Delta r_{\text{orbit}i,j}\ge \Delta r_{\text{min,}}$$(11)

where Δr_orbiti,j is the orbital distance of an embryo at r_i to an embryo at r_j. Δr_min was chosen from the Gaussian. As our analytic model does not track the subsequent dynamical evolution of a placed embryo, mergers and scattering were not accounted for. The embryos that are formed with the one-dimensional analytic model are compared to the results of the N-body simulations in Sect. 4.3.

3 LIPAD and the growth of planetesimals

3.1 LIPAD

The Lagrangian integrator for planetary accretion and dynamics (LIPAD, see Levison et al. 2012), is a particle-based (i.e., Lagrangian)code. LIPAD was developed to follow the collisional, accretional, and dynamical evolution of a large number of meter- to kilometer-sized objects throughout the entire growth process of becoming planets, making it ideal for our study. A detailed description, as well as an extensive suite of tests of LIPAD can be found in Levison et al. (2012). In addition, LIPAD has been successfully employed in previous studies of planet formation, as well as in the collisional evolution of meter- to kilometer-sized planetesimals interacting with planets or protoplanets (Kretke & Levison 2014; Levison et al. 2015; Walsh & Levison 2016, 2019; Deienno et al. 2019, 2020).

LIPAD uses the concept of tracer particles to represent a large number of small bodies with roughly the same orbit and size. Tracers are characterized by three numbers: the physical radius, the bulk density, and the constant total mass of the disk particles that it represents.

Collisional routines are employed to determine when collisions between tracers occur. In this event, following a probabilistic outcome based on a fragmentation law by Benz & Asphaug (1999), tracers can be assigned new physical radii. Therefore a distribution of tracers in LIPAD represents the size distribution of the evolving planetesimal population. The interaction among tracers results from statistical algorithms for viscous stirring, dynamical friction, and collisional damping.

Tracers that are large enough can become planetary embryos. Planetary embryos interact among themselves, as well as with tracers, through normal N-body routines (Duncan et al. 1998).

LIPAD also contains a prescription of the gaseous nebula from Hayashi (1981). This gas disk provides aerodynamic drag, eccentricity, and inclination damping on every object.

3.2 Planetesimal formation in LIPAD

We investigated different total masses of planeteimals and surface density profiles while taking their formation timescales into account. Forthis purpose, we applied the formation rate from our one-dimensional model and scaled it to the total masses after 10⁶ yr between 0.5 and 5 au. The formation of planetesimals as described in Sect. 2 scales linearly with the planetesimal formation efficiency ϵ, therefore we chose the same qualitative formation rate for our various setups. The normalized disk mass change is shown in Fig. 1. Our planetesimal formation model uses a surface density distribution to describe planetesimals in the disk, whereas LIPAD uses tracer particles. We therefore transformed our surface density into a discrete number of tracer particles. We initially defined a total number of tracer particles $N_{\text{Tracer}}\mathcal{O}(\approx {10}^4)$ to be generated in the simulation within 10⁶ yr. To obtain the mass of the individual tracers M_Tracer, we used the final mass that is in planetesimals M_Pts after 10⁶ yr: $$M_{\text{Tracer}}=\frac{M_{\text{Pts}}}{N_{\text{Tracer}}}.$$(12)

The domain of the one-dimensional surface density is split into individual rings of mass. Every 10⁴ yr we added new planetesimal tracers to the LIPAD simulation according to the formation of the planetesimal surface density ΔM_disk. Each of these newly formed tracers was assigned a heliocentric distance that was chosen randomly between the inner and outer edge of the ring in which it formed. By repeating this at every time step and in every ring, we ensured that the overall heliocentric distribution of planetesimal tracers in the LIPAD simulations matched the density slope and planetesimal distribution of the one-dimensional model. LIPAD then continues with the newly included planetesimal tracers in addition to the previously included objects that by then have grown and fragmented until the next group of tracers is included. Using this setup, we connected the timescales of pebble growth and drift, the formation of planetesimals, and their simultaneous growth. The qualitative mass change in individual setups is shown in Fig. 1. The peak of the planetesimal formation rate occurs at $T_{\Delta M_{\text{max}}}~$ 115 ky. Because the formation ofplanetesimals requires a radial pebble flux, we can reach conclusions from the planetesimal formation rate onto the remaining pebbles. About 90% of planetesimals have formed within 400 ky of our setup.

Fig. 1 Qualitative change in planetesimal disk mass Mdisk (red dots), normalized by the total disk mass after 106 yr that we use in the analytic setups. The green dots indicate the disk mass increase every 104 yr ΔMdisk, normalized by the maximum mass change ΔMmax.

4 Numerical results

In the following we present the results of nine different setups in which we varied the total mass within 0.5 to 5 au and the surface density slope with which planetesimals enter the simulation. The total masses after 10⁶ yr were 6 M_⊕, 13 M_⊕, and 27 M_⊕ and for each we varied the density slope with Σ_p ∝ r^−1.0, Σ_p ∝ r^−1.5, and Σ ∝r^−2.0, respectively. The planetesimal formation rate for these analytic setups is shown in Fig. 1. We focused on the mass and semimajor axis evolution of planetary embryos in LIPAD (Figs. 2– 4, Sect. 4.1). The embryo mass occurrences are shown in Fig. 5, Sect. 4.2. The LIPAD results are compared to the analytic model for embryo placement in Fig. 6, Sect. 4.3. Following up on this, we display the cumulative number of embryos formed (Fig. 7, Sect. 4.4), their orbital separation (Fig. 8, Sect. 4.5), and the active embryo number as well as the total mass in embryos (Fig. 9, Sect. 4.6).

4.1 Massand semimajor axis evolution

Figures 2–4 show the evolution of the N-body system within 1 Myr. We show the time and semimajor axis evolution of objects that were classified as planetary embryos in the LIPAD simulation. The classification of an embryo occurs after a tracer particle represents a single object of lunar mass. The tracer is then promoted to a planetary embryo and is treated as a single N-body object with an initial lunar mass. The subsequent growth of a given embryo is represented by the color bar and its semimajor axis evolution by a gray line. The time at which a tracer is promoted to an embryo is shown as black dots.

Because embryos can collide and eventually merge during their evolution, we distinguished active embryos and initial embryos. The number of initial embryos are the events in which tracers have been promoted to planetary embryos (number of black dots), and the number of active embryos is the number of embryos at a given time t. The red line in the plots refers to the analytic model. It indicates where M_P has surpassed a lunar mass (criterion I) when the same analytic planetesimal surface density evolution is assumed as in the N-body simulation. The red line is shown only for reference, comparing the N-body simulation with the analytical result. Even though all planetesimals that enter the LIPAD simulation have an initial semimajor axis larger than 0.5 au, we also find embryo formation within 0.5 au. This is due to dynamical interactions or scattering of the LIPAD tracer particles that lead to a nonzero planetetsimal distribution within the edge of its original formation. Because this effect is not taken into account in the analytical model density distribution, we cannot see a change in the red line within 0.5 au. This effect also has to be considered when the cumulative number of initial embryos is compared (see Fig. 7). Finally, our results show that the more massive disks (Fig. 4) form embryos earlier at close distance. Additionally, embryos in massive disks can form at larger heliocentric distances than in their less massive counterparts (Fig. 2).

4.2 Embryo mass occurrences

Figure 5 shows the number of embryo masses at $T_{M_{\text{disk}}>90\%}$(400 ky) and at 1 Myr from the simulations of Figs. 2–4. Most embryos in each simulation are found at the higher mass end of their simulation. Embryos that have low masses (≈ 0.0123 M_⊕) are less abundant than embryos that share the highest possible masses in the system. No single embryo grows substantially larger than the others in the system, in agreement with standard oligarchic growth models (Kokubo & Ida 1998).

4.3 Comparison with the analytical model

Figure 6 shows the time and location at which an object has reached the mass of a planetary embryo for both the LIPAD simulation and the analytical model from Sect. 2.2(Eqs. (9)–(11)). We show the inner edge of planetesimal formation in LIPAD at 0.5 au and give the time by which the planetesimaldisk mass has reached 90% of its total value ($T_{M_{\text{disk}}>90\%}=400$ ky). The randomization of the semimajor axis in our analytical model is given by 2.5R_Hill, as explained in Sect. 3.2. We find that the overall time and semimajor axis distribution of the analytical model agree well with the larger N-body simulations from LIPAD. The randomization of the semimajor axis reproduces the stochastic nature of the N-body process well, as does the analytic growth equation for the time it takes until embryo formation (based on the local planetesimal surface density evolution at a given distance from the star).

The innermost embryos (0.5 to 1 au) in every setup form well below $T_{M_{\text{disk}}>90\%}$, but no embryo outside 2 au forms within $T_{M_{\text{disk}}>90\%}$. The implications for possible pebble accretion from this behavior are discussed in Sect. 5.2.

4.4 Cumulative number

Figure 7 shows the cumulative number of initial embryos from the LIPAD simulation and the analytic model. The orbital separation of the embryos in the analytic model scales linearly with their Hill radius, which again scales linearly with their distance to the star. The cumulative number of planetary embryos in the analytic model therefore scales logarithmically with distance. Because the orbital separation of embryos in the N-body simulation converges to the same amount of Hill radii, we also find a logarithmic trend in the cumulative number of embryos formed in LIPAD. The total number of initial embryos is related to the total mass in planetesimals after 1 Myr. The reason for this is that embryos can form at larger distances in more massive disks.

The N-body simulations show embryo formation within 0.5 au, which is not possible in the analytic model because the planetesimal formation within 0.5 au is neglected. The innermost embryos that form in the N-body simulation are therefore due to planetesimals that moved within 0.5 au due to their dynamical interactions.

This spatial area of embryo formation is clearly defined by criterion I, see Figs. 2–4. The numberof embryos within this area can be determined using criterion II by setting their orbital separations.

Fig. 2 Time over semimajor axis evolution of the N-body simulation in LIPAD. The time and location at which an object has first reached lunar mass is indicated by the black dots in the plot. The subsequent growth of the embryo is tracked and connected with the gray lines (its mass is given by the color bar). The mass after one million years in planetesimals is 6 M⊕ in these runs. The planetesimal surface density slope is varied (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0). The red line indicates where Mprot surpasses the mass of a lunar-mass planetary embryo in the analytical model from Sect. 2.2, assuming the same evolution of the planetesimal surface density that is given to the N-body simulation.

Fig. 3 Time over semimajor axis evolution of the N-body simulation in LIPAD. The time and location at which an object has first reached lunar mass is indicated by the black dots in the plot. The subsequent growth of the embryo is tracked and connected with the gray lines (its mass is given by the color bar). The mass after one million years in planetesimals is 13 M⊕ in these runs. The planetesimal surface density slope is varied (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0 ). The red line indicates where Mprot surpasses the mass of a lunar-mass planetary embryo in the analytical model from Sect. 2.2, assuming the same evolution of the planetesimal surface density that is given to the N-body simulation.

Fig. 4 Time over semimajor axis evolution of the N-body simulation in LIPAD. The time and location at which an object has first reached lunar mass is indicated by the black dots in the plot. The subsequent growth of the embryo is tracked and connected with the gray lines (its mass is given by the color bar). The mass after one million years in planetesimals is 27 M⊕ in these runs. The planetesimal surface density slope is varied (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0 ). The red line indicates where Mprot surpasses the mass of a lunar-mass planetary embryo in the analytical model from Sect. 2.2, assuming the same evolution of the planetesimal surface density that is given to the N-body simulation.

Fig. 5 Number of embryos within a given mass bin at T = 400 ky ($T_{M_{\text{disk}}>90\%}$) and T = 1 Myr from the simulations of Figs. 2–4.

4.5 Orbital separation

Figure 8 shows the time evolution of the average orbital separation of initial embryos for the LIPAD simulation and the analytical model, see Fig. 6. The mean orbital separation of all systems converges to a value around 10–15 R_Hill after 200–400 ky. The orbital separation is a free parameter from criterion II of the analytical model that we chose to fit the numerical results from our N-body simulations. In combination with criterion I, this allows us to predict the number, spatial distribution, and formation time of planetary embryos for a specific planetesimal surface density evolution. The total number of embryos is given as the number of orbital separations (criterion II) within the possible area of embryo formation (criterion I). Their spatial distribution is determined by their orbital separation, which is a function of the mutual Hill radii. In this way, the absolute orbital separation between embryos increases linearly with increasing distance to the star, leading to a logarithmic cumulative number of initial embryos (see Fig. 7).

Because of the low number of embryos for early times, the mutual distance can differ strongly between the analytical model and the LIPAD runs. This behavior however would also occur if we were to attempt to compare two LIPAD runs with similar initial conditions. This is due to the chaotic nature of the N-body evolution. We can show that for a larger number of embryos, the orbital separation in the analytical model shows the same behavior as in the N-body simulations.

Fig. 6 Analytical model for embryo formation and embryo formation in LIPAD. The black dots indicate the time and location at which an object reached planetary embryo mass in the LIPAD simulations. The red dots indicate the time for each distance from the star at which a planetary embryo is placed using ouranalytical model. The orbital separation as input to the analytical model is given as 17 RHill, with a randomization of 2.5 RHill The inner edge of planetesimal formation in the LIPAD runs was chosen to be at 0.5 au for numerical performance. We vary the total mass in planetesimals after 106 yr from 0.5 to 5 au (6 M⊕,13 M⊕, and 27 M⊕) and the planetesimal surface density slope (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0).

4.6 Active number of embryos and mass in embryos

Figure 9 shows the number of active embryos over time, the total mass that is given in these embryos, and the fraction of the total mass in embryos after 1 Myr (M_Emb) over the total mass in the system after 1 Myr (M_D). The numberof active embryos after 1 Myr is between 30 and 40 embryos for eight out of our nine runs. Only the 6 M_⊕ and Σ_P ∝ r^−1.0 run contains fewer embryos (N_active = 22) after 1 Myr. While the total number of active embryos seems insensitive to the total planetesimal mass or the planetesimal surface density profile, the same is not true for the total mass that is in planetary embryos after one million years.

The total mass in embryos increases for steeper planetesimal surface density profiles and higher total masses after 1 Myr. The fraction of mass M_Emb∕M_Disk that is transformed into embryos increases for both higher masses and the slope of the planetesimal surface density. The numberof embryos does not simply increase for more massive planetesimal disks in our runs because the embryos that form grow larger in more massive disks. They thereby increase their orbital separation to that of the other embryos again. While higher planetesimal disk masses allow for a larger zone in which embryo formation is possible (criterion I), the present embryos increase their orbital spacing due to their higher masses as well. In the case of 27 M_⊕ in planetesimals after 1 Myr, the number of embryos decreases slightly (N_active = 38 for Σ_P ∝ r^−1.0, N_active = 36 for Σ_P ∝ r^−1.5, N_active = 32 for Σ_P ∝ r^−2.0) for the steeper planetesimal surface density profiles, but their mass increases drastically (M_Emb ≈ 5 M_⊕ for Σ_P ∝ r^−1.0, M_Emb ≈ 10 M_⊕ for Σ_P ∝ r^−1.5, M_Emb ≈ 13 M_⊕ for Σ_P ∝ r^−2.0).

Fig. 7 Cumulative number of embryos formed during the LIPAD runs from Figs. 2–4 (black dots). The red dots show the cumulative number of embryos that would be placed according to the analytical model from Sect. 2.2.

5 Discussion

5.1 Embryo formation - LIPAD

Figures 2–4 clearly show that embryo formation for every power-law planetesimal surface density profile occurs from the inside out. This is an expected result because the growth-timescales in the inner disk are shorter and the densities in planetesimals are correspondingly higher. Even though the individual moment and location at which an embryo forms (black dots) appears to be stochastic, there is a pattern to be found in the embryo formation of the system. The red curve that marks criterion I is well within the area of the initial embryos. The individual locations of the embryos, even though they follow the trend of the red line, appear chaotic. The exact location and time at which an object reaches the size of a planetary embryo appears stochastic due to the stochastic behavior of the N-body, but the analytic growth equations constrain the zone of their individual formation well.

Anothereffect that can be found is that embryos increase their orbital distance to other embryos when they grow in mass. This effect has previously been found and discussed by Kokubo & Ida (1998) and Kobayashi et al. (2011) when they studied the oligarchic growth of massive objects. In the general picture, initial embryos begin to form earliest at closer distance to the star.

Furthermore, the orbital separation of planetary embryos when expressed in terms of their Hill radii converges to a similar value in every setup we studied, as shown in Fig. 8. This directly results in a cumulative number of embryos that scales logarithmic with distance, as shown in Fig. 7. The stochastic behavior of the orbital separation in Fig. 8 stems from the randomization, as explained in Sect. 2.2. With the stochasticity in our approach, the cumulative number of embryos in Fig. 7 would vary for the analytic and the N-body run. Overall, however, they agree very well. Additionally, the more embryos form, the more reliable the statistics. We see this as the reason for the better agreement within the setups that result in a higher number of embryos. As a comparison we show the cumulative number of embryos that would form with the analytic model, in which the orbital separation is always expressed in terms of the Hill radius of the previously placed embryos. The cumulative number of embryos and the number of active embryos does not vary sensitively with the initial parameters (total mass and planetesimal surface density slope). The total mass that is converted into embryos, however, does depend strongly on the planetesimal surface density slope and the total mass in planetesimals, as shown in Fig. 9. The number of active embryos even decreases slightly for higher disk masses and steeper planetesimal surface density profiles. As Figs. 2–4 show, the area in which planetary embryos form becomes larger for higher masses and steeper density profiles. Because their orbital separation increases for higher masses and because the mean orbital distance converges to the same number of Hill radii (Fig. 8), the total number of embryos within 1 Myr is not sensitive to different input parameters either.

Fig. 8 Mean orbital separation of the initial embryos from the LIPAD runs and the analytical model embryos over time. The input parameters for the analytical model are given as 17 RHill with a randomization of 2.5 RHill.

5.2 Implications for pebble accretion

While the effect of pebble accretion on the formation of planetary embryos will be the main subject of our companion paper, we can already discuss some viable constraints here. It is notable to mention that the formation timescale of planetesimals is well within the formation timescales of the planetary embryos. This states that the formation of planetesimals continues to occur after planetary embryos have already formed from previously formed planetesimals. Because the growth rate of planetary embryos depends linearly on the local planetesimal surface density (Eq. (9)), the local formation of planetesimals has to be taken into account to model the growth-timescales consistently. We conclude that because the formation of planetesimals requires a radial pebble flux, we can estimate first constraints on this pebble flux using our setup, and subsequently, we can estimate first constraints on the possibility of continuous pebble accretion.

Even though we did not take the accretion of pebbles onto planetesimals or planetary embryos into account in our simulations, we wish to highlight their importance in the general context of planetary growth, as has been displayed by several studies such as Ormel & Klahr (2010), Bitsch et al. (2015), and Ndugu et al. (2017). The efficiency of pebble accretion directly depends on the local pebble flux at the location of an accreting body of sufficient mass. Becausee the formation of planetesimals ΔM_disk scales linearly with the local pebble flux, we can also derive from Fig. 1 that the pebble flux decreases drastically within the first 10⁶ yr of the systems evolution. However, because we continue to form planetsimals well after the first embryos have formed, these embryos could grow by the remaining pebble flux that continues to form planetesimals as well. This indicates that the growth time scales of planetary embryos is a determining factor in defining the global efficiency of pebble accretion. A certain embryo size has to be reached at first to effectively accrete pebbles.

Another crucial effect on the pebble flux evolution is the formation of planetesimals itself because they form based on the disk evolution. The more planetesimals form, the earlier we also form planetary embryos that could accrete pebbles. However, the more planetesimals are formed, the lower the pebble flux because mass is transferred into planetesimals. Even though the exact evolution of the pebble flux differs for every disk, the results of our study can already be used to apply first constraints on the magnitude of the pebble flux based on the formation timescales of planetary embryos. Because ~ 90% of our planetesimals form within 400 ky of our setup, we conclude that the magnitude of the pebble flux has decreased significantly before that time. Embryos that form after 400 ky would therefore not be able to undergo significant pebble accretion in our model.

In our setup most embryos that form within 400 ky also form within 1 au. The formation of farther out embryos around 1.5–2.0 au occurs well after 400 ky. In conclusion, it is not possible for farther out planetary embryos to undergo pebble accretion in our setup. This statement holds true when we assume a power-law distribution for the planetesimal surface density such as the minimum mass solar nebula hypothesis, the dust profile of a viscous disk, or the pebble-flux-regulated planetesimal formation surface density profile. As a first estimate on the timescales necessary for an embryo to form at larger distances, we extrapolated our analytic model for 10 Myr. The times for an embryo to form at 2, 3, and 4 au are shown in Table 1.

Fig. 9 Active number of planetary embryos and mass in planetary embryos over time for the LIPAD runs from Figs. 2–4. We also show the fraction (MEmb∕MD) of mass in embryos over the final planetesimal mass that entered the disk after 1 Myr.

5.3 Architecture of planetary systems

Following our findings from Sect. 5.2, it is not too far-fetched to state that the architecture of planetary systems might very well be determined within the first few 100 ky of their formation in terms of pebble accretion. Our study assumed power-law density profiles for the planetesimal surface density, and our results of inside-out planetary embryo formation are a direct consequence of this. If we were to assume deviation from the power-law profile due to local substructures in the disk, for instance, around the ice line (Drążkowska & Alibert 2017), this picture might change.

The early formation of planetary embryos around the ice line might lead to the formation of cold giant planets through pebble accretion. The formation of these planets can then have major consequences for the subsequent evolution of the inner system. Assuming that outer planetary embryos form early enough to undergo significant pebble accretion, they could alter the evolution of the inner system drastically because they would reduce the pebble flux that reaches the terrestrial planet region. The additional planetesimal formation itself will also have strong consequences for the interior pebble flux. An early decrease in the pebble flux would also lead to a decrease in the formation of planetesimals in the terrestrial planet region. This would again affect the formation of planetary embryos and planetary growth. It becomes clear that the formation of planetesimals, the formation of planetary embryos, and the evolution of the pebble flux are tightly connected within the first few 100 ky of a circumstellar disk.

Anotherscenario that might change the evolution of the system would be the stochastic formation of a planetesimal with an initial size much larger than 100 km (Johansen et al. 2007). The formation of a significantly larger planetesimal in a reservoir of 100 km planetesimals and pebbles could reduce the timescales of planetary embryo formation significantly. This could lead to the presence of planetaryembryos at much larger distances within the lifetime of the pebble flux.

5.4 Embryo formation: analytic model

The one-dimensional analytic parameterized approach agrees well with the sophisticated N-body simulations in terms of the formation timescales of a lunar-mass object and the total number of objects that reach this given size. In two out of our nine runs, the deviation of the total number of embryos is below 5%, in four out of nine runs it is below 10%, and in eight out of nine runs it is below 25%. Only the 6 M_⊕, Σ_P ∝ r^−1.0 run deviates more strongly (≈40%). The Hill criterion for the orbital separation of planetary embryos completely determines the number of planetary embryos without additional assumptions. Considering the time and location at which an object reaches the mass of a planetary embryo, we show that the analytic prescription represents the analytic planetesimal surface density evolution well (Sect. 4.3).

The N-body simulations require weeks (sometimes months) of computation time with the same planetesimal input, whereas the parameterized model takes mere seconds. While the N-body simulations clearly involve more complexity that allows for a more complete picture of the problem, the question of the manner, location and time at which many initial planetary embryos form is well reproduced with the analytic model. This makes the analytic approach well suited for other studies that aim for statistical properties in which computational time is a limiting factor, such as planet population synthesis.

Even though our study focused on an area from 0.5 to 5 au, the analytical model should also hold true at farther locations and might be a valuable asset in considering planetary embryo formation in far out ring-like structures of circumstellar disks, as seen in ALMA observations. Other studies regarding planet formation through pebble accretion may use our findings to modify their initial conditions in terms of the available pebble flux, as explained in greater detail in Sect. 5.2.

6 Summary and outlook

We studied the spatial distribution and formation timescales of planetary embryos from an initial disk of gas and dust. For this purpose, we coupled a one-dimensional model for viscous disk evolution and planetesimal formation to the LIPAD code that studies the dynamical N-body evolution of the evolving planetesimal system. The size of an initial planetesimal was given as 100 km in diameter, and it grew dynamically through collisions with other planetesimals. We analyzed the first million years of nine different systems in which we varied the total mass in planetesimals and their surface density profile.

In combination with analytic estimates of growth rates of planetesimals based on their local surface density, we derived an analytical model for planetary embryo formation. Our model reproduces the spatial distribution and formation time of planetary embryos well. We used their orbital separation as a free parameter that can be fit to match the N-body simulations. The model can be used in further studies (e.g., global models of planet formation and population synthesis) that use a planetesimal surface density description to consistently model the spatial distribution and formation time of planetary embryos. The main findings of planetary embryo formation based on pebble flux regulated planetesimal formation are listed below:

• embryos form first in the innermost regions of planetesimal formation due to shorter growth time scales close to the star and higher planetesimal surface densities;

• the innermost embryos (<1 au) form well within the presence of an active pebble flux for most planetesimal disks, whereas the outer embryos (>2 au) fail to do so in any disk we studied;

• higher planetesimal disk masses or steeper planetesimal surface density profiles do not result in a higher number of active embryos, but in more massive embryos within a larger area.

We linked the formation timescale of planetesimal formation and the evolution of the radial pebble flux to the formation timescale of lunar-mass objects that formed by planetesimal collisions. In doing so, we determined crucial constraints for the possibility of pebble accretion as a planet formation process. These constraints need to be considered in studies that involve pebble accretion on planetary embryos because we showed that a planetary embryo and a pebble flux strongly depend on the radial distance of the embryo to the star.

It is shown that a power-law planetesimal surface density profile cannot build planetary embryos at larger distances within the timescale of a radial pebble flux. This consequence arises from the interplay of pebble-flux-regulated planetesimal formation and the timescales involved to form planetary embryos from 100 km sized bodies. The more planetesimals are formed, the earlier planetary embryos are formed, but the more planetesimals are formed, the lower the mass that remains in pebbles. Vice versa, when we decrease the formation of planetesimals to maintain a higher pebble flux, the growth time scales for planetary embryos increase as a result of lower planetesimal surface densities.

Future studies will include disk-consistent pebble accretion in the N-body simulation to study the effect of an active pebble flux on the formation of planetary embryos. Another study that will follow the approach we presented here will study the formation of embryos in far-out planetesimal rings that could result from pressure bumps during the disks evolution.

Acknowledgments

We thank Gonzalo Carlos de Elía for his very constructive and careful review that helped improving this work. R.D. and K.K. acknowledge support from NASA’s Emerging Worlds program. We also wish to thank Kevin Walsh, Hal Levison and Christian Lenz for many fruitful discussions on this project. This research has been supported by the Deutsche Forschungsgemeinschaft Schwerpunktprogramm (DFG SPP) SPP 1992 “Exploring the diversity of extrasolar planets” contract : KL 1469/17-1 Consistent Planetesimal Formation from Pebbles for Synthetic Population Syntheses of Exo-Planets.

References

All Tables

All Figures

	Fig. 1 Qualitative change in planetesimal disk mass Mdisk (red dots), normalized by the total disk mass after 106 yr that we use in the analytic setups. The green dots indicate the disk mass increase every 104 yr ΔMdisk, normalized by the maximum mass change ΔMmax.
In the text

Fig. 2 Time over semimajor axis evolution of the N-body simulation in LIPAD. The time and location at which an object has first reached lunar mass is indicated by the black dots in the plot. The subsequent growth of the embryo is tracked and connected with the gray lines (its mass is given by the color bar). The mass after one million years in planetesimals is 6 M⊕ in these runs. The planetesimal surface density slope is varied (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0). The red line indicates where Mprot surpasses the mass of a lunar-mass planetary embryo in the analytical model from Sect. 2.2, assuming the same evolution of the planetesimal surface density that is given to the N-body simulation.

In the text

Fig. 3 Time over semimajor axis evolution of the N-body simulation in LIPAD. The time and location at which an object has first reached lunar mass is indicated by the black dots in the plot. The subsequent growth of the embryo is tracked and connected with the gray lines (its mass is given by the color bar). The mass after one million years in planetesimals is 13 M⊕ in these runs. The planetesimal surface density slope is varied (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0 ). The red line indicates where Mprot surpasses the mass of a lunar-mass planetary embryo in the analytical model from Sect. 2.2, assuming the same evolution of the planetesimal surface density that is given to the N-body simulation.

In the text

Fig. 4 Time over semimajor axis evolution of the N-body simulation in LIPAD. The time and location at which an object has first reached lunar mass is indicated by the black dots in the plot. The subsequent growth of the embryo is tracked and connected with the gray lines (its mass is given by the color bar). The mass after one million years in planetesimals is 27 M⊕ in these runs. The planetesimal surface density slope is varied (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0 ). The red line indicates where Mprot surpasses the mass of a lunar-mass planetary embryo in the analytical model from Sect. 2.2, assuming the same evolution of the planetesimal surface density that is given to the N-body simulation.

In the text

	Fig. 5 Number of embryos within a given mass bin at T = 400 ky ($T_{M_{\text{disk}}>90\%}$) and T = 1 Myr from the simulations of Figs. 2–4.
In the text

Fig. 6 Analytical model for embryo formation and embryo formation in LIPAD. The black dots indicate the time and location at which an object reached planetary embryo mass in the LIPAD simulations. The red dots indicate the time for each distance from the star at which a planetary embryo is placed using ouranalytical model. The orbital separation as input to the analytical model is given as 17 RHill, with a randomization of 2.5 RHill The inner edge of planetesimal formation in the LIPAD runs was chosen to be at 0.5 au for numerical performance. We vary the total mass in planetesimals after 106 yr from 0.5 to 5 au (6 M⊕,13 M⊕, and 27 M⊕) and the planetesimal surface density slope (ΣP ∝ r−1.0, ΣP ∝ r−1.5, and ΣP∝ r−2.0).

In the text

	Fig. 7 Cumulative number of embryos formed during the LIPAD runs from Figs. 2–4 (black dots). The red dots show the cumulative number of embryos that would be placed according to the analytical model from Sect. 2.2.
In the text

	Fig. 8 Mean orbital separation of the initial embryos from the LIPAD runs and the analytical model embryos over time. The input parameters for the analytical model are given as 17 RHill with a randomization of 2.5 RHill.
In the text

	Fig. 9 Active number of planetary embryos and mass in planetary embryos over time for the LIPAD runs from Figs. 2–4. We also show the fraction (MEmb∕MD) of mass in embryos over the final planetesimal mass that entered the disk after 1 Myr.
In the text