© The Authors 2024

1 Introduction

In contrast to the Solar System, a large fraction of exoplanetary systems harbour super-Earth/mini-Neptune planets in orbits interior to that of Mercury (e.g. Mayor et al. 2011; Fressin et al. 2013; Mulders et al. 2018). These systems of close-in planets normally consist of several planets that can be in or close to mean-motion resonances (MMRs) (e.g. Fabrycky et al. 2014; Pichierri et al. 2019). Furthermore, it seems that the majority of these systems behave like ‘peas in a pod’, where the planets have roughly equal sizes and there is an equal distribution of period ratios within each system (Millholland et al. 2017; Weiss et al. 2018; Millholland & Winn 2021). However, their formation pathway is still under debate.

While it is clear that the formation of systems in resonance requires planetary migration (e.g. Ogihara et al. 2015; Izidoro et al. 2017, 2021; Hühn et al. 2021), it is unclear whether the planets originate from beyond the water-ice line or formed in the inner regions of the disc. This is independent of whether or not the planetary systems are in resonance. For selected systems, the chemical composition of the planets could be used to decipher the location of their formation (e.g. Benneke et al. 2019; Sun et al. 2024), as the migration behaviour sets the composition of these planets (e.g. Bitsch et al. 2019b, 2021; Schoonenberg et al. 2019; Izidoro et al. 2022; Chatziastros et al. 2024). However, even with JWST operating and the future ARIEL mission, we are still far away from a bulk characterization of planets in these systems. In addition, a spectroscopic analysis of the planetary atmospheres requires that the planets be transiting, and that they orbit stars close enough to allow a sufficient signal-to-noise ratio. However, for most of the Kepler sample, this remains a challenge due to the large distance of the Kepler field from Earth, leaving it unclear as to whether the majority of the close-in planets formed interior or exterior to the water ice line from a detailed compositional point of view.

On the other hand, the analysis of the Kepler sample revealed a valley in the radius distribution of these close-in planets (Fulton et al. 2017; Van Eylen et al. 2018), where peaks at ≈1.5 and ≈2.0 Earth radii and a dip at ≈1.8 Earth radii were found. The interpretation of this radius valley is that it reflects a dichotomy in planetary properties: while the close-in planets with radii of less than 1.5 Earth radii are generally thought to be atmosphereless rocky planets, planets with radii larger than 2.0 Earth radii are thought to be mini-Neptunes with rocky cores and hydrogen–helium envelopes or ice-dominated planets (e.g. Owen & Wu 2013; Lopez & Fortney 2014; Jin & Mordasini 2018). This indicates a compositional difference between these types of planets and may suggest a different origin. The radius valley in itself seems to be universal in the sense that it appears at all stellar masses (Ho et al. 2024).

The radius valley can be explained by photoevaporation of close-in planetary atmospheres (e.g. Owen & Wu 2013; Lopez & Fortney 2014) or by core-powered mass loss (e.g. Gupta & Schlichting 2019). Alternatively, atmospheres of planets could be lost during collisions with other planets (e.g. Liu et al. 2015; Inamdar & Schlichting 2016). The current data support all of these scenarios and are in particular consistent with a combined scenario, where all these effects can act at the same time (Izidoro et al. 2022).

As the planets migrate inwards, they will stop at the inner edge of the disc (Masset et al. 2006; Flock et al. 2019), where they can then build up resonance chains (Izidoro et al. 2017, 2021; Ogihara et al. 2018). As the gas disc disappears, the damping forces acting on the planets diminish and the resonance chains can be broken by instabilities, which rearranges the planetary orbits and increases the eccentricities and inclinations of the close-in planets. These instabilities can either be self-triggered (Izidoro et al. 2017, 2021) or can be triggered due to the presence of outer giant planets (e.g. Bitsch & Izidoro 2023). In addition, the instabilities are affected by the proximity of the planets to each other (e.g. Gladman 1993) as well as by their masses. The rearrangement of the planets within these systems can naturally explain the observed number of transiting planets, as well as the general period distribution of these planets (Izidoro et al. 2021), even though the simulations have trouble explaining systems in the 2:1 resonance configuration.

The trapping of planets into MMRs depends on the relative migration speed of the planets (e.g. Huang & Ormel 2023). Faster migrating planets can result in more compact systems than slower migrating planets. The migration velocity in itself, however, is set by the properties of the disc, such as its viscosity. At low viscosity, even mini-Neptunes can start to open partial gaps that can reduce their migration speed (e.g. Paardekooper et al. 2011; Kley & Nelson 2012; Bitsch et al. 2013). Consequently, systems forming in low-viscosity environments could harbour wider resonances, with consequences for the long-term evolution of these planetary systems. On the other hand, a lower migration speed might imply that the planets are not trapped in resonance at all at the end of the gas–disc phase, because they migrate too slowly to reach the inner edge of the disc (e.g. Brasser et al. 2017).

Past simulations of the ‘breaking the chains’ model focused on systems without giant planets in high-viscosity environments (Izidoro et al. 2017, 2021, 2022; Esteves et al. 2020, 2022, 2023) or on systems with outer giants, but in low-viscosity environments (e.g. Bitsch et al. 2020; Bitsch & Izidoro 2023). The low-viscosity simulations with giant planets already seem to suggest that systems with a 2:1 resonance configuration could be more common in these environments. Here we want to bridge this gap by analysing simulations within the breaking the chains model at low disc viscosity without outer giant planets.

This paper is structured as follows. In Sect. 2, we give a short overview of our model. We show the results of our simulations in Sect. 3 and compare these with the results of the high-viscosity simulations of Izidoro et al. (2021). We then discuss our results in the context of other works in Sect. 4 and summarise our conclusions in Sect. 5.

2 Methods

N-body simulations are needed to study the dynamics and evolution of planetary systems. To this end, we employ FLINT-STONE, an N-body framework based on the MERCURY integrator (Chambers 1999) that includes planetary growth and migration (Bitsch et al. 2019a; Izidoro et al. 2021). We use the same code here as that described in detail in Izidoro et al. (2021) and Bitsch et al. (2019a), and only summarise the main ingredients of our simulations.

We initialise our simulations with 30 planetary embryos separated by 0.25 AU starting at 2.75 AU in a low-viscosity environment (α = 10⁻⁴ for all our simulations). In contrast, the simulations by Izidoro et al. (2021) used a viscosity of α = 5.4 × 10⁻³, which is about a factor of 50 more viscous. We note that the change of viscosity only influences the migration velocity of planets (see below and Appendix A), while the other disc parameters (e.g. gas surface density and temperature) are kept the same for simplicity. The gas surface density and temperature thus evolve in time exactly in the same way for all our simulations. This approach, while artificial, allows a clean approach to understand how the viscosity influences the formation of planetary chains. This approach is the same as in our previous works (Bitsch et al. 2020; Bitsch & Izidoro 2023). The planetary embryos of around Moon mass start in an already evolved protoplanetary disc that is 2 Myr old in the framework of Bitsch et al. (2015a), as in our previous work (Bitsch et al. 2019a, 2020; Bitsch & Izidoro 2023). The embryos have an initial eccentricity and inclination that is randomly assigned (e < 0.01 and i < 0.5 degrees).

