© The Authors 2025

1 Introduction

Planet formation occurs in discs around young stars and two main pathways for this process have been proposed: core accretion (CA, Safronov 1972; Pollack et al. 1996; Hubickyj et al. 2005; Alibert et al. 2005) and disc instability (DI, Kuiper 1951; Mizuno et al. 1978; Cameron 1978; Boss 1997, 2003). In the case of DI (also known as gravitational instability, GI), a cir-cumstellar disc fragments under its own gravity to form one or several gravitationally bound objects. These so called ‘clumps’ have masses of the order of 1 Jovian mass, sizes of ~ 1 au (astronomical unit), and they are close to hydrostatic equilibrium. The subsequent evolution of the clumps determines whether they survive and become compact companions. Companions formed this way can lie in the planetary, brown dwarf, or stellar mass regime.

This includes a number of physical processes that are discussed in detail below.

Originally, DI was proposed as a model for the formation of the Solar System (Kuiper 1951). It was later superseded by CA as the standard model for planet formation, due to a number of (potential) limitations. We review these limitations in the following and discuss their relevance.

One of these limitations is very rapid inward migration. Baruteau et al. (2011) studied the migration of massive planets in gravito-turbulent discs in 2D hydrodynamic simulations, including self-gravity. They found that these objects migrate to the inner disc on a type I timescale, without opening a gap. It was therefore concluded that these planets are unlikely to remain at large separations, at least for the underlying model assumptions (simplified cooling prescription and only a single planet per system). More recently, Rowther & Meru (2020) performed 3D self-gravitating smoothed particle hydrodynamics (SPH) simulations of single migrating planets in gravitationally unstable discs. They applied a variable cooling parameter to account for the expected longer cooling time in the inner disc. The authors found that this effect enables planets to open a gap as soon as they reach the gravitationally stable region in the inner disc; thus, they are able to survive on long timescales. Both studies, however, neglected the effect of additional companions. Gravitational interactions between these companions could also facilitate the survival of at least some objects (Emsenhuber et al. 2021a).

Another related limitation is tidal destruction by Roche lobe overflow (e.g. Kratter & Murray-Clay 2011). When a clump moves close to the primary due to migration or it approaches is periapsis on an eccentric orbit, it can be tidally disrupted (see also Sect. 3.6). Whether this occurs depends on the clump evolution (size and mass). As a clump contracts (Sect. 3.4, its temperature rises. Once the central temperature exceeds ≈2000Κ, the clump undergoes a dynamical collapse (corresponding to the second collapse in star formation), which makes it stable against tidal disruption (Bodenheimer 1974). The pre-collapse timescale depends strongly on the clump’s mass, with massive clumps collapsing sooner (Helled & Bodenheimer 2011). Gas accretion can also shorten this timescale (Vazan & Helled 2012). Furthermore, tidal destruction can be avoided if the planet does not move too close to the primary, which can happen if it opens a gap (e.g. Rowther & Meru 2020). Similarly, a clump could become unbound and dissolve again due to irradiation from its environment (i.e. the thermal bath effect, Vazan & Helled 2012). This fate could be avoided if the clump contracts fast enough.

An additional limitation of DI sometimes brought forward is that it would not be able to form a heavy-element core. While this is not specifically required, core formation is certainly a possibility in a clump (Decampli & Cameron 1979; Boss 1998; Baehr & Klahr 2019). Helled & Schubert (2008) investigated the sedimentation of silicate grains inside isolated clumps of different masses. They found that protoplanets with masses below 5M₄ (Jovian masses) allow for grain sedimentation and core formation, while more massive ones are hotter and in such cases, core formation via grain settling would not be possible. However, cores in more massive clumps could be a result of fragmentation near spiral arms where dust is accumulated (Baehr & Klahr 2019). The accretion of heavy elements post formation could further support this process (Helled et al. 2006).

Another important challenge in studying DI is the computational cost. The huge density contrast between the circumstellar disc and the clump’s interior is very challenging numerically when studying fragmenting discs (Boley et al. 2010; Galvagni et al. 2012). Advanced hydrodynamical codes are expected to allow for progress in this area in the future (see e.g. Matzkevich et al. 2024).

However, none of these limitations alone are expected to prevent the formation of companions in DI in a general way. Still, at present, these limitations in our understanding of several physical processes occurring during DI are significant. One way to make progress in to create more realistic model companion formation in the DI scenarios in the future by combining models of the relevant physical processes and including their interplay into a global model. This output can then be confronted with observational results. Differences between these global model predictions and the observations can then be used to constrain the underlying specialized models and eventually improve our understanding. This is our aim in constructing the model presented in this work.

In recent years, DI has again received growing interest as the number of planetary mass companions on wide orbits has increased (Kratter et al. 2010b; Nielsen et al. 2019). Other classes of planets that are difficult to explain in the classical CA paradigm have also been identified. These include giant planets around very low-mass stars, such as GJ 3512 b (Morales et al. 2019; Kurtovic et al. 2021; Mercer & Stamatellos 2020), and massive planets in very young systems, possibly including AB Aur b (Currie et al. 2022). Overall, CA has been studied extensively, also on a population level (Ida & Lin 2004; Emsenhuber et al. 2021b, see also the reviews in Benz et al. 2014; Mordasini 2018; Burn & Mordasini 2024). Population synthesis projects also exist in the disc instability paradigm (Forgan et al. 2018; Nayakshin & Fletcher 2015). However, these projects are often less mature and/or do not include all the relevant physical processes. Disc instability is expected to mainly populate the brown dwarf (BD) and stellar mass regime (e.g. Rafikov 2005; Kratter et al. 2010a,b; Rice et al. 2015; Xu et al. 2025). The exact shape of the population formed in DI is still unknown. In particular, it is unclear to what extent DI has the capacity to form planetary-mass objects.

A numerical model that includes as many of the relevant physical processes as possible is needed to answer these questions. In particular, it is important to include the early phase of disc formation in a collapsing molecular cloud core (MCC). Here, we present a new global model for Disc Instability Population SYnthesis (DIPSY) that includes a number of improvements over previous global models. While being rich in physical mechanism, it is low-dimensional numerically, so that it can still be used to make quantitative statistical predictions that can be compared to observations. We built on previous works, where an earlier version of the model was used in population syntheses of circumstellar discs (Schib et al. 2021, 2023). The current large model update introduces clumps into the disc if the conditions for fragmentation are satisfied and follows their evolution in detail.

This is the first of a series of papers. In Schib et al. (2025) (Paper II), we applied this model to perform a large-scale population synthesis. The current paper is organised as follows. In Sect. 2, we describe the model for the primary star and the cir-cumstellar disc, including the formation stage through infall, the transport of angular momentum, the treatment of self-gravity, the evolution of the central star, and the temperature model. Section 3 details the treatment of fragmentation as well as the evolution of the clumps. This includes the initial clump mass, mass growth by gas accretion, clump contraction, disc irradiation of the clumps, and mass loss. In Sect. 4, we discuss the exchange of angular momentum with the disc (leading to orbital migration and eccentricity+inclination damping) as well as gravitational interactions among clumps (N-body). Section 5 demonstrates the capabilities of the model by presenting and discussing several systems of varying complexity. In Sect. 6, we discuss the assumptions and limitations of the model. We present a summary and our conclusions in Sect. 7.

2 Disc model

The disc model yields a description for the formation and evolution of the circumstellar disc and forms a central part of the overall global model. Figure 1 gives a schematic overview of the processes occurring in the disc. As illustrated in the figure, the model includes the time evolution of the disc and the primary star. The star and the disc interact via accretion, irradiation, and internal photoevaporation. The disc further interacts with the environment by infall from the MCC and external photoevaporation.

Adaptations of the disc model used in this work were already employed in previous work: Schib et al. (2021, 2022, 2023). We refer to these papers as S21, S22, and S23. Here, we use a modified and extended version, as described below. While we focus on the new aspects of the model, for completeness, we describe all the relevant elements of the model. We used cylindrical polar coordinates, with r denoting the radial direction (distance from the star). We considered a single primary star of (variable) mass M_* at the origin. The heliocentric disc’s mid-plane is located at z = 0au. The disc is assumed to be rotationally symmetric and is described via the vertically integrated surface density, Σ ≡ Σ(r, t), which evolves with time t (Sect. 2.1). The simulations do not start from a fully formed star-and-disc system. Instead, the star and disc are initialised with a negligible mass and are fed by infalling material from the MCC (Sect. 2.5). The star is evolved according to pre-calculated stellar evolution tracks (Sect. 2.4).

The disc is assumed to be subject to internal EUV as well as external FUV photoevaporation. Together with accretion, this ultimately leads to the disc’s dispersal (Sect. 2.7).

Fig. 1 Schematic representation of the star and disc system, including its physical processes.

2.1 Disc evolution

The evolution of the disc’s surface density Σ is calculated by solving the diffusion equation (Lüst 1952 and Lynden-Bell & Pringle 1974): $$\frac{\partial \Sigma }{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left[3r^{1/2}\frac{\partial }{\partial r}\left(v\Sigma r^{1/2}\right)-\frac{2\Lambda \Sigma }{\Omega }\right]+S.$$(1)

The quantity ν denotes the kinematic viscosity. We applied the alpha-parametrisation (Shakura & Sunyaev 1973, Sect. 2.2). Then, Λ is the torque density distribution used to describe the two-way exchange of angular momentum between the disc and companions (Sect. 4.1). The angular frequency, Ω ≡ Q(r, t), includes a contribution from the disc’s auto-gravitation (Sect. 2.3). S ≡ S (r, t) stands for a source and sink term: $$S\left(r,t\right)=S_{\inf }-S_{\text{evap}}-S_{\text{acc}}+S_{\text{loss}}.$$(2)

The contributions to S are: S_inf the source term for the infalling material from the MCC (Sect. 2.5), S_evap the sink term for photoevaporation (Sect. 2.7), S_acc the sink term for mass accreted by companions (Sect. 3.3) and S_loss the source term for mass returned to the disc by companions undergoing mass loss (Sect. 3.6). Equation (1) is solved using the implicit donor cell advection-diffusion scheme with piecewise constant values from Birnstiel et al. (2010). We used a grid with 2800 logarithmically space grid points extending from 3 × 10⁻² au to 3 × 10⁴ au. At the inner edge, the disc is truncated at 0.05 au. The mass flowing across the truncation radius is added to the star, with 10 % considered lost in outflows (Vorobyov 2010).

2.2 Viscosity

The kinematic viscosity ν is parametrised as (Shakura & Sunyaev 1973): $$v=\alpha \frac{c_{\text{s}}^2}{\Omega },$$(3)

where c_s is the isothermal sound speed: $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu \text{u}}}.$$(4)

Here, k_B denotes the Boltzmann constant, T is the disc’s mid-plane temperature (Sect. 2.6), μ ≡ 2.3 is the mean molecular weight and u is the atomic mass unit. For the calculation of α, we considered different regimes as described below.

2.2.1 Global instability of the disc

At early times, the disc is supplied by infalling material from the MCC and can become comparable in mass to the star. In this regime, a global gravitational instability can occur (e.g. Harsono et al. 2011). Kratter et al. (2010a)¹ performed a parameter study of rapid accretion in a series of grid-based 3D hydrodynamic simulations. Their results indicated that the disc’s behaviour can be characterised by two dimensionless parameters, ξ and Γ. They parametrised α in the regime in question as follows: $${\alpha }_{\text{d}}=\frac{1}{18{\left(2-k_{\Sigma }\right)}^2{\left(1+l_j\right)}^2}\frac{{\xi }^{3/2}}{{\Gamma }^{1/3}}.$$(5)

We applied this prescription during the infall phase and follow the authors by setting k_Σ = 3/2 and l_j = 1. The parameters ξ (relating the mass accretion rate of infall to the disc’s sound speed) and Γ (relating the accretion timescale of the infalling gas to its orbital timescale) are given by: $$\xi \frac{{\dot{M}}_{\text{in}}G}{c_{\text{s}}^3}\text{and\hspace{0.17em}}\Gamma =\frac{{\dot{M}}_{\text{in}}}{m_{*d}{\Omega }_{\text{in}}}.$$(6)

Here, Ṁ_in is the infall rate, M_*d the combined mass of the star and the disc, and Ω_in the angular frequency of the infalling gas at the location where it hits the disc. The sound speed, c_s, is also evaluated at the same location; ξ and Γ are calculated at each time step in our model. This makes α_d constant in r but varying in time.

2.2.2 Transport of angular momentum by spiral arms

Once the infall phase has ended, the prescription from Sect. 2.2.1 is no longer applicable (ξ = 0). However, the disc can still be sufficiently massive to be self-gravitating. In this regime, the transport of angular momentum is governed by spiral arms. We applied the following parametrisation (Zhu et al. 2010; Kimura et al. 2016): $${\alpha }_{\text{GI}}=\exp \left(-Q_{\text{Toomre}}^4\right),$$(7)

