© The Authors 2023

1 Introduction

As of February 2023, more than 5200 exoplanets have been confirmed from stellar observations¹, and several thousand planetary candidates discovered by missions such as Kepler, K2, and TESS are awaiting confirmation. The majority of these planets are so-called super-Earths and mini- or sub-Neptunes, which fall between Earth and Neptune in radius and/or mass; such planets do not have an equivalent in our Solar System. Almost all known exoplanets orbit main-sequence stars of spectral categories F, G, K, or M, and the current programs searching for exoplanets tend to concentrate on such stars. Generally, Earth-mass planets are difficult to detect because of their low mass compared to the star, and Earth-mass planets on Earth-like orbits are even harder to detect because it requires long pointing periods from telescopes (e.g., Petigura et al. 2013b) and extremely high-precision measurements. However, this situation is slowly changing. New space telescopes, such as the upcoming ESA PLATO mission, aim to detect and study thousands of exoplanets. The ESA-stated primary goal of the PLATO mission is the “detection and characterization of terrestrial exoplanets around bright solar-type stars, with emphasis on planets orbiting in the habitable zone”. F, G, and K dwarfs are considered to be “solar-like stars”. This planned mission, among others, will provide exciting new opportunities for discovering exoplanets with Earth-like orbits and sizes; we must prepare our science accordingly.

While the new space telescopes allow us to detect and study exoplanets, the ground-based telescope Atacama Large Millimeter/sub-millimeter Array (ALMA), which is currently the largest ground-based astronomical project on Earth, contributes to our understanding of the formation of planets by imaging the gas and dust in their protoplanetary disks (ALMA Partnership 2015). The initial conditions for planet formation depend on the properties of the protoplanetary disk the planets are forming in. These properties vary with the spectral type of the central star, amongst other things, which indicates that planet formation likely also varies across the stellar spectral types (Andrews et al. 2013). This motivated us to study planet formation around different types of main-sequence stars. In this paper, we aim to explore the planet formation process and the resulting most common trends of the exoplanet population orbiting K-dwarf stars. Specifically, we model terrestrial planet formation in the inner disk of this type of star. Additionally, we tested the validity of assumptions commonly made for conditions in the early Solar System.

Several modeling studies have focused on planet formation around M dwarfs (e.g., Laughlin et al. 2004; Lissauer 2007; Raymond et al. 2007; Miguel et al. 2020; Zawadzki et al. 2021), since they are the most common stars in our Galaxy (e.g., Henry et al. 2006; Winters et al. 2014). For this study we focus on K-dwarf stars, which are main-sequence stars with masses and radii between those of red M dwarfs and yellow G-type stars, such as our Sun. They have masses between 0.5 and 0.8 times the mass of the Sun and surface temperatures between 3900 and 5200 K (Pecaut & Mamajek 2013)². K dwarfs are not as common as M dwarfs (which comprise 75% of the main-sequence stellar population in our Galaxy); they comprise about 15% of the main-sequence population, meaning they are more abundant than G and F stars, which comprise 7% and 2%, respectively (Safonova et al. 2021). Their lifetimes are longer than those of F and G stars; therefore, they are of particular interest in the search for extraterrestrial life. Using N-body simulations, we examined how the properties of planetary systems around K dwarfs are established from planetesimal accretion, that is, during the late stages of assembly when planets accrete from planetesimals and planetary embryos.

The currently known population of planets around this type of star is characterized by compact multi-planet systems of mostly small, dense planets with short periods and generally low orbital inclinations (e.g., Fang & Margot 2012). This characterization could be partially due to observational biases, since these planets have been mostly found via the transit method. Additionally, it can be explained by the fact that the rate of solid accretion is proportional to the surface density of planetesimals; therefore, the planetary formation and evolution mostly happens inside and around the snow line (i.e., within a few AU) because at this location there could be a sudden increase in solids (e.g., Stevenson & Lunine 1988; Ciesla & Cuzzi 2006; Schoonenberg & Ormel 2017). K dwarfs are cooler than G-dwarf stars, so their snow line is closer to the star. Therefore, we tested the hypotheses that the most massive planetary embryos and protoplanets form there and that the final systems tend to be compact with short orbital periods. This is the reason our simulations are focused on the inner part of the disk around the star.

The 46 observed multi-planet systems around K dwarfs are shown in Fig. 1 in comparison to the terrestrial Solar System planets. The sample was retrieved from the NASA Exoplanet Archive³. The systems contain at least two planets each (single planets orbiting K dwarfs are not included in the sample), and there are 139 planets in total (only planets with known masses are presented). Out of these 46 systems, 15 contain at least one giant planet (shown in orange in the figure). Here, we use the definition by Clanton & Gaudi (2014), who suggest that “giant planets” need to have masses of ~30 M_⊕ or more. In our sample, about 30% of the systems contain one or more giants, but since giant planets are typically easier to detect, we assume that this percentage may actually be lower for typical systems around K dwarfs. This assumption is supported by studies that show that giant planets with periods shorter than a few years detected around Sun-like stars have an occurrence rate of only around 10% (Cumming et al. 2008; Winn & Fabrycky 2015). Other studies estimate that giant planets specifically orbiting late K-dwarf stars with masses of ~0.5–0.75 M_⊕ have an occurrence rate of under 1% (Gaidos et al. 2013; Bryant et al. 2023). Even though not all studies use exactly the same classification for giant planets, it is clear that these planets are not frequent around this type of star.

It is likely that at least some of the systems in the retrieved sample are not complete, as suggested by studies of completeness of Kepler systems. For example, Zink et al. (2019) estimate the average number of planets to be ~6 per G or K dwarf within the radius and period parameter space of Kepler detections, whereas the systems in our sample contains on average only three planets (3.0 ± 1.3). Giant planets in such systems would probably have already been discovered, so there is a higher probability that they are mostly made up of less massive distant planets that have not yet been detected. Earth-mass planets are currently only detectable when they are on very short orbits around lower-mass stars. In our sample, more than 60% of the planets are “super-Earths” or “sub-Neptunes” with semimajor axes mostly between 0.04 and 0.2 AU. Almost all the discovered planets are much more massive than Earth and orbit their stars more closely than Mercury. In fact, the vast majority of the known systems are remarkably different from the Solar System (e.g., Petigura et al. 2013a; Winn & Fabrycky 2015; Bryan et al. 2019). The type of planet with the highest occurrence rate in K-dwarf systems is absent from the Solar System, but it is currently the most abundant class of exoplanet. Planets with radii between 1 and 4 R_⊕ and orbital periods shorter than 100 days are present around at least 30%–50% of main-sequence stars of spectral types G and K (Mayor et al. 2009, 2011; Howard et al. 2010, 2012).

A typical known K-dwarf planetary system does not contain giant planets, but it comprises on average three planets, mostly super-Earths or sub-Neptunes with very short semima-jor axes. This characteristic is in large part due to observational biases, and we expect the existence of thus far undetected, less massive planets to mainly orbit farther from the host stars. In the next sections, our simulated systems are compared to this retrieved sample in order to evaluate our results. Our study is focused on the formation of terrestrial-type planets (Earths and super-Earths). As such, we do not analyze the distribution of planets with masses above 10 M_⊕, the upper mass limit for super-Earths (Stevens & Gaudi 2013), in the observed systems when comparing them with our simulated systems.

The interaction between the growing planetary embryos and the gas in the protoplanetary disk induces a torque that causes these embryos to migrate toward the star. The onset of migration occurs when the planets reach approximately the mass of Mars (i.e., about 0.1 M_⊕; Ward 1986; Tanaka et al. 2002; Paardekooper et al. 2011). The observed planets on very short orbits described earlier suggest that orbital migration (Type I migration in the case of smaller planets) was important during the evolution of these planetary systems, as it is not likely that these planets formed at their current orbits (e.g., Izidoro et al. 2014). After the gas disk dissipates, the migration stops and rocky planets may continue growing via collisions between embryos (Agnor et al. 1999) because their gravitational interactions might lead to an excitation of their eccentricities, making collisions possible. The gas disk damps their eccentricities and inclinations, so when the gas is gone, the evolution of embryos gradually becomes chaotic.

Our understanding of terrestrial planet formation has improved over the last few decades, but the formation of multiple close-in super-Earths around stars with lower masses than the Sun is still an open issue. Specifically, the mechanism that prevents these planets from migrating too close to their host star is not yet satisfactorily known due to the uncertain effects of torques near the disk’s inner edge (e.g., Brasser et al. 2018). The gas disk does not extend all the way to the star, probably due to the star’s magnetosphere. The magnetic field disrupts the disk out to a radius at which the magnetic energy density and the kinetic energy density of the gas are equal (e.g., Long et al. 2005). At this radius the Type I migration slowly stops, so this gas-free cavity works as a migration trap where there is no migration, nor eccentricity or inclination damping. For a typical T Tauri star, the magnetospheric boundary is located at around 0.05–0.1 AU (Romanova & Lovelace 2006; Romanova et al. 2019). Migrating protoplanets usually end up captured in mean motion resonance (MMR; Terquem & Papaloizou 2007). In multiple N-body simulation studies, several planets migrate inward together as a resonant chain of low-mass protoplanets. When the innermost of them reaches the inner disk, their migration slows until they eventually stall near the disk’s inner edge (Ogihara & Ida 2009; Cossou et al. 2014; Coleman & Nelson 2014, 2016). Simulations also show that either the further evolution of these planets will keep them in resonances and dynamically stable for billions of years (Esteves et al. 2020) or they will become dynamically unstable, typically at the end of or shortly after the gas disk dispersal. Late dynamical instabilities break resonances and often lead to orbital crossings and giant impacts (Izidoro et al. 2017, 2021). This scenario is commonly referred to as “breaking the chains” (e.g., Terquem & Papaloizou 2007; Ogihara & Ida 2009; McNeil & Nelson 2010; Cossou et al. 2014; Izidoro et al. 2017, 2021). The final configurations of systems that underwent late dynamical instabilities are expected to be nonresonant. Most currently known super-Earths are not found in resonant systems (Lissauer et al. 2011; Fabrycky et al. 2014). On the other hand, systems such as TRAPPIST-1 and Kepler-223 have several planets in resonant chains, and, according to multiple studies (Cossou et al. 2014; Izidoro et al. 2017, 2021; Ogihara et al. 2018), they represent a small fraction of systems that did not become unstable after the dispersal of the gas disk. Resonant systems are naturally produced by migration; therefore, this is an important effect to consider when it comes to terrestrial planet formation.

In the next sections, we first explain the methods used in this study and introduce our model. We describe the protoplanetary disk model and density profile together with the initial conditions for each of our setups. Then we present and discuss the results of our simulations and the dynamical evolution, final architecture, and characteristics of the systems we have built, as well as their significance and implications. We also compare the simulated systems to observations. Finally, in the last section, we present the conclusions of our study.

Fig. 1 Currently known multi-planet systems around K-dwarf stars. Only planets with known masses are presented. The sample was retrieved from https://exoplanetarchive.ipac.caltech.edu/ on December 4, 2022. Giant planets with masses > 30 M⊕ are shown in orange, and all planets with lower masses in blue. Terrestrial Solar System bodies are included for comparison (in red). The point size indicates the masses in M⊕ of the observed planets. Note the different size scale of the different colored bodies (blue and orange). Additionally, the markers representing the Solar System bodies are enhanced by a factor of 10 with regard to the planets indicated in blue for better visualization.

2 Methods

For the N-body simulations we employed the GENGA software (Grimm & Stadel 2014; Grimm et al. 2022). GENGA is a GPU implementation of a hybrid symplectic N-body integrator, developed based on another N-body code MERCURY (Chambers 1999), for simulating planet and planetesimal dynamics in the final stages of planet formation, and the evolution of planetary systems. It integrates planetary orbits over long timescales with excellent energy conservation. Gravitational interactions between planetary bodies are computed as perturbations of the Keplerian orbits (Wisdom & Holman 1991). The GENGA code has been successfully used for simulations of the terrestrial planets formation in the Solar System (e.g., Clement et al. 2020; Woo et al. 2021, 2022). Planet formation simulations usually start from the runaway growth phase when planetesimals collide and form bigger objects called planetary embryos (Kokubo & Ida 1995). Then they continue with oligarchic growth (Kokubo & Ida 1998), when a small number of largest embryos begin to dominate the dynamics of the disk until they reach their isolation mass by accreting all planetesimals in their feeding zones, and eventually protoplanets and planets form (e.g., Kokubo & Ida 2000; Chambers 2006).

We model the formation of terrestrial planets around a 0.6 M_⊙ and a 0.8 M_⊙ star from planetesimals and planetary embryos to evolved planets using planetesimal accretion. We do not actually expect any substantial differences in the simulated outcomes for the different stellar masses, but we have chosen these values to cover most of the K-dwarf stars mass range, which will be potentially useful for our future studies. We have not adapted the structure or characteristics of the disk to account for the lower or higher stellar mass as we assume these to be just second-order differences. In the end, we do not draw any conclusions for the outcomes (and their comparison to the observed systems) based on the different stellar masses as the uncertainties on exoplanet radii, masses, and orbital parameters do not allow for it.

In this study we do not attempt to reproduce any specific observed planetary systems. Rather, we consider the generic outcome of planet formation for this type of star and the results of our simulations are then compared to the known systems (as described in Sect. 1). All together we ran 48 high-resolution N-body simulations that varied the initial protoplanetary disk mass, and solid and gas surface density profiles. We investigated whether the planetary systems that were formed allow us to decide which initial conditions and other parameters are responsible for their final characteristics. Simulations start with planetesimals and planetary embryos in their predefined initial locations. Then the bodies begin interacting gravitationally around the star, which results in perturbations of their orbits, their close encounters and collisions, and often ejections of some of the bodies. We followed their collisional growth for up to 20 Myr, until a system of planets or protoplanets (and some leftover planetesimals) is formed, where the planets generally do not collide with each other anymore.