At the start of the simulations, the planetary embryos start to grow via pebble accretion (see Johansen & Lambrechts 2017 for a review), where we follow the recipe of Johansen et al. (2015). We use an exponentially decaying pebble flux that corresponds to a total of 350 Earth masses over the 3 Myr lifetime of the gas disc. The pebble sizes are limited by drift (Lambrechts & Johansen 2014) and we reduce the pebble flux by 50% interior of the water ice line to account for water-ice evaporation. Furthermore, we then also reduce the pebbles sizes to millimetre (mm) size, in agreement with chondrules in the Solar System, as done in our previous work (Bitsch et al. 2019a; Izidoro et al. 2021). Once the planets have reached the pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018b; Ataiee et al. 2018), they stop growing by pebble accretion.

Once planets have reached the pebble isolation mass, where we ignore the dependency on viscosity for simplicity, they should start to accrete gas (Lambrechts et al. 2014; Bitsch et al. 2015b). The exact efficiency of gas accretion is influenced by the planetary core mass as well as by the envelope opacity (Ikoma et al. 2000; Schulik et al. 2019; Bitsch & Savvidou 2021; Brouwers et al. 2021; Moldenhauer et al. 2021). In addition, efficient atmospheric accretion can be hindered by recycling flows (Cimerman et al. 2017; Lambrechts & Lega 2017; Moldenhauer et al. 2021; Wang et al. 2023), which transport material that entered the Hill sphere away from the planet before it can be accreted, which is especially efficient in the inner disc regions. As we are primarily interested in the formation of super-Earth/mini-Neptune systems, we chose to not model the gas accretion process onto the planets as in Izidoro et al. (2021). In principle, one can interpret this approach as a limit of very high envelope opacity that prevents efficient contraction of the planetary atmospheres, resulting in the formation of planetary cores with a small hydrogen/helium envelope.

As the planets start to grow, they eventually become big enough to start to migrate efficiently via type-I migration, where we follow the formulae from Paardekooper et al. (2011). When the planets start to become even larger, they can open gaps in the protoplanetary disc, reducing their migration velocity, where we follow the approach of Kanagawa et al. (2018), as discussed in Appendix B. In addition, we apply forces to mimic the eccentricity and inclination damping onto the growing planets, following our previous approaches (Bitsch et al. 2019a; Izidoro et al. 2021). In principle, the damping of eccentricity and inclination is affected by partial gap opening (Pichierri et al. 2023, 2024), which results in a deviation from the damping formulae employed in this work, which are from Cresswell & Nelson (2008). We nevertheless use the same damping formulae as used in previous works in the interest of consistency, and will investigate the influence of the damping formulae in future works. Planets in low-viscosity environments can open partial gaps, resulting in slower migration compared to planets in high-viscosity environments (see Appendix A).

After 3 Myr, our protoplanetary disc dissipates on an exponential timescale, where we reduce the gas surface density uniformly in the disc with an exponential function on a timescale of 12 kyr within the last 100 kyr. Once the protoplanetary disc has dissipated, we integrate the simulations up to 100 Myr to allow sufficient time for instabilities to occur. During the instabilities, the planets can also collide. We treat collisions as perfect mergers between the planets, due to the fact that even more detailed collisional treatments that include fragmentation (Leinhardt & Stewart 2012) give the same results in our simulations (Esteves et al. 2022).

At the end of the simulations, we investigate the structure of the planetary systems and also run synthetic observations of the formed systems. The exact method for this is described in Izidoro et al. (2021). These synthetic observations are necessary to effectively account for observational bias and compare our results with observations.

The initial conditions (eccentricity, inclination, and planetary mass) are all randomised. In the interest of statistical robustness, we ran 50 simulations with slightly varying initial conditions. We show the resulting planetary systems after 3 (end of the gas–disc phase) and 100 Myr of evolution in Appendix C.

3 Planetary systems from low-viscosity environments

In this section, we discuss the outcome of our results. We first investigate the general outcomes of our simulations and compare these to systems formed in the simulations by Izidoro et al. (2021). We then investigate how synthetic observations of these formed systems can explain the period distribution.

3.1 Individual system evolution

Figure 1 shows the time evolution of a system with respect to its evolution in semi-major axis, planetary masses, eccentricity, and inclination. The planets initially grow by accretion of pebbles, resulting in a smooth increase in the planetary masses. As the planets grow above 1 Earth mass, they start to migrate into the inner disc regions and get caught in resonance configurations. As the pebble-isolation mass increases with orbital distance (see Appendix E), the planets forming further out in the disc become slightly more massive compared to the inner planets. During their migration, smaller planets can already be scattered.

Once the gas disc dissipates, the planetary eccentricities start to increase to a few percent, while the inclinations increase to around 1 degree for the massive inner planets. The eccentricities of the larger planets remain quite low, resulting in stable planetary systems for 100 Myr of evolution. Consequently a chain of resonant planets (5:4, 2:1 and 3:2 MMRs) survives in the inner disc.

Figure 1 shows the evolution of a stable system, which survived until 100 Myr of evolution. Our simulations show two other typical outcomes: Fig. 2 shows a system that underwent instabilities after the gas–disc phase that sculpted the inner system, resulting in only two remaining inner planets, while Fig. 3 shows the evolution of a system where the inner system is completely removed by scattering caused by the outer planets.

The systems shown in Figs. 2 and 3 experience strong scattering events after the removal of the gas discs, which leads to large increases in the planetary eccentricities to values of around 10% (i.e. around 10 times larger than the stable system in Fig. 1). This is a consequence of the dynamical instabilities in the outer parts of the system (r > 1 AU), which eventually affect the inner system (< 0.5 AU). Consequently, only two planets survive outside of any resonance configuration interior to 1 AU (Fig. 2). On the other hand, the instabilities for the system shown in Fig. 3 are so strong that the whole inner system is removed. However, the cause of the instabilities is the same: an increase in the orbital eccentricities in the outer system after the gas disc dispersal, leading to an instability that affects and removes the inner system.

3.2 Planetary systems directly after the gas disc phase

During the first 3 Myr of our simulations, the planets can grow via pebble accretion (and collisions) and migrate, while simultaneously their eccentricity and inclination are damped by gas–disc interactions. Once the gas disc dissipates, only gravitational interactions remain, leading to an increase in eccentricity and inclination with the consequence that instabilities will occur. However, the structure of the systems after the gas–disc phase can already provide some initial hints as to the differences between low- and high-viscosity planet formation environments.

Fig. 4 shows the number of planets that have formed at the end of the gas–disc phase (left), as well as their period ratios (right). We also show the systems that could be identified by observations as well as the synthetic observations of the simulations of Izidoro et al. (2021). In addition, we show the observed Kepler sample (thin grey lines) as well as the planetary chains observed by Kepler (thick grey lines). We note that we only consider planets interior to 0.7 AU as observable, even though more planets can be hidden exterior to that (see Appendix C). In addition, we only include planets above 1 Earth mass in our statistical analysis, because smaller planets are very likely to be missed in transit observations (e.g. Mulders et al. 2018)

Within 0.7 AU, our low-viscosity simulations show a maximum of six planets, while the high-viscosity simulations of Izidoro et al. (2021) show systems with up to ten planets. This difference can be explained by the slightly different setups used here. In particular, Izidoro et al. (2021) used 60 planetary embryos per simulation, allowing more planets to be present in the first place. Another difference is the migration velocity of the growing planets, which is determined by the disc viscosity. Discs with lower viscosity allow early (partial) gap opening, which results in a transition of the planets into a slower migration regime compared to their high-viscosity counterparts (see Appendix A). Consequently, not as many planets from the outer disc migrate interior to 0.7 AU, resulting in different systems.

