# Planet Formation Regulated by Galactic-scale Interstellar Turbulence

Andrew J. Winter, Myriam Benisty, Sean M. Andrews

► To cite this version:

Andrew J. Winter, Myriam Benisty, Sean M. Andrews. Planet Formation Regulated by Galactic-scale Interstellar Turbulence. The Astrophysical Journal Letters, 2024, 972, 10.3847/2041-8213/ad6d5d. insu-04836808

HAL Id: insu-04836808

https://insu.hal.science/insu-04836808v1

Submitted on 15 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Planet Formation Regulated by Galactic-scale Interstellar Turbulence

Andrew J. Winter<fnref n="1" /> ID, Myriam Benisty<fnref n="1" />,<fnref n="2" /> ID, and Sean M. Andrews<fnref n="3" /> ID

Received 2024 May 15; revised 2024 August 7; accepted 2024 August 7; published 2024 August 23

## Abstract

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 late-stage infall from the interstellar medium (ISM). To quantify the prevalence of late-stage 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 Bondi–Hoyle–Lyttleton (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 20%–70% 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 $\alpha_{SS} \sim 10^{-5}$ to $10^{-1}$, 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 large-scale turbulence, with broad consequences for planet formation theory. This possibility urgently motivates deep observational surveys to confirm or refute our findings.

*Unified Astronomy Thesaurus concepts: Protoplanetary discs (1300); Star forming regions (1565); Planet formation (1241)*

## 1. Introduction

Planets form and begin their evolution within disks of gas and dust that persist for a few Myr around young stars (e.g., Haisch et al. 2001). The properties of these protoplanetary disks, as well as the physics that shapes them, are 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–disk 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 disks include irradiation from nearby massive stars (Winter & Haworth 2022, and references therein), neighbor star–disk encounters (Cuello et al. 2023, and references therein), and the late-stage infall 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.

Recent 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 disks, on spatial scales several times the disk outer radius. Infall 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 (Hales et al. 2020; Pérez et al. 2020; Dong et al. 2022), and misaligned inner disks (Ginski et al. 2021). Evidence that is suggestive of infall events has been found not only during the earliest stages of disk 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 all mature star–disk systems of age $\sim 1–3$ Myr. A recent estimate from a sample size of 43 disk-hosting stars suggests that approximately 16% are interacting with the ambient interstellar medium (ISM; Garufi et al. 2024). This material should fall onto the disk on a timescale much shorter than the disk 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 $\sim 1$ Myr of the disk lifetime is hinted at by the fact that disks younger than 1 Myr in $\rho$ Ophiuchus and Corona Australis are somewhat more compact and lower in mass than observed in $\sim 1–2$ Myr old regions (Testi et al. 2016, 2022; Cazzoletti et al. 2019; Cieza et al. 2019; Williams et al. 2019). Infall may also help to resolve an apparent mass budget problem for protoplanetary disk populations compared with observed exoplanets (Manara et al. 2018). It is possible that disks contain a substantial fraction of recently accreted material for much of their lifetime, which would have significant consequences for disk evolution models.

Infalling gas may be the consequence of Bondi–Hoyle–Lyttleton (BHL) accretion, where low-velocity gas is captured owing to the gravitational potential of the star (Bondi 1952;

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 06300 Nice, France</footnote>

<footnote n="2">² Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="3">³ Center for Astrophysics—Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>
Edgar 2004). Theoretically, the rate of capture of material from the 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 has shown that late-stage accretion may be a significant source of mass and angular momentum for ~1 Myr old disks in high gas density environments (Kuffmeier et al. 2023; Padoan et al. 2024; Pelkonen et al. 2024). Infall has also been suggested as a mechanism to explain examples of substantial protoplanetary disks 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 disks 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.

In this work, we present a semianalytic approach for tracking the BHL accretion rate over the disk 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 disk lifetime for a population of young stars that resembles those surveyed for disk properties in the solar neighborhood. This includes both the period in which the star is bound to the gaseous overdensity in which it forms (Throop & Bally 2008) and the period once the bound gas has dispersed. Statistically interpreting the role of BHL accretion on a disk population level requires simulating the turbulent cascade from galactic to protoplanetary disk size scales for the several Myr over which disks persist (Klessen & Hennebelle 2010; Colman et al. 2022; Nusser & Silk 2022). To do so, we apply an excursion set formalism that 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 semianalytic planet formation models (e.g., Weder et al. 2023; Mordasini & Burn 2024).

In 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 disk model and compare to observed disk properties in Section 3. We discuss our results in the broader context of planet formation in Section 4, and we draw conclusions in Section 5.

## 2. Modeling the Interstellar Medium

### 2.1. Overview

In this section, our aim is to estimate the distribution of accretion histories for a population of young stars that resemble those surveyed for disk properties in the solar neighborhood. 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 (see Throop & Bally 2008) and the period 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.

Simulations of sustained, supersonic turbulent flows that are not self-gravitating have clearly shown that the volume-weighted 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 that meet this 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_{\text{cloud}}$ and velocity dispersion $\sigma_v$ of molecular clouds in the galactic ISM (Larson 1981). Since the original relations, improved data have established that generally the velocity dispersion $\sigma_v \propto R_{\text{cloud}}^{1/2}$ (Heyer & Brunt 2004), as expected from supersonic, rapidly cooling turbulence (Burgers 1974), although the normalization 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.

To track turbulent fluctuations in the ISM, we adopt the excursion set approach described by Hopkins (2012), with the key equations summarized 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 timescale, based on which we can compute the temporal evolution of the turbulent fluctuations (Appendix A.3). The last two aspects allow us to estimate local density and velocity fluctuations during and after the star formation process.

In the remainder of this section, we explain the following steps in our procedure:

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 disk surveys (within ~200 pc).

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.

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.

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).

5. Compute accretion rate onto the star-disk 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-disk system over time.

In 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

