© The Authors 2025

1 Introduction

Planetary population synthesis is a method used to put different theories of planet formation to the quantitative observational test by statistically comparing their predictions with the observed exoplanet population (e.g. Ida & Lin 2004; Mordasini et al. 2009a). It has become possible thanks to the fact that over the last three decades since Mayor & Queloz (1995), large surveys employing different observational techniques have discovered thousands of exoplanets as well as sub-stellar companions. Important examples of radial velocity surveys include the HARPS survey (Mayor et al. 2011) and the California Legacy survey (Rosenthal et al. 2021); of the transit method, the Kepler satellite (Borucki et al. 2010; Petigura et al. 2018); direct imaging surveys such as GPIES (Nielsen et al. 2019), SPHERE SHINE (Vigan et al. 2021), and BEAST (Janson et al. 2021b); and finally, microlensing (Suzuki et al. 2016).

Observations of planets as they form remain challenging, and to date, only a handful of bona fide examples are available (Keppler et al. 2018). As a result, the process of planet formation cannot be directly constrained from observations. The limited knowledge about most individual exoplanets (such as some orbital elements and a mass or radius) further affects our ability to construct planet formation models, because they cannot be easily confirmed or rejected through these approaches.

On the other hand, statistical approaches can be used to better constrain planet formation theory. They make use of the large number of observations of young stellar objects (YSOs) and, in particular, their discs (Andrews 2020) in which planets and other companions form. The observations of discs and YSOs mark the initial and boundary conditions for formation. Observed (exo)planetary systems with their statistical properties (Udry & Santos 2007) represent the final outcome of the formation process. This is the principle of planetary population synthesis (Benz et al. 2014; Burn & Mordasini 2024), which tries to find a theory that can best reproduce the observed demographics of planetary systems.

The most important element of the population synthesis approach is a global formation model (Mordasini et al. 2015; Raymond & Morbidelli 2022), which can produce synthetic populations. It is based on an underlying formation paradigm and combines the results of many models of specific physical processes such as disc formation and evolution, gas accretion, collisions, and orbital migration.

Two fundamentally different planet formation paradigms have been proposed. Disc instability (DI), also known as gravitational instability (GI), was the first proposed for the formation of the Solar System (Kuiper 1951). It was later mostly superseded by the core accretion (CA) model (Safronov 1972; Pollack et al. 1996; Alibert et al. 2005) as the main formation pathway of Solar System planets as well as observed exoplanets. However, the DI model is still a valid alternative in some cases and has remained the favoured explanation for the formation of massive giant planets at large separation. The formation of these systems is hard to explain by CA. A famous example of such a system might be the HR 8799 system (Marois et al. 2008, Marois et al. 2010).

While not studied as much in detail as CA, a number of important results concerning the outcome of DI exist. A fundamental prediction is that DI should, in particular, form brown dwarfs (BDs). For example, Stamatellos & Whitworth (2009) performed radiation hydrodynamics simulations of fragmenting discs and found that the most likely outcome is companions in the BD mass range with high eccentricity and inclination, although stellar and planetary-mass companions were also formed. Kratter et al. (2010) studied fragmentation and the subsequent evolution of fragments and investigated whether wide companions such as those in the HR 8799 system could form through DI. They explored the influence of the initial fragment mass, gas accretion, and gap formation on the final masses and found that the formation of an object with a final mass in the planetary regime is rather unlikely. Instead, masses in the BD or stellar regime were found to be the most probable outcome.

Helled & Bodenheimer (2010) studied the composition of companions formed in DI. The influence of a number of factors, most notably the pre-collapse timescale, was considered in planetesimal capture through clumps. The authors found that, if formed in DI, the HR 8799 planets would have low metallicities.

Zhu et al. (2012) performed 2D hydrodynamical simulations of self-gravitating discs including infall. They studied the accretion of gas and the migration of the clumps formed, and concluded that many are lost by Type I migration to primary or tidal destruction, although gap formation is possible in some cases due to gas accretion.

Boley et al. (2010) studied the initial conditions for clumps formed in DI. They provided an analytic estimate for the initial clump mass that places this quantity in the gas giant regime. This estimate was confirmed in radiation hydrodynamics simulations. It was found that clumps could form a heavy-element core and have a tidally stripped envelope, thus potentially leading to the formation of rocky planets.

Baehr & Klahr (2019) performed 3D self-gravitating shearing sheet simulations in order to study solid accretion by clumps. They find that the formation of cores of up to several Earth masses is possible. The process will depend on the metallicity of the disc. Core formation leads to a solid depletion of the clump’s atmosphere, whereas if no core forms, a high metallicity atmosphere is expected. For a more thorough overview, we refer to the reviews in Durisen et al. (2007); Mayer (2010); Helled et al. (2014); Kratter & Lodato (2016).

The population of companions formed via DI is currently unknown, and it remains unclear whether some observed companions formed in this way. In order to answer this question, it is necessary to perform a population synthesis study that links observationally constrained initial conditions to a population of surviving companions, considering all relevant physical processes. Various previous studies have aimed at establishing a population synthesis model for DI, but they have limitations (see Sect. 2.3 for discussion). A key limitation is that previous studies did not include the earliest formation phase where the star-and-disc system forms during infall. This phase is crucial because this is when the disc is most massive and prone to fragmentation. Another key limitation of existing studies is that the initial conditions are not statistically linked to observations. Especially, the fraction of systems that fragment is important to determine. This fragmentation fraction directly translates into the fraction of systems that can potentially host a companion formed through DI. In this paper, we present a comprehensive population synthesis model that aims to overcome these issues. We applied the model described in Schib et al. (2025) (Paper I) with initial conditions informed by observed quantities such as stellar masses and disc sizes, as well as results from hydrodynamic simulations. This allows us to make predictions about the properties of formed systems and their occurrence rates. For example, the fraction of companions in a given interval of mass and separation can be predicted as a function of the primary mass. Such predictions enable a quantitative comparison to observational surveys (e.g. Nielsen et al. 2019). Furthermore, it is possible to estimate the contribution of DI to the observed exoplanet population and guide future observations as facilities improve (e.g. Lagrange et al. 2025).

This paper is organised as follows: In Sect. 2, we give a brief overview of the population synthesis method with a focus on existing projects in the DI paradigm. Section 3 briefly summarises the model. The initial conditions for the population synthesis are described in Sect. 4. In Sect. 5, we present the results. Section 6 is dedicated to the variant populations we performed in order to study a number of parameters. Section 7 compares results from previous studies. In Sect. 8, we compare with observations. Section 9 contains a discussion of our results and limitations of the model. We summarise and conclude in Sect. 10.

2 Planetary population synthesis

2.1 The population synthesis method

Planetary population synthesis requires several ingredients (Benz et al. 2014; Emsenhuber et al. 2023; Burn & Mordasini 2024), including demographics of YSOs and their discs for the initial conditions, a formation model to predict an exoplanet population, and exoplanet demographics to compare with the model results. The model must meet specific requirements: it must be global (i.e. including all the known relevant physical processes and them being coupled), end-to-end (i.e. going from YSOs to mature systems, predicting as many directly observable properties as possible), but still remain relatively fast computationally such that a full population with thousands of systems can be computed.

In the CA paradigm, these criteria are met by a number of existing studies, where models have reached a significant level of complexity. Thus, there are detailed predictions about the expected population of CA planets. Conversely, only a few GI population synthesis projects exist. They are less comprehensive in comparison to their Ca counterparts and have received less attention in the literature. This reflects that global models can only be constructed if sufficient specialised models on the individual relevant physical processes exist. Below, we briefly review some of the existing CA projects, followed by a more detailed description of published GI-population studies.

2.2 CA population synthesis models

The era of planetary population synthesis began with the work of Ida & Lin (2004). The authors developed a simple global model of planet formation and used it to study the mass-semi-major axis distribution of planets predicted by CA. They included stars of different masses and metallicities, but neglected Type I orbital migration. Since then, CA population synthesis and global models have been improved and include more processes such as the interior structure, N-body interactions between planets, and the long-term evolution of planets. In addition, more detailed comparisons with the increasingly large observational surveys have been performed. Key studies include Ida & Lin (2008b, a); Thommes et al. (2008); Mordasini et al. (2009a, b); Miguel et al. (2011b, a); Alibert et al. (2011); Hellary & Nelson (2012); Alibert et al. (2013); Ida & Lin (2010); Ida et al. (2013); Coleman & Nelson (2014, 2016) and references therein. The most extensive population synthesis project to date is the Bern Model (Mordasini et al. 2009a, b; Alibert et al. 2013; Emsenhuber et al. 2021; Schlecker et al. 2021a; Burn et al. 2021; Schlecker et al. 2021b; Mishra et al. 2021; Emsenhuber et al. 2023). Further details on planetary population synthesis models in the CA paradigm can be found in several review papers, including Benz et al. (2014); Mordasini (2018); Drazkowska et al. (2023); Emsenhuber et al. (2023); Burn & Mordasini (2024).

2.3 DI population synthesis models

Although GI has not been studied as extensively as CA, a few DI population synthesis models have been published. Two key models are the ‘Forgan & Rice model’ (Forgan & Rice 2013; Forgan et al. 2018) and the ‘Nayakshin model’ (Nayakshin 2015a; Nayakshin & Fletcher 2015; Nayakshin 2016). The existing studies have rather different foci, and they also differ significantly in the choice of initial conditions and the physical processes that are included, as we discuss below. A table that compares the most important features of the different projects is presented in Appendix A.

2.3.1 Forgan and Rice model

D. Forgan and K. Rice published a series of papers on selfgravitating discs and fragmentation. Their first population synthesis model was presented in Forgan & Rice (2013). Their population focused on solar mass stars and used pre-evolved discs to study the outcome of disc fragmentation. Fragments with an initial fragment mass M_J,FR (Sect. 6.2) were inserted at random locations within the disc and were allowed to migrate. The fragments could migrate and open gaps. Grain growth and settling were also included, while gas accretion was not considered. The authors found that the majority of surviving companions are BDs, and about 40% form heavy-element cores. The formation of terrestrial planets was found to be very unlikely. In Forgan et al. (2018), several adjustments to the migration prescription were made, and more importantly, gravitational interactions between fragments were included using N-body simulations (though mergers were not considered). It was found that gravitational interactions dominate the orbital evolution, leading again to a population dominated by BD, but with a much wider spread in semi-major axes. Some of the fragments with an evolved grain content were found to be disrupted, possibly leading to the formation of planetesimal belts. The authors also predicted a population of free-floating objects due to ejections.

2.3.2 Nayakshin model

S. Nayakshin and collaborators presented a population synthesis study focused on tidal downsizing, which is the fragmentation of the gaseous disc to form bound clumps followed by solid accretion by the clumps, core formation, and tidal stripping of the gas envelope (Nayakshin 2015a; Nayakshin & Fletcher 2015; Nayakshin 2016). These models focused on solar mass stars and used one companion per system. Gas accretion was not considered, but instead, fragments could accrete pebbles. Planet-disc interaction was modelled by means of a simple Type I migration timescale, and by the impulse approximation in the gap-opening regime. The authors found a population of companions with masses from 10^–3 M_♃ to 2 × 10¹ M occupying the region from 10^–1 au to 10² au. The population includes rocky planets as well as gas giants with a range of bulk metallicities.

2.3.3 Other studies

Here we briefly mention a few other studies that aimed to understand the planetary population predicted by the DI paradigm. Baruteau et al. (2011) performed 2D hydrodynamic simulations of gravito-turbulent discs with prescribed cooling. They studied the migration of a single clump around a 1 M_⊙ star and found it to migrate inwards very quickly (~ 10 kyr) without forming a gap, independent of clump mass. Rowther & Meru (2020) studied a very similar setting, but included a prescription for cooling that has a longer timescale in the inner disc. It was found that this more realistic cooling allows the clumps to open a gap in the disc. As a result, clumps are saved from catastrophic inward migration.

Müller et al. (2018) explored the migration of single clumps in static (after fragmentation) discs around predominantly solar mass stars. They investigated different gap opening criteria, initial clump masses, and clump density profiles. Gas accretion and tidal mass loss were also considered. They find that surviving clumps on wide orbits are rare, and that the inferred population depends strongly on the model assumptions.

Humphries et al. (2019) performed simulations of migrating clumps formed in DI, aiming at constraining their migration using constraints from radial velocity surveys of massive giants inside 5 au. They simulated single evolving clumps around solar type stars, neglecting accretion of gas or solids. Thermal irradiation and tidal disruption were considered. The authors found that tens of percent of systems need to fragment, and clumps need to migrate quickly to the inner disc, in order for their constraints to be met.

