Linking planetary embryo formation to planetesimal formation I: The impact of the planetesimal surface density in the terrestrial planet zone

The growth time scales of planetary embryos and their formation process are imperative for our understanding on how planetary systems form and develop. They determine the subsequent growth mechanisms during the life stages of a circumstellar disk. We quantify the timescales and spatial distribution of planetary embryos via collisional growth and fragmentation of dynamically forming 100km sized planetesimals. In our study, the formation timescales of viscous disk evolution and planetesimal formation are linked to the formation of planetary embryos in the terrestrial planet zone. We connect a one dimensional model for viscous gas evolution, dust and pebble dynamics and pebble flux regulated planetesimal formation to the N-body code LIPAD. Our framework enables us to study the formation, growth, fragmentation and evolution of planetesimals with an initial size of 100km in diameter for the first million years of a viscous disk. Our study shows the effect of the planetesimal surface density evolution on the preferential location and timescales of planetary embryo formation. A one dimensional analytically derived model for embryo formation based on the local planetesimal surface density evolution is presented. This model manages to reproduce the spatial distribution, formation rate and total number of planetary embryos at a fraction of the computational cost of the N-body simulations. The formation of planetary embryos in the terrestrial planet zone occurs simultaneously to the formation of planetesimals. The local planetesimal surface density evolution and the orbital spacing of planetary embryos in the oligarchic regime serve well as constraints to model planetary embryo formation analytically. Our embryo formation model will be a valuable asset in future studies regarding planet formation.


Physical motivation
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). It goes without saying that the formation of these planetary cores shapes the general picture of planet formation. To fully model the process of planet formation, one needs 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 to a global model of planet formation  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 di-mensional, parameterized prescription of planetary embryo formation, that can be included into a global model of planet formation. In our companion paper, we add 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), Ida & Lin (2004) to mention just a few) generally begin with the initial presence of massive objects in the circumstellar disk, that are mostly referred to as embryos. Once an embryo has formed, it can grow by the accretion of solids and eventually the accretion of gas. While the accretion of gas on planets begins to be important at larger masses of around 10 M ⊕ (Pollack et al. 1996), the presence of these 10 M ⊕ objects in the disk is the consequence of a previous phase of solid accretion on smaller embryos. For clarity, we define planetary embryos as objects of at least the mass of the earths moon (M = 0.0123 M ⊕ ). These objects are massive enough to accrete planetesimals and pebbles from their surrounding orbits, but far from massive enough to effectively accrete gas. The growth of planetary embryos depends on the local disk environment, like e.g. the availability of planetesimals and pebbles. These quantities change over the course of the disks evolution and depend on the global evolution of the disk. Understanding A&A proofs: manuscript no. main where and when planetary embryos form, based on the circumstellar disks evolution is of vital relevance, as 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 over roughly 4 orders of magnitude in mass, one has to keep in mind that these lunar mass embryos themselves are the product of long term planetesimal growth (Kokubo & Ida 1998;Kobayashi et al. 2011;. How, where and when these planetary embryos form out of much smaller planetesimals will be the main subject of this paper. Despite all the uncertainties regarding the initial sizes that planetesimals were formed (Schlichting et al. 2013;Schäfer et al. 2017;Walsh et al. 2017;Morbidelli et al. 2009); based in observational or theoretical arguments), here we will for simplicity assume planetesimals all formed with a diameter of 100 km (Morbidelli et al. (2009)). Even though this planetesimal size is much larger than that 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). 100 km 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 large in comparison to 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. Forming massive planetary cores of 10 M ⊕ at larger distances to the star within the lifetime of a gaseous disk is currently challenging planetesimal accretion models (Johansen & Bitsch 2019). A solution to this conundrum has appeared in the form of pebble accretion on 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 towards the star and is shown to be an 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 believed 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 on 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 presence of the gas disk due to radial drift (Birnstiel et al. 2012), we face a similar conundrum as before. Under which circumstances can a planetary embryo at a given radial distance form within the lifetime of a radial pebble flux? To answer this question, one needs a global study that models the formation of planetesimals from pebbles and track their following growth up to the size of lunar mass objects.