We 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_{\text{traj}}$, of which some number $N_{\text{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 normalization constant. The absolute frequency of star-forming regions within a certain age range in the solar neighborhood depends on the turbulent timescale. This is itself a function of spatial scale:

$$ \tau_t(R) = \eta t_{\text{cross}} = \eta R / \langle v_t^2(R) \rangle^{1/2}, \quad (1) $$

where $t_{\text{cross}}$ is the crossing timescale and $\eta \approx 1$ (Pan & Scannapieco 2010). The rms turbulent velocity $\langle v_t^2(R) \rangle^{1/2}$ depends on the turbulent energy spectrum, as described in Appendix A. An appropriate timescale $\tau_{t,0}$ over which the trajectory “resets” corresponds to the turbulent timescale at a given scale. For simplicity, here we adopt $\tau_{t,0} \approx \tau_t(h)$, reasoning that fluctuations at a scale $R \sim h$ are approximately the highest amplitude (we discuss this simplification in Section 2.2.3). We can then estimate the differential number $dN$ of star-forming regions younger than a certain age $\Delta\tau$ within a given volume $\Delta V$:

$$ \left\langle \frac{dN}{dM} \right\rangle \approx \frac{1}{2^3} \frac{\Delta\tau}{\tau_{t,0}} \frac{\Delta V \rho}{M} \cdot \frac{N_{\text{coll}}}{N_{\text{traj}}}, \quad (2) $$

where $\rho$ 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_{\text{SFR}}$, any other region with $R_{\text{SFR}}$ is forbidden inside a distance of $2R_{\text{SFR}}$ (such a region would merge with the original region).

<figure><img src="image_4.png" /><figcaption>Figure 1. Sampling of the mass function via the excursion set formalism for GMCs (blue circles) based on $10^5$ randomly drawn density perturbation trajectories in the solar neighborhood. Sampling errors are smaller than the marker size. We have normalized 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 that is well sampled, corresponding to typical local star-forming regions that are surveyed for protoplanetary disk properties. In cyan, we show the range of initial GMC masses we consider to generate our disk population synthesis. In orange, we show the same mass function assuming a fixed star formation efficiency $\epsilon = 0.5$.</figcaption></figure>

### 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_{\text{obs}} = 200$ pc ($\sim h$) younger than 10 Myr. This is approximately the radius within which the majority of protoplanetary disk statistics have been compiled (for a review, see Manara et al. 2023) and the dispersal timescale for protoplanetary disks (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_{\text{local}} \approx 500 M_{\odot}$. Note that 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_{\odot}$ (Goldsmith et al. 2008). However, the complex as a whole is not collapsing to form stars. Instead, collapse and star formation have occurred in several smaller subregions. 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 $\sim 1400 M_{\odot}$ (Preibisch et al. 2002), but it is also composed of several individual populations (Miret-Roig et al. 2022). The nearest individual star-forming region with a total mass of several $\times 10^3 M_{\odot}$ is the Orion Nebula Cluster at a distance of $\sim 400$ pc (Hillenbrand & Hartmann 1998; Reid et al. 2009). 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.

Motivated by this result, hereafter we will truncate the GMC mass function above $500 M_\odot$. We will also exclude regions with masses $<10 M_\odot$, for which we expect the turbulent velocity to approach the sound speed ($M \sim 1$). In principle, they may still be accurately sampled by the excursion set formalism if the density distribution does not diverge greatly from lognormal. However, we cannot choose an arbitrarily 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} \ll 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 that using the average density is sufficient for $R_{obs} \gtrsim h$.

We 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 ISM—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 timescale to strongly influence our results, primarily because the turbulent timescale is a shallow function of $R$ ($\tau_t \propto R^{1/2}$ for $R \ll h$), and the mass of the star-forming region is a steep function of $R$ ($M \propto R^3$ for $R \ll h$), so that the turbulent timescale is a very weak function of the mass of the star-forming region ($\tau_t \propto M^{1/6}$ for $R \ll h$). Indeed, we are vindicated in our assumptions in that we approximately recover $\langle dN/dM \rangle \propto M^{-2}$, as is found observationally across a wide range of galactic environments (e.g., Mok et al. 2020).

## 2.3. Collapse of the Giant Molecular Cloud

### 2.3.1. Approach Summary

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. 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 freefall 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.

We 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 freefall. Physically, it does not do so monolithically, but because smaller subregions of the cloud may be higher in density (shorter freefall timescales), they fragment first and continue to accrete from their parent structures (Hoyle 1953). By assuming that each part of the cloud is in freefall, we can estimate the fraction of the cloud that 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 freefall time of 0.1 Myr, although our results are not sensitive to this choice providing that it is much shorter than the global collapse timescale (~1–3 Myr for our star-forming regions—see Figure B1). We emphasize that this does require knowledge of a geometrical density structure during the cloud collapse. We need only assume a geometry when we compute the evolution of the local density and velocity for individual stars (Section 2.4). We assume that the collapse proceeds until the cloud is instantly dispersed when a star formation efficiency $\epsilon = 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 semianalytic approaches have been developed to describe this process (e.g., Vázquez-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. Stellar Dynamics in the Bound Cloud

### 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 freefall 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_* = \epsilon M_{tot}$, for $\epsilon = 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 $\epsilon = 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 to 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 it 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 timescale in which stars remain in the bound region.

Another related simplification is that we also do not capture the protostellar/Class 0 collapse phase, where for the first ~$10^5$ yr 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 semianalytic model that follows gravitational collapse and the influence on the local density distribution is possible (e.g., Vázquez-Semadeni et al. 2024). Such a treatment of the formation stage may be able both to self-consistently capture the initial mass function and to produce a physically motivated initial disk 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 disk evolution model in Section 3.4.1.

Finally, 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

Once 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 that evolves from the initial conditions determined at the moment 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 timescale. This density structure is not used to compute the ambient environment of the star until the collapsing star-forming region is dispersed.

We 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. Accretion Rate

### 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:

$$R_{BHL} = \frac{2Gm_*}{(c_s^2 + \Delta v_{gas}^2)} \quad (3)$$

Streamlines 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 a 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):

$$\dot{M}_{BHL} \approx \pi R_{acc}^2 \rho (\lambda c_s^2 + \Delta v_{gas}^2)^{1/2}, \quad (4)$$

where $\lambda = e^{3/2}/4 \approx 1.12$ in the isothermal limit. Material is accreted on the turbulent scale that maximizes the accretion rate according to Equation (4). We therefore define the ambient gas density $\rho_{gas}$ and relative velocity $\Delta 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.
In 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:

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

$$R_T = r_b \left( \frac{m_*}{3M_b} \right)^{1/3}. \qquad (5)$$

Initially, this radius is determined by the collapsing cloud, where $r_b$ is the instantaneous location of the star with respect to the center of the region. Once the cloud has dispersed, we take the galactic values $r_b = 8$ kpc and $M_b = 10^{12} M_\odot$.

2. *Turbulent scale:* Clearly the accretion radius cannot exceed the turbulent size scale under consideration, $R_T$.

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 disk surface (e.g., Ercolano & Pascucci 2017).

4. *External irradiation:* Stellar feedback from neighboring stars on scales much larger than the disk may heat the ambient gas, increasing the sound speed and reducing the local density.

In our fiducial model, we do not include photoevaporative suppression or large-scale feedback, but we discuss them in greater detail in Sections 3.4.2 and 3.4.3. Based on the above considerations, we define the accretion radius as

$$R_{acc} = \min \{R_{BHL}, R_T, R_t\}. \qquad (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 proportional to the total stellar mass of the star-forming region, at a time with probability proportional to the instantaneous star formation rate within that 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 and $3 M_\odot$. 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 disk 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.

We show the outcome for the evolution of the local density, relative velocity, and BHL accretion rate (normalized 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

$$\langle \rho_{enc} \rangle_V \sim \frac{3M_{enc}}{4\pi(A_{surf}/\pi)^{3/2}}. \qquad (7)$$

However, 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 $\langle \rho_{enc} \rangle_V$ for the collapsing regions. Nonetheless, the typical $\langle \rho_{enc} \rangle_V \sim 3 \times 10^{-21} \text{ 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 those observed in local star-forming regions.

The resultant BHL accretion rates, normalized for 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 $\dot{M}_{BHL, 1/2} \sim 10^{-8} M_\odot \text{ yr}^{-1}$ (Figure 2, bottom). This median rate then quickly drops to $\dot{M}_{BHL, 1/2} \sim 10^{-10} \text{ to } 10^{-11} M_\odot \text{ 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, $\dot{M}_{acc} \sim 10^{-8} \text{ to } 10^{-9} M_\odot \text{ yr}^{-1}$ (Manara et al. 2023). One might therefore assume that BHL accretion no longer plays an important role for disk evolution. However, in a turbulent medium, vorticity of the infalling material means that accretion is mediated via a (protoplanetary) disk (Krumholz et al. 2006; Kuffmeier et al. 2020). In this case, the change in accretion rate can be averaged over the timescale $\tau_{acc}$ at which angular momentum is transported within the disk. Empirically, this timescale can be a considerable fraction of the disk lifetime (Almendros-Abad et al. 2024). In this case the median underestimates the true contribution of BHL accretion rate onto the disk to the accretion rate onto the star. By contrast, the mean accretion rate when averaged over all stars in our model is $\langle \dot{M}_{BHL} \rangle \sim 10^{-8} \text{ to } 10^{-9} M_\odot \text{ yr}^{-1}$ at 1–3 Myr, as shown with the dashed red line in Figure 2. However, the mean overestimates the typical BHL accretion rate contribution for individual disks, 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–disk systems is therefore not a simple task. Compounding this problem further, disks that survive to compose the observed sample may also preferentially be those that have been replenished by the ISM. Interpreting the role of BHL accretion for disk evolution therefore requires further calculation, which is the subject of Section 3.
<figure><img src="image_5.png" /><figcaption>**Figure 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 normalized to be independent of mass $m_*$, by which lines are colored. 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) normalized BHL accretion rate for the entire sample and a normalized two-dimensional histogram.</figcaption></figure>
# 3. Protoplanetary Disk Properties

## 3.1. Simple Disk Model

In a homogeneous medium the velocities of streamlines 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 disk 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 disk (Moeckel & Throop 2009; Wijnen et al. 2017).

We now aim to generate a simplified “population synthesis” of disks 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 $\dot{M}_{\text{acc}} = M_{\text{disk}}/\tau_{\text{acc}}$, for disk mass $M_{\text{disk}}$ and fixed accretion timescale $\tau_{\text{acc}}$. For each of the disks in our sample discussed, we solve the initial value problem:

$$ \dot{M}_{\text{disk}} = \dot{M}_{\text{BHL}} - \dot{M}_{\text{acc}}. \qquad (8) $$

To evolve the disk mass, in our fiducial model we fix the initial $M_{\text{disk}}(t=0)=0$. This is equivalent to the assumption that the initial disk that forms as a result of core collapse is low in mass, or rapidly depleted. This runs counter to the conventional assumption that the protoplanetary disk represents an evolved state of the primordial protostellar disk, 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 disks arising from BHL accretion alone. Such a model requires no additional free parameters. We emphasize that no specific observational constraints contradict a low initial disk mass. Any such evidence would need to connect the primordial circumstellar material, which forms over $\sim 10^2$ yr, to the protoplanetary disk, which in the context of this work forms over $\sim 10^5$ yr and evolves over $\sim 10^6$ yr. Although (massive) disks 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 $\sim 10^4$ yr), the nature of the disk, and the timescale over which any such disk should survive is not simple. Since few, if any, unambiguous observations exist constraining the properties of primordial disks that are categorically around stars of age $\ll 10^5$ yr, there is no reason to assume that such a primordial disk 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 disk population. Nonetheless, we also discuss the impact of this assumption on our results in Section 3.4.1.

For $\tau_{\text{acc}} = M/\dot{M}_{\text{acc}}$ we assume a log-uniform distribution $\pm 1$ dex around $\log \tau_{\text{acc}}/1$ Myr = -0.5, based on visual inspection of the observed distribution (Almendros-Abad et al. 2024). Empirically, the timescale $\tau_{\text{acc}}$ is not strongly correlated with stellar mass but may be a complex function of the disk properties or stellar age. For example, visually interpreting Figure 12 of Almendros-Abad et al. (2024), the median accretion timescale appears to be $\tau_{\text{acc},1/2} \sim 0.1$ Myr for the young region $\rho$ Ophiuchus, while $\tau_{\text{acc},1/2} \sim 1$ Myr for the older Lupus region. Indeed, we expect rapid BHL accretion to instigate turbulence in the disk, temporarily decreasing $\tau_{\text{acc}}$. However, empirically it is not clear the degree to which any variation of $\tau_{\text{acc}}$ is a survival bias, disks with shorter $\tau_{\text{acc}}$

having been depleted (we will show a similar evolution in our model that has fixed $\tau_{\text{acc}}$ for each star in Section 3.2.2). While exploring time- or infall-dependent $\tau_{\text{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).

Given the initial condition and accretion timescale, we evolve each disk mass according to Equation (8) over the entire time period under consideration. However, we assume a simple criterion for disk dispersal: $M_{\text{disk}} < 3 \times 10^{-5} M_{\odot}$. We adopt this criterion when showing disk properties in subsequent analysis. This threshold is motivated empirically by the fact that few disks have masses below this value; our results are not strongly sensitive to this choice within a sensible range (see Appendix C).

## 3.2. Observational Comparison

We wish to compare our model to observed disks with comparable stellar ages. The ages of star-forming regions both can 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 $\rho$ Ophiuchus complex are $\sim 0.1-1$ Myr (Luhman & Rieke 1999) in Taurus are $\sim 1-3$ Myr (Krolikowski et al. 2021), those in Lupus are $\lesssim 3$ Myr (Galli et al. 2020), those in Chameleon I have been estimated at $\sim 1-2$ Myr (Galli et al. 2021a) and $\sim 3-6$ Myr (Luhman 2007), and those in Upper Sco are $\sim 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 disk masses we take $M_{\text{disk}} = 10^2 M_{\text{dust}}$, where $M_{\text{dust}}$ is as quoted by Manara et al. (2023), inferred from the millimeter continuum flux with standard assumptions for opacity and temperature. In order to draw comparisons with each, we perform a Gaussian kernel density estimate (KDE) 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 data set. 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_{\text{obs},i}/N_{\text{mod},i}$, where $N_{\text{obs},i}$ and $N_{\text{mod},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 disk mass, stellar accretion rate, and outer disk radius distributions. We explore these findings in greater detail as follows.

### 3.2.1. Disk Masses

Disk masses are possibly the best metric of success of a disk 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 millimeter flux are also relatively complete for several local star-forming regions. In the top row of Figure 3, disk masses from our model show remarkable agreement with the observed data at all ages. Disk masses from our model fill out the observed “wedge” in the $m_*-M_{\text{disk}}$ plane.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Linear KDEs for the disk mass (top), stellar accretion rate (middle), and accretion radii (bottom) vs. 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, 3, and 5 Myr. The red points represent the distribution of 100× dust disk masses (top) and 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 disk radii as a red dashed line in the bottom panels (Andrews 2020).</figcaption></figure>

We show the time evolution of the disk mass distribution in our model more quantitatively in Figure 4, where we plot the cumulative distribution function for disk masses around stars $0.5 M_{\odot} < m_{*} < 1 M_{\odot}$. 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 disk 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 disk mass distributions. More speculatively, we note that the disk takes ~0.3 Myr to reach the maximum mass distribution. This is interesting in the context of ρ Ophiuchus, which has disks 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 disk mass distribution, since we have imposed $\dot{M}_{\text{acc}} \propto M/\tau_{\text{acc}}$. They are therefore not an independent test of our model but a validation or repudiation of this assumption. Given that we have anchored $\tau_{\text{acc}}$ to observations and have found observational agreement with the disk mass distribution, we also find, as expected, excellent agreement for accretion rates for the majority of regions (middle row of 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. (2021a) argue that the older estimate obtained by Luhman (2007) may be explained by a systematic
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Cumulative distribution of disk 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 disk populations for clarity. Shaded region represents the statistical uncertainty interval for observational data sets.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Cumulative distribution of the accretion timescales $\tau_{\text{acc}}$ at different times for surviving disks in our model. This evolution is purely due to disks with short accretion timescales dropping out of the sample.</figcaption></figure>

distance underestimate (then Chameleon I may be ~1–2 Myr old). Overall, we find excellent agreement between stellar accretion rates in our model and those in the observed samples.

We may also ask whether the data support our assumption of a fixed $\tau_{\text{acc}}$ distribution with stellar mass and time. It is possible, even probable, that this accretion timescale does vary over time, particularly if internal disk turbulence is driven to some degree by the material inherited from the ISM (which is plausible; see Section 4.1.1). As previously discussed, empirically there is some variation in $\tau_{\text{acc}}$ between regions (Almendros-Abad et al.

2024). In this case, the true distribution of $\tau_{\text{acc}}$ may be a function of the state of the ISM in each specific region. However, even for fixed $\tau_{\text{acc}}$ for individual stars, the fact that star–disk systems with shorter $\tau_{\text{acc}}$ lose their disks earlier means that they drop out of the sample of surviving disks. For our model, this is shown in Figure 5, where we find that once disks are ~3 Myr old, the median $\tau_{\text{acc}} \sim 1$ Myr, a factor of three larger than for the initial disk population. This moderate increase is similar to the median for the observed distribution for disks of this age (Almendros-Abad et al. 2024). We conclude that while $\tau_{\text{acc}}$
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Fraction of surviving disks in our model (solid black line), defined to be disks that have a mass $M_{\text{disk}} > 3 \times 10^{-5} M_{\odot}$ for stars with masses in the range 0.1–3 $M_{\odot}$. We also vary this threshold, shown by the dashed and dotted lines (factor 1/3 and 3, respectively). We show the evolution of disk fraction given by $\exp(-t/\tau_{\text{disk}})$ for $\tau_{\text{disk}} = 3$ Myr (Haisch et al. 2001) and 5 Myr (Ribas et al. 2015) as blue and red lines, respectively.</figcaption></figure>

may vary over time for individual stars, such temporal variation
is not required to produce the observed disk population.

the scaling and approximate normalization of the $m^*-R_{\text{disk}}$
relationship are predicted directly by BHL accretion.

3.2.3. Outer Disk Radii

We can also consider the distribution of gaseous disk radii $R_{\text{disk}}$ we infer from our model. However, the appropriate definition of $R_{\text{disk}}$ is less clear than for disk 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_{\text{disk}} \sim R_{\text{acc}} \propto m_*$ when this disk is accreted from a turbulent medium, with the disk radius typically a factor of order unity smaller than the accretion radius. With the above caveats, we therefore attempt to compare $R_{\text{acc}}$ with the observed outer disk radii. While there are fewer measurements of gaseous disk 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 disk radius as inferred from CO molecular line measurements: $R_{\text{CO}} = 250(m_*/1 M_{\odot})^{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 disk properties, with a tail of large $R_{\text{acc}} > R_{\text{CO}}$. The agreement between the normalization in our model and the data may be somewhat coincidental, since $^{12}$CO is usually optically thick so that $R_{\text{CO}}$ probably overestimates where the bulk of the mass is, while $R_{\text{acc}}$ also somewhat overestimates $R_{\text{disk}}$. It is nonetheless clear that both

For a fair comparison of disk property distributions, we must also confirm that the total fraction of disks included in our sample is similar in our model to those in star-forming regions of a given age. We therefore show the disk fraction as a function of time in Figure 6, by the definition $M_{\text{disk}} > 3 \times 10^{-5} M_{\odot}$. We confirm that our model yields reasonable disk fractions over the range of times we consider in our disk property comparisons. We find a very similar disk fraction evolution to the exponential decay relationships suggested in the literature, typically with dispersal timescale $\tau_{\text{disk}} \sim 3-5$ Myr (Haisch et al. 2001; Ribas et al. 2015). In Figure 6, we also show the disk fractions we obtain when multiplying the threshold mass by 1/3 and 3, finding only a moderate difference of ~10% between the most extreme choices.

Interestingly, we find that the disk fraction never reaches 100% in our model but is maximal at ~90%. Observations indicate that, even for very young disk populations, the disk fraction never exceeds ~80% (Michel et al. 2021). However, when accounting for close binaries, this fraction could be closer to 100% for single stars (Kraus et al. 2012). The model prediction of a somewhat suppressed early disk 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 disk. We also find a somewhat flatter disk fraction as a function of time than a simple exponential decay with linear time exponent. Some authors suggest such a flattened functional form for the disk fraction (Pfalzner et al. 2022), supported by a nonnegligible disk fraction in clusters that are tens of Myr old (Galli et al. 2021b). Such surviving disks can be easily understood in the context of the stochastic accretion process that defines our model.

3.2.4. Disk Lifetimes
<figure><img src="image_10.png" /><figcaption>Figure 7. The fraction of mass that has been accreted in the last half of the disk lifetime, $f_{M,1/2}$, as a function of time for 500 disks in our model. On the left, we show the time evolution of $f_{M,1/2}$ for all 500 disks (faint lines), with a random subset of 10 disks shown boldly for clarity. These lines are colored by the ratio of the disk mass to stellar mass $M_{\text{disk}}/m_*$, shown in the color bar along the top. On the right we show a cumulative fraction for $f_{M,1/2}$ at 1, 2, 3, and 5 Myr. Only disks with a total mass $M_{\text{disk}} > 3 \times 10^{-5} M_\odot$ are included in this cumulative fraction; lower-mass disks are assumed to be dispersed.</figcaption></figure>

A caveat for the above discussion in that we have not included any end-stage disk 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 disks, 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 disk is essentially set by the accretion process in isolation. Once a disk is “dispersed” in our model, if at some later stage the disk replenishment occurs more rapidly than the true photoevaporation rate, then our model remains valid. Nonetheless, an investigation of how photoevaporation may interact with BHL accretion for disk populations is warranted in future work.

### 3.3. Replenished Mass Fraction

It is not clear immediately from Figure 3 whether disk properties are set by a continuous process of BHL accretion. Does BHL just set the initial disk properties that then simply evolve over our imposed timescale $\tau_{\text{acc}}$, or is disk material constantly replenished over the disk lifetime? This is an important distinction because if disk material is constantly being replenished, then, even with well-determined initial conditions, the distribution of disk properties can never be a direct probe of isolated disk evolution.

To understand what fraction of material in observed disks 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 disk that has been accreted within the most recent half of the disk lifetime. In the absence of BHL accretion, Equation (8) has a trivial analytic solution that we label $M_{\text{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 disk, then at time $t_2$ the fraction of mass that has been added to the disk since time $t_1$ is

$$f_M(t_1|t_2, M_{\text{disk}}(t_2)) = 1 - M_{\text{isol}}(t_1|t_2, M_{\text{disk}}(t_2))/M_{\text{disk}}(t_1). \quad (9)$$

From 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 disk dispersal processes. Processes such as planet formation or disk winds should be expected to reduce $M_{\text{isol}}$, thus increasing $f_M$.

Applying this metric in Figure 7, we find a substantial fraction of disks 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 disks. Typically, $f_{M,1/2} \approx 0.5$ early on, and then the fraction drops once the phase in which stars occupy the bound star-forming region finishes. As the disk reduces in mass, disks 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 disk mass. Not all disks follow this pattern, however, with some disks having both high replenishment fractions (≥0.5) and considerable disk masses ($M_{disk}/m_* \gtrsim 10^{-3}$) at relatively old ages (≥5 Myr). The right panel shows that the fraction of disks composed of mostly replenished material drops from ~70% at 1 Myr to ~30% at 2 Myr and ~20% for disks older than 3 Myr. This fraction therefore remains substantial throughout the disk lifetime.

We conclude that disk properties are constantly “reset” for many disks in local star-forming regions during their first few Myr of evolution. This would suggest that the distributions of observed disk properties should not be used as a test of theoretical models for isolated disk evolution.

## 3.4. Caveats and Tests

### 3.4.1. Class 0 Phase and the Initial Disk 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 disk mass is negligible compared to the subsequently accreted material. As discussed, this assumption is justified in the desire for a “clean” test of the disk 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 disk replenishment over time holds in the event of the coincidence that the primordial disk mass is comparable to the subsequently captured mass.

We explore this assumption in Appendix D, repeating our experiment but assuming massive initial disks. While we still do not capture the Class 0 disk assembly phase, this is equivalent to the opposite assumption to that for our fiducial model, that the majority of the young disks' mass is assembled within the first ~10⁵ yr. From Figure D2, we conclude that the replenishment fraction is only influenced for young disks (≤1 Myr), after which the distributions remain similar (see Figure 7). Our conclusion is that disk replenishment remains an important process throughout the disk lifetime and therefore does not depend on this assumption.

### 3.4.2. Internal Winds

As 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.

We discuss the results of including the internally driven wind in detail in Appendix E, but we 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. Second, the fraction of the disk 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% of disks 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 $\gtrsim 2 M_{\odot}$ (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 neighboring 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 accretion, 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 disk sweep-up (Appendix G.3). While we do not expect this to strongly influence our results for the local regions that are the focus of this work (although perhaps more so in Upper Sco), this is clearly a factor that must be considered when interpreting disk 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 timescale, 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_{\odot}$, 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 disk 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 disk evolution calculations can be summarized concisely: BHL accretion is capable of reproducing observed disk masses, accretion rates, radii, and lifetimes as a function of stellar mass and time. The BHL accretion rate is sufficient to have recently resupplied a considerable fraction (~20%–70%) of disks of age 1–3 Myr with a high fraction (≥50%) of replenished material.
# 4. Discussion

## 4.1. Disk Turbulence and Stellar Accretion

In our model, we have assumed that disk mass is accreted onto the star with a fixed timescale $\tau_{acc}$ via an undefined angular momentum transport process. This process could be internal to the star-disk 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 infall more generally, could drive turbulence and/or stellar accretion in protoplanetary disks. In terms of the latter, the case against BHL accretion as the origin for stellar accretion (i.e., the accretion of the disk 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 Disk

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 disk 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 $\lesssim 10$ au (see Section 4.1 of Pascucci et al. 2023), and other internal processes may also drive turbulence in the inner disk (e.g., Carr et al. 2004; Najita et al. 2009; Ilee et al. 2014; Bosman et al. 2023). Then, infalling 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 disk may be driven by infall in two ways.

The 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 timescales $\lesssim 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 presumably lead to shrinking of the disk over time and may therefore leave a similar imprint in disk populations to MHD-driven winds (Tabone et al. 2022; Manara et al. 2023). It would also result in misalignment between inner and outer disks, which is indeed evident in many disks (Bohn et al. 2022; Villenave et al. 2024). Whether or not this results in a continuous or bursty supply of material to the inner disk would depend on whether warps instigated by infalling material decay on comparatively long timescales (e.g., the $\sim 10^5$ yr estimated by Kimmig & Dullemond 2024), or whether associated instabilities drive high levels of turbulence that rapidly decays (Deng et al. 2020; Fairbairn & Ogilvie 2023).

A second way in which BHL accretion may contribute to stellar accretion throughout the disk lifetime is via injecting turbulent energy directly into the disk. We can estimate how

much turbulent energy can be supplied to the outer disk by the ISM during BHL accretion by computing the balance between the rate at which energy is deposited and the decay of that turbulent energy. The rate at which turbulent energy decays is (Klessen & Hennebelle 2010)

$$ \dot{E}_{\text{decay}} = -\frac{3^{3/2} M_{\text{disk}} v_{\text{t,disk}}^3}{2 H_{\text{disk}}}, \quad (10) $$

where $H_{disk} \sim 0.1 R_{disk}$ is the disk scale height and the factor $3^{3/2}$ is included because $v_{t,disk}$ is defined to be the turbulent velocity in a single dimension. The in-flowing turbulent energy from the ISM is

$$ \dot{E}_{\text{in}} = \frac{1}{2} \dot{M}_{\text{BHL}} v_{\text{in}}^2. \quad (11) $$

If the velocity of the flow is dominated by the gravitational acceleration from the central star, then

$$ v_{\text{in}}^2 \approx \frac{2Gm_*}{R_{\text{disk}}}. \quad (12) $$

Then, if the fraction of the energy that goes into turbulent motion in the disk is $\epsilon_t = |\dot{E}_{\text{decay}}/\dot{E}_{\text{in}}| \le 1$ and the disk is in a steady state, then (rewriting Equation (30) of Klessen & Hennebelle 2010)

$$ \left(\frac{v_{t,disk}}{c_s}\right)^3 \approx \frac{2\epsilon_t}{3^{3/2}} \left(\frac{c_s}{0.2 \text{ km s}^{-1}}\right)^{-3} \left(\frac{\dot{M}_{\text{BHL}}}{1.9 \times 10^{-8} M_{\odot} \text{ yr}^{-1}}\right) \times \left(\frac{M_{\text{disk}}}{10^{-2} m_*}\right)^{-1} \left(\frac{H_{\text{disk}}}{R_{\text{disk}}}\right). \quad (13) $$

Instantaneous steady state is probably a good approximation, because the decay timescale is

$$ \tau_{\text{decay}} = |E_{t,disk}/\dot{E}_{\text{decay}}| = H_{\text{disk}}/v_{t,disk} \ll \tau_{\text{disk}}, \quad (14) $$

where

$$ E_{t,disk} = \frac{1}{2} M_{\text{disk}} v_{t,disk}^2. \quad (15) $$

Therefore, the BHL accretion rate $\dot{M}_{BHL}$ needs to remain continuously sufficient to drive any observed turbulence throughout the disk lifetime. However, this requirement is probably met for most disks in the context of our model.

We can estimate the distribution of turbulent velocities statistically by directly computing the energy balance in our model. We can write

$$ \dot{E}_{t,disk} = \epsilon_t \dot{E}_{\text{in}} + \dot{E}_{\text{decay}} + \dot{E}_{\text{acc}}, \quad (16) $$

where we include the accretion term

$$ \dot{E}_{\text{acc}} = -\frac{1}{2} \dot{M}_{\text{acc}} v_{t,disk}^2, \quad (17) $$

although in practice $|\dot{E}_{acc}| \ll |\dot{E}_{decay}|$. Then, we can solve the initial value problem, including the mass evolution of the disk as before. In principle, we could then compute the disk evolution fully consistently assuming only the ISM as the driver of internal physics, since we then have $\alpha_{SS}$ or equivalently $\tau_{acc}$. However, in practice, this requires introducing a free parameter $\epsilon_t$, which may not be a constant. Strictly,
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Evolution of the cumulative distribution function for the turbulent Mach number (or viscous $\alpha_{SS}$; top axis) for the disks in our model, assuming a maximal energy input efficiency $\epsilon_t = 1$. Each line color represents a different disk age in our model, from 0.1 to 9 Myr. We show as colored vertical lines some examples of disks with constraints on their turbulent velocity from the literature, several of which are upper limits (Rosotti 2023).</figcaption></figure>

it also means self-consistently evolving the disk surface density (or outer radius). For these reasons, self-consistently solving the disk evolution would considerably complicate our model. We therefore leave this to future work. Nonetheless, we can adopt sensible parameters, $\epsilon_t = 1$, $R_{disk} = 250(m_*/1 M_\odot)$ au, and $H_{disk} = 0.1 R_{disk}$ with our fixed $\tau_{acc}$, to obtain an estimate for the maximal turbulent disk velocity $v_t$, over time.

The outcome of this exercise is shown in Figure 8, where we find a median $\alpha_{SS} \sim 3 \times 10^{-3}$ for $\sim 3$ Myr old disks. We can compare this to the $\tau_{acc}$ distribution (Figure 5), assuming that

$$ \tau_{acc} \approx 3.8 \left( \frac{\alpha_{SS}}{10^{-3}} \right)^{-1} \left( \frac{H_s/R_s}{0.1} \right)^{-2} \times \left( \frac{m_*}{1 M_\odot} \right)^{-1/2} \left( \frac{R_s}{50 \text{ au}} \right)^{3/2} \text{ Myr,} \quad (18) $$

where we define $H_s$ as the scale height at the viscous scale radius $R_s$, which is some fraction of the outer disk radius. Given that the median $\tau_{acc}$ in our model at 3 Myr is approximately 1 Myr, the median $\alpha_{SS}$ is approximately sufficient to sustain accretion self-consistently in our model by direct kinetic energy injection alone (depending on disk geometry and injection efficiency). In Figure 8 we also compare to observational constraints for turbulence in the outer regions of some disks. With the exception of DM Tau, which has turbulent velocity $v_{t,disk} \sim 0.25c_s - 0.33c_s$, and IM Lup, which has $v_{t,disk} \sim 0.4c_s - 0.6c_s$ (Paneque-Carreño et al. 2024), attempts to measure turbulence in protoplanetary disks have resulted in upper limits. For example, for the disk around HD 163296 this upper limit is $v_{t,disk} \lesssim 0.05c_s$ (Flaherty et al. 2015, 2017), and for TW Hydra it is $v_{t,disk} \lesssim 0.08c_s$ (Flaherty et al. 2018; Teague

et al. 2018). For $M_{disk}/m_* = 0.01$, $H_{disk}/R_{disk} = 0.1$ and $\epsilon_t = 1$, and then even substituting a modest $\dot{M}_{BHL} \sim 10^{-11} M_\odot \text{yr}^{-1}$ yields $v_{t,disk} \sim 0.03c_s$, similar to the empirical upper limit on turbulence in some disks. This is approximately the median BHL accretion rate for a solar-mass star in our model at an age of $\sim 3$ Myr (see the bottom panel of Figure 2). Typically values of $M_{disk}$ are even smaller than $10^{-2}m_*$, so even greater turbulence may be imparted by the ISM. Broadly, for disks older than $\sim 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 disks, 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 disks from the ISM to sustain the observational constraints on turbulence in the outer disk.

### 4.1.2. Stellar Accretion in Regions of Low Gas Density

The second objection raised by Hartmann et al. (2006) is that stars in ionized H II regions still accrete. The typical density in such regions is far lower, and the sound speed and/or velocity dispersion are much larger, such that the BHL accretion rate may be suppressed by several orders of magnitude. The examples quoted by Hartmann et al. (2006) are some stars in Trumpler 37 that were reported by Sicilia-Aguilar et al. (2005) to be accreting substantially at a rate $\sim 10^{-8} M_\odot \text{yr}^{-1}$. This is despite a sound speed $c_{s,ext} \sim 10 \text{ km s}^{-1}$ (which we now distinguish from the sound speed in the disk $c_{s,disk}$) and an electron density $N_e \sim 3-15 \text{ cm}^{-3}$ or $\rho_{gas} \sim 2 \times 10^{-23} \text{ g cm}^{-3}$. The corresponding BHL accretion rate would therefore be $\dot{M}_{BHL} \sim 10^{-13} M_\odot \text{yr}^{-1}$ for a solar-mass star, far too small to drive observed accretion rates. This would appear to rule out infalling gas as the origin of angular momentum transport in at
least these disks. However, this conclusion requires further investigation for the following reasons.

First, as discussed above, turbulence injected into the outer disk may only serve to feed internally mediated transport in the inner disk. Therefore, turning off the supply from the reservoir of material accreted from the ISM should not immediately end accretion onto the star.

Second, the small BHL radius resulting from large $c_{s,ext}$ means that, at least initially after the gas is ionized, accretion onto the disk will be facilitated by sweep-up and not BHL accretion (see Appendix G.3). If we recalculate the sweep-up rate for a disk outer radius 250 au (expected for a solar-mass star), this increases the gas capture rate by two orders of magnitude to $\dot{M}_{su} \sim 10^{-11} M_{\odot} \text{yr}^{-1}$. While this remains much smaller than the accretion rates reported by Sicilia-Aguilar et al. (2005) in Trumpler 37, if the H II region is young ($\lesssim 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.

Finally, as the mode of mass accretion changes, so does the rate of turbulent energy injection. This is because the hot turbulent medium will typically have larger relative velocity $\Delta v_{gas} \sim c_{s,ext}$ and the material also injects thermal energy. In addition, if turbulent fluctuations are suppressed in the ionized medium, then accreted material contains little angular momentum and can fall directly into the inner disk regions. To estimate the resultant turbulent velocity within the disk, we can assume that $v_{in} = \zeta c_{s,ext}$, where $\zeta$ is a factor of order unity. Then,

$$ \dot{E}_{in} = \frac{1}{2} \dot{M}_{su} \zeta^2 c_{s,ext}^2, \quad (19) $$

where $\dot{M}_{su}$ is the sweep-up rate. The resultant steady-state Mach number is

$$ \frac{v_{t,disk}}{c_{s,disk}} \approx 0.016 \epsilon_t^{1/3} \left(\frac{H_{disk}}{10 \text{ au}}\right)^{1/3} \left(\frac{M_{disk}}{10^{-2} M_{\odot}}\right)^{-1/3} \times \left(\frac{\dot{M}_{su}}{10^{-11} M_{\odot} \text{yr}^{-1}}\right)^{1/3} \left(\frac{\zeta c_{s,ext}}{10 \text{ km s}^{-1}}\right)^{2/3} \left(\frac{c_{s,disk}}{0.2 \text{ km s}^{-1}}\right)^{-1}. \quad (20) $$

While our chosen normalization yields $\alpha_{SS} \sim 10^{-4}$, some combination of a factor three increase in the infall velocity or a factor eight increase in the sweep-up accretion rate would be sufficient to drive more substantial $\alpha_{SS} \sim 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 disk for long in an ionized 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 disks. We consider the interesting case of DM Tau, which is a star of mass $\sim 0.4-0.5 M_{\odot}$ hosting a disk with a dust cavity and accreting at a rate $\sim 10^{-8} M_{\odot} \text{yr}^{-1}$ (Manara et al. 2014; Francis & van der Marel 2020). It has a turbulent velocity $v_{t,disk} \sim 0.25c_s - 0.33c_s$,

or viscous $\alpha_{SS} \sim 0.1$ (Flaherty et al. 2020). For this disk, a recent and relatively rapid period of accretion would be needed to drive the observed turbulence. Estimating $M_{disk}/m_* \approx 0.03$ (Andrews et al. 2013), this would imply recent accretion from the ISM at a rate $\dot{M}_{BHL} \gtrsim 2 \times 10^{-8} M_{\odot} \text{yr}^{-1}$ that ended in the past $\sim 10^4$ yr. This is similar to the current stellar accretion rate (Manara et al. 2014), suggesting that the majority of the mass of the present-day disk could have originated from this putative infall 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 infall. This underlines the importance of a systematic search for observational signatures of late-stage infall events.

### 4.1.4. To What Degree Does Infall Mediate Angular Momentum Transport?

In this section, we have considered whether infall can drive the separate but related phenomena of turbulence and stellar accretion in protoplanetary disks. It is clear that there is sufficient energy in the infalling material to drive turbulence compatible with what is observed in the surface layers of the outer disk in the few cases where such constraints exist. This is not the same as the statement that this infall *does* drive observed turbulence; to conclude as such will require careful case studies of observational examples. Similarly, understanding whether infall can regulate stellar accretion, either by angular momentum extraction or by turbulent energy injection, requires more detailed theoretical and observational analysis. Internal processes are probably responsible for angular momentum transport in at least some regions of the disk, but infall could play a role in regulating this process—for example, by supplying material from the outer disk. We conclude that, although it is far from clear that these fundamental problems in disk evolution are solved by infall, ruling out its importance is not currently justified.

## 4.2. Appearance of BHL-accreting Disks

A pertinent question is what we expect the appearance of a disk undergoing BHL accretion to be. As previously discussed, when accreting material with a high vorticity, the resultant disk 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 disk 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” disks, or misalignment between the inner and outer regions (Kuffmeier et al. 2021). By contrast, at least the millimeter dust component of disks is frequently observed to be extremely settled (Duchêne et al. 2003; Pinte et al. 2007, 2016; Villenave et al. 2020, 2022; Pizzati et al. 2023). However, (sub)millimeter emission is also commonly more centrally concentrated than the gaseous disk by a factor $\sim 2.5$ (Ansdell et al. 2018; Facchini et al. 2019; Andrews 2020), which may be due to radial drift or optical effects (Rosotti et al. 2019; Trapman et al. 2019). It is not clear whether turbulence driven in the outer disk should be reflected in the dust distribution, which may have drifted inward and/or settled vertically, even in the presence of the warping induced by infall (Aly et al. 2024). In fact, we can consider one of the best-known disks, HL Tau (ALMA Partnership et al. 2015), which
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The column number density of material assuming a mean molecular mass $\mu = 2.3$ and $\epsilon_{geo} = 1$ following Equation (22) for the surviving disks in our fiducial model. The colors of the lines represent the age of the disks. 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).</figcaption></figure>

has an extremely settled dust disk, evidenced by both dust height and polarization that suggest low $\alpha_{SS} \sim 10^{-5}$ (Pinte et al. 2016; Ueda et al. 2021). It is also one of the earliest known examples of a disk being fed by infalling gas (Hayashi et al. 1993; Welch et al. 2000; Yen et al. 2019; Garufi et al. 2022) at a substantial rate ($\gtrsim 5 \times 10^{-9} M_{\odot} \text{ yr}^{-1}$; Gupta et al. 2024). This example therefore makes a compelling argument that the millimeter dust component of the disk is not necessarily strongly influenced by infall. It also suggests that turbulence may be substantially lower in the midplane with respect to the surface layers. This represents an important topic for future study.

### 4.3. Observability of Infall 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 infall events from the data that are currently available. In fact, growing evidence suggests that a considerable fraction of disks are undergoing some kind of interaction with their environment. For example, from their recent scattered light survey of near-IR polarimetric imaging for 43 disks in Taurus, Garufi et al. (2024) reported that 16% of those disks showed evidence of interaction with ambient material (although this sample is not necessarily representative). While inferring observational signatures of late-stage infall 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 disks in our model. For accretion onto a disk from a radius $R_{acc}$, the maximum freefall time onto a star of mass $m_{*}$ is

$$ \tau_{\text{ff,max}} \approx 2 \times 10^3 \left( \frac{R_{\text{acc}}}{250 \text{ au}} \right)^{3/2} \left( \frac{m_{*}}{1 M_{\odot}} \right)^{-1/2} \text{ yr.} \quad (21) $$

Given a BHL accretion rate $\dot{M}_{BHL}$, we can then estimate the column density of material external to the disk:

$$ N_{\text{ext}} \sim 2 \times 10^{19} \epsilon_{\text{geo}}^{-1} \left(\frac{\mu}{2.3}\right)^{-1} \left(\frac{\dot{M}_{\text{BHL}}}{10^{-9} M_{\odot} \text{ yr}^{-1}}\right)^{-1/2} \times \left(\frac{R_{\text{acc}}}{250 \text{ au}}\right)^{-1/2} \left(\frac{m_{*}}{1 M_{\odot}}\right)^{-1/2} \text{cm}^{-2}, \quad (22) $$

where $\epsilon_{geo} \le 1$ is a geometric factor accounting for the nonuniform distribution of material and $\mu$ is the mean molecular mass. We will assume $\epsilon_{geo} = 1$ here, although, based on the diverse ambient structures seen in observations, a wide range of $\epsilon_{geo}$ may be plausible (e.g., Garufi et al. 2024; Zurlo et al. 2024).

The typical column density implied by Equation (22) is relatively low, such that we would expect that infall 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 disks with ALMA. These isotopologues have a self-shielding column density of $\sim 10^{15} \text{cm}^{-2}$ (van Dishoeck & Black 1988; Visser et al. 2009) and typical abundance relative to molecular hydrogen $\sim 10^{-5}$ (Trapman et al. 2023). As an order-of-magnitude estimate, we may therefore expect column densities $\lesssim 10^{20} \text{cm}^{-2}$ to have the CO photodissociated. From Equation (22), we plot the $N_{\text{ext}}$ distribution over time resulting from the accretion rates in our model in Figure 9. We see that immediately infalling CO would be photodissociated for the majority of disks in the age range $\sim 1-3 \text{ Myr}$. Interestingly, in this age range, the fraction of disks above the threshold $N_{\text{ext}} = 10^{20} \text{cm}^{-2}$, approximately the CO self-shielding column, drops from approximately 30% to close to 0%. Given that this is the range of ages in Taurus, this would be commensurate with the inferred $\sim 16\%$ of disks 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 disks. Such data could be applied to make statistical inferences about the prevalence of late-stage infall, correlating it with disk properties and testing the theoretical predictions we have made in this work.

### 4.4. Implications

Our results strongly suggest that disk 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 summarize some examples as follows:

1. *Disk substructures:* A diverse range of substructures have been observed in protoplanetary disks (Andrews 2020). We have demonstrated that the disk is undergoing constant perturbation by infalling material. We would therefore expect some disks to exhibit structure owing to this process (Kuffmeier et al. 2020). For example, by accreting gas from the ISM with a misaligned angular momentum vector, we expect disk 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 disk warping are evident in some disks (Walsh et al. 2017), and it seems probable that they originate from late-stage infall at least in some instances.

2. *Accretion outbursts and short accretion times:* FU Orionis events are short accretion bursts that last ~100 yr (Hartmann & Kenyon 1996), named for the archetypal example (Herbig 1966, 1977). Infall 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 whether infall events could produce short rise times, but they may indirectly contribute by instigating instabilities in the disk (e.g., Deng et al. 2020; Speedie et al. 2023). Extended nebulosity surrounding several stars undergoing outburst supports the role of infall in at least some cases (Hales et al. 2024; Zurlo et al. 2024). Outbursting star-disk systems may also be more broadly categorized as a subset of systems with short accretion timescales (Manara et al. 2023; Almendros-Abad et al. 2024). Among such disks, a lifetime problem appears particularly pronounced for disks around some Herbig stars (Grant et al. 2023). The assumptions of our model mean that we cannot directly explain these short accretion timescales; we *impose* this timescale for our model. Nonetheless, given a fixed timescale $\tau_{\text{acc}}$, disks with short $\tau_{\text{acc}}$ may be episodically replenished in our model. Thus, to a degree the inferred accretion timescale is decoupled from the expected frequency of observed disks, which no longer needs to be ~$\tau_{\text{acc}}/\tau_{\text{disk}}$, for disk age $\tau_{\text{disk}}$. Instead, the frequency of disks with a given $\tau_{\text{acc}}$ would be dictated by the frequency of infall events that lead to rapid mass gain. Indeed, in our model a small number of disks with very

short accretion timescales (~0.01 Myr) are retained for ~3 Myr (Figure 5).

3. *Long-lived disks:* Stochastic replenishment of disks represents a simple explanation for so-called “Peter Pan” disks, which are (low-mass) stars that are tens of Myr old but still exhibit a gas-rich accretion disk (Lee et al. 2020; Silverberg et al. 2020; Galli et al. 2021b). It is not clear whether Peter Pan disks are exclusively around low-mass stars, or whether finding disks 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 disks appear across a number of regions of different ages, even where the majority of disks have been depleted (Ansdell et al. 2020). Late-stage infall seems to be a viable way to explain some of these long-lived disks. This further motivates deep observational programs to characterize the super-disk-scale gas surrounding protoplanetary disks.

4. *Planet properties:* Possibly the most exciting consequence of protoplanetary disks 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 disks 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 the rate at which protoplanetary disks capture material from their surroundings via 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 disk. Our three most important conclusions from this experiment are as follows:

1. Disk masses, radii, and stellar accretion rates as a function of stellar mass and age can result directly from BHL accretion.

2. Material in the disk is constantly replenished over its lifetime. We estimate that 20%–70% of disks 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 disk.

3. The energy influx due to infalling material onto the disk is sufficient to drive empirically inferred protoplanetary disk turbulence. If this infall produces turbulence in the outer disk, then we expect a broad distribution for $\alpha_{\text{SS}}$, with typical value $\alpha_{\text{SS}} \sim 10^{-3}$, while some disks with much higher $\alpha_{\text{SS}} \sim 0.1$ (as for DM Tau) and lower
$\alpha_{SS} \sim 10^{-5}$ are expected. We highlight that stellar accretion is not necessarily driven by turbulence in the outer disk, and whether infall can contribute to global angular momentum transport to the inner disk remains an open question.

Our approach is highly simplified, representing a novel method that can be built on and adapted to semianalytic models in order to include the role of BHL accretion in planet formation. Several processes may alter the rate of replenishment in the disk, but in particular stellar feedback may curtail late-stage accretion for relatively massive stars ($\gtrsim 2 M_\odot$) and stars born in regions of many OB stars. These feedback processes must be considered in future works.

If the ISM plays an important role in disk evolution, this represents a considerable change of direction for protoplanetary disk theory, which has largely focused on processes operating on an isolated star-disk system. Consequences of this new picture of disk 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 infall during formation influences protoplanetary disk structures, angular momentum transport, and the observed exoplanet population. In particular, deep observational surveys aiming to correlate extended (infalling) structures with disk properties and/or stellar accretion rates would test our findings and are an urgent goal for the future.

## Acknowledgments

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. A.J.W. has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-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 program (PROTOPLANETS, grant agreement No. 101002188).

## Appendix A

### Modeling 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 $\rho_0 = \Sigma_0/2h$ (where we hereafter fix $\Sigma_0 = 12 M_\odot \text{ pc}^{-2}$) and Mach number $M \gg 1$ on a spatial scale $R \sim 1/k$, the volume-weighted distribution of gas density is well described by a lognormal distribution:

$$dp(\delta | k) = \frac{1}{\sigma_k \sqrt{2\pi}} \exp\left(-\frac{\delta^2}{2\sigma_k^2}\right) d\delta \quad (A1)$$

<figure><img src="image_13.png" /><figcaption>Figure A1. The fractional change in the variance $\Delta S/S$ of the Gaussian random field at each spatial scale with varying spatial resolution $\Delta R/R$. We mark the scale height $h$ as a vertical dotted line. The locations at which each $\Delta S$ is computed are shown as scatter points.</figcaption></figure>

$$\delta \equiv \ln\left(\frac{\rho}{\rho_0}\right) - \left\langle\ln\left(\frac{\rho}{\rho_0}\right)\right\rangle, \quad (A2)$$

where

$$\left\langle \ln \left( \frac{\rho}{\rho_0} \right) \right\rangle = - \frac{\sigma_k^2}{2}. \quad (A3)$$

The density dispersion for $k$ can be estimated as

$$\sigma_k \approx \left( \ln \left[ 1 + \frac{3}{4} \frac{\langle v_t^2(k) \rangle}{c_s^2 + \kappa^2 k^{-2}} \right] \right)^{1/2}, \quad (A4)$$

where we will assume throughout this work that the sound speed $c_s = 0.2 \text{ km s}^{-1}$. Here $v_t$ is the turbulent velocity, for which the average $\langle v_t^2(k) \rangle \approx \sigma_g^2(k) \sim kE(k) \propto k^{1-p}$, where $E(k)$ is 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 $\langle v_t^2(k) \rangle^{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 to the square of the Mach number. The additional term in the denominator is proportional to the square of the galactic epicyclic frequency $\kappa$, for which we will assume a flat rotation curve with $\kappa = \sqrt{2}\Omega$, adopting the orbital frequency $\Omega = 2.6 \times 10^{-2} \text{ Myr}^{-1}$ appropriate for the solar neighborhood. This modification is necessary because at large $R$ the differential rotation $\kappa/k$ plays a similar role to thermal pressure on small scales.

Assuming that the turbulent properties of the medium depend only on local gas properties at that scale, the distribution of densities averaged over any spatial scale is also
lognormal with variance:

$$\sigma^2(R) = \int d\ln(k) \sigma_k^2(k) |\tilde{W}(k, R)|^2, \qquad (A5)$$

where $\tilde{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 $\tilde{W}(k, R)$ to be a Fourier-space top-hat function: $\tilde{W}(k | R) = 1$ if $k \le R^{-1}$ and $\tilde{W}(k | R) = 0$ if $k > R^{-1}$.

With 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 $\delta_1$ to $\delta_2 \equiv \delta_1 + \Delta\delta$ given a change from scale $k_1$ to $k_2$ is given by

$$p(\delta_1 + \Delta \delta) d \Delta \delta = \frac{1}{\sqrt{2\pi \Delta S}} \exp\left(-\frac{(\Delta \delta)^2}{2 \Delta S}\right) d(\Delta \delta) \qquad (A6)$$

$$\Delta S \equiv S_2 - S_1 \equiv \sigma^2(R_2) - \sigma^2(R_1), \qquad (A7)$$

where we define the variance

$$S(R) \equiv \sigma^2(R). \qquad (A8)$$

With 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 \to \infty$, we have $\delta \to 0$ and $\sigma(R) \to 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$, $\Delta R_i$. A single trajectory is then

$$\delta(R_i) \equiv \sum_{j}^{R_i > R_j} \Delta \delta_j, \qquad (A9)$$

where each $\Delta \delta_j$ is drawn from the probability density function defined by Equation (A6). We will hereafter define a grid with a maximum scale $R_{\max} = 5h$ down to a minimum scale $R_{\min} < 0.01$ pc, each successive turbulent scale in the trajectory being 95% 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 ($\Delta S/S < 0.1$). In other words, any fluctuations in density or velocity on scales between those for which we calculate the structure of the medium will not contribute significantly to the local field. If, instead, we use $\Delta R/R \lesssim 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 \sim h$ (see red lines in Figure A1).

We have defined $S(R)$ as the variance in the logarithmic density, but the same argument applies to any property that 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 scale. The only difference is that the variance is determined directly by the turbulent spectrum (in our case, $\sigma_g^2(R) \propto 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 disk of finite thickness $h$ (e.g., Begelman & Shlosman 2009), it can be shown that unstable collapse occurs if

$$\frac{\rho}{\rho_0} \ge \frac{\rho_c}{\rho_0} = \frac{Q_0(h)}{2 \tilde{\kappa}} (1 + \tilde{k}) \left[ \frac{\sigma_g(k)^2}{\sigma_g(h^{-1})^2} \tilde{k} + \tilde{\kappa}^2 \tilde{k}^{-1} \right], \qquad (A10)$$

where we have defined the dimensionless versions of the epicyclic frequency $\tilde{\kappa} \equiv \kappa/\Omega = \sqrt{2}$ and wavenumber $\tilde{k} \equiv |k| h$. For a disk in vertical equilibrium, $\sigma_g(h^{-1}) = h\Omega$, and we will fix $\sigma_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-disk scale height (e.g., Vieira et al. 2023). Then, at any given scale, we can define an analogous critical logarithmic density perturbation:

$$\delta_c(R) = \ln\left(\frac{\rho_c}{\rho_0}\right) - \left\langle \ln\left(\frac{\rho}{\rho_0}\right) \right\rangle. \qquad (A11)$$

The unstable, collapsing star-forming region is defined at the largest scale $R_{SFR}$ for which $\delta \ge \delta_c$. Interpolating to the cylindrical geometry, if the $R_{SFR} \gtrsim h$, the gas mass of the collapsing cloud is defined as

$$M(\rho | R) \equiv 4 \pi \rho(R) h^3 \left[ \frac{R^2}{2h^2} + \left(1 + \frac{R}{h}\right) \exp\left(-\frac{R}{h}\right) - 1 \right], \qquad (A12)$$

where $\rho = \rho_c$ and $R = R_{SFR}$.

## A.3. Evolution of the Turbulent Medium

We can compute the temporal evolution of the ISM by assuming that the turbulence is globally in a steady state, such that the turbulent velocity cascade is maintained outside of collapsing regions. The probability density function for any Gaussian random variable $\delta$ under these conditions obeys a generalized Fokker-Planck equation. We can therefore redefine our trajectory at $t + \Delta t$ by computing a random walk step:

$$\Delta \delta_j(t + \Delta t) = \Delta \delta_j(t) \exp(-\Delta t / \tau_t) + \mathcal{R} \sqrt{\Delta S(1 - \exp(-2 \Delta t / \tau_t))}, \qquad (A13)$$

where $\mathcal{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% of the shortest turbulent timescale in our grid. Equation (A13) applies both to the logarithmic density perturbations and to each component of the interstellar gas velocity.

# Appendix 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
<figure><img src="image_14.png" /><figcaption>**Figure B1.** Evolution of the mean density (black lines; left axis) and fraction of star formation in three star-forming regions (red lines; right axis) of different mass in our model. We show the results for three different initial masses of the star-forming regions, $M_{SFR}$.</figcaption></figure>

follow the collapse of the GMC and formation of stars. Once the GMC becomes gravitationally unstable, it collapses on a freefall timescale. It does not do so monolithically, but because smaller subregions of the cloud may be higher in density (and therefore have shorter freefall timescales), 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 semianalytic approaches have been developed to describe this process (see Vázquez-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.

We adopt the approach of Girichidis et al. (2014) in modeling the rate of star formation during the collapse of the GMC. The approach of those authors is attractive both for its simplicity and for 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 freefall collapse from rest, which may be contained within a larger gravitationally unstable cloud. The outer radius $r$ obeys the equation

$$ \frac{d^2 r}{dt^2} = - \frac{GM}{r^2} = - \frac{4\pi G \rho_1 r_1^3}{3r^2}. \qquad (B1) $$

Here $M = \rho V = (4\pi/3)\rho r^3$ is the fixed mass of the region with volume $V$, initial density $\rho_1 = \rho(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

$$ \theta = t / \tau_{ff}, \qquad (B2) $$

where

$$ \tau_{ff} = \sqrt{\frac{3\pi}{32G\rho_1}} = (\gamma\rho_1)^{-1/2} \qquad (B3) $$

is the freefall time, and the second equality defines $\gamma$ for convenience. While there is no exact solution to Equation (B1),

Girichidis et al. (2014) point out that

$$ r(\theta) = r_1(1 - \theta^2)^{p/3}, \qquad (B4) $$

with $p=2$, is a good approximation. If the mass is constant, then the density evolves as

$$ \rho(\theta) = \rho_1(1 - \theta^2)^{-p}. \qquad (B5) $$

To 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 $\rho_a$. We will choose $\rho_a$ to yield a freefall time of 0.1 Myr, although this choice does not strongly influence the star formation rate as long as this timescale is much shorter than the global freefall timescale. The problem is now reduced to finding the initial density $\rho_*$ that evolves to $\rho_a$ at time $t$. If $p=2$, then, inverting Equation (B5),

$$ \rho_*(\rho_a, t) = \frac{2\gamma\rho_a t^2 + 1}{2\gamma^2\rho_a t^4} + \sqrt{\left(\frac{2\gamma\rho_a t^2 + 1}{2\gamma^2\rho_a t^4}\right)^2 - \frac{1}{\gamma^2 t^4}}. \qquad (B6) $$

Now we just need to calculate the fraction of the total cloud mass that has $\rho > \rho_*$. Given that we have an initial probability density function in terms of volume (as in Equation (A1)),

$$ P_V = \frac{1}{V_{tot}} \frac{dV}{d\rho}, \qquad (B7) $$

we can convert this into the mass-weighted version,

$$ P_M = \frac{1}{M_{tot}} \frac{dM}{d\rho}, \qquad (B8) $$

via the equation

$$ \frac{dV}{d\rho} = \frac{dV}{dM} \frac{dM}{d\rho} = \frac{1}{\rho} \frac{dM}{d\rho}, \qquad (B9) $$

for some total volume $V_{tot}$ or mass $M_{tot}$. Here we define $M_{tot} = M(\rho_c|R_{SFR})$ from Equation (A12). Now the total mass in stars is just

$$ M_* = M_{tot} \int_{\rho_*}^{\infty} P_M(\rho) d\rho = \frac{M_{tot}}{2} \left[ 1 + \operatorname{erf} \left( \frac{\ln \frac{\rho_c}{\rho_*} + \sigma_k^2}{\sqrt{2} \sigma_k} \right) \right], \qquad (B10) $$

where the last equality assumes an initially lognormal probability density distribution in volume. As in Section 2.2, the critical density of the region is $\rho_c(R_{SFR})$, and $\sigma_k$ is the dispersion given $k \sim 1/R_{SFR}$, for $R_{SFR}$ the largest scale of the collapsing star-forming region. We stop star formation and assume that the cloud is immediately dispersed once $M_* = \epsilon M$, where we adopt $\epsilon = 0.5$, consistent with calculations for core formation (Matzner & McKee 2000).

We 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_\odot$ 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 $\times 10^5$ yr for $M_{SFR} \sim 10 M_\odot$ to $\sim 1$ Myr at $M_{SFR} \sim 500 M_\odot$. 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 disks during the bound, star-forming phase.
## Appendix C
Sensitivity to Disk Dispersal Criterion

We have used a fiducial criterion for disk dispersal (when $M_{\text{disk}} < 3 \times 10^{-5} M_{\odot}$). While few disk-hosting stars have been found with nominal mass substantially 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 disk properties and replenishment fractions.

We show a comparison to Figure 3 with $M_{\text{disk}} > 10^{-5} M_{\odot}$ in Figure C1. There are some slight differences in the distributions of disk 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 for the very low mass disks 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 with those in Figure 7. We find a small increase in the replenishment fraction for disks, particularly at late times. This is because disks with a higher fraction of replenished material are often lower in mass. In summary, we find only minor differences in our model outcomes if we reduce the disk dispersal criterion by a factor of a few.

<figure><img src="image_15.png" /><figcaption>Figure C1. Same as Figure 3, but with a reduced dispersal disk mass threshold, $M_{\text{disk}} < 10^{-5} M_{\odot}$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure C2. Same as Figure 7, but with a reduced dispersal disk mass threshold, $M_{\text{disk}} < 10^{-5} M_{\odot}$.</figcaption></figure>

## Appendix D
Sensitivity to the Initial Disk Mass

We consider how our conclusions on disk replenishment as a function of time change when we do not assume a negligible initial disk mass, as we assumed in our fiducial model. We therefore repeat our calculations with disk masses initially similar to the observed distribution in young star-forming regions. In practice, we adopt a lognormal distribution for the disk mass with 1 dex of scatter around the relationship

$$M_{\text{disk}}(t=0) = 0.01 \left( \frac{m_{*}}{1 M_{\odot}} \right)^2 M_{\odot}. \quad (\text{D1})$$

We show the resultant disk properties in Figure D1, which again are in good agreement with the observed disk properties. Our adopted initial disk masses are approximately the

maximum we could adopt while still remaining consistent with the observed distribution of disk masses. The outcome for the half-life disk replenishment factor is shown in Figure D2. Although in this case $f_{M,1/2}$ is somewhat reduced early on in the disk evolution, overall differences in the distribution over time are minor. The fraction of disks with at least some recently accreted material remains similar, while we find that $f_{M,1/2} \sim 0.2-0.45$ throughout the disk lifetime. Therefore, even under the assumption that the primordial disk is coincidentally similar in mass to the subsequently accreted one, and ignoring all internal depletion processes, we still expect a substantial fraction of disks to be composed of recently captured disk material. Our conclusion that the present-day disk properties are not a good probe of internal disk physics therefore appears robust.
<figure><img src="image_17.png" /><figcaption>**Figure D1.** Same as Figure 3, but for an initial disk mass $M_{\text{disk}}(t=0) = 0.01(m_*/1 M_{\odot})^2 M_{\odot}$ with 1 dex of scatter.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure D2. Same as Figure 7, but for an initial disk mass $M_{\text{disk}}(t=0) = 0.01(m_*/1 M_{\odot})^2 M_{\odot}$ with 1 dex of scatter.</figcaption></figure>

## Appendix E
### Stellar Feedback and BHL Accretion

#### 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 2003, 2004). For massive stars, when dust 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 (Appendix 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 Appendix 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 the evolution model that we use to compute the evolution of the stellar luminosity.

#### 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;<fnref n="4" /> Paxton et al. 2011; Choi et al. 2016; Dotter 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 (2003). 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

$$L_{\text{acc}} = \frac{Gm_*\dot{M}_{\text{acc}}}{R_*}, \qquad (\text{E1})$$

where $\dot{M}_{\text{acc}}$ is the accretion rate onto the star (not to be confused with $\dot{M}_{\text{BHL}}$) and $R_*$ is the stellar radius. We assume that the accretion flow emits as a blackbody with temperature

$$T_{\text{acc}} = \left( \frac{L_{\text{acc}}}{4\pi R_*^2 \sigma_{\text{SB}}} \right)^{1/4}. \qquad (\text{E2})$$

When we come to consider the short-wavelength ($\lambda < 2070$ Å) photons that drive stellar winds, we will integrate over the

<footnote n="4">⁴ https://waps.cfa.harvard.edu/MIST/</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure E1.** The thermal radius $R_{\text{therm}}$, defined to be the radius at which the bolometric stellar luminosity heats the gas sufficiently to 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 colored by the age of the star.</figcaption></figure>

stellar spectrum:

$$L_{*,<2070} = \int_{0}^{2070} L_{*,\lambda} d\lambda, \qquad (E3)$$

where $L_{*,\lambda}$ is the wavelength-specific stellar luminosity.

### E.3. Thermal Radius

We now consider whether, given the thermal structure of the flow imposed by stellar irradiation, the temperature is ever sufficient to raise the sound speed $c_s > c_{s,crit} \sim v_{\text{esc}}$, where $v_{\text{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_{\text{acc}}$. The problem is then very similar to a passively irradiated protoplanetary disk (Chiang & Goldreich 1997). The surface of the column is heated at a rate $Q_{\text{irr}}$, while the disk cools at a rate $Q_{\text{cool}}$ at the surface of the accretion column. The different geometry of the accretion flow (e.g., a cone that widens far from the star) compared to the disk just means that both the heating and cooling areas at a given radius are increased by the same factor. Then, the center of the column has temperature $T_{\text{col}}$ such that

$$\sigma_{\text{SB}} T_{\text{col}}^4 = \frac{\phi L_{*}}{8 \pi R_{\text{acc}}^2}. \qquad (E4)$$

Here $L_{*}$ is the bolometric luminosity, and $\phi$ in the context of a protoplanetary disk is the flaring angle—i.e., the angle at which

the stellar radiation meets the surface structure. The largest value this can be is $\phi = \pi/2$, which would correspond to a column where the ISM interior to $R_{\text{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.

It is clear from the above discussion that $c_s/c_{s,crit} \propto R_{\text{acc}}^{1/4}$, 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 that the central star does not set the ISM temperature). Therefore, the accretion radius $R_{\text{acc}}$ can always therefore be limited by $R_{\text{therm}}$. The outcomes of these calculations are shown in Figure E1. We find that $R_{\text{therm}}$ is typically much smaller for low-mass stars but that even for $m_{*} = 0.1 M_{\odot}$ we have $R_{\text{therm}} \sim 5 \times 10^4$ au $\sim R_T$, the tidal radius. This means that $R_{\text{therm}} \gtrsim R_T$ in most scenarios, and we do not expect the passive stellar irradiation to limit BHL accretion.

### E.4. Internally Driven Wind

#### E.4.1. Wind Radius

A second feedback effect is the photoevaporative wind that can be driven by ionizing 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 $\dot{M}_{\text{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
<figure><img src="image_20.png" /><figcaption>**Figure 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 $\lambda < 2070$ Å at a fixed ambient gas density $\rho = 10^{-20}$ g cm$^{-3}$. Points represent the outcome of stellar models with no stellar accretion (circles), accretion rates following $\dot{M}_{acc} = 10^{-8}(m_{*}/1 M_{\odot})^2 M_{\odot} yr^{-1}$ (triangles), and accretion rates following $\dot{M}_{acc} = 10^{-7}(m_{*}/1 M_{\odot})^2 M_{\odot} yr^{-1}$ (squares). These points are colored by the age of the star.</figcaption></figure>

to the momentum injected from within. A proper prescription of the role of feedback on the BHL accretion flow for low-mass (\lesssim 3 $M_{\odot}$) stars will require a devoted study. We nonetheless present a crude approximation as follows.

We base our estimate on injecting a fixed fraction $\epsilon_{wind}$ of energy from high-energy photons (wavelength $\lambda < 2070$ Å) 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} \sim v_{esc} = \sqrt{2Gm_{*}/R_{acc}}$. In this case,

$$ \begin{aligned} \frac{1}{2} \dot{M}_{\text{wind}} v_{\text{launch}}^2 &= \epsilon_{\text{wind}} L_{*} < 2070 \Leftrightarrow \dot{M}_{\text{wind}} \\ &= \epsilon_{\text{wind}} \frac{L_{*} < 2070 R_{\text{acc}}}{Gm_{*}}. \end{aligned} \quad (E5) $$

BHL accretion can only proceed if $\dot{M}_{BHL} > \dot{M}_{wind}$, or

$$ R_{\text{acc}} \gtrsim (2Gm_{*})^{-3} \left( \frac{\epsilon_{\text{wind}} L_{*} < 2070}{\pi \rho} \right)^2 = R_{\text{wind}}. \quad (E6) $$

The factor $\epsilon_{wind}$ folds in the fractional solid angle covered by the accretion flow and the complex wind physics. Hereafter we will assume $\epsilon_{wind} = 0.1$. We have also assumed that the wind is launched at the minimum possible velocity $v_{esc}$, so that 10% 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

maximizing replenishment. In practice, many photons will escape along low-density lines of sight (reducing $\epsilon_{wind}$), and higher-energy photons may also launch more rapid winds with velocity $\sim 10$ km s$^{-1}$. In making this choice, we therefore explore the opposite extreme compared to the “no wind” scenario of our fiducial model. Examples of the $R_{wind}$ values obtained from this calculation for stars with varying masses and accretion rates are shown in Figure E2.

### E.4.2. BHL Accretion Rates with Wind Suppression

In order to estimate the role of internally driven photo-evaporative winds in limiting the BHL accretion flow, we make the assumption that no accretion onto the disk can proceed if $R_{acc} < R_{wind}$. We assume $\dot{M}_{BHL} = 0$ if this condition is met. In practice, for efficiency we precompute a table of $R_{wind}$ in stellar age-mass space and then interpolate and rescale 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 $\dot{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 disk radius, by the sweep-up of gas (see Appendix G.3). However, quantifying this process in a star-forming region that may be partially ionized on large
<figure><img src="image_21.png" /><figcaption>**Figure E3.** Same as Figure 2, but including a prescription for the suppression of the BHL accretion by internally driven photoevaporative winds.</figcaption></figure>

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.

For 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 timescale 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 disk evolution calculations to understand the degree to which BHL shapes the disk properties over its lifetime.

### E.4.3. Disk Evolution Including Winds

We calculate the disk 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_{\text{wind}}^2/R_{\text{acc}}^2)$, to avoid discontinuities in the second derivative. The disk fraction resulting from our calculation, analogous to that in Figure 6, is shown in Figure E4. Predictably, the overall disk fraction is somewhat reduced by a moderate factor. However, the influence of winds depends strongly on stellar mass. Dividing the disk lifetimes into different mass bins, we find that disks with host stars of mass $0.3 M_{\odot} < m_{*} < 1 M_{\odot}$, which make up a large fraction of the observed samples, show a similar disk fraction evolution to the full sample shown in Figure 6. Disks around higher-mass stars are depleted more rapidly, since BHL accretion is suppressed at a lower ambient density. Stars with $m_{*} < 0.3 M_{\odot}$ are again depleted more rapidly, but this is due to our fixed mass threshold for a disk to be considered dispersed. Low-mass stars host lower-mass disks, with masses that are more frequently below this threshold.

The disk 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 disks (which have typically undergone faster BHL accretion) are also similar to the negligible wind case. However, we can ask whether these disks 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 ($\gtrsim 8$ Myr), only ~10% of disks have accreted more than half of their mass recently. Nonetheless, based on our simple estimate, a substantial fraction of disks (~20%–60% at ages ~1–3 Myr) have undergone recent replenishment when photoevaporative winds are considered.
<figure><img src="image_22.png" /><figcaption>**Figure E4.** Same as Figure 6, but including a prescription for the suppression of the BHL accretion by internally driven photoevaporative winds. We have also separated the disk survival fraction into different mass bins, shown by different line colors. Black represents all stars in our sample with $m_{*} > 0.1 M_{\odot}$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure E5.** Same as Figure 3, but including a prescription for the suppression of the BHL accretion by internally driven photoevaporative winds.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure E6.** Same as Figure 7, but including a prescription for the suppression of the BHL accretion by internally driven photoevaporative winds.</figcaption></figure>

## Appendix 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_{\odot}$. We therefore repeat the disk evolution experiment, but with a lower cutoff of $2 M_{\odot}$ in the masses of the star-forming region, below which there is not sufficient mass to produce a solar-mass star at a star formation efficiency of 50%.

We show the outcome of our calculations in terms of the disk property distribution in Figure F1 and the replenishment fraction in Figure F2. The distribution of properties in the model remains 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 disks with a substantial quantity of replenished material remains comparable to our fiducial model. We conclude that our findings are not strongly sensitive to our choice of minimum mass of the star-forming regions.
<figure><img src="image_25.png" /><figcaption>Figure F1. Same as Figure 3, but with a lower mass limit of the collapsing star-forming region of 2 M<sub>☉</sub>.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure F2. Same as Figure 7, but with a lower mass limit of the collapsing star-forming region of 2 M⊙.</figcaption></figure>

## Appendix 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 we explore some open questions about the process of capture of interstellar material.

#### G.1. BHL Accretion Mediated by a Disk

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., $\Delta v \sim c_s$), which is rarely the case in our numerical calculations. In a magnetized medium, the BHL accretion rate may also be reduced by a factor $\sim 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 disk 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 nonradial 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 and star–star cases, capture requires a periastron $a_p$ during an encounter between a cloud with radius $R_{cloud}$ and mass $M_{cloud}$ and a star with mass $m_*$ such that (Fabian et al. 1975)

$$a_p \lesssim a_{\text{capt}} \approx R_{\text{cloud}} \left[ \frac{Gm_*}{R_{\text{cloud}} v_{\infty}^2} q_{\text{cloud}}^{-1} (1 + q_{\text{cloud}}) \right]^{1/6} . \quad (G1)$$

Here $q_{cloud} = M_{cloud}/m_*$, and $v_{\infty}$ is the relative velocity between the star and cloud at infinite separation. Then, if $a_{capt} > R_{cloud}$, the cross section is

$$\sigma_{\text{capt}} \approx \frac{2\pi Gm_*(1 + q_{\text{cloud}})(a_{\text{capt}} - R_{\text{cloud}})}{v_{\infty}^2}. \quad (G2)$$

We 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.

The 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} \sim R_{BHL}$. Therefore, we can estimate the total rate of capture for a cloud
<figure><img src="image_27.png" /><figcaption>**Figure G1.** The expected mass accretion rate $\dot{M}_{\text{capt}}$ onto a disk around a star of mass $m_* = 1 M_{\odot}$ due to gravitationally focused cloudlet capture in a gaseous medium of density $\rho = 100 M_{\odot} \text{ pc}^{-3} = 6.77 \times 10^{-21} \text{ g cm}^{-3}$. The color bar shows the averaged mass accretion rate, and the rate $\Gamma_{\text{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_{\text{cloud}} \approx 0.1 R_{\text{BHL}}$ (dashed cyan line).</figcaption></figure>

of radius $R_{\text{cloud}}$ as

$$ \Gamma_{\text{capt}}(R_{\text{cloud}}|\rho_{\text{cloud}}, \sigma_v, \nu) = \int_{0}^{\nu_{\text{BHL}}} n_{\text{clouds}}(R_{\text{cloud}}) \sigma_{\text{capt}} \tilde{v} g_3(\tilde{v}|\sigma_v) d\tilde{v}, \quad (G3) $$

where

$$ g_{N_D}(\Delta v|\sigma_v) = \frac{2\Delta v^{N_D-1}}{\Gamma(N_D/2)} \left(\frac{1}{4\sigma_v^2}\right)^{N_D/2} \exp\left(-\frac{\Delta v^2}{4\sigma_v^2}\right) \quad (G4) $$

is the Maxwell–Boltzmann distributions in $N_D$ dimensions. We have cast Equation (G3) in terms of the cloud density $\rho_{\text{cloud}} = 3M_{\text{cloud}}/4\pi R_{\text{cloud}}^3$. If large-scale turbulence dominates density fluctuations, we can estimate that a typical cloud density will be $\rho_{\text{cloud}} \sim \rho_{\text{gas}}$. We have introduced $n_{\text{clouds}} = 1/(2R_{\text{cloud}})^3$, where $2R_{\text{cloud}}$ is the wavelength of the turbulent oscillations in the ISM that produce a cloud of radius $R_{\text{cloud}}$. Finally, we have defined the BHL velocity $v_{\text{BHL}}$, which is the velocity that yields a particular $R_{\text{BHL}} = R_{\text{cloud}}$ when substituted into Equation (3). The overall accretion rate is then $\dot{M}_{\text{capt}} = \Gamma_{\text{capt}}M_{\text{cloud}}$.

In Figure G1, we show the results of these calculations for a solar-mass star in an ISM of density $\rho_{\text{gas}} = 100 M_{\odot} \text{ pc}^{-3}$. We find that capture rates are typically maximized for cloud sizes $R_{\text{cloud}} \sim 0.1 R_{\text{BHL}} - 1 R_{\text{BHL}}$. Assuming $R_{\text{acc}} = R_{\text{BHL}}$ and substituting $\Delta v_{\text{gas}} = 1 \text{ km s}^{-1}$ for a solar-mass star into Equation (4) yields $\dot{M}_{\text{BHL}} \approx 2.4 \times 10^{-8} M_{\odot} \text{ 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 disk assembly, alongside BHL accretion. This may increase the rate of interstellar gas capture by a factor of order unity, and it may be the subject of future numerical experiments.

### G.3. Face-on Accretion

A third mechanism for ISM accretion is face-on disk accretion, where an existing disk produces a large physical cross section (Moeckel & Throop 2009; Wijnen et al. 2017). In this case, the ram pressure between a preexisting disk 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 disk shrinks as a result of deposition of low angular momentum material onto the disk surface (Wijnen et al. 2017). There is no steady state for this process unless material can be constantly replenished. This replenishment of outer disk material may be viscous, but then the quantity of material accreted is limited by the initial disk mass. Replenishment via BHL or cloud capture would keep the disk at a radius $R_{\text{disk}} \sim R_{\text{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_{\text{disk}} > R_{\text{BHL}}$. Face-on accretion would then temporarily dominate the overall rate, shrinking the disk to maintain a balance between shrinkage and 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 disk evolution.

### ORCID iDs

Andrew J. Winter https://orcid.org/0000-0002-7501-9801

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

### References

Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 216
ALMA Partnership, Vlahakis, C., Hunter, T. R., et al. 2015, ApJL, 808, L4
Almendros-Abad, V., Manara, C. F., Testi, L., et al. 2024, A&A, 685, A118
Aly, H., Nealon, R., & Gonzalez, J.-F. 2024, MNRAS, 527, 4777
Andrews, S. M. 2020, ARA&A, 58, 483

Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129

Ansdell, M., Haworth, T. J., Williams, J. P., et al. 2020, AJ, 160, 248

Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

Begelman, M. C., & Shlosman, I. 2009, ApJL, 702, L5

Benisty, M., Dominik, C., Follette, K., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 605

Boccaletti, A., Di Folco, E., Pantin, E., et al. 2020, A&A, 637, L5

Bohn, A. J., Benisty, M., Perraut, K., et al. 2022, A&A, 658, A183

Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440

Bondi, H. 1952, MNRAS, 112, 195

Bonnell, I., & Bastien, P. 1992, ApJL, 401, L31

Borchert, E. M. A., Price, D. J., Pinte, C., & Cuello, N. 2022, MNRAS, 510, L37

Bosman, A. D., Appelgren, J., Bergin, E. A., Lambrechts, M., & Johansen, A. 2023, ApJL, 944, L53

Burgers, J. 1974, The Nonlinear Diffusion Equation: Asymptotic Solutions and Statistical Problems (Dordrecht: Reidel) https://books.google.de/books?id=zK6lswEACAAJ

Burleigh, K. J., McKee, C. F., Cunningham, A. J., Lee, A. T., & Klein, R. I. 2017, MNRAS, 468, 717

Carr, J. S., Tokunaga, A. T., & Najita, J. 2004, ApJ, 603, 213

Castelli, F., & Kurucz, R. L. 2003, in IAU Symp. 210, Modelling of Stellar Atmospheres, Poster Contributions, ed. N. Piskunov, W. W. Weiss, & D. F. Gray (San Francisco, CA: ASP), A20

Cazzoletti, P., Manara, C. F., Liu, H. B., et al. 2019, A&A, 626, A11

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698

Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485

Clarke, C. J., & Syer, D. 1996, MNRAS, 278, L23

Codella, C., Podio, L., Simone, M. D., et al. 2024, MNRAS, 528, 7383

Colman, T., Robitaille, J.-F., Hennebelle, P., et al. 2022, MNRAS, 514, 3670

Cuello, N., Ménard, F., & Price, D. J. 2023, EPJP, 138, 11

Deng, H., Ogilvie, G. I., & Mayer, L. 2020, MNRAS, 500, 4248

Derkink, A., Ginski, C., Pinilla, P., et al. 2024, A&A, 688, A149

Dong, R., Liu, H. B., Cuello, N., et al. 2022, NatAs, 6, 331

Dotter, A. 2016, ApJS, 222, 8

Duchêne, G., Ménard, F., Stapelfeldt, K., & Duvert, G. 2003, A&A, 400, 559

Dullemond, C. P., Kimmig, C. N., & Zanazzi, J. J. 2022, MNRAS, 511, 2925

Dullemond, C. P., Küffmeier, M., Goicovic, F., et al. 2019, A&A, 628, A20

Edgar, R. 2004, NewAR, 48, 843

Edgar, R., & Clarke, C. 2003, MNRAS, 338, 962

Edgar, R., & Clarke, C. 2004, MNRAS, 349, 678

El Mellah, I., & Casse, F. 2015, MNRAS, 454, 2657

Ercolano, B., & Pascucci, I. 2017, RSOS, 4, 170114

Evans, N. J., II 1999, ARA&A, 37, 311

Fabian, A. C., Pringle, J. E., & Rees, M. J. 1975, MNRAS, 172, 15p

Facchini, S., van Dishoeck, E. F., Manara, C. F., et al. 2019, A&A, 626, L2

Fairbairn, C. W., & Ogilvie, G. I. 2023, MNRAS, 520, 1022

Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A., et al. 2015, ApJ, 813, 99

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

Francis, L., & van der Marel, N. 2020, ApJ, 892, 111

Fukagawa, M., Hayashi, M., Tamura, M., et al. 2004, ApJL, 605, L53

Galli, P. A. B., Bouy, H., Olivares, J., et al. 2020, A&A, 643, A148

Galli, P. A. B., Bouy, H., Olivares, J., et al. 2021a, A&A, 646, A46

Galli, P. A. B., Bouy, H., Olivares, J., et al. 2021b, A&A, 654, A122

Garufi, A., Avenhaus, H., Pérez, S., et al. 2020, A&A, 633, A82

Garufi, A., Ginski, C., van Holstein, R. G., et al. 2024, A&A, 685, A53

Garufi, A., Podio, L., Codella, C., et al. 2022, A&A, 658, A104

Ginski, C., Facchini, S., Huang, J., et al. 2021, ApJL, 908, L25

Girichidis, P., Konstandin, L., Whitworth, A. P., & Klessen, R. S. 2014, ApJ, 781, 91

Goldsmith, P. F., Heyer, M., Narayanan, G., et al. 2008, ApJ, 680, 428

Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 539

Grady, C. A., Woodgate, B., Bruhweiler, F. C., et al. 1999, ApJL, 523, L151

Grant, S. L., Stapper, L. M., Hogerheijde, M. R., et al. 2023, AJ, 166, 147

Green, A., Wong, T., Indebetouw, M., et al. 2024, ApJ, 966, 51

Gupta, A., Miotello, A., Manara, C. F., et al. 2023, A&A, 670, L8

Gupta, A., Miotello, A., Williams, J. P., et al. 2024, A&A, 683, A133

Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJL, 553, L153

Hales, A. S., Gupta, A., Ruíz-Rodríguez, D., et al. 2024, ApJ, 966, 96

Hales, A. S., Pérez, S., Gonzalez-Ruilova, C., et al. 2020, ApJ, 900, 7

Hamer, J. H., & Schlaufman, K. C. 2019, AJ, 158, 190

Hanawa, T., Garufi, A., Podio, L., Codella, C., & Segura-Cox, D. 2024, MNRAS, 528, 6581

Hartmann, L., D'Alessio, P., Calvet, N., & Muzerolle, J. 2006, ApJ, 648, 484

Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207

Haworth, T. J., Clarke, C. J., Rahman, W., Winter, A. J., & Facchini, S. 2018, MNRAS, 481, 452

Haworth, T. J., Coleman, G. A. L., Qiao, L., Sellek, A. D., & Askari, K. 2023, MNRAS, 526, 4315

Hayashi, M., Ohashi, N., & Miyama, S. M. 1993, ApJL, 418, L71

Herbig, G. H. 1966, VA, 8, 109

Herbig, G. H. 1977, ApJ, 217, 693

Heyer, M. H., & Brunt, C. M. 2004, ApJL, 615, L45

Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540

Hopkins, P. F. 2012, MNRAS, 423, 2016

Hoyle, F. 1953, ApJ, 118, 513

Huang, J., Bergin, E. A., Bae, J., Benisty, M., & Andrews, S. M. 2023, ApJ, 943, 107

Huang, J., Bergin, E. A., Öberg, K. I., et al. 2021, ApJS, 257, 19

Hunt, R. 1971, MNRAS, 154, 141

Hunt, R. 1979, MNRAS, 188, 83

Ilee, J. D., Fairland, J., Oudmaijer, R. D., et al. 2014, MNRAS, 445, 3723

Joncour, I., Duchêne, G., Moraux, E., & Motte, F. 2018, A&A, 620, A27

Jørgensen, J. K., van Dishoeck, E. F., Visser, R., et al. 2009, A&A, 507, 861

Kimmig, C. N., & Dullemond, C. P. 2024, arXiv:2406.02754

Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17

Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2012, ApJ, 745, 19

Krieger, A., Kuffmeier, M., Reissl, S., et al. 2024, A&A, 686, A111

Krolikowski, D. M., Kraus, A. L., & Rizzuto, A. C. 2021, AJ, 162, 110

Krumholz, M. R., McKee, C. F., & Klein, R. I. 2006, ApJ, 638, 369

Kuffmeier, M., Dullemond, C. P., Reissl, S., & Goicovic, F. G. 2021, A&A, 656, A161

Kuffmeier, M., Goicovic, F. G., & Dullemond, C. P. 2020, A&A, 633, A3

Kuffmeier, M., Jensen, S. S., & Haugbølle, T. 2023, EPJP, 138, 272

Kuffmeier, M., Pineda, J. E., Segura-Cox, D., & Haugbølle, T. 2024, arXiv:2405.12670

Larson, R. B. 1981, MNRAS, 194, 809

Lee, A. T., Cunningham, A. J., McKee, C. F., & Klein, R. I. 2014, ApJ, 783, 50

Lee, H. M., & Ostriker, J. P. 1986, ApJ, 310, 176

Lee, J., Song, I., & Murphy, S. 2020, MNRAS, 494, 62

Lesur, G., Flock, M., Ercolano, B., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 465

Lin, D. N. C., Papaloizou, J., & Faulkner, J. 1985, MNRAS, 212, 105

Lodato, G., & Clarke, C. J. 2004, MNRAS, 353, 841

Longmore, S. N., Chevance, M., & Kruijssen, J. M. D. 2021, ApJL, 911, L16

Luhman, K. L. 2007, ApJS, 173, 104

Luhman, K. L., & Esplin, T. L. 2020, AJ, 160, 44

Luhman, K. L., & Rieke, G. H. 1999, ApJ, 525, 440

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 539

Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3

Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18

Marino, S., Perez, S., & Casassus, S. 2015, ApJL, 798, L44

Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364

Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76

Mawet, D., Absil, O., Montagnier, G., et al. 2012, A&A, 544, A131

Mellema, G., Arthur, S. J., Henney, W. J., Iliev, I. T., & Shapiro, P. R. 2006, ApJ, 647, 397

Mesa, D., Ginski, C., Gratton, R., et al. 2022, A&A, 658, A63

Michel, A., van der Marel, N., & Matthews, B. C. 2021, ApJ, 921, 72

Miret-Roig, N., Galli, P. A. B., Olivares, J., et al. 2022, A&A, 667, A163

Miyazaki, S., & Masuda, K. 2023, AJ, 166, 209

Moeckel, N., & Throop, H. B. 2009, ApJ, 707, 268

Mok, A., Chandar, R., & Fall, S. M. 2020, ApJ, 893, 135

Mordasini, C., & Burn, R. 2024, RvMG, 90, 55

Najita, J. R., Doppmann, G. W., Carr, J. S., Graham, J. R., & Eisner, J. A. 2009, ApJ, 691, 738

Nakajima, T., & Golimowski, D. A. 1995, AJ, 109, 1181

Nusser, A., & Silk, J. 2022, MNRAS, 509, 2979
Ostriker, E. C., Gammie, C. F., & Stone, J. M. 1999, ApJ, 513, 259

Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415

Padoan, P., Kritsuk, A., Norman, M. L., & Nordlund, Å. 2005, ApJL, 622, L61

Padoan, P., Pan, L., Pelkonen, V.-M., Haugboelle, T., & Nordlund, A. 2024, arXiv:2405.07334

Pan, L., & Scannapieco, E. 2010, ApJ, 721, 1765

Paneque-Carreño, T., Izquierdo, A. F., Teague, R., et al. 2024, A&A, 684, A174

Pascucci, I., Cabrit, S., Edwards, S., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 567

Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3

Pecaut, M. J., & Mamajek, E. E. 2016, MNRAS, 461, 794

Pelkonen, V.-M., Padoan, P., Juvela, M., Haugbølle, T., & Nordlund, Å. 2024, arXiv:2405.06520

Pérez, S., Hales, A., Liu, H. B., et al. 2020, ApJ, 889, 59

Pfalzner, S., Dehghani, S., & Michel, A. 2022, ApJL, 939, L10

Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691

Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S. Inutsuka et al. (San Francisco, CA: ASP), 233

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Pinte, C., Fouchet, L., Ménard, F., Gonzalez, J. F., & Duchêne, G. 2007, A&A, 469, 963

Pizzati, E., Rosotti, G. P., & Tabone, B. 2023, MNRAS, 524, 3184

Preibisch, T., Brown, A. G. A., Bridges, T., Guenther, E., & Zinnecker, H. 2002, AJ, 124, 404

Raycheva, N. C., Haverkorn, M., Ideguchi, S., et al. 2022, A&A, 663, A170

Reid, M. J., Menten, K. M., Zheng, X. W., et al. 2009, ApJ, 700, 137

Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52

Robe, H. 1968, AnAp, 31, 475

Rosotti, G. P. 2023, NewAR, 96, 101674

Rosotti, G. P., Tazzari, M., Booth, R. A., et al. 2019, MNRAS, 486, 4829

Ruffert, M. 1994, A&AS, 106, 505

Ruffert, M. 1996, A&A, 311, 817

Sanchis, E., Testi, L., Natta, A., et al. 2021, A&A, 649, A19

Scicluna, P., Rosotti, G., Dale, J. E., & Testi, L. 2014, A&A, 566, L3

Scott, D. W. 1992, Multivariate Density Estimation (New York: John Wiley & Sons)

Shima, E., Matsuda, T., Takeda, H., & Sawada, K. 1985, MNRAS, 217, 367

Sicilia-Aguilar, A., Hartmann, L. W., Hernández, J., Briceño, C., & Calvet, N. 2005, AJ, 130, 188

Silverberg, S. M., Wisniewski, J. P., Kuchner, M. J., et al. 2020, ApJ, 890, 106

Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson: Univ. of Arizona Press), 219

Speedie, J., Dong, R., Hall, C., et al. 2023, ALMA at 10 years: Past, Present, and Future, 80

Tabone, B., Rosotti, G. P., Cridland, A. J., Armitage, P. J., & Lodato, G. 2022, MNRAS, 512, 2290

Teague, R., Henning, T., Guilloteau, S., et al. 2018, ApJ, 864, 133

Testi, L., Natta, A., Manara, C. F., et al. 2022, A&A, 663, A98

Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111

Throop, H. B., & Bally, J. 2008, AJ, 135, 2380

Tobin, J. J., Sheehan, P. D., Reynolds, N., et al. 2020, ApJ, 905, 162

Trapman, L., Facchini, S., Hogerheijde, M. R., van Dishoeck, E. F., & Bruderer, S. 2019, A&A, 629, A79

Trapman, L., Rosotti, G., Zhang, K., & Tabone, B. 2023, ApJ, 954, 41

Ueda, T., Kataoka, A., Zhang, S., et al. 2021, ApJ, 913, 117

van Dishoeck, E. F., & Black, J. H. 1988, ApJ, 334, 771

van Terwisga, S. E., & Hacar, A. 2023, A&A, 673, L2

Vazquez-Semadeni, E. 1994, ApJ, 423, 681

Vázquez-Semadeni, E., Gómez, G. C., & González-Samaniego, A. 2024, MNRAS, 530, 3445

Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés, M. 2019, MNRAS, 490, 3061

Vieira, K., Korchagin, V., Carraro, G., & Lutsenko, A. 2023, Galax, 11, 77

Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164

Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11

Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2024, ApJ, 961, 95

Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323

Vorobyov, E. I., & Basu, S. 2005, ApJL, 633, L137

Walsh, C., Daley, C., Facchini, S., & Juhász, A. 2017, A&A, 607, A114

Weder, J., Mordasini, C., & Emsenhuber, A. 2023, A&A, 674, A165

Welch, W. J., Hartmann, L., Helfer, T., & Briceño, C. 2000, ApJ, 540, 362

Wijnen, T. P. G., Pols, O. R., Pelupessy, F. I., & Portegies Zwart, S. 2017, A&A, 602, A52

Williams, J. P., Cieza, L., Hales, A., et al. 2019, ApJL, 875, L9

Winter, A. J., & Alexander, R. 2021, MNRAS, 505, 869

Winter, A. J., & Haworth, T. J. 2022, EPJP, 137, 1132

Winter, A. J., Kruijssen, J. M. D., Longmore, S. N., & Chevance, M. 2020, Natur, 586, 528

Wolfire, M. G., & Cassinelli, J. P. 1986, ApJ, 310, 207

Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, ApJ, 880, 69

Zhang, K., Schwarz, K. R., & Bergin, E. A. 2020, ApJL, 891, L17

Zink, J. K., Hardegree-Ullman, K. K., Christiansen, J. L., et al. 2023, AJ, 165, 262

Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Natur, 601, 334

Zurlo, A., Weber, P., Pérez, S., et al. 2024, A&A, 686, A309