© The Authors 2025

1 Introduction

Theoretical studies of planet formation assume idealized disk models, in line with a traditional picture of star formation where the stellar mass reservoir is fully contained in a protostellar core: when the collapse of the protostellar core is finished, the circumstellar disk is fully formed and its total mass can only decrease over time. This traditional picture of star formation, known as core collapse, has been recently challenged. Simulations of star formation under physical conditions characteristic of molecular clouds (MCs), such as Larson’s relations (e.g., Larson 1981; Solomon et al. 1987; Heyer & Brunt 2004), demonstrate that the timescale of star formation is longer than the free-fall time of prestellar cores, with a significant fraction of the stellar mass coming from larger scales through converging flows (Padoan et al. 2020; Pelkonen et al. 2021; Kuffmeier et al. 2023). Because these converging flows are randomly occurring in supersonic turbulence without the aid of stellar gravity, we refer to them as “inertial inflows” (Padoan et al. 2020). The formation time of a star (the duration of the inertial inflows) increases with the final stellar mass, and can exceed 1 Myr for massive stars. Being characterized by infall from random dense filaments (Kuffmeier et al. 2017, 2018), this formation phase is quite different from the spherical collapse of a protostellar core, impacting the formation and evolution of protoplanetary disks.

The same simulations have also shown that once the star has reached its final mass and emerges from the native dense gas, as a pre-main sequence (PMS) class II star, further accretion from larger scales is still possible at a much lower rate through Bondi- Hoyle (BH) accretion. Although too low to increase the mass of the star significantly, early BH accretion may be important for the mass budget of the disk (Padoan et al. 2005; Throop & Bally 2008; Padoan et al. 2014, 2025), so it should be included in models of disk evolution. As for the inertial inflows in the formation phase, this mass infall deviates strongly from spherical symmetry, because the relative velocity between the captured gas and the star is much larger than the sound speed, so the accretion proceeds through dense filaments in the BH trail of the star (Padoan et al. 2025). The idea of a very dynamic evolution of disks, characterized by destabilizing filamentary infall well after the initial protostellar collapse, has also been supported by simulations following the early evolution of disks in realistic large-scale environments (e.g., Kuffmeier et al. 2017, 2018), in idealized setups of late disks (Moeckel & Throop 2009; Kuffmeier et al. 2020, 2021), and by using large-scale models that capture the infall at late times (Kuffmeier et al. 2023, 2024).

Padoan et al. (2025) have shown that BH infall in class II stars can account not only for the observed disk masses, but also for the disk specific angular momentum, or disk sizes, which implies that BH infall, rather than the initial protostellar collapse, may be the main formation mechanism of class II disks. They found that, in supersonic turbulence, the specific angular momentum increases linearly with scale, j(R) ∼ R, and they derived the angular momentum of gas captured by BH infall, j_BH, as the characteristic value of j(R) at a scale equal to the BH radius, R_BH, of the PMS star, j_BH = j(R_BH), using the following expression for the BH radius: $$R_{\text{BH}}=\frac{2G{M}_{\text{star}}}{c_{\text{s}}^2+v_{\text{rel}}^2},$$(1)

where c_s is the sound speed and v_rel the relative velocity between the gas and the star. The expression reduces to the Hoyle- Lyttleton radius, $R_{\text{HL}}=2G{M}_{\text{star}}/v_{\text{rel}}^2$, in the pressureless case where c_s = 0 (Hoyle & Lyttleton 1939), and to the Bondi radius for spherical accretion, $R_{\text{B}}=G{M}_{\text{star}}/c_{\text{s}}^2$, in the limit of v_rel = c_s (Bondi 1952). The characteristic specific angular momentum of the gas captured by the star at the BH radius is $$j_{\text{BH}}=4.3\times {10}^{20} {\text{cm}}^2 {\text{s}}^{-1}{\left({\sigma }_{\text{v},\text{rel}}/2 \text{km} {\text{s}}^{-1}\right)}^{-1}\left(M_{\text{star }}/1M_{\odot }\right)$$(2)

and the corresponding disk radius (assuming a simple Keplerian disk model), $$R_{\text{d}}=3.6\times {10}^2\text{au}{\left({\sigma }_{\text{v},\text{rel}}/2 \text{km} {\text{s}}^{-1}\right)}^{-2}\left(M_{\text{star}}/1M_{\odot }\right)=R_{\text{BH}}/4.1$$(3)

where σ_v,rel, the standard deviation of the gas velocity relative to the star’s velocity, has been defined to include the sound speed, ${\sigma }_{\text{v},\text{rel}}^2=c_{\text{s}}^2+v_{\text{rel}}^2$, and we have adopted the same normalization of 2 km s⁻¹ as in Padoan et al. (2025), because it is a typical value for the stars in the simulation in the early BH phase.

Eqs. (2) and (3) should be considered as upper limits for the average values of j_BH and R_d, because the angular momentum of the gas captured by a star in its trajectory is generally not constant over a ∼1 Myr timescale, causing partial cancelation in the average angular momentum vector, as mass accumulates in the disk. In addition, the scatter in the turbulence scaling of j(R) and σ_v(R) is expected to cause a significant scatter of M_d, j_BH, and R_d around the mean values predicted by the above relations. The time evolution of the angular momentum of the accreting gas was not addressed in Padoan et al. (2025), because of the Eulerian nature of the numerical analysis, where the angular momentum was computed within a sphere of radius equal to R_BH around each star. In this work, we continue the analysis of the same simulation, but with a Lagrangian perspective, using tracer particles to track the time evolution of the angular momentum of the gas accreted by the PMS stars in the simulation. We also address the potential for observational studies of BH infall in PMS stars, by carrying out synthetic James Webb Space Telescope (JWST) observations of BH trails from the simulation.

The paper is structured as follows. In Sect. 2 we briefly describe the simulation and the BH infall rates in comparison with observed accretion rates of PMS stars (Sect. 2.1); we also present the method of tracking the angular momentum evolution using the tracer particles in the simulation (Sect. 2.1). Our results on the angular momentum of the gas captured by BH infall and its time evolution are presented in Sect. 3. In Sect. 4 we describe the synthetic near-infrared (NIR) scattering observations around selected PMS stars from the simulation. Results are further discussed in Sect. 5 in the context of actual observations of streamers, and our conclusions are summarized in Sect. 6.

Fig. 1 Column density maps of the simulation from three cardinal directions (x, y, and z). The sink particles are shown with star symbols.

2 Numerical methods

2.1 Simulation and accretion rates

This work is based on the analysis of a numerical simulation of randomly driven, supersonic, magneto-hydrodynamic (MHD) turbulence, including self-gravity and sink particles to capture the star-formation process. It is the same global simulation used in Jørgensen et al. (2022), Jensen et al. (2023), Kuffmeier et al. (2023, 2024), and Padoan et al. (2025), and it is equivalent to the phhigh reference simulation in Haugbølle et al. (2018), except that the numerical resolution (the root grid) is larger by a factor of two, corresponding to a factor of 8 increase in the number of computational cells. Details of the numerical method can be found in Haugbølle et al. (2018) and are only briefly summarized here. The simulation solves the three-dimensional MHD equations with the adaptive-mesh-refinement (AMR) code RAMSES (Teyssier 2007), with a root grid of 512³ cells, representing a region with a size of 4 pc, and six levels of refinement, corresponding to a smallest cell size ∆x = 25 au. The total mass, mean density, and mean magnetic field strength are M_box = 3000 M_⊙, ${\overline{n}}_{\text{H}}=1897$ cm⁻³, and $\overline{B}=7.2$ µG. We use an isothermal equation of state with a gas temperature T = 10 K and periodic boundary conditions. The turbulence is first driven, without self-gravity, for ≈20 dynamical times, with a random solenoidal acceleration giving an rms sonic Mach number of approximately 10. The simulation is then continued for ∼2 Myr (snapshot number 448) with self-gravity and sink particles (maintaining the random driving force) to capture the ongoing formation of individual stars, yielding 317 stars with a mass distribution consistent with the observed stellar initial mass function (IMF) (Salpeter 1955; Chabrier 2005). Figure 1 shows the column density maps of the simulation from three lines-of sight, with the stars overplotted on the maps.