Jin et al. (2020) conducted a population synthesis study using a 1D model with a simple infall model. They varied the properties of the parent molecular cloud core (MCC; temperature, mass, and rotation) and investigated the influence on the initial fragment mass as well as the location and time of fragmentation. Applying the Toomre mass (Nelson 2006), they inferred initial masses typically 3 M_♃ to 35 M_♃ and fragmentation locations between 20 au to 200 au.

3 Model

The model used in this work was developed, applied, and extended in Schib et al. (2021) (S21), Schib et al. (2022) (S22), and Schib et al. (2023) (S23). The updated version used in the present work is described in detail in Paper I. Thus, we only give a brief summary here.

The model describes the fragmentation of discs around single stars. In each simulation, the star-and-disc system is initialised at a very small mass and then grows from infalling material from the MCC. To this end, the 1D viscous evolution equation for the disc’s radial structure is solved numerically, yielding the disc’s surface density (Lüst 1952; Lynden-Bell & Pringle 1974). The disc’s vertical structure and temperature are calculated based on a thermal balance criterion, including the effects of accretion heating, viscous heating, stellar irradiation, and radiative cooling (Nakamoto & Nakagawa 1994; Hueso & Guillot 2005,S23).

The infalling material is added to the disc, realised as a source term in the equation. Matter accreting across the fixed inner edge of the disc is added to the star. Clumps (young fragments, bound blobs of gas) emerge when the conditions for fragmentation are satisfied. Their initial mass is assumed to be M_F (Boley et al. 2010, Eq. (35) in Paper I). The corresponding gas is removed from the disc. Clumps evolve according to pre-calculated evolution tracks of isolated 1D objects, corrected by disc irradiation (Helled et al. 2006, 2008; Vazan & Helled 2012; Vazan et al. 2013, 2015; Humphries et al. 2019). The evolution tracks include the second collapse of the clumps due to the dissociation of molecular hydrogen (Bodenheimer 1974). If this happens, the clump sizes shrink by several orders of magnitude, and the objects become compact gas spheres. These are gravitationally bound and can end up as giant planets, BDs, or stars.

Clumps can gain mass by gas accretion and collisions and lose mass due to stellar tides. They can also be completely disintegrated by tidal or thermal disruption. Any mass that is lost by the clump is returned to the disc, where the exchange of mass between companions and the disc is modelled with a source/sink term in the evolution equation to ensure the conservation of mass.

Companions migrate and open gaps via two-way exchange of angular momentum, which is realised through torque densities in the disc surface density evolution equation (S22). Damping is calculated based on dynamical friction (Ida et al. 2020), and the gravitational interaction between companions is modelled using N-body integration (Chambers 1999). Clumps can become unbound and leave the planetary system, or be accreted by the central star. Companions can also collide, and in this work, we assume that collisions lead to perfect merging (see discussion in Sect 9.4).

The systems were evolved until the disc disappeared. If any companions remained at this point, we continued the N-body integration up to 100 Myr.

The model thus includes a significant number of physical processes that interact and feed back on each other in multiple ways. To still be able to simulate a high number of systems from t=0 to 100 Myr, the model is low-dimensional, as typical for global models: the disc is 1+1D (radial and vertical) while the clumps/companions are 1D spherically symmetric. Only the N-body is fully 3D.

4 Initial conditions

To perform a population synthesis, the initial conditions for the simulations must be defined. This includes, on the one hand, parameters that are fixed for all simulations, and on the other hand, Monte Carlo variables that are drawn from given probability distributions. Table 1 lists the most important fixed parameters (top) as well as the Monte Carlo variables (bottom). For the latter, the mean value is given. The quantities are explained in the following.

4.1 Overall setup

In planetary population synthesis for the CA paradigm, setting initial conditions typically means defining the initial (gas) disc mass and the (fixed) stellar mass. The outer radius and profile of the gas disc are then defined based on observational constraints and theoretical considerations. During the simulation, the disc’s mass can only decrease, and the stellar mass remains constant (for example, Sect. 3 in Emsenhuber et al. 2021).

Since for DI we want to study the earliest stages of planet formation when discs are most massive, we needed a different approach. Our star-and-disc systems are initialised very early, at 1 kyr after the formation of the central protostar. The masses of the star and the disc at this time are of the order of 0.01 M_⊙. This is very low compared to their maximum value reached during the evolution. The large majority of the mass is added to the system by infall, where all the infalling material is added to the disc. The star’s mass is increasing as the disc gas accretes onto the star across the inner edge. During the simulation, the masses of the star and disc increase significantly, in contrast to CA simulations. Knowledge of infall rates at the earliest moments of star formation is therefore required. In addition, we required some information about the sizes of young discs in order to identify the locations where gas infalls onto the disc.

4.2 Initial star and disc masses and infall rates

The stellar masses, disc masses, and disc sizes are determined by the source term for the infall. Unfortunately, these quantities and their time evolution are difficult to observe directly and therefore we built upon the numerical results from Bate (2018) where a hydrodynamic simulation of star formation (Bate 2012) was analysed in detail with respect to the discs formed around the forming stars. The original calculation consists of a 500 M_⊙ molecular cloud with a radius of ≈8 × 10⁴ au and uses 3.5 × 10⁷ smoothed particle hydrodynamics (SPH) particles. During the simulation, 183 protostars (represented by sink particles) are formed, and the discs around these protostars are analysed. An important result of the study is the chaotic nature of the star formation process. It shows that the ideal picture of a disc around a single star in isolation is over-simplified and rather rare. Instead, most discs are truncated in star-disc interactions, multiple systems form, and discs can be destroyed and formed again through accretion. It is therefore challenging to extract the required data. For simplicity, we selected 35 single (i.e. not bound to another star) protostars from the SPH simulation. The selection criterion was that these systems are not affected by other protostars for a time span of ≈10kyr.

We used the growth rates of these systems (protostar + disc) in this early phase as infall rates for our simulations. The indices of the systems we chose are given in Appendix E of S21. We assumed that the stellar mass, disc mass, and infall rate extracted this way are correlated, and used them to construct a probability distribution from which we drew the initial stellar mass, disc mass, and infall rate. More details on this procedure are given in Sect. 3 of S21. We also assumed that the infall rates are constant throughout the infall phase. This is a strong simplification, and we discuss it further in Sect. 9.

4.3 Infall locations

A disc is formed because the angular momentum initially present in the MCC must be conserved. Given some assumptions about the initial angular momentum distribution (e.g. a spherical MCC in uniform rotation), the source term for the infall onto the disc can be obtained (semi-)analytically. Such recipes have been applied in the past and there is a wealth of literature on them (e.g. Ebert 1955; Bonnor 1956; Shu 1977; Terebey et al. 1984; Visser & Dullemond 2010). However, the star formation process is more complex than suggested by such simple models. The disc sizes from the SPH simulation can, in principle, be used to determine where the infalling material needs to be added, as was done in S21 (see Sect. 2.10.1 for details). In particular, this implies that the location where the infalling material is added to the disc (R_i) is increasing with time. In S21, we find that this way of determining the infall location leads to discs that are much larger than observed Class 0 disc sizes (Tobin et al. 2020). In order to make this comparison, we compared the gas disc sizes at the end of the infall phase to the observed Class 0 dust disc sizes. Due to their large sizes, almost half of these discs were found to fragment. In S21, we also investigated cases with much smaller infall radii, as suggested by some non-ideal magneto-hydrodynamic (MHD) simulations (Hennebelle et al. 2016). In this case, we found that all discs remain stable against fragmentation, and that the discs are much smaller than observed. In S23, we identified the initial conditions that lead to early disc sizes compatible with observations. We iteratively determined infall locations that lead to early disc sizes compatible with the Class 0 disc sizes from Tobin et al. (2020) and found the efficiency of accretion heating that simultaneously reproduces the observed luminosities. We used infall locations constructed in this way in the present work.

The infall location is of the order of 10 au, with a dependency on the stellar mass, so a different value is used in each of the 100 stellar mass bins (see Sect. 5.2.2). The values used in the baseline population are given in Appendix B.

4.4 Initial orbits

When a clump is introduced into the disc (Sect. 3), its initial orbital elements must be defined. This should be done in a way that reflects the stochastic nature of fragmentation (Paardekooper 2012; Young & Clarke 2015; Xu et al. 2025). We assigned an initial separation at a random location in the region of the disc where Q_Toomre < 2. Since clumps are formed in spiral arms (e.g. Mayer et al. 2004), they are expected to be on initially eccentric orbits. We therefore assigned them a random eccentricity between 0 and 0.3 (strictly speaking, we assigned an angular momentum deficit, as explained in Sect. 2.3 of Paper I). Clumps are expected to have nearly flat orbits (e.g. Kubli et al. 2025), and we therefore assigned them a very small random inclination from 0 to 0.03°. The value of the inclination is not very important, as long as it is not zero (as this would confine all orbits to a plane). Clumps start at apoapsis, with a random argument of periapsis and longitude of the ascending node.

5 Results for the baseline population DIPSY-0

In this section, we present and discuss the main results of our study, which are the fundamental statistical properties of the baseline population DIPSY-0. In this population, we used M_F as the initial fragment mass (Sect. 3). The gas accretion rate onto companions is limited to 10^–3 M_♃/yr (see Sect. 9.5). A random number (from one to three) of clumps are inserted. Early disc sizes are set to be consistent with observed Class 0 discs (Tobin et al. 2020, S23), and the values for the background viscosity and the rate of external photoevaporation are according to Table 1. We discuss several variant populations, where we studied deviations from these assumptions, in Sect. 6.

The present section is organised as follows: in Sect. 5.1, we discuss the possible fates of fragments. In Sect. 5.2, we show the mass-semi-major axis diagram of surviving companions, in Sect. 5.3 the companion mass function (CMF), and in Sect. 5.4 the distribution of semi-major axes. Section 5.5 focuses on the eccentricities and inclinations of the companions, Sect. 5.6 on their multiplicity. In Sect. 5.8, we discuss the dependency of our results on the host stellar mass. Section 5.9 finally looks at collisions and close encounters.

5.1 Fates of fragments

Over the 100 000 systems of the baseline population DIPSY- 0, a total of ≈430 000 fragments form. Only a small fraction of the fragments (7400, 1.7%) survive for 100 kyr.

Figure 1 gives an overview of the possible fragment fates. We find that the vast majority of objects (more than 60%) are lost in collisions (see Sect. 5.9 for a detailed discussion). Another large group of objects (26%) is found to be ejected. In this study, ejected objects correspond to those that move beyond 10000au from their host star. This typically occurs via gravitational interaction with other objects. Another 2% of objects are found to merge with the central star. This occurs either because they migrate all the way across the inner edge of the grid or because they collide with the star via a highly elliptical orbit. The fraction of clumps that are thermally disrupted is found to be 9%. Only a very small fraction of clumps (0.5%) are tidally disrupted.

Fig. 1 Occurrence of the different possible fragment fates for all fragments formed in the baseline population. From left to right: Collision and merger with another clump, ejection from the system, thermal disruption, accretion on the primary, survival of the fragment, tidal disruption.

5.2 Mass-semi-major axis diagram

The 100 000 systems were divided into 100 logarithmically spaced bins in final primary mass, ranging from 0.05 M_⊙ to 5 M_⊙. This is necessary due to the strong mass dependence of the initial mass function (IMF) of stars, as we discuss in Sect. 5.2.2. A key result of population syntheses is the prediction of the mass-distance diagram, which might be of similar importance as the stellar HR-diagram (Benz et al. 2014). Figure 2 shows the masses of all companions in the baseline population DIPSY-0 against their semi-major axis. The semi-major axis is calculated in Jacobi coordinates. The final mass of the primary (the host star) is indicated by the colour code in logarithmically spaced bins. For these 100 000 host stars, we find that 7400 companions survive. They cover a large region in the mass-semi-major axis space, spanning five orders of magnitude in semi-major axis and two and a half orders of magnitude in mass. Figure 2 shows all the surviving companions of the baseline population. The distribution of the underlying primary masses is approximately flat in log M_*.final, so does not follow the IMF. In general, we find that the properties of surviving companions depend on primary mass. Due to the binning, this dependence is not directly visible in Fig. 2. We discuss the host stellar mass dependence in Sect. 5.8.

5.2.1 Properties of the mass-semi-major axis distribution

Overall, the envelope of the points covers a region that has roughly the shape of an inverted and rotated ‘L’ with a remarkable paucity of lower mass (planetary-mass) companions inside of 100–1000 AU. Within the region covered, there are no very strong pile-ups or sharp features. In orbital distance, the distribution is, roughly speaking, uniform in log. The mass has a more complex distribution that depends on orbital distance. We discuss the origin of the ‘L’ shape, which is a key result of this study, in Sect. 9.2.

Although we find that all the companions share the same formation mechanism, it is still interesting to classify them according to their final mass into planetary-mass, BD, and stellar companions, because this can be observationally constrained. For this, we used the same mass boundaries between the planetary, BD, and/or stellar regime as Delorme et al. (2024).