by setting α globally and using the minimum of Q_Toomre (given in Eq. (28)) across the disc. This treatment may somewhat over-predict the transport of angular momentum through spiral arms. We found that using a radially varying prescription (see also Kratter et al. 2008) leads to pile-ups of material outside of the unstable regions, which can cause an extreme number of fragmentation event (given that we already observed ~100 events in some cases). Therefore, in this work, we kept the assumption we used in Schib (2021); Schib et al. (2023) However, this is an important topic that needs further investigation.

2.2.3 Background viscosity

In the absence of global instabilities or spiral arms, we applied a background viscosity, α_BG. In principle, the model can handle any reasonable value for α_BG. In the simulations presented in Sect. 5 we applied a value of α_BG = 10⁻² throughout the disc in order to reproduce observed disc lifetimes and for consistency with S21; S23 and Paper II. This choice should integrate all mechanisms for the transport of angular momentum, such as effects of disc winds (Bai & Stone 2013; Turner et al. 2014; Suzuki et al. 2016; Weder et al. 2023) or hydrodynamic instabilities like the vertical shear instability (see e.g. Urpin & Brandenburg 1998; Klahr et al. 2023). The value of α_BG may seem on the higher side compared to values on the order of 10⁻³ seen in the literature (Flaherty et al. 2017, 2018). We discuss this topic further in Sect. 6.1. In Appendix A, we study two additional calculations with a lower background viscosity.

In summary, we set α as the maximum of a gravitational α_G and α_BG, $$\alpha =\max \left({\alpha }_{\text{G}},{\alpha }_{\text{BG}}\right),$$(8)

where α_G is equal to α_d during the infall and to α_GI afterwards.

2.3 Auto-gravitation

Since our discs are comparable in mass to their host star at early times, the disc’s self-gravity (auto-gravitation) must be considered. We followed the approach of Hueso & Guillot (2005) and assumed that the disc is in vertical hydrostatic equilibrium, $$\frac{1}{\rho }\frac{\text{d}p}{\text{d}z}=-\frac{G{M}_{*}}{r^3}z-4\pi G\Sigma .$$(9)

The right-most term is the contribution from the disc’s gravity. Then, ρ ≡ ρ(r, z) and p ≡ p(r, z) are the volume density and pressure, respectively. The disc self-gravity leads to a modification of the angular frequency relative to the Keplerian case, $$\Omega \left(r,t\right)={\left[\frac{G{M}_{*}}{r^3}+\frac{1}{r}\frac{\text{d}{V}_{\text{d}}}{\text{d}r}\right]}^{1/2}.$$(10)

In Eq. (10), V_d is the gravitational potential of the disc given by (Hueso & Guillot 2005): $$V_{\text{d}}\left(r\right)={\displaystyle {\int }_{R_{*}}^{\infty }-G\frac{4\text{K}\left[-\frac{4r/r_1}{{\left(r/r_1-1\right)}^2}\right]}{\left|\left(r/r_1-1\right)\right|}}\Sigma \left(r_1\right)\text{d}{r}_1,$$(11)

where R_* is the stellar radius, K the elliptic integral of the first kind. The solution of Eq. (9) becomes $$\rho \left(r,z\right)={\rho }_0\left(r\right)\exp \left(-\left(\frac{\left|z\right|}{H_0}+{\left(\frac{z}{H_1}\right)}^2\right)\right),$$(12)

with ρ₀(r) ≡ ρ(r, 0), where $$H_0=\frac{c_{\text{s}}^2}{4\pi G\Sigma }$$(13)

is the contribution of the disc’s self-gravity and $$H_1=\frac{\sqrt{2}{c}_{\text{s}}}{\sqrt{\frac{\text{G}{M}_{*}}{r^3}}}$$(14)

is that of the star. Requiring $$\Sigma \left(r\right)={\displaystyle {\int }_{-\infty }^{\infty }\rho \left(r,z\right)\text{d}z}.$$(15)

Thus, in our convention, we have $${\rho }_0\left(r\right)=\frac{\Sigma \left(r\right)}{\sqrt{2\pi }H},$$(16)

which leads to the expression of the vertical pressure scale height, $$H=\frac{H_1}{\sqrt{2}}\exp \left(\frac{H_1^2}{4H_0^2}\right)\left(1-\text{erf}\left(\frac{H_1}{2H_0}\right)\right).$$(17)

We refer to Sect. 2.2 of S21, Sect. 3.4.1 of Hueso & Guillot (2005) as well as Huré (2000) and Paczynski (1978) for additional comments and subtleties of this treatment of autogravitation.

The disc’s auto-gravitation also has important implications for embedded companions. They are also subject to the additional acceleration resulting from V_d and will no longer orbit on Keplerian trajectories (see Sect. 4.2). Therefore, it would not make sense to use the usual orbital elements, since these assume Keplerian orbits. Instead, we usually use the separation, r, from the origin to describe a companion’s position. A companion on an eccentric Keplerian orbit (with semi-major axis a and eccentricity e) will have a deficit in angular momentum with respect to an equal mass companion with the same semi-major axis on a circular orbit. The angular momentum is reduced by a factor $\sqrt{1-e^2}$. We can thus define an ‘eccentricity’ that works in the presence of auto-gravitation by demanding that the angular momentum deficit relative to a circular² orbit remains the same.

2.4 Stellar evolution

We applied the stellar evolution tables from Yorke & Bodenheimer (2008) where the radius and luminosity of an isolated protostar are tabulated. We interpolated these tables linearly for a given mass and age. In addition to the star’s intrinsic luminosity, L_int, we included a contribution, L_acc, from accretion of disc material onto the star to calculate the total luminosity, $$L_{*}=L_{\mathrm{int}}+L_{\text{acc}},$$(18)

where L_acc is given by $$L_{\text{acc}}=f_{\text{acc}}\frac{G{M}_{*}{\dot{M}}_{*}}{2R_{*}},$$(19)

with the stellar accretion rate, Ṁ_*, that is given by the disc model. The efficiency factor of stellar accretion heating, f_acc, is set to 1/12 in order to reproduce on average the distribution of observed luminosities of class 0 systems (Tobin et al. 2020, S23).

2.5 Infall

Our simulations were initialised with a seed star and disc. Then the disc was fed by gas infalling from a collapsing MCC. A number of simple prescription how to calculate the source term of infalling material in a spherically symmetric setting exist. Examples are the classical Shu collapse model (Shu 1977) or the Bonnor-Ebert sphere model (Ebert 1955; Bonnor 1956). More complex infall models have been developed on the basis of these ideas, such as the TSC model (Cassen & Moosman 1981; Terebey et al. 1984). While these models provide a relatively simple and intuitive description of the disc formation process, they largely neglect its complex and turbulent nature. Hydrodynamic simulations of star formation show that high angular momentum impacts the disc early, something that is not seen in inside-out collapse. Magnetic fields add another level of complexity that we discuss in the following.

To take into account the chaotic nature of the star formation process, we based our source term on the hydrodynamic population synthesis study of discs by Bate (2018). The paper analyses the evolution and properties of discs formed in a 3D hydrodynamic star cluster simulation. Stellar masses, disc masses, and disc radii are given as a function of time for 183 synthetic protostars. In S21 we discuss in detail how we select a sample of these systems in order to construct probability distributions for initial stellar mass, disc mass, infall rate and disc size. The first three of these quantities are considered correlated: we used multi-variate distributions to obtain the initial values for these quantities. This data was then used as input quantity to perform a population synthesis of forming discs and analyse their properties. As found in S21 this baseline run ‘hydro’ leads to massive (~0.3 M_⊙) and large (~200au) discs. Almost half of these discs were found to fragment. However, magnetic fields, which were not considered in Bate (2018), are thought to alter the disc formation process. Discs are found to be smaller in magnetohydrodynamic (MHD) collapse simulations than in pure hydrodynamic simulations, though the specifics depend on which aspects of magnetic fields are considered. In order to investigate this effect, another run based on MHD was performed in S21. There, we imposed much smaller infall radii (the location where the infalling material hits the disc) in such a way that the distribution of early disc sizes agrees with the prediction of an MHD simulation of disc formation (Hennebelle et al. 2016). The discs formed in this simulation are less massive (~0.1 M_⊙), much smaller (~40au) and none of them fragment. Interestingly, the sizes of observed Class 0 discs (Tobin et al. 2020) lies between the early disc size in ‘hydro’ and ‘MHD’ runs. This led us to perform an additional set of simulations with an additional modification of the infall radii set to fit the observed Class 0 disc sizes. Together with the choice of f_ace (Eq. (19)), this leads to the run ‘OBS_REDIRR’ presented in S23, which exhibits a distribution of both early disc radii and luminosities in agreement with observations. In other words, in this work we use infall radii that lead to discs with sizes that agree with the observations of Tobin et al. (2020). The infall radii found in this way are discussed in Paper II. The source term for the infall S_inf is given by (Eq. (32) of S21): $$S_{\inf }\left(r,t\right)=\frac{{\dot{M}}_{\text{in}}}{2\sqrt{2}{\pi }^{3/2}{R}_{\text{i}}{\sigma }_{\text{i}}}\exp \left[-\left(\frac{r-R_{\text{i}}\left(t\right)}{\sqrt{2}{\sigma }_{\text{i}}}\right)\right].$$(20)

Here, R_i is the infall location and the width σ_i is chosen as R_i/3.

2.6 Temperature model

Our temperature model is based on the vertical energy conservation of the disc which is in hydrostatic equilibrium (Eq. (9)). We considered the effects of the following processes: viscous heating, irradiationfromthe central star(including anaccretionterm, Eq. (18)), shock heating from the gas falling from the MCC onto the disc and a constant background irradiation from the environment. We follow Nakamoto & Nakagawa (1994); Hueso & Guillot (2005) and consider both optically thick and an optically thin regimes and assume an energy balance at the disc’s surface, $$\sigma T_{\text{S}}^4=\frac{1}{2}\left(1+\frac{1}{2\Sigma {\kappa }_{\text{p}}}\right)\left({\dot{E}}_v+{\dot{E}}_{\text{s}}\right)+\sigma T_{\text{env}}^4+\sigma T_{\text{irr}}^4,$$(21)

where T_S is the disc’s surface temperature. This leads to an expression for the disc’s mid-plane temperature: $$\begin{array}{lll}\sigma T_{\text{mid}}^4\hfill & =\hfill & \frac{1}{2}\left[\left(\frac{3\Sigma {\kappa }_{\text{R}}}{8}+\frac{1}{2\Sigma {\kappa }_{\text{P}}}\right){\dot{E}}_v+\left(1+\frac{1}{2\Sigma {\kappa }_{\text{P}}}\right){\dot{E}}_{\text{s}}\right]\hfill \\ \hfill & \hfill & +\sigma T_{\text{env}}^4\hfill \\ \hfill & \hfill & +\sigma T_{\text{irr}}^4.\hfill \end{array}$$(22)

In Eq. (22), κ_R ≡ κ_R(ρ₀, T_mid) and κ_P ≡ κ_Ρ(ρ₀, T_mid) are the Rosseland and Planck mean opacities. Then, κ_P and κ_R are calculated based on Malygin et al. (2014) for the gas and Semenov et al. (2003) for the dust. The calculation is done analogous to Marleau et al. (2017, 2019), assuming ‘normal silicate’ dust grains made of homogeneous spheres (per the NRM model in Semenov et al. 2003). We assume a dust-to-gas ratio of 0.01. In Eq. (22), ${\dot{E}}_v=\Sigma v\left(T_{\text{mid}}\right){\left(r\frac{\text{d}\Omega }{\text{d}r}\right)}^2$ is the viscous heating term, ${\dot{E}}_{\text{s}}=S_{\inf }{\left(r\Omega \right)}^2/2$ the shock heating term due to infall (Kimura et al. 2016). The temperature contribution due to stellar irradiation, T_irr, is calculated as in Hueso & Guillot (2005): $$T_{\text{irr}}=T_{*}{\left[\frac{2}{3\pi }{\left(\frac{R_{*}}{r}\right)}^3+\frac{1}{2}{\left(\frac{R_{*}}{r}\right)}^2\left(\frac{\text{d}\ln \left(H\right)}{\text{d}ln\left(r\right)}-1\right)\right]}^{1/4},$$(23)

where the effective stellar temperature T_* is calculated as $$T_{*}={\left(\frac{L_{*}}{4\pi R_{*}^2\sigma }\right)}^{1/4},$$(24)

with L_* from Eq. (18). In Eq. (23), we set d ln(H)/d ln(r) ≡ 9/7 (Chiang & Goldreich 1997; Fouchet et al. 2012).

2.7 Photoevaporation and disc dispersal