Current computational hardware including recently improved GPU hardware is increasingly powerful and allows us to run the necessary calculations for planet formation processes faster than ever before. Nevertheless, the high-resolution long-term simulations still take a relatively long time to run; therefore, we limited the simulations to 20 Myr. This time period is long enough for our purposes as previous studies already show that this type of planet forms quickly. For example, Ogihara et al. (2015) examine the formation of close-in super-Earths around G dwarfs using N-body simulations and find that Type I migration can cause planets to form in only a few megayears. Our N-body simulations include the gas disk effects (Type I migration as well) on the formation of this kind of planet (close-in super-Earths), but around a bit less massive stars. We therefore expect that planets around K-dwarf stars may form within a similar time frame.

2.1 Initial conditions and parameters

All the simulations begin with 12 000 bodies grouped in radius between 200–2000 km depending on the initial mass of the disk of solids, which we vary between 5 M_⊕ and 40 M_⊕. Even though this disk mass may seem excessive, recent studies have been performed with similar initial masses. For example, Morbidelli et al. (2022) manage to produce the total mass of planetesimals exceeding 40 M_⊕ at the silicate sublimation line in the forming Solar System (at ~1 AU, but closer to the star for K dwarfs) by reducing the viscosity parameter a and suggested that this large amount of mass could explain the formation of the frequently observed super-Earths. In our simulations all planetesimals are fully self-gravitating. The higher the initial disk mass, the larger (and therefore more massive) bodies we start with to preserve the number of 12 000 particles. The planetesimal number was chosen based on a series of simulations testing the speed of the code. Since more bodies result in more calculations and longer computational time scaling roughly as N² (Grimm et al. 2022), using 12 000 particles achieves a high-enough resolution, but still within a reasonable time limit of around 2 months and available computation time. The fact that simulations with higher initial disk masses start with larger bodies should not affect the simulation outcome. According to Miguel et al. (2020) the initial size of planetesimals does not affect the final population, at least for bodies between 100 and 500 km, and our test-runs showed the same for larger sizes. With larger initial bodies the planetary growth happens faster (e.g., Woo et al. 2021), but the formation timescale is not the main focus of this study. At the beginning of a simulation, the masses of the planetesimals and embryos are calculated from the initial radii and density; the latter is chosen to be 3 g cm⁻³. This value is commonly assumed for rocky planetesimals formed in the inner part of a protoplanetary disk. As planetesimals collide and grow, the radius of the new planetary body is set by conserving the mass and mixing the densities of the merged bodies. Since the density is the same for all the particles, the final planets will all have the same densities if any compressional effects are ignored; hence, these simulations do not allow us to properly assess the radii of the final bodies. This is a limitation that will be improved in the future studies.

The planetesimals are distributed between 0.2 and 2 AU with a surface density (Σ_P) profile Σ_P ∝ r⁻¹ for the initial disk. Models for the viscous evolution of the gas in the disk suggest a shallow planetesimal surface density slope of Σ_P ∝ r⁻⁰⁹ (Shakura & Sunyaev 1973), whereas the minimum mass solar nebula (MMSN) hypothesis provides a steeper profile of Σ_P ∝ r⁻¹⁵ (Weidenschilling 1977; Hayashi 1981). Other similar studies use these surface density profiles or even a steeper slope of Σ_P ∝ r⁻²¹ from Lenz et al. (2019). Zawadzki et al. (2021) run their simulations with Σ_P ∝ r⁻⁰⁶ and Σ_P ∝ r⁻¹⁵, and show that the solid surface density of the final planetary systems appear to be independent from the initial distribution of embryos. Therefore, we decided to use a mean value of the last two slopes. Figure 2 shows the initial solid surface density distribution for simulations with the starting planetesimal locations between 0.2 and 2 AU, and different initial disk masses. The values of surface density at 1 AU are in the range from 11.8 (for the initial disk mass 5 M_⊕) to 94.4 g cm⁻² (for the initial disk mass of 40 M_⊕). The higher-end values are significantly higher than in other studies, for example Mulders et al. (2020) do not exceed 50.0 g cm⁻² in their N-body simulations around a solar mass central star. We extended the range to these hypothetical values to study whether massive planets can form in such conditions.

During the simulations, bodies are removed when they come closer to the central mass than the inner truncation radius (they are assumed to collide with the star), or move farther out than the outer truncation radius of the disk, or when they collide and the less massive body merges into the more massive one. We assume perfect accretion during collisions, ignoring fragmentation, since it seems that this effect does not radically alter the outcome of the simulations (Kokubo & Genda 2010; Chambers 2013; Quintana et al. 2016). The formation of the planetesimals themselves is suggested to happen very early in the lifetime of a protoplanetary disk, already after a few thousand years at 2– 3 AU (Lenz et al. 2019). In the case of large planetesimals or planetary embryos, age determination of iron meteorites, representing the metal cores of differentiated planetesimals, suggests that these were formed in the inner Solar System within the first megayear (Halliday & Kleine 2006; Kruijer et al. 2014; Schiller et al. 2015). This time period is still short, so for our purposes we set the start time of the simulations to 0 and assume that planetesimals and embryos have not moved much while they were forming. In this manner their initial distribution of orbital elements reflects that of a young disk. The planetesimals had initial values of the inclination and orbital eccentricity randomly ranging between 0 to 0.1 degree and 0 to 0.01, respectively, but these values increase quickly after the simulations start due to particle interactions. An example of the initial planetesimal/embryo distribution with their masses and radii is shown in Fig. 3, where the different colors indicate mass ranges of the different objects (i.e., planetesimals, protoplanetary embryos, and embryos). Planetary embryos are often defined as objects of at least the mass of the Earth’s Moon, M_Moon = 0.012 M_⊕ (e.g., Woo et al. 2021, 2022; Voelkel et al. 2021), we also use this definition. There is no difference in treatment of these objects during the N-body simulations.

Our simulations represent a disk model with the particles initially distributed between 0.2 and 2 AU and the inner truncation radius of 0.05 AU; in other words, the inner edge of the simulations is set to quarter of the semimajor axis of the innermost body at the start of the simulations. All other parameters are then varied in the model. The time step necessary for the N-body integration runs depends on the orbital period of the possible innermost object. It is recommended to use 1/20 of the orbital period at the cutoff distance for sufficient accuracy (Wisdom & Holman 1991). However, a shorter time step results in a longer computation time and such N-body simulations already require many weeks, even months of computation time. Therefore, we used a longer time step than recommended. Results that are presented in the following sections do not seem affected by this compromise. This issue is discussed more in Sect. 4.3, where we present results of a run with the recommended time step. Table 1 summarizes the main parameters and initial conditions of our model. The inner and outer truncation radius are cutoff distances for the simulations; bodies with a separation to the central mass smaller or larger than these values are taken out of the simulation. Initial semimajor axes a_min and a_max limit the planetesimal spatial distribution at the beginning of a simulation. We varied the initial mass of the solid disk and gas surface density profile, presented in Table 2, and simulated 20 Myr of planetary evolution around two different central masses. We ran all combinations of these parameters.

The initial conditions for planet formation depend on the properties of the protoplanetary disk. For the estimates of the initial mass of the protoplanetary disk we followed Manara et al. (2018). In their study, they focus only on the solid part of the disks and present recalculated disk dust masses (the original disk masses are from Pascucci et al. 2016 and Ansdell et al. 2016) around various spectral types of stars, including K dwarfs with different masses, located in two young star-forming regions. The dispersion of the disk mass estimates is large, and for K-dwarf stars this ranges from less than 5 M_⊕ to more than 70 M_⊕, for one of the star-forming regions and to more than 110 M_⊕ for the other one. But for the majority of the disks the dust masses are less than 10 M_⊕. Of course, planetary systems are concentrated within a few AU of the central star, whereas most of the disk mass is farther out; therefore, we should consider only a small portion of this estimate for our simulations. On the other hand, these estimated disk dust masses might be highly underestimated due to multiple reasons, as several studies have pointed out. For example, they are derived from disk measurements that are only sensitive to certain grain sizes, so it is possible that a significant fraction of the dust is undetected (e.g., Williams 2012). Also, the estimates depend on the values of dust opacity and disk temperature, which are still debatable (e.g., Andrews et al. 2013; Pascucci et al. 2016). Some dust mass might be confined to optically very thick inner regions of disks (e.g., Tripathi et al. 2017); however, the measurements are based on the emission from the outer regions (R > 10 AU) of disks (e.g. Bergin & Williams 2017). After taking all these uncertainties into account, we chose initial solid mass values in a range of 5 to 40 M_⊕ for the simulations.

Fig. 2 Initial solid surface density distribution for simulations with the starting planetesimal locations between 0.2 and 2 AU and different initial disk masses (IDM). The increasing linewidth with the values on the y-axis represents the growing fraction of larger bodies in the disk with increasing initial disk mass.

Fig. 3 Initial distribution, masses, and radii of 12 000 bodies for a simulation with an initial solid disk mass of 10 M⊕. The x-axis represents the distance from the central star in AU, and the y-axis the radii (left panel) and masses (right panel). Red, blue, and black indicate different mass ranges of the objects, i.e., mass M < 10−3 M⊕ (considered as planetesimals), 10−3M⊕ < M < 10−2 M⊕ (proto-embryos), and M > 10−2 M⊕ (embryos), respectively. There is no difference in the treatment of these objects during the N-body simulations. The isolation mass, which is the mass of embryos that accreted all the bodies within their feeding zone, is represented by the dashed line.

2.2 Gas surface density profile

The gas disk implementation in GENGA follows that of Morishima et al. (2010). GENGA supports all important gas disk effects, such as aerodynamic gas drag, disk-planet interaction including the gas disk tidal damping, Type I migration, and the global disk potential. All these effects are included in the simulations. The gas surface density dissipates exponentially in time (t) and uniformly in space (r) following $${\displaystyle {\sum }_{\text{gas}}\left(r,t\right)=}{\displaystyle {\sum }_{\text{gas,0}}{\left(\frac{r}{1\text{\hspace{0.17em}AU}}\right)}^{-1}\exp \left(-\frac{1}{{\tau }_{\text{decay}}}\right)}.$$(1)

In Eq. (1), Σ_gas,0 is a gas surface density (Σ_gas) at 1 AU and t = 0. The MMSN, which is used to define the initial conditions of the protoplanetary disk required to make the planets around the Sun, has Σ_gas,0 = 1700 g cm⁻² based on Hayashi (1981), whereas Morishima et al. (2010) assume Σ_gas,0 = 2000 g cm⁻². We decided to test four different values: 1000, 1250, 1500, and 1750 g cm⁻², which are mostly lower than assumed MMSN values and reflect the lower masses of our stars compared to the Sun. The quantity τ_decay represents the gas disk decay timescale, which we fixed at 1 Myr. Woo et al. (2021, 2022) examine decay timescales of 1 Myr and 2 Myr for simulations of the Solar System formation, and show that decay timescale ≤ 1 Myr can better explain the specific characteristics of our planetary system. They also show that the mass loss to the Sun due to Type I migration is lower in the case of the shorter decay time. Since both our stars are smaller than the Sun, and we are trying to form more massive terrestrial planets than exist in the Solar System, we chose the shorter decay time of 1 Myr. However, it is important to mention that real disks may not dissipate smoothly at all orbital radii, instead inside-out dissipation can alter the orbital architectures of planetary systems (Liu et al. 2017) and even trigger instability (Liu et al. 2022). After the gaseous part of the protoplanetary disk dissipates, gas drag and gas dynamical friction disappear as well. This causes the accretion rate to slow down and eccentricity and inclination damping of planetary embryos will cease too (e.g., Bitsch & Kley 2010). In the version of GENGA employed here, the gas disk inner edge is hard-coded to 0.1 AU. As mentioned in Sect. 1, the magnetospheric boundary is located at around 0.05–0.1 AU. The gas disk inner edge represents the boundary in our simulations and we chose to keep this value of 0.1 AU for our runs.

This gas surface density profile is used in our model, the inner and outer truncation radius are at 0.05 AU and 3 AU, respectively. That means that, since the gas disk inner edge is set to 0.1 AU, there is a cavity without gas between the inner truncation radius and the inner edge of the gas disk. Many known planetary systems around K-dwarf stars contain planets with orbital distances shorter than 0.1 AU (see Fig. 1); therefore, we decided to set the inner truncation radius to 0.05 AU. Figure 1 shows that some planets are observed even closer to their host stars than 0.05 AU; however, due to limited computation time we did not use an even lower value of the inner radius.

3 Results

In this section, we present our simulated planetary systems around the K-dwarf stars after 20 Myr of planet formation. In total, we performed 48 runs with different combinations of parameters with the target of finding parameters that reproduce the known K-dwarf system characteristics. We varied the initial protoplanetary disk mass, and solid and gas surface density profiles. But first we start by discussing the dynamical evolution of the simulations.

3.1 Dynamical evolution of the simulated systems

The dynamical evolution (20 Myr) of a typical simulation is presented in Fig. 4. The panels on the left show the evolution of semimajor axes and masses of all bodies with a mass > 0.1 M_⊕, and the panels on the right present eccentricities and orbital inclinations, and show only bodies with masses > 0.5 M_⊕ to better visualize these orbital characteristics. The top left plot shows a typical evolution of such a simulation based on the migration model as described in Sect. 1. During the gas disk phase, planetary embryos grow and migrate inward, occasionally end up falling onto the central star as we see just before 2 Myr. The planets quickly settle into chains of MMRs. According to Izidoro et al. (2017), the chains are typically established within 1.5 Myr; this is the case in our simulations as well. As the gas in our disk decays exponentially with the decay time of 1 Myr, it will never completely dissipate, but we see that at ~5 Myr the amount of gas in the disk must be relatively low as migration mostly stops. Additionally, eccentricity and inclination damping disappears together with the gas and this is clearly demonstrated in both eccentricity and inclination plots at around 7 or 8 Myr when both properties, but inclination particularly, increase. Therefore, we estimated the time of the (most) gas dispersal at ~8 Myr. After this time, the plots display several dynamical instabilities occurring at around 11 Myr and later. They lead to orbit crossings, collisions and merging of some of the planets, as displayed in both left plots. Almost all our systems undergo such instability phases, but usually evolve to a stable configuration fairly quickly (in the order of 10⁵ yr). Eventually, the planetary system becomes more dynamically excited and less compact; however it still contains eight planets, six of them interior to 1 AU. In this simulation, the three most massive planets have several M_⊕ and the rest are under 1 M_⊕. Both the first and the third (from the central star outward) subsequent planet pairs are in 2:1 resonance. The final orbital eccentricities are typically around 0.05 (median = 0.05 ± 0.02) and inclinations around 1.42 deg (median = 1.42 ± 0.74 deg). These values are within the current eccentricities and inclinations of the three largest terrestrial planets of the Solar System.