While the systems in our simulations always have at least two planets with 1 Earth mass or more, observationally, the picture is different. Due to the mutual inclinations between the planets, nearly 30% of all systems would only show one transiting planet, while the number of systems with five transiting planets is significantly reduced compared to the pure simulation data. This is in line with the results of Izidoro et al. (2021), who observed the same behaviour in their simulations. Comparing to the observational constraints from Kepler clearly shows that the number of single transiting planets is not reproduced at the stage directly after dissipation of the gas disc; this requires some instability to increase the level of mutual inclination and to reduce the number of planets (see below).

The period ratios on the other hand show that most of the formed planets are in resonant configurations at the end of the gas–disc lifetime, which is caused by inward migration and trapping at the inner edge of the disc. This is, of course, in contrast to the observed period ratios by Kepler, where only a small fraction of systems are expected to be in a resonant configuration. This result is in line with the findings of Izidoro et al. (2021), where instabilities are necessary to break the resonant configuration; we study this aspect further below. However, there is a significant difference in the observed distribution of the period ratios between the simulations presented here and those of Izidoro et al. (2021): the resonant chains in Izidoro et al. (2021) are much more compact than those in the simulations presented here.

The difference in the period ratios of the systems can again be explained by planetary migration. As the planets migrate faster in the Izidoro et al. (2021) simulations, they can naturally come closer to each other than they are able to in simulations in a low-viscosity environment. In particular, the simulations by Izidoro et al. (2021) seem to overestimate planetary systems in the 7:6 and 5:4 MMRs, while simultaneously underestimating the planets in the 3:2 and 2:1 MMRs. In contrast, the low-viscosity simulations show a much better match at these MMRs. We show in Appendix D that this difference is truly related to disc viscosity rather than the number of embryos. On the other hand, here we only discuss the simulations at the end of the gas–disc lifetime, and instabilities are expected to change these configurations. We discuss this aspect further in the following section.

Fig. 1 Evolution of a single planetary system as a function of time regarding semi-major axis (top left), mass (top right), eccentricity (bottom left), and inclination (bottom right). The coloured lines correspond to the surviving inner (r < 4 AU) planets, while the grey lines depict the evolution of the outer (also less massive) planets. The vertical line at 3 Myr marks the dissipation of the gas disc phase, where the planets lose the damping forces on eccentricity and inclination from the disc. The inner system remains stable until 100 Myr. The system corresponds to system number 3 in Fig. C.1.

Fig. 2 Same as Fig. 1, but for a system that became unstable after the gas–disc phase with two surviving inner planets. We mark all surviving planets in colour. The planet marked in yellow stops growing at around 700 kyr, when it acquired an eccentricity too large to allow the accretion of pebbles. Once its eccentricity is damped below 1%, it continues to grow via pebble accretion. The system corresponds to system number 21 in Fig. C.1.

3.3 Planetary systems after the instability phase

After the gas phase, the planets can only interact gravitationally and are no longer influenced by the damping effects of the disc that previously acted on their eccentricity and inclination. Consequently, the planetary systems become unstable and their orbital configurations change. In the following, we classify systems that have not undergone a major rearrangement after the gas–disc phase as ‘chains’.

Fig. 5 shows the number of detected planets (left) and their period ratios (right) in our simulations after 100 Myr of system evolution. As before, we also show the distribution of the detected systems as well as the results from Izidoro et al. (2021) in comparison to the Kepler data.

As a consequence of the instability, many planetary systems have no inner system left and therefore show no inner planets (see Appendix C). Consequently, the number of systems with three or more inner planets is also greatly reduced. However, the reduced number and increased mutual inclination of the systems allow a much better match to the Kepler observations regarding the systems with only one transiting planet, in line with the results of Izidoro et al. (2021). The systems with multiple transiting planets also match much more closely after the instability phase.

The resulting period ratios are determined by the systems with at least two surviving planets (limiting our statistics to some level; see Appendix C). The general results of our low-viscosity simulations do not match the period distribution from the Kepler systems, as the systems in our simulations remain too tightly packed compared to the Kepler systems, even when applying synthetic observations. In contrast, the simulations of Izidoro et al. (2021) are a much better match to the Kepler systems in terms of their period ratio distribution.

The differences in these results could be explained by the different strengths of the instabilities between the different system configurations. While the instabilities in our simulations mostly result in a removal of planets from the system, the instabilities in the simulations of Izidoro et al. (2021) mostly show a rearrangement of the system, with more survivors, which leads to period ratio distributions that are in better agreement with observations. This is caused by the difference in the outer system architecture. While the systems formed in the high-viscosity environments show almost no massive planets exterior to ~200 days due to the fast inward migration of planets of a few Earth masses, the opposite is true at low viscosity, where even the more massive planets cannot migrate to the inner disc (see Appendix A), resulting in outer systems with several planets above a few Earth masses. In combination with the lower planetary masses of the inner systems in the low-viscosity case (see Appendix E), the instabilities triggered by the outer systems are more likely to result in the removal of the inner planets rather than simply a dynamical rearrangement. Using a mixture of 2% stable and 98% unstable systems (as done in Izidoro et al. 2021 to achieve the best match to observations) does not allow a good match to the overall period ratios (see Appendix F) of the low-viscosity simulations as well.

In Fig. 6, we plot the number of observed planets and their period distributions only for the chains formed in our simulations, as well as the chains from Izidoro et al. (2021) after 100 Myr of evolution. The chains from our simulations harbour only four or more planets, in contrast to the Kepler observations; this is also different from the simulations of Izidoro et al. (2021), who also show chains with a smaller number of planets. In addition, the high-viscosity simulations predict chains with seven or more planets, which have not yet been detected in such abundance, but should exist in theory (e.g. Trappist-1).

On the other hand, the surviving chains show a wider period ratio distribution in the low-viscosity simulations than in the high-viscosity simulations. In particular, the low-viscosity simulations show more planets in the 3:2 and 2:1 configurations compared to the high-viscosity simulations, which leads to better agreement with the observational constraints from the chains alone.

The eccentricities of the planets are shaped by the gravitational interactions between the bodies, where more violent instabilities can cause larger eccentricities. Consequently, the eccentricity distributions of the two sets of simulations (Fig. 7) are different. The low-viscosity simulations experience more violent instabilities, which are caused by the outer planets, while the higher viscosity simulations do not feature outer planets that are both sufficiently massive and sufficiently close to the inner system to affect it. Consequently, they experience only ‘mild’ instabilities, resulting in planets with larger eccentricities in the low-viscosity simulations.

In contrast, the eccentricity distribution of the surviving planets in chains (right in Fig. 7) is very similar. This is caused by the fact that chains, if they are to survive, must also some-how avoid being strongly affected by external planets and must therefore keep their eccentricities low. In the planetary systems shown in Figs. 3 and 2, the eccentricities increase above 10%, resulting in violent instabilities that break the chains. In contrast, systems that remain in a resonant configuration have low eccentricities (Fig. 1), as the low eccentricities prevent orbit crossing and instabilities in general.

