{'fraction_non_alphanumeric': 0.03210260972716489, 'fraction_numerical': 0.019770660340055358, 'mean_word_length': 4.531027884089666, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The link between infall location, early disc size, and the fraction of self-gravitationally fragmenting discs

Many protoplanetary discs are self-gravitating early in their lives. If they fragment under their own gravity, they form bound gaseous clumps which may evolve to become giant planets. Today, the fraction of discs that undergo fragmentation, and the frequency of conditions that may lead to giant planet formation via gravitational instability, is still unknown. We perform a population synthesis of discs from formation to dispersal. In varying the infall radius, we study the relationship of the early disc size with fragmentation. Furthermore, we investigate how stellar accretion heating affects the fragmentation fraction. We find that discs fragment only if they become sufficiently large early in their lives. This size depends sensitively on where mass is added to the discs during the collapse of their parent molecular cloud core. By choosing intermediate infall locations leading to a synthetic disc size distribution that is in agreement with the observed one, we find a fragmentation fraction between 0.1 and 11 %, depending on the efficiency of stellar accretion heating of the discs. We conclude that the early disc size is mainly determined by the infall location during the collapse of the molecular cloud core and controls the population-wide frequency of fragmentation. Stellar accretion heating plays an important role for fragmentation and must be studied further. Our work is an observationally-informed step towards a prediction of the frequency of giant planet formation by gravitational instability. Upcoming observations and theoretical studies will deepen our understanding of the formation and early evolution of discs, eventually allowing to understand how infall, disc morphology, giant planet formation via gravitational instability, and the observed extrasolar planet population are linked.


Introduction
Protoplanetary discs are of central importance in astrophysics: they are the birth places of planets. Two main pathways for planet formation exist: core accretion (Safronov 1972;Pollack et al. 1996;Ida & Lin 2004;Mordasini et al. 2012) and gravitational instability (GI, Kuiper 1951;Cameron 1978;Boss 1997;Kratter & Lodato 2016). GI remains the leading formation mechanism for some observed giant planets on wide orbits (Marois et al. 2008(Marois et al. , 2010Teague et al. 2018) and may explain the formation of giants around low-mass stars (Morales et al. 2019). GI has also been proposed for the formation of intermediatemass planets (Deng et al. 2021) and very young giants (Currie et al. 2022;Cadman et al. 2021). Fragmentation, the collapse of parts of a protoplanetary disc to form a bound object, is a necessary condition for GI. It is still debated Whether a bound clump formed by GI can survive to become a giant planet. The interplay among several physical processes like migration, accretion of gas and solids, dynamical collapse, grain sedimentation and core formation needs to be studied in order to answer this question (Bodenheimer et al. 1980;Helled & Schubert 2008;Forgan & Rice 2013;Nayakshin 2017;Forgan & Rice 2013). While many questions around GI are still open, the first step of fragmentation has been studied extensively using 2D and 3D hydrodynamic simulations (Toomre 1964;Gammie 2001;Mayer et al. 2002;Vorobyov & Basu 2010;Baehr et al. 2017;Deng et al. 2017;Jin et al. 2020). Discs are found to be prone to fragmentation in the outer regions (tens of au) if they are cold and massive enough (Boss 1998(Boss , 2000Gammie 2001;Clarke 2009;Rice et al. 2004Rice et al. , 2005. However, while many discs are expected to undergo a phase of gravitational instability early in their evolution, (Durisen et al. 2006) it is uncertain how many of them fragment. Fragmentation could also be triggered externally by binary companion in discs that are not massive enough to fragment on their own (Cadman et al. 2022). Gravitational instability also drives transport of angular momentum via spiral arms and discs can remain in a quasi-steady state of self-regulation Cossins et al. 2009;Rice 2016). A better understanding of the frequency of fragmentation is important for population synthesis studies of GI as it determines the expected number of clumps. The number of clumps is in turn an upper bound for the number of planets that may from via gravitational instability. A theoretical prediction of the fragmentation fraction is thus a necessary step towards a prediction of the frequency of planet formation via GI. This is of high interest for comparisons with the extrasolar planet population.
Even if GI does not produce a significant number of planets, the formed clumps may still be responsible for free-floating planetary mass objects (Forgan et al. 2018) or episodic accretion (Audard et al. 2014).
In previous work (Schib et al. 2021, hereafter Paper I) we performed a population synthesis of protoplanetary discs with a focus on fragmentation. We included the formation of the starand-disc system by infall from the molecular cloud core (MCC), since fragmentation is most likely to occur during or shortly after disc formation where the disc-to-star mass ratio is highest. It was found that fragmentation depends sensitively on the disc formation process. If young discs are compact, as expected by some non-ideal magnetohydrodynamic (MHD) simulations of cloud collapse (Hennebelle et al. 2016), fragmentation is suppressed completely. If discs are large already early on, as expected from radiation hydrodynamic simulations (Bate 2018), fragmentation may be prevalent. Our comparison with observed sizes of Class 0 discs suggested that real disc sizes lie somewhere between the two extremes. This raises the question of the likelihood of discs to fragment if they have sizes consistent with observations. While we studied a broad parameter space in Paper I, in this paper we focus on a specific aspect: how the disc size distribution from infall sets the fraction of self-gravitationally fragmenting discs.
First, we briefly review the model. Then we describe our analysis and present the results. A short discussion is followed by the summary and conclusions.

Model and investigated parameter space
Our model is described in detail in Paper I. We apply it without changes. The most important aspects of the model are summarised here for convenience.
The model consists of a 1D, vertically integrated gas disc around a single star. We use cylindrical polar coordinates, with r denoting the radial direction. The disc's mid-plane is located at z = 0 au. The surface density Σ evolves in time t according to the viscous evolution equation (Lüst 1952 andLynden-Bell &Pringle 1974): S ≡ S (r, t) is a source/sink term: The summands describe, from left to right, the infall from the MCC, internal (Clarke et al. 2001) and external (Matsuyama et al. 2003) photoevaporation, and the removal of mass due to fragmentation. Detailed descriptions and analytical expressions are given in Paper I. We assume a short phase of infall at the beginning of our simulations. The precise implications of this assumption are difficult to assess. We discuss this further in Sect. 4. The constant infall rate is chosen based on a selection of systems from the radiation hydrodynamic simulation of star formation in Bate (2018) (B18). Our calculation of the viscosity is based on the αprescription (Shakura & Sunyaev 1973): ν = αc s 2 /Ω, with c s the isothermal sound speed and Ω the disc's angular frequency. We use different values for α to describe different physical processes. During the infall phase, the disc-to-star mass ratio can reach values of 0.5 or higher, and the transport of angular momentum is dominated by gravitational torques (Shu et al. 1990;Harsono et al. 2011). We use the prescription given in Eq. (32) of Kratter et al. (2010) to account for these effects. After the infall phase, the transport of angular momentum can remain very high as long as the disc is self-gravitating (Kratter et al. 2008). We use the corresponding parametrisation for the transport through spiral arms given in Zhu et al. (2010) in this case. In the absence of gravitational torques we apply a minimal background value of 0.01 (Kimura et al. 2016) that describes any other physical process responsible for the transport of angular momentum, like the magneto-rotational instability (MRI, Balbus & Hawley 1991), hydrodynamic instabilities or MHD winds (Turner et al. 2014). We note that MRI is no longer considered the main source of turbulence in protoplanetary discs due to their weak ionisation (e.g. Lesur et al. 2022;Cui & Bai 2021). If the long-term evolution is indeed better described by a value of α between 10 −4 and 10 −3 , our model would predict disc lifetimes much longer than observed, as discussed in Paper I, unless photoevaporation is much stronger than we assume.
The disc's midplane temperature is calculated assuming an energy balance on the surface. The following processes are considered: viscous heating, shock heating from the infall, irradiation from a constant background and irradiation by the star. The irradiation temperature is calculated as (Hueso & Guillot 2005): In Eq. 3, T * and R * are the star's temperature and radius, respectively. We set d ln(H)/d ln(r) = 9/7 for numerical reasons (Hueso & Guillot 2005). The star's intrinsic temperature and radius are interpolated from stellar evolution tables (Yorke & Bodenheimer 2008). As the disc evolves, some of its mass crosses the inner truncation radius assumed to be constant at 0.05 au. It is accreted onto the central star. For the majority of the simulations, we assume that half of the gravitational binding energy is radiated away, thermalised and heating the disc. For the total stellar luminosity L * we write (Baraffe et al. 2009;Vorobyov & Elbakyan 2018;Oliva & Kuiper 2020): where L int is the star's intrinsic luminosity and the second term is half the stellar accretion (shock) luminosity, with f acc an efficiency factor, G the gravitational constant, M * the stellar mass andṀ * the accretion rate of disc material onto the star. The effective stellar temperature in Eq. 3 is then calculated as: with σ the Stefan-Boltzmann constant. Nominally, we assume, f acc =1 ( Baraffe et al. 2009;Kimura et al. 2016;Oliva & Kuiper 2020;Kuffmeier et al. 2018). Treating stellar irradiation this way is a strong assumption and likely an overestimate, since the accretion luminosity is absorbed/thermalised in the infalling material. We discuss this further in Sect. 4. There is also some shock heating related to material from the MCC as it reaches the disc. Paper I suggest that its effect is small compared to the accretion on the star and we include it in all runs.
During the disc's evolution, the criteria for fragmentation are checked in every time step. We use the same criteria as in Paper I.
The main (necessary) condition is the Toomre criterion (Toomre 1964): where Q Toomre is the Toomre parameter and κ the epicyclic frequency. If the disc is supplied by infalling material faster than it can be transported away, it will invariably fragment when Q Toomre drops below unity (Boley 2009). This is called the "infall dominated regime". The condition for this regime is given in Paper I (Eq. (25). If the infall rates are not high enough, the disc transitions into the "cooling dominated regime". In that case it can only fragment with sufficient cooling, which can be stated as t cool < β c (Gammie 2001), where we set the critical cooling β c to 3 (Deng et al. 2017). The cooling timescale t cool is given in Eq. (27) in Paper I. When the disc fragments, a bound clump composed predominantly of gas forms. Such a clump may undergo a number of processes like migration, accretion and disruption. The combined outcome of these is uncertain. In this work we avoid this complexity and follow the approach described in Sect. 2 of Paper I: we assume the clump migrates to the inner disc quickly and accretes on the central star. The initial fragment mass of (Boley et al. 2010, tough see also Forgan & Rice 2011) is removed from the disc an added to the star. In Paper I, we performed a disc population synthesis with five runs, each of them consisting of 10 000 simulations. We varied several parameters to investigate how they affect the masses, sizes and fragmentation of the discs. The run "hydro" (the baseline case) was set up with initial conditions extracted from results presented in the hydrodynamic simulation of star formation by Bate (2018). This includes the infall location, which was chosen to approximately have the same angular momentum of the discs compared to Bate (2018). Based on the "hydro" run, we performed three additional runs where we varied the background viscosity, the accretion heating efficiency and the initial fragment mass. We also performed a run where we chose the infall location differently. The run "MHD" is described in more detail below. The simulations we conducted in the present work are executed in an analogous fashion as the ones in Paper I: Every simulation starts with a very low stellar mass and disc mass (∼ 0.01 M ), is supplied with infalling material for a given length of time and then evolves further until the disc dissolves. First, we investigated the dependency of the fragmentation fraction on the the early disc size distribution. For this, we performed a number of disc population syntheses where we vary the infall radii. A prediction for the early disc size is given in Hennebelle et al. (2016). The authors perform 3D non-ideal MHD simulations and show that the disc radii at early stages of disc formation agree within a factor of two with their analytic expression: In Eq. 8, A is a measure of the ambipolar diffusivity, B z denotes the magnetic field in the inner part of the core and M * d the combined mass of star and disc. In Paper I we performed a run "MHD" in which the infalling material was deposited at an infall radius constant in time that lies close to the star (a few au, where the specific value depends on stellar mass). We demonstrate that the disc radii at the end of the infall phase agree well with Eq. 8 for A = 0.1 s and B z = 0.1 G. We find that none of these discs fragmented due to their small sizes. In Paper I we Run Infall radius Stellar accretion heating efficiency f acc 1 "5× MHD" 1 2 "8× MHD" 1 3 "12× MHD" 1 4 OBS IRR 1 5 OBS NOIRR 0 6 OBS REDIRR 1/12 Table 1: Overview of the runs. The infall radius is the location at which the infalling material is deposited in the disc (see Sect. 2). In runs OBS IRR, OBS NOIRR and OBS REDIRR, the infall radii were chosen to match observed disc sizes. In run OBS REDIRR the accretion heating was reduced by a factor of twelve relative to runs 1-4. also performed a "hydro" run, where the infall radii were taken from the (non-magnetized) radiation-hydrodynamic simulations of B18. This run produced discs of a much larger size of which a significant fraction (22%) fragment 1 . Interestingly, the resulting synthetic disc radius distribution in these two runs brackets the observed disc radius distribution (Tobin et al. 2020).
For the present work we performed in total six new runs. In the first three new runs ("5x MHD", "8X MHD" and "12x MHD") we subsequently increased the infall radii. We did this by multiplying the (small) infall radii from the "MHD" run from Paper I by a constant factor that differs in each of the new runs: 5, 8 and 12. A different value for the infall radius was used in the "MHD" run in each mass bin in order for the early disc radius to agree with Eq. 8. In all but the least massive system, this radius is close to 2 au. Therefore, typical values for the new runs are 10, 16 and 24 au.
In a second step, we modified the infall radii more finely in order to reproduce observed radii of Class 0 discs (run OBS IRR). All these runs assume a full efficiency of stellar accretion heating ( f acc = 1). In a last step, we studied the effect of accretion heating. We performed two additional runs ("OBS NOIRR" and "OBS REDIRR" ) with inactive ( f acc =0) or reduced accretion heating ( f acc =1/12). Table 1 gives an overview of the six runs performed in this work. As in Paper I, we chose the duration of the infall (the length of the constant infall phase) in such a way, that the distribution of stellar masses at the end of the simulations agrees with the observed initial mass function (IMF) (Chabrier 2005). A figure demonstrating this can be found in App. A. It shows the stellar number density in pc −3 per logarithmic interval of mass for all the runs performed in this work. The agreement is reasonable for all runs across the parameter space from 0.05 M to 5 M we studied. In order to compensate for the decreasing number density towards high masses, we divided this parameter space into 100 logarithmically spaced mass bins and set up our simulations in such a way that we obtain approximately 100 systems with a final stellar mass lying in each bin. This ensures our statistics are reasonably robust across the entire mass range. O. Schib et al.: The link between infall location, early disc size, and the fraction of self-gravitationally fragmenting discs  Table 2: An overview of the results. Listed are the disc's radii and masses at the end of the infall phase, the fraction of discs that fragment, and the reduced disc lifetimes (see Sect. 3 and Appendix A for further details).

Results
In this section, we present the results of our simulations. First, we discuss how different infall radii influence the early disc sizes, then we concentrate on the influence of stellar accretion heating on discs with similar early disc sizes. An overview of the results for all six runs is given in Table 2. The numbers shown represent weighted mean values calculated assuming that the final distribution of stellar masses agrees with the IMF. The disc sizes are increasing from run "5x MHD" to "12x MHD" as expected, while they are the same for the remaining runs as further discussed in Sect. 3.1. The fragmentation fractions change significantly for different runs as discussed in Sect. 3.1 to 3.3. On the contrary, disc masses and lifetimes are very similar. We discuss these in Appendix A.
3.1. Infall radii and early disc radii Figure 1 depicts the disc radii of the first three runs, measured at the end of the infall phase, as cumulative distributions. The disc radius is defined as the radius containing 63.2 % of the disc's mass (B18). The figure also shows the observational Result from Tobin et al. (2020). The authors perform a multi-wavelength survey of hundreds of protostars and use dust continuum emission to measure Class 0 dust disc radii. During this early phase the gas disc radii should not yet differ very much from those obtained from continuum emission (Birnstiel & Andrews 2014, B18). Nevertheless, this is a clear limitation of our study and we discuss this further in Sect. 4. Fig. 1 shows that runs 1 to 3 cover a region of parameter space much closer to observed disc radii than the "hydro" and "MHD" runs we performed in Paper I. It is therefore interesting to check which of these systems fragment. The region in a disc prone to fragmentation is typically outside of ∼ 50 au. There, run "5x MHD" has radii at the lower end of the observed population, run "8x MHD" is roughly in agreement, and run "12x MHD" exhibits larger radii. This increase in disc sizes is reflected in an increasing fragmentation fraction: We find that in run "5x MHD", none of the systems fragment. For run "8x MHD" we find 0.19 % of systems fragment, while in run "12x MHD" the fraction of fragmenting systems is 3.5 %. These values are computed by weighting the different mass bins such that the final stellar masses are consistent with the IMF.
Taken at face value, this puts the observed disc radii right at the edge of fragmentation, with only about 2 out of 1000 discs fragmenting. This would imply that a necessary condition for giant planet formation via gravitational instability is only very rarely met (but also not never). However, the observed distribu-Tobin + '20 "MHD" "5x MHD" "8x MHD" "12x MHD" "hydro"  Fig. 1: Cumulative distributions of the disc radii at the end of the infall phase for runs "5x MHD", "8x MHD" and "12x MHD". The runs "MHD" and "hydro" from Paper I are shown for reference. The blue shaded region shows the sizes of observed Class 0 discs for systems in isolation from Tobin et al. (2020). tion of radii has a different, broader shape in comparison to our syntheses.

Agreement with observed radii
The different shapes of the size distributions discussed raise the question whether a synthetic population of discs with a distribution of radii in (better) agreement with observations can also be simulated. This is not straight-forward, since there is no simple relationship between the infall radius and the early disc radius. The disc's radius changes with time through the evolution of the disc and this effect depends on the disc's mass. In Sect. 2 we describe how our simulations are distributed in different mass bins. More massive systems tend to produce larger discs in our model. Such a correlation is also seen in the hydrodynamic simulations of Bate (2018) as well as in observed discs (Tobin et al. 2020). It is also expected from Eq. 8. Increasing the infall radii in the bins corresponding to higher masses would therefore lead to a broader distribution of radii. Of course, for such a modified set of initial conditions, the infall times need to be adapted again in order to satisfy the requirement to fit the IMF (Sect. 2) which would affect the distribution of radii. After iterations, we infer the distribution of infall radii that satisfies both constraints at the same time (run OBS IRR). The resulting distribution of radii is shown in Fig. 2 (orange solid line). The fraction of discs fragmenting in run OBS IRR is 0.12 %, similar to run "8x MHD". This is expected as these two runs have a very similar distribution of radii above 80 au, where fragmentation is most likely. If this result is robust, it means fragmentation is very rare indeed, independently of subtleties in the disc radius distribution.

Irradiation from stellar accretion heating
Next we revisit our assumption concerning stellar accretion (shock) heating. If this effect is indeed overestimated by our treatment as mentionned in Sect. 2 and discussed below, this would lead to too high temperatures in the outer discs and inhibit fragmentation through Eq. 6. We tested this hypothesis by setting as a limiting case the second summand in Eq. 4 to zero ( f acc =0), effectively reducing the irradiation to the star's intrin- sic luminosity. We then constructed a new set of initial conditions that again satisfies the constraints in stellar mass and disc radius: run OBS NOIRR. This run is also shown in Fig. 2. In run OBS NOIRR, a high fraction of 25 % of the systems fragment, demonstrating the strong influence of accretion heating on fragmentation besides disc size. The assumptions about accretion heating made in runs OBS IRR ( f acc =1) and OBS NOIRR ( f acc =0) are likely extremes. The former case appears, however, clearly more likely than the latter, because there are direct observations (Cieza et al. 2016) that the high stellar luminosity during episodes of strong accretion shifts the water iceline to large orbital distance, implying that f acc must be larger than zero. We discuss this in more detail in Sect. 4. Fig. 3 shows the stellar luminosities at the end of the infall phase together with the observed luminosities of isolated Class 0 protostars (Tobin et al. 2020). The figure shows first that a very large range in luminosities is covered. Second, we note that the synthetic luminosity distributions when including only the intrinsic luminosity ( f acc =0) or when considering the the sum of intrinsic and full accretion luminosity f acc =1) are bracketing the observed luminosity distribution. This suggest to consider an intermediate value for f acc . In run OBS REDIRR we thus assess a  scenario with such reduced accretion heating. To roughly reproduce the observed luminosity distribution, we find that the accretion term (second summand in Eq. 4) needs to be reduced by about a factor of 12, i.e., f acc =1/12. Again we chose infall radii and infall times in such a way that the constraints from both disc sizes and final stellar mass distribution are satisfied. The fraction of systems that fragment is 11 % in run OBS REDIRR. This fragmentation fraction lies between what we found for runs OBS IRR ( f acc =1, 0.12% fragmentation) and OBS NOIRR ( f acc =0, 25% fragmentation) as expected: fragmentation is reduced less than in run OBS IRR because the discs are not as strongly heated and thus stabilised. The number is also interesting in another context: the frequency of substellar companions on wide orbits, which might have formed via gravitational instability, is on the order of a few to a few ten percent (Vigan et al. 2021). Even though our approach to fit the luminosity distribution is simplistic, run OBS REDIRR is still our most realistic, observationally-informed scenario. It satisfies the constraints from observed stellar masses and disc sizes while also exhibiting luminosities in the same range as is observed. We discuss some limitations of our study in the following section. In Appendix A we show additional results for runs OBS IRR, OBS NOIRR and OBS REDIRR: early disc masses, final stellar masses and disc lifetimes of these systems. The infall radii used in these runs are given in Appendix B, the infall times in Appendix C.

Fragmentation
The likelihood of the disc to fragment and the location of the instability depend on the evolution of the disc's surface density and temperature. Table 3 gives an overview of the most important properties for our runs OBS NOIRR and OBS REDIRR. Fig. 4 shows the fraction of fragmenting discs as a function of final stellar mass as well as the mean values of the number, location and mass of the fragments. The top left panel reveals that fragmentation only becomes important for systems with a final stellar mass 0.3 M , with typical values of ∼ 40 % in run OBS REDIRR and ∼ 90 % in run OBSNOIRR. The fraction decreases slightly for final stellar masses above 3 M . This was not seen in Paper I and is a consequence of the different choice of infall radii. Interestingly, a similar trend is found in radial velocity surveys. An increase in giant planet occurrence rate with host star mass is found up to a host star mass of 1.9 M with a decline at higher host masses (e.g. Reffert et al. 2015). It is still unclear whether the trend in fragmentation fractions persists in fully formed giant planets. Also, it is still unknown whether a similar decrease at higher stellar masses exists in directly imaged planets due to insufficient detections at higher stellar masses (Vigan et al. 2021, see also Janson et al. 2021;Wagner et al. 2022). The top right panel of Figl 4 shows that the number of fragments de-pends strongly on irradiation. While the runs with accretion heating of the outer disc (OBS IRR and OBS REDIRR) typically have around ten fragments, this number is larger without accretion heating and can lead to more than a hundred fragments. The number of fragments in this case is lower than what we found in Paper I, less than half in the OBS IRR and OBS REDIRR runs compared to the "hydro" run. This is a consequence of the more compact discs.
The location where the discs fragment (bottom left panel in Fig. 4) depends mainly on the accretion heating: The hotter the outer discs are, the further out they fragment, in agreement with Paper I.
The initial fragment mass (bottom right panel) is influenced both by the location of fragmentation and the accretion heating. More accretion heating means higher temperatures in the outer disc which leads to more massive fragments through Eq. 7 directly. It also does so indirectly by moving the location of fragmentation further out. This is in good agreement with what was found in Paper I. Run OBS REDIRR features quite low initial fragment masses, typically 0.5 M J to 1 M J .

Discussion
In our analysis we had to make several strong assumptions due to theoretical and observational uncertainties. Here we discuss how these may influence our results.
The most important assumption is the comparability of gas and dust disc radii. For older discs like Class II this would certainly not be justified, as dust would have had time to grow and drift (Birnstiel et al. 2010;Testi et al. 2014;Tobin et al. 2020). Since we compare the sizes of young discs (Class 0), the differences are expected to be smaller. However, it is still possible that there is a mismatch between observed Class 0 dust discs sizes and the sizes of the corresponding gas discs in both directions.
The definition of the disc radius used in this work (63.2 % of the mass) is identical to that used in Bate (2018). It is well suited for comparisons with observations of discs with truncated power-law surface density profiles Σ ∝ r −γ independent of the index γ (as long as γ < 2, see the discussion in Sect. 2.3.1 in Bate 2018). For our comparison we evaluate the disc radii at a welldefined point in time: the end of the infall phase. If this coincides with the end of the Class 0 phase, this means we overestimate the disc sizes, since the observed discs are expected to to represent an (unknown) distribution of Class 0 ages. If on the contrary, the Class 0 phase lasts longer than our infall phase and the discs would have more time to spread, or if the size of the observed gas disc is smaller due to growth/drift of dust particles this would make our discs too small. Based on our discussion in Sect. 3, an overestimate of the disc sizes would lead to an overestimate of the fraction of discs that fragment. Conversely if our discs are too small, we are likely to underestimate the occurrence rate of fragmentation.
Our results indicate that only a minority of discs undergo fragmentation. This does not mean, however, that gravitational instabilities are not important in the non-fragmenting bulk of the population. In fact, most of the discs in run OBS REDIRR have Q Toomre < 2 somewhere, at some point during their evolution. This suggests that these discs were self-gravitating for at least a short period during their early lives.
In comparing the size distribution of our full sample with the sample from Tobin et al. (2020) we also implicitly assume their sample will be representative of the IMF. Unfortunately, it is very difficult to measure stellar masses of such young systems, and the precise distribution of stellar masses of the ob-served sample is unknown. The observed luminosities cannot be translated into masses directly due to accretion. Our comparison of luminosities with and without accretion (Fig. 3) demonstrate the potential extent of this discrepancy: the total luminosities in our run OBS IRR are higher by more than an order of magnitude compared to the observations, which are again much higher than our intrinsic luminosities. This could be explained in part by a different distribution of masses in the observed sample. Our populations have a stellar mass distribution as in the Chabrier (2005)-IMF. However, a different IMF dominated may be more adequate for the Orion region (Da Rio et al. 2012, see also Luhman 2000). The difference could also be caused by the way the stars accrete. The difference between observed and theoretically expected luminosities has long been discussed and is historically known as the "protostellar luminosity problem" (Kenyon et al. 1990). This apparent problem has been solved by demonstrating that both episodic accretion and a longer accretion phase can produce luminosities that are in line with observations (e.g. Offner & McKee 2011;Dunham & Vorobyov 2012).
The precise nature in which accretion happens is still not fully understood. Also, it is still unknown which of the two proposed mechanisms described above contributes to what degree to the explanation of the "protostellar luminosity spread" (Sect. 3 in Fischer et al. 2022). It is therefore unclear, how our assumption of a phase of constant accretion should be modified in order to account for the observed luminosity spread. This is an important topic that needs to be studied further. The exact nature of accretion affects the stellar evolution, and hence disc fragmentation. A study that assesses the interplay of episodic accretion, and stellar evolution in a 3D MHD simulation is presented by Kuffmeier et al. (2018).
Finally, we also assumed that once the disc fragments, a clump of a given mass M F (Eq. 7) forms, migrates inwards quickly and accretes on the star. While this work does not focus on the fate of the clumps, it should be noted that our simplifying assumptions could influence the disc's evolution and its fragmentation at later stages. While M F is the initial clump's mass, it could be significantly different from the mass accreted onto the star and/or the final mass of the surviving clump. The clump's mass is expected to decrease via mass loss (e.g., tidal disruption in the inner disc, Boley et al. 2010;Nayakshin 2010). The forming clump could also grow in mass via further accretion of disc gas bringing it into the Brown Dwarf or stellar mass regime. If gas accretion is rapid, a gap may form which would slow down migration (D'Angelo & Lubow 2008;Rowther & Meru 2020;Zhu et al. 2012;Oliva & Kuiper 2020;Schib et al. 2022). A deep gap could reduce gas flow to the inner disc and make the outer disc prone to further fragmentation.

Summary, conclusions, and outlook
We performed a population synthesis of protoplanetary discs focusing on disc sizes and the fragmentation likelihood. In three sets of initial simulations (runs "5x MHD", "8x MHD" and "12x MHD") we demonstrated how increasing the infall radii controls the fragmentation of protoplanetary discs. We then constructed initial conditions (infall locations) in such a way that the synthetic disc radius distribution at the end of the infall phase agrees with the observed radius distribution of Class 0 discs, while the distribution of stellar masses at the end of our simulations agrees with the IMF (runs OBS IRR, OBS NOIRR and OBS REDIRR). For these three runs, we also investigated how fragmentation is influenced by accretion heating. Our most im- -Protoplanetary discs need to reach a certain size while they are massive in order to fragment. Observed discs appear to have radii just large enough to fragment in a few cases. -The early disc size distribution is mainly determined by the infall location. To get a synthetic size distribution that is in agreement with the observed one, we find infall locations that are between those predicted by some magnetized and pure hydrodynamic star formation simulations. Specifically, the infall locations need to be increased by factors 5-8 relative to the MHD simulations of Hennebelle et al. (2016) which corresponds to a reduction by factors 2-3 relative to pure HD simulations B18. -The fraction of discs that actually fragments does, however, not only depend on their early sizes, but also on the fraction f acc with which stellar accretional (shock) luminosity contributes to heating the outer disc. -In disc population syntheses fitting the observed Class 0 disc radius distribution, we find fragmentation fractions of 0.12% and 25% for the limiting cases of f acc =1 and 0, respectively. Here, the former scenario is more likely than the latter, since heating of discs by stellar accretion luminosity is an observed phenomenon Cieza et al. (2016). -Finally, in our most realistic case, we empirically determine f acc ≈1/12 by requiring that our synthetic luminosity distribution is in approximate agreement with the observed luminosity distribution of Class 0 protostars (while still also fit-ting the observed disc radius distribution). In this case, the fragmentation fraction is 11%, which is our observationallyinformed best estimate of the fraction. -A better understanding of accretion is crucial for the understanding of the evolution of protoplanetary discs and the formation of planets. It needs to be investigated further using simulations of discs with an evolving dust component, including radiative transport, and evolution models of accreting stars.
The discovery of AB Aurigae b (Currie et al. 2022) as a massive, still actively accreting giant planet at a large orbital distance of 93 au has recently given observational support to gravitational instability as a giant planet formation mechanism. However, it is still unknown at which frequency GI contributes to giant planet formation. Qualitatively, it is often assumed to be a rather rare process (Vigan et al. 2021), but quantitative theoretical predictions are currently sparse.
A necessary condition for giant planet formation is disc fragmentation. The fragmentation fractions of about 0.1 to 11% that we find in this work are thus an upper bound for the frequency of giant planet formation via GI. Our work represents a step towards a quantitative theoretical prediction of the importance of GI as a giant planet formation mechanism.
Future work promises to help in our understanding of the formation and evolution of stars, discs, and planets. Simulations of star formation are becoming more and more detailed, allowing for a more precise treatment of the interaction between stars and discs Kuiper & Hosokawa 2018). The study of disc kinematics continues to improve our understanding of the mass distribution and motion of the gas in the outer discs and the presence of structures and embedded protoplanets (Teague et al. 2019;Pinte et al. 2020;Terry et al. 2022). Observations of the gas emission of discs may help constraining young disc sizes in the near future (e.g. Rota et al. 2022). Finally, exoplanet surveys continue to improve our understanding of the demographics of planets, including in particular distant and forming planets which are prime candidates for a formation via GI (Marois et al. 2008;Chauvin et al. 2017;Currie et al. 2022). This will eventually allow to understand how star formation, infall, disc morphology, giant planet formation via gravitational instability, and the observed extrasolar planet population are connected.

Appendix A: Properties of runs OBS IRR, OBS NOIRR and OBS REDIRR
Here we show some additional results for runsOBS IRR, OBS NOIRR and OBS REDIRR. This can be useful when comparing these runs to the ones from Paper I or other studies. An overview of these results is given in Table 2. Fig. A.1 depicts the distributions of stellar masses at the end of the simulations (i.e. when the disc is gone). It shows that the agreement with the IMF is reasonable for all runs. In Fig. A.2 we show the distribution of disc masses at the end of the infall phase. Also shown is the observational result from Tychoniec et al. (2018). The masses are similar for runs OBS IRR, OBS NOIRR and OBS REDIRR. They are comparable, tough somewhat larger than the observed Class 0 masses. Systems with less accretion heating are slightly less massive. This result may seem counter-intuitive, since lower disc temperatures should lead to less accretion to the star and therfore more massive discs. However the infall radii in the less massive (and more abundant) systems are lower for systems with less accretion heating. herefore more gas is accreted on the star early, which explains the lower masses. Figure A.3 depicts the disc fractions based on the reduced near-infrared lifetimes t life . It shows excellent agreement with the disc fractions based on the fits from Richert et al. (2018), when the "magnetic" pre main  sequence (PMS) model from Feiden (2016) is used to determine cluster ages. The reduction is based on the beginning of PMS (Kimura et al. 2016). We note that there is substantial uncertainty involved when comparing disc lifetimes from simulations to observed disc fractions. A detailed discussion of this is found in Sect. 6.3 of Paper I. The same is done for the initial values for the stellar mass, disc mass and infall rate. In contrast to the infall times, the distributions for the other initial values are identical to what is used in Paper I and are given there (Appendix E).