While Fig. 4 displays the dynamical evolution of a typical simulation that underwent several dynamical instabilities after the gas disk dissipation, Fig. 5 shows a less common evolution of a system that did not do through any instabilities. The panels again show 20 Myr of evolution in semimajor axes, masses, eccentricities and inclinations, both before and after the gas dispersal at ~8 Myr. The plots show some collisions after the gas disk phase as well, but they are very rare and involve only small leftover bodies farther from the central star. No dynamical instabilities that would affect the architecture of the system are visible, but we still see some planet growth in the bottom left plot, which displays bodies with masses > 0.1 M_⊕. After 20 Myr, the final planetary system consists of 11 planets (> 0.5 M_⊕), with the most massive planets being half the mass of the most massive planets in the previous system (Fig. 4). This system is also more compact than the previous example. The inner planets are in a 4:3-4:3-4:3—5:4 resonant chain, that is, the three innermost subsequent planet pairs are all in 4:3 period ratio (P_out/P_in) and the fifth pair is in 5:4 resonance. Multiple other pairs are very close to being in resonance. In this study, we allowed for a 1% deviation of a resonant loci away from the exact commensurability based on Brasser et al. (2022, Fig. 2), independently of the eccentricities. The final orbital eccentricities are typically around 0.02 (median = 0.02 ± 0.03) and the inclinations around 0.84 deg (median = 0.84 ± 1.15 deg). These values are a bit lower than the typical values for the system that underwent several dynamical instabilities, which is expected as the instabilities excite the planets. But the values are still very similar, which can be explained by the fact that the previous system had enough time to stabilize again after the last instability phase. The dynamical evolution examples presented in Figs. 4 and 5 are representative of all our simulations. Although we present the most contrasting examples, we see that their evolutionary paths show many similarities, hinting at a generic pathway.

To demonstrate that a dynamical instability can also occur much later in the simulation time, even in a simulation that did not undergo an instability after the gas dispersal during the standard time of our simulations (20 Myr), we present one example of a few simulations that were extended to 100 Myr (see Fig. 6). The plots show several collisions resulting in planetary growth at approximately 21 Myr and then later at around 28 Myr. The final outcome of the simulation is a less compact, dynamically stable system (for 70 Myr) with 5 super-Earths and 3 Earth-sized planets with orbital eccentricities reaching values of ~ 10% and inclinations up to a few degrees, so within the current eccentricities and inclinations of the three largest terrestrial planets of the Solar System.

To determine which systems underwent a late instability phase with giant impacts and which did not, Izidoro et al. (2017, 2021) use the distributions of the last collision epochs in their simulated systems. They subsequently divide the systems into “stable”, with no collisions occurring after the gas dispersal, and “unstable” with some collisions happening afterward. We cannot use this criterion for our simulations as basically all our systems experienced collisions after the gas disk phase, and often long after. Most such collision events happened at much lower collision frequency than a dynamical instability that results in disruption of resonant chains, and thereby significantly changing the architecture of a system. Also, we cannot state when exactly the gas disk phase ends, whereas in their simulations it ends abruptly. In addition, their simulations start with only embryos, no planetesimals that exert dynamical friction on protoplanets and can damp their eccentricities and inclinations (Wetherill & Stewart 1993; Ida & Makino 1993; O’Brien et al. 2006). Therefore, it is difficult to compare the evolution of our versus their simulated systems. According the collision criterion of Izidoro et al. (2017, 2021), all our systems would be classified as unstable, even though some of them remained in resonance and clearly appear to be long-term stable. The histogram in Fig. 7 shows a number of embryo-embryo collisions in all our simulations during their 20 Myr of evolution. The number of collisions decreases mostly smoothly from the start of the simulations until ~8 Myr. Then the number stays at a continuous low level of collisions and the collision rate is roughly constant, with occasional potential dynamical instabilities until the end of the 20 Myr period. We say “potential” because it is questionable whether the spikes in the number of collisions are actually statistically significant. Our systems that get disturbed by the instabilities then stabilize relatively quickly and readjust, and those with a lower number of planets sometimes become resonant again (some gas might still be present in the system). Additionally, we extended one-third of all our runs to 40 Myr and one-fifth to 100 Myr. At 40 Myr, we see that in around one-third of the extended simulations one or several instabilities occurred after 20 Myr of evolution (see again Fig. 6); while between 40 Myr and 100 Myr, we see no more instabilities. Izidoro et al. (2017) show that at least 75% (and probably 90–95%) of their simulated systems must be unstable to match the period ratio distribution of observations. By investigating each simulation and its dynamical evolution individually, we recognized at least ~85% of our systems that underwent at least one late dynamical instability, which resulted in the resonant chains becoming unstable. However, due to the reasons mentioned previously, the number is a bit uncertain.

Fig. 4 Dynamical evolution of planets in one of the typical systems that underwent several late dynamical instabilities after the gas dispersal. The plots show the temporal evolution (20 Myr) of planetary bodies, specifically their semimajor axes, masses, eccentricities, and orbital inclinations. Each line color represents an individual object; however, the color is random in each plot. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. For better visualization, the panels on the left (semimajor axes and masses) show the evolution of all bodies with masses > 0.1 M⊕, and the panels on the right (eccentricities and orbital inclinations) show only bodies with masses > 0.5 M⊕.

Fig. 5 Dynamical evolution of planets in a less common system that did not undergo any late dynamical instabilities after the gas dispersal. The plots show the temporal evolution (20 Myr) of planetary bodies, specifically their semimajor axes, masses, eccentricities, and orbital inclinations. Each line color represents an individual object; however, the color is random in each plot. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. For better visualization, the panels on the left (semimajor axes and masses) show the evolution of all bodies with masses > 0.1 M⊕, and the panels on the right (eccentricities and orbital inclinations) show only bodies with masses > 0.5 M⊕.

3.2 Effects of disk parameters

In this subsection, we briefly present and compare the various generated systems, and discuss the effects of the disk parameters and initial conditions on the final systems before exploring their architecture more in the next subsection. The simulated planetary systems are displayed in Fig. 8. The architectures of the systems around a 0.6 M_⊙ and 0.8 M_⊙ central star are plotted for different masses of the initial disk, and two (for 0.6 M_⊙ star) and four (for 0.8 M_⊙ star) gas surface density values. In this paper, a “planet” is a body with M_P ≥ 0.5 M_⊕. We used this cutoff to facilitate comparison with observations: at present smaller planets are very rarely detected, and our retrieved observational sample contains only a few of them. In Fig. 8, the symbol size is proportional to the mass of the formed planets. Many of our planets are similar in mass to the close-in super-Earths observed around these stars. Specifically, the simulations with higher initial disk masses reproduce the observed sample quite well (this will be discussed in detail in Sect. 4), whereas the lowest disk mass (5 M_⊕) forms only one or two very low-mass planets (or none at all) barely reaching the “planet cutoff”; hence, they are not displayed in Fig. 8, but will be discussed a bit later. The increasing initial disk mass clearly results in generally more massive planets in a system, particularly in the case of planets very close to the host star. In Fig. 9, the ratios of final masses of the systems M_tot,final to initial disk masses M_tot,initial are plotted against M_tot,initial to show the planet formation efficiency for each simulation. We can see that the simulations with the lowest initial disk mass are not very successful in forming planets. On the other hand, the simulations with higher initial disk masses are quite efficient in producing planets. The ratio M_tot,final/M_tot,initial ranges from approximately 0.6 to 0.84 (when not considering the initial disk mass of 5 M_⊕). In other words, 60 to 84% of the initial mass will end up in the planets. Figure 9 also shows that higher central mass and higher gas density usually produce higher total mass of the simulated systems. This general trend is mostly followed by all the parameter combinations except for 0.8 M_⊙ & 1250 g cm⁻² (darker blue line with circles), which displays somewhat erratic behavior. Particularly, when initial disk mass = 25 M_⊕ the total final mass is considerably lower than expected according to the trend observed in the figure. This can be explained by the fundamentally chaotic nature of the formation process. We discuss this more in Sect. 3.4, where we look closer at the evolutionary track of this simulated system.

In addition, we examined the scaling between disk mass and planet mass in our systems by plotting the average masses of the largest 〈M₁〉 and second-largest 〈M₂〉 planets against the corresponding initial disk masses and their empirical fits (see Fig. 10). This was studied by Kokubo et al. (2006) in a migration-free setting. They showed that the planet mass scales roughly linearly with the available mass in the disk (i.e., $M_{\text{pl}}\propto M_{\text{disk}}^{0.97-1.1}$ for a fixed stellar mass). Raymond et al. (2007) independently found the same result for a star with 0.4 M_⊙. In our simulations, 〈M₁〉 and 〈M₂〉 also typically increase with the disk mass, but predominantly not linearly. Using the least-squares fit method we obtained $M_{\text{pl}}\propto M_{\text{disk}}^{0.6}$ for 0.6 M_⊙ star and $M_{\text{pl}}\propto M_{\text{disk}}^{0.1-1.0}$ for 0. 8 M_⊙ star, so the growth mostly follows a power-law function rather than linear. However, the values are quite scattered and for the 0.8 M_⊙ star, and the fit is significantly different for the largest versus second-largest planets. The comparison with non-migrating simulations shows that migration results in larger variations in planet masses (of the largest planets) and generally a slower increase in planet mass with the initial disk mass. The mass of the most massive planets in the systems will be discussed a bit further in this subsection.

From a theoretical perspective, most models proposed to explain the formation of exoplanets are based on processes that are quite inefficient. Some of the initial solid mass is often lost to the star due to rapid Type I migration. Migration is happening also in our simulations, but the gas-free cavity between the inner cutoff of the simulation range and the inner edge of the gas disk creates a planet trap, which prevents a significant mass loss to the star. The bodies keep piling up around the inner edge of the gas disk until they form planets similar to the close-in super-Earths with orbital distances mostly within ~0.3 AU and masses of several M_⊕ (e.g., Ogihara et al. 2015). Figure 8 shows that many planets migrated farther inward than the inner edge of the gas disk. As mentioned in Sect. 1 when discussing the breaking-the-chains scenario, planets often become locked into resonant chains while migrating. The innermost planet enters the cavity and stops being affected by the disk torques, that is to say, it slowly stops migrating. But if the planet is in resonance with one or several planets still inside the gas disk experiencing the inward Type I migration, the torques induced on those planets might push the inner planet farther inward through the resonance lock (e.g., Terquem & Papaloizou 2007; Brasser et al. 2018). This way one or several planets might end up located very close to the central star (as observed) without colliding with the star. At the same time, planets that formed later or farther away might be prevented from migrating too close, either by joining the resonant chain or if migration simply stops due to the dissipation of the gas disk, so not all planets in a system end up too close to the star (as displayed in Fig. 8).

Tables 3–5 present some of the final characteristics of the simulated systems grouped by the value of the gas surface density and central mass. Average planet mass 〈M_p〉 in an individual system mostly grows systematically with increasing initial disk mass. The same applies to the mass of the most massive planet M_max in each system. In the case of the gas surface density values higher than 1000 g cm⁻², the increase is not that systematic. Both the average mass and the mass of the most massive planet reach a maximum when the disk mass is 30–35 M_⊕, and then often drop. Figure 11 displays this tendency for M_max. Interestingly, in most cases it is actually not the highest initial disk mass that produces the highest M_max. It seems that there is an upper limit for the mass of the most massive planet in a system and that adding more disk mass will not result in a much more massive planet. Generally, both higher central mass and gas density result in more massive planets. At the same time, with more massive planets in a system the dispersion in mass is also higher as all systems contain small planets as well. In addition, it seems that more gas in the disk and higher central mass results in more irregularity for this typically increasing trend. The highest 〈M_p〉 is 4.97 ± 3.58 M_⊕ and the highest M_max is 13.34 M_⊕, so we managed to form some more massive planets than generally considered super-Earths. The average planet in the observed K-dwarf sample (when considering only the super-Earths) has 〈M_p〉 = 5.9 ± 2.5 M_⊕ (planet mass is calculated as an average value from the whole sample, not individually per system). In order to be able to compare the calculated values for the simulated systems presented in Tables 3–5 with the average observed mass, we have to consider observational biases and take into account the fact that 〈M_p〉 values in the table also include small planets located farther from the star than the currently used methods are able to detect (as discussed in Sect. 1). This is done in Sect. 4.

Fig. 6 Dynamical evolution of planets in an extended simulation (100 Myr) that underwent multiple dynamical instabilities much later after the gas dispersal, at approximately 21 and 28 Myr. The plots show the temporal evolution of planetary bodies, specifically their semimajor axes, masses, eccentricities, and orbital inclinations. Each line color represents an individual object; however, the color is random in each plot. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. For better visualization, the panels on the left (semimajor axes and masses) show the evolution of all bodies with masses > 0.1 M⊕, and the panels on the right (eccentricities and orbital inclinations) show only bodies with masses > 0.5 M⊕.

Fig. 7 Histogram showing a number of collisions during the 20 Myr of evolution in all our simulations. Only collisions where both bodies are above the lunar mass (MMoon = 0.012 M⊕) are presented.

Fig. 8 Architecture of our simulated planetary systems around a 0.6 M⊙ and a 0.8 M⊙ K-dwarf star, for two or four different gas surface density values and different initial disk masses (IDM) after 20 Myr of N-body integration using GENGA. Outcomes for the initial disk mass 5 M⊕ are not displayed. The point size indicates the masses in M⊕ of the formed planets. The gray zone is the region inside the inner truncation radius, and the dashed gray line shows the inner edge of the gas disk.

Fig. 9 Mtot,final/Mtot,initial ratio versus the initial solid disk mass, Mtot,initial. The total mass of a system, Mtot,final, is calculated by summing up the masses of all individual planets in each system. Values for 0.6 M⊙ and 0.8 M⊙ star masses and various gas surface densities are presented.