For such a classification, our results suggest that only about 11% of all companions are planets (or planetary-mass companions), i.e. objects with masses < 13 M_♃. The vast majority (75%) of the objects are found to be BDs (13 M_♃ to 73.3 M_♃), in agreement with previous studies (e.g. Kratter et al. 2010; Zhu et al. 2012). Finally, we find that about 14% of the fragments end up as stars. The three classes are indicated through the different background colours in Fig. 2. The least massive companion is a 1M_t giant planet, the most massive one a 225 M_♃ (0.21 M_⊙) M-star. The companion closest to the host star has a semi-major axis of 0.15au. The furthest bound companions would be beyond our cutoff of 10⁴ au. We note that companions at very large separations (beyond ~ 1000 au as indicated by the grey shaded region in Fig. 2) might become unbound by stellar flybys or galactic tides (Nordlander et al. 2017; Portegies Zwart et al. 2018; Bancelin et al. 2019; Portegies Zwart 2021; Portegies Zwart et al. 2021; Gratton et al. 2025). This effect is not easy to quantify in our simulations where each star-disc system forms and evolves in the absence of any dynamical interactions with the surrounding environment. The outcome will clearly depend on the age, surrounding stellar environment, and mass of the system. We encourage future studies to investigate this topic.

5.2.2 The initial mass function of the host stars

In this study, we simulated the formation of companions via DI for a wide range of host star masses ranging from 0.05 to 5 M_⊙. This leads to special requirements when presenting the results, as stars of different masses have different occurrence rates. It has long been known that stars with masses >1M_⊙ are much less abundant than their lower mass counterparts (e.g. Salpeter 1955). The observationally determined IMF of the host stars is used to describe this dependency. To have statistical significance (i.e. a high enough number of systems) also at higher host star masses, we divided the parameter space of final host star mass into 100 logarithmic bins in mass, and constructed our simulations such that each bin contains roughly 1000 systems. When presenting the combined results obtained this way, we implicitly present them as if each stellar mass were equally likely in log, by discussing the population as a whole. For example, the companions shown in Fig. 2 only show a preference for the primary mass to the degree that the fragmentation fraction has such a preference, not because certain primary masses are more or less abundant according to the IMF. The advantage is, however, that we can also see the outcome for massive stars, which are of particular interest but are much less frequent.