The spatial resolution of the simulation is insufficient to resolve the accretion disks around the sink particles, and the mass and angular momentum loss by winds and jets. Thus, the sink particle should be viewed as containing both the star and its accretion disk. To mimic the mass loss by protostellar winds and jets, the sub-grid model for mass accretion onto the sink particles (see Haugbølle et al. 2018, for details) prescribes that only 50% of the gravitationally captured mass (the infall onto the disk) is accreted onto the sink particles; the other 50% of the captured mass is removed from the simulation. Although we use the terms “accretion” and “infall” interchangeably, with regard to the simulation, we always refer to the mass infalling toward the (unresolved) disks, but whenever we give a numerical value for the infall (e.g., in Fig. 2), it is already halved to account for the mass loss. We assume that the infalling mass will settle onto a disk according to its angular momentum, but the accretion from the disk to the star (i.e., the depletion of the disk mass), which is measured by the observers, is not captured by the simulation. However, by comparing the computed infall rates (from the interstellar medium (ISM) to the disks) with the observed accretion rates (from the disks to the stars) of real stars we can evaluate the importance of the infall for the mass budget of disks.

Fig. 2 shows the infall rates, Ṁ, versus the stellar mass, M, for a subset of the stars found in the final six snapshots of the simulation (a period of 30 kyr). Because this work and the observational samples we compare with focus on class II stars, we have selected for this figure the sink particles that are representative of class II stars according to the criterion discussed in Padoan et al. (2025), that is the gas density at the sink particle location must be ≤5 × 10⁴ cm⁻³. This criterion yields a total of 672 stars (from the six snapshots), with an upper limit on the accretion rate of order 10⁻⁷ M_⊙ yr⁻¹. The median and maximum ages of the stars are 1.1 and 2.0 Myr respectively, consistent with estimated ages of class II stars (Evans et al. 2009). Fig. 2 also shows the accretion rates of 207 class II stars in four star-forming regions with estimated ages ≲3 Myr (Testi et al. 2022; Gangi et al. 2022). The infall rates from the simulation are comparable to the observed accretion rates, except for the lack of infall rates above ∼10⁻⁷ M_⊙ yr⁻¹, due to our conservative classification of class II sink particles, and the lack of low accretion rates in the observations, because of the sensitivity limit arising from chromospheric noise in the spectra (Manara et al. 2013, 2017). As the computed infall rates onto the disks are comparable to the observed accretion rates from the disks to the stars, this late infall must play an important role in the disk mass budget.

Fig. 2 Mass accretion rates as a function of stellar mass. Red diamonds are Class II YSOs in Chamaeleon I, Lupus, and L1688 from Testi et al. (2022), and in Taurus from Gangi et al. (2022). Blue points are large-scale infall rates of star particles from the MHD simulation (see Section 2.1) with local gas density lower than 5 × 104 cm−3. The dashed and dotted black lines are the expected chromospheric noise levels (Manara et al. 2013, 2017) using the 3 Myr isochrones of Baraffe et al. (2015) and Feiden (2016), respectively.

2.2 Tracer particles and angular momentum evolution

The simulation includes 512³ passively advected tracer particles, each representing a fluid element with a characteristic mass of approximately 2.2 × 10⁻⁵ M_⊙. The tracer particles are advected with a symplectic kick-drift-kick scheme, using the fluid velocities sampled with cloud-in-cell interpolation. They record all the hydrodynamic variables and are tagged once they accrete onto a sink particle. By selecting all the tracer particles that eventually accrete onto a given star, we can track the Lagrangian history of fluid elements contributing to the mass of that star, with a time resolution of 5 kyr, that is the time interval between the simulation snapshots that were saved. Because observed disk masses in class II stars are of order 1% of the stellar mass on average (Testi et al. 2022), we study the time evolution of the accreting gas in the last ten epochs, where each epoch corresponds to the accretion of 1% of the stellar mass, as representative of the accretion process that is most likely responsible for the current dynamical state of class II disks.

Although the disk is not resolved in the simulation, we aim at measuring the angular momentum supplied to the disk by the accreting gas. To that purpose, following the reasoning in Padoan et al. (2025), which is summarized in Sect. 1, we assume that the disk receives the angular momentum of the gas gravitationally captured by the star, hence the angular momentum is measured at a scale equal to the Bondi-Hoyle radius of the star, R_BH. This is done by computing the angular momentum vector of each accreting tracer particle, J_tr, at a distance R_BH from the star: $$J_{\text{tr}}=m_{\text{tr}}{r}_{\text{tr}}\times {\upsilon }_{\text{tr}},$$(4)

where m_tr is the mass of the individual tracer particle, and r_tr and v_tr are the position and velocity vectors of the individual tracer particle with respect to the accreting star when the tracer particle crosses R_BH. Due to the tracer particle data being recorded every 5 kyr, the last position of a tracer is often outside (low- mass stars) or inside (high-mass stars) of R_BH. Thus, we need to interpolate or extrapolate J_tr at the distance R_BH based on its values at the available distances. For that purpose, we examine the dependence of the specific angular momentum of each tracer particle, j_tr, on the distance, r_tr, between the tracer and the star accreting it. We expect that j_tr scales with r_tr according to different laws outside and inside R_BH, hence we consider the two cases separately, as explained in the following paragraph.

We first take all the tracers contributing to the last 10% of the mass of each star, splitting them into ten 1% mass bins (epochs) in chronological order of their accretion time. The chronological order among the tracers accreted in the same snapshot is not available, so such tracers are sorted by their index number, which, by this point in the simulation, is effectively a random order. For every tracer particle in each epoch and each star, we take the tracer’s last ten positions (50 kyr) before crossing R_BH (when the tracer is outside R_BH), and derive the least-square fit slope of j_tr versus r in log-log space, $j_{\text{tr}}\propto r_{\text{tr}}^{a_{\text{out},\text{ in }}}$. We then compute the mean value of the slope, a_out, by averaging the slopes of all the tracer particles in each mass bin of each star. The procedure is then repeated for the 10 positions after the tracer particle has crossed R_BH (when the tracer is inside R_BH), yielding the mean slopes a_in. It is possible that the tracer particles spend less than 50 kyr between crossing R_BH and accreting onto the star, in which case the least-square fit to derive the slope relies on less than 10 positions. When less than three positions are available, the slope is not derived for that tracer particle. If no slope is derived for any of the tracers within a single epoch (1% mass bin) of a specific star, no a_in value is derived for that epoch. This happens only for low-mass stars. Moreover, in the case of low-mass stars, the value of R_BH is usually low, corresponding to scales that are not well resolved in the simulation. Hence, we prefer not to rely on the values of a_in derived for low-mass stars, but only on those for stars with mass ≥0.5 M_⊙.