Fig. 10 Average masses of the largest 〈M1〉 (filled circles) and second-largest 〈M2〉 (open circles) planets against the corresponding initial disk masses (IDM) for the two central star masses 0.6 M⊙ (in green, the top plot) and 0.8 M⊙ (in blue, the bottom plot), and the empirical fits for 〈M1〉 (solid line) and 〈M2〉 (dotted line). Since the systems starting with 5 M⊕ typically contain only one planet, they are not included in the plots.

Fig. 11 Mass of the most massive planet, Mmax (M⊕), in a system for each initial disk mass value (IDM) and both star masses, 0.6 and 0.8 M⊙, shown in green and blue, respectively. Mean values and their standard errors are included (green squares with the dotted line and blue circles with the dashed line for 0.6 and 0.8 M⊙, respectively).

3.3 Architecture of the simulated systems

Knowledge of architecture of planetary systems is important for constraining the theories of planet formation and evolution. Systems of planets with nearly coplanar and circular orbits, such as the Solar System, are consistent with the standard model of planet formation via planetesimal accretion in a gaseous protoplanetary disk. While planetary systems with differently inclined, noncircular orbits suggest past events that increased eccentricities and inclinations in the system (e.g., Jurić & Tremaine 2008; Ford & Rasio 2008; Chatterjee et al. 2008), such as resonant encounters between planets, or planet-planet scattering. Figure 8 presents the architecture of the 42 simulated systems after 20 Myr of planetary evolution. We ran 48 different combinations of parameters, but we omit the runs with the initial disk mass of 5 M_⊕. In the end, our simulated population contains 335 planets in 42 systems. We examine their architecture but do not describe each individual system, instead focusing on the general trends. At the end of this subsection, we look specifically at the final eccentricities and inclinations of the simulated planets. All simulations started with planetesimals and planetary embryos distributed between 0.2 and 2 AU. Planets in the final systems are located between 0.06–0.09 AU to 2 AU (just a reminder, during the planet formation if bodies cross the inner radius of 0.05 AU or outer radius of 3 AU, they are taken out of the simulation). Basically, in almost all systems the innermost planet, in rare cases two planets, is found inside the cavity.

The planets can be divided into two loose categories: more massive planets of several M_⊕ located close to the star almost exclusively within 0.3 AU, and lower-mass planets of at most 2 M_⊕, but often much lower mass, which can mostly be found farther from the star (see Fig. 12). We find many lower-mass planets close to the star as well, but no really massive planets with masses above 2 M_⊕ beyond 0.5 AU. Stevens & Gaudi (2013) proposed 2 M_⊕ as the lower limit for super-Earths. This is not a universally accepted definition but if we adopt it, then we basically did not manage to form any super-Earths beyond 0.5 AU. In many cases, there is actually something like a spatial division between the two categories of planets with a region containing no planets (see again Fig. 8). The “zone of more massive planets” typically extends farther out for larger initial disk masses. Generally, in this zone a system contains either a few very massive planets that are far from each other or several planets with lower masses packed more tightly together. This seems to be independent of the gas surface density, central mass or initial disk mass. Relative spacing between the planets is then approximately the same, although the lowest values of disk masses show generally a bit larger distances between their less massive planets (the dynamical spacing between planets is examined in Sect. 3.7).

In the simulated population, we have four systems with a single planet. This is an interesting result as single planet systems are quite common around K-dwarf stars, and not only around them. The large number of Kepler systems with single transiting planets versus multiple transiting planets is known the as Kepler dichotomy (Johansen et al. 2012). Of course, these systems might just be incomplete, even though some studies claim that the dichotomy is real (e.g., Fang & Margot 2012; Johansen et al. 2012; Moriarty & Ballard 2016). On the other hand, other studies demonstrate that super-Earth systems are naturally multiple and the Kepler dichotomy is just an observational effect, a consequence of the mutual inclinations becoming excited due to the dynamical instabilities the systems experience as they evolve (e.g., Izidoro et al. 2017, 2021). These studies predict an insignificant number of real single-planet systems in the Kepler sample, which agrees with the outcomes of our simulations. However, our single-planet systems are formed by the lowest initial disk masses (5 M_⊕) and have very low masses (0.7 M_⊕ at most), which suggests that their initial disks may not have contained enough material to form multiple planets. The dynamical evolution of one these simulations is presented in Fig. 13, and it is representative of all the simulations that end with one planet. As after 20 Myr of evolution there were still many planetesimals available in the simulation, we extended the simulation time to 30 Myr, but not longer due to the extremely long computation time necessary for running the simulations starting with lowest disk masses (i.e., a larger fraction of smaller bodies). In this system, the evolution follows the general trend with embryos migrating inward, but there is not much growing happening. Most of the closest (to the central star) and heaviest embryos do not experience any collisions after 4 Myr. Embryos orbiting farther from the star keep growing until the end of the simulation time, though not at a significant rate. Basically, the embryos grew too slow and then it became too late, because the planet growth phase mostly stopped. Nevertheless, the material is still available in the disk. The simulation started with a disk mass of 5 M_⊕, and during the 30 Myr of planet formation, only ~0.4 M_⊕ is ejected from the system. The rest of the initial material remains in the system as leftover embryos and planetesimals, mostly contained in bodies of the mass of Mars and above. At the end of the simulation, several of the bodies are piled up close to the inner edge, but only one of them actually classifies as a planet in our study, with a mass of ~0.6 M_⊕. A longer simulation time would probably not improve the chances of additional planet growing much; however, a higher surface density or a longer gas disk decay time could increase the number of planets in the system as it would increase the planetesimal accretion or make the planet growth phase longer. Of course, if we defined a planet differently, then we would get a different number of planets in the system (and in all the other “single-planet” systems). It seems that not enough initial building material to form larger bodies causes that all the formed planetary bodies tend to be small and similarly sized, not even reaching the mass of Earth. Our results suggest that the large number of detected single-planet systems may simply be an observational effect as these relatively low-mass planets should have even lower-mass companions.

We also explore the number of planets in the individual systems. The majority of our systems have a relatively high number of planets (up to 12 in one case) because they have a high number of small planets farther from the star. Extending the simulation time would slightly reduce their number in some cases as collisions would still occasionally occur, as discussed in Sect. 3.1. Since the main purpose of this study is to reproduce the observed planet population around K dwarfs, we mostly focus on the region close to the star. The typical number of planets in a system (when not considering systems with the initial disk mass of 5 M_⊕) is between 7 and 9. For a 0.6 M_⊙ star the average number of planets is 7.57 ± 1.81 and 7.43 ± 1.72 for 1000 and 1250 g cm⁻², respectively. For a 0.8 M_⊙ star the average number of planets is 8.86 ± 1.07, 8.71 ± 2.14, 7.57 ± 1.99 and 6.71 ± 1.98 for 1000, 1250, 1500 and 1750 g cm⁻², respectively. It appears that the higher gas surface density the lower the final number of planets, which is the same trend as found in a non-migration case of similar N-body simulations (Kokubo et al. 2006; Raymond et al. 2007), and in the case of the higher star mass for the same gas surface density value more planets form (or survive). Generally, systems with massive planets harbor fewer planets than systems with less massive planets, and these fewer massive planets are usually spaced farther apart.

Figure 14 displays the final eccentricities and inclinations of planets around 0.6 M_⊙ star (in green) and 0.8 M_⊙ star (in blue) plotted against their semimajor axes. Current eccentricities and inclinations (with respect to the invariable plane) of the Solar System terrestrial planets are displayed for comparison. Essentially, all the planets in our sample are within the range of the Solar System planets for both orbital elements (except for a few exceptions). As Mercury’s orbit has the highest eccentricity and inclination of all the eight planets, we can use it as the upper limit. The 0. 8 M_⊙ star plots (the bottom row) display several very massive planets with quite high inclinations and, particularly, eccentricities (compared to the Earth and Venus), which suggest more violent evolution of their planetary systems with possibly recent events of dynamical instabilities. We do not want to investigate these parameters in detail, only show that they have reasonable values based on the studies of similar known exoplanets and the Solar System values. Median value for eccentricities in the simulated population is 0.03 ± 0.03 with a maximum of 0.19, and for the observed sample around K dwarfs it is ${0.08}_{-0.01}^{+0.08}$ (1σ ranges), but the observed sample contains only a small number of planets with calculated eccentricity. Nevertheless, studies show that most Kepler planets tend to have relatively low eccentricities (e.g., Fabrycky et al. 2014; Hadden & Lithwick 2014; Mills et al. 2019; Van Eylen et al. 2019). In addition, for observed planets the inclination is defined as the angle of the plane of the orbit relative to the plane perpendicular to the line-of-sight from Earth to the object, and this makes the comparison with the simulated population a bit difficult. Regardless, as discussed in Sect. 1, the currently known population of planets around similar stars generally have low orbital inclinations, mostly less than 3 deg relative to a common reference plane (Fang & Margot 2012). Median value for inclinations in the simulated population is 1.11 ± 1.38 deg with a maximum of 6.57 (if we do not consider the single outlier inclination in the top right plot with a value of 9.58 deg), which agrees with the study. When comparing values of these parameters specifically for different central masses and gas surface densities, we see that both eccentricities and tend to be a bit lower for the runs with initially higher gas surface density and also for higher central mass, which at least in the case of gas surface density can be explained by the fact that more gas in the disk should result in more eccentricity and inclination damping (as discussed in Sect. 1). Inclinations show this trend as well, but not as systematically. Either way, orbits of the simulated planets are similar to the coplanar and circular orbits of the Solar System planets. This is consistent with planets forming in a protoplanetary disk, followed by evolution mostly without significant or lasting perturbations from other bodies. Figure 14 also shows how much closer to the star the majority of the simulated planets orbit, compared to, for example, Mercury; only lower-mass planets are located farther out.

Fig. 12 Mass-distance relationship for the planetary systems. Magenta and blue circles respectively represent our simulated sample and the observed sample around K dwarfs. The gray zone is the region inside the inner truncation radius, and the dashed gray line shows the inner edge of the gas disk.

Fig. 13 Dynamical evolution of one of the single planet systems. The plots show the temporal evolution (30 Myr) of planetary bodies, specifically their semimajor axes and masses. Only bodies with masses > 0.1 M⊕ are displayed. Each line color represents an individual object; the same object is indicated by the same color in both plots. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. Only one object (in light purple) has a mass above 0. 5 M⊕ and is classified as a planet in this study.

Fig. 14 Eccentricities (on the left) and inclinations (on the right) of the simulated planets around a 0.6 M⊙ star (in green; the top row) and a 0.8 M⊙ star (in blue; the bottom row) versus their semimajor axes. Different shades of green and blue denote the different initial values of the gas surface densities at 1 AU, i.e., 1000, 1250, 1500, and 1750 g cm−2; the higher initial density, the darker color of the circles. The point size indicates the masses in M⊕ of the formed planets. Eccentricities and inclinations of terrestrial Solar System planets are displayed for comparison (in red); their sizes are enlarged by a factor of 2 for better visualization.

3.4 Reproducibility of the outcomes

Now we briefly explore the reproducibility of our simulated outcomes. Even though we might expect that simulations with similar initial conditions produce planetary systems with similar characteristics, we have to take into account the chaotic nature of this stage of planetary evolution. We examine an evolutionary track of the simulated system with the starting parameters of 25 M_⊕, 0.8 M_⊙, and 1250 g cm⁻², whose total final mass is much lower than expected according to the general trend observed in Fig. 9. This overall increasing trend of M_tot,final with the initial disk mass is followed relatively well by all the parameter combinations except for 0.8 M_⊙ and 1250 g cm⁻² (darker blue line), which displays quite irregular behavior. Particularly, in the case of initial disk mass = 25 M_⊕ the total final mass is considerably lower than expected; therefore, we focused on this simulation. We ran another three simulations with exactly the same parameters (also planetesimals and planetary embryos have the same sizes and are distributed in the same way) in order to reproduce the first outcome.

Figure 15 displays the temporal total planet mass evolution of these four runs in total. We plot their individual evolutionary tracks against time; the original simulation is 73_1 Sim, represented by the blue line. Immediately we see that this track behaves very differently compared to the additional three runs. The sudden decrease in the total mass at approximately 8.6 Myr of evolution is followed by a gradual increase, but the final mass of the system is still much lower than for the other simulations, where the final masses are very similar. Each of their masses would actually neatly followed the increasing trend in Fig. 9. After investigating GENGA output files, it is clear that this 73_1 Sim feature in Fig. 15 is caused by a sudden orbital instability that happens when several relatively massive planets, piled-up very tightly close to the inner edge of the disk and locked in a resonance chain, are suddenly disrupted by another planet that gets too close and excites the group. The resonant chain breaks and two planets (with ~1 M_⊕ and ~3 M_⊕) fall into the star by venturing closer than 0.05 AU from the central mass (see Fig. 16). Planets that form later and migrate toward the inner edge might get caught up in resonance with the planets already anchored there, and pump up their eccentricities. Since the innermost planets have no eccentricity damping, it is expected that their eccentricities grow large and they merge. This is what happens in this simulation: the planets merge, but the excitation causes two of the planets to be lost to the star. Some of the smaller bodies then continue colliding and growing for some time, which is represented by the step-wise increase in the total mass following the sudden drop, but too much mass is already lost. The other three simulations do not experience such an instability (when so much mass is removed from the system) and show regular increase in the total mass all the way. This is in line with the chaotic nature of the accretion process. As 73_1 Sim behaves differently than the other three simulations, we treated it as a special case and used 73_3 Sim in the final statistics and comparison of the simulations to the observations. Overall, we see that most significant collisions and merging happen within the first ~ 11 Myr in all four runs, but some small events happen later as well. The various outcomes of the simulations are displayed in Fig. 17. It is obvious that simulations with the same initial conditions and parameters can produce planetary systems with differing characteristics. In two cases, the final systems include fewer massive planets in the region close to the inner edge of the disk; in two other cases, a higher number of smaller planets is located in this region. These two types of systems can still have the same total mass. The erratic behavior makes is difficult to decide which initial conditions and other parameters might be responsible for the final architectures of the formed planetary systems. On the other hand, random variations between the outcomes of successive simulations can help explain the diversity of exoplanet systems, as discussed in detail in Raymond et al. (2020).