In the following, we discuss how the host star IMF influences our results, i.e. when the host stars are weighted according to the IMF. We used the IMF given by Chabrier (2005): $$\frac{\text{d}n}{\hspace{0.17em}\text{d}\log m}=\{\begin{array}{ll}0.093\times \exp \left\{-\frac{1}{2}\frac{{(\log m-\log 0.2)}^2}{{0.55}^2}\right\},\hfill & m\le 1{\text{M}}_{\odot }\hfill \\ 0.041m^{-1.35},\hfill & m\ge 1{\text{M}}_{\odot }\hfill \end{array}$$(1)

In Eq. (1), dn/dlog m ξ ξ(log(m)) denotes the stellar number density in pc⁻³ per logarithmic interval of mass and m is the stellar mass in M_⊙. For unresolved systems, the IMF ξ_sys is slightly different below 1M_⊙: $$\frac{\text{d}{n}_{\text{sys\hspace{0.17em}}}}{\text{d}\log m}=0.076\times \exp \left\{-\frac{1}{2}\frac{{(\log m-\log 0.25)}^2}{{0.55}^2}\right\},m\le 1{\text{M}}_{\odot }$$(2)

Instead of weighting each stellar mass bin equally as in Fig. 2, we can also plot the mass-semi-major axis diagram inferred from our simulations for a region in space with a distribution of stellar masses according to the IMF, by weighting the systems with ξ(log m).

This is done by plotting a random selection of companions for each primary mass, where the fraction to be plotted depends on the stellar mass. This results in only a fraction of companions (≈1000) being shown, with almost all coming from the low mass bins. The result is shown in Fig. 4. We find that the mass- semi-major axis distribution is qualitatively similar to the one discussed in Sect. 5.2. The strong preference for M-stars means the properties of this sub-population (discussed in Sect. 5.8) are dominant

The different weights need to be kept in mind when discussing statistical properties of the simulated population as a whole, except for properties with a weak primary mass dependence. Using equal weights for each mass range results in averaged quantities that are not representative of observed distributions of stars. On the other hand, quantities weighted as ξ(log m) are completely dominated by the systems of the lowest masses. Therefore, in the following sections (Sect. 5.3 to Sect. 5.7) we individually analyse narrower host star mass bins for the fundamental companion properties such as the mass or orbital distance distributions.

Fig. 2 Mass vs. semi-major axis for all the surviving companions of the baseline population. The colour code of each companion gives the (final) mass of the primary (host star). Differences in occurrence rate due to the different weights of stellar masses according to the IMF are not reflected in this figure; instead, for each of the 100 stellar mass bins, ≈1000 systems containing a total of 7400 companions are shown. The number of systems (Nsyst), the number of surviving companions (Ncomp), the number of objects accreted on the primary (Nacc), and the number of ejected objects (Nejec) are given in the bottom left corner. The horizontal background colours indicate how the companions might be classified as stars, BDs, or planets according to their mass. Clearly, the distinction between these classes is not as clear-cut as these limits may suggest. Objects at very large separations (⪆ 1000 au) may not survive long enough to be observed, as the grey area suggests, as they can become unbound.

5.3 The companion mass function

The top left panel of Fig. 5 shows the CMF for the baseline population. The distribution of the masses is a central quantity that carries important clues about the formation (e.g. Butler et al. 2006; Nielsen & Close 2010; Ma & Ge 2014). The CMF is given separately for M/K-, G/F- and A/B-stars, respectively. Within the three bins, contributions from stellar hosts of different masses are weighted according to the host star IMF ξ(log m) as explained in Sect. 5.2.2. The shaded area indicates the type of companion according to mass, equivalent to Fig. 2.

The CMF peaks at ≈30 M_♃ to 70 M_♃ in the massive BD range. The distribution is not symmetric, with a sharp dropoff towards the higher masses (above about 100–200 Jovian masses), and approximately a power law towards the smaller masses. In the range from 2 M_♃ to 73.3 M_♃ we find for M/K-stars: N_M/K ∝ M_C^1.1 and for more massive stars: N_G/F/A/B ∝ M_C^1.2.

The inferred CMF is found to be very similar across spectral types, indicating that dynamical processes (collisions and ejections) set it. The only notable difference is at the high-mass end, where M- and K-stars have fewer companions.

Also shown in the figure is ξ_sys, the IMF for unresolved systems (see Sect. 5.2.2). ξ_sys is only constrained down to ≈90M_♃. The extrapolation below this value is shown by the dashed line. At lower masses, the CMF seems to follow ξ_sys approximately. This result suggests that if DI is responsible for the formation of a fraction of the observed low-mass stars, then the BDs and planetary mass objects are just the low-mass tail of this star formation mechanism. We note that our CMF may underestimate the number of stars formed by DI due to limitations in our model, which assumes that the companions formed always remain of lower mass than the primary (see Sect. 6.1 and 9.5).

The CMF shows many objects in the brown dwarf desert (around 30 M_♃ e.g. Grether & Lineweaver 2006). However, the brown dwarf desert is observed at short orbital periods, while most companions in our population are on orbits >> 10 au. Our occurrence rates agree well with observations of distant companions (Sect. 8).

5.4 Distribution of semi-major axes

Similarly to the CMF, we can also investigate how the surviving companions are distributed in distance from the primary, which is also a fundamental quantity, as different formation mechanisms populate different regions (e.g. Fernandes et al. 2019; Fulton et al. 2021). The top right panel of Fig 5 displays the number of companions binned in semi-major axis, again in groups of spectral types.

For G-/F-stars and A-/B-stars, the distribution is similar. Very few companions are found inside of 10–20 AU. Outside, the distribution is approximately uniform in log, with a local maximum at several 100–1000 AU. Companions around M- and K-stars seem to be more evenly distributed in log distances, with additional companions closer to the primary than for more massive host stars. As mentioned, companions outside of roughly 1000 AU might become unbound due to gravitational interactions with other stars.

5.5 Eccentricities and inclinations

Next, we investigated two other important orbital elements (e.g. Bowler et al. 2020; Bryan et al. 2020; Bowler et al. 2023), which are the orbital eccentricities e and inclinations i of the DIPSY-0 population. For multiple systems, a massive inner companion can cause the orbital elements of the outer companions to become meaningless in primary-centric coordinates. Therefore, they are calculated in Jacobi coordinates. The bottom panels of Fig. 5 show histograms of the eccentricities (left) and inclinations (right), respectively, of all companions.

As for the masses and distances, the eccentricities and inclinations shown in Fig. 5 are given as separate histograms for groups of primary spectral types: M- and K-stars, G- and F-stars, as well as A- and B-stars. Within the three bins the eccentricities and inclinations are weighted according to the IMF (ξ(log m)) as described in Sect. 5.2.2, and only systems with a final stellar mass M_*,final > 0.2 are shown. Below this value, very few systems have surviving companions.

We find that the eccentricities are substantial, with typical values of 0.4–0.6. They are found to be very similar across the large stellar mass range, with mean values of 0.41–0.44 and a standard deviation of ≈0.2. We find additional companions with 0.0 < e ≲ 0.2 around M- and K-stars. These high eccentricities are caused by the many scattering events occurring while many objects are in the system.

Inclinations are found to range between 20 and 60°, with mean values rising slightly with stellar mass from 35 to 45°. There is a small number of companions with i up to 180°. This is another indication that the final population is dominated by strong gravitational interactions between the (initially much more numerous) companions. Again, there is little difference across spectral types except for a tendency of more low-inclination orbits in M- and K-stars. This is because in some lower mass systems, only one or a few fragments form, leading to less gravitational interaction. Mean eccentricities and inclinations are given in Table 2.

5.6 Multiplicity of the companions

An important observational constraint for any formation model is the multiplicity of the resulting systems. Here, we only consider systems with at least one surviving companion. The maximum number of surviving companions is three, and Fig. 6 shows the fractions of systems with one, two, and three companions, again separately for M/K-, G/F-, and A/B-stars. As seen in the figure, about 80% of systems have only a single companion, while about 20% have two. Thus, one surviving companion is the most likely outcome. Three surviving companions in a system are very rare (fraction of a percent). Furthermore, the results are again similar for all the spectral types. We discuss the physical processes driving these results further in Sect. 9.1.

5.7 Mass semi-major axis diagram as a function of primary spectral type

In Sect. 5.2, we discussed the mass semi-major axis diagram for the surviving companions of all host star spectral types. Due to the mass dependence of the IMF, it is difficult to compare the inferred populations for different primary spectral types when plotting them all together. Thus, here we individually discuss the resulting populations for the various spectral types. Figure 7 presents the inferred populations for six different spectral types from M to B stars. The spectral types differ in the widths of the corresponding mass bins (Fig. 3). We therefore selected 5000 systems from each spectral type (without scaling according to the IMF inside each type) to be plotted in Fig. 7. The different numbers of companions shown in the different panels reflect the dependence of the fragmentation fraction on the primary mass.

First, one notes that the rotated and inverted ‘L’ of the envelope of points covered in the mass semi-major axis diagram exists for all spectral types in a similar fashion, but that its exact locus and density shifts as a function of host star type. The top left panel shows companions around M-stars (0.08 M_⊙ to 0.45 M_⊙). We only selected systems with M_*,final > 0.2 M_⊙ because very few systems with lower stellar mass fragment. In the selected mass range, the mean number of surviving companions per system, n_c, is 0.04. The number of companions per system accreted on the primary, n_acc, is roughly twice as high, and the number of ejected companions per system, n_ejec is almost 0.6. Among systems with at least one surviving companion, the mean number of companions (i.e. the companion multiplicity), n_mult, is 1.23. Most companions around M-stars are found to be in the BD mass range, some planets can also form, they end up mainly at large distances, but almost none of the companions are in the stellar regime.

For K-stars (0.45 M_⊙ to 0.8 M_⊙), n_c is about four times larger (0.17), resulting in additional points in the top-right panel. The general distribution of the companions in the mass semi-major axis diagram is similar to that for M-stars, with slightly more massive companions at larger radial distances. n_acc is comparable to nc here, n_ejec is larger than for M-stars (≈2).

The behaviour for the populations around G-, F-, A-, and B- stars is also found to be similar. There is some variation in the fragmentation fraction, as seen in n_c, with a corresponding variation in n_acc and n_ejec. The population shifts to more massive, more distant companions as the primary mass increases, leading to more stellar mass companions. We discuss this in Sect. 5.8.

A small fraction of planetary mass companions is inferred around primaries of all spectral types, especially at large separations. We discuss this further in Sect. 9.2. The quantities n_c, n_acc, n_ejec and n_mult forthe various spectral types are summarised in Table 3.

Fig. 3 Initial mass function (Chabrier 2005) for the range of host star masses studied. The objects at the lower mass end are BDs, and the stellar spectral types from M to B are coloured1. Bins in final host star mass used in our study are shown as vertical dashed lines (only every tenth bin is shown for visibility). This demonstrates the range in primary mass covered in our study as well as the logarithmic binning.

Fig. 4 Mass-semi-major axis diagram, weighted according to ξ(log m) (Chabrier-2005 IMF, see text Sect. 5.2.2). With this weighting, the population is completely dominated by M-dwarfs (note the difference in Fig. 2).

5.8 Hoststarmass dependence

So far, we discussed the properties of the baseline population as a whole or in rather large classes according to the spectral type. Here, we systematically discuss the direct dependency of the inferred properties on the final stellar mass M_*,final and its physical reasons. A larger M_*,final means a higher infall rate and/or a longer infall phase. This typically leads to the formation of more fragments. The larger mass reservoir also leads to more gas accretion. Furthermore, discs around more massive stars tend to be larger (S21), and fragments formed at larger distances are more massive (Eq. (35) in Paper I). We can therefore expect that the inferred populations depend on M_*,final, shifting to more massive and more distant companions with increasing mass. Below, we discuss this dependency with respect to the masses, semimajor axes, eccentricities, and multiplicities of the companions. But first, we investigate how the prerequisites for companion formation (disc gravito-turbulence and fragmentation depend on host star mass).

5.8.1 Gravito-turbulence, fragmentation, and survival

Gravitational instabilities can play a role in systems even if they do not fragment. In fact, almost all discs/systems in our baseline population have at least a short phase of gravito-turbulence (Q_Toomre < 2) early in their evolution. Only around the least massive primaries (less than 0.1 M_⊙), 10% to 20% of systems never reach this state.

The fraction of fragmenting discs is much lower. Spiral arms driven by gravito-turbulence are very efficient in transporting angular momentum and mass, stabilising the discs against fragmentation (e.g. Lodato & Rice 2004; Forgan et al. 2011). The least massive systems do not fragment at all, only systems with M_*,final ≳ 0.15M_⊙ will fragment. The fragmentation fraction rises with primary mass up to 0.6 M_⊙ and then decreases again. The increase is caused by the increasing disc mass around more massive primaries. The decrease towards even more massive primaries is due to a combination of several effects. An important factor is the early size of the disc. However, it is unclear whether the observed sample of young discs we used to constrain the early disc sizes is indeed representative for stars with a final mass ≳ 1 M_⊙. If massive stars are absent from the sample, we may under-predict the fragmentation fraction at high masses. The ‘survival fraction’, i.e. the fraction of systems with at least one surviving companion, is again smaller, slightly more than half of the fragmentation fraction or between 5% to 15% with a weak dependency on M_*,final.

In summary, above afinal stellar mass of 0.3 M_≳, all the systems undergo a phase of GI, where 10% to 20% fragment, and about half as many have at least one surviving companion. In Fig. 8 we visualise these results, showing these three quantities as fractions of the total number of systems in bins of spectral type.

5.8.2 Companion mass

The top left panel of Fig. 9 displays the mean final companion mass as a function of final stellar mass. The spread (1σ) in the fragment mass distribution is shown as shaded region. The mean companion mass is increasing with primary mass. Even around 0.2 M_⊙ hosts, the mean companion mass is in the BD regime (≈30M_♃). It increases to 50M_♃ for 5M_⊙ hosts. This means that there is a positive correlation of the mean companion mass with the primary mass although the correlation is rather weak (increase by a factor of about 1.7 for a 20-fold increase in the host star mass). The best-fit power law M_c/M_♃ = 40(M_*,final/M_⊙)^0.2 is shown in the top left panel of Fig. 9. Interestingly, the distribution of masses becomes somewhat wider around more massive hosts (≳1 M_⊙). This happens because more massive systems tend to produce more fragments, and fragments formed later grow less in mass since the bulk of the available gas is accreted by their older, more massive counterparts. The effect is also seen in Fig. 2, where many of the lower mass companions are seen to orbit around the most massive stars.

Fig. 5 Distributions of four key properties for the companions of the baseline population DIPSY-0, separately for M/K-, G/F-, and A/B host stars, weighted according to ξ(log m) within the three bins, for systems with M*,final > 0.2M⊙. Top left: CMF. For comparison, ξsys(log m) is shown as a thin black line (dashed part indicates extrapolation). Top right: Distribution of semi-major axes. Bottom left: Distribution of eccentricities. Bottom right: Orbital inclinations.

Fig. 6 Companion multiplicity for different spectral types (for systems with at least one surviving companion).

5.8.3 Semi-major axes

In the top right panel of Fig. 9, the median semi-major axes S_c are shown vs. final stellar mass, with the dashed lines indicating the first and third quartile. The figure shows that, on average, companions are further from the primary with increasing primary mass. The median semi-major axis is ≈150 au for the least massive primaries and increases to ≈1200au for 5 M_⊙ primaries. The top right panel of Fig. 9 also shows the best-fit power law S_c/au = 550(M_*,final/M_⊙)^0.6.

Fig. 7 Mass vs semi-major axis for surviving companions of the baseline DIPSY-0 population for different spectral types of the primary. The spectral type, from M to B, is indicated in the title of each panel. Each panel contains the companions from a random selection of 5000 systems from the corresponding mass range. In the bottom left corner, the following quantities are given as mean numbers for every sub-population: surviving companions per system nc, objects accreted on the primary per system nacc, ejected objects per system nejec, and number of surviving companions per system with at least one surviving companion, ncs (companion multiplicity). The envelope of points has a similar shape for all spectral classes but systematically shifts in position and density.

Fig. 8 Fractions of systems that are gravito-turbulent. These develop spirals (‘GI fraction’) that fragment and have at least one surviving companion. The fractions are given for groups of spectral types according to the final host star mass M*,final. All values of M*,final 0.05 M⊙ to 5 M⊙ are included.

5.8.4 Eccentricities

It is also interesting to investigate the dependency between the inferred eccentricities and the primary mass. The bottom left panel of Fig. 9 depicts the mean companion eccentricities, including spread, after 100Myr. The figure shows that there is only a very weak dependence, with the mean eccentricity growing from ≈0.34 to ≈0.44. This is caused by the increasing gravitational interaction between companions as the number of fragments increases with host mass.

5.8.5 Multiplicities

Since moremass is added to the disc during a longer infall phase in systems with a higher final stellar mass, such systems typically fragment more often. Fragmentation is also possible further out in the disc if the infall lasts longer and the disc can spread outwards while it is still self-gravitating.

In the bottom right panel of Fig 9, the dashed orange line (using the right y-axis) shows the mean initial number of fragments per fragmenting system (i.e. the mean number of clumps formed in total during the disc’s lifetime). This demonstrates that there is an increase from about 10 to 70 fragments with stellar mass increasing from 0.25 to 3 M_⊙.

Given this trend, one wonders if more massive host stars also have additional surviving fragments in the final state. For higher primary masses, more fragments form, and also of higher mass. However, both of these aspects mean that there are also more gravitational interactions. Typically, in a system where strong gravitational interactions are important (like here because massive clumps emerge relatively close to each other), the system in the end only has a few (one or two, rarely three) companions, independently of the initial number, because of dynamical stability requirements (Ford et al. 2001; Matsumura et al. 2010). Interestingly, this mechanism is also observed in Class IV systems of dynamically active giant planets in a population synthesis model for CA (Emsenhuber et al. 2023).

This self-limiting effect is indeed visible in the bottom right panel of Fig. 9 with the thick red line, which is the mean number of surviving companions per fragmenting system (left y-axis). It is essentially constant as a function of host star mass. The value is ≈1.2 around 0.2 M_⊙ stars and ≈1.3 around 1 M_⊙ stars. The slight decrease towards the most massive stars is probably caused by our treatment where companions beyond 10000au are considered as ejected objects. The comparison of the two curves also shows that only avery small fraction ofall fragments survive.

This result on the final companion number (and the generally weak dependencies on host star mass) suggests that the outcome is not mainly driven by external processes and boundary conditions (like the stellar mass), but by physical processes happening inside of the system itself. As discussed further in Sect. 9.1, the insight that dynamical processes (collisions, scatterings, and ejections) are very important for the final outcome of DI is an important result. It is driven by the high number of fragments that form overall and then interact. This is difficult to fully capture for radiation-hydrodynamic simulations as they typically simulate shorter time spans (Boley 2009;B 18;Boss 2021; Xu et al. 2025).

5.9 Collisions and close encounters

The gravitational interaction between clumps is modelled using the mercury N-body integrator (Chambers 1999, Sect. 4.2 in Paper I). We assume that all companions are point masses, as long as their closest approach exceeds the sum of their outer radii. This excludes the possibility of tidal interactions in close encounters. In addition, the point-mass assumption is clearly an over-simplification as clumps are very extended, have a lot of mass far from their centre, and are likely not spherically symmetric. Therefore, in reality, close encounters could strip some of the clumps’ mass, or even unbind them (see Sect. 9.4 for further details). Given that the radius of each companion is known at every point in time, mercury’s collision handling algorithm can be used to detect when two objects collide. We found that this happens very often (Sect. 5.1). In the ≈12000 systems that fragment with at least two fragments formed, ≈260000 collisions are registered. Among these, ≈231 000 (90%) happen between pairs ofclumps, i.e. beforeeitherofthetwo has undergone a second collapse. This is expected, given that the radii of the clumps before the second collapse are very large (~au), so that their size is not negligible compared to the size of their orbit.

All the collisions are assumed to lead to perfect merging. If two companions collide, they are combined into one object, which has the sum of the masses and momenta (inelastic collision). This is a strong assumption that might not always be appropriate (see further discussion in Sect. 9.4).

Fig. 9 Quantities as a function of M*,final. Top left: Mean final companion mass, Mc. The thin black line shows the best-fit power law. Top right: Median final semi-major axis, Sc. The thin black line shows the best-fitpower law. Bottom left: Mean final eccentricity. Bottom right: Mean number of surviving fragments (left y-axis) per fragmenting system. The mean number of fragments formed initially (i.e. emerging in total during the disc’s lifetime) is shown in orange (right y-axis), also per fragmenting system.

5.9.1 Impact velocities

We found that the majority of collisions occur at low velocity. The left panel of Fig. 10 shows a histogram of the quantity v_c/v_e„ the ratio of relative velocity at impact to mutual escape velocity, for all clump-clump collisions. It reveals that most collisions (70%) actually happen at or slightly below escape velocity. Only 4% have v_c/v_e > 1.5. This happens due to the large sizes of the clumps and the strong damping that occurs on timescales comparable to or shorter than orbital timescales.

5.9.2 Impact angles

Another important quantity is the impact angle, as it can greatly influence the outcome of a collision. The right panel of Fig. 10 shows a cumulative histogram of impact angles for all clumpclump collisions. As seen in the figure, more than 75% of the collisions occur at impact angles less than 45^°. This makes the assumption of perfect merging more reasonable. Also for collisions where at least one body has undergone a second collapse, the perfect merging assumption seems to be reasonable (e.g. Leinhardt & Stewart 2012; Chau et al. 2021). In the figure, the expected distribution of impact angles for a system of objects with random orientations (‘3D’) and for a flat system (‘2D’) is also displayed (Shoemaker 1962). The clumps are initially close to the disc’s midplane (though with small random inclinations). We could therefore expect the distribution of impact angles to lie between the ‘2D’ and ‘3D’ curves. We note that the statistics basedonvc/ve and θ is usually used for compact objects on Keplerian orbits. Both assumptions are not necessarily appropriate for clumps, and therefore, the results may not agree with what is expected from collisions in planetary systems. We discuss collisions further in Sect. 9.4.

6 Variant populations

Several of the underlying physical processes occurring during DI have not yet been studied in great detail with dedicated models. Thus, to calculate the baseline population, we had to make several uncertain assumptions that could significantly affect our results. In the following, we investigate the effect of some of these assumptions using variant populations. We simulated a total of six additional variant populations as summarised in Table 4.

Each of the populations listed in Table 4 consists of 10000 systems, constructed as the baseline population, except for the specific changes we are going to discuss. The mass-semi-major axis diagrams are shown in Fig. 11.

Fig. 10 Left: Distribution of the ratio, vc/ve, for all clump-clump (non-collapsed objects only) collisions (see text Sect. 5.9.1) Right: Distribution of the impact angle, θ, for all clump–clump collisions. The solid orange line shows only collisions between clumps, while the thin dotted black lines shows those between collapsed objects. The expected distribution for random orientations (dotted red line) as well as for a flat system (dashed green line) are also shown. Further explanations are given in Sect. 5.9.2.

6.1 DIPSY-1: No limit on gas accretion

For the baseline population, we limited the rate of gas accretion onto companions to 10^–3M_♃yr^–1. To study the effect of higher gas accretion rates, we removed the limit in population DIPSY-1. The results of this calculation need to be interpreted with caution. Some of the simulations are not physical due to the underlying limitations of our model. In assuming that the disc is always primary-centric, we implicitly assume that companion masses are much smaller than the primary. With unlimited gas accretion, this is no longer valid. In some cases, companions can become even more massive than the primary. While the N-body interactions of said companion with the primary as well as other companions would still be correct, results related to migration, gas accretion, and the disc’s gravitational potential will be unphysical. The top left panel of Fig. 11 shows the mass semi-major axis diagram for population DIPSY-1, with the baseline population given as shaded dots in the background for comparison.

Without limiting the gas accretion rate, we see the emergence of a sub-population of stars with masses of 0.5 M_⊙ to 4M_⊙ at distances of 10 to 100 au that is absent in any of the other calculations. Such stellar companions begin to dominate accretion when reaching a certain mass, and accretion on the primary is cut off. However, as mentioned, our model cannot be used to realistically study such systems. This experimental population is nevertheless helpful. It indicates that, depending on the precise nature of gas accretion, DI might form systems of binary or multiple stars of comparable masses. It is recalled that for stars, multiple systems are very common, with a ratio of single:double:triple:quadruple of 1.50:1.0:0.105:0.026 (Duquennoy & Mayor 1991).

Interestingly, population DIPSY-1 also shows a fairly high number of planetary mass companions. It seems thus possible that this fraction would be higher in systems characterised by two or more stars of similar mass.

6.2 DIPSY-2: No gas accretion and higher initial fragment mass

In this section, we explore the other limiting case of no gas accretion whatsoever. Gas accretion was not included in Forgan et al. (2018) and Nayakshin (2015a). Additionally, we applied a higher initial mass for the fragments. We used the local Jeans mass from Forgan & Rice (2011): $$M_{\text{J},\text{FR}}=\frac{4}{3}\frac{2^{1/4}{\pi }^{11/4}{c}_{\text{s}}^3{H}^{1/2}}{G^{3/2}{\Sigma }^{1/2}\sqrt{1+{\beta }^{-1}}}$$(3)

In Eq. (3), β = t_coolΩ is the local cooling time in units of the orbital timescale (Eq. (27) in S21). M_J,FR is calculated in Forgan & Rice (2011) as a measure for the initial fragment mass in a self-gravitating disc and used as fragment mass in Forgan et al. (2018), see also Sect. 2.9.2 in S21. In a Keplerian disc, it holds: M_F ≈ 0.04 M_J,FR. We use the higher initial mass here for comparability with Forgan et al. (2018), but also for computational reasons: no gas accretion means much more gas is available for later fragmentation and the number of fragments is increased a lot. The top right panel of Fig. 11 shows mass versus semi-major axis of this population.

The similarity of populations DIPSY-2 and DIPSY-0 is remarkable. While there is a stronger concentration at about 100 au, the population without gas accretion overall occupies a similar mass semi-major axis space as the baseline population. This is also seen in the CMF. Figure 12 shows the CMF for population DIPSY-2 for G- and F-stars, the CMF from the baseline population is also shown for comparison. The mass functions of the two populations are similar except for lower mass planets, which would not form due to the higher initial fragment mass. The 2 M_♃ companion at 140 au is a rare case of a clump that survived tidal mass loss. It appears that the higher initial fragment mass compensates for the mass gain through gas accretion. The actual effect is more complex, though, since companions also grow through collisions in our work.

Fig. 11 Mass vs semi-major axis (like Fig. 2) for the variant populations DIPSY-1 to 6 as indicated in the panel titles. The baseline population DIPSY-0 is shown in grey in the background. Please note that the top left panel shows the result of an experimental calculation in which some systems are unphysical because the secondary mass is large compared to the primary mass – a situation not correctly captured by our model. Apart from this aspect, the variant populations are remarkably similar to the baseline one, indicating that internal processes dominate the outcome.

Fig. 12 CMF for population DIPSY-2 (similar to the top left panel of Fig. 5), with no gas accretion and a higher initial fragment mass. The thin black line shows the baseline-CMF for G- and F-stars for comparison.

6.3 DIPSY-3: No initial eccentricity

In population DIPSY-3, we studied the impact of assuming that clumps are formed on circular orbits instead of having a random eccentricity from 0 to 0.3 at birth, as in the baseline population. The middle left panel of Fig. 11 shows this population, again with the baseline population shaded in the back.

The figure does not reveal any clear difference from the baseline population. The three - 1M companions inside of 1000 au may indicate a tendency for additional closer-in low mass companions. Compared to the baseline population, the fraction of surviving companions is slightly smaller, in favour of the fraction accreted onto the host star. This may be because the migration rate tends to be higher on circular orbits. We checked that the distribution of final eccentricities is also very similar to the baseline population, meaning that they are dominated by dynamical interactions rather than by the initial conditions. This is again a sign that the outcome is driven by internal processes among a highly dynamically and collisionally evolved population.

6.4 DIPSY-4: Larger initial disc sizes

The next variant population, DIPSY-4, has an infall location of the gas added onto the discs that is a factor of two larger compared to the baseline population. This leads to much more fragmentation, because more mass is available in the gravitationally unstable region of the disc. The middle right panel of Fig. 11 demonstrates this, showing the more than 3700 companions forming around 10000 stars. For comparison, there were 7400 companions forming around 100 000 stars in the baseline population.

The mean number of surviving companions per system is thus five times higher than in the baseline population, demonstrating once again the importance of the disc size for fragmentation (S23). On the other hand, the distribution in mass semi-major axis space is very similar to DIPSY-0.

6.5 DIPSY-5: Only one clump per fragmentation event

In population DIPSY-5, we only allowed the formation of one fragment in each fragmentation event instead of a random number from one to three in DIPSY-0. The envelope of points in the mass-distance diagram, shown in the bottom left panel of Fig. 11, is quite similar, but in terms of the number of fragments, only roughly half as many fragments per system survive. This outcome is not fully obvious: while only half as many fragments per fragmentation event form on average, the mass not removed will be available for future fragmentation. The reason why fewer fragments survive is the reduced N-body interaction between them. Having several fragments form at the same time makes their dynamical interaction more pronounced, and this tends to prevent some of these from undergoing rapid inward migration. This again demonstrates the importance of N-body interactions.

6.6 DIPSY-6: Lower viscosity and stronger disc photoevaporation

The next variant population, DIPSY-6, is used to probe the influence of the background (minimal) α-viscosity assumed in the simulations. The baseline population uses a rather high value of the background α_BG = 10^–2. It is recalled that the actual (total) α is calculated based on the stability of the disc (see Sect. 2.2 of Paper I) and can be higher in the early phases. At later times, this value might, however, be too high, and if accretion in protoplane- tary discs is mainly driven by magnetohydrodynamic disc winds (see e.g. Bai & Stone 2013; Weder et al. 2023), the remaining ‘background viscosity’ would be lower. We test this by imposing α_BG = 10^–3 in population DIPSY-6. Additionally, we increased the rate of external photoevaporation by a factor of ≈200 (see Sect. 2.7 of Paper I) to 5 × 10^–6 g cm^–1 yr^–1 in order to reproduce the disc lifetimes seen in the baseline population. The resulting population is shown in the bottom right panel of Fig. 11.

The resulting population is again similar to DIPSY-0, in terms of the mass semi-major axis distribution as well as in the number of surviving, accreted, and ejected companions per system. This similarity may seem surprising, as such a different viscosity could influence the disc evolution significantly. However, the early evolution of the fragmenting population is completely dominated by self-gravity. In this phase, the viscosity is given by Eq. (8) in Paper I. The background viscosity is only important for the longer-term evolution of the discs.

In summary, we find that none of the assumptions we tested in this section has a major impact on the resulting population. The exception is the limit on gas accretion discussed in Sect. 6.1, which potentially restricts our study to the formation of systems with a primary that dominates in mass. Of course, there could be other assumptions that would have a larger impact (Sect. 9), and we hope to address this in future studies.

7 Comparison to previous studies

As mentioned in Sect. 1, various population synthesis projects in the DI paradigm have been presented. Below, we compare our results with previous studies.

7.1 Forgan, Rice, and collaborators

Forgan & Rice (2013) performed a population synthesis of selfgravitating disc fragmentation and tidal downsizing (Boley et al. 2010; Nayakshin 2010). Their model used pre-evolved discs around 1 M_⊙ stars and clumps are inserted at random locations in the outer disc. This first study did not include N-body interactions, but this was added later in Forgan et al. (2018). Their Fig. 2 shows a comparison of the inferred population with and without N-body interactions in terms of the distribution of final mass versus semi-major axis. They found that the majority of the surviving objects have masses between 1 M_♃ to 200 M_♃. The masses can at most be as large as the initial fragment mass M_J,FR, since no gas accretion or merges were considered.

Without N-body interactions, most surviving companions have semi-major axes of 10 au to 100 au, close to the birth locations but with possible inward migration. The 2018-update demonstrates the importance of gravitational interaction, with an inferred population ranging from 10^–1 au to 10⁵ au. The bulk of the population lies in a similar mass-semi-major axis space as our baseline population, i.e. their population is also dominated by massive giant planets and BDs at large semi-major axes. There is even a hint of a similar rotated ‘L’- of the envelope of points in their mass semi-major axis diagram that is characteristic for our results, i.e. a presence of lower/planetary mass objects at distances outside of about 100 AU and an absence of such objects between 1 and 100 AU.

The distribution of multiplicities they inferred (right panel in their Fig. 4) is similar to our result (Sect. 5.6). They find a 5:1 ratio of single:double companions, while we obtain approximately a 4:1 ratio.

However, there are also significant differences: their simulations predicted the presence of a significant number of objects inside of 1 AU covering a large mass range, which is absent in our population. These companions are formed by inward scattering. Their mass function is bimodal with peaks at about 25 and 80 Jovian masses, whereas we obtain an unimodal distribution. The eccentricity distribution is also different: while ours peaks at 0.4 with a fall-off towards both lower and higher values, their distribution peaks at about 0.05 and is flat between 0.1 and 0.95.

7.2 Nayakshin models

Nayakshin & Fletcher (2015) performed a population synthesis study of planet formation around solar mass stars in the DI paradigm including tidal downsizing. Their surviving population mainly consists of two groups. A group of rocky planets with masses between approximately 0.1 M_⊕ to 10 M_⊕ and separations from 0.1 au to 10 au. This group results from inward migration and tidal downsizing stripping the gaseous envelope of a solid core that forms via the accretion of solids by the gaseous clumps. Solid accretion is not included in our model, meaning that such a scenario cannot happen by construction (tidal downsizing in the sense of loss of gas only is included in our model, but not the formation of solid cores). The second group consists of giant planets (0.5 M_♃ to 10 M_♃) that did not suffer such catastrophic mass loss of the gaseous envelope. The separations of these objects lie predominantly between 10 au to 90 au, with a fraction residing closer to the star (down to 0.1 au). This is a region of parameter space that is almost not populated in our study. Their masses are also clearly lower (about 0.4–10 Jovian masses, with no objects more massive than 20 M_♃), in strong contrast to our findings. This is mainly due to the absence of gas accretion and mergers. Initial fragment masses are typically lower than in Forgan et al. (2018). Furthermore, while our companions often have large semi-major axes (typically at least tens of AUs or even hundreds of AUs), no objects are found in their simulations outside of 90 AU. Thus, there are significant differences between the predictions made by these two models. The main reason for these differences is clearly that Nayakshin & Fletcher (2015) consider one fragment per system and consequently no N-body interactions.

8 Comparison to observations

The majority of companions formed in our population are on wide orbits, with semi-major axes between 10¹ au to 10⁴ au. Direct imaging surveys are therefore best suited to compare our results against. A quantitative statistical comparison of our populations to observational surveys is not straightforward due to observational biases, the underlying assumptions such as the mass-luminosity relation, different sample properties, etc. We thus perform an in-depth comparison with observations in a dedicated publication. In the following, we compare our results only with a limited number of observational results, concentrating on the frequency of companions, a question which has been at the centre of attention in past surveys. In Sect. 8.1 we discuss the conclusions from a review of older surveys. Sections 8.2 and 8.3 look at results from two specific, more recent surveys.

8.1 Review by Bowler & Nielsen 2018

A number of large imaging surveys have been conducted in the past. One of the first is Biller et al. (2013). It is possible to perform meta-analyses of such surveys to constrain the occurrence rate of wide companions. Bowler & Nielsen (2018) review the analyses available at the date of writing. The authors conclude that the occurrence rate of massive giant planets (from 5 M_♃ to 13 M_♃) at separations between 5au to 500 au lies close to 1%. Our baseline population predicts a lower number in comparison. There is a population of planetary mass objects. However, as is obvious from Fig. 2, most of these companions lie beyond 500 au. For the limits in mass and semi-major axis given above, we find companion fractions of 0.4‰ to 0.6‰ around all primary spectral types, so approximately 20 times less. Bowler & Nielsen (2018) also give an occurrence rate between 1% to 4% of BDs (with masses from 13 M_♃ to 75 M_♃) at the same separations. Here, the agreement with our baseline population is very good: we find 2% to 3%. This is a non-trivial result given that the model was not optimised to reproduce any observational constraints.

The only population we synthesised that has more planetary mass companions inside 500 au is population DIPSY-4 with a larger early disc size and thus more fragmentation. In this population, ≈0.5% of systems have a companion with a mass of 5 M_♃ to 13 M_♃ between 5au to 500 au. However, the fraction of systems with a BD companion is also larger ≈(10%) than what is found observationally by Bowler & Nielsen (2018), showing that in our model, we cannot simply adjust one quantity without (undesired) effects for others. An overview of these numbers is given in Table 5.

8.2 SHINE project

The SHINE (SPHERE Infrared Survey For Exoplanets) project is a large survey that uses the SPHERE instrument at the Very Large Telescope (Desidera et al. 2021; Langlois et al. 2021). Vigan et al. (2021) perform a statistical analysis of a subsample of 150 stars with spectral types from B to M. The authors use the BEX-COND-hot evolutionary tracks (Marleau et al. 2019) to infer the mass of the observed companions. Using a parametric model, they predict the frequencies of systems with at least one companion with a mass between 1M to 75 M and a separation from 5 au to 300 au. Their parametric model also attributes a frequency to a ‘planet-like’ (i.e. CA) formation pathway called ‘PPL/LN’ and a ‘binary star-like’ formation pathway called ‘BDB’, which might also be fuelled by DI. The results from Vigan et al. (2021), together with the results from our model in the same mass/semi-major axis limits, are summarised in Table 6. The result from the SHINE survey predicts companion frequencies from binary-like formation at the percent level. This agrees qualitatively with the prediction from DIPSY. There is quantitative agreement within a factor of 2.5, which is a fair agreement given the many rough model approximations involved, such as the low dimensionality of the model. The trend with stellar mass is at this point unclear. Vigan et al. (2021) state that their results seem to show a local minimum in the frequency of substellar companions around FGK stars. They caution, however, that this result should not be overinterpreted because their sample only contains 20 M-stars, and that the analysis of the full SHINE sample should yield more robust results. Also, within their 1-sigma error bars, the observed BDB- frequencies are compatible with a BDB-frequency independent of spectral type.

At the present time, the preliminary SHINE results are compatible with an increase in companion frequency from K/G/F to A/B stars, while we find a slight opposite trend (decrease by a factor of 0.68). A possible explanation lies in our initial conditions. The sizes of young discs are quite uncertain, especially as a function of stellar mass. In our simulations, the infall location is chosen such that (S23) the observed early disc sizes as found by Tobin et al. (2020) are reproduced. It is not known if the observed sample contains A/B stars.

8.3 BEAST survey

The BEAST (B-star Exoplanet Abundance Study, Janson et al. 2021b) survey is a direct imaging survey also using the SPHERE instrument, but targeting B-stars and focusing on the detection of exoplanets. Delorme et al. (2024) study the frequency of giant exoplanets by statistically analysing the first half of the BEAST survey. The sample consists of 34 young B-stars with a median mass of 4.8 M_⊙. In the selection of this sample, stellar companions between ≈12au to 720 au are excluded. When comparing to our population, we therefore also exclude such systems. The occurrence rate of substellar companions around B-stars is given as ${11}_{-5}^{+7\%}$ by Delorme et al. (2024). Using the aforementioned constraints, the corresponding number from DIPSY-0 is ≈7%. Again, given the many underlying simplifications, this seems a fair agreement. However, while ≈60% of the companions from DIPSY are massive BDs (30 M_♃ to 73.3 M_♃), the two detected companions that are part of the BEAST survey are both at most low-mass BDs: b Centauri b (Janson et al. 2021a) has a mass of (11 ± 2) M_♃, µ² Sco B (Squicciarini et al. 2022) a mass of (14 ± 1) M_♃. In terms of the companion-to-star mass ratio, the discrepancy becomes even bigger. It therefore appears that b Centauri b and µ² Sco B are not typical outcomes of DI as represented in our current model². This will be studied further in future work. Future comparisons will also include other quantities, such as the orbital distance, multiplicity, or the eccentricity distributions (Wahhaj et al. 2013; Bowler et al. 2020).

9 Discussion and limitations

We performed a population synthesis with a model that includes a significant number of self-consistently linked physical effects compared to some previous studies. While doing so, we still had to make a number of assumptions and simplifications. In this section, we discuss some of the physical processes driving our results but also how the most important assumptions and simplifications could affect our findings.

9.1 High number of clumps formed and consequences: Scatterings and collisions as internal processes

Our simulations agree with previous radiation-hydrodynamic studies in the finding that mass loading onto the protoplanetary disc from infall can drive the disc into multiple episodes of fragmentation (Boley 2009; Zhu et al. 2012; Oliva & Kuiper 2020). In such hydrodynamic simulations, which simulate (in comparison to our infall duration) often only shorter time intervals (several 10³ to a few 10⁴ years), the formation of a significant number of fragments can be observed especially if the mass loading is high and the infall radius is large. For example, in one simulation of Boley (2009), 11 fragments form, while in the simulation of Zhu et al. (2012) with the highest mass loading rate and infall radius, half a dozen fragments emerge rapidly. In our simulations, the durations of the infall phase are chosen to lead to a stellar mass distribution following the IMF. While typical infall durations are a few 10⁴ years, there is also a non- negligible fraction of longer durations of 10⁵ years (S23). These longer timescales are especially prevalent for more massive stars, because in our setup, in order to form more massive stars, the infall duration needs to be longer on average.

This explains the (at first sight) non-trivial result, that in the baseline population, among 100 000 stars, a huge number of about 430 000 fragments form in total. Given that only a minority of systems fragment (Sect. 5.8.1), fragmenting systems will have tens of fragments formed over their lifetime on average. In Fig. 9, the mean number of fragments is shown as a function of the final host star mass. One can see a clear increase from about 10 at 0.2 M_⊙ to 70 at 30 at 3 M_⊙. This reflects the longer infall duration, allowing more fragments to form in total.

This very high number of initially (i.e. in total) formed clumps needs to be sharply distinguished from the final number of companions per fragmenting disc, which is only about 1.2 with nearly no dependency on host star mass. This implies that the formation process is very inefficient at turning initial fragments into bound companions. Instead, almost all fragments are either destroyed in collisions, or (with a factor of ~2 lower probability) ejected. With these results at hand, a new vision of DI arises where the final outcome is a strongly collisionally and scattering-dominated process. These ‘internal’ processes explain key outcomes such as the number of final companions or the eccentricities, and the weak dependency of many results on host star mass.

Fig. 13 Example of the formation of a system with an inner more massive and an outer, less massive companion. The lines show the growth tracks of the clumps. Arrows to the right denote ejections, and the ‘T’ stands for the (ultimate) thermal destruction.

9.2 Origin of the ‘L’ shape in the mass-distance plane

One of the key findings of this study is the prediction of the rotated and inverted ‘L’ envelope that is populated by the companions in the mass semi-major axis diagram with a paucity of planetary mass objects inside of 10 au to 100 au. Figure 13 shows the formation tracks in an example system with two surviving companions of about 40 M_♃ at 100 AU and a 13 M_♃ at 2300 AU. The final primary mass is 1.06 M_⊙. Initially, 23 clumps emerge within the first 30 kyr between 50 and 200 AU with masses of 1 to 2 Jovian masses. They then grow by rapid gas accretion and collisions. Typically at 0M_♃, the second collapse occurs. One clump forming inside of the others undergoes inward migration and is destroyed by the combined action of disc irradiation and tidal mass loss. After 100 kyr, the system contains 9 collapsed clumps, which interact gravitationally. By 320 kyr, only the two final companions are left. From the growth track of the outer lower mass companion, one can see that outward scattering out of the gas-rich part of the disc leads to the formation of the lower mass object further out, while the clump scattered outside grows by gas accretion and, initially, collisions. Thus, in order to not grow to masses higher than planetary mass, the orbital distance needs to be larger than the discs’ size of 100–1000 AU. Otherwise, they grow into the brown dwarf or even stellar regime. Thus, the combination of the emergence of many interacting clumps, scatterings, and gas accretion inside and the absence of it outside leads to the ‘L’ shape.

We note that there is an interesting analogy in low-mass star formation: One hypothesis suggested that low-mass stars might form due to ejection from denser regions of a collapsing molecular cloud. While current understanding particularly emphasises that low-mass stars form from the turbulent compression and fragmentation in molecular clouds producing collapsing cores over a wide range of masses, it is still possible that ejection also plays a role depending on the formation environment (see e.g. Luhman 2012).

9.3 Infall model

For the formation of the star-and-disc systems, we assumed a relatively short infall phase characterised by a high and constant infall rate at the beginning of the simulations (Sect. 4, Sect. 2.5 of Paper I). While short episodes with high infall rates are also seen in hydrodynamic simulations (e.g. Bate 2018), a number of observations show that accretion varies strongly during stellar mass assembly (Fischer et al. 2023). Such accretion variability could modify our results. For example, if the mass supply to the disc is dominated by episodes of high accretion rates that are separated in time, this could lead to fewer fragments being in the disc at the same time, which could reduce the number of surviving fragments. At the same time, sudden drops in accretion rate will favour fragmentation since they drive the disc towards instability through the decreasing temperature. The net effect of accretion variability on the outcome of DI is difficult to predict, since it depends on the precise nature of the variability, which also has an influence on the host’s properties (Kuffmeier et al. 2018). We plan to investigate this topic further in the future.

9.4 Collisions

Collisions occur very often in our simulations and they play a key role in setting the final companion masses. Clump-clump collisions are particularly frequent due to the very large sizes of the bodies (typically -au, which is orders of magnitude larger than ‘normal’ (collapsed) planets and BDs). In assuming perfect merging, we made a strong assumption for the outcome of collisions. This base assumption is to some degree justified, as we observed low impact velocities around or even below the mutual escape velocity because of the very strong damping (see Sect. 5.9). This is a very particular collisional regime, different from what is typically seen in collisions between, for example, rocky planets. Interestingly, a similar regime exists for mergers of asteroids that are also large (with respect to their semi-major axes) and damped by their debris disc (Wimarsson et al. 2025). A very recent study (Matzkevich et al. 2024) focusing specifically on collisions between clumps finds a number of different possible outcomes. Clump-clump collisions can, for example, result in disruption, erosion, or hit-and-run, depending on their ages, masses, and impact angles. In some cases, the collision can even trigger the dynamical collapse of a clump. This outcome may be especially relevant for our study, since the collapse makes the clumps robust against tidal disruptions. However, Matzkevich et al. (2024) find this phenomenon is most likely for evolved clumps, some of which would soon undergo dynamical collapse even without a collision. Thus, we do not expect this outcome to change our results very much. Furthermore, the collision outcomes for clump-clump collisions from Matzkevich et al. (2024) (their Figs. 9–11) suggest that in some cases perfect merging is a reasonable assumption for collisions below mutual escape velocity at intermediate impact parameter. Such collisions are found to be rather common in this work. Clearly, further investigations of this topic are required. It is difficult to judge in what way a more elaborate collision model would influence our results. Less perfect mergers may decrease the clumps’ growth, possibly keeping them at a lower final mass and making them more prone to tidal disruption. However, regardless of the outcome of collision, the total mass involved needs to be conserved. So if a lot of mass is returned to the disc after an erosive or disruptive collision, this gas would again be available for accretion or further fragmentation.

Another simplification related to the collision model are close encounters. As long as the centres of two bodies do not approach closer than the sum of their radii given by the interior structure model, they interact gravitationally as point masses. For clumps, this may not be a good approximation, since a large fraction of their mass is relatively far from their centre (Galvagni et al. 2012). and may be stripped/unbound in a close encounter. To our knowledge, no hydrodynamic simulations exist at the moment that specifically study close encounters of clumps (though see Fig. 6 in Xu et al. 2025). It will also be important to address this topic in the future.

9.5 Primary-centric disc and gas accretion

The accretion of disc gas is among the major factors influencing the final population of companions. At the same time, it is a complex, three-dimensional process, which is still not very well understood. Hydrodynamical simulations of gas accretion onto giant planets do not agree very well on the accretion rates. For example, Fletcher et al. (2019) perform a comparison of seven different hydrodynamical codes simulating the migration and gas accretion of clumps at wide separation. They find a factor of several difference for the accretion rate of an initially 2M clump at 120 au with identical initial conditions. The differences are lower in the case of clumps with an initial mass of 12 M, and accretion rates are of the order of 10⁻³ M yr⁻¹. We used this value as an upper limit for the gas accretion rate. This seems reasonable also when comparing to Shabram & Boley (2013). The authors perform radiation-hydrodynamic simulations to study the early evolution of circumplanetary discs (CPD) around companions formed through disc fragmentation. They find that long-term accretion rates from the CPD to the companion reach 3 × 10⁻⁴ M yr⁻¹.

Without the limit on the accretion rate, our accretion model would allow some companions to approach the mass of the primary (Sect. 6.1). While this may be a pathway for the formation of binary- or multiple-star systems, it is problematic in the context of our model. The assumption of a primary-centric disc only makes sense if the primary is much more massive than any of the companions. If much higher accretion rates are possible (as seen e.g. in Zhu et al. (2012) at much higher infall rate or in Oliva & Kuiper (2020) in the case of massive star formation) this would lead to a shift of the barycentre, clearly invalidating the primarycentric disc assumption. Our limit on the accretion rate therefore effectively restricts the probed parameter space to systems with a massive primary.

It is also possible that accretion onto clumps is inhibited due to effects not included in our model. Nayakshin & Cha (2013) study the migration and accretion of clumps and use analytical and numerical arguments to show that radiative feedback should reduce the accretion of disc gas onto low-mass clumps. A similar effect is seen in Stamatellos (2015). Deng et al. (2021) investigate disc fragmentation in the presence of the spiral-driven dynamo (Riols & Latter 2019; Deng et al. 2020) and find that initial fragment masses are much lower than in the absence of the dynamo, and that gas accretion is stifled due to magnetic shielding in the self-gravitating phase. The combination of lower (approximately an order of magnitude) fragment masses and suppressed gas accretion would likely influence our results. However, it is not straightforward to say how, as any mass not bound in clumps would be available for further fragmentation or future accretion.

9.6 Magnetic fiields

Magnetic fields do not only affect fragmentation and the accretion of gas, but also play an important role in the formation of the star-and-disc system. We discuss this in detail in S21. The effect of the magnetised collapse is indirectly taken into account in our model by setting the infall location in such a way that we retrieve early disc sizes in agreement with observations (S23). Clearly, a more physically motivated approach that quantitatively captures the physics of all non-ideal MHD effects and its impact on the infall location as well as the angular momentum would be desirable and opens a perspective for future work.

Additionally, magnetically driven disc winds are thought to play a key role in driving accretion and removing mass in pro- toplanetary discs (Bai & Stone 2013; Turner et al. 2014). MHD winds are also likely to influence gas-disc migration, though this is an active field of research and the details are not yet clear (Ogihara et al. 2015; Lega et al. 2021; Paardekooper et al. 2023; Wu & Chen 2025; Hu et al. 2025). Including the effects of MHD winds in the future would be desirable (see e.g. Weder et al. 2023; Suzuki et al. 2016).

Furthermore, it was shown when properly accounting for the interplay of self-gravity and magnetic fields in young discs (Riols & Latter 2019; Deng et al. 2020), GI may also lead to the formation of intermediate-mass planets (Deng et al. 2021). Accretion and migration in such systems may behave differently in such systems (Kubli et al. 2023; Kubli et al. 2025).

9.7 Inner disc edge

We did not study in detail what happens at the inner edge of the disc. The inner truncation is fixed at 0.05 au, and whenever a companion crosses it, it is considered accreted onto the primary. This assumption is adequate because we are mostly interested in the surviving population at wide orbits. At the same time, our results show that the number of companions accreted on the central star is comparable to the number of surviving objects. If even a fraction of the accreted population was able to survive close to the primary, it could explain the formation of some observed companions, for example relatively close giant planets around M-dwarfs (Morales et al. 2019). A more detailed treatment of the inner edge, including a dependency on the stellar magnetic field and rotation rate would be desirable for this. Additionally, the companions that do accrete on the primary also need to be considered. Such events will change the composition of the star and will lead to outbursts that may be observable (Vorobyov & Basu 2010; Soares et al. 2025).

9.8 Metallicity correlations

Giant planets have long been known to be more frequent around metal-rich stars (Gonzalez 1997; Santos et al. 2001). This giant planet-metallicity correlation has become well established (e.g. Fischer & Valenti 2005, see Adibekyan 2019 for a review). It is considered strong support for CA being responsible for the formation of the majority of giant planets (Ida & Lin 2004; Mordasini et al. 2009b; Wang et al. 2018; Emsenhuber et al. 2021). There is some discussion in the literature about this trend only being valid up to a ‘transition mass’ of a few M_♃. Above this mass, planets and BDs would form through DI (Schlaufman 2018; Santos et al. 2017). Adibekyan (2019) analyses and discusses these claims in some detail based on a large, homogeneous dataset from the SWEET-Cat catalogue (Santos et al. 2013). They find that the correlation of planet occurrence with host metallicity continues above 4M, though it is less strong. They conclude that the claim about the different formation mechanism is not justified. Instead, they argue that the formation environment is different. This finding seems to be compatible with the results of our study. As discussed in Sect. 8.1, our populations seem to produce fewer giant planets between 5 au to 500 au than observed by direct imaging.

In its current state, our model would likely predict an anticorrelation of host star metallicity with companion frequency. This is because in our clump evolution model, clumps of lower metallicity have lower opacities and shorter pre-collapse timescales and are therefore less likely to be tidally or thermally disrupted (see also Helled & Bodenheimer 2010; Humphries et al. 2019).

However, this simplistic argument overlooks important physics. A higher metallicity could also lead to a lower clump opacity through grain growth and settling. This is discussed in Helled & Bodenheimer (2011). These authors model isolated clumps by solving the stellar structure equations. They study the effect of grain growth and settling on the pre-collapse timescale of the clumps and find it to be much shorter than without this effect, and insensitive to clump mass. Additionally, the accretion of pebbles and planetesimals should also be considered (Baehr & Klahr 2019). Nayakshin & Fletcher (2015) include accretion of pebbles and find a positive correlation of giant planet occurrence with metallicity for planets inside of 5 au, and no correlation outside of 10 au. Since our model considers neither accretion of solids nor solid dynamics inside clumps, we cannot make any clear predictions about metallicity correlations at this time. We plan to address this important topic in a subsequent paper of this series.

9.9 Clump geometry and clump rotation

Throughout this work, we assumed that clumps are spherically symmetric objects that do not rotate. While necessary at the current stage of the model, both assumptions are not satisfied for clumps that form in a disc. In a more realistic scenario, a clump forming from a collapsing region of a disc will inherit angular momentum from the Keplerian shear and is likely to have an elongated shape. This is seen clearly in hydrodynamic simulations of disc fragmentation (see e.g. Fig. 2 in Boley et al. 2010, Fig. 5 in Zhu et al. 2012, Fig. 7 in Vorobyov & Elbakyan 2018, or Fig. 6 in Xu et al. 2025). The clump can also gain angular momentum through gas accretion and will spin up as it contracts (Zhu et al. 2012; Galvagni et al. 2012). Clump rotation could be expected to increase the clump’s pre-collapse timescale, as it works against the contraction. But the effect is more complicated, and the development of non-axisymmetric modes can shorten the pre-collapse (Galvagni et al. 2012).

At this time, it is difficult to judge how these effects will influence the inferred population. Future studies should include and address these topics.

10 Summary and conclusions

We performed a large-scale population synthesis study of companion formation in the DI paradigm, applying the model described in Paper I. By combining descriptions of physical processes in star and disc formation and evolution (e.g. infall, disc structure, and stellar evolution), disc GI and fragmentation, multiple clump formation, evolution through gas accretion, interior structure evolution with second collapse, clump orbital migration, clump thermal and tidal destruction, as well as N-body interactions with collisions and many other effects, we were able to globally follow companion formation by DI and quantitatively predicted occurrence rates and physical properties of companions around stars from 0.05 M_⊙ to 5M_⊙. We find DI to be a mechanism that forms companions with masses ranging from the low-mass star down to the planetary regime, although it very rarely produces planetary mass objects inside of about 100 AU (at a per-mill level). The main results from our baseline population DISPY-0 are summarised below:

• For final stellar masses above 0.1–0.2 M_⊙, virtually all discs undergo an early self-gravitating phase where they form spiral arms;

• Fragmentation (i.e. the formation of bound gas clumps) occurs in - 10% to 20% of systems with final stellar masses M_*,final ≳ 0.1 M_⊙,;

• When fragmentation occurs, the total number of fragments formed is large, with tens of fragments (clumps) forming on average per fragmenting system and up to hundreds in extreme cases. This high number is explained by the fact that continuing infall drives the discs through several fragmentation events;

• However, only a small fraction (2%) of the clumps survive as bound companions on a Myr timescale;

• Fragments are lost due to (approximately) collision (60%), ejection (25%), thermal disruption from disc irradiation (10%), accretion on star (2%), and tidal disruption (0.5%);

• Based on their mass, surviving companions are grouped into low-mass stars (15%), BDs (75%), and planetary mass objects (10%);

• The CMF is unimodal and peaks at about 50 Jovian masses. Below this value, the frequency of companions scales with mass approximately as M_c¹. Our model might underestimate the number of stellar companions formed;

• Brown dwarfs and stars occupy a large range of semi-major axes (10¹ au to 10⁴ au), while most planetary mass objects lie beyond 500 au. This means that in the mass-distance plane, the envelope of points covers a characteristic inverted and rotated ‘L’ with a paucity of planetary-mass companions inside of 100–1000 AU;

• In terms of frequencies, DI produces planetary-mass companions on a low percent level at distances larger than ~100 au. Inside of 100 AU, planetary-mass companions are, in contrast, very rare (low per-mill level);

• This key outcome is due to the following physical process based on dynamical clump-clump interactions (scattering) and growth: When clumps remain in the disc, there is typically enough gas for them to grow into the BD or stellar mass regime by both gas accretion (Kratter et al. 2010) and clumpclump collisions. Only ifthey are scattered outwards to large distances beyond the discs’ outer radii, and thus deprived of gas to accrete, can they more commonly remain in the planetary mass regime;

• Compared to results of direct imaging surveys (Bowler & Nielsen 2018), the synthetic occurrence rate of BDs is in good agreement with observations, while the number of planetary mass objects is about a factor of 20 lower, pointing at CA as the dominant formation mechanism (Janson et al. 2012; Kraus & Ireland 2012);

• Surviving companions have large eccentricities (e - 0.4 ± 0.2) and typical inclinations around 40-, both with a large spread;

• Surviving companions around more massive stars are on average more massive and on wider orbits, though the dependencies are weak;

• Among systems with at least one surviving companion, four out of five only have one single companion; one out of five systems has two; and very few (<1%) have three companions, with almost no primary mass dependency;

• Companions are typically −25 × as massive after 100Myr than at fragmentation, due to gas accretion and mergers;

• Clump-clump collisions are very frequent and also contribute considerably to the final masses. We thus find that collisions and gravitational interactions (scatterings) in general are key in shaping the emerging companion population formed by DI. The fact that the outcome is driven by clumpclump interactions as an internal process explains why we observe remarkably weak dependencies of our results on many external (boundary) conditions such as stellar properties but also several important model assumptions;

• Ejections, as a related phenomenon, are very frequent and lead to £1 ejected object per star for all primaries.

Figure 14 provides a graphical summary ofmany ofthese results, with a flow chart splitting up into different outcomes from the top to the bottom. For a sample of 1000 F/G stars, it first displays the number of discs undergoing GI. Given that we only used a sub-range of (particularly relevant) stellar masses in this study, the exact fractions differ to some extent from the ones given in the textual summary. As discussed, virtually all discs are selfgravitating at some time during the infall and early evolutionary phase. For those undergoing such an instability, we can then see how many fragment (about 19%). The majority of systems do not fragment due to the fast transport of mass and angular momentum by the spiral arms (Sect. 5.8.1). Onthenextlevel, the clump level, the total number formed is given, as well as their fundamental fate (ending up as bound companions, ejected, or destroyed). It is again remarkable to see how many initial clumps form (5581, i.e. on average about 29 per fragmenting disc), but how few (128) remain in the end as bound companions. Ejection as the fate is already more than ten times more common (1567 clumps), but destruction (3886 clumps) is clearly the most likely outcome via one of the four listed mechanisms. Among these, collision with another clump is by far the most important channel. Finally, for the bound companions that remain at the end of the simulation, we classify the bodies according to semimajor axis- smaller (76) or larger (52) than 1000 au-and further into planetary mass, BD, and stellar companion based on their mass. We see that inside 1000 AU, only three planetary-mass objects eventually result (about 3%o), compared to 55 BD and 18 stars.

As a global model, our framework depends on a non- negligible number of parameters and base assumptions about the governing physical processes. We thus also quantified a number of variant populations to study the sensitivity on different parameters and assumptions. From these we find:

• Relinquishing the baseline limit on the gas accretion rate of 10⁻³ Jovian masses per year, often seen in hydrodynamic simulations (Shabram & Boley 2013; Fletcher et al. 2019), can lead to the formation of very massive companions, including equal mass binaries. However, these simulations are not quantitatively physically realistic, as they are not consistent with our base assumption that the central object (primary star) always remains more massive than any companion formed;

• Disabling gas accretion altogether and, instead, using a higher initial fragment mass (as in Forgan et al. 2018) leads to a population surprisingly similar to the baseline;

• Whether or not clumps have initial eccentricity does not have a significant influence on their survival in our model, a testament to the strong internal driving of the system by gravitational interactions;

• A larger, early disc size (larger than observed) in contrast to the calibrated sizes in the baseline population) leads to a disproportionately higher number of fragments, as well as surviving companions. The fraction of planetary mass objects is higher in this case, but the fraction of BDs is too high compared to observations;