Figure 3 shows the values of a_out (upper panel) and a_in (lower panel) as a function of the star’s final mass. The horizontal red lines correspond to the median values of the slopes, measured over the whole mass range for a_out, and for M ≥ 0.5 M_⊙ in the case of a_in, as explained above. These median values are a_out = 0.82 and a_in = 0.53. Interestingly, the median value of a_in shows that the tracer particles inside R_BH have velocities close to Keplerian, which is consistent with the prediction that the disk size should be just a few times smaller than R_BH (see Eq. (3)). We apply these median scaling exponents to rescale the angular momentum of each tracer particle to a distance equal to R_BH. For each tracer with at least one recorded position inside R_BH, we scale its first inner J_tr value back to R_BH using the a_in slope. If no position inside R_BH is recorded, we take the tracer’s last position before accretion and scale its J_tr value to R_BH using the median slope a_out. Finally, we scale the angular momentum of each tracer to the distance of the expected disk radius from Eq. (3), R_d = R_BH /4.1, using the median slope a_in. All rescaling applies to the modulus of J_tr, the direction is left unchanged. The scaling procedure is summarized in the following formulas, $$J_{\text{tr},\text{disk}}=\{\begin{array}{ll}{J}_{\text{tr}}\times (R_{\text{BH}}/r_{\text{tr}}{)}^{a_{\text{in}}}/{4.1}^{a_{\text{in}}},\hfill & \text{ if }{r}_{\text{tr}}<R_{\text{BH}},\hfill \\ J_{\text{tr}}\times {\left(R_{\text{BH}}/r_{\text{tr}}\right)}^{a_{\text{out}}}/{4.1}^{a_{\text{in}}},\hfill & \text{ if }{r}_{\text{tr}}\ge R_{\text{BH}},\hfill \end{array}$$(5)

where J_tr,disk is the rescaled angular momentum vector of an individual tracer particle at the distance of the expected outer edge of the disk.

The total angular momentum vector, J, in each epoch of each star, is then calculated by summing over the angular momentum vectors of the tracers accreted during that epoch: $$J={\displaystyle \sum _{i=0}^n{J}_{\text{tr},\text{disk},\text{i}}}.$$(6)

The corresponding specific angular momentum is derived from the length of the angular momentum vector divided by the mass accreted during that epoch (1% of the final stellar mass): $$j=\frac{\left|\mathit{J}\right|}{{{\displaystyle \sum }}_{i=0}^n{m}_{tr,\text{i}}},$$(7)

where m_tr,i is the mass of the individual tracer particles that are accreted during that particular epoch.

In addition, we calculate the mass-weighted average angle, δϕ, between the angular momentum vector of individual mass tracer particles, J_tr,disk,i, and their summed up angular momentum vector of each epoch, J. This is a measure of how ordered the accretion was within an epoch (“tracers-to-disk” dispersion). We also calculate the angle ∆ϕ between the angular momentum vectors of consecutive epochs (J_k vs. J_k−1 ), measuring the misalignment between the epochs (“disk-to-disk” misalignment). We use the following equations: $$\delta \varphi =\frac{{\displaystyle {\sum }_{i=0}^n{m}_{\text{tr},\text{i}}}\arccos \left(\frac{J_{\text{t},\text{disk},\text{i}}\cdot J}{\left|J_{\text{tr},\text{disk},\text{i}}\right||J\mid }\right)}{{\displaystyle {\sum }_{i=0}^n{m}_{\text{tr},\text{i}}}},$$(8)

and $$\text{Δ}\varphi =\arccos \left(\frac{J_{\text{k}}\cdot J_{\text{k}-1}}{\left|J_{\text{k}}\right|\left|J_{\text{k}-1}\right|}\right).$$(9)

Figure 4 shows four example stars and the accretion of the last 10% of their mass, as well as local column density images of the middle point in time of the accretion of the last epoch (the last 1% of the accreted mass), which clearly shows the BH accretion wakes. They all share some general characteristics: the accretion rate starts high and then drops quite rapidly toward the end, as the star is no longer accreting from converging flows and instead accretes via BH accretion. The mass of the star does seem to influence in which epoch this switch happens. The least massive star (top-most) accretes very quickly at first, almost 1% of its mass per kyr, before its accretion rate quickly drops around 97% of its final mass. In the case of the second-least massive star, this transition is more gradual, and the last two epochs are accreted significantly more slowly. By contrast, in the two most massive stars, the accretion drops significantly only toward the last epoch. Since all stars that are fully formed in the simulation have both the accretion rate and the flow morphology consistent with BH accretion in the epoch when they receive the last 1% of their mass, we will use this last epoch (the final mass bin plotted in Figure 4) to study the contribution of BH accretion to the disk formation process.

In addition to the accretion rate, Figure 4 also shows the evolution of the magnitude of the specific angular momentum, j, and of the two angles δϕ and ∆ϕ defined above. There is a general trend of j to slightly decrease toward the last epoch, as might be expected when the stars decouple from their surrounding gas, causing a reduction of their Bondi-Hoyle radius (from the increased relative velocity between the star and the gas), with an exception to that trend shown in the third panel from the top. It is also notable that the angle δϕ is rather large, indicating the accretion of material with different angular-momentum directions, and that the angle ∆ϕ is occasionally very large, signaling episodes of sudden change in the average direction of the accreted angular momentum, which is discussed further in Sect. 3.2.

Fig. 3 Power-law slope of jtr versus rtr, $j_{\text{tr}}\propto r_{\text{tr}}^{a_{\text{out},\text{in}}}$, averaged for all the tracer particles within a single epoch, versus the final stellar mass (blue dots), as explained in the main text. The upper panel is for the scaling outside of RBH, the lower panel for the scaling inside RBH. The red lines show the median slope values aout = 0.82 (upper panel) and ain = 0.53 (lower panel, computed over masses ≥0.5 M⊙), and the shaded regions the rms value around the median.

Fig. 4 Four example stars (one per row) showing the accretion of the last 10% of their mass on the left (the top x-axis is the mass fraction, fM, and the bottom x-axis is the corresponding age of the star), and their late-accretion flows on the right. Left column: Angle ∆ϕ (dashed blue line, “disk-to-disk” misalignment), angle δϕ (dotted blue line, “tracers-to-disk” dispersion), specific angular momentum of the epoch (solid red line), j27 = j × 10−27, and the median accretion rate within 0.5% massbins (dash-dotted red line), Ṁ. Right column: Column density images of log($N_{{\text{H}}_2}$) in a 6400 au box, from three different cardinal directions (xyz), using the middle snapshot of the last epoch for each of the stars. The text in the images list the final mass, the age and the average accretion rate at this snapshot, and the minimum and the maximum limits used for the colorscale of the logarithm images, in all directions.

3 Disk angular momentum and alignment

3.1 Extraction of class II sources

Because we seek to characterize the contribution to the disk formation process by BH accretion in class II PMS stars, we need to target our analysis to stars that reach the class II phase before the end of the simulation. A simple way to identify stars that have reached nearly 100% of their real final mass before the end of the simulation is to consider their average accretion rate during the time they accumulate the last 1% of their total mass at the end of the simulation. If the accretion rate is low enough, the star is essentially fully formed, though the remaining BH accretion may still be significant for the mass balance of the disk. The value of the accretion rate that separates the formation phase from the BH phase (essentially class I or earlier from class II), is expected to be mass dependent, as we have previously demonstrated that the formation time (for example to reach 95% of the final mass), scales approximately as the square root of the final stellar mass, ∼M^1/2 (Padoan et al. 2014; Haugbølle et al. 2018; Padoan et al. 2020; Pelkonen et al. 2021), hence the average accretion rate scales approximately as ∼M^1/2 (though with a large scatter, increasing toward lower masses).

We express the accretion rate in terms of the time needed to accrete the last 1% of the stellar mass, t₉₉₉ ¹, which is plotted as a function of the stellar mass at the end of the simulation in the left panel of Figure 5. As apparent from the distribution of the points in the scatter plot, and confirmed by the histograms in the right panel of the same figure, the stars are separated into two populations, where the transition value of t₉₉₉ has a mass dependence, t₉₉₉ ∼ M^0.7, qualitatively consistent with the above argument based on the star-formation time (which would predict t₉₉₉ ∼ M^0.5). Because the transition between the two populations is well-defined in the scatter plot, we identify as class I (or earlier embedded phases) the stars with t₉₉₉ values below the transition (red dots in the left panel of Figure 5), and as class II all the stars above the transition (blue dots in the figure). This selection returns 160 class II stars, out of a total of 317 stars. The fraction of class II stars is lower than that of nearby star-forming regions because the median age of the stars is only 0.79 Myr (less than half the duration of the simulation with self-gravity, 2 Myr, as the star-formation rate increases with time).

Fig. 5 Scatter plot of the formation time of 90% of the last epoch (t999) versus mass on the left, and kernel density estimator (KDE) probability density plots of the formation time in four bins of final stellar mass on the right. We find the local minima of the KDE functions, which are then plotted on the scatter plot as black squares. The black dashed line is a fit to the black squares, log y = 0.72 ± 0.03 log x − 1.54 ± 0.03. This line is used to separate the sinks fed by slower Bondi-Hoyle accretion (blue points) and sinks accreting at a much higher rate, likely still in the class I or even earlier protostellar phases (red points).

3.2 Specific angular momentum and alignment of class II disks

For each class II star in the simulation, we compute the average specific angular momentum, j, of the last ten 1% mass bins (epochs) as explained in Sect. 2.2. We consider the last epoch as the most representative contribution of BH accretion to the formation of class II disks, and plot the absolute value of its specific angular momentum vector as a function of the final stellar mass, M, in Figure 6 (blue dots). There is a clear positive correlation between j and M, with a least-squares fit giving the relation: $$j=4\times {10}^{20} {\text{cm}}^2 {\text{s}}^{-1}{\left(M/1M_{\odot }\right)}^{1.1}.$$(10)

This result is in striking agreement with the analytical prediction (Eq. (2)), as well as the result of the Eulerian analysis in our previous paper (Padoan et al. 2025), which is shown by the black line in Figure 6. We also find a significant overlap between our predicted j values and those derived from observations of disk sizes (Hendler et al. 2020; Cieza et al. 2021; Long et al. 2022; Stapper et al. 2022), shown by the smaller red dots in Figure 6. Our predicted j values are on average a factor of two higher than those of the observed disks, which is to be expected because there are multiple processes that could contribute to angular momentum transport when the actual disk is assembled from the gravitationally captured gas. Such processes are not included here, as the simulation does not resolve the disk formation, and we predict the maximum angular momentum by sampling it at the disk edge.

We also performed the above study using disk mass proportions that scaled with the square root of the stellar mass: m_d/m_s = 0.01 × (m_s/1 M_⊙)^0.5, motivated by Testi et al. (2022). The results were very similar to the 1% case, but with lower statistics for the low-mass stars which did not have enough mass tracer particles.

An important consequence of BH accretion in class II disks is the possibility of disk misalignment. The mass accretion in the BH phase is a stochastic process that depends on the gas mass distribution encountered by the class II star in its trajectory within the parent cloud. Both the mass accretion rate and the accreted angular momentum can have significant time fluctuations. The angular momentum arises primarily from the relative velocity between the gas and the star, and the offset between the center of mass of the captured gas and the position of the star, rather than from local gas vorticity (see discussion in Padoan et al. 2025). Thus, the direction of the angular momentum vector of the accreted gas can change over time. To quantify this effect in the simulation, we compute the angle between the average J vectors in the last two epochs of every class II star, as defined in Eq. (9). The result is shown as a histogram in Figure 7. The ∆ϕ angle distribution peaks at low values, but cover a broad range from 0° to 160°, with a median value of 47°. This suggests that, besides a significant contribution to the disk mass budget, BH accretion in class II stars may also cause disk misalignment or warping, and strongly affect the angular momentum evolution of the disks. Thanks to the tracer particles, besides the change in direction of the mean J vector between different epochs, ∆ϕ, we can also compute the internal dispersion in the direction of the angular momentum of each tracer particle contributing to a single epoch (a 1% mass bin), δϕ, as defined in Eq. (8). The probability distribution of δϕ is shown by the shaded histogram in Figure 7, and has a median value of 70°. This corresponds to a cancelation factor of approximately three for the modulus of J obtained from the vector sum of individual tracer particles. So, the cancelation effect is significant, but the accreted angular momentum is far from being completely random to the point of resulting in nearly total cancelation. This relatively modest cancelation effect explains the good agreement of our result, expressed by Eq. (10), with the analytical prediction from Eq. (2) and the Eulerian analysis in Padoan et al. (2025).

Fig. 6 Scatter plot of specific angular momentum versus mass in the last epoch of the class II stars in the simulation (blue dots), and for observations of spatially resolved disks (Hendler et al. 2020; Cieza et al. 2021; Long et al. 2022; Stapper et al. 2022) (smaller red dots). The blue line is a least-squares fit (log y = (1.14 ± 0.07) log x + (20.7 ± 0.06)) to the blue data points, while the solid black line is the fit (log y = (0.86 ± 0.03) log x + (20.6 ± 0.03)) from the Eulerian analysis in Padoan et al. (2025) and the dashed black line a fit (log y = (0.71 ± 0.04) log x + (20.3 ± 0.02)) to the observations. The gray dotted lines are the specific angular momentum for a given disk radius (1, 100, and 10 000 au, with increasing j) as a function of stellar mass, using Eq. (29) in Padoan et al. (2025) for Keplerian disks.

Fig. 7 Histograms of the last ∆ϕ angle (blue step histogram), and of the last δϕ angle (blue-shaded bar histogram), for the class II stars in the simulation (see the main text for the definition of the angles). The median values of ∆ϕ and δϕ are shown by the black dashed and dotted vertical lines, respectively.

Fig. 8 Scatter plot of the specific angular momentum versus mass in the last epoch of the class I and earlier protostars. The red line is the leastsquares fit (log y = (0.90 ± 0.06) log x+ (20.8 ± 0.04)) to the data points, while the blue line is the fit to the simulated class II stars from Figure 6. The gray dotted lines are the specific angular momentum for a given disk radius (1, 100, and 10 000 au, with increasing j) as a function of stellar mass, using Eq. (29) in Padoan et al. (2025) for Keplerian disks.

3.3 Embedded phase versus class II

This work stresses the importance of late BH accretion in class II disks. However, the same analysis can be applied to protostars. For that purpose, we select the stars that never reach the class II phase in the simulation, based on the criterion presented in Sect. 3.1 and illustrated in Figure 5 (red dots in the figure). Similarly to Section 3.2, we compute j in the last 10 epochs, and plot the value of j of the last epoch versus the stellar mass in Figure 8 for the embedded stars. The plot is similar to the case of class II stars. The main differences are the presence of more stars with masses around 10 M_⊙ (most of the most massive stars are still forming at the end of the simulation), and a slightly shallower slope of the least-squares fit. The latter is primarily due to the presence of a significant number of low-mass stars with relatively large values of j compared to the class II stars. Most of the low-mass class II stars have spent significant time in the BH phase, during which time their v_rel (velocity difference to the surrounding gas) has increased, as the stars have gradually decoupled from the surrounding gas. Increased v_rel reduces R_BH (see Eq. (1)), hence decreasing the accreted j (see discussion of the time dependence in Padoan et al. 2025). The low-mass class I stars, instead, being still fed by inertial inflows (Padoan et al. 2020; Pelkonen et al. 2021), and even more so in earlier proto- stellar phases, are more tightly coupled with the surrounding gas (smaller v_rel ), hence their R_BH is comparatively larger.

A more striking difference between the accretion process during earlier protostellar phases (protostellar collapse and inertial inflows) and the class II phase (BH accretion) is illustrated in Figure 9. The unshaded histograms compare the last ∆ϕ distribution for the two classes, where ∆ϕ, as in Figure 7, is the angle between the J vectors in the last two epochs. The histogram peaks at a much lower angle for the protostars compared to the class II stars, with median values of 14° and 47° respectively. This difference is to be expected, because the star-formation process is characterized by the protostellar collapse first, and by inertial inflows of spatially coherent dense filaments later, while BH accretion, after the star has decoupled from those converging flows (or said flows have ended), is a more stochastic process. As a result, disk misalignment is expected to be more likely in the later class II phase than during the star-formation phase.

Figure 9 also shows the comparison between the δϕ distributions, where δϕ is a measure of the dispersion of the J vector orientation of individual tracer particles inside the last epoch. The difference between the two histograms is less pronounced, but in the same direction, as that between the ∆ϕ histograms: class I stars have slightly smaller δϕ values than class II stars, with median values of 59° and 70° respectively. The existence of a significant difference between the median values of δϕ and ∆ϕ also in the class I phase is expected, because the converging inertial flows forming the star come from different directions, each one bringing a different J orientation (as measured by individual tracer particles), even if they are rather coherent in time (hence the lower values of ∆ϕ).

Fig. 9 Histograms of the last ∆ϕ (red step histogram) and of the last δϕ (red-shaded bar histogram), for class I and earlier protostars. The median values of ∆ϕ and δϕ are shown by the black dashed and dotted vertical lines, respectively. For comparison, the corresponding histograms for class II stars from Fig. 7 have been replotted here in blue.

4 Synthetic observations

We are proposing that late BH accretion plays a major role in the mass and angular momentum budget of class II disks. This scenario leads to disk infall rates on the order of the observed accretion rates of class II stars, and disk specific angular momentum and sizes consistent with the observations (see also Padoan et al. 2025). In this section we show that our main hypothesis can be directly tested observationally. We demonstrate that observations with the James Webb Space Telescope (JWST) can detect the predicted BH trails and also help to constrain their corresponding infall rates.

4.1 Radiative transfer: SOC

We used the radiative transfer code SOC (Juvela 2019) to calculate images of NIR scattered light at 1.50 µm and 2.77 µm. The dust properties were based on the Weingartner & Draine (2001) model with the selective extinction of R_V = 5.5. The observed level of scattered light depends mainly on the albedo and the extinction cross section. The albedo, A, of the R_V = 5.5 model is higher than in the corresponding model for normal Milky Way dust with R_V = 3.1 (e.g. 0.68 vs. 0.54 at 1.5 µm). Such increase is expected inside molecular clouds as a result of grain evolution (ice deposition and grain coagulation), the effect increasing toward mid-infrared wavelengths (Steinacker et al. 2010). If the medium is optically thin at the wavelength of interest, a larger cross section (per gas mass) results in a higher intensity for the observed scattered light, but the effect would be the opposite for optically thick structures (Juvela et al. 2020). The NIR opacities of the R_V = 3.1 model would be some 30% higher than in the adopted R_V = 5.5 model, where τ_ext(1.5 µm) = 1 is reached at N(H) = 1.5 × 10²² cm⁻² and τ_ext (2.77 µm) = 1 at N(H) = 3.8 × 10²² cm⁻². Typical streamers have lower column densities, and observed class II YSOs in Chamaeleon I (Manara et al. 2017), Upper Scorpius (Manara et al. 2020), and Taurus (Gangi et al. 2022) have extinctions less than A_V = 10 mag, typically from 1 to 5 mag. Thus, although dust evolution may increase scattering via higher albedos, the negative effect of smaller absolute cross sections tends to be more significant.

To estimate the surface brightness due to the full line-of- sight, we started by calculating the scattered surface brightness for the full AMR model, using all AMR levels (smallest cell size of 25 au). All stars were treated as separate point sources, and the model was further illuminated from the outside by the normal interstellar radiation field (Mathis et al. 1983). Maps of scattered light were registered separately for photons originating from the diffuse background and from the point sources. These showed that the surface brightness is smooth at the small scales relevant for the present study, apart from the signal produced by local embedded stars.

For more detailed modeling of the environment around the stars (and to reach higher S/N for the resulting maps), we extracted from the full AMR grid smaller unigrid cubes that were each centered on a selected source (see below). The cubes had a size of 12 800 au and a resolution of 50 au. The scattered light due to all embedded sources inside the cube was simulated using 8 × 10⁸ photon packages. The numerical Monte Carlo noise in the resulting maps is mostly below 1% of the mean surface brightness of the map, and the relative noise is lower for the detected streamers.

The luminosity and the temperature of the stars were calculated with the Modules for Experiments in Stellar Astrophysics (MESA) code (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). The code has been augmented with modules for accretion of mass with a variable entropy and diagnostics to calculate the accretion luminosity, see Jensen & Haugbølle (2018) for details. In the following we will refer to the total stellar luminosity, the sum of photospheric and accretion luminosities, as the “luminosity” or “stellar luminosity”.

The Hertzsprung-Russell diagram for all stars in the last snapshot are shown in Fig. 10, where class I and class II sources are classified based on their local gas density: class I have local gas densities that are higher than 5 × 10⁴ cm⁻³, and class II have local gas densities that are lower (see Section 2.1). At the time of the snapshot, only stars above 3.5 solar masses have reached the zero-age main sequence and commenced burning hydrogen, even while most of them are still deeply embedded. In lightgreen the evolutionary tracks, including self-consistent account for mass and entropy accretion, highlighting the protostellar evolution. For 23% (15%) of the stars the accretion luminosity, L_acc, contribute more than 20% (50%) of the total luminosity, as can be seen by considering offsets between the magenta and green points for class I stars (the red and blue points for class II) in Fig. 10.

The stars of interest were selected from the last snapshot of the simulation, with a luminosity cut (L > 0.1 L_⊙) and an accretion-rate cut (10⁻¹¹ < Ṁ < 10⁻⁶ M_⊙/yr), irrespective of their classification as class I or class II stars. Based on an initial set of calculations, stars less luminous than this are not able to illuminate their surroundings sufficiently, and those with lower accretion rates do not have enough dust to scatter their light. All the other stars within the computational volume were included regardless of their luminosity and accretion rate, as they may still contribute to the surface brightness of the scattered light. However, this also means that we need to be careful about attributing a detection when it is due to a nearby star. To alleviate this concern, in the cases where we have more than one star in the subcube, we consider only the brightest one (called the dominant star). This results in 48 single stars and 21 dominant stars, observed from three directions with two wavelengths, for a total of 207 maps (69 stars times three directions) for each of the two wavelengths and for the column density.

Fig. 10 Hertzsprung-Russell diagram for all stars in the last snapshot (top panel) and for the 69 stars in the JWST synthetic observations (bottom panel). The small symbols (magenta and red) are for photospheric luminosities, while the large symbols (green and blue) are for stellar luminosities which include the accretion luminosity. Magenta and green points are class I stars, while red and blue points are class II. When Lacc is small (especially class II stars), the small points overlap with the larger ones, resulting in the appearance of red/magenta points with blue/green borders for class II and class I stars, respectively. The zeroage main sequence is indicated by the dashed black line, while the thin transparent light-green lines show the evolutionary tracks of all stars from an age of 104 yr to the time of the snapshot, and the thin transparent light-blue lines are the evolutionary tracks with episodic accretion included. Stellar evolution is computed following the methodology of Jensen & Haugbølle (2018).

4.2 Scattered NIR light with JWST

The 414 surface brightness maps are included with the corresponding 207 column density maps for the analysis. We assume that the distance to the target cloud is 200 pc, close to the distance of Chamaeleon I. We then perform a very simple ring subtraction to get rid of symmetrical features: we divide the map into 1-pixel thick rings around the star and subtract the median of the surface brightness within each of the rings from the pixels in the corresponding ring. This highlights the asymmetric structures that are expected from BH accretion in class II stars, or from accreting filaments in class I stars. As a very simple way to analyze the success of the detection, we calculate the signal-to-noise ratio (S/N) using the Exposure Time Calculator of James Webb Space Telescope (JWST). We assume a rapid readout, leading to an exposure time of 107.4 seconds (about two minutes), and an aperture size of 0.125” to match with our map pixel size, which results in 1σ noise levels of 0.043 MJy/sr and 0.0074 MJy/sr for filters F150W and F277W, respectively. We also test longer exposure times of about 8 minutes and 50 minutes. If we have traced a continuous structure from the center to at least 1000 au at S/N ≥5σ, then we consider it a successful detection of an asymmetric structure, hence a candidate for an accretion flow. In addition, we measure the surface brightness within a ring between 250 au and 1000 au, in order to compare it with the accretion rate and luminosity of the star.

As the MHD model did not have disks, we are unable to model disk shadows. We also did not consider the effect of the occulting masks for coronagraphy in order to block the central star, which will also cause some limits on the usable filters. However, as we are looking at larger distances from the central star, the results are still valid.

Figure 11 shows a relatively high-mass star (4.7 M_⊙), accreting at a rate of 3.2 × 10⁻¹⁰ M_⊙/yr. Despite its modest accretion rate, its high luminosity illuminates the surrounding dust, but also the streamer seen in the column density map, as the NIR photons scatter from the dust particles. The resulting surface brightness is enough that even after we perform a very rudimentary ring subtraction, the remaining surface brightness structure aligns almost perfectly with the BH trails seen in the column density map, leading us to conclude that this is a successful detection. The white contour shows that even a two-minutes exposure with JWST would be enough for this detection.

Figure 12 shows a lower-mass star (0.64 M_⊙), accreting at a rate of 3.1 × 10⁻⁹ M_⊙/yr. Its luminosity is only slightly higher than that of the Sun, but it is enough for the BH accretion feature to be detected. The opening angle of the BH wake is larger than in the previous case, indicating a smaller relative velocity between the star and the gas.

Figures 13 and 14 show the total number of maps in bins of the accretion rate and the luminosity of the sources, respectively, and the fraction of the successful detections of accretion features for both filters, using 2 min and 50 min exposure times. Obviously, the detection fraction is higher for longer exposure time and higher luminosity. It is also higher for the shorter wavelength, due to the increased stellar flux, and higher accretion rate, because that usually corresponds to higher column density in the accreting gas. The drop in the detection fraction at 100 L_⊙ is due to two sources with very low accretion rates, as shown in Figure 15. Even with just 2min exposure, JWST’s detection fraction rises above 50% for Ṁ ≥ 10⁻⁹ M_⊙ /yr, or for L ≥ 1 L_⊙. The longer exposures can push the 50% success rate to smaller values of Ṁ and L. Reducing the resolution of the map by averaging over several pixels would improve the signal-to-noise as well, at the cost of the resolution, as long as the resolution is still high enough to resolve the asymmetric structures.

Figure 15 shows scatter plots of the accretion rate and luminosity for all sources. The colors and the symbols refer to the detection rate per source with a 2 min exposure time and using either F150W (left panel) or F277W filter (right panel): a small black dot corresponds to a non-detection in all three maps (x, y, and z directions), a magenta triangle to a detection in one map out of three, a green square to a detection in two maps out of three, and a large blue dot to a detection in all three maps. The protostars (class I and earlier) are highlighted with red circles around the markers. The figure shows that a selection based on both accretion rate and luminosity thresholds improves the detection fraction. If we require both Ṁ > 5 × 10⁻¹⁰ M_⊙ /yr and L > 0.5 L_⊙ (red dotted lines), the detection fraction is 100% for F150W, regardless of further cuts by protostellar classification or maximum luminosity. For F277W filter with the same requirements, the detection fraction is approximately 90%, 74 detections out of 81 maps. If we consider only the class II stars, the detection fraction is still 84%, 38 detections out of 45 maps. If we also discard the three class II stars with L > 100 L_⊙, to reflect the scarcity of such sources in the solar neighborhood, we get a detection fraction of 81%, 29 detections out of 36 maps.

In Figure 16, we compare the NIR flux density in an aperture from 250 to 1000 au with both filters to the accretion rate. We do not find a significant difference between the distribution of single and dominant stars, so we do not differentiate between them in the plot. Instead, we distinguish between class I and class II, and also split the class II stars into high and low accretion rates, based on the criterion of Ṁ > 5 × 10⁻¹⁰ M_⊙/yr found in the previous paragraph. We find that there is a large scatter in flux density between stars of a similar accretion rate, which is partially explained by the luminosity differences between those stars. Thus, we divide the flux density with the square root of the luminosity, which helps to reduce the scatter and separate the different accretion regimes, especially when using F277W. Figure 16 shows a clear relationship between the accretion rate and the measured flux density, which is to be expected, as the higher accretion rate implies a denser medium that the star is traveling through and a denser accretion flow, meaning that there is more dust to scatter the NIR photons. We fit a power-law to each panel, which could be used as a rough order-of-magnitude estimate of the accretion rate using the two observables, the stellar luminosity and the F277W flux density, even if the surface brightness of the BH wake itself were too weak for a high signal- to-noise detection. However, more complicated radiative transfer modeling would be needed to accurately derive the accretion rate from the surface brightness, possibly using the morphology of the surface brightness structure rather than simple aperture photometry.

Fig. 11 JWST synthetic NIR observations (see Section 4) around a relatively high-mass star (located at the center of the map) in a 6400 × 6400 au frame. The title lists the sink particle number and direction of projection, the age, the mass, the luminosity and the accretion rate. The top row shows, in log scale the surface brightness in the two selected wavelengths (left and central images respectively, with the JWST filter as a label), and the local column density integrated across the 12 800 au depth of the selected volume (right image, with minimum and maximum log value of the color scale as a label). The contours (white, red, and yellow) are 5σ limits for 2, 8 and 50 min exposure time, respectively (the column density image shows the F277W contours). The bottom row shows the corresponding surface brightness images after the ring subtraction described in the text (the column density map is the same as in the upper row). Negative and zero pixel values are masked out (white background). Contours show the 5-sigma limits of the remaining surface brightness, after the ring subtraction.

Fig. 12 As Fig. 11, but for a lower-mass star.

Fig. 13 Total number of maps in bins of the accretion rate, and the detection fractions of the maps where accretion features were detected with S/N ≥5σ at least 1000 au from the central star. The black line is the number of all the maps, non-detections included, with y-axis on the left. The colored lines (blue and red for F150W and F277W filters, and dashed and dotted for 2 min and 50 min exposure times, respectively) and the y-axis on the right show the successful detection fraction of the maps.

Fig. 14 As Fig. 13, but in bins of the stellar luminosity instead of the accretion rate.

Fig. 15 Scatter plots of stellar luminosity versus accretion rate for all the synthetic observations, with symbols and colors referring to the detection rate after 2 min of exposure time with the F150W (left panel) and F277W (right panel) filters. The small black dots correspond to a non-detection of an accretion feature in all three directions (x, y, and z), the magenta triangles detection in only one direction, the green squares detection in two directions, and the large blue dots detection in all three directions. The protostars (class I or earlier) are highlighted by red circles around the markers. The red dotted vertical and horizontal lines mark the thresholds of Ṁ > 5 × 10−10 M⊙/yr and L > 0.5 L⊙, respectively.

Fig. 16 Comparison of the NIR flux density (F150W (left panel) or F277W (right panel) filter, aperture from 250 to 1000 au), divided by the square root of the luminosity of the star in solar luminosities, to the accretion rate. Blue dots (Ṁ ≥ 5 × 10−10 M⊙/yr) and cyan triangles (Ṁ < 5 × 10−10 M⊙/yr) are for the class II stars (local density ≤ 5 × 104 cm−3), and the small red squares are for the class I stars (local density > 5 × 104 cm−3). The black solid line is a least-squares fit giving log y = (1.24 ± 0.07) log x − (11.9 ± 0.18) (F150W, left panel) and log y = (1.25 ± 0.05) log x −(11.1 ± 0.11) (F277W, right panel).

5 Discussion

This work provides further evidence that late BH accretion on class II stars can fully replenish their disks. More importantly, thanks to its Lagrangian analysis with tracer particles, it shows that the random nature of BH accretion results in significant time variations of the angular momentum direction of the accreted gas. Such variations may cause angular momentum cancelation in the disk, instabilities and pressure traps, affecting the disk evolution and its accretion rate onto the star. They may also result in observable evidence of disk misalignment, such as misalignment of inner and outer disks (Marino et al. 2015; Pinilla et al. 2018; Ansdell et al. 2020), misaligned single and multi-planet systems (Huber et al. 2013; Sanchis-Ojeda et al. 2013; Kamiaka et al. 2019; Hjorth et al. 2021; Albrecht et al. 2022), distributions of projected stellar rotation velocities of stars with transiting planets (Louden et al. 2021; Esteves et al. 2023), hot Jupiters with orbit obliquity beyond the Kozai angle (Kawai et al. 2024), and possibly misalignment between the disk and the stellar spin in class II stars (too high disk extinction of class II stars with rotation-axis inclination larger than ∼60°, Kounkel et al. 2023). Evidence for the warped or misaligned disks (Sai et al. 2020; Yamato et al. 2023) and misaligned outflows (Okoda et al. 2021; Sai et al. 2024) have also been found in earlier star formation phases.

This scenario of late-time mass infall is consistent with recent discoveries of large-scale flows (streamers) around class II disks (Ginski et al. 2021; Huang et al. 2021; Valdivia-Mena et al. 2024), as well as the general detection of reflection nebulae around class II stars (Gupta et al. 2023). This suggests that targeted observations, such as the JWST observations discussed here, may result in a high detection rate of streamers around class II stars as well. Follow-up ALMA observations may then probe the kinematics of class II streamers, to be compared with accreting filaments in BH trails seen in the simulations. Such large-scale flows have also been found around class 0 stars (e.g., Thieme et al. 2022; Kido et al. 2023; Aso et al. 2023), and around class I stars or stars in transition between class I and class II (e.g., Yen et al. 2019; Pineda et al. 2020; Alves et al. 2020; Grant et al. 2021; Valdivia-Mena et al. 2022; Cacciapuoti et al. 2024; Valdivia-Mena et al. 2024).

We have derived a dependence of j on the stellar mass, as shown in Figure 6 and Eq. (10). As commented in Sect. 3.2, the derived values of j are close to those from our Eulerian analysis (Padoan et al. 2024) and only a factor of two larger than those inferred from the observed disk sizes, which is to be expected because j is most likely not perfectly conserved as material settles to the disk from the BH radius (a factor of four contraction – see Eq. (3)). A correlation between observed disk sizes and stellar mass, for class II stars, qualitatively consistent with our result, was first suggested by Andrews et al. (2018), and later confirmed in other works (e.g., Hendler et al. 2020; Andrews 2020; Long et al. 2022). The reported slope of the dependence of disk size on stellar mass is different in different observational samples, varying from 0.3 in Long et al. (2022) to 0.9 in Andrews (2020). In our sample that combines results from different surveys we find a slope of ∼0.4 ± 0.1.

Because $j~M_{\text{star}}^{1/2}{R}_{\text{d}}^{1/2}$ in a Keplerian disk, our model prediction of Eq. (2), j ∼ M_star, corresponds to R_d ∼ M_star. The result of the Lagrangian analysis of this work corresponds to a slightly steeper relation, with an exponent of 1.28, while the Eulerian analysis in Padoan et al. (2024) results in a slightly shallower exponent of 0.72. At face value, the model prediction is steeper than the observed relation (though comparable in the case of Andrews (2020)). However, the observed relation is affected by a significant observational bias: the fraction of unresolved disks with sizes smaller than ∼20 au in the surveys is very large, ∼70% according to Manara et al. (2023), and clearly skewed toward lower stellar masses (expected to host the smaller disk sizes), so the real mass dependence of the disk sizes should be steeper than what is currently derived from the observations.

We note that a relation between disk size and stellar mass has been recently found also from observations of class 0 and I disks, with a reported slope of 1.1 in Yen et al. (2024), and 0.7 in Yen & Lee (2024). Our fit for embedded stars in Figure 8 corresponds to the relation $R_{\text{d}}~M_{\text{star}}^{0.8}$, consistent with the observations. However, we stress that we have not attempted to predict the value of j for actual class 0 and I disks, because j (or the size) of very young disks must have a significant contribution from the initial collapse of the prestellar core, so it should be primarily determined not by the mass infall at a specific time, which is what we compute here, but by the cumulative infall up to that time. The disk size is in fact expected to grow during the collapse of a prestellar core, though with a shallower mass dependence than observed (an exponent ≤0.5), depending on the initial orientation of the magnetic field (e.g., Hennebelle et al. 2016; Lee et al. 2021). For class II stars, over Myr timescales the Bondi- Hoyle infall (from the ISM to the disks), comparable to the observed accretion rates (from the disks to the stars) as shown in Fig. 2, can fully rejuvenate the disks, so the disk properties inherited from the initial prestellar-core collapse are less relevant.

We find that the angular momentum is not conserved from larger to smaller scales. We should in general expect magnetic breaking to play a role in this MHD simulation, as well as in nature, causing a transfer of angular momentum to larger scales (see for example the MHD disk-formation models by Hennebelle et al. (2016) and Lee et al. (2021)). However, observations of L1527 IRS, a class 0 protostar, by Ohashi et al. (2014), and of L1489 IRS, a class I protostar, by Sai et al. (2020) suggest that the specific angular momentum of infalling materials may be conserved within infalling envelopes (∼1000 au scale). Even if we assumed that angular momentum were conserved from the Bondi-Hoyle radius to the disk radius, this would only scale up the specific angular momenta found in this study by a factor of two, meaning that there would be an even larger angular momentum budget from late-stage infall to explain the observed class II disk sizes.

6 Conclusions

Using a star-formation simulation with a spatial resolution of 25 au, where the full stellar IMF is successfully reproduced, but stellar disks are not resolved, we have studied the angular momentum of the gas gravitationally captured by the stars, with a Lagrangian analysis based on tracer particles. Our main results are summarized in the following:

• BH infall rates of class II stars in the simulation are of the same order as the observed accretion rates of class II stars, suggesting that class II disks may be replenished from larger scale at a similar mass rate as they accrete onto the central star;

• The specific angular momentum of the gas gravitationally captured by class II stars scales approximately linearly with the stellar mass, in agreement with the analytical prediction and numerical analysis in Padoan et al. (2025), and with the observational j − M relation for spatially resolved disks;

• In class II stars, the direction of the accreted angular momentum may vary strongly with time, with a median angle of 47° between consecutive epochs where the accreted mass is of the order of one disk mass. The time evolution is more coherent in the class I phase, with a median angle of 14°;

• Synthetic observations of class II stars in the simulation with accretion rates Ṁ > 5 × 10⁻¹⁰ M_⊙/yr show that their BH accretion flows can be detected by JWST in NIR scattered light. Imposing an additional selection constraint for the stellar luminosity of L > 0.5 L_⊙, the estimated detection success rate with JWST is over 80% with both filters (100% in F150W) and only 2 minutes of exposure time.

With regard to the mass budget, BH accretion in class II stars may alleviate the problems of the short disk accretion times relative to the disk ages (Fedele et al. 2010; Manara et al. 2020; Maucó et al. 2023), and of the insufficient disk masses relative to the masses of exoplanetary systems (e.g., Ansdell et al. 2016; Manara et al. 2018; Mulders et al. 2021; Stefánsson et al. 2023). As for the angular momentum evolution, the contribution of BH accretion in the class II phase cannot be neglected, as it provides a large enough angular momentum to explain the observed size of class II disks, as well as significant time variations to induce strong perturbations in the disks.

Our scenario for class II disks leads to the prediction that, for a high enough accretion rate and stellar luminosity, accretion streamers should be readily detected by targeted JWST observations. If confirmed by future observations, this scenario would have important implications for models of planet formation and disk evolution.

Acknowledgements

We thank the anonymous referee for their helpful comments, leading to an improved and more clear presentation of our work. VMP and PP acknowledge financial support by the grant PID2020-115892GB-I00, funded by MCIN/AEI/10.13039/501100011033 and by the grant CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033. VMP gratefully acknowledges also financial support from the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130). MJ acknowledges support from the Research Council of Finland grant No. 348342. The research leading to these results has received funding from the Independent Research Fund Denmark through grant No. DFF 802100350B (TH). We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France. The astrophysics HPC facility at the University of Copenhagen, supported by research grants from the Carlsberg, Novo, and Villum foundations, was used for carrying out the postprocessing, analysis, and long-term storage of the results.