Fig. 15 Temporal evolution (20 Myr) of the total planet mass (sum of all planetary bodies with MP > 0.5 M⊕ in a system) in the simulated system called 73 Sim (initial disk mass = 25 M⊕, central mass = 0.8 M⊙, gas surface density = 1250 g cm−2). Shown are the evolutionary tracks of four simulations with exactly the same parameters. The total mass included in Fig. 9 is the final state of 73_1 Sim, represented by the blue line.

Fig. 16 Snapshots of the 73_1 Sim simulation showing eccentricity versus the semimajor axis of the planetary bodies. Panel A shows simulation at ~8.6 Myr of evolution, and each subsequent panel displays the simulation a few thousand years later (B–F). Snapshot A presents a resonant chain of planets located very close to the star, which gets disrupted by another planet that gets too close and excites the group (B). This dynamical instability results in a ~1 M⊕ planet at first (C), and later on a ~3 M⊕ planet (F) falling onto the central star.

Fig. 17 Outcomes of four 73 Sim simulations with initial disk mass = 25 M⊕, central mass = 0.8 M⊙, and gas surface density = 1250 g cm−2. The original simulation is 73_1 Sim. The size of the circles indicates the mass of the formed planets.

3.5 Time needed to grow an Earth-mass planet

In this subsection, we examine the time needed to grow an Earth-mass planet in our simulations. Our findings are presented in Fig. 18, which displays the mass of the most massive body in each system during the 20 Myr of simulation time. The runs are grouped by their initial disk mass and it is clear that higher disk mass generally results in shorter time necessary for growing a planet with the mass of Earth (or higher). Figure 19 then shows the actual time it takes to grow an Earth-mass planet depending on the initial disk mass for the individual planetary systems (circles) and for the average values of each disk mass bin (purple squares with error bars). The light blue and dark blue circles represent the times from the simulations starting with the gas surface density = 1000 g cm⁻² and 1750 g cm⁻², respectively. This clearly shows that the time is typically shorter for the higher gas surface density. The average time ranges from less than 0.1 Myr for the initial disk mass of 40 M_⊕ to around 4 Myr for the disk mass of 10 M_⊕, and displays the generally decreasing trend with the increasing initial disk mass. This is expected as higher surface density of solids as well as larger initial bodies lead to the planetary growth happening faster as discussed already in Sect. 2.1. Embryos grow faster because the number of collisions required to reach a certain mass is lower when the simulations start with larger planetesimals. We start all our simulations with the same number of particle (i.e., the same resolution), but since we increase the initial mass of the disk, we ran the simulations with larger and larger planetesimals. This causes the simulations with higher initial mass of the disk running much faster than simulations with lower disk masses. Even though, there is no commonly accepted formation time of an Earth-mass planet, the shortest times (~0.1 Myr and less) are probably not realistic. The simulations starting with the disk mass 10 M_⊕ show the largest dispersion of the time needed to grow an Earth-mass planet from ~400 000 yr to even close to 10 Myr in two cases. However, the vast majority of the simulations form their Earth-mass planet(s) within the lifetime of the gas disk, and the planetary bodies then often keep colliding and growing.

In Fig. 18, we also see that most collisions happen within a few million years (~2 Myr), which is in agreement with other similar studies (e.g., Ogihara et al. 2015; Zawadzki et al. 2021). In some runs, occasional collisions occur even after 10 Myr. Generally in our simulations, ~80 to 95% of collisions happen during the first 2 Myr, and only ~3 to less than 1% occur after the first 10 Myr (for the runs starting with 10 M_⊕ and 40 M_⊕, respectively). So, particularly for the higher initial disk masses, almost all collisions happen during the first 10 Myr of the evolution. According to the core accretion scenario (Perri & Cameron 1974; Mizuno et al. 1978) a protoplanet starts efficiently capturing a massive gas envelope from the protoplanetary disk to become a gas giant when its mass is ~10 M_⊕ (e.g., Mizuno et al. 1978; Stevenson 1982; Bodenheimer & Pollack 1986; Hubickyj et al. 2005). In our simulated sample, there are only about eight planets that actually reach the mass necessary to get into this runaway gas accretion phase, and this happens in all cases after more than ~6 Myr of evolution when most of the gas is already gone. This could explain the fact that giant planets are quite rare around this kind of star, as discussed in Sect. 1. Smaller planets of only several M_⊕ are not assumed to accrete substantial amounts of gas.

3.6 Period ratio distribution and logarithmic spacing

Now, we discuss the simulated systems and further compare them to the observed systems. The primary goal of this study is to reproduce the known exoplanets around K-dwarf stars with their characteristics and by exploring how the initial conditions of the late stage of planet formation affect the properties of the simulated planetary systems. To evaluate the results we used the sample of the observed K-dwarf systems described in Sect. 1 for comparison. The sample is actually quite small; therefore, in some cases we extend the sample by known systems with at least two planets with known masses discovered mostly by the Kepler telescope using transit method around M-, K-, and G-dwarf stars. At the present time, transiting exoplanets provide the most informative and complete collection of data on exo-planet characteristics, but only within ~0.2 AU of their host stars. This sample (hereafter the MKG sample) contains 384 planets from 143 planetary systems and it was also retrieved from NASA Exoplanet Archive⁴. In this section, we compare the results of our simulations directly to the observed samples. However, since all observational data are affected by observational biases and most of the known exoplanets were discovered by Kepler mission, we apply the detection bias of this survey to our simulated data and then compare them to the Kepler observations again in the following section.

First we examine period ratios and spacing between planets in the simulated population. These are two of the most important characteristics of planetary systems. We compute period ratios P_out/P_in of all subsequent pairs of planets in the individual systems for our simulated sample and compare them to the period ratios of the MKG pairs. Figure 20 shows that many planetary pairs in both samples are located in or near orbital resonances, indicated by the dashed lines. Orbital migration often moves planets into resonances (Terquem & Papaloizou 2007), which can then either stabilize or destabilize a system if the resonant chain becomes too long (Matsumoto et al. 2012; Goldberg et al. 2022). We see that the MKG sample has a quite broad distribution of period ratios quantitatively similar to our sample that is, however, a little narrower. The regions of the most occupied ratios partially overlap, but in the simulated population we do not see the decrease in occurrences at around 1.8, which is quite significant in the MKG sample. In both samples, the majority of the pairs are in 3:2 resonance, but this peak is much more prominent in our simulated sample compared to the MKG systems. Many exoplanets have been found to be in resonances, specifically, 3:2 resonance has been confirmed to be well populated by multiple studies (Lithwick & Wu 2012; Batygin & Morbidelli 2013; Fabrycky et al. 2014). In our sample, the 5:3, 4:3, and much less pronounced 2:1 resonances, and in the observed sample 2:1 and 5:4/4:3 resonances, are also very common. A weaker peak near the 2:1 is present in the Kepler distribution according to Lissauer et al. (2011) as well. Other studies of the Kepler sample also confirmed weaker peaks near the 5:3, 4:3 and 5:4 resonances (e.g., Aschwanden & Scholkmann 2017). Nevertheless, as already discussed in Sect. 1, most currently known super-Earths are not found in resonant systems (Lissauer et al. 2011; Fabrycky et al. 2014). They spend time trapped in resonances, but most are removed by late instabilities according to breaking-the-chains scenario (e.g., Terquem & Papaloizou 2007; Ogihara & Ida 2009; McNeil & Nelson 2010; Cossou et al. 2014; Izidoro et al. 2017, 2021). In our simulations, almost all systems go through one or several instabilities, but the instabilities often disrupt only part of the chain, or the planets manage to readjust back to a resonant chain after the instability. At 20 Myr, around 75% of the planet pairs in our systems are not in resonances. Also, around 75% of the systems do not contain multiple-planet resonant chains, even though some of them still contain individual pairs in resonances. As shown in Fig. 4, even a system that underwent several instabilities may have some planets in resonances. If we define systems with surviving resonant chains as the stable ones, then we get around 75% of unstable systems. Previously we showed that 85% of our systems potentially underwent at least one late instability (Sect. 3.1). These numbers are consistent with the estimates in Izidoro et al. (2017), where they presented that at least 75% of their simulated systems must be unstable to match the observations; and probably even 90–95%. We also showed that extending the simulation time would probably increase the numbers (Sect. 3.1).

Planetary systems generally follow a simple logarithmic spacing rule (e.g., Bovaird & Lineweaver 2013; Mousavi-Sadr et al. 2021). The relation can be written as follows: log P_n = c₁ + c₂(n − 1), where n is the planet number from the inside out, and c₁ and c₂ are (fitting) constants, where c₁ = log P₁ (i.e., the log of the period of the innermost planet). We compute logarithmic spacing between pairs of neighboring planets in each planetary system (see Fig. 21) for both the simulated and the observed systems, and plot the obtained values of c₂ versus c₁. The figure shows our simulated sample together with the sample around MKG dwarfs. The observed sample around K dwarfs is included in the calculations, but is not displayed as it looks similar to the MKG sample (with fewer planets though). For the simulated systems, c₁ ranges from 0.8 to about 1.3 and the typical value is close to 1 (mean = 0.98 ± 0.10) − that is, P₁ ~ 9.5 days - while c₂ ranges from 0.1 to 0.5 and is typically around 0.24 (mean = 0.24 ± 0.08) − that is, $P_{\text{out}}/P_{\text{in}}={10}^{c_2}\sim 1.7$. For known systems around K and MKG dwarfs, the values of c₁ and c₂ show a large range: for the K-dwarf sample, the mean of c₁ is 0.68 ± 0.42 and of c₂ it is 0.27 ± 0.33, and for the MKG sample, the mean of c₁ is 0.75 ± 0.46, and of c₂ it is 0.26 ± 0.33. This means that the typical period of the innermost planet is ~5 days and ~6 days, and the typical P_out/P_in ratio is ~1.9 and ~1.8 for the K and MKG samples, respectively, whereas P₁ ~ 9.5 days and P_out/P_in ~ 1.7 for the simulated population. We note that the innermost planets in the simulated systems are farther from their host stars than in both observed samples, which results in generally higher values of c₁. This can be explained by the fact that the simulations are limited by the inner truncation radius. Many of the known planets are located closer to the star than the inner truncation radius of our simulations. Also, there is a small difference in the period ratios: the orbits of two adjacent planets are generally a bit closer to each other in the simulated systems compared to the known systems (P_out/P_in ~ 1.7 versus ~1.8/1.9). Small planets are usually packed tighter than larger planets and the simulated planets are generally less massive than the planets in the observed samples, so this behavior might be expected. In addition, some of the currently known systems are possibly incomplete and these “missing” planets can affect the typical period ratio of the whole system. Nevertheless, the difference is so small that we did not explore it more. We also clearly see (in Fig. 21) that the simulated systems are much more similar to each other than the observed ones; the latter are very diverse. This is also expected since the initial conditions and parameters of the simulations are limited to a several specific values; therefore, they might not be able to produce very diverse planetary systems.

Fig. 18 Temporal evolution (20 Myr) of the mass of the most massive surviving body in the system. It shows the time it takes to grow an Earth-mass planet for each system (where the solid lines representing each simulation meet the dashed line). The simulations are grouped by their initial disk mass (IDM), indicated by the same color. Bodies in the systems starting with the initial disk mass 5 M⊕ never reach 1 M⊕ and are therefore not included in the plot.

Fig. 19 Time it takes to grow an Earth-mass planet depending on the initial disk mass (IDM) for the individual planetary systems (circles) and for the average values of each mass bin (purple squares with error bars). The dashed purple line indicates the descending trend of needed time with increasing disk mass. The light blue and dark blue circles represent the times from the simulations starting with the gas surface densities = 1000 g cm−2 and 1750 g cm−2, respectively.

Fig. 20 Period ratios Pout/Pin of all pairs of subsequent planets for our simulated population of planets (in magenta) compared to the period ratios of the pairs from our MKG sample (see the hatched area of the histogram). Main orbital resonances between the planet pairs are indicated by the dashed gray lines and identified by the fractions at the top of the plot.

Fig. 21 Logarithmic spacing between pairs of neighboring planets in each planetary system. The constant c2 is plotted against the values of constant c1. Magenta and blue circles represent our simulated sample and the observed sample around MKG-dwarf stars, respectively.

3.7 Long-term stability of the simulated systems

To assess a long-term orbital stability of our simulated objects in their positions in multi-planet systems without resorting to computationally expensive simulations, usually simple, approximate conditions can be used, such as a minimum separation between neighboring planets in mutual Hill radii or dynamical spacing Δ (Chambers et al. 1996; Pu & Wu 2015). Following the Hill stability criterion (Chambers et al. 1996), which defines the critical value of the mutual separation as $2\sqrt{3}$ Hill radii, we can access the Hill-stability of the systems. Planetary orbits in Hill-stable systems should be unable to cross each other. None of the systems, neither observed nor simulated, contains a pair of planets with a mutual distance below this critical value (see Fig. 22), so none are supposedly Hill-unstable. Still, in any system, two adjacent planets are less likely to be long-term (gigayear-scale) stable when their dynamical spacing is smaller than Δ = 10 (Pu & Wu 2015). Only ~5% (13 out of 288) of all simulated pairs are less likely to be stable compared to almost a quarter of the K-dwarf systems and almost a third of the MKG systems. Yet, if they are long-term stable, these planets are likely to be in orbital resonance with each other, as is the case for many exoplanetary systems. In general, the notion of Hill stability should be interpreted with caution. For example, Jupiter and Saturn have dynamical spacing Δ ~ 8 and still they are very dynamically stable as a pair (Laskar 2008).