• The possibility for several fragments to form per fragmentation event increases their chance of survival;

• The background alpha-viscosity has a minor influence on the population as long as the disc lifetimes remain similar. The underlying reason is that in the most important early phase, self-gravity determines the actual (high) alpha.

We conclude that DI is a robust mechanism for the formation of low-mass stars, BDs and planets, although it largely avoids populating the inner 100-1000 au with planetary mass objects.

The number of planets formed this way may be low (by a factor of ~20 lower when compared to some direct imaging results), in agreement with predictions from the literature (Rice et al. 2015; Kratter et al. 2010). However, important physics is still missing in our model (and sometimes not even well understood), so it is too early to firmly draw this conclusion (Sect. 9). The number of BDs produced by our model agrees well with occurrence rates from direct imaging surveys. A more thorough comparison with observations, including data from upcoming larger surveys, will be conducted in future studies to put the wealth of quantitative predictions provided by the synthetic systems to the observational test. Such comparisons will reveal differences between synthesis and observations; however, the differences will precisely allow us to better understand the formation process and make it possible to improve the models.

In agreement with previous studies (Papaloizou & Terquem 2001; Forgan et al. 2015; Forgan et al. 2018), we predict DI to produce a high number of ejected objects.

We identified a number of processes that have a key influence on the outcome of DI. The most important are the following: the formation of a very high number of clumps due to continuous fragmentation in infalling discs; dynamical interactions (scatterings and ejections) among the clumps concurrently present in a disc; gas accretion by the clumps leading to significant mass growth; and a very high number of clump-clump collisions. The characteristics of disc formation (in particular, the early disc size) are also important. The precise nature of these processes should be studied further by dedicated models in the future. The results of such studies can then be included in the next generation of population synthesis models.