References

All Figures

	Fig. 1 Column density maps of the simulation from three cardinal directions (x, y, and z). The sink particles are shown with star symbols.
In the text

Fig. 2 Mass accretion rates as a function of stellar mass. Red diamonds are Class II YSOs in Chamaeleon I, Lupus, and L1688 from Testi et al. (2022), and in Taurus from Gangi et al. (2022). Blue points are large-scale infall rates of star particles from the MHD simulation (see Section 2.1) with local gas density lower than 5 × 104 cm−3. The dashed and dotted black lines are the expected chromospheric noise levels (Manara et al. 2013, 2017) using the 3 Myr isochrones of Baraffe et al. (2015) and Feiden (2016), respectively.

In the text

Fig. 3 Power-law slope of jtr versus rtr, $j_{\text{tr}}\propto r_{\text{tr}}^{a_{\text{out},\text{in}}}$, averaged for all the tracer particles within a single epoch, versus the final stellar mass (blue dots), as explained in the main text. The upper panel is for the scaling outside of RBH, the lower panel for the scaling inside RBH. The red lines show the median slope values aout = 0.82 (upper panel) and ain = 0.53 (lower panel, computed over masses ≥0.5 M⊙), and the shaded regions the rms value around the median.

In the text

Fig. 4 Four example stars (one per row) showing the accretion of the last 10% of their mass on the left (the top x-axis is the mass fraction, fM, and the bottom x-axis is the corresponding age of the star), and their late-accretion flows on the right. Left column: Angle ∆ϕ (dashed blue line, “disk-to-disk” misalignment), angle δϕ (dotted blue line, “tracers-to-disk” dispersion), specific angular momentum of the epoch (solid red line), j27 = j × 10−27, and the median accretion rate within 0.5% massbins (dash-dotted red line), Ṁ. Right column: Column density images of log($N_{{\text{H}}_2}$) in a 6400 au box, from three different cardinal directions (xyz), using the middle snapshot of the last epoch for each of the stars. The text in the images list the final mass, the age and the average accretion rate at this snapshot, and the minimum and the maximum limits used for the colorscale of the logarithm images, in all directions.