Figure 22 displays the dynamical spacing Δ computed between all pairs of adjacent planets in each system. Planets in Kepler multi-planet systems tend to have a spacing of around 20 Hill radii, and are usually not closer together than 10 Hill radii (e.g., Lissauer et al. 2011; Fabrycky et al. 2014; Weiss et al. 2018). For the simulated systems, only ~5% of the adjacent pairs have Δ < 10 and the typical Δ value is ~21 (median = 20.9 ± 15.9), in agreement with previous studies. For the observed K-dwarf sample, ~24% of the adjacent pairs have Δ < 10 and the typical value is ~15 (median = 15.2 ± 9.0), and for the observed MKG-dwarf systems, ~31% of the adjacent pairs have Δ < 10 and the typical value is ~13 (median = 12.8 ± 12.1). This is similar to the findings of Fang & Margot (2012). Even so, the simulated systems tend to be wider spaced than the observed ones. The dispersion is very large for all three typical Δ. This is because all three populations include neighboring pairs of very massive planets orbiting very close to each other, as well as low-mass planets located far away from each other. Again, the possibility that at least some of the currently known systems might be incomplete can also play a role here. Additionally, if we focus only on the zone of more massive planets within 0.3 AU of the star (which is where basically all observed planets around K dwarfs are found), then we get a lower typical value of Δ ~ 15 (median = 15.2 ± 5.6) for the simulated systems. This value agrees very well with Δ ~ 15 from the K-dwarf sample. Planetary orbits are usually spaced farther apart as distance from the star increases, our results support this as well. As the observed samples contain mostly only planets very close to the star while the simulated systems contain both very close but also distant planets, this can explain the generally larger dynamical separation of the complete simulated population.

The Hill stability criterion is intended for two-planet cases, not for multi-planet systems, as it does not consider mutual inclinations between orbits. Mutual orbital inclinations are known to influence the evolution of planetary systems and their long-term stability strongly. Therefore, we also use the angular momentum deficit (AMD; Laskar 1997, 2000; Laskar & Petit 2017; Petit et al. 2017) to predict the long-term stability of our simulated systems. The AMD is a measure of the deviation of the system from being perfectly circular and coplanar. This is a more sophisticated approach that accounts for both eccentricity and mutual inclination, and it is applicable to multi-planet systems. A planetary system is AMD-stable if the AMD in the system is not sufficient to allow for collisions between its planets. We used the AMD stability criterion (Laskar & Petit 2017) for 42 simulations (systems starting with a disk mass of 5 M_⊕ contain only one or two very small planets and are therefore not included), and in Fig. 23 we plot the values of the AMD stability coefficient,β (in log scale), which is defined as $$\beta =\frac{C}{{\text{Λ}}^{\prime }{C}_c},$$(2)

where C is the total AMD of a system, Λ′ is the circular momentum of the outer planet and C_c is the critical AMD, such that for smaller AMD, collisions are forbidden. A pair is considered as AMD-stable if its AMD stability coefficient β < 1 or logβ < 0; thus, collisions are not possible. A system is AMD-stable if every adjacent pair of planets is AMD-stable. In the case of the innermost planet, the pair is formed with the star. The majority of the simulated systems (~55%) meets the criterion for AMD stability at the end of 20 Myr, which means that these systems quite quickly settled into a long-term orbital arrangement. The systems that are AMD-unstable are often very close to stability with only one or two pairs with logβ positive, but very close to 0. The closest pair (the innermost planet and the star) is in all cases very AMD-stable; hence, we do not expect the planet to be in danger of colliding with the star, unless the AMD is transferred inward from a more distant AMD-unstable pair. However, AMD-unstable systems are not necessarily unstable. They might be stabilized by the presence of MMRs, which can prevent pairs from colliding or prevent tidal evolution. In fact, around 80% of all pairs in the AMD-unstable systems are in or very near MMR. It seems that our simulations have formed mostly long-term stable planetary systems already after 20 Myr. Additionally, we extended the simulation time for some of the AMD unstable systems; one-third of all runs was extended to 40 Myr and one-fifth to 100 Myr. At 40 Myr, we see changes in AMD stability coefficient values in every system. AMD has been transferred between the planet pairs, so that some of them have become more AMD-stable and some less. At 100 Myr, we see almost no changes in AMD of the systems compared to the values at 40 Myr. Generally, we do not see much improvement in the AMD stability of the systems after extending the simulation time.

Fig. 22 Dynamical spacing, Δ, between all pairs of adjacent planets in each system in units of mutual Hill radius plotted against their period ratios, Pout/Pin. The figure shows our simulated sample (in magenta) together with the observed sample around MKG dwarfs (in blue). Two important values of dynamical spacing are represented by the dashed lines: the Hill stability criterion $\Delta =2\sqrt{3}$ (Chambers et al. 1996) in orange and Δ = 10 in gray.

3.8 Peas in a pod

Exoplanets within a system tend to be of similar sizes and masses, often ordered in size and/or mass and evenly spaced. In addition, smaller planets are frequently packed in tight configurations, while large planets often have wider orbital spacing (relative spacing between the planets then being approximately the same). This “peas in a pod” pattern, proposed for the first time by Weiss et al. (2018), Millholland et al. (2017), is not universally accepted but recent observations indicate the presence of several of these correlations in the architecture of exoplanetary systems (Millholland & Winn 2021; Otegi et al. 2022). Some correlations have also been reproduced in planet formation simulations (Mishra et al. 2021; Mulders et al. 2020). Since our simulations do not allow us to assess radii of the final bodies, we examine masses of adjacent planets in the simulated systems and their spacing to evaluate whether our data follow at least some of these similarity trends. Figure 24 shows the mass ratios of all outer/inner neighboring planets (the top plot) and the period ratios of the outer/inner planet pairs (the bottom plot). According to the peas in a pod pattern, M_i+1/M_i ≈ 1 and (P_i+2/P_i+1)/(P_i+1/P_i) ≈ 1 are expected for these ratios. For mass ratios, we analyzed all 288 planet pairs in all simulated systems that contain more than one planet, and using the Pearson correlation test, we calculated that there is a moderate positive correlation for masses of adjacent planets with the R-value of 0.58 and P-value of 4.27 × 10⁻²⁷ for all planets in the data set. For period ratios, we have found no possible correlation for the entire sample with R-value of 0.04 and P-value of 0.53. This is likely due to the spatial division between the two categories of planets, the close-in more massive planets and the less massive distant planets with a region often containing no planets in the middle (as discussed in Sect. 3.3). However, if we consider only the more massive planets within 0.5 AU of the star (i.e., 127 planet pairs), then we again find a moderate positive correlation for periods of adjacent planet pairs with R-value of 0.54 and P-value of 4.48 × 10⁻¹¹. Thus, we have shown that the architectures of our simulated systems often follow the peas in a pod pattern, at least when it comes to planetary mass and period ratios (only within 0.5 AU) of adjacent planets/planet pairs. Mamonova et al. (2023) analyzed an observed sample containing planets around M- and K-dwarf stars, and found a strong correlation of the R-value of 0.896 and P-value of 5.6 × 10⁻¹⁵ in mass ratios but no correlation in period ratios. This agrees with our findings as we find no correlation in period ratios, if we do not restrict the population.

Fig. 23 AMD stability of the simulated systems. Each planet is represented by a circle, the size of which is proportional to the mass of the planet. The color represents the AMD stability coefficient, β (in log scale), of the inner pair associated with the planet. The innermost planet is represented by the AMD stability coefficient associated with the star. Numbers in the simulation names (y-axis) are random, but the systems are ordered by their initial disk mass value, which decreases from top to bottom. Systems starting with the disk mass = 5 M⊕ contain only one or two very small planets and are therefore not included in the plot (but they are AMD-stable).

Fig. 24 Mass ratios of the outer/inner neighboring planet of the entire sample (the left plot) and period ratios of the outer/inner planet pair within 0.5 AU of the star (the right plot) in the simulated population. The dashed gray line is the 1:1 line.

4 Discussion

In this section, we present the final comparison of the simulated and observed populations and analyze whether our effort to reproduce the known exoplanet sample around K-dwarf stars was successful and to what extent, and which initial conditions seem to be the most favorable. Limitations of the study and potential future work are discussed at the end.

4.1 Comparison to the known systems

The final visual comparison of the simulated versus observed planet population around K dwarfs is shown in Fig. 25. Our focus is mainly on the known planets with a super-Earth mass (~ 10 M_⊕) or less, although the simulations produced also several planets with slightly higher masses. There is a significant region where the two samples overlap, which shows that we managed to reproduce the observed population at least partially. The fact that simulations are limited by the inner truncation radius can explain why there are no simulated planets with semimajor axes below 0.05 AU. Our simulated systems also contain many small planets mostly at larger distances, which are basically “undetectable” by the current methods (as discussed in Sect. 1). These planets are located in the bottom right corner of the figure, where no observed planets can be, and has not been found.

Our simulated sample clearly differs from the Solar System planets and contains many closely-packed hot super-Earths; this is expected based on the observed sample and from our initial conditions, and it also confirms the hypothesis that the final systems tend to be compact with short orbital periods (as discussed in Sect. 1). According to Migaszewski (2016), if a system ends up in a compact configuration, it might be attributed to multiple chains of MMRs between neighboring planets, arising from the planet-disk interactions causing inward migration at early evolutionary stages of planet formation. This is something what we observe in our simulations. The simulated sample does not consist of a large number of systems (42), but still a basic comparison to the observations is possible. In order to ensure that we make meaningful comparisons, we perform our statistics in several ways. We focus only on planets closer than 0.3 AU, as the K-dwarf sample basically does not include planets located farther from the star, except for a few giant planets (Sect. 1). We also focus only on the super-Earths, so that the larger observed planets are not considered. Additionally, we omit all the planets with masses below 2 M_⊕ as the known population around K-dwarf stars contains almost no planets with such low masses. At last, we need to perform a slight modification to our data set. Due to the limitation of our model, our population is farther from the central star than the known population. To be able to assess whether we managed to reproduce the observations in term of semimajor axes, we need to eliminate this issue. As the inner truncation radius of the simulations is located at 0.05 AU from the star, we move all the simulated planets closer to the star by this distance. Finally, we recalculate the values from Tables 3–5 based on the new cutoffs and present the updated values in Tables 6–8. As mentioned earlier, the average planet mass in the observed K-dwarf sample (when considering only the super-Earths, not the giant planets) is 〈M_p〉 = 5.9 + 2.5 M_⊕. When analyzing our data, specifically 〈M_p〉 and M_max, we see that the simulations starting with < 15 M_⊕ do not reach the masses necessary to reproduce the average mass of the observed planets. It seems that the model needs a higher initial disk mass. On the other hand, as already shown in Fig. 11, the most massive planets are generated with the disk mass of 30 M_⊕ and above, but typically not with the highest disk mass. It appears as if the chosen initial disk mass values cover the whole range: the lowest do not reproduce the expected outcome, whereas the highest disk masses might be a bit too high. The typical number of planets in a system (when not considering disk mass = 5 M_⊕ systems) is now 3.1 ± 1.22, which agrees very well with the value for the observed sample around K dwarfs, which is 3.0 ± 1.3.

The chaotic nature of the formation process presents a challenge in comparing the outcomes of the simulations to the known population. It also makes it difficult to determine what initial conditions best correspond to reality as similar conditions can produce quite different final architectures of the systems (already shown in Sect. 3.4). Now we describe the model results in the context of our target expectations and variability of the parameter space we explored, and explain why not all models “match” the observed systems and examine how well they reproduce general trends in the currently known population. Investigating the cumulative distribution functions (see Fig. 26) of the semimajor axes in the simulated systems (the top row), we see that the cumulative distribution of the whole population (on the left) is very close to the observed samples. We used the two-sample Kolmogorov Smirnov (K S) test to determine whether the two sets of samples (simulated + K-dwarf sample and simulated + MKG-dwarf sample) come from the same distribution. The p-values for the respective pairs of samples are presented in Tables A.1 and A.2. In the case of the whole simulated sample versus K-dwarf sample and MKG-dwarf sample, we find no significant difference between the distributions for either of the two sets. The simulated distribution actually follows the MKG sample better than the K-dwarf sample. This just might be due to the fact that the K-dwarf sample is smaller. Distributions for the different gas density values (the top middle plot) are also very close to the observed samples. The K-S test showing no significant difference between the distributions. The simulated distribution with 1250+ g cm⁻² (i.e., simulations starting with the gas density 1250 g cm⁻² and above) is actually closest to the K-dwarf sample. Distributions for the different initial disk masses (the top right plot; IDM = 10+ means that only the systems starting with the disk mass of 10 M_⊕ or more are considered) are still mostly very close to the observed samples with K–S test showing no significant differences, but we see that the higher disk masses, particularly 40 M_⊕, are not that similar anymore based on the commonly used 5% significance level. This might be partially caused by the small size of this particular sample. The simulated distribution for the disk masses of 10+ M_⊕ (the same as for the whole simulated sample) is the most similar to the K-dwarf sample. When we examine the cumulative distributions for the planetary masses in the simulated systems (the bottom row), we again look at the whole simulated sample at first (the bottom left plot). Here, the simulated distribution looks very different compared to the observed ones and K–S test actually rejects that the two data sets are coming from the same distribution. This suggests that for these parameters our model is incapable to reproduce the observed population. The other two plots for the masses look more promising, as with the increasing gas surface density and initial disk mass (or more precisely, with removing their lowest values) the simulated distributions are getting closer to the observed distributions. Specifically, the distributions for the different gas densities (the bottom middle plot) show that the curves for the higher values actually resemble the observed distributions much more. Both samples 1500+ and 1750 g cm⁻² possibly come from the same distribution as the K-dwarf sample. Finally, when we focus on the initial disk masses (the bottom right plot), we see that the higher disk mass values improve the distribution for values up to IDM = 30+ M_⊕, above that higher disk mass brings no improvement. Distributions for masses < 5 M_⊕ are different from each other and are mostly getting closer to the observed ones with the increasing disk mass, while above 5 M_⊕ the cumulative distributions are quite similar for all disk masses, and are different from the observed ones. The distributions are converging, but not to the observations. The cumulative mass distributions for both observed samples increase almost linearly, implying that the probability distribution of the planetary masses is flat; in other words, finding planets with masses between 2 and 12 M_⊕ is all equally probable (uniform distribution). This is not the case for the simulated population. Our cumulative mass distributions show a strong preference for planets with masses <6 M_⊕, and fewer planets with larger masses, suggesting that we are missing a mechanism in our simulations for creating more massive planets. The K–S test also shows that our samples do not come from the same distribution as the K-dwarf population.