Previous work
Since lunar mass objects are not believed to form from the spontaneous collapse of a pebble cloud, there have been numerous studies that investigate the growth from planetesimals to planetary embryos in a circumstellar disk. Estimating the timescales of planet formation from a disk of planetesimals go back to Safronov & Zvjagina (1969) and Lissauer (1987). Following up, it was shown that the growth of planetary embryos can be split in different growth phases, like runaway growth (Kokubo & Ida 1996) and eventually oligarchic (Kokubo & Ida 1998) once the embryo enhances the eccentricity of its surrounding planetesimals, effectively decreasing the accretion on the planet. Not only is the accretion of planetesimals suppressed in the runaway regime, but also arrange the embryos themselves around stable orbital separations when expressed in their mutual Hill radii (Kokubo & Ida (1998), ). While more has been done on the formation of planetary embryos, their growth timescales and orbital separation will be main subject to our study. Planetary embryo formation depends on the spatial distribution of planetesimals within the circumstellar disk, as they 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, if considering 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 and suggests that the global planetesimal surface density distribution has major consequences for planet formation (Voelkel et al. 2020). Therefore, studying the formation of planetary embryos based on the planetesimal surface density slope is the next logical step.

The goal of this study
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. For that, we conduct N-body simulations and model the dynamical evolution, growth and fragmentation of planetesimals with an initial size of d = 100km. Our study ranges from the initial gas and dust distribution, over pebble and planetesimal formation up to the finally formed planetary embryos within 0.5au and 5au 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 6 years at a 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 will explain the physical models that we use in our study and our prescription on planetary embryo formation. (Sect. 2). The connection between the one dimensional planetesmal formation model and LIPAD, as well as their explanation can be found in Sect. 3. Results and their dis-cussion can be found in Sect. 4 and Sect. 5. Sect. 6 contains a brief summary of our study and an outlook on how to proceed with the obtained results.

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

Disk evolution and planetesimal formation
We have chosen to use a one dimensional viscous disk with an α prescription for turbulence (Shakura & Sunyaev 1973), in which we have added a two population model for solids (Birnstiel et al. 2012). Based on the radial drift of the solids we form planetesimals with a parameterized efficiency. An exact description of the two population model can be found in Birnstiel et al. (2012). Here we will 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 on whether particle growth is limited by radial drift or fragmentation respectively. Each time step solves one advection diffusion equation of the combined solid density with Σ s the solid density, Σ g the gas density, D g the diffusion coefficient andū the weighted velocity of the two populations. The weighted velocity is given as where f m (r) is given as the aforementioned mass relation that separates the two populations with their corresponding velocities u 0 and u 1 . The individual populations are then given as: 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 m = 0.97, drift limted case 0.75, fragmentation limited case The decision on whether a particle is within Σ 0 or Σ 1 is done by its Stokes number (Birnstiel et al. 2012). Σ 0 contains particles with a small Stokes number (S t 1). The motion of these particles is coupled to the motion of the gas. Σ 1 contains larger particles with S t ≥ 1, which are no longer coupled to the gas. In the following we will refer to Σ 0 as dust and Σ 1 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 assume 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 observations of 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 aṡ with the formation efficiency, d(r) the radial separation of pebble traps and r the radial distance to the star. We have chosen d(r) to be 5 gas pressure scale heights and = 0.05 in our simulations.Ṁ peb is the radial pebble flux, which in our model is defined aṡ 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 ). Due to the fact we do not specify the physical process that form planetesimals (e.g., streaming instability, Kelvin Helmholts, etc.), this one dimensional planetesimal formation description can be considered model independent. The formation of planetesimals is regulated by the local pebble flux. The latter 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, as it connects the viscous timescales of the disk with the formation of planetesimals.