In the text

Fig. 5 Scatter plot of the formation time of 90% of the last epoch (t999) versus mass on the left, and kernel density estimator (KDE) probability density plots of the formation time in four bins of final stellar mass on the right. We find the local minima of the KDE functions, which are then plotted on the scatter plot as black squares. The black dashed line is a fit to the black squares, log y = 0.72 ± 0.03 log x − 1.54 ± 0.03. This line is used to separate the sinks fed by slower Bondi-Hoyle accretion (blue points) and sinks accreting at a much higher rate, likely still in the class I or even earlier protostellar phases (red points).

In the text

Fig. 6 Scatter plot of specific angular momentum versus mass in the last epoch of the class II stars in the simulation (blue dots), and for observations of spatially resolved disks (Hendler et al. 2020; Cieza et al. 2021; Long et al. 2022; Stapper et al. 2022) (smaller red dots). The blue line is a least-squares fit (log y = (1.14 ± 0.07) log x + (20.7 ± 0.06)) to the blue data points, while the solid black line is the fit (log y = (0.86 ± 0.03) log x + (20.6 ± 0.03)) from the Eulerian analysis in Padoan et al. (2025) and the dashed black line a fit (log y = (0.71 ± 0.04) log x + (20.3 ± 0.02)) to the observations. The gray dotted lines are the specific angular momentum for a given disk radius (1, 100, and 10 000 au, with increasing j) as a function of stellar mass, using Eq. (29) in Padoan et al. (2025) for Keplerian disks.