We note that, in principle, more compact resonance configurations (e.g. 7:6, 6:5, 5:4) should result in more violent instabilities compared to wider resonance configurations (e.g. 3:2, 2:1). An example of a well-described less violent instability would be the dynamical instability of the Solar System in the Nice model, although it features more massive planets on much larger orbits (e.g. Tsiganis et al. 2005; Nesvorný 2018). In the cases analysed here, instabilities tend to be more violent, not necessarily due to the compactness of the resonant configurations themselves, but instead due to the effect of external perturbers (see Figs. 3 and 2). We will explore this mechanism in more detail in future work.

Fig. 3 Same as Fig. 1, but for a system that became unstable after the gas–disc phase with no surviving inner planets. We mark all surviving planets in colour. The planet marked in orange acquires some eccentricity during the gas–disc phase, preventing pebble accretion, and it only continues to grow once its eccentricity is damped below 1%. The system corresponds to system number 1 in Fig. C.1.

Fig. 4 Number of synthetically observed planets (left) and period ratios of adjacent planets (right) for the planetary systems directly at the very end of the gas–disc phase, before the systems undergo instabilities. The thin grey line marks all the observations, while the thick grey line marks the observed chains of planets. The red colour corresponds to the direct results of our simulations, while the black line corresponds to the synthetical observations of our simulations. The purple line represents the synthetic observations from Izidoro et al. (2021). A clear difference is noted in the compactness of the chains between the high-viscosity simulations of Izidoro et al. (2021) and the low-viscosity simulations used here.

Fig. 5 Number of synthetically observed planets (left) and period ratios of adjacent planets (right) for the planetary systems after 100 Myr of integration; i.e. after the instability phase. The colours are the same as in Fig. 4. The high-viscosity simulations from Izidoro et al. (2021) produce wider systems after the instabilities compared to the low-viscosity simulations, which is also a consequence of the real number of planets in the systems before instabilities (left in Fig. 4), where more planets are initially present at high viscosities.

Fig. 6 Number of synthetically observed planets (left) and period ratios of adjacent planets (right) for chains of planets that are stable after 100 Myr. The colours are the same as in Fig. 4.

Fig. 7 Eccentricity distribution of our low-viscosity simulations and of the high-viscosity simulations of Izidoro et al. (2021) for all planets (left) and only for the surviving chains (right).

4 Discussion

We show that the outcome of our simulations depends crucially on the viscosity of the disc (see Appendix D). We now discuss the influence of the different (changed) parameters on the outcomes of our simulations.

4.1 Planetary migration

The planetary migration rates are not only set by the viscosity of the disc (e.g. Paardekooper et al. 2011), but also by its other properties, and in particular its surface density and temperature (which in turn set the entropy). Here we use the migration formulae by Paardekooper et al. (2011), which allow reduced (or even outward) migration rates at higher viscosities, depending on the disc profile. We also use the disc profile of Bitsch et al. (2015a), where the disc has already evolved for 2 Myr. At this stage, the aspect ratio of the inner disc is relatively small (H/r ≈ 0.025–0.03), resulting in an easier partial gap opening by small planets (Crida & Morbidelli 2007; Bitsch et al. 2013; Kanagawa et al. 2018; Pichierri et al. 2023), which reduces their inward migration speed compared to simulations at higher viscosity (see Appendix A), where gap opening is hindered. Consequently, the exact outcomes of our simulations are influenced by our choice of disc profile.

We also note that other migration formulae (e.g. Jiménez & Masset 2017) give different migration rates, where the heating torque in particular can change the migration direction of fast-growing small-mass planets (e.g. Benítez-Llambay et al. 2015; Masset 2017). However, the heating torque diminishes once the planets reach the pebble-isolation mass, resulting in nominal inward migration at low viscosity. The overall outcome of planet-formation simulations is not greatly influenced by the heating torque or different migration prescriptions (Baumann & Bitsch 2020); however, the detailed migration rates and disc parameters are important in the formation of specific systems (e.g. Hühn et al. 2021).

The relative migration velocities of the planets set the resonances in which planets can be trapped (e.g. Huang & Ormel 2023). This is independent of the factor that sets the migration rates, be it disc profile or viscosity. Here we explore how low viscosities influence the migration rates. A similar effect could also be achieved via different disc parameters, such as a lower overall disc surface density. However, this would also influence the rate at which planets grow via pebble accretion, as well as the damping rates of eccentricity and inclination, making it more difficult to make direct comparisons and draw robust conclusions.

In order to test the influence of disc viscosity on resonance trapping, we show in Appendix G the outcomes of simulations of two migrating equal-mass planets in different disc–viscosity environments. While the exact outcomes of the simulations depend on the disc parameters, it is clear that the viscosity – which sets the migration rate – determines the resonance that the planets are trapped in. Planets migrating in high-viscosity environments are trapped in closer resonances compared to planets migrating in low-viscosity environments.

4.2 Planetary growth

The growth of the planets is initially driven by pebble accretion (e.g. Johansen & Lacerda 2010; Ormel & Klahr 2010; Lambrechts & Johansen 2012). While pebble accretion is an efficient mechanism of planet growth, individual pebble-accreting planets are actually very inefficient in reducing the overall inward flux of pebbles (Lambrechts et al. 2014; Bitsch et al 2019a), as long as they do not reach the pebble isolation mass As the pebble flux is not reduced significantly by one growing planet, systems of multiple planets can efficiently form (Lambrechts et al. 2019; Bitsch et al. 2019a; Izidoro et al. 2021; Matsumura et al. 2021). In addition, the growth in the outer regions of the disc happens on longer timescales, meaning that the planets forming in the inner system reach the pebble isolation mass more quickly than those forming in the outer regions, diminishing the influence of more outer embryos on the growth of inner embryos (Bitsch et al. 2019a, 2020; Bitsch & Izidoro 2023).

In our simulations, we used 30 planetary embryos, while the simulations of Izidoro et al. (2021) used initially 60 planetary embryos. Consequently, more planets can initially grow in the simulations of Izidoro et al. (2021), but they are also distributed in a wider radial space compared to our simulations Furthermore, planets forming in the outer regions (e.g. exterior to 10 AU) will not migrate interior to 0.7 AU in low-viscosity environments (see Appendix A), and remain in the outer disc, causing late instabilities of the inner systems. Consequently, we think that a greater number of embryos in the low-viscosity simulations would not significantly change the outcome of our simulations with respect to the inner system structure (see also Appendix D).

On the other hand, in the high-viscosity simulations of Izidoro et al. (2021), planets migrate inwards faster and can thus reach the inner edge of the disc even if they originate from beyond 10 AU. Consequently, more planetary embryos can migrate into the inner region of the disc. This can result in a squeeze of the formed planetary chains, where incoming planets can push already trapped planets into tighter resonances. In combination with the faster migration speed, which can allow tighter resonances as well, the formed systems are more tightly packed, with more planets in the high-viscosity simulations.

In Appendix D, we show a comparison to the simulations of Bitsch & Izidoro (2023), who also performed simulations with 15 and 60 initial embryos. These simulations support the idea that the migration velocity – set by disc viscosity – is responsible for the different chain configurations rather than the number of planetary embryos. In order to verify this result, we conducted simulations with only two planets in discs with different viscosities in Appendix G, confirming this phenomenon in simulations very clearly.

4.3 Formation of Kepler-223 and TOI-178