Planetesimal growth and embryo formation
In the following we will describe the one dimensional analytical model that determines where and when lunar mass planetary embryos will be formed (based on the local planetesimal surface density evolution). The model connects analytic growth rates with the orbital seperation 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). Once planetesimals have fomred at a time t 0 at a distance r, we introduce We set the initial mass to that of a 100 km in diameter planetesimal with a solid density of ρ s = 1.0 g/cm 3 . During the evolution of the planetesimal disk, we integrate the mass growth rate of M p within a swarm of planetesimals in every timestep. The local mass growth rate is then given as (Lissauer 1993) with Ω as the orbital keppler frequency, v esc the escape velocity of an object with mass M p and v ∞ the mean dispersion velocity within the swarm of planetesimals. We choose 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 kepplerian velocity at an orbital distance r. Eq. 9 is integrated in every timestep with the updated values for Σ P , v esc and v ∞ , hence new planetesimals form over time. Once 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 determine this to be the location at which a lunar mass planetary embryo can be formed. We do not track the subsequent evolution of the embryo. Our approach is solely designed to estimate the local timescales involved to form an embryo mass object within an evolving planetesimal disk. The eccentricity for the analytical model computation is given as e(r, t) = 5 · 10 −4 (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 considers all planetesimals, including M p (r, t = 0) to be 100 km in diameter. Eq. 9 however is still valid for different planetesimal sizes, by adapting v ∞ , v esc and M p (r, t = 0). 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 seperation to already present embryos. Criterion I for the embryo formation model can be written as: 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 by numerical studies by Kokubo & Ida (1998) (2019) we find an orbital separation of planetary embryos in the oligarchic growth regime of ∆r orbit ∼10-20R Hill . We choose a randomized Gaussian distribution for the orbital separation around 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: where ∆r orbiti, j is the orbital distance of an embryo at r i to an embryo at r j . ∆r min is chosen from the Gaussian. 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.

LIPAD and the growth of planetesimals
LIPAD (Lagrangian Integrator for Planetary Accretion and Dynamics; 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 meterto kilometer-sized objects through the entire growth process to become 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 succesfuly employed in previous studies of planet formation, as well as collisional evolution of meter-to kilometer-sized planetesimals interacting with planet/protoplanets (Kretke & Levison 2014;Levison et al. 2015;Walsh & Levison 2016Deienno et al. 2019Deienno et al. , 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 will represent.
Collisional routines are employed to determine when collisions between tracers will happen. 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 will represent the size distribution of the evolving planetesimal population. The interaction among tracers is resultant from statistical algorithms for viscous stirring, dynamical friction, and collisional damping.
Large enough tracers can be promoted to become planetary embryos. Planetary embryos interact among themselves, as well as with Tracers, via normal N-body routines (Duncan et al. 1998).
LIPAD also has a prescription of the gaseous nebula from Hayashi (1981). This gas disk provides aerodynamic drag, eccentricity, and inclination damping on every object.

Planetesimal formation in LIPAD
We aim to investigate different total masses of planeteimals and surface density profiles, while taking their formation timescales into account. For that purpose we apply the formation rate from our one dimensional model and scale it to the total masses after 10 6 years between 0.5 au and 5 au. The formation of planetesimals as described in Sect. 2 scales linearly with the planetesimal formation efficiency , which is why we choose the same qualitative formation rate for our various setups. The normalized disk mass change can be seen in Fig. 1. Our planetesimal formation model uses a surface density distribution to describe planetesimals in the disk, whereas LIPAD uses tracer particles. For that matter we transform our surface density into a discrete number of tracer particles. We initially define a total number of tracer particles N T racer O(≈ 10 4 ) to be generated in the simulation within 10 6 years. To get the mass of the individual tracers M T racer we use the final mass that is in planetesimals M Pts after 10 6 years: The domain of the one dimensional surface density is split into individual rings of mass. Every 10 4 years we add new planetesimal tracers to the LIPAD simulation according to the formation of the planetesimal surface density ∆M disk . Each of those newly formed tracers is assigned a heliocentric distance that is chosen randomly between the inner and outer edge of the ring in which it formed. Doing this in every time step and every ring, we ensure that the overall heliocentric distribution of planetesimal tracers in the LIPAD simulations will match the density slope and planetesimal distribution of the one dimensional model. LIPAD then continues with the newly included planetesimal tracers additional to the previously included objects that by then have grown and fragmented until the next group of tracers is included. Using this setup we connect the timescales of pebble growth and drift, the formation of planetesimals and their simultaneous growth. The qualitative mass change of the individual setups can be seen in Fig. 1

Mdisk/Mtot
Mdisk/ Mmax T Mmax 115ky TM disk > 90% 400ky Fig. 1: Qualitative change of the planetesimal disk mass M disk (red dots), normalized by the total disk mass after 10 6 years that we use in the analytic setups. The green dots indicate the disk mass increase every 10 4 years ∆M disk , normalized by the maximum mass change ∆M max .

Numerical results
In the following we will present the results of nine different setups in which we vary the total mass within 0.5 au to 5 au and the surface density slope with which planetesimals enter the simulation. The total masses after 10 6 years are 6 M ⊕ , 13 M ⊕ and 27 M ⊕ and for each we vary the density slope with Σ p ∝ r −1.0 , Σ p ∝ r −1.5 and Σ ∝ r −2.0 respectively. The planeteismal formation rate for these analytic setups is shown in Fig. 1. We focus on the mass and semimajor axis evolution of planetary embryos in LIPAD ( Fig. 2 -Fig. 4 Fig. 4 show the evolution of the N-body system within 1 Myrs. 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 grey line. The occurrence of when a tracer is promoted to an embryo is shown as black dots. Since embryos can collide and eventually merge during their evolution we make a distinction between 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 assuming the same analytic planetesimal surface density evolution 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 en-ter the LIPAD simulation have an initial semimajor axis of above 0.5 au, we find embryo formation within 0.5 au as well. This is due to dynamical interactions/scattering of the LIPAD tracer particles that lead to a nonzero planetetsimal distribution wihtin the edge of its original formation. Since this effect is not taken into account in the analytical model density distribution, we cannot see a change of the red line within 0.5 au. This effect also has to be considered when comparing the cumulative number of initial embryos (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).  Fig. 2 -Fig. 4. Most embryos in each simulation are found in 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. There is no single embryo growing substantially larger than the others in the system, in agreement with standard oligarchic growth models (Kokubo & Ida 1998). Fig. 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 (Eq. 9, 10, 11). The red dots refer to analytical model with the same planetesimal surface density increase as in the LIPAD simulations, whereas the black dots are those from the N-body simulation shown in Fig. 2 -Fig.  4. We show the inner edge of planetesimal formation in LIPAD at 0.5 au and give the time by which the planetesmal disk mass has reached 90% of its total value (T M disk >90% = 400ky). The randomization of the semimajor axis in our analytical model is given by 2.5R Hill as explained in Sec. 3.1. We find that the overall time and semimajor axis distribution of the analytical model is in well agreement with the larger N-body simulations from LIPAD. The randomization of the semimajor axis does well in reproducing the stochastic nature of the N-body process, as well as the analytic growth equation does for the time it takes until embryo formation (based on the local planetesimal surface density evolution at a given distance from the star) is possible. The innermost embryos (0.5 au to 1 au) in every setup form well below T M disk >90% , but no embryo outside 2 au forms within T M disk >90% . The implications on possible pebble accretion from this behavior are discussed in Sect. 5.2. Fig. 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 logarithmic with distance. Since 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 Myrs. The reason for this is that in more massive disks, embryos can form at larger distances. The N-body simulations show embryo formation within 0.5 au, which is not possible in the analytic model, since 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 dynamcial interactions. This spatial area of embryo formation is well defined by Criterion I, see Fig. 2 -Fig. 4. The number of embryos within this area can be determined using Criterion II by setting their orbital separations. Fig. 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 sys-tems converges to a value around 13-15 R Hill after 200-400 ky. The orbital separation is a free parameter from Criterion II of the analytical model that we have chosen 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. 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). Due to 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 one The red dots indicate the time for each distance from the star at which a planetary embryo is placed using our analytical model. The orbital seperation as input to the analytical model is given as 17 R Hill with a randomization of 2.5 R Hill 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 10 6 years from 0.5 au to 5 au (6M ⊕ ,13M ⊕ ,27M ⊕ ) and the planetesimal surface density slope (Σ P ∝ r −1.0 , Σ P ∝ r −1.5 , Σ P ∝ r −2.0 ). attempts to compare two LIPAD runs with similar initial conditions, 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. 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 number of active embryos after 1 Myr is between 30 and 40 embryos for 8 out of our 9 runs. Only the 6 M ⊕ and Σ P ∝ r −1.0 run contains less 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 number of embryos does not simply increase for more massive planetesimal disks in our runs. The reason being that the embryos that form grow larger in more massive disks. They thereby increase their orbital separation to the other embryos again. While larger 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 Myrs we can see that 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. 2 -Fig 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. Fig. 2 -Fig. 4 clearly show that embryo formation for every power law planetesimal surface density profile occurs from the inside out. This is an expected result due to the shorter growth time scales in the inner disk and the correspondingly higher densities in planetesimals. 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 embryos individual locations, even though following 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 do well in constraining the zone of their individual formation. Another effect that can be found is that embryos increase their orbital distance to other embryos when they grow in mass. This effect has already been found and discussed by Kokubo & Ida (1998) and Kobayashi et al. (2011) when studying the oligarchic growth of massive objects. In the general picture, initial embryos begin to form the 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 studied, as can be seen in Fig 8. This directly results in a cumulative number of embryos that scales logarithmic with distance, as can be see in Fig. 7. As 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 it it shown in Fig. 9. The number of active embryos even decreases slightly for higher disk masses and steeper planetesimal surface density profiles. As Fig. 2 -Fig. 4 show, the area in which planetary embryos form becomes larger for higher masses and steeper density profiles. Since their orbital separation increases for higher masses and since the mean orbital distance converges to the same number of Hill radii (Fig. 8), we see that the total number of embryos within 1 Myrs also does not sensitively abbreviate for different input parameters.

Implications for pebble accretion
While the effect of pebble accretion on the formation of planetarey 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. Since the growth rate of planetary em- bryos 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 since the formation of planetesimals requires a radial pebble flux, using our setup we can estimate first constraints on said pebble flux and subsequently on the possibility of continuous pebble accretion. Even though we do 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 already displayed by several studies like Ormel & Klahr (2010), Bitsch et al. (2015) and Ndugu et al. (2017) to mention just a few. The efficiency of pebble accretion directly depends on the local pebble flux at the location of an accreting body of sufficient mass. Since 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 6 years of the systems evolution. However, since 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. There has to be a certain embryo size reached at first to effectively accrete pebbles. Another crucial impact on the pebble flux evolution is the formation of planetesimals itself, since they form based on the disks evolution. The more planetesimals form, the earlier we also form planetary embryos, that could accrete pebbles. However, the more planetesimals one forms, the lower the pebble flux would become, due to the mass transfer 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. Since ∼ 90% of our planetesimals form within 400ky of our setup, we conclude that the magnitude of the pebble flux has decreased significantly before that time. Embryos that form after 400ky would therefore not be able to undergo significant pebble accretion in our model. We find that in our setup most embryos that form within 400 ky also form within 1 au. The formation of further out embryos around 1.5-2.0 au occurs well after 400 ky. In conclusion, it is not possible for further out planetary embryos to undergo pebble accretion in our setup. This statement yields true if one assumes a power law distribution for the planetesimal surface density like the minimum mass solar nebula hypothesis, the dust profile of a viscous disk, or the pebble flux regulated planetesimal formation surface density profile.

On the architecture of planetary systems
Following up on our findings from Sect. 5.2 it is not too farfetched 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 assumes power law density profiles for the planetesimal surface density and our results of inside out planetary embryo formation is a direct consequence of this. If one would assume abbreviations from the power law profile due to local substructures in the the disk, like e.g. around the iceline (Drążkowska & Alibert 2017), this picture might change. The early formation of planetary embryos around the iceline could lead to the formation of cold giant planets via pebble accretion. The formation of those planets can then have major consequences to 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, as they would reduce the pebble flux that reaches the terrestrial planet region. Also the additional planetesimal formation itself will 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 effect 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. Another scenario 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 planetary embryos at much larger distances within the lifetime of the pebble flux.

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 2 out of our 9 runs, the deviation of the total number of embryos is below 5 %, in 4 out of 9 it is below 10 % and in 8 out of 9 it is below 25 %. Only the 6M ⊕ , Σ P ∝ r −1.0 run deviates stronger (≈ 40%). The Hill crite-rion 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 does well in handling the analytic planetesimal surface density evolution (Sect. 4.3). It is worth mentioning that the N-body simulations require weeks (sometimes months) of computation time with the same planetesimal input, whereas the parameterized model takes merely seconds. While the N-body simulations clearly involve more complexity that allow for a more complete picture of the problem, the question on where, when and how 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, like e.g. planet population synthesis. Even though our study focused on an area from 0.5 au to 5 au, the analytical model should also yield true at further locations and could be a valuable asset in considering planetary embryo formation in far out ring-like structures of circumstelar disks, as seen in ALMA observations. Other studies regarding planet formation via 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.

Summary & Outlook
We study the spatial distribution and formation timescales of planetary embryos from an initial disk of gas and dust. For this purpose, we couple 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 is given as 100 km in diameter and dynamically grows due to collisions with other planetesimals. We analyze the first million years of nine different systems in which we vary the total mass in planetesimals and their surface density profile. In combination with analytic estimates on growth rates of planetesimals based on their local surface density, we derive an analytical model for planetary embryo formation. Our model does well in reproducing the spatial distribution and formation time of planetary embryos. We use 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, population synthesis etc. ) 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: -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(>2au) fail to do so in any disk 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 link 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 find 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, since as we show, the presence of a planetary embryo and the presence of 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 one forms, the earlier one forms planetary embryos, but the more planetesimals one forms, the less mass remains in pebbles. Vice versa, if one decreases 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 our presented approach will study the formation of embryos in far out planetesimal rings that could result from pressure bumps during the disks evolution.