In the text

	Fig. 7 Histograms of the last ∆ϕ angle (blue step histogram), and of the last δϕ angle (blue-shaded bar histogram), for the class II stars in the simulation (see the main text for the definition of the angles). The median values of ∆ϕ and δϕ are shown by the black dashed and dotted vertical lines, respectively.
In the text

Fig. 8 Scatter plot of the specific angular momentum versus mass in the last epoch of the class I and earlier protostars. The red line is the leastsquares fit (log y = (0.90 ± 0.06) log x+ (20.8 ± 0.04)) to the data points, while the blue line is the fit to the simulated class II stars from Figure 6. The gray dotted lines are the specific angular momentum for a given disk radius (1, 100, and 10 000 au, with increasing j) as a function of stellar mass, using Eq. (29) in Padoan et al. (2025) for Keplerian disks.

In the text

	Fig. 9 Histograms of the last ∆ϕ (red step histogram) and of the last δϕ (red-shaded bar histogram), for class I and earlier protostars. The median values of ∆ϕ and δϕ are shown by the black dashed and dotted vertical lines, respectively. For comparison, the corresponding histograms for class II stars from Fig. 7 have been replotted here in blue.
In the text

Fig. 10 Hertzsprung-Russell diagram for all stars in the last snapshot (top panel) and for the 69 stars in the JWST synthetic observations (bottom panel). The small symbols (magenta and red) are for photospheric luminosities, while the large symbols (green and blue) are for stellar luminosities which include the accretion luminosity. Magenta and green points are class I stars, while red and blue points are class II. When Lacc is small (especially class II stars), the small points overlap with the larger ones, resulting in the appearance of red/magenta points with blue/green borders for class II and class I stars, respectively. The zeroage main sequence is indicated by the dashed black line, while the thin transparent light-green lines show the evolutionary tracks of all stars from an age of 104 yr to the time of the snapshot, and the thin transparent light-blue lines are the evolutionary tracks with episodic accretion included. Stellar evolution is computed following the methodology of Jensen & Haugbølle (2018).