The Kepler-223 system consists of four super-Earths locked in a resonance chain with 4:3, 3:2, and 4:3 resonances (Mills et al. 2016). Understanding the formation of a particular resonance chain requires detailed studies of the migration history of the planets. Hühn et al. (2021) carried out an extensive parameter study using N-body simulations with prescribed migration rates to explain the migration history of Kepler-223. In particular, these authors varied the gas surface density value and the disc viscosity in order to find a parameter combination that allows the formation of that system.

Hühn et al. (2021) find that surface densities higher than the Minimum Mass Solar Nebular, or MMSN (Weidenschilling 1977; Hayashi 1981) (which are also higher than the gas surface densities used here) and α-viscosities with α ≤ 10⁻³ favour the formation of the orbital configuration of Kepler-223. This result is in line with our low-viscosity simulations, which can reproduce a significant fraction of systems with 4:3 and 3:2 resonances. This result is encouraging considering the completely different setup of the simulations of Hühn et al. (2021) and ours, indicating that lower viscosities could favour the production of systems in the 4:3 and 3:2 resonances.

The TOI-178 system consists of six transiting planets following a 2:4:6:9:12 chain of Laplace resonances (Leleu et al. 2021). Our simulations indicate that such wide chains could be better explained by low-viscosity environments, as the chains built in high-viscosity environments do not produce planets in the 2:1 resonance configuration.

4.4 Consequences for inner and outer system structure

Our simulations show that the structure of the inner system is influenced by the outer regions (P>200 days). While the systems that form in the low-viscosity environments experience dramatic instabilities that can result in a large fraction of completely destroyed systems, the instabilities in the high-viscosity systems mostly simply break the chains but leave the overall systems intact. We note that these instabilities do not necessarily eject the planets from the system; indeed, a large fraction of the planets are thrown into the host star, where they can influence the stellar atmospheric abundances (e.g. Liu et al. 2024). This difference is caused by the greater number and larger masses of planets in the outer disc regions (exterior to 1 AU; see Fig. C.1), which trigger the violent instabilities in the systems formed in the low-viscosity environments. This is in line with the simulations of Bitsch & Izidoro (2023), where outer giant planets destroy the inner systems, with the difference that the systems in Bitsch & Izidoro (2023) already become unstable during the gas–disc phase due to the greater impact of the giants compared to that of the outer sub-Neptunes in the simulations presented here.

This difference has two important implications for the occurrence rates of planets in the inner systems as well as in the outer systems. Even though all our simulations allow the formation of inner super-Earths during the gas–disc stage, only around 50% of our final systems harbour inner planets after 100 Myr of integration. On the other hand, all our systems still harbour planets exterior to 1 AU after 100 Myr of integration. This result is broadly in line with the occurrence rate of sub-Neptunes exterior to the water-ice line from microlensing surveys (e.g. Suzuki et al. 2016), which predict that these planets should be more common than hot inner super-Earths.

We note that the presence of the outer systems (r > 1 AU) is influenced by our choice of the initial embryo distribution. If we were to include only a small number of embryos in our simulations (e.g. less than 8) within a few astronomical units (AU), all embryos would migrate into the inner disc regions and no outer system would exist. In this case, the instabilities would be less pronounced and the inner systems would survive, as they do in the high-viscosity simulations. However, the lower viscosity would still determine the formed resonance chains, because the lower viscosity allows slower migration speeds, trapping the planets in wider configurations. This can be seen in Fig. 4, where the period ratios of planets clearly show wider systems for the low-viscosity simulations compared to the high-viscosity simulations. On the other hand, if no outer embryos were to exist initially, it is unclear how the sub-Neptunes exterior to the water-ice line, as found in microlensing surveys (e.g. Suzuki et al. 2016), would come to exist.

4.5 Comparison to other planet-formation studies

The formation of super-Earth and mini-Neptune systems has been studied by various groups in general (e.g. Ogihara et al. 2018; Emsenhuber et al. 2021; Ogihara et al. 2024; Batygin & Morbidelli 2023) and specific systems have also been examined (e.g. Coleman et al. 2019). These studies differ in their methods regarding accretion (pebbles or planetesimals), disc evolution (viscous or disc winds), and initial embryo formation (in ring-like structures or distributed all over the disc), but they all include planetary migration, which is necessary for the build up of resonance chains. More notably, most studies use the migration prescription by Paardekooper et al. (2011), as we do in the work presented here. Furthermore, all simulations include an inner edge, where inward migration can stop and planets can pile up (e.g. Masset et al. 2006; Flock et al. 2019).

However, independently of how the planets form (pebbles or planetesimals) or the configuration in which they form (in rings or widely spread), they have to migrate to the inner edge of the protoplanetary disc. We believe migration speed, which is determined by the viscosity of the gas disc (as well as by the gas surface density and temperature; see Appendix B), to be the crucial component that determines the eventual resonance configuration of the planets. This belief is motivated by the fact that capture into resonances is determined by the relative migration velocity and not by the initial distance between the embryos.

On the other hand, as larger viscosities allow the inward migration of planets from beyond the water-ice line, a mixture in the composition of the systems is expected (e.g. Izidoro et al. 2022). At low viscosity, this may not be the case, unless the discs are very cold and the water-ice line is initially close to 1 AU (e.g. Bitsch et al. 2019b), allowing planets to reach the inner disc even in low-viscosity environments (see also Appendix A). However, how much of the water ice in pebbles can be retained by growing planets is still debated, as the pebbles might evaporate in the hot atmosphere of the planets, preventing the efficient accretion of water ice (Johansen et al. 2021; Müller et al. 2024).

5 Summary

We conducted simulations similar to those of Izidoro et al. (2021), but using a low-viscosity environment. This low-viscosity environment allows planets of a few Earth masses to already open partial gaps, reducing their rate of inward migration. Consequently, the planetary chains that form at the inner edge of the disc are in wider configurations compared to planets that form in a high-viscosity environment. This has consequences for the stability of the systems and for their observability. Our findings are as follows:

• The lower viscosity causes slower migration speeds and consequently planets remain in wider configurations after the gas–disc phase, where planets are mostly trapped in the 4:3, 3:2, and 2:1 resonances. The higher-viscosity simulations, on the other hand, do not populate the 2:1 resonance directly after the gas–disc phase and only marginally populate the 3:2 resonance. Most planets are in the 4:3 resonance or exhibit an even tighter configuration;

• Following dissipation of the gas disc, the planetary systems experience instabilities, which break the resonance chains. Consequently, the number of observed planets can be matched, as also shown in our previous works (Izidoro et al. 2021; Bitsch & Izidoro 2023);

• The overall period ratios of the low-viscosity simulations, after instabilities, remain too narrow to explain the observed period ratios. In contrast, the high-viscosity simulations, in general, provide a better match to the overall period ratio distribution (Izidoro et al. 2021);

• On the other hand, the simulated multi-planet chains in the low-viscosity environment very closely match observations. In particular, our low-viscosity simulations can reproduce the large number of planets in the 3:2 and 2:1 configurations. This result is also not affected by outer giant planets, which can destabilise the inner systems (Bitsch & Izidoro 2023), indicating that low-viscosity migration might be the major driver of wider resonances than the 5:4 (or tighter) resonances, which should therefore be more likely to form in discs with high viscosities.