Some physical mechanisms, such as solid accretion and MHD-effects, are still partly or fully missing from the model (Nayakshin & Fletcher 2015; Baehr & Klahr 2019; Bai & Stone 2013; Riols & Latter 2019; Deng et al. 2021, 2020; Kubli et al. 2023; Kubli et al. 2025). Including them in the future will be important as they could influence our results in a relevant way.

Gravitational instabilities (not necessarily including fragmentation) appear to be ubiquitous in protoplanetary discs. Their effect should be considered in studies of planet formation beyond the DI framework as well, as it influences the disc’s structure of both the gas and the dust component (Boss 2015; Rowther et al. 2024). Such disc inhomogeneities might play a role in the dynamics of pebbles as well as the formation of planetesimals and of (solid) protoplanets. Disc instability and Ca might thus not be considered fully distinct formation pathways. Hybrid models have been proposed (e.g. Durisen et al. 2007), and they may be necessary to explain some of the recent observations of young protostellar and protoplanetary systems. The quantitative statistical comparison to observations, especially with large surveys such as SHINE (Desidera et al. 2021; Langlois et al. 2021; Vigan et al. 2021), GPIES (Nielsen et al. 2019), BEAST (Janson et al. 2021b; Delorme et al. 2024), or COBREX (Squicciarini et al. 2025), will help untangle the mysteries of star and planet formation.