In the text

Fig. 11 JWST synthetic NIR observations (see Section 4) around a relatively high-mass star (located at the center of the map) in a 6400 × 6400 au frame. The title lists the sink particle number and direction of projection, the age, the mass, the luminosity and the accretion rate. The top row shows, in log scale the surface brightness in the two selected wavelengths (left and central images respectively, with the JWST filter as a label), and the local column density integrated across the 12 800 au depth of the selected volume (right image, with minimum and maximum log value of the color scale as a label). The contours (white, red, and yellow) are 5σ limits for 2, 8 and 50 min exposure time, respectively (the column density image shows the F277W contours). The bottom row shows the corresponding surface brightness images after the ring subtraction described in the text (the column density map is the same as in the upper row). Negative and zero pixel values are masked out (white background). Contours show the 5-sigma limits of the remaining surface brightness, after the ring subtraction.

In the text

	Fig. 12 As Fig. 11, but for a lower-mass star.
In the text

Fig. 13 Total number of maps in bins of the accretion rate, and the detection fractions of the maps where accretion features were detected with S/N ≥5σ at least 1000 au from the central star. The black line is the number of all the maps, non-detections included, with y-axis on the left. The colored lines (blue and red for F150W and F277W filters, and dashed and dotted for 2 min and 50 min exposure times, respectively) and the y-axis on the right show the successful detection fraction of the maps.