For our simulations, we chose two different values of the density for the less massive star with 0.6 M_⊙ and four values for the more massive star with 0.8 M_⊙. As our test-runs showed no major differences between the outcomes from the two different star masses, we limited the 0.6 M_⊙ star simulations to the two density values to save the computation time. Also, we assumed lower gas density for the lower mass star. However, the unknown mechanism for creating more massive planets might actually be the higher gas density, as we actually manage to reproduce the masses of the observed sample with the values of 1500 g cm⁻² and above. Both distributions (Fig. 26, the bottom middle plot) show a trend similar to the uniform distributions of the observed populations. Additionally, extending the simulation running time could further increase the number of more massive planets and reduce the number of less massive planets, and result in a more uniform distribution of the simulated population. Our simulations also produce many low-mass planets with M < 2 M_⊕, both close to the star and farther away, which are not found in the observed sample. This appears to be due to the observational biases as we show in the following Sect. 4.2. Whether these planets are indeed there and not yet discovered, or whether they are not there and our models are insufficient, at the moment we cannot say.

Models presented in this study do not include gas accretion onto the planets or atmosphere formation, which is a current limitation of GENGA (gas accretion is not yet implemented in the code). Even though the simulated planet population had mostly formed before the gas disk dissipated (see Sect. 3.1), the few planets with the potential of becoming gas giants reached the necessary mass to accrete significant atmospheres only after the gas was mostly gone. So since significant gas accretion is not a process that is expected to take place for the planets formed in our simulations, we assume that the gas then disappeared mainly through the photo-evaporation by the stellar radiation.

Fig. 25 Masses versus semimajor axes of planets of our simulated planetary systems (in magenta) compared to the currently known systems around K dwarfs (in blue). Only planets with known masses are presented. Known planets with masses above the super-Earths mass range (above 10 M⊕) are indicated in light gray. Here, we use the classification of planets by mass proposed by Stevens & Gaudi (2013). The point size indicates the masses in M⊕ of the planets. The sample of the known systems around K-dwarf stars was retrieved from https://exoplanetarchive.ipac.caltech.edu/ on December 4, 2022.

Fig. 26 Cumulative distributions of the simulated semimajor axes (the top row) and planetary masses (the bottom row) for the whole population (the left column), for different gas surface densities (the middle column), and for different initial disk masses (IDM, the right column). IDM = 10 + M⊕ means that only systems starting with a disk mass of 10 M⊕ or more are considered for that distribution. The same is true for the other IDMs as well as for the gas surface density values, e.g., 1000+ g cm−2. IDM = 5+ M⊕ contains only a few very small planets compared to IDM = 10 + M⊕; therefore, we omitted that distribution. Cumulative distributions for the K-dwarf (blue dotted line) and MKG (dashed black line) samples are displayed for comparison.

4.2 Synthetic observations and their comparison to the Kepler population

To account for Kepler observational biases, we performed synthetic observations of our systems using the Exoplanet Population Observation Simulator (EPOS; Mulders et al. 2018, 2019), designed to simulate survey observations of synthetic exoplanet populations. The survey detection efficiency is the average detection efficiency of all main-sequence stars of spectral types F, G, and K. Figure 27 presents the synthetic observations of our simulated population and the associated statistics of the systems. The code uses planetary radii for the statistics. As our model does not simulate radii, we used a mass-radius relation based on Chen & Kipping (2016) to calculate them. All statistics are computed for planets with orbital periods between 1 and 400 days and radii between 0.5 and 5 R_⊕. Panel A shows the orbital periods and planet radii of detected (magenta) and undetected (gray) planets. We see that the undetected ones are mostly smaller planets of radii less than ~1 R_⊕ (or a bit above that) and therefore masses less than ~1 M_⊕, and with orbital periods longer than ~100 days or semimajor axes more than ~0.4 AU. These cutoffs are similar to the cutoffs used in Sect. 4.1, and actually exclude the majority of the planets with masses below 2 M_⊕ from our simulated population as most of these low-mass planets are located farther than 0.4 AU from the central star (see again Fig. 12). The statistics in panels B D present the intrinsic distribution of our systems’ properties, the synthetic observable distribution, and the distribution observed with Kepler. Plot B displays the distribution of period ratios between adjacent planet pairs. The observable distribution (in magenta) is slightly moved to the higher values compared to the intrinsic distribution (gray dotted line). Debiasing typically shifts the intrinsic period ratio distribution to larger values because planet pairs with one of the planets at larger orbital periods are less likely to be detected (e.g., Brakensiek & Ragozzine 2016). The opposite is happening here, since the intrinsic distribution is the “debiased” one. The reason why the data exhibit the opposite trend is the fact that many simulated planets at larger distances from the star have small period ratios (see Fig. 28), and since they also have lower detection probabilities, they are removed from the observable population, which causes the shift toward higher ratios in the observable distribution. This was observed and explained in Mulders et al. (2019). The observable distribution follows the Kepler data (red dashed line) a bit better than the intrinsic distribution. Plot C shows the frequency of multiplanet systems, or the number of detectable planets per star. This statistic traces mainly the mutual inclination distribution together with the spacing between planets, their sizes, and orbital periods. Since usually only a few planets in each system are transiting, the frequency presented here (in magenta) does not really reflect the intrinsic multi-planet frequency (gray dotted line) of our population. However, in Sect. 4.1 after applying the cutoffs to the population, we determined that the typical number of planets in our system is 3.1 ± 1.22, which seems to be in good agreement with the data in plot C. It is clear that the frequency does not show Kepler dichotomy (Johansen et al. 2012), the large number of systems with single transiting planets versus multiple transiting planets, which is clearly visible in the Kepler distribution (red dashed line). Plot D displays the distribution of radius ratios of the adjacent planet pairs. The intrinsic distribution (gray dotted line) peaks at values slightly lower than 1, whereas the observable distribution (in magenta) peaks at higher values, slightly above 1, corresponding to somewhat larger outer planets compared to the inner ones in the planet pairs. The observable distribution follows the Kepler data (red dashed line) much better than the intrinsic distribution. We will explore planetary radii more in future studies, when we can actually simulate them in our model. Figure 28 shows, apart from orbital period ratios versus semi-major axis, also the innermost planet locations. The histogram on top shows the marginalized distribution of the innermost planets compared to the distribution derived from Kepler, and the histogram on the right shows the marginalized period ratio distribution compared again to that of Kepler. This plot is generated by EPOS, but it displays unprocessed data produced by our model.

Fig. 27 Synthetic EPOS (Mulders et al. 2018, 2019) observations of our simulated population (A) and the associated statistics of the systems (B–D). Panel A shows the orbital periods and planet radii (calculated using a mass-radius relation from Chen & Kipping 2016) of detected planets (in magenta) and undetected planets (in gray). The statistics in panels B–D show the intrinsic distribution of our systems’ properties with dotted gray lines, the synthetic observable distribution in magenta, and the distribution observed with Kepler in dashed red lines. Plot B displays the distribution of period ratios between adjacent planet pairs. Plot C shows the frequency of multi-planet systems, and plot D the distribution of radius ratios of the adjacent planet pairs. All statistics are calculated for planets with orbital periods between 1 and 400 days and radii between 0.5 and 5 R⊕.

Fig. 28 Orbital period ratios (magenta) and innermost planet locations (orange). This plot was generated using EPOS (Mulders et al. 2018, 2019). The period ratios of adjacent planet pairs are plotted at the semimajor axis of the outer planet. The histogram on the right shows the marginalized period ratio distribution compared to that of Kepler. The semimajor axis of the innermost planet in each system is displayed in orange. The histogram on top shows the marginalized distribution of the innermost planets compared to the distribution derived from Kepler (red line). Both Kepler distributions are derived by Mulders et al. (2018).

4.3 Study limitations

One of the main limitations of this study are the planetary radii not being realistic. As discussed in Sect. 2.1, all the simulated planets have the same density as they formed by accreting planetesimals with the same density, and gravitational compression effects are not taken into account. Therefore, in this paper we do not draw any conclusions from the radii of our population, but we are planning to address them in the future studies. In Sect. 2.1, we also mentioned using a slightly longer time step than is recommended due to extremely long computation time (several months) necessary for our simulations. Here, we compare the results of a run with the recommended time step with a run using the longer time step. Our simulations run with a time step of 0.7305 days. Additionally, we ran one of the simulations with a four-times-shorter time step of 0.182625 days (the numbers are chosen so that we get an integer number for the total simulation time), to show that the outcomes of the simulations do not seem affected by the length of the time step. Figure 29 shows the planetary systems formed after 10 Myr of simulation with the parameters: central mass = 0.8 M_⊙, initial disk mass = 20 M_⊕, gas surface density = 1250 g cm⁻², and a time step of 0.7305 days (simulation 61) and 0.182625 days (simulation 61_st, shorter time step). We see that the architectures of the simulated systems and planet characteristics are similar. Simulations contain almost the same number of planets, 11 and 12, and the total mass of the planets is in both cases very similar, 14.4 and 14.3 M_⊕. We see that the mass ratio of adjacent planets is somewhat different, with simulation 61 containing close-in planets with mostly similar masses, while the simulation 61_st planets are quite different in mass. However, this is a common feature, which we can see in other simulations as well, either in Sect. 3.2 or in the various outcomes of 73 Sim simulations in Sect. 3.4. We used the K–S test to compare these two systems and get a p-value of 0.94 and 0.99 for masses and semimajor axes, respectively, which confirms that they are sampled from the same distribution. The simulation with the shorter time step needs a very long computation time (4 times shorter time step results in approximately 4 times longer execution time); therefore, we compared results after 10 Myr, even though in the rest of the paper we use the results after 20 Myr of simulation for simulation 61 (with the longer time step).

Another limitation is connected to the length of the time step, and it is the inner truncation radius set at 0.05 AU. Many of the known exoplanets are located closer to the star than the inner truncation radius of our runs; therefore, in order to reproduce the observed population, the simulations must be able to generate planets even closer to the central mass. However, setting the radius closer to the star will require an even shorter time step, and this will in turn increase the necessary computation time. Our simulations were carried out for 20 Myr due to the project time limit. It would be interesting to explore whether and how the planetary systems would change over even longer timescales and subjected to stellar and planetary tides, since in Sect. 3.1 we see that some of the systems are still evolving after 20 Myr. We ran a third of our simulations (the most AMD-unstable ones) for another 20 Myr, and for the majority of the runs the extended time resulted in no significant changes in the architecture of the systems, their long-term stability, etc. However, some of the simulations show substantial changes in their architecture; therefore, in future studies, we should simulate at least 40–50 Myr of evolution. So this limitation can also be potentially addressed and maybe eliminated in the future work, if we have sufficient computational resources.

Fig. 29 Planetary systems formed after 10 Myr of simulation with the parameters central mass = 0.8 M⊙, initial disk mass = 20 M⊕, and gas surface density = 1250 g cm−2, with a time step of 0.7305 days (simulation 61) and 0.182625 days (simulation 61_st).

5 Conclusions

We have explored the effects that various initial model conditions and configurations have on the final outcomes of 20 Myr of planetary accretion around two K-dwarf stars with different masses. We find that:

• The more massive the solid disk we start the simulation with, the more mass will be projected into the total mass of the planets in each simulated system. Generally, our simulations have a high planet formation efficiency as 60% to 84% of the initial mass will end up in the planets. Additionally, the scaling between the disk mass and the planet mass (of the largest and second-largest planet) in our systems shows a mostly logarithmic increase with empirical fits of $M_{\text{pl}}\propto M_{\text{disk}}^{0.6}$ for a 0.6 M_⊙ star and $M_{\text{pl}}\propto M_{\text{disk}}^{0.1-1.0}$ for a 0.8 M_⊙ star. The comparison with non-migrating simulations shows that migration results in larger variations in planet masses (of the largest planets) and generally a slower increase in planet mass with the initial disk mass.

• We manage to reproduce the main characteristics and architectures of the known systems, and produce mostly long-term stable systems after 20 Myr of evolution with an initial disk mass of 10 M_⊕ and above and a gas surface density value of 1500 g cm⁻² and above. Our simulations also produce many low-mass planets with M_p < 2 M_⊕, both close to the star and farther away, which are not found in the observed sample. This appears to be due to the observational biases as shown by the performed synthetic observations of our systems.

• The average planet mass in the observed K-dwarf sample (when not considering the giant planets) is 〈M_p〉 = 5.9 ± 2.5 M_⊕. Our data, specifically the average planet mass 〈M_p〉 and the mass of the most massive planet M_max, show that simulations that start with a disk mass of < 15 M_⊕ do not reach the masses necessary to reproduce the average mass of the observed planets. It seems that the model needs a higher initial disk mass.

• Our cumulative mass distributions for planets with masses between 2 and 12 M_⊕ show a strong preference for planets with masses M_p < 6 M_⊕ and a lesser preferences for planets with larger masses. The cumulative mass distributions for the observed samples increase almost linearly. This suggests that we are missing a mechanism for creating more massive planets.

• Around 75% of our systems do not contain multiple-planet resonant chains, and approximately 85% of the systems underwent at least one late instability. These numbers are consistent with the lower estimate determined based on the observations and would probably increase if the simulation times were extended.

• Earth-mass planets form quickly around 0.6 and 0.8 M_⊙ stars, mostly before the gas disk dissipates. The final systems after 20 Myr of evolution contain only a small number of planets with masses M_p > 10 M_⊕, and these formed after the gas was mostly gone.

Future models that span a larger range of stellar and planetary masses and characteristics as well as properties ofprotoplanetary disks, and which use longer simulation times combined with new observations, will help us further improve our understanding of planet formation around K-dwarf stars.

Acknowledgments