Fig. 14 Graphical summary of the results of the baseline population DIPSY-0 for F/G host stars in terms of the frequency of different outcomes. On the left, the basic level (star, disc, clump, and companion) is shown in brown. The rest shows a decision-tree-like flow chart with the different outcomes. The colour coding is taken from Figs. 2 and 1 for the mass ranges of the surviving companions and the fates of the destroyed clumps, respectively. A central outcome, the frequency of bound companions inside of 1000 AU, is shown in blue bold face. Please note that the numbers are subject to model assumptions.

Data availability

The companion data for the populations is available for download on Zenodo.

Acknowledgements

We thank Michael Meyer and Lucio Mayer for the insightful discussions. OS and CM acknowledge the support from the Swiss National Science Foundation under grant 200021_204847 ‘PlanetsInTime’. R.H. acknowledges support from the Swiss National Science Foundation under grant 200020_215634. Part of this work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. Calculations were performed on the Horus cluster of the Division of Space Research and Planetary Sciences at the University of Bern.

References

Appendix A Comparison of existing GI population syntheses

The table in this appendix (adapted from Schib (2021)) gives a comparison of several GI population syntheses conducted to date in terms of the physics included in the models. The different studies use different prescriptions for migration. The following symbols will be used in the table: $$t_{\text{I},\text{FR}}={\left(\frac{M_{\text{p}}}{M_{*}}\Omega \right)}^{-1}\frac{H}{r},$$(A.1) $$t_{\text{I},\hspace{0.17em}\text{N}}=\frac{M_{*}^2}{M_{\text{p}}{M}_{\text{d},\text{loc}}}{\Omega }_{\text{p}}^{-1}$$(A.2) where $$M_{\text{d},\text{loc}}=\pi \Sigma r^2$$(A.3) $$t_{\text{II},\text{FR}}=\frac{1}{\alpha \Omega }{\left(\frac{H}{r}\right)}^{-2},$$(A.4) and M_p is the companion mass. The updated model from Forgan et al. (2018) is not included explicitly in the table. It is the only model that includes gravitational interaction between the clumps besides the model presented here. Other than that, the only difference to Forgan & Rice (2013) is the gap opening criterion given in the explanations below.