In the text

	Fig. 14 As Fig. 13, but in bins of the stellar luminosity instead of the accretion rate.
In the text

Fig. 15 Scatter plots of stellar luminosity versus accretion rate for all the synthetic observations, with symbols and colors referring to the detection rate after 2 min of exposure time with the F150W (left panel) and F277W (right panel) filters. The small black dots correspond to a non-detection of an accretion feature in all three directions (x, y, and z), the magenta triangles detection in only one direction, the green squares detection in two directions, and the large blue dots detection in all three directions. The protostars (class I or earlier) are highlighted by red circles around the markers. The red dotted vertical and horizontal lines mark the thresholds of Ṁ > 5 × 10−10 M⊙/yr and L > 0.5 L⊙, respectively.

In the text

Fig. 16 Comparison of the NIR flux density (F150W (left panel) or F277W (right panel) filter, aperture from 250 to 1000 au), divided by the square root of the luminosity of the star in solar luminosities, to the accretion rate. Blue dots (Ṁ ≥ 5 × 10−10 M⊙/yr) and cyan triangles (Ṁ < 5 × 10−10 M⊙/yr) are for the class II stars (local density ≤ 5 × 104 cm−3), and the small red squares are for the class I stars (local density > 5 × 104 cm−3). The black solid line is a least-squares fit giving log y = (1.24 ± 0.07) log x − (11.9 ± 0.18) (F150W, left panel) and log y = (1.25 ± 0.05) log x −(11.1 ± 0.11) (F277W, right panel).

In the text