Clearly, a large variety of parameters influence the outcome of planet-formation simulations. Here we have identified that different migration rates (caused by different viscosities) have consequences on the chains that are formed in these simulations. As the low-viscosity simulations better match the chains than the high-viscosity simulations, and as the high-viscosity simulations better match the overall observed period ratios than the low-viscosity simulations, we suggest that the majority of the resonance chains are most probably formed at low viscosity, most likely without nearby external perturbers, as these tend to destroy the resonance chains effectively. In contrast, the majority of the systems in non-resonant configurations are likely to have formed in high-viscosity environments with a small number of external perturbers that allow a breaking of the chains rather than their total destruction. This interpretation indicates that a variety of viscosities (and outer companions) is probably needed to explain the different observed systems. Furthermore, this result is independent of the exact formation mechanism of the planets (e.g. pebbles or planetesimals), and of the initial embryo configuration, as long as the embryos can migrate and form resonance chains at the inner edge of the disc.

Acknowledgements

We thank the referee for their comments that helped us to improve our manuscript.

Appendix A Migration of single planets

We show in Fig. A.1 the semi-major axis evolution of 10 Earth mass planets as a function of time. The planets start at either 1 or 10 AU and are embedded in discs with different viscosities. The planets at low viscosity migrate the slowest, while planets at higher viscosity migrate fastest, due to partial gap opening at low viscosity. Consequently, we expect the planets in the high-viscosity regime to result in tighter configurations due to the larger mutual migration velocities (see Appendix G).

Fig. A.1 Single migrating planets starting at 1 AU (left) or 10 AU (right) in discs with different viscosities. The viscosity of α = 10−4 corresponds to our nominal simulations, while the viscosity of α = 5 × 10−3 corresponds to the viscosity used in Izidoro et al. (2021).

Appendix B Migration rates

We show in Fig. B.1 the migration rate of 5 and 10 Earth mass planets in a ficticious disc model with $${\mathrm{\Sigma }}_{\mathrm{g}}={\mathrm{\Sigma }}_0{\left(\frac{r}{\mathrm{A}\mathrm{U}}\right)}^{-3/2},T=170{\left(\frac{r}{\mathrm{A}\mathrm{U}}\right)}^{-1/2}\mathrm{K},$$(B.1)

where Σ₀ varies from 100 to $3400\frac{\mathrm{g}}{{\mathrm{c}\mathrm{m}}^2}$. We additionally vary α from 10⁻⁴ to 8 × 10⁻³. The migration rates are influenced by the gas surface density as well as by the viscosity.

Fig. B.1 Migration rates in AU per Myr of a 5 Earth mass planets (left) and 10 Earth mass planets (right) at 1 AU in a protoplanetary disc following the migration rate formula of Paardekooper et al. (2011) including a reduction of the migration rate by gap opening via the formalism of Kanagawa et al. (2018).

As the gas surface density decreases, the migration rate decreases, because the torque acting on the planet is proportional to the gas surface density. A reduction by a factor of 3 in gas surface density therefore reduces the migration speed of the planets by a factor of 3 as well.

In contrast, as the viscosity decreases, the planet can start to open a gap in the protoplanetary disc, which influences the migration velocity of the planet, where we follow the approach by Kanagawa et al. (2018). Kanagawa et al. (2018) relate the type-II migration time-scale to the type-I migration time-scale (which we calculate as explained above) in the following way $${\tau }_{\text{migII }}=\frac{{\mathrm{\Sigma }}_{\text{up }}}{{\mathrm{\Sigma }}_{\text{min }}}{\tau }_{\text{migI}},$$(B.2)

where Σ_up corresponds to the unperturbed gas surface density and Σ_min to the minimal gas surface density at the bottom of the gap generated by the planet. The ratio Σ_up/Σ_min can be expressed through $$\frac{{\mathrm{\Sigma }}_{\mathrm{u}\mathrm{p}}}{{\mathrm{\Sigma }}_{\text{min }}}=1+0.04K_{\text{mig }},$$(B.3)

where $$K_{\text{mig }}={\left(\frac{M_{\mathrm{P}}}{{\mathrm{M}}_{\odot }}\right)}^2{\left(\frac{H}{r}\right)}^{-5}{\alpha }^{-1}.$$(B.4)

Clearly, a smaller α values allows a faster opening of a gap and thus a reduction in the migration speed, as long as the gas surface density remains constant.

Low gas surface densities reduce the migration speeds and prevent the build up of chains, which is commonly assumed in the in-situ formation scenarios (e.g Hansen & Murray 2013; Dawson et al. 2015; Lee & Chiang 2016). In these scenarios, the planets are normally fully formed before they start to migrate in the low gas surface density environments towards the end of the disc’s lifetime. However, if the discs were very low in surface density at all time, it would be difficult to form planets in the first place due to the lack of planetary building blocks. If the disc had larger gas surface densities initially, the planets would form efficiently, but would then also start to migrate efficiently, resulting in the formation of planetary chains as in Izidoro et al. (2021).

On the other hand, a reduction in viscosity avoids this issue. In fact a lower viscosity is even benefitial for the growth of planets due to the lower pebble scale height (e.g. Lambrechts & Johansen 2014; Bitsch et al. 2018a). In the viscous accretion scenario, a lower viscosity would result in lower accretion rates of the disc, in potential conflict with observations (e.g. Hartmann & Bae 2018; Somigliana et al. 2024). On the other hand, the accretion rates of discs could be driven by disc winds (e.g. Lesur et al. 2023), which still allow a low mid-plane viscosity while maintaining large disc accretion rates.

We note that the opening of a gap also depends on the disc’s aspect ratio, where lower aspect ratios allow an easier opening of a gap. While the disc’s aspect ratio scales with the power of −5, it’s value normally does not change by large quantities. In our used disc model (Bitsch et al. 2015a), the aspect ratio changes by about a factor of 2 within the first Myr of disc evolution. This would then amount to a similar change of K_mig as reducing the viscosity by a factor of 50. However, we implant the embryos into a disc that has already evolved for 2 Myr. At this stage the changes of the disc’s aspect ratio (especially in the inner disc regions) are minimal. This is caused by the fact that the temperature of the disc is determined at these late stages by stellar heating, which also only changes minimally during these stages (Baraffe et al. 2015), resulting in a nearly constant aspect ratio in our model during planet formation (Bitsch et al. 2015a). We thus think that evolution of migrating planets is mainly influenced by changes in the viscosity rather than the disc’s surface or temperature.

Appendix C Low viscosity systems

We show in Fig. C.1 the configurations of the planetary systems after 3 Myr (end of the gas–disc phase) as well as after 100 Myr. The planetary systems become unstable after the gas–disc phase, which reduces the number of planets. We note that for the analysis presented in the main paper only planets more massive than 1 Earth mass are taking into account, because smaller planets are not reliably detectable.

Fig. C.1 Systems formed in our low-viscosity simulations after the gas–disc phase (left) and after 100 Myr of integration (right). The size of the dots is proportional to the planetary mass, while the arrow marks the aphelion and perihelion of the corresponding planet. We note that the synthetic observations only take planets up until 0.7 AU into account.

Appendix D Influence of embryo number and disc viscosity

The total number of planetary embryos can – of course – influence the outcome of the simulations, especially when looking at the limiting case of just a few planetary embryos, which might even be placed too far apart to gravitationally influence each other. While it is beyond the scope of this work to run an extensive suite of simulations with different number of initial embryos, we relate back to our previous work (Bitsch & Izidoro 2023). In this work, we studied the formation of inner super-Earth systems with outer giant planets in the same set-up as presented here – with the exception that gas accretion onto planetary embryos was allowed so that gas giants can form.