Here we give some additional explanations for the following Table B.2:

• 0

  The updated model in Forgan et al. (2018) uses instead the gap opening condition from Crida et al. (2006) with an additional condition for gap crossing from Malik et al. (2015).

• 1

  In addition to Crida et al. (2006), the authors investigate two additional conditions related to the gap crossing and gap formation timescale.

• 2

  with α = 0.005

• 3

  fragments are randomly space within the unstable region

• 4

  disc is unstable at the beginning of the simulation

• 5

  fragment is inserted at the initial fragment location

• 6

  Alexander & Armitage (2007) for UV and Owen et al. (2012) for X-ray

• 7

  External photoevaporation as in Clarke (2007)

• 8

  All models use the ‘first core/follow the adiabats’ approach described in Sect. 4 of Nayakshin (2015a).

• 9

  An N-body integrator is used in the updated model in Forgan et al. (2018).

Appendix B Infall locations

Table B.1 lists the infall location Ri of the baseline population DIPSY-1 for each of the 100 logarithmically spaced bins. The bins are for the final stellar mass M_*,final.

All Tables

All Figures

	Fig. 1 Occurrence of the different possible fragment fates for all fragments formed in the baseline population. From left to right: Collision and merger with another clump, ejection from the system, thermal disruption, accretion on the primary, survival of the fragment, tidal disruption.
In the text

Fig. 2 Mass vs. semi-major axis for all the surviving companions of the baseline population. The colour code of each companion gives the (final) mass of the primary (host star). Differences in occurrence rate due to the different weights of stellar masses according to the IMF are not reflected in this figure; instead, for each of the 100 stellar mass bins, ≈1000 systems containing a total of 7400 companions are shown. The number of systems (Nsyst), the number of surviving companions (Ncomp), the number of objects accreted on the primary (Nacc), and the number of ejected objects (Nejec) are given in the bottom left corner. The horizontal background colours indicate how the companions might be classified as stars, BDs, or planets according to their mass. Clearly, the distinction between these classes is not as clear-cut as these limits may suggest. Objects at very large separations (⪆ 1000 au) may not survive long enough to be observed, as the grey area suggests, as they can become unbound.

In the text

Fig. 3 Initial mass function (Chabrier 2005) for the range of host star masses studied. The objects at the lower mass end are BDs, and the stellar spectral types from M to B are coloured1. Bins in final host star mass used in our study are shown as vertical dashed lines (only every tenth bin is shown for visibility). This demonstrates the range in primary mass covered in our study as well as the logarithmic binning.

In the text

	Fig. 4 Mass-semi-major axis diagram, weighted according to ξ(log m) (Chabrier-2005 IMF, see text Sect. 5.2.2). With this weighting, the population is completely dominated by M-dwarfs (note the difference in Fig. 2).
In the text

Fig. 5 Distributions of four key properties for the companions of the baseline population DIPSY-0, separately for M/K-, G/F-, and A/B host stars, weighted according to ξ(log m) within the three bins, for systems with M*,final > 0.2M⊙. Top left: CMF. For comparison, ξsys(log m) is shown as a thin black line (dashed part indicates extrapolation). Top right: Distribution of semi-major axes. Bottom left: Distribution of eccentricities. Bottom right: Orbital inclinations.

In the text

	Fig. 6 Companion multiplicity for different spectral types (for systems with at least one surviving companion).
In the text

Fig. 7 Mass vs semi-major axis for surviving companions of the baseline DIPSY-0 population for different spectral types of the primary. The spectral type, from M to B, is indicated in the title of each panel. Each panel contains the companions from a random selection of 5000 systems from the corresponding mass range. In the bottom left corner, the following quantities are given as mean numbers for every sub-population: surviving companions per system nc, objects accreted on the primary per system nacc, ejected objects per system nejec, and number of surviving companions per system with at least one surviving companion, ncs (companion multiplicity). The envelope of points has a similar shape for all spectral classes but systematically shifts in position and density.

In the text

	Fig. 8 Fractions of systems that are gravito-turbulent. These develop spirals (‘GI fraction’) that fragment and have at least one surviving companion. The fractions are given for groups of spectral types according to the final host star mass M*,final. All values of M*,final 0.05 M⊙ to 5 M⊙ are included.
In the text

Fig. 9 Quantities as a function of M*,final. Top left: Mean final companion mass, Mc. The thin black line shows the best-fit power law. Top right: Median final semi-major axis, Sc. The thin black line shows the best-fitpower law. Bottom left: Mean final eccentricity. Bottom right: Mean number of surviving fragments (left y-axis) per fragmenting system. The mean number of fragments formed initially (i.e. emerging in total during the disc’s lifetime) is shown in orange (right y-axis), also per fragmenting system.

In the text

Fig. 10 Left: Distribution of the ratio, vc/ve, for all clump-clump (non-collapsed objects only) collisions (see text Sect. 5.9.1) Right: Distribution of the impact angle, θ, for all clump–clump collisions. The solid orange line shows only collisions between clumps, while the thin dotted black lines shows those between collapsed objects. The expected distribution for random orientations (dotted red line) as well as for a flat system (dashed green line) are also shown. Further explanations are given in Sect. 5.9.2.

In the text

Fig. 11 Mass vs semi-major axis (like Fig. 2) for the variant populations DIPSY-1 to 6 as indicated in the panel titles. The baseline population DIPSY-0 is shown in grey in the background. Please note that the top left panel shows the result of an experimental calculation in which some systems are unphysical because the secondary mass is large compared to the primary mass – a situation not correctly captured by our model. Apart from this aspect, the variant populations are remarkably similar to the baseline one, indicating that internal processes dominate the outcome.

In the text

	Fig. 12 CMF for population DIPSY-2 (similar to the top left panel of Fig. 5), with no gas accretion and a higher initial fragment mass. The thin black line shows the baseline-CMF for G- and F-stars for comparison.
In the text

	Fig. 13 Example of the formation of a system with an inner more massive and an outer, less massive companion. The lines show the growth tracks of the clumps. Arrows to the right denote ejections, and the ‘T’ stands for the (ultimate) thermal destruction.
In the text

Fig. 14 Graphical summary of the results of the baseline population DIPSY-0 for F/G host stars in terms of the frequency of different outcomes. On the left, the basic level (star, disc, clump, and companion) is shown in brown. The rest shows a decision-tree-like flow chart with the different outcomes. The colour coding is taken from Figs. 2 and 1 for the mass ranges of the surviving companions and the fates of the destroyed clumps, respectively. A central outcome, the frequency of bound companions inside of 1000 AU, is shown in blue bold face. Please note that the numbers are subject to model assumptions.

In the text