We assumed that the discs are subject to external and internal photoevaporation. Our model of external photoevaporation is based on Matsuyama et al. (2003), which considers FUV irradiation by massive stars. We applied the following sink term, $$S_{\text{ext}}\left(r,t\right)=\{\begin{array}{ll}{S}_{\text{wind}}\left(1-\frac{1}{1+s{m}_{\text{ext}}^{20}}\right),\hfill & \text{if\hspace{0.17em}}r>0.1{\beta }_{\text{M}}{r}_{\text{mI}}\hfill \\ 0,\hfill & \text{otherwise,}\hfill \end{array}$$(25)

with a smoothing term $s{m}_{\text{ext}}=\frac{r}{{\beta }_{\text{M}}{r}_{\text{g,ext}}}$ and $r_{\text{g,ext}}\left(t\right)=G{M}_{*}\left(t\right)/c_{\text{s,ext}}^2$ the gravitational radius. We set c_s,ext = 2.5 km s⁻¹, β_M = 0.14 and S_wind = 2.8 × 10⁻⁸ g cm⁻¹ yr⁻¹, giving an evaporation rate of 10⁻⁸ M_⊙ yr⁻¹ for a disc that extends to 1000 au. This value is relatively low, though comparable with previous work (Armitage et al. 2003; Mordasini et al. 2009).

For the internal photoevaporation, we closely follow Clarke et al. (2001), an EUV model based on the ‘weak stellar wind’ case studied in Hollenbach et al. (1994). We used a sink term given by $$S_{\mathrm{int}}\left(r,t\right)=2c_{s,\mathrm{int}}{d}_{\text{r}}\text{u}s{m}_{\mathrm{int}},$$(26)

with a sound speed c_s,int = 11.1 km s⁻¹, $$d_{\text{r}}=1.8\times {10}^4{\left(M_{*}\left(t\right)/M_{\odot }\right)}^{-0.25}{\left(\frac{r_{\text{g,int}}}{{10}^{14}}\right)}^{-1.5}{\left(\frac{r}{r_{\text{g,int}}}\right)}^{-2.5}.$$

sm_int is a smoothing factor of $$s{m}_{\mathrm{int}}=\{\begin{array}{ll}1-{\left(1+{\left(\frac{r}{0.14r_{\text{g,int}}}\right)}^{20}\right)}^{-1},\hfill & \text{if\hspace{0.17em}}r>r_{\mathrm{int}}\hfill \\ 0,\hfill & \text{otherwise,}\hfill \end{array}$$(27)

where r_int = 0.07 r_g,int, $r_{\text{g,int}}\left(t\right)=G{M}_{*}\left(t\right)/c_{\text{s,int}}^2$.

Fig. 2 Schematic overview of the physical processes related to fragmentation and the evolution of companions that are included in the global model.

3 Fragmentation and clump evolution

In this section, we discuss the topic of fragmentation, and how we handled it in the model. This includes the formation of a clump and its subsequent evolution. Figure 2 gives a schematic overview of the relevant processes. When the disc fragments, one or several clumps are inserted. The clumps can accrete gas (Sect. 3.3) and evolve in their internal structure according to pre-calculated clump evolution tracks (Sect. 3.4) under the influence of disc irradiation (Sect. 3.5). Mass loss of clumps is discussed in Sect. 3.6. Clumps can undergo orbital migration (Sect. 4.1) and interact gravitationally with each other (Sect. 4.2).

3.1 Fragmentation

Fragmentation is the breaking up, due to self-gravity, of parts of the disc into bound objects which we call clumps or fragments. The degree to which a disc is self-gravitating can be quantified by the Q_Toomre parameter (Toomre 1964): $$Q_{\text{Toomre}}=\frac{c_{\text{s}}\kappa }{\pi G\Sigma }.$$(28)

It results from the stability analysis of a thin, differentially rotating fluid sheet. A Q_Toomre < Q_crit ~ 1 permits exponentially growing modes. In Eq. (28), κ denotes the epicyclic frequency: $$\kappa =\sqrt{4{\Omega }^2+2r\Omega \frac{\text{d}\Omega }{\text{d}r}},$$(29)

equal to Ω only in Keplerian discs. For fragmentation, Q_Toomre < Q_crit is not sufficient. Such discs can remain in a state of marginal stability, where the transport of angular momentum is dominated by spiral arms (Hohl 1971; Bertin & Romeo 1988; Lodato & Rice 2004; Cossins et al. 2009; Forgan et al. 2011; Rice 2016). A disc fragments only if saturation mechanisms fail to prevent the growth of perturbations. There is a wealth of literature on fragmentation, including Goldreich & Lynden-Bell (1965); Gammie (2001); Matzner & Levin (2005); Lodato & Rice (2005); Kratter & Matzner (2006). An extensive review can be found in Kratter & Lodato (2016). In our context, it is important to consider two regimes of fragmentation.

3.1.1 Infall-dominated regime of fragmentation

During the infall phase, material from the MCC can reach the disc at a rate greater than it can be transported away. If this occurs, the disc invariably fragments (Boley 2009). We call this the ‘infall-dominated’ regime. The ξ and Γ parameters introduced in Sect. 2.2.1 can be employed to check if we are in this regime. The condition is given in Kratter et al. (2010a) (see also Sect. 3.6 of Kratter & Lodato 2016): $$\Gamma <\frac{{\xi }^{2.5}}{850}.$$(30)

In our model, we assume the disc fragments if Q_Toomre < 1 in the infall-dominated regime.

3.1.2 Cooling-dominated regime of fragmentation

If the condition from Eq. (30) is not satisfied (in particular after the infall phase), the disc fragments only if it can cool efficiently. Gammie (2001) performed shearing sheet simulations of a razor thin disc to demonstrate that this is the case, if the cooling timescale is less than a few times the orbital timescale, following the expression $$t_{\text{cool}}\lesssim {\beta }_{\text{c}}{\Omega }^{-1},$$(31)

a condition known as the Gammie criterion. β_c ≈ 3 is the critical cooling parameter as determined by Deng et al. (2017), see also Baehr et al. (2017). We employ the following cooling timescale (Mordasini et al. 2012): $$t_{\text{cool}}=\frac{3\gamma \Sigma c_{\text{s}}^2}{32\left(\gamma -1\right)\sigma T^4}{\tau }_{\text{eff}},$$(32)

where τ_eff = κ_RΣ/2 + 2/(κ_RΣ) is an effective optical depth, γ ≡ 1.45 is the adiabatic index. We therefore assume that the disc fragments if Q_Toomre < 1 and Eq. (31) are satisfied simultaneously.

In summary, our condition for fragmentation is as follows: Q_Toomre < 1 and either Eq. (30) or (31) must be satisfied. These conditions are checked throughout the disc at every time step.

3.2 Initial fragment mass and insertion of fragments

If the disc fragments, the initial fragment mass M_F (given in Eq. (36)) is removed from the disc at the location (annulus), where Q_Toomre is minimal at the time of fragmentation. (For population synthesis calculations, the procedure is slightly different, as explained in Sect. 4.4 of Paper II). The clump’s properties (aside its mass and age) are determined by means of evolution tracks, as described in Sect. 3.4. A clump of the same mass is inserted at this location r_c. In order to avoid unphysical effects by removing mass from the disc instantly, we instead removed the initial fragment mass on a fiftieth of an orbital timescale. We found that varying this value over reasonable ranges does not strongly affect the outcomes. We apply the following sink term in Eq. (2): $$S_{\text{acc}}\left(r\right)=\frac{{\dot{M}}_{\text{insert}}}{\sqrt{2\pi \delta }}\exp \left(-\frac{{\left(r-r_{\text{c}}\right)}^2}{2\delta }\right),$$(33)

with Ṁ_insert = 50M_FO(r_c)/(2π). For $\sqrt{\delta }$, we chose $2\sqrt[3]{M_{\text{F}}/3M_{*}}$ (the exact value of this width is not very important).

How massive emerging clumps are is important, since the pre-collapse timescale has a strong mass dependency (Sect. 3.6). A rough estimate for the initial fragment mass is the Toomre mass (Nelson 2006): $$M_{\text{Toomre}}=\frac{\pi c_{\text{s}}^4}{G^2\Sigma }.$$(34)

Forgan & Rice (2011) performed more detailed calculations and derived a measure for the local Jeans mass inside the spiral structure of a self-gravitating disc. In our conventions (note the comments in Sect. 2.9.2 in S21) it is equivalent to $$M_{\text{J,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 }_{\text{cool}}^{-1}}},$$(35)

where β_cool = t_coolΩ. The cooling timescale t_cool is given in Eq. (32). A different estimate of the initial fragment mass is given in Boley et al. (2010). The authors performed a calculation that considers the density perturbation near the corotation of a spiral arm. They initial clump mass in the context of a fragmenting spiral arm was found to be $$M_{\text{F}}=\frac{1.6c_{\text{s}}^3}{G\Omega }.$$(36)

We observe that there is a substantial difference between these three estimates, with M_F being by far the lowest. For a Keplerian disc with Q_Toomre = 1, we have $$M_{\text{F}}\approx 0.2M_{\text{J}}^{\text{Sphere}}\approx 0.16M_{\text{Toomre}}\approx 0.04M_{\text{J,FR}},$$

where $M_{\text{J}}^{\text{Sphere}}$ is the Jeans mass given in Eq. (13) of Nelson (2006). Boley (2009) performed 3D SPH simulations with radiative cooling and found good agreement of their analytic expression (M_F) with the simulation results. This has also been investigated by Tamburello et al. (2015), albeit in a different context. These authors performed simulations of isolated galaxies using an N-body + SPH code and also found initial fragment masses comparable to or even lower than M_F. A possible interpretation of the much higher value of M_J,FR is that it does not represent the mass of a clump at the time of fragmentation, but instead, after a few orbits, where it could have accreted a substantial amount of gas. When including magnetic fields, we note that the gravitational instability dynamo may lead to even lower fragment masses (Deng et al. 2021; Kubli et al. 2023). Given these results, we use M_F as the nominal initial mass for our clumps. The influence of a higher initial fragment mass is investigated in Paper II.

During the short period where the fragment is set up, it is kept at a constant separation and is not allowed to interact with the disc (except for the mass removal). No other fragments are allowed to form inside of a Hill radius (for M_F) during this time.

3.3 Gas accretion

We applied the model for gas accretion as presented in S22. This model is based on the Bondi and Hill accretion regimes investigated in D’Angelo & Lubow (2008) (DL08). These authors performed 3D nested grid hydrodynamic simulations of migrating planet undergoing rapid gas accretion. They estimated the accretion rate onto the companion in their Eq. (9), expressed as $${\dot{M}}_{\text{c}}~\frac{\Sigma }{H}\Omega R_{\text{f}}^3.$$(37)

Here, R_f is the radius of the feeding zone of a companion of mass M_c at a separation r_c. It is taken to be the smaller of either the Bondi radius $R_{\text{B}}=G{M}_{\text{c}}/c_{\text{s}}^2$ or the Hill radius³ $R_{\text{H}}=r_{\text{c}}\sqrt[3]{M_{\text{c}}/\left(3M_{*}\right)}$. The accretion rates in the Bondi and Hill regime then become: $${\dot{M}}_{\text{B}}=C_{\text{B}}\Omega \frac{\Sigma r_{\text{c}}^2}{M_{*}}{\left(\frac{r_{\text{c}}}{H}\right)}^7{\left(\frac{M_{\text{c}}}{M_{*}}\right)}^2{M}_{\text{c}},$$(38) $${\dot{M}}_{\text{H}}=\frac{1}{3}{C}_{\text{H}}\Omega \frac{\Sigma r_{\text{c}}^2}{M_{*}}\left(\frac{r_{\text{c}}}{H}\right),$$(39)

where C_B and C_H are dimensionless coefficients of order unity. DL08 found that the accretion rate onto the protoplanet agrees well with min (Ṁ_B, Ṁ_H) as long as the disc is able to supply enough gas. Once the local reservoir is depleted, the accretion rate drops.

In S22, we derived an accretion model based on Eq. (37). Instead of using global values of Σ and Ω, we calculate the contributions from each grid cell inside the feeding zone separately. The gas is then removed self-consistently from the disc at the location from where it was accreted. We obtain the following accretion rate (Eq. (D.5) in S22): $$\begin{array}{l}{\dot{M}}_{\text{B,H}}=\sqrt{\pi }{C}_{\text{B,H}}{\displaystyle {\int }_{r_{\text{c}}-R_{\text{f}}}^{r_{\text{c}}+R_{\text{f}}}{\rho }_0\left(r\right)H_1\exp \left(\frac{H_1^2}{H_0^2}\right)}\\ \times \left(\text{erf}\left[\frac{\sqrt{R_{\text{f}}^2-{\left(r-r_{\text{c}}\right)}^2}}{H_1}+\frac{H_1}{2H_0}\right]-\text{erf}\left[\frac{H_1}{2H_0}\right]\right){\upsilon }_{\text{rel}}\text{\hspace{0.17em}d}r.\end{array}$$(40)

In Eq. (40), C_B,H denotes the numerical factor, calibrated with the results from DL08 to be C_B = 10 and C_H = 0.19, respectively in S22. The quantities H₀ and H₁ are related to the disc’s auto-gravitation and are described in Sect. 2.3. υ_rel = |rΩ – r_cΩ_c| is the velocity of the gas relative to the companion.

In S22, we compare the accretion model, together with our migration model with the results from two different hydro-dynamic simulations of accreting, migrating companions. We typically find a reasonable agreement for a range of parameters.

3.4 Clump evolution tracks

Just after formation, clumps are extended, au-sized objects of low density. They will contract and heat up on a timescale that strongly depends on their mass. During the time that a clump is young and extended, it is prone to tidal destruction. However, if the interior reaches temperatures of ~2000 K, the clump undergoes a dynamical collapse due to the dissociation of molecular hydrogen (Bodenheimer 1974). The radius then shrinks by about three orders of magnitude to a few Jovian radii (R₄). In order to account for this evolution, we employed interior evolution tracks of isolated clumps. These tracks are the result of 1D calculations of clumps published in Humphries et al. (2019). The code used solves the standard stellar structure equations on a fully implicit, adaptive grid (Helled et al. 2006, 2008; Vazan & Helled 2012; Vazan et al. 2013, 2015). It assumes a protosolar hydrogen-helium composition and the Saumon et al. (1995) equation of state. Heat transport inside the clump is modelled by convection, conduction or radiation, depending on the local conditions. In the radiative regions, opacities from Pollack et al. (1985) are used. The clumps evolve in isolation (i.e. with a zero pressure boundary condition), which means that no exchange of mass or irradiation are considered. Tracks are available for a number of masses, ranging from 0.5 M₄ to 12M₄. Since the clump masses in our simulations evolve with time, we linearly interpolate the quantities from the evolution tracks. This allows us to determine the clumps’ (and later compact objects’) properties, such as radius and temperature, at any time during their evolution. Figure 3 shows the time evolution of the radius, central density and central temperature for a selection of clump masses. Also shown is the pre-collapse time for all masses. We define the pre-collapse time as the age of the clump when its central temperature reaches 2000 K. At this time, the clump radius plummets, while the central density and temperature increase.

Fig. 3 Evolution tracks of isolated clumps. Top left: clump radius. Top right: central density. Bottom left: central temperature. Bottom right: precollapse time for all clump masses. These tracks were first published in Humphries et al. (2019).

3.5 Clump irradiation

Our evolution tracks were calculated for isolated clumps. However, when clumps form, they are embedded in the surrounding circumstellar disc. The clumps are therefore irradiated by the surrounding disc material. This influences their evolution, and this effect is not considered in the tracks.

In the absence of tabulated tracks including irradiation we follow an approximate approach inspired by Humphries et al. (2019)⁴. The authors assume that, because of the irradiation, the clump’s evolution is slowed down: The clump’s ‘internal time’, τ, passes more slowly than the physical time t depending on the incident radiation. For the clump’s total energy Ε we have: $$\frac{\text{d}E\left(\tau \right)}{\text{d}\tau }=-L_0\left(\tau \right),$$(41)

where L₀(τ) is the clumps radiative luminosity in isolation. If the clump is embedded in the disc, its energy balance is modified by the incident irradiation from the disc, and at time t it must hold: $$\frac{\text{d}E\left(t\right)}{\text{d}t}=-L\left(t\right)=-L_0\left(t\right)+4\pi R_{\text{c}}^2\left(t\right){\sigma }_{\text{B}}{T}_{\text{mid}}^4\left(t\right),$$(42)

where the right-most term denotes the radiation incident on the clump from the disc (R_c is the clump’s radius and T_mid is the disc’s mid-plane temperature evaluated at the clump’s location).

The clump’s radius and luminosity for a given clump mass and age are linearly interpolated from the evolution tables (Sect. 3.4). Equation (42) expresses that, in order to infer the clump’s effective luminosity L, its luminosity in isolation must be reduced by the incident luminosity from the disc. In order to find the relationship between t and τ, we divide Eq. (42) by Eq. (41) and obtain $$\frac{\text{d}E\left(t\right)}{\text{d}t}{\left(\frac{\text{d}E\left(\tau \right)}{\text{d}\tau }\right)}^{-1}=\frac{-L\left(t\right)}{-L_0\left(\tau \right)}=\frac{\text{d}\tau }{\text{d}t}=\frac{L_0\left(t\right)-4\pi R_{\text{c}}^2\left(t\right){\sigma }_{\text{B}}{T}_{\text{mid}}^4\left(t\right)}{L_0\left(\tau \right)}.$$(43)

In the second equality we required dE(t) = dE(τ) since we want to evaluate the evolution tables at a time where the change in total energy of the isolated case is the same as in the irradiated case. Equation (43) allows us to infer τ for a given time t. For example for τ₀ = t₀ = 0 s, at t₁, we get $${\tau }_1={\tau }_0+\frac{\text{d}\tau }{\text{d}t}\left(t_1-t_0\right)=\frac{\text{d}\tau }{\text{d}t}{t}_1.$$(44)

Knowing the relationship between t and τ allows us to calculate the properties of the irradiated clump: for a clump with internal age τ₁ the evolution tables are evaluated at $\frac{\text{d}\tau }{\text{d}t}{t}_1$.

This approach is a reasonable approximation when T_mid is much lower than the clump’s effective temperature T_eff; however, it becomes problematic if the two reach comparable values. For example, if they are the same, Eq. (43) gives τ = const. In other words, the clump stops evolving. In a realistic scenario, we would expect the clump to adapt to a different interior structure, possibly with a larger radius or losing some of its mass. The situation becomes even more problematic if T_mid exceeds T_eff. In practice, we kept τ constant in this case, and instead compared T_mid with the clump’s central temperature. If the disc’s temperature exceeds even that, we assumed that the clump is thermally destroyed (Cameron et al. 1982; Vazan & Helled 2012). Once available from dedicated evolutionary simulations, we will use evolution tracks including irradiation in future iterations of the global model.

3.6 Mass loss

When a clump moves closer to the star, its Hill radius R_H decreases. At some point, the clump’s physical radius given by the tracks can become larger than its Hill radius. We assumed that in this case the clump loses the mass outside of its Hill sphere (Roche lobe overflow). In our model, R_c is compared with R_H at each time step. If R_c > R_H, the clump’s mass is reduced by a small amount and R_c is redetermined from the evolution tables. The new radius is again compared with R_H for the reduced mass. This process is repeated until R_c ≤ R_H. In some cases, no stable configuration is found, which implies that the clump was entirely disrupted (tidal destruction). Any mass lost in this process is returned to the disc at the clump’s location via a source term (S_loss in Eq. (1)). The same happens if a clump is destroyed due to thermal disruption (Sect. 3.5).

4 Clump–disc and clump–clump interactions

4.1 Gas disc-driven orbital migration

A companion embedded in a circumstellar disc interacts gravitationally with the disc gas. The exchange of angular momentum leads to orbital migration (Goldreich & Tremaine 1980). Classically, gas disc migration has been divided into two main regimes. In the type I regime, the surface density remains largely unperturbed by the presence of the companion. In a non-self-gravitating disc, this is the case for planets up to a few to a few tens of Earth masses (Μ_⊕). The exact value depends on the stellar mass and disc properties. Type I migration is complex and is strongly affected by the disc’s thermodynamics. Some of the contribution to the total torque comes from Lind-blad resonances and is typically negative (inwards) (Goldreich & Ward 1973). Additional contributions come from the co-rotation region and these can be both negative and positive, such that the net migration can be both inwards and outwards (Paardekooper & Mellema 2006). An important aspect in determining whether outward migration is possible is the saturation of the co-rotation torque at higher masses, since the saturation works against outward migration (Paardekooper et al. 2011). More massive companions start to perturb the surface density significantly until they clear a gap in the vicinity of their orbit. They are then said to migrate in type II. The simple picture consisting of type I and type II migration is helpful when studying planets formed in CA. Such planets are assembled from initial masses of ≪ 1 Μ_⊕ that need to accrete mass for many orbits to eventually open a gap. Prescriptions for the type I migration (e.g. Tanaka et al. 2002; Paardekooper et al. 2010, 2011) can then be applied until some criterion for gap opening is fulfilled (Crida et al. 2006; Kanagawa et al. 2018). Above gap-opening masses, type II prescriptions can then be applied (Dittkrist et al. 2014; Kanagawa et al. 2018).

In DI, the situation is different. Fragments are born with high masses ~M₄ and they may migrate too fast to carve a deep gap despite their high masses. Furthermore, migration is expected to be to some degree stochastic (Baruteau et al. 2011; Kubli et al. 2025). The conventional formulae for type I migration assuming a steady state cannot be applied, and the conditions for gap opening are still uncertain (Malik et al. 2015; Müller et al. 2018; Rowther & Meru 2020). Furthermore, while the transfer of angular momentum from a ~Μ_⊕ planet on the disc gas may be negligible, this is not the case for a clump of several M₄.

4.1.1 Torque densities

In order to overcome these limitations, we applied the migration model from S22. The migration model uses torque densities to account for the two-way exchange of angular momentum between the disc and the companions. The torque density distribution per unit disc mass, $\frac{dT}{dm}{\left(r\right)}_i$, for companion i is defined such that $${\mathcal{T}}_i=2\pi {\displaystyle {\int }_0^{\infty }{\left(\frac{\text{d}\mathcal{T}}{\text{d}m}\right)}_i\left(r\right)\Sigma \left(r\right)r\text{\hspace{0.17em}d}r},$$(45)

where T_i is the total torque on companion i: the sum of the contributions from all radial locations in the disc. In the evolution equation (Eq. (1)) we thus have for n companions: $$\Lambda ={{\displaystyle \sum _{i=1}^n-\left(\frac{\text{d}\mathcal{T}}{\text{d}m}\right)}}_i\left(r\right).$$(46)

In principle, this formalism fully integrates the exchange of angular momentum between the disc and the companions in our model. Hereafter, we drop the subscript i and consider the case of an individual planet for readability.

We note, however, that two important questions remain unanswered regarding (1) the shape of $\frac{dT}{dm}\left(r\right)$ and (2) what happens with planets on eccentric orbits.

For the first question, possible choices are the classical impulse approximation (Lin & Papaloizou 1979a,b, 1986) or an improved formalism (Armitage et al. 2002). These approaches can work reasonably well in a pure type II regime where a gap has already been fully opened. However, they do not work well in type I or in the transition between type I and type II because the approach explicitly excludes co-rotation torques. As discussed above, co-rotation torques are important there. Fortunately, torque densities are well suited for the type I regime as well – provided they include co-rotation torques. A model that accomplishes this is presented in D’Angelo & Lubow (2010). The authors performed 3D nested grid hydrodynamic simulations of locally isothermal discs. They provided a functional form of the torque density distribution parametrised by the disc’s surface density and temperature gradients, $$\frac{\text{d}\mathcal{T}}{\text{d}m}\left(r\right)=\mathcal{F}\left(x,\beta ,\zeta \right)\Omega {\left(r_{\text{c}}\right)}^2{r}_{\text{c}}^2{\left(\frac{M_{\text{c}}}{M_{*}}\right)}^2{\left(\frac{r_{\text{c}}}{H\left(r_{\text{c}}\right)}\right)}^4,$$(47)

where $$\mathcal{F}\left(x,\beta ,\zeta \right)=\left[p_1{e}^{-{\left(x+p_2\right)}^2{p}_3^2}+p_4{e}^{-{\left(x-p_5\right)}^2{p}_6^2}\right]\tanh \left(p_7-p_{8}x\right)$$(48)

is a dimensionless function describing the shape of the torque density. In Eq. (48), β = –d ln Σ/d ln r and ζ = –d in T/d ln r, $x=\frac{r-r_{\text{c}}}{\max \left(H,R_{\text{H}}\right)}$ is a scaling factor. The parameters p₁ to p₈ depend on β and ζ and set the amplitude and width of the torque density. They are determined by D’Angelo & Lubow (2010) via fits to their simulations. This form of the torque density is adequate to the type I regime, D’Angelo & Lubow (2010) study 1 Μ_⊕-planets in most of their simulations. However, they also explored what happens in the case of gap-opening planets. They showed that when the planetary mass increases, the amplitude (scaled) torque density decreases, although only by a factor of order unity.

4.1.2 High mass torque

In addition to the decreasing amplitude, as the planetary mass approaches 1 Μ₄, an inversion of the torque density starts to appear near the position of the planet (their Fig. 15). The exact shape of this inversion, which depends strongly on the planet’s mass is, however, not very important. What is important is that the torque is small near the planet. In S22 we studied this effect by applying a torque model named the ‘high-mass torque model’ resulting from an interpolation of Fig. 15 in D’Angelo & Lubow (2010) to accreting, migrating planets. We argued that the inversion of the torque density prevents premature gap opening and compared the migration tracks of these 1D simulations with the 3D nested grid hydrodynamic simulations of DL08, finding reasonable agreement for a range of parameters. We also showed that without the modifications related to the planet’s mass, a gap opens too fast and migration slow down too quickly. While the torque model performs well in comparison with locally isothermal 3D simulations, it cannot reproduce features such as the outward migration that depend on the (non-)saturation of the co-rotation torque. We propose a simple method to overcome this limitation in Appendix F of S22. For the present study, these effects should not be important, since we are dealing with massive objects for which the co-rotation torque saturates. Therefore, we applied the high mass torque from S22 in Eq. (45). Since we also applied this shape of the torque density to higher-mass companions, we prevented the inversion close to the companion by setting $\frac{\text{d}T}{\text{d}m}=0$ near the companion (Hallam & Paardekooper 2017). We note that despite our treatment of the torque densities at high masses, their application to 1D discs leads to gap shapes that are different from what is seen in hydrodynamic simulations. This is seen in simulations running for hundreds of orbits and the gaps in 1D tend to become too deep and too narrow. We discuss this issue and its possible consequences in Sect. 6.2 and will address this topic in future research.

4.1.3 Eccentric and inclined orbits

The second question raised in Sect. 4.1.1 about companions on eccentric orbits is particularly important for multiple systems. Even if all companions (like newly formed clumps) were initially on (approximately) flat, circular orbits, the mutual gravitational interaction of the companions tends to increase their eccentricity and inclinations. Furthermore, in the context of DI, we expect the clumps to be on initially eccentric orbits, since they are born near the corotation of spiral arms.

Strictly speaking, the torque density formulation discussed in Sect. 4.1.1 is only applicable for nearly circular, flat orbits. The same is true for the type I migration formulae mentioned in Sect. 4.1. These must be modified (Cresswell & Nelson 2008; Bitsch & Kley 2010; Fendyke & Nelson 2014; Coleman & Nelson 2014). The treatment in the context of torque densities is not trivial.

To our knowledge, no published torque densities for eccentric, inclined orbits exist (but see Fairbairn & Dittmann 2025). However, simply applying, for example, the torque density corresponding to a companion’s semi-major axis to the disc when the companion is on a strongly eccentric orbit violates the conservation of angular momentum. Numerical experiments show that doing so tends to add angular momentum to the system, which can lead to unphysical fragmentation events and other complications. We stress that such difficulties arise only if the companion’s torque is applied on the disc through Eq. (1). An additional complication arises in relation to the damping of eccentricities and inclinations. This damping is expected as long as the companions are embedded in a gas disc. For massive companions in type II migration, hydrodynamic simulations show that damping occurs on timescales shorter than the migration timescale (Kley et al. 2004; Kley 2019). In this situation, the so-called K-damping is often applied (Lee & Peale 2002), where the damping timescale is set to 1/K times the migration timescale. K is typically of the order of ten Emsenhuber et al. (2021a). The damping is then applied to the companion as an additional acceleration in the N-body integrator. In the context of our model, this becomes problematic, since it means applying a torque or force to a companion without applying the opposite equal to the disc. Again, this leads to unphysical effects in the disc evolution.

For this reason, we developed a more robust way to model migration and damping for massive companions. It is based on the formulation of planet-disc interaction in the framework of dynamical friction presented in Ida et al. (2020). The authors studied planet-disc interactions in the subsonic and supersonic regimes and compared and discussed in detail various approaches to calculate damping found in the literature. Based on this discussion and on the results of a hydrodynamic simulation, they provided a simple prescription for the damping timescales of semi-major axis, eccentricity and inclination for low-mass planets. In order for this approach to be applicable in our model, we need to demonstrate that, in the limit or circular orbits, we retrieve the same torque as the one obtained from the high mass torque (Sect. 4.1.2). Furthermore, it also needs to be applicable for massive companions.

Ida et al. (2020) assumed discs with smooth surface densities in their Sect. 4 to arrive at equations of motion that can be expressed in terms of the migration timescale (their Eq. (46)). In order to apply the dynamical friction formalism to gap-opening companions, we need to relax this assumption. Instead, we start from their Eqs. (29) and (32). Their Eq. (29) gives the (specific) force from dynamic friction in the subsonic limit (Tanaka & Ward 2004): $${\mathit{F}}_{\text{DF,sub}}=-0.780\frac{\Delta \mathit{\upsilon }}{t_{\text{wave}}},$$(49)

with the characteristic (inverse) timescale of $$t_{\text{wave}}^{-1}=q\frac{{\Sigma }_{\text{c}}{r}_{\text{c}}^2}{M_{*}}{h}^{-4}{\Omega }_{\text{c}},$$(50)

we have the mass ratio of the companion to the star, $$q=\frac{M_{\text{c}}}{M_{*}},$$(51)

the disc aspect ratio, $$h=\frac{H\left(r_{\text{c}}\right)}{r_{\text{c}}},$$(52)

the torque $${\Gamma }_{\text{o}}={\left(\frac{q}{h}\right)}^2{\Sigma }_{\text{c}}{r}_{\text{c}}^4{\Omega }_{\text{c}}^2,$$(53)

and the velocity of the companion relative to the gas, $$\Delta \mathit{\upsilon }={\mathit{\upsilon }}_{\text{c}}-{\mathit{\upsilon }}_{\text{g}},$$(54)

where: $${\upsilon }_{\text{c}}=\left(\begin{array}{c}{\upsilon }_{\text{c,r}}\\ {\upsilon }_{\text{c,}\theta }\\ {\upsilon }_{\text{c,z}},\end{array}\right),$$(55)

and, assuming an unpertured disc velocity field, $${\mathit{\upsilon }}_{\text{g}}=\left(\begin{array}{c}{\upsilon }_{\text{g,r}}\left(r_{\text{c}}\right)\\ {\upsilon }_{\text{c}}{\Omega }_{\text{g}}\left(r_{\text{c}}\right)\\ 0\end{array}\right).$$(56)

These are the velocities of the companion and the gas (at the companion’s separation), respectively, in cylindrical coordinates. For Δυ ≡ |Δυ| we get $$\Delta \upsilon =\sqrt{{\left({\upsilon }_{\text{c,r}}-{\upsilon }_{\text{g,r}}\right)}^2+{\left({\upsilon }_{\text{c},\theta }-r{\Omega }_{\text{g}}\right)}^2+{\upsilon }_{\text{c,z}}^2}.$$(57)

The angular frequency of the gas, Ω_g is given by $${\Omega }_{\text{g}}=\Omega \left(1-\eta \right),$$(58)

with Ω from Eq. (10). The factor $$\eta \equiv -\frac{h^2}{2}\frac{\text{d}\ln p}{\text{d}\ln r}$$(59)

accounts for the deviation of the gas velocity from the purely gravitational value due to the pressure gradient. For smooth discs this effect is very small, but it is important here because otherwise the θ-component of Δυ vanishes for a companion on a circular orbit. For the supersonic limit, Ida et al. (2020) apply the expression from Muto et al. (2011) (their Eq. (32), simplified): $${\mathit{F}}_{\text{DF,sup}}\simeq -2\pi \frac{\Delta \mathit{\upsilon }}{t_{\text{wave}}}{\left(\frac{\Delta \upsilon }{c_{\text{s}}}\right)}^3.$$(60)

Summation of timescales then leads to (as in their Eq. (35)): $${\mathit{F}}_{\text{DF}}=\frac{\Delta \mathit{\upsilon }}{\Delta \upsilon }{\left[\frac{1}{-0.78\Delta \upsilon t_{\text{wave}}^{-1}}+\frac{1}{-2\pi \Delta \upsilon t_{\text{wave}}^{-1}}{\left(\frac{\Delta \upsilon }{c_{\text{s}}}\right)}^3\right]}^{-1}.$$(61)

As discussed in Ida et al. (2020), the dynamical friction formalism can recover the migration rate in the subsonic case, but only to an order of magnitude. This is because in the subsonic regime, the contribution by density waves is important. A way to include this contribution without assuming low masses is to modify Eq. (61). To do this, a relationship between the total torque calculated through Eq. (45) and $t_{\text{wave}}^{-1}$ is needed. We write the total torque on a companion as: $$T=\tilde{K}{\Gamma }_{\text{o}},$$(62)

where $\tilde{K}$ is a dimensionless prefactor. With the migration timescale, τ: $${\tau }^{-1}=\frac{2T}{L_{\text{c}}}=\frac{2\tilde{K}{\Gamma }_{\text{o}}}{M_{\text{c}}{r}_{\text{c}}^2{\Omega }_{\text{c}}},$$(63)

where $$L_{\text{c}}=M_{\text{c}}{r}_{\text{c}}^2{\Omega }_{\text{c}},$$(64)

is the orbital angular momentum of the companion, which can be expressed as $$t_{\text{wave}}^{-1}=\frac{1}{2\tilde{K}{h}^2}{\tau }^{-1}=\frac{{\Gamma }_{\text{o}}}{L_{\text{c}}{h}^2}.$$(65)

If we assume Σ(r) ~ Σ_c in Eq. (45) and substitute from Eq. (47), then formally we can write: $$\tilde{K}=2\pi H^{-2}{\displaystyle {\int }_0^{\infty }F\left(\frac{r-a}{H},\beta ,\zeta \right)rdr.}$$(66)

In summary, we have $$\begin{array}{lll}T\hfill & =\hfill & \tilde{K}{\Gamma }_{\text{o}}=\tilde{K}{L}_{\text{c}}{h}^2{t}_{\text{wave}}^{-1}\hfill \\ \hfill & =\hfill & 2\pi H^{-2}{L}_{\text{c}}{h}^2{t}_{\text{wave}}^{-1}{\displaystyle {\int }_0^{\infty }F\left(\frac{r-a}{H},\beta ,\zeta \right)rdr.}\hfill \end{array}$$(67)

Now Eq. (61) can be modified. We replace the factor –0.78 by $\tilde{K}{h}^2{\eta }^{-1}$. In order to see that this choice is reasonable, we consider the θ-component for a circular orbit. Then we have Δυ/c_s << 1 and Δυ = υ_c,θ – r_cΩ_g = r_cηΩ, and thus: $$\begin{array}{lll}{F}_{\text{DF},\theta }\hfill & =\hfill & {\left[\frac{\eta }{\tilde{K}{h}^2{r}_{\text{c}}\eta \Omega t_{\text{wave}}^{-1}}\right]}^{-1}\hfill \\ \hfill & =\hfill & \frac{\tilde{K}{h}^2{r}_{\text{c}}^2{M}_{\text{c}}{\Omega }_{\text{c}}{t}_{\text{wave}}^{-1}}{M_{\text{c}}{r}_{\text{c}}}\hfill \\ \hfill & =\hfill & \frac{\tilde{K}{L}_{\text{c}}{h}^2{t}_{\text{wave}}^{-1}}{M_{\text{c}}{r}_{\text{c}}}\hfill \\ \hfill & =\hfill & \frac{T}{M_{\text{c}}{r}_{\text{c}}}.\hfill \end{array}$$(68)

In the first step, we use Ω = Ω for a circular orbit. In the second step, we substitute the expression from Eq. (64) and in the third, we drew from Eq. (67). Thus, we retrieve the acceleration of the companion when subject to the torque from Eq. (45). In the general case, we therefore apply the following acceleration (see Sect. 4.2 for the sign of T): $$F_{\text{TD,DF}}=\frac{\Delta \upsilon }{\Delta \upsilon }{\left[\frac{M_{\text{c}}{r}_{\text{c}}}{T}-\frac{1}{2\pi \Delta \upsilon t_{\text{wave}}^{-1}}{\left(\frac{\Delta \upsilon }{c_{\text{s}}}\right)}^3\right]}^{-1}.$$(69)

It includes the contribution to migration by density waves, is not limited to smooth discs and reproduces the supersonic limit.

Equation (69) also illustrates what we have discussed above: in case of circular orbits, the torques on companion and disc are equal but opposite. Conversely, on eccentric orbits, υ_c,θ is typically different from rΩ_g. Therefore, it would not be correct to apply the torque from Eq. (45) to the disc. Instead, when applying the torque to the disc, it is scaled by the sine of the angle between the position vector r_c and the velocity vector, υ_g, $$T_{\text{corr}}=T\frac{r_{\text{c}}\times {\upsilon }_{\text{g}}}{\left|r_{\text{c}}\right|\left|{\upsilon }_{\text{g}}\right|}.$$(70)

This ensures that the torque on the disc is reduced for eccentric orbits and gives the correct limit for radial motion (zero torque on the disc).

We note that our treatment of synchronising the torques between the companion and the disc fully accounts for accelerations only in θ direction. Applying torques in r- and z- direction would make the disc locally eccentric/inclined. This is not possible in our model which corresponds to a rotationally symmetric disc centred in the z = 0 plane.

4.2 Gravitational interaction

In addition to the tidal interaction with the gas disc, companions also interact with each other gravitationally. We employ the mercury N-body code (Chambers 1999) to model such interactions. It is a symplectic integrator that applies a hybrid method in order to correctly handle close encounters. mercury is called in each main time step. Since the companion positions are updated in the N-body code, the result of the companions’ interaction with the gas disc must also be included there. This is done by means of additional (specific) forces. The form of these forces related to gas disc migration as well as i- and e-damping is discussed in Sect. 4.1. Here we give the final expressions that are passed to the N-body code. There is a subtlety to the application of these forces. It is done in cylindric, heliocentric coordinates. This means that the companions’ positions and velocities first need to be transformed from democratic to heliocentric coordinates, and then decomposed into tangential and radial components. After the application, the transformation needs to be inverted. The force in r-direction is $$F_{\text{r}}=\frac{{\upsilon }_{\text{c,r}}+\frac{\dot{m}}{2\pi r{\Sigma }_{\text{c}}}}{\Delta \upsilon }{\left[-\frac{M_{\text{c}}{r}_{\text{c}}}{\left|T\right|}-\frac{1}{2\pi \Delta \upsilon t_{\text{wave}}^{-1}}{\left(\frac{\Delta \upsilon }{c_{\text{s}}}\right)}^3\right]}^{-1}-\frac{d}{dr}{V}_{\text{d}}\left(r\right).$$(71)

The first summand in Eq. (71) is almost as in the r-component of Eq. (69), (we substituted $-\frac{\dot{m}}{2\pi r{\Sigma }_{\text{c}}}$ for υ_g,r, with ṁ ≡ ṁ(r) the gas accretion rate in the disc)⁵. The difference is the sign of the first summand in the square brackets. This term is always negative since it corresponds to the damping, which needs to be in the direction opposite to the companion’s relative motion. The second summand in Eq. (71) is the additional acceleration resulting from the disc’s gravitational potential. In θ-direction we have the force: $$F_{\theta }=\frac{{\upsilon }_{\text{c},\theta }-r_{\text{c}}{\Omega }_{\text{g}}\left(r_{\text{c}}\right)}{\Delta u}{\left[\frac{M_{\text{c}}{r}_{\text{c}}}{T}-\frac{1}{2\pi \Delta \upsilon t_{\text{wave}}^{-1}}{\left(\frac{\Delta \upsilon }{c_{\text{s}}}\right)}^3\right]}^{-1},$$(72)

exactly as in the θ-component of Eq. (69). Migration is typically directed inward (T < 0), although it can also become positive. For the z-direction, we obtained $$F_{\text{z}}=\frac{{\upsilon }_{\text{c,z}}}{\Delta u}{\left[-\frac{M_{\text{c}}{r}_{\text{c}}}{\left|T\right|}-\frac{1}{2\pi \Delta \upsilon t_{\text{wave}}^{-1}}{\left(\frac{\Delta \upsilon }{c_{\text{s}}}\right)}^3\right]}^{-1}-sign\left(z_{\text{c}}\right)2\pi G{\Sigma }_{\text{enc}},$$(73)

with the same sign in the T-term as in Eq. (71). z_c is the companion’s z-coordinate. In Eq. (73), the last summand is the acceleration an inclined companion experiences due to the disc’s gravity. We assumed that this acceleration results from an infinitely extended, thin disc with the surface density Σ_enc enclosed by the companion. We note that this additional term does not damp the inclination of the companion, since it is directed against υ_c,z only if z_c and υ_c,z have the same sign.

4.3 Collisions

When the companions interact gravitationally they can experience close encounters, leading to scatterings, or collisions. In this study, we treated all collisions as perfect mergers: when two companions come closer than the sum of their radii, a perfectly inelastic collision is assumed to occur. This may appear as an extreme assumption since it neglects the possibility of eruptive or hit-and-run collisions that would change the orbits of the bodies in question without completely destroying one of them. This is especially true for clumps due to their large extension. Indeed, it was shown that collisions between clumps can lead to various outcomes depending on the assumed conditions (see Matzkevich et al. 2024 for details). However, we found that most clump-clump collision typically occur below the mutual escape speed. In this regime (see also Wimarsson et al. (2025) for a similar collisional regime in a different context), perfect merging may be a reasonable assumption (e.g. Leinhardt & Stewart 2012; Chau et al. 2021). We study collisions further in Paper II.

Overall, the global disc instability formation model presented here contains – in the low dimensional approximation – a significant number of physical processes that are coupled in a self-consistent way. At the same, it is still a strong simplification of the actual companion formation and evolution process. We discuss important limitations in Sect. 6.

5 Formation of individual systems

Thanks to its low-dimensional character, it is possible to employ the global model for population syntheses. However, before studying the result of the model in a statistical sense in Paper II, it is important to discuss the evolution of several individual systems, studying outcomes of increasing complexity. Thus, we began by looking at the simple case of a system that does not fragment. Then we studied a system with a single fragmentation event that only produces one clump. Next, we looked at a system with two fragmentation events. Finally, we investigated a more complex system evolution with multiple fragmentation events and several interacting clumps.

Fig. 4 Evolution of a system that becomes self-gravitating but that does not fragment. Left panel: stellar mass, disc mass and disc size as a function of time. Right panel: minimum QToomre and α as a function of time. The black vertical line shows the end of infall. It should be noted that there two distinct y-axes in both panels to which the different quantities belong to.

5.1 A non-fragmenting system

Figure 4 shows the time evolution of a system where the disc becomes self-gravitating (i.e. where Q_Toomre approaches unity and spiral waves would form) but where no fragmentation occurs. As we shall see in Paper II, such an outcome is the most probable one: while almost all discs are found to have early on a self-gravitating phase, only 10–20% fragment. Thus, only the disc part of the model is used. This system is initialised at 1 kyr with a stellar seed of ≈6.5 × 10⁻² M_⊙ and a disc of (≈2 × 10⁻² M_⊙). It experiences infall at ≈3 × 10⁻⁵ M_⊙yr⁻¹ until t ≃ 36 kyr. The left panel of Fig. 4 displays the stellar and disc masses as a function of time.

The disc mass is initially less than a third of the stellar mass. Due to the rapid infall, the disc mass increases quickly to two thirds the stellar mass. This leads to a very fast transport of angular momentum via the global instability of the disc (Sect. 2.2.1). The associated α parameter is shown in the right panel of the figure. It reaches almost 0.3 in the early phase. The Q_Toomre parameter, also shown in the right panel, drops to low values, but remains just above unity. After the infall phase terminates (indicated with a thin vertical line), α drops quickly, although it remains elevated above background α_BG for more than 100 kyr. This coincides well with the overall minimum value of Q_Toomre, which is lower than about 1.5 in the same period. We expect such a system to exhibit spiral arms for about 100 kyr. Here, we see the typical auto-regulation of self-gravitating discs: as Q_Toomre approaches unity, spiral waves develop which efficiently transport away mass and angular momentum. In our model, this is expressed by the increase of α reaching a high value of almost 0.3. This increases in turn Q_Toomre via a higher disc temperature and reduced surface density. In this system, fragmentation is prevented this way. It is also a typical behaviour that the overall lowest Q_Toomre is reached directly after the end of infall, because at this moment, the heating (and thus stabilisation) of the disc by shock heating of the disc’s surface by infalling gas vanishes.

5.2 A simple fragmenting system

The next system we investigate has a lower total mass than the one discussed above, although the behaviour of stellar mass and disc mass is qualitatively similar (top left panel of Fig. 5). The initial values are 0.02 M_⊙ and 0.04 M_⊙ for stellar mass and disc mass, respectively, with infall at 4.9 × 10⁻⁵ M_⊙ yr⁻¹ until 5.5 kyr. In this case (with a very high disc-to-star mass ratio) the unstable initial state leads to accretion bursts. This is seen as an oscillatory feature in Q_Toomre and α. We checked that these oscillations are not numerical in nature by rerunning the simulation at a much lower time step. The thin black dash-dotted line in the top right panel of the figure, shows the accretion rate of disc material on the star mdot (in units of 10⁻⁴ M_⊙ yr⁻¹ vs the left y-axis). The accretion rate also shows strong oscillations between 1 and 2 kyr. The final stellar mass is ≈0.2 M_⊙.

In Fig. 5, the end of the infall phase is again marked with a black vertical line. In this case, given the combination of quantities setting Q_Toomre, a value slightly less than unity is reached and a fragment is formed with an initial mass (set according to Eq. (36) based on the local disc conditions) of 0.9 M₄ at 44 au. The middle left panel shows the time evolution of mass (solid line) and separation (dashed line) of the fragment. The red vertical line at ≈6 kyr marks the fragmentation event. Initially, the fragment migrates inwards rapidly, while accreting gas from the disc. At approximately 10 kyr (or after about six initial orbits), the migration is slowed down. At this time, the forming companion has already accreted almost 4 M₄ of gas. The middle right panel shows a zoom-in on the same data with a linear x-axis, revealing more details of what is happening in the early phase. The accretion rate is at its maximum of 10³ M₄ yr⁻¹ during the early phase due to the massive disc and the fast transport. The rapid inward migration is slowed down after ≈ 10 kyr with the onset of the formation of a gap.

The fragment in Fig. 5 starts with a moderate initial angular momentum deficit (corresponding to an eccentricity of ≈0.08, see Sect. 2.3) and does not exhibit strong oscillations in its separation at 6 kyr (middle right panel). To understand the role of different initial orbital eccentricities, we run two variant simulations which only differs by the initial eccentricity. They are shown only in the middle right panel. In the first variant simulation, we studied how the fragment would behave with a higher initial ‘eccentricity’ (see Sect. 2.3) of 0.25, shown by the thin dotted line in the middle right panel. The fragment starts from apocentre and initially moves inwards quickly, but is essentially circularised through damping in a single orbit. The difference in separation after 5 kyr is less than 1 au. In the second variant simulation (dashed black line), a fragment with double that eccentricity (corresponding to 0.5, higher than the maximum used for the population synthesis) is inserted. In this case, cir-cularisation happens in about three orbits and the fragment has moved ≈ 2 au further inwards (at 25 kyr) compared the one with lowest angular momentum deficit. This indicates that the initial eccentricity does not have a major impact on companions formed by DI, at least in the single fragment case. We discuss this further in Paper II.

During the first 3 kyr after fragmentation, the disc receives only little perturbation from the fragment, which migrates quickly. Later, with increasing fragment mass, the disc material is pushed away more vigorously. After ≈4kyr to 6kyr, the combination of tidal interaction and gas accretion has led to a deep gap. This process is shown in the bottom left panel of Fig. 5. It shows the surface density centred on r – r_c for several snapshots in the early migration. With the formation of the gap, the accretion rate starts to decrease due to the depletion of the disc near the companion.

The bottom-right panel in Fig. 5 finally displays the fragment’s internal properties, namely its radius, central temperature, T_c, and the effective temperature, T_eff,c. The radius of the clump increases from an initial 0.6 au to 1 au under the combined influence of accretion and contraction. The central temperature riises from ≈300 K to ≈2000 K. At this point, the clump undergoes a dynamical collapse at 16 kyr. Then, radius drops to 1.7 × 10⁻³ au or about 3.6 R₄ (Jovian radii).

This is an example of a clump that survives the early inward migration due to the fast accretion. If it had stayed at its initial mass, it would have migrated into the hot inner disc long before its pre-collapse time, and would have been tidally or thermally disrupted. Due to the increasing mass, the clump’s Hill radius remains larger than 3 au before 16 kyr. The disc’s temperature is always lower than the clump’s effective temperature, and the clump is not thermally disrupted. However, the time until the second collapse occurs is slightly delayed through the irradiation from the disc (Sect. 3.5). This outcome is aided by the fact that the fragment is born after the end of the infall phase, where the outer disc is colder. Also, the fast accretion is possibly facilitated by the absence of other fragments.

The collapsed companion continues to migrate inwards on a timescale of 10² kyr as seen in the middle left panel. It crosses the disc’s inner edge at 235 kyr at 0.05 au and is considered in the model to be accreted by the primary. It could, however, in principle also become a Hot Jupiter, but we currently do not model the further evolution once a planet enters the stars magnetospheric cavity. The corresponding increase of the primary mass is visible in the top left panel. Observationally speaking, this system contains after the end of (main) infall phase for about 0.2 Myr a massive hot (2000 K) companion migrating through the inner 10 au of the disc.

Fig. 5 Evolution of a system forming a single (transient) fragment which is eventually accreted into the host star because of orbital migration. Top left: stellar mass, disc mass and disc radius. Top right: minimum QToomre and α. The stellar accretion rate is shown in black (see text Sect. 5.2). Middle left: time evolution of fragment orbital separation rc and mass Mc. Middle right: same with a zoom-in on the early phase, two variant calculations are shown additionally in black (see text Sect 5.2). Bottom left: surface density, centred on the fragment (and thus at different absolute orbital distances), for a selection of times. Bottom right: time evolution of the fragment radius, central and effective temperature as well as the disc temperature at the fragment’s location.

5.3 A system forming two fragments

Increasing the complexity further, we next discuss a system where two fragments form. Here, N-body interactions and collisions start to play a role. This system is similar to the one discussed above and is initialised with a 0.03 M_⊙ star and a 0.026 M_⊙ disc. The disc starts less massive than the star, but accretes quickly and its mass reaches that of the star at 2 kyr. The star accretes faster after that, and reaches 0.2 M_⊙ by the end of the simulation. The infall continues up to 5.1 kyr, and at this time the stellar mass is 0.17 M_⊙ and the disc’s mass is 0.11 M_⊙. The evolution of disc and stellar mass is shown in the top left panel of Fig. 6. The disc fragments at 5.5 kyr, i.e. again shortly after the end of infall. In contrast to the system discussed above, two clumps form. The evolution of Q_Toomre and α (top right panel) is qualitatively similar to that discussed before. However, there is more oscillation in both parameters due to the more complex dynamics of the two fragments.

The middle panels of Fig. 6 show the evolution of mass and semi-major axis of both fragments. These clumps do not migrate inwards. Instead, the combination of their mutual gravitational interaction and the interaction with the disc leads to a slow outward motion until the two clumps collide at 8 kyr, i.e. shortly after their emergence. The merged fragment with a combined mass of ≈6 M₄ migrates outwards in the increasing disc with an increasing eccentricity. The magnitude of the torque decreases until it becomes negative and the forming companion begins migrating inwards again starting at ≈ 50 kyr. The eccentricity is damped slowly. The bottom left panel of Fig. 6 shows the surface density, centred on the first fragment for different times. The gap opens less quickly compared to the system from Fig. 5 due to the orbit’s eccentricity. The bottom right panel depicts the evolution of temperature and radius of both fragments. Until the collision, R_c, T_c, and T_eff are increasing slowly due to the increase in the clumps’ mass. The curves corresponding to the second clump are in thin solid line and are hard to see since they are partly covered by the first clump’s lines. The evolution is almost the same for both clumps due to their very similar masses and ages. After the collision, the surviving clump contracts slowly, increasing its temperature, until it undergoes second collapse at 12.3 kyr. While the dynamical evolution in this system is complicated due to the many factors involved, the overall behaviour of the surviving companion is characterised by an early outward migration followed by a slow inward migration. Outward migration occurs because this companion is located at the outer edge of the disc after the interaction with the other clump. There, the positive torque from the inner disc dominates. The first companion continues to accrete gas, and enters the brown dwarf (BD) regime. It carves a deep gap despite being on an eccentric orbit. This prevents material from the outer disc from crossing its orbit and the inner disc depletes. Once the inner disc is depleted, it no longer exerts a torque on the companion and the migration direction is inward. The inward migration becoming slower with decreasing disc mass and stops after 6 Myr as the disc dissolves. The companion survives as 43 M₄ BD with a semi-major axis of 5 au.

5.4 A complex system with multiple fragmentation events

After a fragmentation event, mass is removed from the disc, which drives it towards stability. If the infall phase is sufficiently long, the disc can become unstable again. Furthermore, the tidal interaction of companions with the disc can also drive it towards instability. In this way, multiple fragmentation events, each yielding one or more fragments, are possible. This leads to situations where many fragments are in the disc at the same time. An example of such a system is shown in Fig. 7. The left panel shows the time evolution of the stellar and disc masses. This is a very massive system, with the primary reaching 4.2 M_⊙, the disc mass is 1 M_⊙ at its highest. The infall phase lasts for 75 kyr, which corresponds to a long duration. The disc fragments dozens of times and forms in total (over its lifetime) ~100 fragments. The fragmentation events (sometimes closely spaced in time) are shown as vertical red lines in the figure. The right panel shows the minimum of Q_Toomre as well as α as a function of time. Both quantities show a lot of variations due to the many fragments and their interaction with the disc. This phase lasts until 100 kyr. Just after the end of the infall phase, a last intense burst of fragmentation occurs. Many clumps are scattered to the outer disc, where the majority of the mass resides. There, they accrete much of the gas disc in a short time, as seen in the sudden drop of M_d in Fig. 7 around 100 kyr. This also drives the disc towards stability, with Q_Toomre increasing sharply and α dropping to its background value.

The behaviour of the companions is very complex. Clumps interact gravitationally with each other or collide, or they might become ejected or thermally disrupted. There are up to 30 objects in the disc at the same time.

Fig. 6 Evolution of a system yielding two fragments in a single fragmentation event. The panels are as in Fig. 6, except here, separation, mass, radius, and temperature for two companions (with subscript 1 and 2) are shown in the middle and the bottom right panels. The surface density shown in the bottom left panel is centered on the first companion. In the end, a massive brown dwarf companion at 5 au remains.

5.5 Diversity of outcomes: system architectures, masses, semi-major axes and eccentricities

The systems’ architecture after 100 Myr can be quite different. We discuss these in detail on a population level in PaperII. Here, we look at a number of individual systems to illustrate possible outcomes of the combined processes described above. Surviving companions have masses from ≈ 1 M₄ to 200 M₄, semi-major axes from ≈ 1au to 10 000 au and appear in single, double or triple systems. Figure 8 shows a collection of example systems from the baseline population DIPSY-0 discussed in Paper II. Each panel depicts mass versus semi-major axis of the companions. Their eccentricity is shown as grey bar. Orbital parameters were calculated in Jacobi coordinates. Companion masses are typically on the order of a few tens of Jovian masses as seen in most panels of the figure. Eccentricities can vary from close to zero up to large values. Small eccentricities are often seen when all but one companion are lost (e.g. through collisions or gravitational interactions) early. Any eccentricity of the remaining companion will then be damped through the interaction with the disc. The top left panel shows such a system. Single companions can also have significant eccentricities though, as seen in the rightmost panels of the third and forth row. This typically happens if a companion is scattered out of the disc or if additional companions are only lost after the disc phase. In both cases, the companion’s eccentricity would not be damped. For systems with two surviving companions, masses can increase or decrease with semi-major axis (as seen e.g. in the second row). The masses of the two companions can be similar (like in the lower two panels in the second column or very different (as in system 8387). In systems with two companions, a typical configuration consists of a massive inner companion on a nearly circular orbit with a lower mass outer companion with a high eccentricity. Examples for this are the top two systems in the rightmost row.

Fig. 7 Evolution of a system with many fragments (same as Fig. 4).

6 Model limitations

6.1 Alpha-viscosity

The simulations presented in this paper were run with thevalue of α = 0⁻² for the background alpha viscosity. We did this to reproduce observed disc lifetimes (S21; S23) and for consistency with Paper II. This also allowed us to avoid the additional complexity that arises from a more detailed treatment of the disc’s viscosity.

Here we consider that the magneto-rotational instability (MRI, Balbus & Hawley 1991) might not be the main driver of angular momentum in circumstellar discs (e.g. Lesur et al. 2023). A single value for α_BG throughout the disc is probably an oversimplification independent of the value. However, our value for α_BG could be reasonable in the outer disc regions we are mostly interested in, if the disc is sufficiently ionised there (e.g. Cilibrasi et al. 2023). Global models of viscous discs typically apply alpha-values between 10⁻³ and 10⁻² (Kimura et al. 2016; Emsenhuber et al. 2021b). Using a lower value of 10⁻³ or 10⁻⁴ for α_BG would lead to disc lifetimes that are longer than the observed ones, as discussed in S21 (see also Haisch, Jr. et al. 2001; Mamajek 2009; Emsenhuber et al. 2023). It may also not reproduce observed accretion rates (see Weder et al., in rev.). The inclusion of magneto-hydrodynamic disc winds (see e.g. Suzuki et al. 2016; Weder et al. 2023) or photoevaporation models with much larger mass loss rates (Picogna et al. 2019; Ercolano et al. 2021; Haworth et al. 2018) may lead to higher overall mass loss rates and would enable the use of a lower α_BG with the same disc lifetimes (i.e. with a similar long-term evolution). We note that the effect of wind-driven accretion on orbital migration, an important ingredient of our model, is not yet well understood (Kimmig et al. 2020; Lega et al. 2022; Wafflard-Fernandez & Lesur 2023). In this work, we are mostly interested in the early phases of disc formation and evolution where self-gravity plays a role. We thus kept the model simple in order to avoid additional complexity and also for consistency with S21 and S23. The choice of the background α_BG has, however, no influence on the early evolution of our discs and their fragmentation, since the transport of angular momentum is dominated by global instabilities and spiral arms at this stage (Sects. 2.2.1 and 2.2.2), which in our model correspond to a (significant) increase of the (total) alpha over the background value. In Appendix A, we tested the effect of using a lower value of 10⁻³ for α_BG. The results shown there indicate that companions may end up on wider orbits in this case due to the slower type II migration. In Sect. 6.6 of Paper II, we performed a similar test on a population level.

6.2 Torque densities

The torque densities we applied (4.1.2) were derived from 3D hydrodynamic simulations of planet-disc interaction, as discussed. Our approach is still new and therefore there are uncertainties regarding the migration rate and gap formation. Here, we applied the torque densities also for companions with masses beyond the range in which they were tested. While the results from D’Angelo & Lubow (2010) indicate that the mass dependence is weak, this regime should be investigated in future studies. Maybe more importantly, we applied the torque densities also when the disc is self-gravitating. To our knowledge, no torque densities for self-gravitating discs have been published yet, and this limitation needs to be overcome in dedicated studies. Another important challenge is that, despite the treatment of the torque densities near the companion, gaps seem to become too deep (tend to zero) over long timescales. While this is probably not crucial for this study (since the torque contribution from the gap bottom is very small anyway), it is of key importance for studies that include solids. We are currently working on an improved prescription that should handle these situations. The simulations on which the high mass torque is based are also locally isothermal. A more comprehensive treatment of thermodynamics might not be very important for massive companions that open deep gaps. However, some recent results (e.g. Ziampras et al. 2025) also indicate that thermodynamics plays a key role in the gap-opening process. The influence of e.g. radiation transport in the parameter range of our study should also be investigated in future research.

Fig. 8 Mass versus semi-major axis at the end of the simulation for a selection of systems. The eccentricity is given as a grey bar. The index of the system as well as the final primary mass is given at the top of each panel. See text Sect. 5.5 for further explanations.

6.3 Heliocentric disc model

The basic type of disc model we use assumes that the primary is at the origin, while the disc is revolving around it. For a disc-star system (only) this assumption is reasonable since the disc is assumed rotationally symmetric. As long as any companion in the system is much less massive than the primary, the barycentre would remain close to the origin. However, when studying massive companions as we do in some cases here, this assumption can become problematic. The barycentre can shift away from the primary which leads to torques on the disc that would make it eccentric. Furthermore, the companions’ orbits would be inconsistent with the disc’s dynamics. For example, a companion with a circular orbit around the barycentre no longer remains at the same location in the disc. Such problems arise when a companion accretes at a very high rate for a prolonged period. However, hydrodynamic simulations indicate that accretion rates onto clumps rarely exceed 10⁻³ M₄/yr for any length of time. Shabram & Boley (2013) performed 3D radiation hydrodynamics simulations of massive companions at wide separation as expected to result from fragmentation. They find long-term accretion rates onto the proto gas-giants of up to 3 × 10⁻⁴ M₄/yr. Fletcher et al. (2019) study the migration of massive, wide companions including gas accretion, comparing seven different hydrodynamic codes. They do not find sustained accretion rates above 10⁻³ M₄/yr also for companions with an initial mass of 12 M₄. We therfore limit the accretion rate onto the companions to 10⁻³ M₄/yr. We note, that accretion rates above this value were reported in Zhu et al. (2012). The accretion limit effectively restricts the applicability of our model to systems with a primary that is dominant in mass. The study of, for example, binaries of comparable mass interacting with the disc is not possible. We discuss this restriction further in Paper II.

6.4 Photoevaporation

The prescriptions for internal and external photoevaporation (Sect. 2.7) are standard models that do not yet include some of the latest results on disc photoevaporation (e.g. Ercolano et al. 2021. The main purpose of including photoevaporation in our work is (together with the background alpha viscosity) to regulate the disc lifetimes to observed lengths (see also S21; S23). The exact nature of mass loss is less important, since we are mostly interested at what happens in the outer disc very early on, before photoevaporation becomes a relevant process. However, in Paper II we also investigate the effect of having a higher rate of external photoevaporation.

6.5 Accretion of solids, grain growth and sedimentation

A detailed treatment of solids (dust, pebbles) is currently missing in our model. Solid grains inside a clump can influence its evolution, for example by changing the opacity (Helled & Bodenheimer 2011; Humphries et al. 2019). A higher opacity would likely lead to a longer pre-collapse timescale, making the clumps’ survival less likely and vice versa (see also Baehr & Klahr 2019). On the other hand, grains could grow and sediment to form a core (Boss 1998; Helled et al. 2006, 2008; Helled & Schubert 2008), and grain growth and settling could also lead to a decrease of the clump’s opacity, as discussed in Helled & Bodenheimer (2011), and thus make it reach the second collapse earlier. These are important aspects that we plan to address in future work.

7 Summary and conclusions

In this work, we present and discuss a new global numerical model that can be used to study the formation of a star-and-disc system, including disc fragmentation, clump formation, and clump evolution as well as various interaction processes (e.g. clump-clump and clump-disc). The model represents an improvement over existing similar models by combining several key physical effects not self-consistently linked previously. We used a precursor of this model in previous studies (S21; S22; S23) to investigate the formation and fragmentation of circum-stellar discs, as well as gas accretion and orbital migration of companions. The new version of the model presented here combines the functionalities of previous work and adds a number of physical effects. It can now be used to study the formation and evolution of systems of companions formed in disc instability ranging in mass from the planetary mass regime over brown dwarfs to (low) stellar masses. By using a low-dimensional approach (axisymmetric discs, spherically symmetric clumps and companions, but with 3D N-body interactions) the model is also suitable to perform population syntheses of objects formed in the context of disc instability. The main addition to the model is the inclusion of one or several clumps after fragmentation and their evolution and interaction with the disc as well as with each other. With this extension, the current model consists of three main parts, described below.

The first main part is the model of the circumstellar disc (Sect. 2). It is characterised by the following aspects:

• The evolution of the gas surface density is modelled in the framework of the viscous disc picture in the 1D axisymmetric vertically integrated approximation (Sect. 2.1). Several source and sink terms are included: infall, disc photoe-vaporation, gas accretion and mass loss by the embedded companions;

• The viscosity of the disc is based on the α-parametrisation (Sect. 2.2) and contains time-variable contributions that depend on the disc’s global (Sect. 2.2.1) and local stability (Q_Toomre) parametrising the transport of angular momentum by spiral arms (Sect. 2.2.2) and a constant background contribution (Sect. 2.2.3);

• As disc masses can be comparable to the stellar mass, the discs’ auto-gravitation is taken into account in the calculation of the angular frequency and the vertical scale height (Sect. 2.3);

• As the primary mass increases strongly and its properties evolves, we use tabulated (proto)stellar evolution tracks (Sect. 2.4), where the luminosity also includes the accretional luminosity;

• Star and disc start from negligible initial masses and are then built up by infall from a molecular cloud core (Sect. 2.5). Infall rates are taken from hydrodynamic simulations, while infall radii are chosen to obtain observed Class 0 disc sizes. Infall durations are set to obtain the given final stellar masses;

• The disc’s temperature structure, which (together with the gas surface density and the angular frequency) is key for the stability, is found from energy conservation in the vertical direction. This includes the contributions from viscous heating, stellar irradiation, shock heating by the infalling gas, background irradiation from the MCC, and radiative cooling at the disc’s surface (Sect. 2.6);

• External FUV-driven and internal EUV-driven disc photoevaporation are included, as they are important for the long-term evolution post-infall and the eventual dissipation of the discs (Sect. 2.7).

The second main part of the global model (Sect. 3) deals with the formation and evolution of gaseous clumps forming from the self-gravitational fragmentation of a disc. Here, the following aspects are included:

• By checking the Toomre parameter and cooling criteria (different during the infall and post-infall phases) at each moment and orbital distance in the disc, we can determine whether and where the disc fragments (Sect. 3.1);

• When a disc is found to fragment, one or several fragments are inserted in the fragmenting region. Their initial masses can be set according to several prescriptions found in the literature (Sect. 3.2);

• After formation, clumps can accrete gas from the disc (Sect. 3.3). Thus, we can apply the Bondi+Hill-accretion model presented in S22. This can lead to significant mass growth. The gas accreted by the clumps is removed from the disc to ensure mass conservation;

• The clumps’ interior evolution (cooling, contraction) is modelled by interpolating in pre-calculated evolutionary tables of isolated clumps (Sect. 3.4). We approximately correct for the delaying effect caused by disc irradiation (Sect. 3.5). The evolution tracks include in particular also the second collapse where the clumps turn from very extended (au-scale) low-density objects into much more compact ones (Jovian radius-scale);

• Clumps can lose mass through tides if the approach close to the primary (Sect. 3.6); namely, if their physical size given by the interior evolution tracks becomes larger than the Hill sphere (Roche lobe overflow). Thermal destruction is included in this step as well.

The third main part of the global model comprises the interaction of clumps and companions with the disc and with each other.

• The gravitational interaction with the gas disc is modelled by using torque densities (S22) to take into account the two-way exchange of angular momentum between companions and the disc (Sect. 4.1). This leads to orbital migration and gap opening, as well as the damping of eccentricities and inclinations of the clumps. We have developed a new approach based on dynamical friction (Ida et al. 2020) that is applicable in the context of gap-opening companions;

• The gravitational interaction between companions is calculated by applying an N-body integrator (Sect. 4.2). The N-body interactions (full 3D, in contrast to the other model parts) lead to scatterings, ejections, and collisions with the host star and among the clumps and companions. Collisions are treated in a simplified way and integration continues up to 100Myr to include potential late dynamical events.

While it is possible to apply this new global model for disc instability population synthesis (see Paper II), aspects of the model, such as the treatment of the infall phase or the companion-disc interaction, can also be used in general studies of companion formation via disc instability. Following the description of the model, we showcase it (Sect. 5) by studying a number of systems with increasing complexity, studying the evolution and impact of parameters such as the primary mass and disc mass as well as mass and semi-major axis of companions.

These simulations demonstrate the complexity of modelling planet and companion formation via DI when the effects included in the model interact. We illustrate this by discussing in detail four simulations of increasing complexity. The first example is a system with a disc that becomes self-gravitating, but does not fragment (Sect. 5.1). Then we look at a system which does fragment, but only a single fragment forms that migrates to the inner edge of the disc (Sect. 5.2). The third example is a system in which two fragments form and later collide, with the merged companion surviving and ending up as a ~40 M₄ brown dwarf (Sect. 5.3. Fourth, we look at a complex system that fragments multiple times with strong interaction between the dozens of fragments that emerge in the disc (Sect. 5.4).

Finally, as an outlook with respect to the second companion paper, we discuss the diversity of outcomes in terms of system architectures in the case where there is at least one surviving companion by discussing the masses, semi-major axes, and eccentricities of sixteen different systems (Sect. 5.5). Despite the complexity of the model and the many coupled physical mechanisms included, our model has some important limitations: A fixed and high value for the background alpha-viscosity (Sect. 6.1), uncertainties related to the application of torque densities (Sect. 6.2), the heliocentric disc model (Sect. 6.3), the prescriptions for photoevaporation (Sect. 6.4) and the absence of potential solid accretion, grain growth, and settling which are still missing in the current model (Sect. 6.5).

We conclude that a variety of interconnected physical processes, such as gas accretion, orbital migration, and N-body gravitational interactions, exert a strong influence on the inferred population of forming objects. Furthermore, our findings highlight that the assumptions underlying theoretical models play a pivotal role in determining the final architecture and properties of the systems that form in such conditions. In Paper II of the series, we used the new model presented here to establish a large-scale population synthesis in the disc instability model that allows us to confront the predictions of this formation paradigm with observations in a quantitative way.

Acknowledgements

We thank Lucio Mayer, Allona Vazan and Simon Grimm for the insightful discussions. O.S. and C.M. acknowledge the support from the Swiss National Science Foundation under grant 200021_204847 ‘PlanetsIn-Time’. 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 Lower background viscosity

For the simulations presented in Sect. 5 we used a value of 10⁻² for the background viscosity α_BG. Here, we test the sensitivity of our results to this choice, by applying a lower value of 10⁻³. Fig. A.1 and A.2 depict the time evolution of mass and separation for the systems discussed in Sect. 5.2 (A.1) and 5.3 (A.2). Fig. A.1 shows the system with a single fragment. As expected, the evolution is identical while the viscosity is determined by spiral arms (until ≈30 kyr, compare to top left panel of Fig. 5). The final masses are the same, since the mass evolution is essentially finished at this point. The subsequent evolution of the separation is much slower at lower viscosity due to the longer type II migration timescale. The outcome is the same, though: the companion migrates all the way to the inner edge of the disc and is accreted on the primary. The evolution of the system with a high viscosity is shown with thin black lines for comparison (the mass is barely visible as it is unchanged).

In Fig. A.2, the analogous evolution for the system with two fragments is shown. Again, the early evolution is the same and the long-term migration is much slower. In this case, the surviving companion ends up at ≈50 au.

Fig. A.1 Mass-separation evolution of the system with one fragment at a lower background viscosity. The results from Sect 5 are shown in black.

Fig. A.2 Mass-separation evolution of the system with two fragments at a lower background viscosity. The results from Sect 5 are shown in black.

Appendix B Hill radius for high mass ratios

The expression for R_H given in Sect. 3.3 was derived assuming q ≡ M_c/M_* << 1. It is inaccurate for q ≳ 0.1 (in fact, R_H is greater than 0.5 for M_c = M_*). This is because the expression for R_H results from using the binomial approximation on a quintic equation in the derivation. In order to avoid this difficulty, we use the following polynomial approximation for q > 0.01: $$\begin{array}{l}{R}_{\text{H}}/r_{\text{c}}\approx 0.1171+3.41867q-31.2062q^2+203.082q^3\\ -879.93q^4+2617.62q^5-5534.93q^6+8550.38q^7\\ -9841.34q^8+8546.69q^9-5634.7q^{1}0\\ +2819.48q^{11}-1062.45q^{12}+296.371q^{13}\\ -59.3077q^{14}+8.04718q^{15}-0.662781q^{16}\\ +0.0250073q^{17},\end{array}$$(B.1)

which is accurate to q ≈ 3.

All Figures

	Fig. 1 Schematic representation of the star and disc system, including its physical processes.
In the text

	Fig. 2 Schematic overview of the physical processes related to fragmentation and the evolution of companions that are included in the global model.
In the text

	Fig. 3 Evolution tracks of isolated clumps. Top left: clump radius. Top right: central density. Bottom left: central temperature. Bottom right: precollapse time for all clump masses. These tracks were first published in Humphries et al. (2019).
In the text

	Fig. 4 Evolution of a system that becomes self-gravitating but that does not fragment. Left panel: stellar mass, disc mass and disc size as a function of time. Right panel: minimum QToomre and α as a function of time. The black vertical line shows the end of infall. It should be noted that there two distinct y-axes in both panels to which the different quantities belong to.
In the text

Fig. 5 Evolution of a system forming a single (transient) fragment which is eventually accreted into the host star because of orbital migration. Top left: stellar mass, disc mass and disc radius. Top right: minimum QToomre and α. The stellar accretion rate is shown in black (see text Sect. 5.2). Middle left: time evolution of fragment orbital separation rc and mass Mc. Middle right: same with a zoom-in on the early phase, two variant calculations are shown additionally in black (see text Sect 5.2). Bottom left: surface density, centred on the fragment (and thus at different absolute orbital distances), for a selection of times. Bottom right: time evolution of the fragment radius, central and effective temperature as well as the disc temperature at the fragment’s location.

In the text

Fig. 6 Evolution of a system yielding two fragments in a single fragmentation event. The panels are as in Fig. 6, except here, separation, mass, radius, and temperature for two companions (with subscript 1 and 2) are shown in the middle and the bottom right panels. The surface density shown in the bottom left panel is centered on the first companion. In the end, a massive brown dwarf companion at 5 au remains.

In the text

	Fig. 7 Evolution of a system with many fragments (same as Fig. 4).
In the text

	Fig. 8 Mass versus semi-major axis at the end of the simulation for a selection of systems. The eccentricity is given as a grey bar. The index of the system as well as the final primary mass is given at the top of each panel. See text Sect. 5.5 for further explanations.
In the text

	Fig. A.1 Mass-separation evolution of the system with one fragment at a lower background viscosity. The results from Sect 5 are shown in black.
In the text

	Fig. A.2 Mass-separation evolution of the system with two fragments at a lower background viscosity. The results from Sect 5 are shown in black.
In the text