Planet formation occurs over a few Myr within protoplanetary disks of dust and gas which are often assumed to evolve in isolation. However extended gaseous structures have been uncovered around many protoplanetary disks suggestive of latestage infall from the interstellar medium ISM. To quantify the prevalence of latestage infall we apply an excursion set formalism to track the local density and relative velocity of the ISM over the disk lifetime. We then combine the theoretical BondiHoyleLyttleton BHL accretion rate with a simple disk evolution model anchoring stellar accretion timescales to observational constraints. Disk lifetimes masses stellar accretion rates and gaseous outer radii as a function of stellar mass and age are remarkably well reproduced by our simple model that includes only ISM accretion. We estimate that 2070 of disks may be mostly composed of material accreted in the most recent half of their lifetime suggesting that disk properties are not a direct test of isolated evolution models. Our calculations indicate that BHL accretion can also supply sufficient energy to drive turbulence in the outer regions of protoplanetary disks with viscous  SUBSSSUB  10SUP5SUP to 10SUP1SUP although we emphasize that angular momentum transport and particularly accretion onto the star may still be driven by internal processes. Our simple approach can be easily applied to semianalytic models. Our results represent a compelling case for regulation of planet formation by largescale turbulence with broad consequences for planet formation theory. This possibility urgently motivates deep observational surveys to confirm or refute our findings.
['Protoplanetary discs', 'Star forming regions', 'Planet formation', '1300', '1565', '1241', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
Planet Formation Regulated by Galacticscale Interstellar Turbulence
{'Planet formation regulated by galactic-scale interstellar turbulence': "<!-- image -->  \n1 Universit'e Cˆote d'Azur, Observatoire de la Cˆote d'Azur, CNRS, Laboratoire Lagrange, 06300 Nice, France 2 Universit'e Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France 3 Center for Astrophysics - Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA", 'ABSTRACT': 'Planet formation occurs over a few Myr within protoplanetary discs of dust and gas, which are often assumed to evolve in isolation. However, extended gaseous structures have been uncovered around many protoplanetary discs, suggestive of late-stage in-fall from the interstellar medium (ISM). To quantify the prevalence of late-stage in-fall, we apply an excursion set formalism to track the local density and relative velocity of the ISM over the disc lifetime. We then combine the theoretical BondiHoyle-Lyttleton (BHL) accretion rate with a simple disc evolution model, anchoring stellar accretion time-scales to observational constraints. Disc lifetimes, masses, stellar accretion rates and gaseous outer radii as a function of stellar mass and age are remarkably well-reproduced by our simple model that includes only ISM accretion. We estimate 20 -70 percent of discs may be mostly composed of material accreted in the most recent half of their lifetime, suggesting disc properties are not a direct test of isolated evolution models. Our calculations indicate that BHL accretion can also supply sufficient energy to drive turbulence in the outer regions of protoplanetary discs with viscous α SS ∼ 10 -5 -10 -1 , although we emphasise that angular momentum transport and particularly accretion onto the star may still be driven by internal processes. Our simple approach can be easily applied to semi-analytic models. Our results represent a compelling case for regulation of planet formation by large-scale turbulence, with broad consequences for planet formation theory. This possibility urgently motivates deep observational surveys to confirm or refute our findings.  \nKeywords: protoplanetary discs - planet formation - star forming regions', '1. INTRODUCTION': "Planets form and begin their evolution within discs of gas and dust that persist for a few Myr around young stars (e.g. Haisch et al. 2001). The properties of these protoplanetary discs, and the physics that shapes them, is critical to our understanding of the origin of the (exo)planets we observe today. For many years, models for planet formation have been based on assumptions that the star-disc system by-and-large evolves in isolation from its surroundings. This assumption is now being challenged on a variety of fronts. External mechanisms that shape protoplanetary discs include irradiation from nearby massive stars (Winter & Haworth 2022, and references therein), neighbour star-disc encounters (Cuello et al. 2023, and references therein) and the late stage in-fall of gas from the surroundings (e.g. Dullemond et al. 2019; Kuffmeier et al. 2020, 2021, 2023). It is this last influence that is the topic of this work.  \nRecent observations in molecular line emission (Huang et al. 2021) and infrared scattered light (Garufi et al. 2020; Benisty et al. 2023) are frequently revealing extended structures associated with protoplanetary discs, on spatial scales several times the disc outer radius. In-fall has been suggested as a possible origin for these structures (Dullemond et al. 2019; Kuffmeier et al. 2020, 2023; Hanawa et al. 2024), which are highly complex, including large scale spirals and streamers around T Tauri (Huang et al. 2023) and Herbig AeBe (Boccaletti et al. 2020) stars, often associated with reflection nebulosity (Gupta et al. 2023), accretion outbursts (P'erez et al. 2020; Hales et al. 2020; Dong et al. 2022) and misaligned inner discs (Ginski et al. 2021). Evidence that is suggestive of in-fall events has been found not only during the earliest stages of disc evolution, when the star is still embedded in its natal environment (Jørgensen et al. 2009; Maury et al. 2019; Codella et al. 2024), but also around stars such as AB Aur (Nakajima & Golimowski 1995; Grady et al. 1999; Fukagawa et al. 2004) DR Tau (Mesa et al. 2022) and SU Aur (Ginski et al. 2021), which are  \nall mature star-disc systems of age ∼ 1 -3 Myr. A recent estimate from a sample size of 43 disc-hosting stars suggests approximately 16 percent are interacting with the ambient ISM (Garufi et al. 2024). This material should fall onto the disc on a time-scale much shorter than the disc lifetime (Dullemond et al. 2019; Kuffmeier et al. 2020), suggesting that either these events occur much more frequently than observed or the external structure is replenished by the ISM for a prolonged period. Indeed, substantial mass accumulation at least over the first ∼ 1 Myr of the disc life-time is hinted at by the fact that discs younger than 1 Myr in ρ Ophiucus and Corona Australis are somewhat more compact and lower mass than observed in ∼ 1 -2 Myr old regions (Testi et al. 2016, 2022; Cieza et al. 2019; Cazzoletti et al. 2019; Williams et al. 2019). It is possible that discs contain a substantial fraction of recently accreted material for much of their lifetime, which would have significant consequences for disc evolution models.  \nIn-falling gas may be the consequence of Bondi-Hoyle-Lyttleton (BHL) accretion, where low velocity gas is captured due to the gravitational potential of the star (Bondi 1952; Edgar 2004). Theoretically, the rate of capture of material from the interstellar medium (ISM) via BHL accretion can be substantial while stars remain bound to a dense gas reservoir. Padoan et al. (2005) proposed this mechanism to explain the observed relationship between stellar accretion rates and stellar mass. This idea was supported by subsequent work (Throop & Bally 2008; Klessen & Hennebelle 2010), while more recent analysis of turbulent hydrodynamic simulations in a periodic box have shown that late stage accretion may be significant source of mass and angular momentum for ∼ 1 Myr old discs in high gas density environments (Kuffmeier et al. 2023; Pelkonen et al. 2024; Padoan et al. 2024). In-fall has also been suggested as a mechanism to explain examples of substantial protoplanetary discs persisting around old stars (Scicluna et al. 2014; Derkink et al. 2024). However, assessing the importance of this process for the observed population of protoplanetary discs throughout their lifetime requires quantifying the temporal evolution of the turbulent ISM surrounding forming stars. This is challenging for hydrodynamic simulations that are both expensive and limited in integration time by the (turbulent) crossing time of the computational domain.  \nIn this work, we present a semi-analytic approach for tracking the BHL accretion rate over the disc lifetime. In a sense, this is a complement and extension of that of Throop & Bally (2008), who tracked BHL accretion rates onto stars that are gravitationally bound to a dense gas reservoir. Our aim here is to estimate the distribution of accretion histories over the whole disc lifetime for a population of young stars that resembles those surveyed for disc properties in the solar neighbourhood. This includes both the period in which the star is bound to the gaseous overdensity in which it forms (Throop & Bally 2008) and once the bound gas has dispersed. Statistically interpreting the role of BHL accretion on a disc population level requires simulating the turbulent cascade from galactic to protoplanetary disc size scales for the several Myr over which discs persist (Klessen & Hennebelle 2010; Colman et al. 2022; Nusser & Silk 2022). To do so, we apply an excursion set formalism which is effectively zero dimensional (Hopkins 2012). This is therefore a highly efficient way of statistically following density fluctuations of the ISM, making it conducive to semi-analytic planet formation models (e.g. Weder et al. 2023; Mordasini & Burn 2024).  \nIn the remainder of this work, we first present a model for tracking the local ISM conditions and resultant BHL accretion rate for stars evolving in a turbulent ISM (Section 2). We then present a simple disc model and compare to observed disc properties in Section 3. We discuss our results in the broader context of planet formation in Section 4 and draw conclusions in Section 5.", '2.1. Overview': "In this section, our aim is to estimate the distribution of accretion histories for a population of young stars that resembles those surveyed for disc properties in the solar neighbourhood. To do so, we must follow the velocity and density distribution of gas surrounding these young stars. This includes both the period in which the star remains within the star forming region (cf. Throop & Bally 2008), and once the bound giant molecular gas has dispersed. The latter, and to some degree the former, requires quantifying the properties of the supersonic, turbulent ISM.  \nSimulations of sustained, supersonic turbulent flows that are not self-gravitating have clearly shown that the volumeweighted probability density function for density fluctuations is lognormal, with a variance that is a weak function of the Mach number (e.g. Vazquez-Semadeni 1994; Ostriker et al. 1999). A small fraction of the ISM in a steady state is therefore at a sufficient density to marginally overcome (turbulent) pressure support, and undergo gravitational collapse (Evans 1999). Regions which meet these criterion are known as giant molecular clouds (GMCs). These structures have long been known to obey so-called Larson relations, which broadly connect the size scale R cloud and velocity dispersion σ v of molecular clouds in the galactic ISM (Larson 1981). Since original relations, improved data has established  \nthat generally the velocity dispersion σ v ∝ R 1 / 2 cloud (Heyer & Brunt 2004), as expected from supersonic, rapidly cooling turbulence (Burgers 1974), although the normalisation of this relation may somewhat vary with environment (Green et al. 2024). Any models aiming to track turbulent fluctuations in the ISM density and velocity, or the resultant properties of collapsing GMCs, must take into account these relationships.  \nTo track turbulent fluctuations in the ISM, we adopt the excursion set approach described by Hopkins (2012), with the key equations summarised in Appendix A for convenience. In brief, this approach involves quantifying the statistical distribution of densities or velocities at a given location in space as the sum of contributions to that density on larger scales (Appendix A.1). A 'trajectory' is defined as the cumulative sum of perturbations from the largest scale to the smallest scales (equation A9). These trajectories define the density and velocity on a single volume element, although the size of the volume element depends on the scale being considered. For any given trajectory, we can determine whether the local density is sufficient to undergo gravitational collapse (Appendix A.2), becoming a star forming region. Each scale also has an associated turbulent time-scale, based on which we can compute the temporal evolution of the turbulent fluctuations (Appendix A.3). These latter two aspects allow us to estimate local density and velocity fluctuations during and after the star formation process.  \nIn the remainder of this section, we explain the following steps in our procedure:  \n- 1. Quantify the properties of local star forming regions (Section 2.2): Apply the excursion set formalism to establish a statistically representative sample of unstable trajectories, and therefore determine the typical properties of star forming regions that are targeted for local disc surveys (within ∼ 200 pc).\n- 2. Follow the collapse of molecular clouds and star formation (Section 2.3): For a representative sample of star forming regions, with associated trajectories, we adopt a simple model for collapse and the star formation rate. This allows us to draw the formation times for a representative population of stars. Note that a star forming region may contain many stars, and is therefore not equivalent to the protostellar collapse or class 0 phase of stellar evolution.\n- 3. Track the dynamical evolution of stars in the bound cloud (Section 2.4): For the representative stellar population, we track the dynamical evolution of stars based on the gravitational potential of the host star forming region.\n- 4. Evolve temporal fluctuations of the turbulent medium (Section 2.5): Once the cloud is exhausted, the star is ejected into the turbulent medium. We continue to evolve the trajectory to track both density and velocity of the ISM during and after the collapse of the cloud according to equation A13.\n- 5. Compute accretion rate onto the star-disc system (Section 2.6): Based on the local density and velocity during and after the collapse of the star forming region, we calculate the BHL accretion rate onto the star-disc system over time.  \nIn the remainder of this section, we describe these steps in more detail. For each stage, we also highlight caveats to our approach, which amount to points for development in future work.", '2.2. Properties of local star forming regions': "2.2.1. Approach summary  \nWe first need to quantify the mass function of gravitationally unstable GMCs. To do so, we follow the approach of Hopkins (2012), as described in Section 3 of that work. In brief, this involves computing the probability of the turbulent medium collapsing at a given scale, based on the turbulent energy spectrum. We randomly sample a number of 'trajectories', N traj , of which some number N coll are unstable at some length scale R corresponding to total mass M (see Appendix A.2). When weighted by the inverse of the relative volume that this star forming region represents, this yields an estimate for the mass function of unstable GMCs. However, to estimate the statistically sampled upper limit of the mass function among local star forming regions we must also estimate a normalisation constant. The absolute frequency of star forming regions within a certain age range in the solar neighbourhood depends on the turbulent time-scale. This is itself a function of spatial scale:  \nτ t ( R ) = η t cross = η R/ ⟨ v 2 t ( R ) ⟩ 1 / 2 , (1)  \nwhere t cross is the crossing time-scale and η ≈ 1 (Pan & Scannapieco 2010). The root mean square turbulent velocity ⟨ v 2 t ( R ) ⟩ 1 / 2 depends on the turbulent energy spectrum, as described in Appendix A. An appropriate time-scale τ t , 0  \n6  \nFigure 1. Sampling of the mass function via the excursion set formalism for GMCs (blue points) based on 10 5 randomly drawn density perturbation trajectories in the solar neighbourhood. Sampling errors are smaller than the marker size. We have normalised the y -axis such that we show the average number of regions younger than 10 Myr that are expected to be found within 200 pc. We highlight in yellow the region which is well-sampled, corresponding to typical local star forming regions which are surveyed for protoplanetary disc properties. In cyan, we show the range of initial GMC masses we consider to generate our disc population synthesis. In orange, we show the same mass function assuming a fixed star formation efficiency ϵ = 0 . 5.  \n<!-- image -->  \n/circledot  \nover which the trajectory 'resets' corresponds to the turbulent time-scale at given scale. For simplicity, here we adopt τ t , 0 ≈ τ t ( h ), reasoning that fluctuations at a scale R ∼ h are approximately the highest amplitude (we discuss this simplification in Section 2.2.3). We can then estimate the differential number d N of star forming region younger than a certain age ∆ τ within a given volume ∆ V :  \n〈 d N d M 〉 ≈ 1 2 3 ∆ τ τ t , 0 ∆ V ρ M · N coll N traj , (2)  \nwhere ρ and M are the density and mass of the star forming region respectively. We include a factor 1 / 2 3 since given a region of radius R SFR , any other region with R SFR is forbidden inside a distance of 2 R SFR (such a region would merge with the original region).", '2.2.2. Mass function': "We show the outcome of this calculation for 10 5 independent trajectories in Figure 1. We plot the number of star forming regions above a given mass M within a spherical volume with radius R obs = 200 pc ( ∼ h ) younger than 10 Myr. This is approximately the radius within which the majority of protoplanetary disc statistics have been compiled (for a review, see Manara et al. 2023) and the dispersal time-scale for protoplanetary discs (Haisch et al. 2001; Pfalzner et al. 2022). We define the mass threshold below which we expect > 10 star forming regions to be the threshold for being well sampled. This yields a mass limit M local ≈ 500 M ⊙ . Note, our definition of a star forming region is a region that is collapsing to form stars. This is in contrast to a broader complex that may not be collapsing globally. As an illustrative example, the Taurus-Auriga association has a total mass much larger than 500 M ⊙ (Goldsmith et al. 2008). However, the complex as a whole is not collapsing to form stars. Instead collapse and star formation has occurred in several smaller sub-regions. It is these regions, or 'NESTS' (Joncour et al. 2018), that are the individual star forming regions in the sense appropriate to this work. In Taurus, these regions are typically a few tens of solar masses. The older region Upper Scorpius has a total stellar mass ∼ 1400 M ⊙ (Preibisch et al. 2002), but is also composed of several individual populations (Miret-Roig et al. 2022). The nearest individual star forming region with total mass of several 10 3 M ⊙ is the Orion Nebula cluster at a distance of ∼ 400 pc (Hillenbrand & Hartmann 1998; Reid et al.  \n2009). We conclude that, although our approach is approximate in that it neglects some physical processes that may influence local star formation, such as possible nearby supernovae (Zucker et al. 2022), the upper mass limit we infer is reasonable. Our subsequent results do not depend strongly on this limit, as long as we exclude much more massive regions that may be dispersed prematurely by stellar feedback.  \nMotivated by this result, hereafter we we will truncate the GMC mass function above 500 M ⊙ . We will also exclude regions with masses < 10 M ⊙ , for which we expect the turbulent velocity to approach the sound speed ( M∼ 1). In principle, they may still be accurately sampled by the excursion set formalism if the density distribution does not diverge greatly from log-normal. However, we cannot choose an arbitrary small lower mass limit, since for much lower mass regions there is insufficient mass to form stellar mass objects. Experimentation shows that our results are not strongly dependent on our choice of lower mass limit (see Section 3.4.4).", '2.2.3. Caveats for the mass function': 'In the averaging for equation 2 we have ignored that similar mass unstable regions may be spatially correlated, since the probability of collapse on a given spatial scale is dependent on the density contribution from larger scales. However, Hopkins (2012) find that this bias is negligible for the low mass star forming regions, with R SFR ≪ h . The shape of the mass function for the local, low mass GMCs relevant in the context of this work should therefore be correct. Clearly, star forming regions should still preferentially be found where the gas density is highest, while we expect using the average density is sufficient for R obs ≳ h .  \nWe also highlight that our approach to computing the mass function of GMCs is not exact because we do not follow the time evolution of the interstellar medium - i.e. we follow Section 3 of Hopkins (2012) rather than Section 7, which tracks turbulent fluctuations and assembly of merger trees. In principle, in order to compute the mass function of observed star forming regions we should take into account both the variations in the turbulent timescale at different spatial scales, and how the cloud evolves after it reaches critical density. To some degree, this can be folded in by evolving the turbulent medium via perturbations to each trajectory, as described in Appendix A.3, and following Section 7 of Hopkins 2012. However, strictly this would also require self-consistently evolving the cloud, including collapse and growth. While this represents an aim for the future, we here wish to keep our model as simple as possible. We do not expect our simplification in assuming a single turbulent time-scale to strongly influence our results, primarily because the turbulent time-scale is a shallow function of R ( τ t ∝ R 1 / 2 for R ≪ h ), and the mass of the star forming region is a steep function of R ( M ∝ R 3 for R ≪ h ), so that the turbulent time-scale is a very weak function of the mass of the star forming region ( τ t ∝ M 1 / 6 for R ≪ h ). Indeed, we are vindicated in our assumptions in that we approximately recover ⟨ d N/ d M ⟩ ∝ M -2 , as is found observationally across a wide range of galactic environments (e.g. Mok et al. 2020).', '2.3.1.': "Approach summary  \nIn the early stage after a star forms, it occupies the local star forming region that is undergoing gravitational collapse. Since we are interested in the local density and velocity evolution throughout a star's lifetime, we require a simple model to follow the collapse of the GMC and formation of stars. To do so, we first draw a representative sample of GMCs following Section 2.2, each with an associated trajectory that we will subsequently use. For each of these regions, we wish to estimate how the region evolves during the collapse and star formation phase. To do so, we adopt the approach of Girichidis et al. (2014) to estimate the free-fall collapse and resultant star formation rate. The mathematics of the approach are described in Appendix B, but we briefly and qualitatively describe the model here.  \nWe start with a gravitationally unstable GMC of radius R SFR , which is initially the largest scale on which a given turbulent trajectory is unstable (Appendix A.2). This region collapses in free-fall. Physically, it does not do so monolithically, but because smaller sub-regions of the cloud may be higher density (shorter free-fall time-scales), they fragment first and continue to accrete from their parent structures (Hoyle 1953). By assuming that each part of the cloud is in free-fall, we can estimate the fraction of the cloud has collapsed beyond some density threshold to form stars from the initial probability density function for density on the scale R SFR . We define this density threshold by a free-fall time of 0 . 1 Myr, although our results are not sensitive to this choice providing it is much shorter than the global collapse time-scale ( ∼ 1 -3 Myr for our star forming regions - see Figure B1). We emphasise that this does require knowledge of a geometrical density structure during the cloud collapse. We need only assume a geometry when when we compute the evolution of the local density and velocity for individual stars (Section 2.4). We assume the  \ncollapse proceeds until the cloud is instantly dispersed when a star formation efficiency ϵ = 0 . 5 is reached, this value being commensurate with simulations of stellar core formation (Matzner & McKee 2000).", '2.3.2. Caveats for the cloud collapse model': "Our philosophy in this work is to maintain a model that is as simple as possible, building a framework that can be developed in subsequent work. As a result, we have neglected numerous processes that may influence star formation. In particular, the cloud is exposed to heating by feedback processes from the stars as they form. A wide range of analytic and semi-analytic approaches have been developed to describe this process (e.g. V'azquez-Semadeni et al. 2019), and these may be adopted in future work. We are helped in this context in our interest in low mass star forming regions that are locally well sampled. We do not expect many massive stars or supernovae among these local regions, and we are therefore somewhat justified in neglecting feedback processes and exclusively considering gravitational collapse.", '2.4.1. Approach summary': 'Having defined a local star formation rate within an individual star forming region (which may form many stars), we can draw formation times for an initial stellar population. Since these stars are initially bound to the cloud from which they form, we integrate the orbit of the star within the gravitational potential of the cloud, ignoring interactions with other stars. For numerical convenience we will calculate the gravitational potential with a smooth Plummer profile with a core radius a = R SFR / 2 . 6, such that the outer radius of the star forming region is twice the half-mass radius. The radius R SFR evolves following free-fall collapse, as described in Section 2.3 - i.e. R SFR (and therefore a ) evolves from its initial value following equation B4, with r 1 being the initial size of the region. We place stars within the core radius, and we draw an isotropic velocity distribution with a dispersion assuming instantaneous virial equilibrium with the cloud. To track local density and relative velocity of the gas, we then assume the same smooth Plummer density profile (as for Throop & Bally 2008, see discussion in Section 2.4.2). We then evolve the star under the time-evolving potential, until the final dispersal of the cloud (when M ∗ = ϵM tot , for ϵ = 0 . 5). We do not modify the total mass in the region as gas is converted into stars, although given that we end the collapse phase at star formation efficiency ϵ = 0 . 5, this should only have an order unity influence on the total potential.', '2.4.2. Caveats for accretion during the bound phase': "During the bound phase, apart from our inclusion in the broader star forming context, this part of our approach makes similar simplifications as those of Throop & Bally (2008). As did Throop & Bally (2008), we neglect turbulent fluctuations of the medium within the collapsing region. Clearly gravitational collapse of the ISM modifies the local turbulent energy spectrum, and therefore a prescription for these fluctuations would not be straightforward. We expect our approach to recover typical accretion rates within a given region, but underestimates the variance in these rates. Unlike Throop & Bally (2008), we consider the increase in overall density resulting from the collapse of the star forming region, an appropriate age dispersion of the stars, and a self-consistent time-scale in which stars remain in the bound region.  \nAnother related simplification is that we also do not capture the protostellar/class 0 collapse phase, where for the first ∼ 10 5 years the young system is clearly being fed by the environment (Pineda et al. 2023). We simply form our star instantaneously at some time. In principle, a semi-analytic model that follows gravitational collapse and the influence on the local density distribution is possible (e.g. V'azquez-Semadeni et al. 2024). Such a treatment of the formation stage may be able to self-consistently capture both the initial mass function and produce a physically motivated initial disc from the residual angular momentum. This would again complicate our model, and we therefore do not present a full prescription for these early stages in this work. However, we explore the consequences of our assumed initial conditions for our disc evolution model in Section 3.4.1.  \nFinally, we do not include the influence of encounters with other stars. This may be included in the future, although for the local star formation regions that are the topic of this work stellar densities are low.", '2.5. Ambient density and velocity in the unbound ISM': '2.5.1. Approach summary  \nOnce the star forming region is dispersed, all the stars it hosts are instantaneously ejected into the background ISM. The background here is defined by the trajectory which evolves from the initial conditions determined at the  \nmoment of collapse of the star forming region. From this moment of collapse, we must then track the density and velocity perturbations at each scale in the trajectory following the procedure of Hopkins (2012), as described in Appendix A.3. At each time-step we (independently) perturb the velocity and logarithmic density perturbation at each scale stochastically, with a magnitude that depends on the turbulent time-scale. This density structure is not used to compute the ambient environment of the star until the collapsing stat forming region is dispersed.  \nWe continue to evolve the ISM after the dispersal of the cloud and the star moves through the unbound ISM. The relative velocity of the star as it moves through the unbound ISM is the gas velocity at the relevant scale, minus the velocity of the star with respect to the bound cloud at the time of dispersal and the velocity of the cloud itself. The velocity of the cloud is fixed as the gas velocity at the time and scale of cloud collapse. Any star therefore has an associated ambient density and relative velocity at each scale.', '2.5.2. Caveats for the turbulent fluctuations of the unbound ISM': 'The main simplifying assumption for this step is that local star formation does not influence the turbulent fluctuations of the ISM. To varying degrees, this is probably not the case. For example, a local increase in the gravitational potential should drive density enhancements to smaller scales, while stellar feedback drives density enhancements to larger scales. We discuss the role of feedback on the ultimate BHL accretion rates in Sections 3.4.2 and 3.4.3. We expect that our focus on lower mass regions at least somewhat mitigates the exclusion of these influences. Incorporating a prescription for how local star formation shapes the turbulent energy spectrum remains a goal for future work.', '2.6.1. Accretion radius': 'BHL accretion has been the topic of many works (for a review, see Edgar 2004). In the purest case, it is the gravitational focusing of ambient gas as the star passes through, resulting in a stagnation point in the wake of the star within which material is accreted. The instantaneous BHL radius within which material is captured is then defined:  \nR BHL = 2 Gm ∗ ( c 2 s +∆ v 2 gas ) . (3)  \nStreamlines with an impact parameter at the BHL radius R BHL hit a stagnation point along the axis coincident with the relative upstream gas velocity vector. In the ideal case, in an homogeneous medium, they then accrete directly onto the star. In general, material within the cross-section carved out by the accretion radius R acc (where typically R acc = R BHL ) is captured by the star. At any turbulent scale, we can estimate the instantaneous rate at which a star captures material in such a scenario (Bondi 1952; Shima et al. 1985):  \n˙ M BHL ≈ πR 2 acc ρ ( λc 2 s +∆ v 2 gas ) 1 / 2 , (4)  \nwhere λ = e 3 / 2 / 4 ≈ 1 . 12 in the isothermal limit. Material is accreted on the turbulent scale that maximises the accretion rate according to equation 4. We therefore define the ambient gas density ρ gas and relative velocity ∆ v gas to correspond to the quantities on this accretion scale. Before the dispersal of the formation cloud, the density and velocity are defined at the location of the star from the appropriate Plummer profile.  \nIn some circumstances the radius within which capture can occur can be smaller than the BHL radius, R acc < R BHL . Broadly, some considerations that may alter the radius from which a star can accrete or suppress accretion entirely are as follows:  \n- 1. Tidal radius: The star is not isolated, but embedded in a larger scale gravitational potential. At the tidal radius R T , the tidal force of this larger scale structure balances the force induced by the stellar potential. If the tidal radius R T < R BHL , then this radius limits the material that can be captured by the star. If a star is bound to a mass M b at a radius r b , then:  \nR T = r b ( m ∗ 3 M b ) 1 / 3 . (5)  \nInitially, this radius is determined by the collapsing cloud, where r b is the instantaneous location of the star with respect to the centre of the region. Once the cloud has dispersed, we take the galactic values r b = 8 kpc and M b = 10 12 M ⊙ .  \n2  \n-  \n)  \n(  \nFigure 2. The relative velocity (top) and density (middle) of the local gas, for a subset of six stars in our sample of 500 trajectories. The BHL accretion rate (bottom) is normalised to be independent of mass m ∗ , by which lines are coloured. The gas density and velocity are shown at the accretion scale (see text for details). In the middle panel, the horizontal dotted lines are the mean gas density for some high density regions in Taurus (Goldsmith et al. 2008). In the bottom panel, the red solid (dashed) line follows the median (mean) normalised BHL accretion rate for the entire sample, and a normalised two-dimensional histogram.  \n<!-- image -->  \n25  \n- 2. Turbulent scale: Clearly the accretion radius cannot exceed the turbulent size scale under consideration, R t .\n- 3. Photoevaporative wind: High energy photons from the central star can heat up the surface layers of the infalling column and drive a thermal wind, similar to the photoevaporative winds that are driven from the protoplanetary disc surface (e.g. Ercolano & Pascucci 2017).\n- 4. External irradiation: Stellar feedback from neighbouring stars on scales much larger than the disc may heat the ambient gas, increasing the sound speed and reducing the local density.  \nIn our fiducial model, we do not include photoevaporative suppression or large scale feedback, but discuss them in greater detail in Sections 3.4.2 and 3.4.3. Based on the above considerations, we define the accretion radius:  \nR acc = min { R BHL , R T , R t } . (6)', '2.6.2. Accretion histories': 'In order to compute accretion histories, we draw a representative sample of stars from the representative sample of star forming regions. This involves drawing a star from a region with a probability that is proportional to both the total stellar mass of the star forming region, and at a time that is proportional to the instantaneous star formation rate within each region. It is also now necessary to define the mass of a star to calculate accretion rates. Hereafter we draw stellar masses from a log-uniform distribution between 0 . 01 -3 M ⊙ . This distribution is not important because we will resample when comparing to observations, except that we need to fully sample the stellar mass range relevant for drawing comparisons to protoplanetary disc surveys. Although the stellar masses do not evolve self-consistently as the stars accrete in our model, the total mass that they would gain is nearly always a small fraction of the initial stellar mass. We then draw 500 accretion histories from a representative sample of synthetic star forming regions, with a draw probability weighted by the stellar mass of the region.  \nWe show the outcome for the evolution of the local density, relative velocity and BHL accretion rate (normalised to show the accretion rate assuming a solar mass star) for a subset of six stars in Figure 2. Given that density in particular can vary over several orders of magnitude, we also draw an approximate comparison to the typical densities in observed star forming regions. To do so, we adopt averaged density estimates from Table 4 of Goldsmith et al. (2008) for a number of high density areas within the Taurus star forming complex. We convert the mass estimates M enc within the quoted surface area A surf to a density estimate by approximating spherical geometry, such that the volume averaged density is:  \n⟨ ρ enc ⟩ V ∼ 3 M enc 4 π ( A surf /π ) 3 / 2 . (7)  \nHowever, this approximation probably yields an underestimate of the typical density within or close to collapsing star forming regions. This is because the area A surf as defined by Goldsmith et al. is not limited to a single, collapsing star forming region. Instead, it contains both high and low density subregions, resulting in an underestimate of ⟨ ρ enc ⟩ V for the collapsing regions. Nonetheless, the typical ⟨ ρ enc ⟩ V ∼ 3 × 10 -21 cm g -3 is commensurate with the high-end for the local gas densities implied by our model (see the middle panel of Figure 2). We conclude that the local densities for stars in our model are broadly consistent with the those observed in local star forming regions.  \nThe resultant BHL accretion rates, normalised for the a solar mass star, are shown in Figure 2, for a subset of six stars in our sample. The phase in which the star remains within the bound gaseous reservoir typically persists for ∼ 1 Myr. During this phase, local low velocity and high density gas results in an early phase of enhanced BHL accretion at a median rate ˙ M BHL , 1 / 2 ∼ 10 -8 M ⊙ yr -1 (Figure 2, bottom). This median rate then quickly drops to ˙ M BHL , 1 / 2 ∼ 10 -10 -10 -11 M ⊙ yr -1 in the age range 1 -3 Myr for a solar mass star. This is considerably lower than the typical rate of accretion onto the central star, ˙ M acc ∼ 10 -8 -10 -9 M ⊙ yr -1 (Manara et al. 2023). One might therefore assume that BHL accretion no longer plays an important role for disc evolution. However, in a turbulent medium, vorticity of the in-falling material means that accretion is mediated via a (protoplanetary) disc (Krumholz et al. 2006; Kuffmeier et al. 2020). In this case, the change in accretion rate can be averaged over the time-scale τ acc at which angular momentum is transported within the disc. Empirically, this time-scale can be a considerable fraction of the disc lifetime (Almendros-Abad et al. 2024). In this case the median underestimates the true contribution of BHL accretion rate onto the disc to the accretion rate onto the star. By contrast, the mean accretion rate when averaged over all stars in our model is ⟨ ˙ M BHL ⟩ ∼ 10 -8 -10 -9 M ⊙ yr -1 at 1 -3 Myr, as shown as the dashed red  \nline in Figure 2. However, the mean overestimates the typical BHL accretion rate contribution for individual discs, since it is dominated by a small number of rapid accretors. Interpreting the contribution from BHL accretion to the stellar accretion rate for typical star-disc systems is therefore not a simple task. Compounding this problem further, discs that survive to comprise the observed sample may also preferentially be those which have been replenished by the ISM. Interpreting the role of BHL accretion for disc evolution therefore requires further calculation, which is the subject of Section 3.', '3.1. Simple disc model': "In an homogeneous medium the velocities of steamlines are zero at the stagnation point. Captured material therefore has zero angular momentum and accretes directly onto the star. However, in an inhomogeneous medium, vorticity results in retention of a bound disc structure that does not instantly accrete onto the star (Krumholz et al. 2006). For a given fluid element, we may then expect accretion onto the star to be mediated via a protoplanetary disc (Moeckel & Throop 2009; Wijnen et al. 2017).  \nWe now aim to generate a simplified 'population synthesis' of discs shaped by BHL accretion. To do so, we wish to initially remain agnostic as to the driver of stellar accretion (Lynden-Bell & Pringle 1974; Tabone et al. 2022). We therefore consider a stellar accretion rate ˙ M acc = M disc /τ acc , for disc mass M disc and fixed accretion time-scale τ acc . For each of the discs in our sample discussed, we solve the initial value problem:  \n˙ M disc = ˙ M BHL -˙ M acc . (8)  \nTo evolve the disc mass, in our fiducial model we fix the initial M disc ( t = 0) = 0. This is equivalent to the assumption that the initial disc that forms due to core collapse is low mass, or rapidly depleted. This runs counter to the conventional assumption that the protoplanetary disc represents an evolved state of the primordial protostellar disc, which formed during the stellar core collapse. We make this decision as a 'clean' test of our model; we wish to explore the properties of discs arising from BHL accretion alone. Such a model requires no additional free parameters. We emphasise that no specific observational constraints contradict a low initial disc mass. Any such evidence would need to connect the primordial circumstellar material, which forms over ∼ 10 2 years, to the protoplanetary disc, which in the context of this work forms over ∼ 10 5 years and evolves over ∼ 10 6 years. Although (massive) discs around Class 0 stars have been inferred (Maury et al. 2019; Tobin et al. 2020) interpreting these observations in terms of the exact age of the (proto-)star (down to ∼ 10 4 years), the nature of the disc and the time-scale over which any such a disc should survive is not simple. Since few, if any, unambiguous observations exist constraining the properties of primordial discs that are categorically around stars of age ≪ 10 5 years, there is no reason to assume that such a primordial disc should be both long-lived and comparable in mass to the subsequently accreted one. Indeed, it would be a remarkable coincidence if BHL accretion and protostar formation contribute an equal fraction of mass to the disc population. Nonetheless, we also discuss the impact of this assumption on our results in Section 3.4.1.  \nFor τ acc = M disc / ˙ M acc we assume a log-uniform distribution ± 1 dex around log τ acc / 1 Myr = -0 . 5, based on visual inspection of the observed distribution (Almendros-Abad et al. 2024). Empirically, the time-scale τ acc is not strongly correlated with stellar mass, but may be a complex function of the disc properties or stellar age. For example, visually interpreting Figure 12 of Almendros-Abad et al. (2024), the median accretion time-scale appears to be τ acc , 1 / 2 ∼ 0 . 1 Myr for the young region ρ Ophiuchus, while τ acc , 1 / 2 ∼ 1 Myr for the older Lupus region. Indeed, we expect rapid BHL accretion to instigate turbulence in the disc, temporarily decreasing τ acc . However, empirically it is not clear the degree to which any variation of τ acc is a survival bias, discs with shorter τ acc having been depleted (we will show a similar evolution in our model which has fixed τ acc for each star in Section 3.2.2). While exploring time- or in-fall-dependent τ acc is an interesting avenue for future exploration, we initially maintain our simplified model, subsequently exploring the potential role of the ISM in driving a time-dependent stellar accretion rate (Section 4.1).  \nGiven the initial condition and accretion time-scale, we evolve each disc mass according to equation 8 over the entire time period under consideration. However, we assume a simple criterion for disc dispersal: M disc < 3 × 10 -5 M ⊙ . We adopt this criterion when showing disc properties in subsequent analysis. This threshold is motivated empirically by the fact that few discs have masses below this value; our results are not strongly sensitive to this choice within a sensible range (see Appendix C).  \nWe wish to compare our model to observed discs with comparable stellar ages. The ages of star forming regions can both be observationally uncertain (e.g. Soderblom et al. 2014) and should exhibit physical scatter given the finite duration of the star formation process. However, nominally ages for the youngest populations of stars in the ρ Ophiuchus complex are ∼ 0 . 1 -1 Myr (Luhman & Rieke 1999) in Taurus are ∼ 1 -3 Myr (Krolikowski et al. 2021), in Lupus are ≲ 3 Myr (Galli et al. 2020) in Chameleon I have been estimated at ∼ 1 -2 Myr (Galli et al. 2021b) and ∼ 3 -6 Myr (Luhman 2007), and Upper Sco are ∼ 5 -10 Myr (Pecaut & Mamajek 2016; Luhman & Esplin 2020). We compare these regions to our model for ages of the individual stars (not the regions), 0 . 5, 1, 2, 3 and 5 Myr respectively; these choices should be considered illustrative, not precise. For total disc masses we take M disc = 10 2 M dust , where M dust is as quoted by Manara et al. (2023), inferred from the mm continuum flux with standard assumptions for opacity and temperature. In order to draw comparisons with each, we perform a Gaussian kernel density estimate with Scott's rule for bandwidth selection (Scott 1992). Where we compare to data directly, we weight the model data points to yield the same distribution in the stellar mass axis as for the observational dataset. In practice, this means dividing the stellar mass distribution into bins, and defining the weight of model data points in each bin i to be N obs ,i /N mod ,i , where N obs ,i and N model ,i are the number of observed and model points in bin i respectively. The results of these comparisons are shown in Figure 3, in which we find remarkable agreement with the observed disc mass, stellar accretion rate, and outer disc radius distributions. We explore these findings in greater detail as follows.", '3.2.1. Disc masses': "Disc masses are possibly the best metric of success of a disc evolution model, since they are time-integrated (i.e. not sensitive to temporal fluctuations) and relatively straightforwardly defined, albeit with some uncertain factors and observational caveats. Surveys inferring dust mass via mm flux are also relatively complete for several local star forming regions. In the top row of Figure 3, disc masses from our model show remarkable agreement with the observed data at all ages. Disc masses from our model fill-out the observed 'wedge' in the m ∗ -M disc plane.  \nWe show the time evolution of the disc mass distribution in our model more quantitatively in Figure 4, where we plot the cumulative distribution function for disc masses around stars 0 . 5 M ⊙ < m ∗ < 1 M ⊙ . We also show observed distributions for Taurus and Upper Sco, leaving out the other regions for clarity. We can see that at ∼ 1 Myr, the disc mass distribution in our model is very similar to the distribution in Taurus. This decreases after a few Myr until reaching a similar distribution to that seen in Upper Sco between 3 and 9 Myr. Within uncertainties, and given the caveat of a range of stellar ages in each region, our model appears fully consistent with the observed disc mass distributions. More speculatively, we note that the disc takes ∼ 0 . 3 Myr to reach the maximum mass distribution. This is interesting in the context of ρ Ophiuchus, which have discs that are more compact (Testi et al. 2016; Cieza et al. 2019) and lower mass (Williams et al. 2019; Testi et al. 2022) than ∼ 1 -2 Myr old regions. A similar result has been found in Corona Australis (Cazzoletti et al. 2019). Our models suggest that this is plausibly the result of ongoing mass assembly proceeding from outside-in, via BHL accretion.", '3.2.2. Stellar accretion rates': 'The accretion rates from our model are entirely determined by the disc mass distribution, since we have imposed ˙ M acc ∝ M disc /τ acc . They are therefore not an independent test of our model, but a validation or repudiation of this assumption. Given that we have anchored τ acc to observations, and we have found observational agreement with the disc mass distribution, as expected we also find excellent agreement for accretion rates for the majority of regions (middle row, Figure 3). We note minor tension particularly in the case of Chameleon I, for which some accretion rates are somewhat larger than predicted by our model. This is plausibly the result of a younger population than the comparison age we have adopted. For example, Galli et al. (2021b) argue that the older esimate obtained by Luhman (2007) may be explained by systematic distance underestimate (then Chameleon I may be ∼ 1 -2 Myr old). Overall, we find excellent agreement between stellar accretion rates in our model and the observed samples.  \nWe may also ask whether the data supports our assumption of a fixed τ acc distribution with stellar mass and time. It is possible, even probable, that this accretion time-scale does vary over time, particularly if internal disc turbulence is driven to some degree by the material inherited from the ISM (which is plausible, see Figure 8). As previously discussed, empirically there is some variation in τ acc between regions (Almendros-Abad et al. 2024). In this case, the true distribution of τ acc may be a function of the state of the ISM in each specific region. However, even for fixed τ acc for individual stars, the fact that star-disc systems with shorter τ acc lose their discs earlier mean that they drop out of the sample of surviving discs. For our model, this is shown in Figure 5, where we find that once discs are ∼ 3 Myr old, the median τ acc ∼ 1 Myr, a factor three larger than for the initial disc population. This moderate increase is similar  \nFigure 3. Linear kernel density estimates (KDEs) for the disc mass (top), stellar accretion rate (middle) and accretion radii (bottom) versus stellar mass in our model compared to observations (red points, red dashed lines) in local star forming regions. We show results at times 0 . 5, 1, 2 and 3 and 5 Myr. The red points represent the distribution of 100 × dust disc masses (top), stellar accretion rate (middle) in ρ Ophiuchus, Taurus, Lupus, Chameleon I and Upper Scorpius (left to right - Manara et al. 2023). When constructing the KDEs, model points are weighted to yield the same distribution along the stellar mass axis as in the observational samples, except in the bottom row. For comparison with the accretion radii, we show the empirically inferred relationship between CO outer disc radii as a red dashed line in the bottom panels (Andrews 2020).  \n<!-- image -->  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \nto the median for the observed distribution for discs of this age (Almendros-Abad et al. 2024). We conclude that, while τ acc may vary over time for individual stars, such temporal variation is not required to produce the observed disc population.', '3.2.3. Outer disc radii': "We can also consider the distribution of gaseous disc radii R disc we infer from our model. However, the appropriate definition of R disc is less clear than for disc masses or stellar accretion rates in terms of the correspondence between the model and observational data. The observed outer radius in dust is dependent on variations in the dust opacity (Rosotti et al. 2019) and in gas molecular tracers, such as CO, outer radii are dependent on thermochemistry (Trapman et al. 2023). Both are dependent on the definition of 'outer radius', often defined by the radius that encloses some fraction of the total flux. In our simple models, we do not explicitly define an outer radius. Simulations (Krumholz et al. 2006; Pelkonen et al. 2024) and analytic estimates (Padoan et al. 2024) indicate that R disc ∼ R acc ∝ m ∗ when this disc is  \n<!-- image -->  \n/circledot  \nFigure 4. Cumulative distribution of disc masses at different times in our model (solid lines) compared to the observed distribution (dashed lines) for Taurus (blue) and Upper Sco (orange). We include only two observed disc populations for clarify. Shaded region represents the statistical uncertainty interval for observational datasets.Figure 5. Cumulative distribution of the accretion time-scales τ acc at different times for surviving discs in our model. This evolution is purely due to discs with short accretion time-scales dropping out of the sample.  \n<!-- image -->  \naccreted from a turbulent medium, with the disc radius typically a factor order unity smaller than the accretion radius. With the above caveats, we therefore attempt to compare R acc with the observed outer disc radii. While there are fewer measurements of gaseous disc outer radii in local star forming regions, in the bottom row of Figure 3 we compare our model to the empirically-inferred relationship for the outer gas disc radius as inferred from CO molecular line measurements: R CO = 250( m ∗ / 1 M ⊙ ) 0 . 9 au (Andrews 2020) typically a factor ∼ 2 -3 larger than the dust component, although with considerable scatter and some examples of much more extended CO radii (Sanchis et al. 2021). We again find a remarkable similarity between our model results and the observed disc properties, with a tail of large  \nFigure 6. Fraction of surviving discs in our model (solid black line), defined to be discs which have a mass M disc > 3 × 10 -5 M ⊙ for stars with masses in the range 0 . 1 -3 M ⊙ . We also vary this threshold, shown by the dashed and dotted lines (factor 1 / 3 and 3 respectively). We show the evolution of disc fraction given by exp ( -t/τ disc ) for τ disc = 3 Myr (Haisch et al. 2001) and 5 Myr (Ribas et al. 2015) as blue and red lines respectively.  \n<!-- image -->  \nR acc > R CO . The agreement between the normalisation in our model and the data may be somewhat coincidental, since 12 CO is usually optically thick so that R CO probably overestimates where the bulk of the mass is, while R acc also somewhat overestimates R disc . It is nonetheless clear that both the scaling and approximate normalisation of the m ∗ -R disc relationship are predicted directly by BHL accretion.", '3.2.4. Disc lifetimes': "For a fair comparison of disc property distributions, we must also confirm that the total fraction of discs included in our sample is similar in our model to those in star forming regions of a given age. We therefore show the disc fraction as a function of time in Figure 6, by the definition M disc > 3 × 10 -5 M ⊙ . We confirm that our model yields reasonable disc fractions over the range of times we consider in our disc property comparisons. We find a very similar disc fraction evolution to the exponential decay relationships suggested in the literature, typically with dispersal time-scale τ disc ∼ 3 -5 Myr (Haisch et al. 2001; Ribas et al. 2015). In Figure 6, we also show the disc fractions we obtain when multiplying the threshold mass by 1 / 3 and 3, finding only a moderate difference of ∼ 10 percent between the most extreme choices.  \nInterestingly, we find that the disc fraction never reaches 100 percent in our model, but is maximal at ∼ 90 percent. Observations indicate that, even for very young disc populations, the disc fraction never exceeds ∼ 80 percent (Michel et al. 2021). However, when accounting for close binaries, this fraction could be closer to 100 percent for single stars (Kraus et al. 2012). The model prediction of a somewhat suppressed early disc fraction should therefore not be interpreted as a unique signature of BHL accretion, and may also depend on the presence or absence of a primordial disc. We also find a somewhat flatter disc fraction as a function of time than a simple exponential decay with linear time exponent. Some authors to suggest such a flattened functional form for the disc fraction (Pfalzner et al. 2022), supported by a non-negligible disc fractions in clusters that are tens of Myr old (Galli et al. 2021a). Such surviving discs can be easily understood in the context of the stochastic accretion process that defines our model.  \nWe caveat this discussion in that we have not included any end-stage disc dispersal processes; in particular the photoevaporative wind driven by the central star (Clarke et al. 2001; Alexander et al. 2006; Gorti & Hollenbach 2009; Owen et al. 2010; Picogna et al. 2019). We would expect that, at least for very low mass discs, mass-loss may be dominated by such a wind. Our assumption is effectively that mass-loss in the wind is inefficient; in this case we can assume that the lifetime of an observable disc is essentially set by the accretion process in isolation. Once a disc is 'dispersed' in our model, if at some later stage the disc replenishment occurs more rapidly than the true  \nphotoevaporation rate, then our model remains valid. Nonetheless, a future investigation of how photoevaporation may interact with BHL accretion for disc populations is warranted in future work.", '3.3. Replenished mass fraction': "It is not clear immediately from Figure 3 whether disc properties are set by a continuous process of BHL accretion. Does BHL just set the initial disc properties that then simply evolve over our imposed time-scale τ acc , or is disc material constantly replenished over the disc lifetime? This is an important distinction because if disc material is constantly being replenished then, even with well-determined initial conditions, the distribution of disc properties can never be a direct probe of isolated disc evolution.  \nTo understand what fraction of material in observed discs should originate from 'fresh' ISM material, we define a metric f M, 1 / 2 to be the 'half-life replenishment factor'. That is, the fraction of material within the disc that has been accreted within the most recent half of the disc life-time. In the absence of BHL accretion, equation 8 has a trivial analytic solution that we label M isol ( t ), given an initial mass. Then, if stellar accretion of the older material occurs before the accretion of younger material, and if stellar accretion is the only mechanism by which mass is removed from the disc, then at time t 2 the fraction of mass that has been added to the disc since time t 1 is  \nf M ( t 1 | t 2 , M disc ( t 2 )) = 1 -M isol ( t 1 | t 2 , M disc ( t 2 )) /M disc ( t 1 ) . (9)  \nFrom equation 9 we define f M, 1 / 2 such that t 2 = t 1 / 2. We expect this estimation of f M, 1 / 2 to be a lower bound because we ignore any other disc dispersal processes. Processes such as planet formation or disc winds should be expected to reduce M isol , thus increasing f M .  \nApplying this metric in Figure 7, we find a substantial fraction of discs composed of a high proportion of recently captured ISM material. On the left panel we show how the replenishment fraction evolves over time for a random subset of discs. Typically, f M , 1 / 2 ≈ 0 . 5 early on, then the fraction drops once the phase in which stars occupy the bound star forming region finishes. As the disc reduces in mass, discs often see an increase in f M , 1 / 2 later in their lifetime, when relatively modest BHL accretion rates can still provide substantial replenished mass compared to the existing disc mass. Not all discs follow this pattern however, with some discs having both high replenishment fractions ( ≳ 0 . 5) and considerable disc masses ( M disc /m ∗ ≳ 10 -3 ) at relatively old ages ( ≳ 5 Myr). The right panel shows that the fraction of discs composed of mostly replenished material drops from ∼ 70 percent at 1 Myr, to ∼ 30 percent at 2 Myr and ∼ 20 percent for discs older than 3 Myr. This fraction therefore remains substantial throughout the disc lifetime.  \nWe conclude that disc properties are constantly 'reset' for many discs in local star forming regions during their first few Myr of evolution. This would suggest that the distributions of observed disc properties should not be used as a test of theoretical models for isolated disc evolution.", '3.4.1. Class 0 phase and the initial disc mass': "To simplify our model, we have assumed that a star initially forms instantaneously, without a period of growth that corresponds to the protostellar and class 0 phase. Indeed, for our fiducial model, we have made the assumption that the initial disc mass is negligible compared to the subsequently accreted material. As discussed, this assumption is justified in the desire for a 'clean' test of the disc properties resulting from BHL accretion, retaining a simple model with few free parameters, and by the coincidence that an alternative would imply. However, given that this assumption is also unconventional, we can ask whether our conclusion regarding the degree of disc replenishment over time holds in the event of the coincidence that the primordial disc mass is comparable to the subsequently captured mass.  \nWe explore this assumption in Appendix D, repeating our experiment but assuming massive initial discs. While we still do not capture the class 0 disc assembly phase, this is equivalent to the opposite assumption as for our fiducial model; that the majority of the young discs mass is assembled within the first ∼ 10 5 yr. From Figure D2, we conclude that the replenishment fraction is only influenced for young discs ( ≲ 1 Myr), after which the distributions remain similar (cf. Figure 7). Our conclusion is that disc replenishment remains an important process throughout the disc lifetime therefore does not depend on this assumption.  \nFigure 7. The fraction of mass that has been accreted in the last half of the disc life-time, f M, 1 / 2 , as a function of time for 500 discs in our model. On the left, we show the time evolution of f M, 1 / 2 for all 500 discs (faint lines), with a random subset of 10 discs shown boldly for clarity. These lines are coloured by the ratio of the disc mass to stellar mass M disc /m ∗ , shown in the colour bar along the top. On the right hand side we show a cumulative fraction for f M, 1 / 2 at 1, 2, 3 and 5 Myr. Only discs with a total mass M disc > 3 × 10 -5 M ⊙ are included in this cumulative fraction, lower mass discs are assumed to be dispersed.  \n<!-- image -->  \nAs discussed in Section 2.6, photoevaporative winds driven by high energy photons from the accreting star may suppress the BHL accretion flow (e.g. Edgar & Clarke 2004). While quantifying the degree to which BHL accretion is suppressed is a complex problem, and should be the subject of a dedicated study, in Appendix E we outline a simplified approach.  \nWe discuss the results of including the internally driven wind in detail in Appendix E, but highlight two key conclusions as follows. The first is that photoevaporation leads to more discrete, episodic periods of accretion, particularly for high mass stars (Figure E3). When stars are in a sufficiently low density environment, the wind completely suppresses accretion in our model, although in reality some accretion may still be possible. Secondly, the fraction of the disc mass that has been recently accreted is reduced by including the wind at late times (Figure E6 compared to Figure 7). Despite including this wind, our estimate suggests that ∼ 20 percent of discs are still composed of mostly replenished material at an age of 2 -3 Myr. This fraction is not very strongly influenced because the periods where the majority of the mass is accumulated are the rapid BHL accretion phases where the dense accretion column is not disrupted by the wind. Nonetheless, clearly it is plausible that internally driven winds can disrupt the BHL accretion flow, at least for stars ≳ 2 M ⊙ (Figure E2), and therefore this represents an important topic for future study.", '3.4.3. External stellar feedback': 'Apart from internal irradiation, we have excluded heating of the ambient medium by neighbouring stars. For older regions, with a lower ambient density, and more massive regions, with a larger number of OB stars, feedback should eventually become important for setting the ambient gas density and sound speed. This will in turn suppress BHL  \naccretion, perhaps in conjunction with the internal wind discussed in Section 3.4.2. Specifically, for feedback-heated regions, the sound speed and relative gas velocity would be enhanced, reducing R BHL and possibly changing the mode of accretion from BHL to disc sweep-up (Section G.3). While we do not expect this to strongly influence our results for the local regions which are the focus of this work (although perhaps more so in Upper Sco), this is clearly a factor that must be considered when interpreting disc properties in different environments (e.g. van Terwisga & Hacar 2023).', '3.4.4. Lower limit of star forming region mass': 'As shown in Figure B1, star forming regions of lower mass collapse over a shorter time-scale, with higher density. We might then expect the mass accumulation for these regions to be more dominated by this early phase. We have chosen a minimum mass of a collapsing GMC to be 10 M ⊙ , but if we reduce this we include more stars in lower mass regions that may then be more dominated by the initial stages of BHL accretion. However, we test the sensitivity of our results to this lower mass threshold in Appendix F, finding that in a statistical sense our results are not strongly sensitive to our choice.', '3.4.5. Other caveats for the BHL accretion rate': 'There remain numerous caveats for our assumed BHL accretion rate. For example, complicating factors include magnetic fields, self-gravity of the gas, and a change in accretion mechanism when the BHL radius is within the disc radius (sweep-up). We discuss these issues in greater detail in Appendix G; generally we expect deviations of order unity from the nominal BHL accretion rate.', '3.5. Results summary': 'The qualitative results of our disc evolution calculations can be summarised concisely: BHL accretion is capable of reproducing observed disc masses, accretion rates, radii and lifetimes as a function of stellar mass and time. The BHL accretion rate is sufficient to have recently re-supplied a considerable fraction ( ∼ 20 -70 percent) of discs of age 1 -3 Myr with a high fraction ( ≳ 50 percent) of replenished material.', '4.1. Disc turbulence and stellar accretion': 'In our model, we have assumed that disc mass is accreted onto the star with a fixed time-scale τ acc via an undefined angular momentum transport process. This process could be internal to the star-disc system, driven by magnetic winds or by viscosity generated by (magneto-)hydrodynamic instability (Tabone et al. 2022; Lesur et al. 2023). Here we assess whether BHL accretion, and late stage in-fall more generally, could drive turbulence and/or stellar accretion in protoplanetary discs. In terms of the latter, the case against BHL accretion as the origin for stellar accretion (i.e. the accretion of the disc material onto the star) was concisely and influentially laid out by Hartmann et al. (2006), based on two arguments. We consider both of these arguments in Sections 4.1.1 and 4.1.2.', '4.1.1. Angular momentum transport to the inner disc': 'The first argument against BHL accretion as a driver of stellar accretion is that accreted material would presumably need to have negligible angular momentum in order to accrete close to the stellar radius. Both simulations of BHL accretion in a turbulent medium (Krumholz et al. 2006; Kuffmeier et al. 2023; Pelkonen et al. 2024) and the fact that disc radii far exceed the stellar radius imply that this is not the case. This would appear to demonstrate that BHL accretion cannot set stellar accretion rates. However, observationally it is clear that, at least for the highest accretors, (magnetically-driven) winds can be launched from radii ≲ 10 au (see Section 4.1 of Pascucci et al. 2023), and other internal processes may also drive turbulence in the inner disc (e.g. Carr et al. 2004; Najita et al. 2009; Ilee et al. 2014; Bosman et al. 2023). Then in-falling material only needs to reach down to these radii in order to eventually accrete onto the host star. Even if the initial angular momentum of material is significantly larger than this, some angular momentum transport in the outer disc may be driven by in-fall in two ways.  \nThe first way is the misalignment of the angular momentum vector of material accreted at different times. For example, Kuffmeier et al. (2024) recently showed that the angular momentum vector of accreted material in a turbulent medium is typically decorrelated on time-scales ≲ 1 Myr. Clearly, adding fluid elements together with the same magnitude of angular momentum but offset angular momentum vectors will reduce the total specific angular momentum of the system. This can simply be understood as a delayed form of BHL directly onto the star. Such a process would  \npresumably lead to shrinking of the disc over time, and may therefore leave a similar imprint in disc populations as MHD driven winds (Tabone et al. 2022; Manara et al. 2023). It would also result in misalignment between inner and outer discs, which is indeed evident in many discs (Bohn et al. 2022; Villenave et al. 2024). Whether or not this results in a continuous or bursty supply of material to the inner disc would depend on whether warps instigated by in-falling material decay on comparatively long time-scales (e.g. the ∼ 10 5 years estimated by Kimmig & Dullemond 2024), or whether associated instabilities drive high levels of turbulence that rapidly decays (Deng et al. 2020; Fairbairn & Ogilvie 2023).  \nA second way in which BHL accretion may contribute to stellar accretion throughout the disc lifetime is via injecting turbulent energy directly into the disc. We can estimate how much turbulent energy can be supplied to the outer disc by the ISM during BHL accretion by computing the balance between the rate at which energy is deposited compared to the decay of that turbulent energy. The rate at which turbulent energy decays is (Klessen & Hennebelle 2010):  \n˙ E decay = -3 3 / 2 2 M disc v 3 t , disc H disc , (10)  \nwhere H disc ∼ 0 . 1 R disc is the disc scale height and the factor 3 3 / 2 is included because v t , disc is defined to be the turbulent velocity in a single dimension. The in-flowing turbulent energy from the ISM is:  \n˙ E in = 1 2 ˙ M BHL v 2 in . (11)  \nIf the velocity of the flow is dominated by the gravitational acceleration from the central star, then:  \nv 2 in ≈ 2 Gm ∗ R disc . (12)  \nThen, if the fraction of the energy that goes into turbulent motion in the disc is ϵ t = | ˙ E decay / ˙ E in | ≤ 1 and the disc is in a steady state, then (rewriting equation 30 of Klessen & Hennebelle 2010):  \n( v t , disc c s ) 3 ≈ 2 ϵ t 3 3 / 2 ( c s 0 . 2 kms -1 ) -3 ( ˙ M BHL 1 . 9 × 10 -8 M ⊙ yr -1 ) ( M disc 10 -2 m ⋆ ) -1 ( H disc R disc ) . (13)  \nInstantaneous steady state is probably a good approximation, because the decay timescale is  \nτ decay = | E t , disc / ˙ E decay | = H disc /v t , disc ≪ τ disc , (14)  \nwhere  \nE t , disc = 1 2 M disc v 2 t , disc . (15)  \nTherefore, the BHL accretion rate ˙ M BHL needs to remain continuously sufficient to drive any observed turbulence throughout the disc lifetime. However, this requirement is probably met for most discs in the context of our model.  \nWe can estimate the distribution of turbulent velocities statistically by directly computing the energy balance in our model. We can write:  \n˙ E t , disc = ϵ t ˙ E in + ˙ E decay + ˙ E acc , (16)  \n˙ E acc = -1 2 ˙ M acc v 2 t , disc , (17)  \nwhere we include the accretion term:  \nalthough in practice | ˙ E acc | ≪ | ˙ E decay | . Then we can solve the initial value problem, including the mass evolution of the disc as before. In principle, we could then compute the disc evolution fully consistently assuming only the ISM as the driver of internal physics, since we then have α SS or equivalently τ acc . However, in practice, this requires introducing a free parameter ϵ t , which may not be a constant. Strictly, it also means self-consistently evolving the disc surface density (or outer radius). For these reasons, self-consistently solving the disc evolution would considerably complicate our model. We therefore leave this to future work. Nonetheless, we can adopt sensible parameters: ϵ t = 1, R disc = 250( m ∗ / 1 M ⊙ ) au and H disc = 0 . 1 R disc with our fixed τ acc to obtain an estimate for the maximal turbulent disc velocity v t , disc over time.  \nFigure 8. Evolution of the cumulative distribution function for the turbulent Mach number (or viscous α SS , top axis) for the discs in our model, assuming a maximal energy input efficiency ϵ t = 1. Each line colour represents a different disc age in our model, from 0 . 1 to 9 Myr. We show as coloured vertical lines some examples of discs with constraints on their turbulent velocity from the literature, several of which are upper limits (Rosotti 2023).  \n<!-- image -->  \nThe outcome of this exercise is shown in Figure 8, where we find a median α SS ∼ 3 × 10 -3 for ∼ 3 Myr old discs. We can compare this to the τ acc distribution (Figure 5), assuming that:  \nτ acc ≈ 3 . 8 ( α SS 10 -3 ) -1 ( H s /R s 0 . 1 ) -2 ( m ∗ 1 M ⊙ ) -1 / 2 ( R s 50 au ) 3 / 2 Myr , (18)  \nwhere we define H s as the scale height at the viscous scale radius R s , which is some fraction of the outer disc radius. Given the median τ acc in our model at 3 Myr is approximately 1 Myr, the median α SS is approximately sufficient to sustain accretion self-consistently in our model by direct kinetic energy injection alone (depending on disc geometry and injection efficiency). In Figure 8 we also compare to observational constraints for turbulence in the outer regions of some discs. With the exception of DM Tau, which has turbulent velocity v t , disc ∼ 0 . 25 -0 . 33 c s , and IM Lup which has v t , disc ∼ 0 . 4 -0 . 6 c s (Paneque-Carre˜no et al. 2024) attempts to measure turbulence in protoplanetary discs have resulted in upper limits. For example, for the disc around HD 163296 this upper limit is v t , disc ≲ 0 . 05 c s (Flaherty et al. 2015, 2017) and for TW Hydra it is v t , disc ≲ 0 . 08 c s (Flaherty et al. 2018; Teague et al. 2018). For M disc /m ∗ = 0 . 01, H disc /R disc = 0 . 1 and ϵ t = 1, then even substituting a modest ˙ M BHL ∼ 10 -11 M ⊙ yr -1 yields v t , disc ∼ 0 . 03 c s , similar to the empirical upper limit on turbulence in some discs. This is approximately the median BHL accretion rate for a solar mass star in our model at an age of ∼ 3 Myr (see the bottom panel of Figure 2). Typically values of M disc are even smaller than 10 -2 m ∗ , so even greater turbulence may be imparted by the ISM. Broadly, for discs older than ∼ 3 Myr (as is the case for all those with upper limits on the turbulence) we expect relatively low levels of turbulence, while for younger discs, like DM Tau (Zhang et al. 2020) and IM Lup (Mawet et al. 2012), turbulence driven by the ISM can also be larger. We conclude that sufficient turbulent energy is supplied to discs from the ISM to sustain the observational constraints on turbulence in the outer disc.', '4.1.2. Stellar accretion in regions of low gas density': 'The second objection raised by Hartmann et al. (2006) is that stars in ionised HII regions still accrete. The typical density in such regions is far lower, and the sound speed and/or velocity dispersion much larger, such that the BHL  \naccretion rate may be suppressed by several orders of magnitude. The example quoted by Hartmann et al. (2006) are some stars in Trumpler 37 which were reported by Sicilia-Aguilar et al. (2005) to be accreting substantially at a rate ∼ 10 -8 M ⊙ yr -1 . This is despite a sound speed c s , ext ∼ 10 km s -1 (which we now distinguish from the sound speed in the disc c s , disc ) and an electron density N e ∼ 3 -15 cm -3 or ρ gas ∼ 2 × 10 -23 g cm -3 . The corresponding BHL accretion rate would therefore be ˙ M BHL ∼ 10 -13 M ⊙ yr -1 for a solar mass star, far too small to drive observed accretion rates. This would appear to rule out in-falling gas as the origin of angular momentum transport in at least these discs. However, this conclusion requires further investigation for the following reasons.  \nFirstly, as discussed above, turbulence injected into the outer disc may only serve to feed internally mediated transport in the inner disc. Therefore turning off the supply from the reservoir of material accreted from the ISM should not immediately end accretion onto the star.  \nSecondly, the small BHL radius resulting from large c s , ext means that, at least initially after the gas is ionised, accretion onto the disc will be facilitated by sweep-up and not BHL accretion (see Appendix G.3). If we re-calculate the sweep-up rate for a disc outer radius 250 au (expected for a solar mass star), this increases the gas capture rate by two orders of magnitude to ˙ M su ∼ 10 -11 M ⊙ yr -1 . While this remains much smaller than the accretion rates reported by Sicilia-Aguilar et al. (2005) in Trumpler 37, if the HII region is young ( ≲ 0 . 1 -1 Myr - Mellema et al. 2006; Raycheva et al. 2022), then there may still be significant density variations that can lead to a spread in BHL capture rates. Indeed, Sicilia-Aguilar et al. (2005) point out a population of embedded stars in Trumpler 37, which are evidence for at least some residual cold, dense gas.  \nFinally, as the mode of mass accretion changes, so does the rate of turbulent energy injection. This is both because the hot turbulent medium will typically have larger relative velocity ∆ v gas ∼ c s , ext , and the material also injects thermal energy. In addition, if turbulent fluctuations are suppressed in the ionised medium, then accreted material contains little angular momentum and can fall directly into the inner disc regions. To estimate the resultant turbulent velocity within the disc, we can assume that v in = ζc s , ext , where ζ is a factor of order unity. Then:  \n˙ E in = 1 2 ˙ M su ζ 2 c 2 s , ext , (19)  \nwhere ˙ M su is the sweep-up rate. The resultant steady state Mach number is:  \nv t , disc c s , disc ≈ 0 . 016 ϵ 1 / 3 t ( H disc 10 au ) 1 / 3 ( M disc 10 -2 M ⊙ ) -1 / 3 ( ˙ M su 10 -11 M ⊙ yr -1 ) 1 / 3 ( ζc s , ext 10 kms -1 ) 2 / 3 ( c s , disc 0 . 2 kms -1 ) -1 (20)  \nWhile for our chosen normalisation yields α SS ∼ 10 -4 , some combination of a factor three increase in the in-fall velocity or a factor eight in the sweep-up accretion rate would be sufficient to drive more substantial α SS ∼ 10 -3 . Turbulence injected by warping could also still play a role in this context (Deng et al. 2020). It is questionable whether either directly injected or warp-induced turbulence can sustain accretion into the inner disc for long in an ionised medium, but the moderate factors involved motivate further numerical experiment.', '4.1.3. A case study: DM Tau': 'We can also consider how the BHL accretion scenario fits with observations of individual discs. We consider the interesting scase of DM Tau, which is a star of mass ∼ 0 . 4 -0 . 5 M ⊙ hosting a disc with a dust cavity and accreting at a rate ∼ 10 -8 M ⊙ yr -1 (Manara et al. 2014; Francis & van der Marel 2020). It has a turbulent velocity v t , disc ∼ 0 . 25 -0 . 33 c s , or viscous α SS ∼ 0 . 1 (Flaherty et al. 2020). For this disc, a recent and relatively rapid period of accretion would be needed to drive the observed turbulence. Estimating M disc /m ∗ ≈ 0 . 03 (Andrews et al. 2013) this would imply recent accretion from the ISM at a rate ˙ M BHL ≳ 2 × 10 -8 M ⊙ yr -1 that ended in the last ∼ 10 4 years. This is similar to the current stellar accretion rate (Manara et al. 2014), suggesting that the majority of the mass of the present day disc could have originated from this putative in-fall event. This scenario is plausible statistically within the context of our model. However, apart from the argument above, there is currently no particular evidence that DM Tau has undergone recent in-fall. This underlines the importance of a systematic search for observational signatures of late stage in-fall events.', '4.1.4. To what degree does in-fall mediate angular momentum transport?': 'In this section, we have considered whether in-fall can drive the separate but related phenomena of turbulence and stellar accretion in protoplanetary discs. It is clear that there is sufficient energy in the in-falling material to  \ndrive turbulence compatible with what is observed in the surface layers of the outer disc in the few cases where such constraints exist. This is not the same as the statement that this in-fall does drive observed turbulence, to conclude as such will require careful case studies of observational examples. Similarly, understanding whether in-fall can regulate stellar accretion, either by angular momentum extraction or by turbulent energy injection, requires both more detailed theoretical and observational analysis. Probably internal processes must be responsible for angular momentum transport in at least some regions of the disc, but in-fall could play a role in regulating this process for example, by supplying material from the outer disc. We conclude that, although it is far from clear that these fundamental problems in disc evolution are solved by in-fall, ruling out its importance is not presently justified.', '4.2. Appearance of BHL-accreting discs': "A pertinent question is what we expect the appearance of a disc undergoing BHL accretion to be. As previously discussed, when accreting material with with a high vorticity, the resultant disc outer radius can be of order the BHL radius (Krumholz et al. 2006; Kuffmeier et al. 2023). What is more uncertain is what should be the structure of a disc if the angular momentum of the accreted material is not aligned with the existing material. Depending on the relative angular momentum vector, this could result in 'puffed up' discs, or misalignment between the inner and outer regions (Kuffmeier et al. 2021). By contrast, at least the mm-dust component of discs is frequently observed to be extremely settled (Duchˆene et al. 2003; Pinte et al. 2007, 2016; Villenave et al. 2020, 2022; Pizzati et al. 2023). However, (sub-)mm emission is also commonly more centrally concentrated than the gaseous disc by a factor ∼ 2 . 5 (Ansdell et al. 2018; Facchini et al. 2019; Andrews 2020), which may be due to radial drift or optical effects (Trapman et al. 2019; Rosotti et al. 2019). It is not clear whether turbulence driven in the outer disc should be reflected in the dust distribution, that may have drifted inwards and/or settled vertically, even in the presence of the warping induced by in-fall (Aly et al. 2024). In fact, we can consider one of the best known discs, HL Tau (ALMA Partnership et al. 2015), which has an extremely settled dust disc, evidenced by both dust height and polarisation that suggest low α SS ∼ 10 -5 (Pinte et al. 2016; Ueda et al. 2021). It is also one of the earliest known examples of a disc being fed by in-falling gas (Hayashi et al. 1993; Welch et al. 2000; Yen et al. 2019; Garufi et al. 2022) at a substantial rate ( ≳ 5 × 10 -9 M ⊙ yr -1 - Gupta et al. 2024). This example therefore makes a compelling argument that the mm dust component of the disc is not necessarily strongly influenced by in-fall. It also suggests that turbulence may be substantially lower in the mid-plane with respect to the surface layers. This represents an important topic for future study.", '4.3. Observability of in-fall events': 'A valid question to ask is the degree to which the model we present in this work predicts more abundant evidence of late-stage in-fall events from the data that is currently available. In fact, growing evidence suggests that a considerable fraction of discs are undergoing some kind of interaction with their environment. For example, from their recent scattered light survey of NIR polarimetric imaging for 43 discs in Taurus, Garufi et al. (2024) reported that 16 percent of those discs showed evidence of interaction with ambient material (although this sample is not necessarily representative). While inferring observational signatures of late stage in-fall requires full hydrodynamic and radiative transfer simulations (Krieger et al. 2024), we can make a simple estimate for the density of gas we expect to be surrounding discs in our model. For accretion onto a disc from a radius R acc , the maximum free-fall time onto a star of mass m ∗ is:  \nτ ff , max ≈ 2 × 10 3 ( R acc 250 au ) 3 / 2 ( m ∗ 1 M ⊙ ) -1 / 2 years . (21)  \nGiven a BHL accretion rate ˙ M BHL , we can then estimate the column density of material external to the disc:  \nN ext ∼ 2 × 10 19 ϵ -1 geo ( µ 2 . 3 ) -1 ( ˙ M BHL 10 -9 M ⊙ yr -1 ) -1 / 2 ( R acc 250 au ) -1 / 2 ( m ∗ 1 M ⊙ ) -1 / 2 cm -2 , (22)  \nwhere ϵ geo ≤ 1 is a geometric factor accounting for the non-uniform distribution of material and µ is the mean molecular mass. We will assume ϵ geo = 1 here, although based on the diverse ambient structures seen in observations, a wide range of ϵ geo may be plausible (e.g. Garufi et al. 2024; Zurlo et al. 2024).  \nThe typical column density implied by equation 22 is relatively low, such that we would expect that in-fall is often challenging to trace observationally. We may consider the example of CO isotopologues, which are often the brightest tracers used to probe the gas content of protoplanetary discs with ALMA. These isotopologues have a self-shielding  \nFigure 9. The column number density of material assuming a mean molecular mass µ = 2 . 3 and ϵ geo = 1 following equation 22 for the surviving discs in our fiducial model. The colours of the lines represent the age of the discs. The red vertical dashed line represents an approximate column density for photodissociation of 12 CO, based on a self-shielding column of 10 15 cm -2 (van Dishoeck & Black 1988) and abundance relative to hydrogen of 10 -5 (Trapman et al. 2023).  \n<!-- image -->  \ncolumn density of ∼ 10 15 cm -2 (van Dishoeck & Black 1988; Visser et al. 2009), and typical abundance relative to molecular hydrogen ∼ 10 -5 (Trapman et al. 2023). As an order of magnitude estimate, we may therefore expect column densities ≲ 10 20 cm -2 to have the CO photodissociated. From equation 22, we plot the N ext distribution over time resulting from the accretion rates in our model in Figure 9. We see immediately in-falling CO would be photodissociated for the majority of discs in the age range ∼ 1 -3 Myr. Interestingly, in this age range, the fraction of discs above the threshold N ext = 10 20 cm -2 , approximately the CO self-shielding column, drops from approximately 30 percent to close to 0 percent. Given that this is the range of ages in Taurus, this would be commensurate with the inferred ∼ 16 percent of discs with evidence of interaction with the ambient medium from scattered light observations (Garufi et al. 2024). However, the physical column density down to which these observations are sensitive is not clear. While these estimates are approximate, they highlight the importance of deep surveys of the environments surrounding large numbers of discs. Such data could be applied to make statistical inferences about the prevalence of late stage in-fall, correlating it with disc properties and testing the theoretical predictions we have made in this work.', '4.4. Implications': "Our results strongly suggest that disc properties are regulated by a continuous process of ISM accretion. This would have numerous consequences, many of which should be investigated in future work. We briefly summarise some examples as follows:  \n- · Disc sub-structures: A diverse range of substructures have been observed in protoplanetary discs (Andrews 2020). We have demonstrated that the disc is undergoing constant perturbation by in-falling material. We would therefore expect some discs to exhibit structure due to this process (Kuffmeier et al. 2020). For example, by accreting gas from the ISM with a misaligned angular momentum vector we expect disc warping (Dullemond et al. 2022; Kimmig & Dullemond 2024), which may also produce observable shadows (Marino et al. 2015; Bohn et al. 2022; Krieger et al. 2024). Both external structures and disc warping are evident in some discs (Walsh et al. 2017), and it seems probable that they originate from late stage in-fall at least in some instances.\n- · Accretion outbursts and short accretion times: FU Orionis events are short accretion bursts that last ∼ 100 years (Hartmann & Kenyon 1996), named for the archetypal example (Herbig 1966, 1977). In-fall of material has been suggested as the origin of these outbursts (e.g. Vorobyov & Basu 2005), alongside other mechanisms such as stellar encounters (Bonnell & Bastien 1992; Borchert et al. 2022) and thermal instability (Lin et al. 1985; Clarke & Syer 1996; Lodato & Clarke 2004). It is unclear if in-fall events could produce short rise times, but may  \nindirectly contribute by instigating instabilities in the disc (e.g. Deng et al. 2020; Speedie et al. 2023). Extended nebulosity surrounding several stars undergoing outburst supports the role of in-fall in at least some cases (Zurlo et al. 2024; Hales et al. 2024). Outbursting star-disc systems may also be more broadly categorised as a subset of systems with short accretions time-scales (Manara et al. 2023; Almendros-Abad et al. 2024). Among such discs, a lifetime problem appears particularly pronounced for discs around some Herbig stars (Grant et al. 2023). The assumptions of our model mean that we cannot directly explain these short accretion time-scales; we impose this time-scale for our model. Nonetheless, given a fixed time-scale τ acc , discs with short τ acc may be episodically replenished in our model. Thus, to a degree the inferred accretion time-scale is decoupled from the expected frequency of observed discs, which no longer needs to be ∼ τ acc /τ disc , for disc age τ disc . Instead, the frequency of discs with a given τ acc would be dictated by the frequency of in-fall events that lead to rapid mass gain. Indeed, in our model a small number of discs with very short accretion time-scales ( ∼ 0 . 01 Myr) are retained for ∼ 3 Myr (Figure 5).  \n- · Long-lived discs: Stochastic replenishment of discs represents a simple explanation for so-called 'Peter Pan' discs, which are (low mass) stars that are tens of Myr old but still exhibit a gas-rich accretion disc (Silverberg et al. 2020; Lee et al. 2020; Galli et al. 2021a). It is not clear if Peter Pan discs are exclusively around low mass stars, or whether finding discs around M dwarf stars is just down to sampling of the mass function. Speculatively, this phenomenon may be related to the fact that individual massive discs appear across a number of regions of different ages, even where the majority of discs have been depleted (Ansdell et al. 2020). Late stage in-fall seems a viable way to explain some of these long-lived discs. This further motivates deep observational programs to characterise the super-disc scale gas surrounding protoplanetary discs.\n- · Planet properties: Possibly the most exciting consequence of protoplanetary discs that are regulated by BHL accretion would be the potential consequences for planets. Multiple studies have shown that mature planet properties are dependent on stellar kinematics (Hamer & Schlaufman 2019; Winter et al. 2020). Some of these differences could be related to covariance with stellar properties, such as age and metallicity (Miyazaki & Masuda 2023). However, it is not clear that all correlations can be explained by stellar properties alone (Longmore et al. 2021; Winter & Alexander 2021; Zink et al. 2023). If discs are constantly replenished by the ISM, stars that form at high velocity relative to the local rest frame should accrete less mass from the ISM after the dispersal of their formation cloud. This represents a direct theoretical link between planet formation and stellar kinematics.", '5. CONCLUSIONS': 'In this work, we have presented a novel method for computing the evolution of rate at which protoplanetary discs capture material from their surroundings via Bondi-Hoyle-Lyttleton (BHL) accretion. To this end, we have adapted the excursion set formalism developed by Hopkins (2012) to track the local density and velocity evolution of stars throughout the lifetime of the protoplanetary disc. Our three most important conclusions from this experiment are:  \n- 1. Disc masses, radii and stellar accretion rates as a function of stellar mass and age can result directly from BHL accretion.\n- 2. Material in the disc is constantly replenished over its lifetime. We estimate that 20 -70 percent of discs are mostly composed of material captured in the most recent half of their lifetime. This estimate is conservative in the sense that we include only stellar accretion, and ignore processes such as planet formation or winds that may further erode the existing disc.\n- 3. The energy influx due to in-falling material onto the disc is sufficient to drive empirically inferred protoplanetary disc turbulence. If this in-fall produces turbulence in the outer disc, then we expect a broad distribution for α SS , with typical value α SS ∼ 10 -3 , while some discs with much higher α SS ∼ 0 . 1 (as for DM Tau) and lower α SS ∼ 10 -5 are expected. We highlight that stellar accretion is not necessarily driven by turbulence in the outer disc, and whether in-fall can contribute to global angular momentum transport to the inner disc remains an open question.  \nOur approach is highly simplified, representing a novel method that can be built upon and adapted to semi-analytic models in order to include the role of BHL accretion in planet formation. Several processes may alter the rate of  \nreplenishment in the disc, but in particular stellar feedback may curtail late stage accretion for relatively massive stars ( ≳ 2 M ⊙ ) and stars born in regions of many OB stars. These feedback processes must be considered in future works.  \nIf the ISM plays an important role in disc evolution, this represents a considerable change of direction for protoplanetary disc theory, which has largely focused on processes operating on an isolated star-disc system. Consequences of this new picture of disc evolution are fundamental and wide-reaching, relating to every stage of planet formation. The inextricable link between planet formation and the ISM underscores the importance of future observational and theoretical studies aiming to quantify how gas in-fall during formation influences protoplanetary disc structures, angular momentum transport and the observed exoplanet population. In particular, deep observational surveys aiming to correlate extended (in-falling) structures with disc properties and/or stellar accretion rates would test our findings, and are an urgent goal for the future.', 'ACKNOWLEDGEMENTS': "We thank the anonymous referee for an extremely thorough and thoughtful report, which helped us to considerably improve both the content and clarity of this work. We also thank Cathie Clarke, Alessandro Morbidelli, Sarah Jeffreson and Carlo Manara for useful discussions. AJW has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skglyph[suppress]lodowska-Curie grant agreement No 101104656. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement No. 101002188).", 'REFERENCES': 'Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 216, doi: 10.1111/j.1365-2966.2006.10293.x ALMA Partnership, Vlahakis, C., Hunter, T. R., et al. 2015, ApJL, 808, L4, doi: 10.1088/2041-8205/808/1/L4 Almendros-Abad, V., Manara, C. F., Testi, L., et al. 2024, Evolution of the relation between the mass accretion rate and the stellar and disk mass from brown dwarfs to stars. https://arxiv.org/abs/2402.10523  \nAly, H., Nealon, R., & Gonzalez, J.-F. 2024, MNRAS, 527, 4777, doi: 10.1093/mnras/stad3494  \nAndrews, S. M. 2020, ARA&A, 58, 483,  \ndoi: 10.1146/annurev-astro-031220-010302  \nAndrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129,  \ndoi: 10.1088/0004-637X/771/2/129  \nAnsdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21, doi: 10.3847/1538-4357/aab890  \nAnsdell, M., Haworth, T. J., Williams, J. P., et al. 2020, AJ, 160, 248, doi: 10.3847/1538-3881/abb9af  \nBegelman, M. C., & Shlosman, I. 2009, ApJL, 702, L5, doi: 10.1088/0004-637X/702/1/L5  \nBenisty, M., Dominik, C., Follette, K., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 605, doi: 10.48550/arXiv.2203.09991  \nBoccaletti, A., Di Folco, E., Pantin, E., et al. 2020, A&A, 637, L5, doi: 10.1051/0004-6361/202038008  \nBohn, A. J., Benisty, M., Perraut, K., et al. 2022, A&A, 658, A183, doi: 10.1051/0004-6361/202142070  \nBond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440, doi: 10.1086/170520  \nBondi, H. 1952, MNRAS, 112, 195,  \ndoi: 10.1093/mnras/112.2.195  \nBonnell, I., & Bastien, P. 1992, ApJL, 401, L31, doi: 10.1086/186663  \nBorchert, E. M. A., Price, D. J., Pinte, C., & Cuello, N. 2022, MNRAS, 510, L37, doi: 10.1093/mnrasl/slab123  \nBosman, A. D., Appelgren, J., Bergin, E. A., Lambrechts, M., & Johansen, A. 2023, ApJL, 944, L53, doi: 10.3847/2041-8213/acb651  \nBurgers, J. 1974, The Nonlinear Diffusion Equation: Asymptotic Solutions and Statistical Problems, University of Maryland Lecture Series (D. Reidel Publishing Company).  \nhttps://books.google.de/books?id=zK6IswEACAAJ  \nBurleigh, K. J., McKee, C. F., Cunningham, A. J., Lee, A. T., & Klein, R. I. 2017, MNRAS, 468, 717,  \ndoi: 10.1093/mnras/stx439  \nCarr, J. S., Tokunaga, A. T., & Najita, J. 2004, ApJ, 603, 213, doi: 10.1086/381356  \nCastelli, F., & Kurucz, R. L. 2004, ArXiv Astrophysics e-prints  \n- Cazzoletti, P., Manara, C. F., Liu, H. B., et al. 2019, A&A, 626, A11, doi: 10.1051/0004-6361/201935273  \nChiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368,  \ndoi: 10.1086/304869  \nChoi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102', 'A. MODELLING TURBULENT FLUCTUATIONS IN THE ISM': 'In order to quantify the mass distribution of star forming regions and the evolution of the ambient gas density and velocity, we apply an excursion set formalism as adapted by Hopkins (2012), which was in turn adapted from the work of Bond et al. (1991). We briefly describe the approach of Hopkins (2012) here, using the same notation for clarity.', 'A.1. Excursion set formalism': "In the first instance, we note that for a region of space with mean density ρ 0 = Σ 0 / 2 h (where we hereafter fix Σ 0 = 12 M ⊙ pc -2 ) and Mach number M≫ 1 on a spatial scale R ∼ 1 /k , the volume weighted distribution of gas density is well-described by a lognormal distribution:  \nd p ( δ | k ) = 1 σ k √ 2 π exp ( -δ 2 2 σ 2 k ) d δ (A1)  \nwhere  \nThe density dispersion for k can be estimated:  \nσ k ≈ ( ln [ 1 + 3 4 ⟨ v 2 t ( k ) ⟩ c 2 s + κ 2 k -2 ]) 1 / 2 , (A4)  \nwhere we will assume throughout this work that the sound speed c s = 0 . 2 km s -1 . Here v t is the turbulent velocity, for which the average ⟨ v 2 t ( k ) ⟩ ≈ σ 2 g ( k ) ∼ kE ( k ) ∝ k 1 -p where E ( k ) the turbulent energy spectrum. Throughout this work we will fix p = 2, appropriate for supersonic, rapidly cooling turbulence (Burgers 1974). On the smallest scales we modify ⟨ v 2 t ( k ) ⟩ 1 / 2 to be the sum in quadrature of the power-spectrum component and the sound speed due to the local thermal pressure. The second term in the logarithm in equation A4 is a modification of the normal term, which would be proportional the square of Mach number. The additional term in the denominator is proportional to the square of the galactic epicyclic frequency κ , for which we will assume a flat rotation curve with κ = √ 2Ω, adopting the orbital frequency Ω = 2 . 6 × 10 -2 Myr -1 appropriate for the solar neighbourhood. This modification is necessary because at large R the differential rotation κ/k plays a similar role to thermal pressure on small scales.  \nAssuming that the turbulent properties of the medium depend only on local gas properties at that scale, then the distribution of densities averaged over any spatial scale is also log-normal with variance:  \nσ 2 ( R ) = ∫ dln( k ) σ 2 k ( k ) | ˜ W ( k, R ) | 2 (A5)  \nwhere ˜ W ( k, R ) is the Fourier transform of the real-space window function W ( x, R ) defined to include the relevant range of spatial scales. Clearly, the density distribution on progressively smaller spatial scales depends on the contribution of all larger scales. Following Hopkins (2012), we define ˜ W ( k, R ) to be a Fourier-space top-hat function: ˜ W ( k | R ) = 1 if k ≤ R -1 and ˜ W ( k | R ) = 0 if k > R -1 .  \nWith these quantities defined, the density field at some infinitesimal volume can be understood as a random walk through Fourier space, which is simple to quantify for a Gaussian random field with the Fourier top-hat window (Bond et al. 1991). In this case, the probability of a change in logarithmic density perturbation from δ 1 to δ 2 ≡ δ 1 +∆ δ given a change from scale k 1 to k 2 is given by  \np ( δ 1 +∆ δ ) d ∆ δ = 1 √ 2 π ∆ S exp ( -(∆ δ ) 2 2 ∆ S ) d(∆ δ ) (A6)  \n∆ S ≡ S 2 -S 1 ≡ σ 2 ( R 2 ) -σ 2 ( R 1 ) (A7)  \nδ ≡ ln ( ρ ρ 0 ) -〈 ln ( ρ ρ 0 )〉 (A2)  \n〈 ln ( ρ ρ 0 )〉 = -σ 2 k 2 . (A3)  \nFigure A1. The fractional change in the variance ∆ S/S of the Gaussian random field at each spatial scale with varying spatial resolution ∆ R/R . We mark the scale height h as a vertical dotted line. The location at which each ∆ S is computed are shown as scatter points.  \n<!-- image -->  \nwhere we define the variance  \nS ( R ) ≡ σ 2 ( R ) . (A8)  \nWith this result, we can therefore define a 'trajectory' for the variation of density across a range of discretely defined spatial scales. Starting on some sufficiently large initial scale R →∞ , we have δ → 0 and σ ( R ) → 0. We then have the well-defined initial conditions for the walk down through scales starting at some large R max , with increments in R , ∆ R i . A single trajectory is then:  \nδ ( R i ) ≡ R j >R i ∑ j ∆ δ j , (A9)  \nwhere each ∆ δ j is drawn from the probability density function defined by equation A6. We will hereafter define a grid with a maximum scale R max = 5 h down to a minimum scale R min < 0 . 01 pc, each successive turbulent scale in the trajectory being 95 percent of the scale above. As shown in Figure A1, choosing this value ensures that the fractional contribution to the overall variance at each successive spatial scale remains small (∆ S/S < 0 . 1). In other words, any fluctuations in density or velocity on scales between the those for which we calculate the structure of the medium will not contribute significantly to the local field. If instead we use ∆ R/R ≲ 0 . 5, this is no longer true; we may miss scales in our calculation that dominate the local density or velocity on an intermediate scale for R ∼ h (see red lines in Figure A1).  \nWe have defined S ( R ) as the variance in the logarithmic density, but the same argument applies to any property which is Gaussian distributed with a known contribution to the variance on each scale. We therefore apply a similar calculation to compute each component of the velocity field at each scales. The only difference is that the variance is determined directly by the turbulent spectrum (in our case, σ 2 g ( R ) ∝ R ), rather than via equation A4.", 'A.2. Critical density and the mass of collapsing GMCs': 'In order to compute the mass function of collapsing GMCs, we must determine the critical density at which a GMC undergoes gravitational collapse. Given the dispersion relation for density fluctuations in a disc of finite thickness h (e.g. Begelman & Shlosman 2009), it can be shown that unstable collapse occurs if:  \nρ ρ 0 ≥ ρ c ρ 0 = Q 0 ( h ) 2 ˜ κ (1 + ˜ k ) [ σ g ( k ) 2 σ g ( h -1 ) 2 ˜ k + ˜ κ 2 ˜ k -1 ] , (A10)  \nwhere we have defined the dimensionless versions of the epicyclic frequency ˜ κ ≡ κ/ Ω = √ 2 and wavenumber ˜ k ≡ | k | h . For a disc in vertical equilibrium σ g ( h -1 ) = h Ω, and we will fix σ g ( h -1 ) = 6 km s -1 . This is the fiducial case explored by Hopkins (2012) and yields a scale height h = 235 pc similar to the observed thin disc scale height (e.g. Vieira et al. 2023). Then at any given scale, we can define an analogous critical logarithmic density perturbation:  \nδ c ( R ) = ln ( ρ c ρ 0 ) -〈 ln ( ρ ρ 0 )〉 . (A11)  \nThe unstable, collapsing star forming region is defined at the largest scale R SFR for which δ ≥ δ c . Interpolating to the cylindrical geometry if the R SFR ≳ h , the gas mass of the collapsing cloud is defined:  \nM ( ρ | R ) ≡ 4 π ρ ( R ) h 3 [ R 2 2 h 2 + ( 1 + R h ) exp ( -R h ) -1 ] , (A12)  \nwhere ρ = ρ c and R = R SFR .', 'A.3. Evolution of the turbulent medium': 'We can compute the temporal evolution of the ISM by assuming the turbulence is globally steady-state, such that the turbulent velocity cascade is maintained outside of collapsing regions. The probability density function for any Gaussian random variable δ under these conditions obeys a generalized Fokker-Planck equation. We can therefore redefine our trajectory at t +∆ t by computing a random walk step:  \n∆ δ j ( t +∆ t ) = ∆ δ j ( t ) exp ( -∆ t/τ t ) + R √ ∆ S (1 -exp( -2 ∆ t/τ t )) , (A13)  \nwhere R is a Gaussian random number with unity variance. The new trajectory is the sum of the components given by equation A13. We define the time-step to be 2 percent of the shortest turbulent time-scale in our grid. Equation A13 applies both to the logarithmic density perturbations and to each component of the interstellar gas velocity.', 'B. CLOUD COLLAPSE AND STAR FORMATION RATE': "In the early stage after a star forms, it occupies the local star forming region that is undergoing gravitational collapse. Since we are interested in the local density and velocity evolution throughout a star's lifetime, we require a simple model to follow the collapse of the GMC and formation of stars. Once the GMC becomes gravitationally unstable, it collapses on a free-fall time-scale. It does not do so monolithically, but because smaller sub-regions of the cloud may be higher density (and therefore have shorter free-fall time-scales), they fragment first and continue to accrete from their parent structures (Hoyle 1953). Simultaneously, the cloud is exposed to heating by feedback processes from the stars as they form. A wide range of analytic and semi-analytic approaches have been developed to describe this process (see V'azquez-Semadeni et al. 2019, and references therein). We are helped in this context in our interest in low mass star forming regions that are locally well sampled. Since we do not expect many massive stars or supernovae among these local regions, we will hereafter ignore feedback processes and exclusively consider gravitational collapse.  \nWe adopt the approach of Girichidis et al. (2014) in modelling the rate of star formation during the collapse of the GMC. The approach of those authors is attractive both for its simplicity and the self-consistent treatment of the gas density distribution as it evolves under gravity. We briefly review the approach here, referring the reader to the original work for details. We consider a spherical region undergoing free-fall collapse from rest, which may be contained within a larger gravitationally unstable cloud. The outer radius radius r , obeys the equation:  \nd 2 r d t 2 = -GM r 2 = -4 πGρ 1 r 3 1 3 r 2 . (B1)  \nHere, M = ρV = (4 π/ 3) ρr 3 is the fixed mass of the region with volume V , initial density ρ 1 = ρ (0) and radius r 1 = r (0). We have dropped the flattened geometry of equation A12 for simplicity, but also because we are considering low mass (small) regions that have r 1 < h . We define the dimensionless time:  \nθ = t/τ ff , (B2)  \nwhere  \nτ ff = √ 3 π 32 Gρ 1 = ( γρ 1 ) -1 / 2 , (B3)  \nis the free-fall time, and the second equality defines γ for convenience. While there is no exact solution to equation B1, Girichidis et al. (2014) point out that  \nwith p = 2, is a good approximation. If the mass is constant, then the density evolves as:  \nr ( θ ) = r 1 ( 1 -θ 2 ) p/ 3 , (B4)  \nρ ( θ ) = ρ 1 ( 1 -θ 2 ) -p . (B5)  \nρ ∗ ( ρ a , t ) = 2 γρ a t 2 +1 2 γ 2 ρ a t 4 + √ ( 2 γρ a t 2 +1 2 γ 2 ρ a t 4 ) 2 -1 γ 2 t 4 . (B6)  \nTo determine the rate at which mass is accreted onto forming stars, we compute the rate at which mass evolves above some threshold density for accretion ρ a . We will choose ρ a to yield a free-fall time 0 . 1 Myr, although this choice does not strongly influence the star formation rate as long as this time-scale is much shorter than the global free-fall time-scale. The problem is now reduced to finding the initial density ρ ∗ that evolves to ρ a at time t . If p =2, then inverting equation B5:  \nNow we just need to calculate the fraction of the total cloud mass that has ρ > ρ ∗ . Given that we have an initial probability density function in terms of volume (as in equation A1):  \nP V = 1 V tot d V d ρ , (B7)  \nwe can convert this into the mass-weighted version:  \nP M = 1 M tot d M d ρ , (B8)  \nvia the equation:  \nd V d ρ = d V d M d M d ρ = 1 ρ d M d ρ , (B9)  \nfor some total volume V tot or mass M tot . Here we define M tot = M ( ρ c | R SFR ) from equation A12. Now the total mass in stars is just:  \nM ∗ = M tot ∫ ∞ ρ ∗ P M ( ρ ) d ρ = M tot 2 [ 1 + erf ( ln ρ c ρ ∗ + σ 2 k √ 2 σ k )] , (B10)  \nwhere the last equality assumes an initially lognormal probability density distribution in volume. As in Section 2.2, the critical density of the region is ρ c ( R SFR ) and σ k is the dispersion given k ∼ 1 /R SFR , for R SFR the largest scale of the collapsing star forming region. We stop star formation and assume the cloud is immediately dispersed once M ∗ = ϵM , where we adopt ϵ = 0 . 5 consistent with calculations for core formation (Matzner & McKee 2000).  \nWe show the evolution of some example star forming regions according to this simple model in Figure B1. Over the typical range of the masses M SFR of the star forming regions that we consider (10 -500 M ⊙ in our fiducial model), the mean density of the regions when the majority of stars form decreases by an order of magnitude. However, the period of time for which stars typically inhabit the collapsing cloud increases from a few 10 5 years for M SFR ∼ 10 M ⊙ to ∼ 1 Myr at M SFR ∼ 500 M ⊙ . Given that the cumulative amount of mass accreted onto stars is proportional to density and time, we then expect a fairly weak dependence on the total mass accumulated by discs during the bound, star forming phase.", 'C. SENSITIVITY TO DISC DISPERSAL CRITERION': 'We have used a fiducial criterion for disc dispersal (when M disc < 3 × 10 -5 M ⊙ ). While few disc-hosting stars have been found with nominal mass subtantially lower than this, the choice is to some degree arbitrary. Here we consider how reducing this threshold by a moderate factor of three alters the distribution of disc properties and replenishment fractions.  \nWe show a comparison to Figure 3 with M disc > 10 -5 M ⊙ in Figure C1. There are some slight differences in the distributions of disc properties, but they remain qualitatively similar to the observed distributions. For example, the typical accretion radius is slightly smaller in the low mass threshold case. This is because the accretion radius  \nFigure B1. Evolution of the mean density (black lines, LHS axis) and fraction of star formation in three star forming regions (red lines, RHS axis) of different mass in our model. We show the results for three different initial masses of the star forming regions, M SFR .  \n<!-- image -->  \nFigure C1. As in Figure 3, but with a reduced dispersal disc mass threshold, M disc < 10 -5 M ⊙  \n<!-- image -->  \nFigure C2. As in Figure 7, but with a reduced dispersal disc mass threshold, M disc < 10 -5 M ⊙ .  \n<!-- image -->  \nfor the very low mass discs may be set by a period of BHL accretion with comparatively low BHL accretion rates (corresponding to a smaller BHL radius). In Figure C2, we also compare the replenishment fractions to those in Figure 7. We find a small increase in the replenishment fraction for discs, particularly at late times. This is because discs with a higher fraction of replenished material are often lower mass. In summary, we find only minor differences in our model outcomes if we reduce the disc dispersal criterion by a factor few.', 'D. SENSITIVITY TO THE INITIAL DISC MASS:': 'We consider how our conclusions on disc replenishment as a function of time change when we do not assume a negligible initial disc mass, as we assumed in our fiducial model. We therefore repeat our calculations with disc masses initially similar to the observed distribution in young star forming regions. In practice, we adopt a log-normal distribution for the disc mass with 1 dex of scatter around the relationship  \nM disc ( t = 0) = 0 . 01 ( m ∗ 1 M ⊙ ) 2 M ⊙ . (D1)  \nWe show the resultant disc properties in D1, which again are in good agreement with the observed disc properties. Our adopted initial disc masses are approximately the maximum we could adopt while still remaining consistent with the observed distribution of disc masses. The outcome for the half-life disc replenishment factor is shown in Figure D2. Although in this case f M , 1 / 2 is somewhat reduced early on in the disc evolution, overall differences in the distribution over time is minor. The fraction of discs with at least some recently accreted material remains similar, while we find that f M , 1 / 2 ∼ 0 . 2 -0 . 45 throughout the disc lifetime. Therefore, even under the assumption that the primordial disc is coincidentally similar in mass to the subsequently accreted one, and ignoring all internal depletion processes, we still expect a substantial fraction of discs to be composed of recently captured disc material. Our conclusion that the present day disc properties are not a good probe of internal disc physics therefore appears robust.', 'E.1. Feedback overview': 'Several authors have explored how feedback during massive star formation can alter the temperature in accretion flows, and therefore reduce the efficiency of mass accretion from the ISM (Wolfire & Cassinelli 1986; Edgar & Clarke  \nFigure D1. As in Figure 3, but for an initial disc mass M disc ( t = 0) = 0 . 01( m ∗ / 1 M ⊙ ) 2 M ⊙ with 1 dex of scatter.  \n<!-- image -->  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n2003, 2004). For massive stars, when dust is sublimates and high temperatures are reached, the accretion flow can be strongly perturbed or suppressed. In our case, we focus on much lower stellar masses and weaker irradiation, such that we do not expect dust sublimation in the flow. However, we are also interested in much lower accretion rates and lower local ISM density than expected during the formation of massive stars. In some cases, we then expect feedback to limit accretion. Understanding the role of feedback requires a devoted study, but here we present a simple estimate for the degree to which feedback limits BHL accretion flow. First (Section E.3), we consider passive heating of the accretion flow, which we show is never a stringent limit on the BHL accretion rate. More importantly, in Section E.4.1, we consider suppression of the accretion flow through photoevaporative winds, which are launched by high energy photons. Both of these mechanisms depend on the evolution of the feedback from the central star, so we first briefly outline our stellar model that use to compute the evolution of the stellar properties.  \nFigure D2. As in Figure 7, but for an initial disc mass M disc ( t = 0) = 0 . 01( m ∗ / 1 M ⊙ ) 2 M ⊙ with 1 dex of scatter.  \n<!-- image -->', 'E.2. Stellar evolution models': 'To compute the feedback efficacy, we require L ∗ given a stellar mass m ∗ of a certain age. For this, we use the precomputed MESA Isochrones and Stellar Tracks (MIST 1 - Paxton et al. 2011; Dotter 2016; Choi et al. 2016). When we consider high energy photons, we also require a stellar spectrum. We therefore couple the MIST isochrones with the stellar atmosphere models of Castelli & Kurucz (2004). However, these models do not account for the accretion luminosity that can dominate energy output at short wavelengths for low mass stars. As a simple estimate of this accretion luminosity, we will assume:  \nL acc = Gm ∗ ˙ M acc R ∗ , (E1)  \nwhere ˙ M acc is the accretion rate onto the star (not to be confused with ˙ M BHL ), and R ∗ is the stellar radius. We assume that the accretion flow emits as a blackbody with temperature:  \nT acc = ( L acc 4 πR 2 ∗ σ SB ) 1 / 4 . (E2)  \nWhen we come to consider the short wavelength ( λ < 2070 ˚ A) photons that drive stellar winds, we will integrate over the stellar spectrum:  \nwhere L ∗ ,λ is the wavelength-specific stellar luminosity.  \nL ∗ ,< 2070 = ∫ 2070 0 L ∗ ,λ d λ, (E3)  \nFigure E1. The thermal radius R therm , defined to be the radius at which the bolometric stellar luminosity heats the gas sufficiently be able to escape the gravity of the central star, as a function of stellar mass and luminosity. We also show the time-evolution for stars of different initial stellar mass in the m ∗ -L ∗ plane. These points are coloured by the age of the star.  \n<!-- image -->  \n∗  \n/circledot', 'E.3. Thermal radius': 'We now consider whether, given the thermal structure of the flow imposed by stellar irradiation, whether the temperature is ever sufficient to raise the sound speed c s > c s , crit ∼ v esc , where v esc is the escape velocity. In such a scenario, the gas at the stagnation point would not be bound to the star, and would therefore escape rather than undergo capture. To estimate the degree to which passive heating by the star can unbind the accretion flow, we will assume that the flow is instantaneously at rest at the accretion radius R acc . The problem is then very similar to a passively irradiated protoplanetary disc (Chiang & Goldreich 1997). The surface of the column is heated at a rate Q irr , while the disc cools at a rate Q cool at the surface of the accretion column. The different geometry of the accretion flow (for example, a cone that widens far from the star) compared to the disc just means that both the heating and cooling areas at a given radius are increased by the same factor. Then the centre of the column has temperature T col such that:  \nσ SB T 4 col = ϕL ∗ 8 πR 2 acc . (E4)  \nHere L ∗ is the bolometric luminosity, and ϕ in the context of a protoplanetary disc is the flaring angle - i.e. the angle at which the stellar radiation meets the surface structure. The largest value this can be is ϕ = π/ 2, which would correspond to a column where the ISM interior to R acc is optically thin. Although in this case our assumptions on the column geometry are broken, we adopt this as the most restrictive possible case.  \nIt is clear from the above discussion that c s /c s , crit ∝ R 1 / 4 acc , so that there always formally exists some radius for which c s > c s , crit (although in practice this radius may be sufficiently far away such that the central star does not set the ISM temperature). Therefore the accretion radius R acc can always therefore be limited by R therm . The outcome of these calculations are shown in Figure E1. We find that R therm is typically much smaller for low mass stars, but that even for m ∗ = 0 . 1 M ⊙ we have R therm ∼ 5 × 10 4 au ∼ R T , the tidal radius. This means that R therm ≳ R T in most scenarios, and we do not expect the passive stellar irradiation to limit BHL accretion.', 'E.4.1. Wind radius': "A second feedback effect is the photoevaporative wind that can be driven by ionising or photodissociating photons. These photons can heat the outer layers of the accretion flow, and drive mass-loss in a thermal wind that may be comparable to the BHL accretion rate. In principle, computing the wind mass-loss rate ˙ M wind requires considering the full microphysics within the photodissociation region (e.g. Haworth et al. 2018, 2023). Perhaps even more importantly, as discussed by Edgar & Clarke (2004), the highly tangential direction of the streamlines as they approach the star makes the flow vulnerable to the momentum injected from within. A proper prescription of the role of feedback on the BHL accretion flow for low mass ( ≲ 3 M ⊙ ) stars will require a devoted study. We nonetheless present a crude approximation as follows.  \nWe base our estimate on injecting a fixed fraction ϵ wind of energy from high energy photons (wavelength λ < 2070 ˚ A) into unbinding mass at a velocity v launch from the BHL radius. As discussed above, this estimate is crude in that it ignores the geometry of the flow. For the purpose of this work, it should be considered an order of magnitude estimate of the conditions under which winds become important. We will assume v launch ∼ v esc = √ 2 Gm ∗ /R acc . In this case:  \n1 2 ˙ M wind v 2 launch = ϵ wind L ∗ ,< 2070 ⇐⇒ ˙ M wind = ϵ wind L ∗ ,< 2070 R acc Gm ∗ . (E5)  \nBHL accretion can only proceed if ˙ M BHL > ˙ M wind , or:  \nR acc ≳ (2 Gm ∗ ) -3 ( ϵ wind L ∗ ,< 2070 πρ ) 2 = R wind . (E6)  \nThe factor ϵ wind folds in the fractional solid angle covered by the accretion flow and the complex wind physics. Hereafter we will assume ϵ wind = 0 . 1. We have also assumed that the wind is launched at the minimum possible velocity v esc , so that 10 percent of the available energy from high energy photons goes into launching the maximal possible mass in the wind. We expect this to be a conservative assumption from the perspective of maximising replenishment. In practice, many photons will escape along low density lines-of-sight (reducing ϵ wind ), and higher energy photons may also launch more rapid winds with velocity ∼ 10 km s -1 . In making this choice, we therefore explore the opposite extreme compared to the 'no wind' scenario of our fiducial model.", 'E.4.2. BHL accretion rates with wind suppression': 'In order to estimate the role of internally driven photoevaporative winds on limiting the BHL accretion flow, we make the assumption that no accretion cannot proceed if R acc < R wind . We assume ˙ M BHL = 0 if this condition is met. In practice, for efficiency we pre-compute a table of R wind in stellar age-mass space, and then interpolate and re-scale for the instantaneous ambient density. We show examples of the evolution of the BHL accretion rate resulting from this prescription in Figure E3, which can be compared to Figure 2. The key difference is that when the ambient gas density is low, instead of experiencing a low ˙ M BHL , the star accretes no mass at all. We highlight that physically this would not necessarily be the case. Material outside of the BHL radius may also be hot and/or dense, such that pressure gradients would not necessarily drive gas away from the star. In this case, the star may continue to accrete, either by BHL accretion or, if the BHL radius is inside the disc radius, then by the sweep-up of gas (see Section G.3). However, quantifying this process in a star forming region that may be partially ionised on large scales is a complex problem. We neglect this consideration here for the sake of simplicity and in the spirit of a conservative estimate for the BHL accretion rate.  \nFor the example cases shown in Figure E3, we get a sense that photoevaporation particularly influences the accretion rates for more massive stars (as expected). For our overall sample (which is not drawn from a physical initial mass function or observationally representative equivalent), the result is that the median accretion rate quickly becomes negligible on a time-scale of ∼ 1 Myr. However, the mean remains relatively unaffected, since the periods of high accretion are not suppressed by the wind. As before, we require disc evolution calculations to understand the degree to which BHL shapes the disc properties over its life-time.', 'E.4.3. Disc evolution including winds': 'We calculate the disc evolution following the prescription described in Section 3.1. For numerical reasons, solving the initial value problem is achieved by weighting the BHL accretion rate by a factor exp( -R 2 wind /R 2 acc ), to avoid  \nFigure E2. The radius R wind inside which the accretion flow may be suppressed by the wind driven by high energy photons, as a function of stellar mass and stellar luminosity at wavelengths λ < 2070 ˚ A. Points represent the outcome of stellar models with no stellar accretion (circles), accretion rates following ˙ M acc = 10 -8 ( m ∗ / 1 M ⊙ ) 2 M ⊙ yr -1 (triangles) and ˙ M acc = 10 -7 ( m ∗ / 1 M ⊙ ) 2 M ⊙ yr -1 (squares). These points are coloured by the age of the star.  \n<!-- image -->  \n∗  \n/circledot  \ndiscontinuities in the second derivative. We show the disc fraction resulting from our calculation, analogous to that in Figure 6, is shown in Figure E4. Predictably, the overall disc fraction is somewhat reduced by a moderate factor. However, the influence of winds depends strongly on stellar mass. Dividing the disc lifetimes into difference mass bins, we find that discs with host stars of mass 0 . 3 M ⊙ < m ∗ < 1 M ⊙ , which make up a large fraction of the observed samples, show a similar disc fraction evolution to the full sample shown in Figure 6. Discs around higher mass stars are depleted more rapidly, since BHL accretion is suppressed at a lower ambient density. Stars with m ∗ < 0 . 3 M ⊙ are again depleted more rapidly, but this is due to our fixed mass threshold for a disc to be considered dispersed. Low mass stars host lower mass discs, with masses that are more frequently below this threshold.  \nThe disc properties we obtain when including photoevaporative winds are shown in Figure E5. When compared to Figure 3, these are not strongly influenced by the wind. Early and rapid phases of accretion are not strongly influenced by the wind, and therefore the resultant masses for surviving discs (which have typically undergone faster BHL accretion) are also similar to the negligible wind case. However, we can ask whether these discs are composed of less recently accreted material, since stars are more likely to inhabit regions of low ambient density later in their evolution. Indeed, from Figure E6 (compared to Figure 7), we see a clear decrease in the fraction of replenished material. In particular, for older stars ( ≳ 8 Myr), only ∼ 10 percent of discs have accreted more than half of their mass recently. Nonetheless, based on our simple estimate, a substantial fraction of discs ( ∼ 20 -60 percent at ages ∼ 1 -3 Myr) have undergone recent replenishment when photoevaporative winds are considered.', 'F. SENSITIVITY TO THE LOWER MASS LIMIT OF STAR FORMING REGIONS': 'Lower mass star forming regions collapse faster to form stars at a higher ambient density, over a shorter period. We can therefore ask how sensitive our results are to our choice to truncate masses of collapsing star forming regions below 10 M ⊙ . We therefore repeat the disc evolution experiment, but with a lower cut-off of 2 M ⊙ in the masses of star forming region, below which there is not sufficient mass to produce a solar mass star at a star formation efficiency of 50 percent.  \n2  \nFigure E3. As in Figure 2, but including a prescription for the suppression of the BHL accretion by internally-driven photoevaporative winds.  \n<!-- image -->  \nWe show the outcome of our calculations in terms of the disc property distribution in Figure F1 and the replenishment fraction F2. The distribution of properties in the model remain similar to the fiducial case (Figure 3). The replenishment fraction drops somewhat for stars of age ∼ 2 Myr, which can be understood in that fewer stars are still embedded in the bound cloud at intermediate ages. However, overall, the fraction of discs with a a substantial quantity of replenished material remains comparable to our fiducial model. We conclude that our findings are not strongly sensitive on our choice of minimum mass of the star forming regions.  \nFigure E4. As in Figure 6, but including a prescription for the suppression of the BHL accretion by internally-driven photoevaporative winds. We have also separated the disc survival fraction up into different mass bins, shown by different line colours. Black represents all stars in our sample with m ∗ > 0 . 1 M ⊙ .  \n<!-- image -->', 'G. CAVEATS FOR THE BHL ACCRETION RATE': 'In this work, we have considered a simple estimate for the BHL accretion rate. Here we discuss some possible reasons why the rate may deviate from this approximation, and discuss some open questions about the process of capture of interstellar material.', 'G.1. BHL accretion mediated by a disc': 'The nature of BHL accretion has a long history of testing with numerical experiments. In a supersonic, uniform medium, these broadly validate the analytic estimate for the BHL accretion rate within a factor of order unity (Hunt 1971, 1979; Ruffert 1994, 1996; Edgar 2004; El Mellah & Casse 2015). In a turbulent medium, the instantaneous accretion rate may be limited by the local vorticity instead of the velocity (Krumholz et al. 2006). In this case there may be a substantial reduction in the accretion rate, but only when the Mach number is low (i.e. ∆ v ∼ c s ), which is rarely the case in our numerical calculations. In a magnetised medium, the BHL accretion rate may also be reduced by a factor ∼ 2, although this is dependent on the ratio of magnetic to thermal pressure (Lee et al. 2014; Burleigh et al. 2017). Overall, we expect our simple estimates to yield a reasonable approximation for the mass accretion rate.', 'G.2. Cloudlet capture': 'Another possible avenue for the capture of interstellar gas is the capture of small cloudlets via tidal energy dissipation. In pure BHL accretion, the self-gravity of the ISM is ignored, but particularly early in the disc evolution, in a dense environment, the gravitational force exerted by the local medium may be substantial. In this case, during the close passage between a star and a cloud, the orbital energy may be dissipated by non-radial oscillations within the cloud, analogous to the formation of tight binaries in globular clusters (Robe 1968; Fabian et al. 1975), in a process that remains similar across a range of internal density profiles (Lee & Ostriker 1986). Assuming identical physics in the star-cloud to the star-star cases, capture requires a periastron a p during an encounter between a cloud with radius R cloud and mass M cloud and star with mass m ∗ such that (Fabian et al. 1975):  \na p ≲ a capt ≈ R cloud [ Gm ∗ R cloud v 2 ∞ q -1 cloud (1 + q cloud ) ] 1 / 6 . (G1)  \nFigure E5. As in Figure 3, but including a prescription for the suppression of the BHL accretion by internally-driven photoevaporative winds.  \n<!-- image -->  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \nHere q cloud = M cloud /m ∗ , and v ∞ is the relative velocity between star and cloud at infinite separation. Then if a capt > R cloud , the cross section is:  \nσ capt ≈ 2 πGm ∗ (1 + q cloud ) ( a capt -R cloud ) v 2 ∞ . (G2)  \nWe have made the conservative assumption that only gravitationally focused encounters may result in capture. We exclude encounters for which the star passes through the cloud, which would result in BHL accretion, although even in this case tidal effects may be relevant for a dense ISM.  \nThe cloud must be compact enough to be captured fully by the star. We can estimate that the cloud can no longer be captured when the differential potential across the cloud is comparable to the kinetic energy that is dissipated during capture. It is easy to show that this yields R cloud ∼ R BHL . Therefore we can estimate the total rate of capture for a cloud of radius R cloud as:  \nΓ capt ( R cloud | ρ cloud , σ v , ) = ∫ v BHL 0 n clouds ( R cloud ) σ capt ˜ vg 3 (˜ v | σ v ) d˜ v, (G3)  \nwhere  \ng N D (∆ v | σ v ) = 2∆ v N D -1 Γ( N D / 2) ( 1 4 σ 2 v ) N D / 2 exp ( -∆ v 2 4 σ 2 v ) (G4)  \nFigure E6. As in Figure 7, but including a prescription for the suppression of the BHL accretion by internally-driven photoevaporative winds.  \n<!-- image -->  \nis the Maxwell-Boltzmann distributions in N D dimensions. We have cast equation G3 in terms of the cloud density ρ cloud = 3 M cloud / 4 πR 3 cloud . If large scale turbulence dominates density fluctuations, we can estimate that a typical cloud density will be ρ cloud ∼ ρ gas . We have introduced n clouds = 1 / (2 R cloud ) 3 , where 2 R cloud is the wavelength of the turbulent oscillations in the ISM that produce a cloud of radius R cloud . Finally, we have defined the BHL velocity v BHL , which is the velocity that yields a particular R BHL = R cloud when substituted into equation 3. The overall accretion rate is then ˙ M capt = Γ capt M cloud .  \nIn Figure G1, we show the results of these calculations for a solar mass star in an ISM of density ρ gas = 100 M ⊙ pc -3 . We find that capture rates are typically maximised for cloud sizes R cloud ∼ 0 . 1 -1 R BHL . Assuming R acc = R BHL and substituting in ∆ v gas = 1 km s -1 for a solar mass star into equation 4 yields ˙ M BHL ≈ 2 . 4 × 10 -8 M ⊙ yr -1 . Inferring the comparable cloud capture accretion rate from Figure G1 yields a very similar rate. In general, we expect that tidal cloud capture may play a significant role in disc assembly, alongside BHL accretion. This may increase the rate of interstellar gas capture by a factor of order unity, and may be the subject of future numerical experiment.', 'G.3. Face-on accretion': 'A third mechanism for ISM accretion is face-on disc accretion, where an existing disc produces a large physical cross-section (Moeckel & Throop 2009; Wijnen et al. 2017). In this case, the ram pressure between a pre-existing disc and the ISM mediates the accretion process. Face-on accretion may be significant when a star moves through a relatively high velocity medium. In this case, the disc shrinks as a result of deposition of low angular momentum material onto the disc surface (Wijnen et al. 2017). There is no steady state for this process unless material can be constantly replenished. This replenishment of outer disc material may be viscous, but then the quantity of material accreted is limited by the initial disc mass. Replenishment via BHL or cloud capture would keep the disc at a radius R disc ∼ R BHL , and would therefore result in a similar accretion rate to BHL accretion (Moeckel & Throop 2009). However, this may significantly enhance accretion rates when the BHL radius decreases; then R disc > R BHL . Face-on accretion would then temporarily dominate the overall rate, shrinking the disc to maintain a balance between shrinkage  \nFigure F1. As in Figure 3, but with a lower mass limit of the collapsing star forming region of 2 M ⊙ .  \n<!-- image -->  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \n∗  \n/circledot  \nand replenishment. Overall, we may expect some moderate enhancement to the rate of capture and stellar accretion during short phases, although this would not be expected to greatly change the overall disc evolution.  \nFigure F2. As in Figure 7, but with a lower mass limit of the collapsing star forming region of 2 M ⊙ .  \n<!-- image -->  \n/circledot  \n/circledot  \nFigure G1. The expected mass accretion rate ˙ M capt onto a disc around a star of mass m ∗ = 1 M ⊙ due to gravitationally focused cloudlet capture in a gaseous medium of density ρ = 100 M ⊙ pc -3 = 6 . 77 × 10 -21 g cm -3 . The colour bar shows the averaged mass accretion rate, and the rate Γ capt of individual cloudlet capture is shown by contours. The BHL radius is shown as a solid cyan line, while the optimum clouds contributing to mass accretion have R cloud ≈ 0 . 1 R BHL (dashed cyan line).  \n<!-- image -->'}