The authors would like to thank the reviewer Sean Raymond for taking the time and effort necessary to review the manuscript. We sincerely appreciate all valuable comments and suggestions, which helped us to improve the quality of the manuscript. This study is supported by the Research Council of Norway through its Centres of Excellence funding scheme, project No. 223272 CEED (P.H., E.M., S.C.W.). It has been done at the Centre for Earth Evolution and Dynamics (CEED) at the University of Oslo. We also acknowledge financial support from the Research Council of Norway through its Centres of Excellence scheme, project number 332523 (PHAB). Our numerical simulations have been performed on Norwegian supercomputers Betzy and Saga operated by Sigma2 as a part of Stephanie C. Werner’s project NN9010K. The authors would like the thank Simon Grimm and Joachim Stadel for allowing them to utilize their software GENGA and for helping with the use of the software. The study has made use of the NASA Exoplanet Archive operated by the California Institute of Technology, and The Extrasolar Planets Encyclopaedia (http://exoplanet.eu). Additionally, the authors acknowledge the help from the members of the Earth and Beyond group at CEED, particularly Yutong Shan for her helpful comments and discussions. R.B. would also like to thank the Research Centre for Astronomy and Earth Sciences (Budapest, Hungary) for financial support.

Appendix A Cumulative distributions of the simulated semimajor axes and planetary masses: Statistics

The two-sample K-S test was used to test whether the two pairs of samples come from the same distribution. The p-values are presented in Tables A.1 and A.2; IDM = 10+ means that only the systems starting with the disk mass of 10 M_⊕ or more are considered. Based on the commonly used 5% significance level, the values above 0.5 confirm that the samples are sampled from the same distribution. In some cases the simulated samples are quite small (the smallest containing only 20–30 planets), but K–S test can be reliably used for the samples of these sizes as well.

References

All Tables

All Figures

Fig. 1 Currently known multi-planet systems around K-dwarf stars. Only planets with known masses are presented. The sample was retrieved from https://exoplanetarchive.ipac.caltech.edu/ on December 4, 2022. Giant planets with masses > 30 M⊕ are shown in orange, and all planets with lower masses in blue. Terrestrial Solar System bodies are included for comparison (in red). The point size indicates the masses in M⊕ of the observed planets. Note the different size scale of the different colored bodies (blue and orange). Additionally, the markers representing the Solar System bodies are enhanced by a factor of 10 with regard to the planets indicated in blue for better visualization.

In the text

	Fig. 2 Initial solid surface density distribution for simulations with the starting planetesimal locations between 0.2 and 2 AU and different initial disk masses (IDM). The increasing linewidth with the values on the y-axis represents the growing fraction of larger bodies in the disk with increasing initial disk mass.
In the text

Fig. 3 Initial distribution, masses, and radii of 12 000 bodies for a simulation with an initial solid disk mass of 10 M⊕. The x-axis represents the distance from the central star in AU, and the y-axis the radii (left panel) and masses (right panel). Red, blue, and black indicate different mass ranges of the objects, i.e., mass M < 10−3 M⊕ (considered as planetesimals), 10−3M⊕ < M < 10−2 M⊕ (proto-embryos), and M > 10−2 M⊕ (embryos), respectively. There is no difference in the treatment of these objects during the N-body simulations. The isolation mass, which is the mass of embryos that accreted all the bodies within their feeding zone, is represented by the dashed line.

In the text

Fig. 4 Dynamical evolution of planets in one of the typical systems that underwent several late dynamical instabilities after the gas dispersal. The plots show the temporal evolution (20 Myr) of planetary bodies, specifically their semimajor axes, masses, eccentricities, and orbital inclinations. Each line color represents an individual object; however, the color is random in each plot. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. For better visualization, the panels on the left (semimajor axes and masses) show the evolution of all bodies with masses > 0.1 M⊕, and the panels on the right (eccentricities and orbital inclinations) show only bodies with masses > 0.5 M⊕.

In the text

Fig. 5 Dynamical evolution of planets in a less common system that did not undergo any late dynamical instabilities after the gas dispersal. The plots show the temporal evolution (20 Myr) of planetary bodies, specifically their semimajor axes, masses, eccentricities, and orbital inclinations. Each line color represents an individual object; however, the color is random in each plot. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. For better visualization, the panels on the left (semimajor axes and masses) show the evolution of all bodies with masses > 0.1 M⊕, and the panels on the right (eccentricities and orbital inclinations) show only bodies with masses > 0.5 M⊕.

In the text

Fig. 6 Dynamical evolution of planets in an extended simulation (100 Myr) that underwent multiple dynamical instabilities much later after the gas dispersal, at approximately 21 and 28 Myr. The plots show the temporal evolution of planetary bodies, specifically their semimajor axes, masses, eccentricities, and orbital inclinations. Each line color represents an individual object; however, the color is random in each plot. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. For better visualization, the panels on the left (semimajor axes and masses) show the evolution of all bodies with masses > 0.1 M⊕, and the panels on the right (eccentricities and orbital inclinations) show only bodies with masses > 0.5 M⊕.

In the text

	Fig. 7 Histogram showing a number of collisions during the 20 Myr of evolution in all our simulations. Only collisions where both bodies are above the lunar mass (MMoon = 0.012 M⊕) are presented.
In the text

Fig. 8 Architecture of our simulated planetary systems around a 0.6 M⊙ and a 0.8 M⊙ K-dwarf star, for two or four different gas surface density values and different initial disk masses (IDM) after 20 Myr of N-body integration using GENGA. Outcomes for the initial disk mass 5 M⊕ are not displayed. The point size indicates the masses in M⊕ of the formed planets. The gray zone is the region inside the inner truncation radius, and the dashed gray line shows the inner edge of the gas disk.

In the text

	Fig. 9 Mtot,final/Mtot,initial ratio versus the initial solid disk mass, Mtot,initial. The total mass of a system, Mtot,final, is calculated by summing up the masses of all individual planets in each system. Values for 0.6 M⊙ and 0.8 M⊙ star masses and various gas surface densities are presented.
In the text

Fig. 10 Average masses of the largest 〈M1〉 (filled circles) and second-largest 〈M2〉 (open circles) planets against the corresponding initial disk masses (IDM) for the two central star masses 0.6 M⊙ (in green, the top plot) and 0.8 M⊙ (in blue, the bottom plot), and the empirical fits for 〈M1〉 (solid line) and 〈M2〉 (dotted line). Since the systems starting with 5 M⊕ typically contain only one planet, they are not included in the plots.

In the text

	Fig. 11 Mass of the most massive planet, Mmax (M⊕), in a system for each initial disk mass value (IDM) and both star masses, 0.6 and 0.8 M⊙, shown in green and blue, respectively. Mean values and their standard errors are included (green squares with the dotted line and blue circles with the dashed line for 0.6 and 0.8 M⊙, respectively).
In the text

	Fig. 12 Mass-distance relationship for the planetary systems. Magenta and blue circles respectively represent our simulated sample and the observed sample around K dwarfs. The gray zone is the region inside the inner truncation radius, and the dashed gray line shows the inner edge of the gas disk.
In the text

Fig. 13 Dynamical evolution of one of the single planet systems. The plots show the temporal evolution (30 Myr) of planetary bodies, specifically their semimajor axes and masses. Only bodies with masses > 0.1 M⊕ are displayed. Each line color represents an individual object; the same object is indicated by the same color in both plots. The dashed black line shows the estimated time of the gas dissipation at ~8 Myr. Only one object (in light purple) has a mass above 0. 5 M⊕ and is classified as a planet in this study.

In the text

Fig. 14 Eccentricities (on the left) and inclinations (on the right) of the simulated planets around a 0.6 M⊙ star (in green; the top row) and a 0.8 M⊙ star (in blue; the bottom row) versus their semimajor axes. Different shades of green and blue denote the different initial values of the gas surface densities at 1 AU, i.e., 1000, 1250, 1500, and 1750 g cm−2; the higher initial density, the darker color of the circles. The point size indicates the masses in M⊕ of the formed planets. Eccentricities and inclinations of terrestrial Solar System planets are displayed for comparison (in red); their sizes are enlarged by a factor of 2 for better visualization.

In the text

Fig. 15 Temporal evolution (20 Myr) of the total planet mass (sum of all planetary bodies with MP > 0.5 M⊕ in a system) in the simulated system called 73 Sim (initial disk mass = 25 M⊕, central mass = 0.8 M⊙, gas surface density = 1250 g cm−2). Shown are the evolutionary tracks of four simulations with exactly the same parameters. The total mass included in Fig. 9 is the final state of 73_1 Sim, represented by the blue line.

In the text

Fig. 16 Snapshots of the 73_1 Sim simulation showing eccentricity versus the semimajor axis of the planetary bodies. Panel A shows simulation at ~8.6 Myr of evolution, and each subsequent panel displays the simulation a few thousand years later (B–F). Snapshot A presents a resonant chain of planets located very close to the star, which gets disrupted by another planet that gets too close and excites the group (B). This dynamical instability results in a ~1 M⊕ planet at first (C), and later on a ~3 M⊕ planet (F) falling onto the central star.

In the text

	Fig. 17 Outcomes of four 73 Sim simulations with initial disk mass = 25 M⊕, central mass = 0.8 M⊙, and gas surface density = 1250 g cm−2. The original simulation is 73_1 Sim. The size of the circles indicates the mass of the formed planets.
In the text

Fig. 18 Temporal evolution (20 Myr) of the mass of the most massive surviving body in the system. It shows the time it takes to grow an Earth-mass planet for each system (where the solid lines representing each simulation meet the dashed line). The simulations are grouped by their initial disk mass (IDM), indicated by the same color. Bodies in the systems starting with the initial disk mass 5 M⊕ never reach 1 M⊕ and are therefore not included in the plot.

In the text

Fig. 19 Time it takes to grow an Earth-mass planet depending on the initial disk mass (IDM) for the individual planetary systems (circles) and for the average values of each mass bin (purple squares with error bars). The dashed purple line indicates the descending trend of needed time with increasing disk mass. The light blue and dark blue circles represent the times from the simulations starting with the gas surface densities = 1000 g cm−2 and 1750 g cm−2, respectively.

In the text

	Fig. 20 Period ratios Pout/Pin of all pairs of subsequent planets for our simulated population of planets (in magenta) compared to the period ratios of the pairs from our MKG sample (see the hatched area of the histogram). Main orbital resonances between the planet pairs are indicated by the dashed gray lines and identified by the fractions at the top of the plot.
In the text

	Fig. 21 Logarithmic spacing between pairs of neighboring planets in each planetary system. The constant c2 is plotted against the values of constant c1. Magenta and blue circles represent our simulated sample and the observed sample around MKG-dwarf stars, respectively.
In the text

Fig. 22 Dynamical spacing, Δ, between all pairs of adjacent planets in each system in units of mutual Hill radius plotted against their period ratios, Pout/Pin. The figure shows our simulated sample (in magenta) together with the observed sample around MKG dwarfs (in blue). Two important values of dynamical spacing are represented by the dashed lines: the Hill stability criterion $\Delta =2\sqrt{3}$ (Chambers et al. 1996) in orange and Δ = 10 in gray.

In the text

Fig. 23 AMD stability of the simulated systems. Each planet is represented by a circle, the size of which is proportional to the mass of the planet. The color represents the AMD stability coefficient, β (in log scale), of the inner pair associated with the planet. The innermost planet is represented by the AMD stability coefficient associated with the star. Numbers in the simulation names (y-axis) are random, but the systems are ordered by their initial disk mass value, which decreases from top to bottom. Systems starting with the disk mass = 5 M⊕ contain only one or two very small planets and are therefore not included in the plot (but they are AMD-stable).

In the text

	Fig. 24 Mass ratios of the outer/inner neighboring planet of the entire sample (the left plot) and period ratios of the outer/inner planet pair within 0.5 AU of the star (the right plot) in the simulated population. The dashed gray line is the 1:1 line.
In the text

Fig. 25 Masses versus semimajor axes of planets of our simulated planetary systems (in magenta) compared to the currently known systems around K dwarfs (in blue). Only planets with known masses are presented. Known planets with masses above the super-Earths mass range (above 10 M⊕) are indicated in light gray. Here, we use the classification of planets by mass proposed by Stevens & Gaudi (2013). The point size indicates the masses in M⊕ of the planets. The sample of the known systems around K-dwarf stars was retrieved from https://exoplanetarchive.ipac.caltech.edu/ on December 4, 2022.

In the text

Fig. 26 Cumulative distributions of the simulated semimajor axes (the top row) and planetary masses (the bottom row) for the whole population (the left column), for different gas surface densities (the middle column), and for different initial disk masses (IDM, the right column). IDM = 10 + M⊕ means that only systems starting with a disk mass of 10 M⊕ or more are considered for that distribution. The same is true for the other IDMs as well as for the gas surface density values, e.g., 1000+ g cm−2. IDM = 5+ M⊕ contains only a few very small planets compared to IDM = 10 + M⊕; therefore, we omitted that distribution. Cumulative distributions for the K-dwarf (blue dotted line) and MKG (dashed black line) samples are displayed for comparison.

In the text

Fig. 27 Synthetic EPOS (Mulders et al. 2018, 2019) observations of our simulated population (A) and the associated statistics of the systems (B–D). Panel A shows the orbital periods and planet radii (calculated using a mass-radius relation from Chen & Kipping 2016) of detected planets (in magenta) and undetected planets (in gray). The statistics in panels B–D show the intrinsic distribution of our systems’ properties with dotted gray lines, the synthetic observable distribution in magenta, and the distribution observed with Kepler in dashed red lines. Plot B displays the distribution of period ratios between adjacent planet pairs. Plot C shows the frequency of multi-planet systems, and plot D the distribution of radius ratios of the adjacent planet pairs. All statistics are calculated for planets with orbital periods between 1 and 400 days and radii between 0.5 and 5 R⊕.

In the text

Fig. 28 Orbital period ratios (magenta) and innermost planet locations (orange). This plot was generated using EPOS (Mulders et al. 2018, 2019). The period ratios of adjacent planet pairs are plotted at the semimajor axis of the outer planet. The histogram on the right shows the marginalized period ratio distribution compared to that of Kepler. The semimajor axis of the innermost planet in each system is displayed in orange. The histogram on top shows the marginalized distribution of the innermost planets compared to the distribution derived from Kepler (red line). Both Kepler distributions are derived by Mulders et al. (2018).

In the text

	Fig. 29 Planetary systems formed after 10 Myr of simulation with the parameters central mass = 0.8 M⊙, initial disk mass = 20 M⊕, and gas surface density = 1250 g cm−2, with a time step of 0.7305 days (simulation 61) and 0.182625 days (simulation 61_st).
In the text