In Fig. D.1 we plot the period ratios of the systems within 1.0 AU for all planets in the simulations after the end of the gas–disc phase at 3 Myrs with initially 15, 30 or 60 planetary embryos and with a K = 5 damping for the giant planets. We use here the full sample of planets, because some planets in the simulations of Bitsch & Izidoro (2023) grew envelopes and became larger than 20 Earth masses – the limit used in the simulations by Izidoro et al. (2021) as a definition of a super-Earth. We are here primarily interested in the influence of the number of embryos on the period ratios and not on their mass distribution.

Fig. D.1 Period ratios of the planets interior to 1.0 AU in the simulations of Bitsch & Izidoro (2023) with K = 5 at the end of the gas–disc phase at 3 Myr. We note that the planets in the 15 embryo simulations mostly turn to gas giants, reducing their migration speeds and thus allow wider period ratios.

Using 30 initial embryos results in a similar period ratio distribution compared to the results presented in Fig. 4, namely a large fraction of systems in the 4:3 and 3:2 configuration. Using 60 initial planetary embryos results in a similar pile-up in the 4:3 and 3:2 configurations. This is in contrast to the simulations by Izidoro et al. (2021), which show a pile-up in tighter configurations. As Izidoro et al. (2021) used 60 initial planetary embryos with the same starting configuration as our simulations, we can conclude that the viscosity is indeed responsible for the differences in the period ratios rather than the number of initial embryos.

On the other hand, if only 15 initial embryos are used, the situation is different. Here we see a pile-up mostly in the 2:1 resonance configuration, rather than the 3:2 and 4:3. This is caused by the wide separation of the planetary embryos, where the embryos are spaced initially around 1.0 AU apart (Bitsch & Izidoro 2023). Consequently only 1–2 planets reach the inner edge of the disc, where the outer planet already transitioned into a gas giant, which reduced their migration speed significantly so that the 2:1 configuration cannot be overcome. Furthermore, the large initial separation prevents the second planet from reaching the inner edge such that in most cases no resonant system is formed.

In Fig. D.2 we show the period ratios of our simulations with α = 10⁻⁴, as discussed in the main part of the paper, as well as from comparison simulations using the same disc and planet parameters, but α = 5 × 10⁻³. Both sets of simulations feature 30 initial embryos. We show the periods at the end of the gas disc lifetime of 3 Myr. Clearly, the planetary systems formed in the high-viscosity simulations result in tighter resonant chains compared to the systems formed in the low-viscosity environments. In particular, the high-viscosity simulations feature a larger fraction of planets in the 4:3 and 5:4 (and tighter) resonance configurations, while the low-viscosity simulations produce many more planets in the 3:2 and 2:1 resonance configurations. We note that the exact period configurations are subject to change after the instability phase, but instabilities break the resonances and do not populate the resonance configurations, indicating that overall trend that higher viscosity systems produce resonant chains in tighter configurations compared to low-viscosity simulations will continue to hold.

Fig. D.2 Period ratios of the planets interior to 1.0 AU for our simulations using α = 10−4 and a comparison sample with the same disc properties, but α = 5 × 10−3. The resulting periods are plotted after the end of the gas disc phase at 3 Myr.

We can conclude from these results that as long as enough planetary embryos are present, high viscosities are needed to make tight chains, while lower viscosities are needed to make wider chains (e.g. for the 3:2 rather than 5:4 configuration). However, if the number of planetary embryos is very limited, even wider configurations are possible. Here, we have to admit the caveat that this might be caused by the larger planetary masses for the simulations with 15 initial embryos that allow deeper gap opening in the simulations of Bitsch & Izidoro (2023) compared to our here presented simulations.

Appendix E Planetary masses

The planetary mass determines not only the migration speed, but are also important for the stability of planetary systems (e.g. Gladman 1993; Chambers et al. 1996). In Fig. E.1 we show the cumulative mass distribution of our planets after 100 Myr of evolution. Clearly, the masses of our low-viscosity planetary systems are too low compared to the inferred masses from the Kepler observations (using the mass-radius relationship of Wolfgang et al. (2016)). The high-viscosity simulation of Izidoro et al. (2021), on the other hand, match very well with the observations, especially for the masses of the chains.

In our model, the only difference between the high viscosity and low-viscosity simulations is the migration velocity. The viscosity enters also into the pebble isolation mass (Bitsch et al. 2018b), but is set to the same value in all our simulations to allow an easier comparison. Thus, the difference in mass only arises from the difference in migration velocity. As shown in Fig. A.1, the planets formed in low-viscosity environment migrate slower and can not reach the inner edge of the disc towards the end of the disc’s lifetime, if they start to migrate from exterior to ≈ 6 − 7 AU (see also Brasser et al. 2017). However, in this region, the pebble isolation mass starts to increase dramatically (see right in Fig. E.1). Therefore the majority of the planets that form via pebble accretion in the low-viscosity environments that migrate to the inner disc only reach up to ≈3 Earth masses. As pebble accretion is very efficient, nearly all planets reach this mass, explaining the lower mass tail in the left of Fig. E.1. The larger planetary masses are caused by collisions after the gas disc phase during the instability phase of the systems.

Fig. E.1 Top: Cumulative distribution of the planetary masses in our simulations (red) and their detections (black) as well as from the simulations of Izidoro et al. (2021) (purple) after 100 Myr of integration time. We also overplot the masses from the Kepler observations using the mass-radius relationship of Wolfgang et al. (2016). Bottom: Time evolution of the pebble isolation mass in our model, following the recipe of Bitsch et al. (2018b). Here t = 0 corresponds to the starting time of the simulations, so to a disc that is already 2 Myr old in the framework of Bitsch et al. (2015a).

In contrast, the planets formed in simulations with higher viscosity can migrate from much further out into the inner disc region due to their faster migration velocity (Fig. A.1). This result in a steady influx of planets that have reached higher masses due to the higher pebble isolation mass exterior to a 6-7 AU. Furthermore, these planets are massive enough that they can collide already during the gas disc-phase (Izidoro et al. 2021), which is not the norm for the planets formed in the low-viscosity simulations. Consequently, the planets formed in the high-viscosity simulations in the inner disc feature higher planetary masses, even though the formula for the pebble isolation mass is the same.

Appendix F Period ratios

We show in Fig. F.1 a mixture of 2% stable chains and 98% systems that underwent instabilities. This mixture was used in the past high-viscosity simulations of Izidoro et al. (2021), as it gave the best fit to the observations (indicated by the purple line). The low-viscosity simulations still fail to reproduce the observed period ratios, as their period ratios remain too tight.

Fig. F.1 Number of synthetically observed planets from our simulations (top) and the period ratios of adjacent planets (bottom) for a mix of planetary systems that includes 2% stable chains and 98% systems that underwent instabilities. The colors are the same as in Fig. 4. Both sets of simulations give a nice match to the number of detected planets, but the period ratios are matched better in the high-viscosity case, compared to the low-viscosity case.

Appendix G Resonant chains of two-planet systems

In order to investigate the influence of viscosity on the period ratio of planetary systems, we conduct a simple experiment. We take systems with 2 equal mass planets and let them migrate through a discs, where only the viscosity varies between the different set-ups. Different viscosities can cause different migration velocities (see Appendix A) due to the change of the strength of the entropy driven corotation torque (e.g. Paardekooper et al. 2011; Bitsch et al. 2013). In particular, high-viscosity environments can lead to outward migration, if the disc’s surface density and temperature gradients are in the right regime. In fact, high viscosities will lead to outward migration in our used disc model (Bitsch et al. 2015a), which complicates our simple experiment. We thus use a simple power law disc model for this experiment: $${\mathrm{\Sigma }}_{\mathrm{g}}=1700{\left(\frac{r}{\mathrm{A}\mathrm{U}}\right)}^{-3/2}\frac{\mathrm{g}}{{\mathrm{c}\mathrm{m}}^2},T=170{\left(\frac{r}{\mathrm{A}\mathrm{U}}\right)}^{-1/2}\mathrm{K}.$$(G.1)

Fig. G.1 Final period ratios of two-planet systems with equal masses in discs with different viscosities.

This simple disc model does not support outward migration for all levels of viscosity, allowing us to study clearly the migration behaviour of pairs of planets of equal mass, which are not growing. We start our planetary systems just outside the 2:1 resonance in order to study how the different migration velocities influence the final resonance configuration. We use here the same migration prescription as in our main paper (Paardekooper et al. 2011).

We show in Fig. G.1 the final period ratios of our planetary systems as function of the disc’s viscosity. An increase in the α-viscosity parameter results in tighter resonance configurations, caused by the larger migration velocities of the planets in the higher viscosity environments. Furthermore, lower planetary masses result in tighter resonances as well, compared to higher mass planets. The reason for this is the strength of the resonance, which scales with the planetary masses, and needs to be overcome by migration.

The here observed trend is reflected in our main simulations, where the simulations at low viscosity with low masses show wider final period ratios compared to the simulations with higher masses and higher viscosities. We note though that the exact period ratio in a multi-planet system is more complicated as in the two-planet case, due to the increase in eccentricities caused by the gravitational interaction of multiple planets that can influence the exact resonance configuration.

References

All Figures

Fig. 1 Evolution of a single planetary system as a function of time regarding semi-major axis (top left), mass (top right), eccentricity (bottom left), and inclination (bottom right). The coloured lines correspond to the surviving inner (r < 4 AU) planets, while the grey lines depict the evolution of the outer (also less massive) planets. The vertical line at 3 Myr marks the dissipation of the gas disc phase, where the planets lose the damping forces on eccentricity and inclination from the disc. The inner system remains stable until 100 Myr. The system corresponds to system number 3 in Fig. C.1.

In the text

Fig. 2 Same as Fig. 1, but for a system that became unstable after the gas–disc phase with two surviving inner planets. We mark all surviving planets in colour. The planet marked in yellow stops growing at around 700 kyr, when it acquired an eccentricity too large to allow the accretion of pebbles. Once its eccentricity is damped below 1%, it continues to grow via pebble accretion. The system corresponds to system number 21 in Fig. C.1.

In the text

	Fig. 3 Same as Fig. 1, but for a system that became unstable after the gas–disc phase with no surviving inner planets. We mark all surviving planets in colour. The planet marked in orange acquires some eccentricity during the gas–disc phase, preventing pebble accretion, and it only continues to grow once its eccentricity is damped below 1%. The system corresponds to system number 1 in Fig. C.1.
In the text

Fig. 4 Number of synthetically observed planets (left) and period ratios of adjacent planets (right) for the planetary systems directly at the very end of the gas–disc phase, before the systems undergo instabilities. The thin grey line marks all the observations, while the thick grey line marks the observed chains of planets. The red colour corresponds to the direct results of our simulations, while the black line corresponds to the synthetical observations of our simulations. The purple line represents the synthetic observations from Izidoro et al. (2021). A clear difference is noted in the compactness of the chains between the high-viscosity simulations of Izidoro et al. (2021) and the low-viscosity simulations used here.

In the text

Fig. 5 Number of synthetically observed planets (left) and period ratios of adjacent planets (right) for the planetary systems after 100 Myr of integration; i.e. after the instability phase. The colours are the same as in Fig. 4. The high-viscosity simulations from Izidoro et al. (2021) produce wider systems after the instabilities compared to the low-viscosity simulations, which is also a consequence of the real number of planets in the systems before instabilities (left in Fig. 4), where more planets are initially present at high viscosities.

In the text

	Fig. 6 Number of synthetically observed planets (left) and period ratios of adjacent planets (right) for chains of planets that are stable after 100 Myr. The colours are the same as in Fig. 4.
In the text

	Fig. 7 Eccentricity distribution of our low-viscosity simulations and of the high-viscosity simulations of Izidoro et al. (2021) for all planets (left) and only for the surviving chains (right).
In the text

	Fig. A.1 Single migrating planets starting at 1 AU (left) or 10 AU (right) in discs with different viscosities. The viscosity of α = 10−4 corresponds to our nominal simulations, while the viscosity of α = 5 × 10−3 corresponds to the viscosity used in Izidoro et al. (2021).
In the text

	Fig. B.1 Migration rates in AU per Myr of a 5 Earth mass planets (left) and 10 Earth mass planets (right) at 1 AU in a protoplanetary disc following the migration rate formula of Paardekooper et al. (2011) including a reduction of the migration rate by gap opening via the formalism of Kanagawa et al. (2018).
In the text

	Fig. C.1 Systems formed in our low-viscosity simulations after the gas–disc phase (left) and after 100 Myr of integration (right). The size of the dots is proportional to the planetary mass, while the arrow marks the aphelion and perihelion of the corresponding planet. We note that the synthetic observations only take planets up until 0.7 AU into account.
In the text

	Fig. D.1 Period ratios of the planets interior to 1.0 AU in the simulations of Bitsch & Izidoro (2023) with K = 5 at the end of the gas–disc phase at 3 Myr. We note that the planets in the 15 embryo simulations mostly turn to gas giants, reducing their migration speeds and thus allow wider period ratios.
In the text

	Fig. D.2 Period ratios of the planets interior to 1.0 AU for our simulations using α = 10−4 and a comparison sample with the same disc properties, but α = 5 × 10−3. The resulting periods are plotted after the end of the gas disc phase at 3 Myr.
In the text

Fig. E.1 Top: Cumulative distribution of the planetary masses in our simulations (red) and their detections (black) as well as from the simulations of Izidoro et al. (2021) (purple) after 100 Myr of integration time. We also overplot the masses from the Kepler observations using the mass-radius relationship of Wolfgang et al. (2016). Bottom: Time evolution of the pebble isolation mass in our model, following the recipe of Bitsch et al. (2018b). Here t = 0 corresponds to the starting time of the simulations, so to a disc that is already 2 Myr old in the framework of Bitsch et al. (2015a).

In the text

Fig. F.1 Number of synthetically observed planets from our simulations (top) and the period ratios of adjacent planets (bottom) for a mix of planetary systems that includes 2% stable chains and 98% systems that underwent instabilities. The colors are the same as in Fig. 4. Both sets of simulations give a nice match to the number of detected planets, but the period ratios are matched better in the high-viscosity case, compared to the low-viscosity case.

In the text

	Fig. G.1 Final period ratios of two-planet systems with equal masses in discs with different viscosities.
In the text