EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) A97 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number A97
Number of page(s) 15
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202558019
Published online 02 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

Recent advances in far-infrared/submillimeter observations (ALMA, NOEMA, and VLA) and magnetohydrodynamic (MHD) simulations have improved our understanding of protoplanetary disks. However, key gaps remain, particularly in the Class 0 phase where disks, which are deeply embedded in their nascent cores, are challenging to observe and simulate due to highly nonlinear processes acting across a large dynamical range (see the reviews by Teyssier & Commerçon 2019; Tsukamoto et al. 2023; Lesur et al. 2023). The magnetically regulated disk-formation scenario (Hennebelle et al. 2016a; Wurster et al. 2018; Vaytet et al. 2018; Wurster & Lewis 2020a,b; Lee et al. 2021; Lebreuilly et al. 2024; Mayer et al. 2025a; Lee et al. 2024), in which the envelope’s magnetic field extracts angular momentum and ambipolar diffusion regulates magnetic braking, is the favored theoretical model. This produces disk radii of ≲50 AU, consistent with the observations of Maury et al. (2019) and Tobin et al. (2020)1. Said models find that ∼0.1 G field strengths are needed to reproduce observations, which aligns with paleomagnetic measurements of Solar System meteorites (Weiss et al. 2021). Once formed, Class 0 disks experience vigorous accretion in the form of infall that shapes their dynamical evolution. Understanding this phase is crucial for planet formation scenarios, as current surveys find Class II disk masses insufficient to explain exoplanetary systems (Greaves & Rice 2010; Williams 2012; Najita & Kenyon 2014; Mulders et al. 2015, 2018). This implies either early planet formation during the Class 0 phase or continuous mass replenishment from the surroundings (Manara et al. 2018; Padoan et al. 2025). In addition, cosmochemical analyses of meteorites have provided constraints on disk kinematics in the early Solar System (Morbidelli et al. 2024), indicating that calcium aluminum inclusions (CAIs) formed near the Sun and were transported to larger distances, where they were incorporated into carbonaceous chondrites. Modeling this early phase by describing accretion processes and their kinematic impacts on the disk is therefore essential to bridging the gap between star and planet formation models.

Several studies have investigated disk formation during gravitational collapse. Lee et al. (2021) simulated the collapse of a 1 M⊙ cloud using radiative MHD (including flux-limited diffusion and ambipolar diffusion); they found that magnetic fields dominate infall dynamics and that the disk midplane contributes less to mass transport than the upper layers. Similarly, Tu et al. (2024) studied a magnetized collapse of a 0.5 M⊙ cloud core with ambipolar diffusion, albeit by adopting a barotropic equation of state. They found that gravo-magnetic sheets feed the disk primarily through its upper and lower surfaces. Xu & Kunz (2021) similarly carried out a collapse with ambipolar and Ohmic diffusion, albeit in a non-turbulent core, and found that angular momentum transport within the disk is driven by gravitational instabilities. On larger scales, Mayer et al. (2025b) and Yang & Federrath (2025) performed zoomed-in barotropic MHD simulations from cloud-scale initial conditions. While Yang & Federrath (2025) used ideal MHD, Mayer et al. (2025b) included ambipolar diffusion. Both found that magnetic fields and inherited cloud-scale turbulence channel accretion into specific regions in the disk. On the high-mass end, He & Ricotti (2023, 2025) carried out similar zoom-in simulations that included ideal MHD and radiative transfer, highlighting how turbulence mitigates magnetic braking. Across both low- and high-mass regimes, anisotropic accretion has been consistently reported and has been further connected to the development of eccentric motions within disks (Commerçon et al. 2024; Calcino et al. 2025). In addition, external accretion can drive turbulence in the disk, thereby transporting angular momentum (Lesur et al. 2015; Hennebelle et al. 2016b, 2017; Kuznetsova et al. 2022).

In this paper we present a comprehensive analysis of the accretion and kinematics of protoplanetary disks formed through self-consistent radiative-MHD simulations of gravitational collapse in an isolated dense core. We characterize how anisotropic accretion shapes both the global disk kinematics and the previously unexplored dynamics of eccentric fluid motions, an important yet overlooked aspect of early disk evolution.

2. Model

We modeled the collapse of an isolated uniform-density core of mass M0 using the RAMSES adaptive mesh refinement code (Teyssier 2002). The core was initialized with a temperature of T0 = 10 K and size R0 (thermal-to-gravitational energy ratio α = 0.4), solid-body rotation along the −z axis (rotational-to-gravitational energy ratio β = 4 × 10−2), and a uniform magnetic field tilted by 10° relative to the rotation axis, with a mass-to-flux ratio μ = (10/3)μcrit. No initial turbulent velocity field was present in this run; instead, we introduced an m = 2 azimuthal density perturbation of amplitude 10%. In this setup, the magnetic interchange instability develops spontaneously during the collapse2.

The collapse proceeds self-consistently until the density exceeds 5 × 1012 cm−3, at which point a sink particle forms (Bleuler & Teyssier 2014), accreting material within a radius of 4Δxmin (where Δxmin is the resolution at the finest adaptive mesh refinement level). The sink particle then accretes 10% of the material whose density exceeds 5 × 1012 cm−3 (see Hennebelle et al. 2020a). The simulation begins with a 643 coarse grid (ℓmin = 6) and is refined to maintain 20 cells per Jeans length, computed with a temperature cap of 100 K. The simulation domain is ≈16 000 AU3 for R1 and ≈48 000 AU3 for R2.

Our radiative-MHD implementation uses a hybrid solver (Mignon-Risse et al. 2021a,b), combining the M1 method for stellar irradiation from sink particles (Rosdahl et al. 2013; Rosdahl & Teyssier 2015) with gray flux-limited diffusion for absorption and radiative transfer (Commerçon et al. 2011, 2014). Nonideal MHD effects are present in the form of ambipolar diffusion (Fromang et al. 2006; Teyssier et al. 2006; Masson et al. 2012).

Two runs are presented in this paper: R1 and R2. R1 has M0 = 1 M⊙ and ℓmax = 14 (Δxmin ≈ 0.97 AU). R2 has M0 = 3 M⊙ and ℓmax = 16 (Δxmin ≈ 0.72 AU). Assuming a ≈30% star formation efficiency (André et al. 2010; Könyves et al. 2015), these should form a 0.3 M⊙ and 1 M⊙ star, respectively. R1 was previously presented in Commerçon et al. (2024) and is similar to Hennebelle et al. (2020a), Lee et al. (2021). The spatial resolution of both runs is traded for extended integration time, enabling us to study the long-term evolution of the disk.

We also accounted for dust dynamics under the monofluid approximation (Lebreuilly et al. 2019; Commerçon et al. 2023). We considered 𝒩d bins of single-size dust species with sizes ranging from 10−3 μm to 20 μm. Dust particles do not undergo coagulation or fragmentation. We used 𝒩d = 10 for R1 and 𝒩d = 5 for R2.

Finally, 104 massless Lagrangian particles were used in run R2 to trace the gas dynamics during the collapse using the particle-mesh method, allowing us to track the temperature-density history of gas parcels during the collapse. These were randomly distributed throughout the volume of the initial dense core. Run R1 does not contain any tracer particles.

3. Global disk and ejecta properties

We began by presenting the qualitative and quantitative properties of the disks formed in runs R1 and R2. Figure 1 shows column density images at t = 58.5 kyr for R1 and t ≈ 27 kyr for R2, with t = 0 corresponding to the epoch of sink formation. In both runs, strong radial magnetic field gradients trigger the interchange instability (Spruit et al. 1995), visible at large scales (Figs. 1a, d) as distinct streamers feeding the disks. This highlights the magnetic field’s dominant role in governing gas kinematics during the collapse and subsequent disk accretion. The instability is noticeably more pronounced in R2 due to its greater mass reservoir. Edge-on views (Figs. 1c, f) reveal numerous streamers connecting to the disk, demonstrating anisotropic accretion. As discussed in Sect. 5.1, this anisotropic accretion directly drives eccentric motions within the disk (see also Commerçon et al. 2024; Calcino et al. 2025).

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Column density images in the z (panels a, b, d, and e) and x (c and f) directions for runs R1 (first row) and R2 (second row), at both large scales (a and d) and small scales (b, c, e, and f). The lime colored contours in panels (b) and (e) are a least-square fit of an ellipse enveloping the disk’s surface. The associated eccentricity is shown in the top-right corner of panels (b) and (e).

Figures 1b and 1e show an ellipse fitted to the disk’s outer edge3, from which we measured the geometric eccentricity via 1 − b 2 a 2 Mathematical equation: $ \sqrt{1-\frac{b^2}{a^2}} $ (a and b being the major and minor axes). The disk exhibits significant eccentricity, with values of 0.33 for R1 and 0.26 for R2, with streamers feeding its outer regions anisotropically4. Figure 2 tracks the temporal evolution post-sink formation (which marks our t = 0)5. The two runs show qualitatively similar protostellar mass growth (M*, Fig. 2a), though R2’s rate is higher. The disk mass (Md, Fig. 2b) differs markedly: R1 declines steadily from ≈2 × 10−1 M⊙ to ≈3 × 10−2 M⊙ over 60 kyr, while R2 shows large fluctuations driven by the interchange instability. The disk radius Rd (90% mass radius as done in Mignon-Risse et al. 2021b, solid lines in Fig. 2b) generally follows Md’s behavior. For R1, Rd increases despite mass loss, which is consistent with viscous disk theory where outward angular momentum transport causes inward mass transport6. The semimajor axis ad obtained from elliptical fitting (Fig. 2b) mirrors the trends in Rd but is systematically larger, as it traces the outermost orbits. The geometric eccentricity (Fig. 2c) fluctuates markedly (shaded regions) due to dynamical instabilities (spiral waves, accretion streamers) altering the disk’s morphology. However, the time-averaged eccentricity (solid lines) remains stable (> 0.2) with smaller variations in R1 than R2.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Global evolution of the protoplanetary disk as a function of time, where t = 0 represents the epoch of sink (i.e., protostellar) formation. The quantities shown are the protostellar mass (solid lines in panel a), disk mass (dotted lines in panel a), disk radius (solid lines in panel b), disk semimajor axis (dotted lines in panel b), and apparent disk eccentricity (c). The disk semimajor axis and apparent eccentricity are inferred from an elliptical fit. The blue (red) curve corresponds to the 1 M⊙ (3 M⊙) run R1 (R2). The shaded regions in panel (d) represent temporal fluctuations in the measurement, and the solid line is an average value.

Overall, R2 seems to display more dynamical behavior over time than R1. We link R2’s large dip in Rd at t ≈ 15 kyr with an influx of material with opposite angular momentum to that of the disk in Appendix B.

4. Disk accretion

The apparent eccentricities, morphological variability, and anisotropic accretion revealed in Figs. 1 and 2 naturally raise the question of how material is actually delivered to and redistributed within the disk. In this section we examine the accretion process in detail by quantifying both its directional dependence and the relative contributions of vertical and radial infall. All analyses here use a reference frame centered on the sink particle, with the disk’s angular momentum defining the z-axis.

4.1. Directional mass fluxes

We quantified the mass carried by the magnetically driven streamers shown in Fig. 1 by analyzing mass fluxes across a fixed Rshell = 30 AU sphere (exceeding both runs’ maximum semimajor axes) around the sink particle. This approach avoids the complexities of tracking the disk’s irregular surface (Fig. 2c). We computed mass fluxes through spherical surface elements via

M ˙ ( Φ , Λ , t ) = − ρ v r d S = − ρ v r R shell 2 cos ( Φ ) d Φ d Λ , Mathematical equation: $$ \begin{aligned} \dot{M}(\Phi , \Lambda , t) = -\rho v_{\mathrm{r} } \mathrm{d}S = -\rho v_{\mathrm{r} } R_{\mathrm{shell} }^{2} \cos (\Phi ) \mathrm{d}\Phi \mathrm{d}\Lambda , \end{aligned} $$(1)

where Φ and Λ are the latitude and longitude, ρ is the gas density, and vr is the radial velocity in spherical coordinates.

The time-integrated Hammer projections of Fig. 3 reveal anisotropic accretion patterns, with concentrated midplane lobes corresponding to streamer-disk connections. The polar regions show low-velocity outflows (negative fluxes) of magnetic origin, carrying ∼100 times less mass than the inflows (positive fluxes; Fig. 4a), predominantly within the first ∼5 kyr. Run R1’s northern outflow exhibits bifurcated lobes due to a reorientation of the outflow during the simulation.

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Hammer projections displaying the time and surface integrated radial mass flux on a sphere of Rshell = 30 AU for runs R1 (left) and R2 (right), displaying both the inflow (top colorbar) and outflow (bottom colorbar) of material throughout the simulation’s duration post-sink formation (tf), computed by integrating Eq. (1) in time. The lime colored contour delimitates the transition from positive to negative values.

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Total inflow (solid lines) and outflow (dotted lines) of mass (panel a) and specific angular momentum (panel b) through a sphere of radius Rshell = 30 AU around the sink particle as a function of time, where t = 0 corresponds to the epoch of sink formation, for runs R1 (blue) and R2 (red). This figure is complementary to Fig. 3.

We also compare in Fig. 4b the specific angular momentum jshell carried by the flow of material, measured as

j shell = − R shell ∫ t ∫ S ρ v r v ϕ d S d t ∫ t ∫ S M ˙ d t , Mathematical equation: $$ \begin{aligned} j_{\mathrm{shell} } = -\frac{R_{\mathrm{shell} }\int _{t} \int _{S} \rho v_{\mathrm{r} }v_{\phi }\mathrm{d}S\mathrm{d}t}{\int _{t}\int _{S}\dot{M}\mathrm{d}t}, \end{aligned} $$(2)

where vϕ is the azimuthal velocity. We see that outflows carry 3.85 times more specific angular momentum than the inflows for R1, and 3.41 times for R2. This illustrates the role of magnetic fields in carrying angular momentum outward. However, other transport processes within the disk must be at play to remove the remaining angular momentum, as the outflow does not have sufficient mass to be efficient in this regard.

4.2. Vertical and radial accretion

Once material crosses the inner 30 AU radius, we examined how it reaches the disk by computing mass fluxes across a fixed cylindrical surface with Rcyl = 30 AU and height hcyl = 16 AU. This enabled a comparison of the radial inflow through the cylinder wall and vertical inflow through its top and bottom surfaces at z = ±8 AU. We emphasize however, that this cylinder does not represent the disk surface, but rather acts as its proxy. The radial mass accretion rate is (locally)

M ˙ cyl ( z , ϕ , t ) = − ρ v R R cyl d ϕ d z , Mathematical equation: $$ \begin{aligned} \dot{M}_{\mathrm{cyl} }(z, \phi , t) = -\rho v_{\mathrm{R} }R_{\mathrm{cyl} }\mathrm{d}\phi \mathrm{d}z, \end{aligned} $$(3)

where vR is the radial velocity in cylindrical coordinates and ϕ the azimuth. The vertical accretion rate through the top and bottom surfaces is (locally)

M ˙ z ( R , ϕ , t ) = − sgn ( z ) ρ v z R d R d ϕ , Mathematical equation: $$ \begin{aligned} \dot{M}_{z}(R, \phi , t) = -\mathrm{sgn}(z)\rho v_{z}R \mathrm{d}R \mathrm{d}\phi , \end{aligned} $$(4)

where sgn is the sign function, vz the vertical velocity, and R = x 2 + y 2 Mathematical equation: $ R=\sqrt{x^2+y^2} $.

Figure 5 shows the time integrated mass fluxes. Radial accretion is anisotropic in ϕ but less so across z. Vertical accretion at z = ±8 AU (second and third columns) reveals negative values in R1 from current sheets that vertically recirculate material7, with mass reinjected elsewhere radially or through the opposite face. Red lines mark the minimum and maximum disk radius over time. Both runs show substantial accretion in the inner disk, with some material landing at radii as small as ∼10 AU.

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Total inflow of material through a cylinder of radius Rshell = 30 AU and height h = 16 AU throughout the simulation’s duration post-sink formation for runs R1 (≈58 kyr, first row) and R2 (≈27 kyr, second row), computed by integrating Eq. (3) (first column) and Eq. (4) (second and third columns) in time. The red circles in the second and third columns correspond to the minimum and maximum semimajor axis of the disk throughout each simulation’s duration, and thus represent the locations where material has landed in the disk.

In Fig. 6, the azimuthally and temporally integrated vertical accretion ∫ t ∫ ϕ M ˙ z d t Mathematical equation: $ \int_{t}\int_{\phi}\dot{M}_{z}\,\mathrm{d}t $ indicates that vertical inflow delivers mass fairly evenly across disk radii8. Figures 7a and 7b compare the total and surface-specific accretion from radial (solid lines) and vertical (dash-dotted lines) directions. In R1, most mass arrives vertically, while in R2, radial and vertical contributions are similar. Nevertheless, radial accretion dominates per unit area over vertical accretion in both runs. R2’s enhanced radial inflow stems from its greater angular momentum budget, which concentrates more material near the midplane (see Figs. 1c and 1f). In contrast with classical 1D hydrodynamic models, where material falls at its centrifugal radius, our simulations reveal that magnetic fields, in concert with turbulence, play an important role by rendering accretion highly anisotropic.

Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Vertical mass accretion onto the disk for runs R1 (blue) and R2 (red), obtained by integrating Eq. (4) over ϕ and time. Dash-dotted and dotted lines show accretion through the top (z = 8 AU) and bottom (z = −8 AU) surfaces, respectively; the solid lines are their sum. Solid vertical lines mark the disk’s minimum and maximum radii during the simulation. The dashed lines represent the cumulative fraction of landed mass (normalized to unity), plotted against the right-hand axis. This figure is complementary to Fig. 5.

Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Comparison of the total (panel a) and surface-specific (panel b) mass accreted through a cylinder of radius Rcyl = 30 AU (solid lines) and the top and bottom disks of the cylinder at, respectively, z = 8 AU and z = −8 AU (dash-dotted lines). These were computed by integrating Eqs. (3) and (4) in space and time.

4.3. Tracer particles

Run R2 employs 104 gas tracer particles, of which 5645 are located within the disk or the sink accretion volume by the final simulation snapshot. Figure 8 shows the trajectories of these particles in the x − y (top) and z − x (bottom) planes. These trajectories are consistent with our earlier results, confirming that the majority of the accreted mass arrives onto the disk near the equatorial plane. In the x − y view, particles clearly trace the prominent interchange instability, with distinct streamers feeding material into the disk. These trajectories arise naturally from the particle–mesh method, which drives tracers toward regions of higher density9. The z − x view highlights the highly irregular and chaotic nature of vertical accretion, with a wide range of particle trajectories in this direction. Unfortunately, once the particles have settled in the disk, our snapshot frequency does not permit us to follow their evolution within it. In addition, most of these particles are ultimately accreted by the sink particle early in the disk’s history, as they migrate inward during phases of inward-driven transport (see Sect. 5.2) and can no longer escape its immediate vicinity.

Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

Trajectories of the 5645 tracer particles that settle into the disk during run R2 in the x − y (top) and x − z (bottom) planes. The color scale reflects the cumulative distance traveled by each particle, with darker shades marking the most recent segments of each path.

Having established the spatial trajectories of the tracers, we turned to their thermodynamic histories during infall. Figure 9 shows the temperature–density histories of randomly selected tracer particles that settle onto the disk either from the polar regions or from near the equatorial plane. Particles accreted through the equatorial plane (Fig. 9b) exhibit a more classical adiabatic heating (gray curve) trajectory prior to disk settling, with a sharp density increase once they settle onto the disk. By contrast, particles arriving from the polar regions (Fig. 9b) follow more complex evolutionary paths, undergoing repeated cooling and heating events prior to disk settling as they traverse different regions of the flow, being influenced by both the magnetic field and the outflow.

Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

Temperature–density histories of randomly selected tracer particles that settle onto the disk from the polar regions (panel a) and from near the equatorial plane (panel b). Each curve corresponds to an individual particle, represented by a different color. The gray curve in panel (b) represents T ∝ ρ5/3 − 1.

5. Disk kinematics

The anisotropic nature of accretion imprints itself on the disk’s velocity field. Rather than settling into purely circular Keplerian orbits, the localized streams of infalling gas produce eccentric motions (Commerçon et al. 2024), and as we show below, drive mass transport in the interior. In this section we explore these kinematic features and how they shape the disk’s evolution.

5.1. Eccentric motions within the disk

As noted earlier, the disks in our simulations exhibit ellipsoidal morphologies, with apparent eccentricities often exceeding ∼0.1 (Fig. 2c). We quantified eccentric motions via the fluid kinematics following Commerçon et al. (2024):

e = v × j G M ∗ − r r , Mathematical equation: $$ \begin{aligned} \boldsymbol{e} = \frac{{\boldsymbol{v}}\times {\boldsymbol{j}}}{GM_{*}} - \frac{{\boldsymbol{r}}}{r}, \end{aligned} $$(5)

where e is the eccentricity vector and r = x 2 + y 2 + z 2 Mathematical equation: $ r=\sqrt{x^{2}+y^{2}+z^{2}} $. M* is the sink (i.e., protostellar) mass. The semimajor axis of each orbit is then obtained as

a = ( 2 r − v 2 G M ∗ ) − 1 . Mathematical equation: $$ \begin{aligned} a = \left(\frac{2}{r}-\frac{v^{2}}{GM_{*}}\right)^{-1}. \end{aligned} $$(6)

To account for pressure contributions on the orbital velocity, we measured the average eccentricity within an orbital bin. The resulting values are presented in Fig. 10. We find that disk eccentricity is progressively damped over time. In run R2 (Fig. 10b), episodes of enhanced accretion (Fig. 2a) coincide with temporary eccentricity growth, followed by rapid damping. Elevated eccentricities persist in the inner disk, probably from the influence of the sink particle. To assess this, we computed the angular momentum deficit (AMD) within the disk, which represents the difference between the circular orbit angular momentum and the actual orbital angular momentum:

AMD = ∫ a ∫ 0 2 π Σ ( a ) a 3 G M ∗ ( 1 − 1 − e 2 ) d a d ϕ , Mathematical equation: $$ \begin{aligned} \mathrm{AMD} = \int _{a}\int _{0}^{2\pi }\Sigma (a)\sqrt{a^{3}GM_{*}}\left(1 - \sqrt{1-e^{2}}\right) \mathrm{d}a \mathrm{d}\phi , \end{aligned} $$(7)

Thumbnail: Fig. 10. Refer to the following caption and surrounding text. Fig. 10.

Distribution of eccentricity, computed from the disk kinematics (Eqs. 5 and 6), for run R1 (a) and run R2 (b). Each curve represents a different time, with t = 0 marking the epoch of sink formation. The dashed line indicates the largest radius enclosing 90% of the disk mass over the entire simulation. The sink accretion radius is indicated by the dotted line.

where Σ(a) is the mean column density within an orbital bin. In the limit of an ideal, inviscid disk, AMD is a conserved quantity. The resulting specific AMD (Fig. 11) shows a general decrease over time for R1, whereas the trend is not as clear for R2. It also shows increases during strong accretion episodes. We attribute the initial (t < 10 kyr) declining trend to grid dissipation, as it is stronger when the disk is small and the orbital resolution is inadequate. The declining trend weakens for R1 when the disk grows in radius and numerical resolution becomes sufficient. In the quiescent regime of run R1 (t > 35 kyr), we can estimate a characteristic damping timescale under the assumption that the disk turbulence acts as an effective bulk viscosity (Goodchild & Ogilvie 2006). This timescale is comparable to the viscous timescale, approximated as (αΩ)−1(R/H)2. Based on the disk properties inferred from Fig. 2b (M* ∼ 0.3 M⊙, Rd ∼ 15 AU, H/R ∼ 0.1) and the effective viscosity α ∼ 0.110, we estimate this timescale to be ∼17 kyr. While the exact dissipative properties of the turbulence remain complex (Ogilvie 2001), this estimate highlights that physical damping timescales are short enough to be relevant within the remaining 25 kyr simulation time. Thus, the damping we observe may not be purely numerical. Consequently, the sustained eccentricity implies that anisotropic accretion actively replenishes AMD against these dissipative processes.

Thumbnail: Fig. 11. Refer to the following caption and surrounding text. Fig. 11.

Specific AMD within the disk as a function of time since sink formation for runs R1 (blue) and R2 (red).

While we observe significant eccentricity in the Class 0 phase driven by anisotropic infall, this state is likely transient. The existence of the Cold Classical population in the Kuiper Belt, which formed during the later Class II phase and exhibits quasi-circular orbits, suggests that the disk eventually circularizes. This provides further evidence that AMD is not conserved; once accretion subsides, damping reduces the eccentricity, leading to the dynamically cold disks in which planetary bodies eventually form.

5.2. Orbital mass transport

Eccentricity is only one aspect of the disks’ dynamics, as the inflow patterns seen in Fig. 6 also suggest vigorous vertical mass accretion across all radii. Such accretion can drive turbulence (Klessen & Hennebelle 2010), which in turn causes angular momentum transport through the turbulent stress tensor (Balbus & Papaloizou 1999).

Assessing mass transport in cylindrical or spherical coordinates is challenging, as eccentric fluid motions bias the measured fluxes. To account for the disk’s nested eccentric orbits, we adopted the orbital formalism of Ogilvie & Barker (2014), which uses (λ, ϕ, z) coordinates, where the semi-latus rectum

λ = a ( 1 − e 2 ) = j 2 G M ∗ Mathematical equation: $$ \begin{aligned} \lambda = a(1-e^2) = \frac{j^2}{GM_{*}} \end{aligned} $$(8)

is conserved within an orbit if angular momentum is conserved, making it a convenient quasi-radial coordinate. We therefore measured mass transport in λ bins via

M ˙ ( λ ) = ∫ ∫ J ρ ( v r R λ − R ϕ R λ Ω ) d ϕ d z , Mathematical equation: $$ \begin{aligned} \dot{M}(\lambda ) = \int \int J\rho \left(\frac{v_{\mathrm{r} }}{R_{\lambda }} - \frac{R_{\phi }}{R_{\lambda }}\Omega \right) \mathrm{d}\phi \mathrm{d}z, \end{aligned} $$(9)

where J = ∂(x, y, z)/∂(λ, ϕ, z) is the Jacobian of the coordinate system (see Ogilvie & Barker 2014)11.

Turbulent stresses will induce transport

F Rey ( λ ) = − ∂ G / ∂ λ d j / d λ , Mathematical equation: $$ \begin{aligned} \mathcal{F} _{\mathrm{Rey} }(\lambda ) = -\frac{\partial \mathcal{G} /\partial \lambda }{\mathrm{d}j/\mathrm{d}\lambda }, \end{aligned} $$(10)

where the subscript stands for “Reynolds” and 𝒢 is the internal torque computed using

G = − ∫ ∫ J R 2 T λ ϕ d ϕ d z , Mathematical equation: $$ \begin{aligned} \mathcal{G} = -\int \int JR^{2}T^{\lambda \phi }\mathrm{d}\phi \mathrm{d}z, \end{aligned} $$(11)

and Tλϕ is a component of the turbulent stress tensor

T λ ϕ = R J T R ϕ − R R ϕ J T ϕ ϕ , Mathematical equation: $$ \begin{aligned} T^{\lambda \phi } = \frac{R}{J} T^{R\phi }-\frac{RR_{\phi }}{J}T^{\phi \phi }, \end{aligned} $$(12)

where

T R ϕ = δ v r δ v ϕ , Mathematical equation: $$ \begin{aligned}&T^{R\phi } = \delta v_{\mathrm{r} }\delta v_{\phi }, \end{aligned} $$(13)

T ϕ ϕ = δ v ϕ 2 . Mathematical equation: $$ \begin{aligned}&T^{\phi \phi } = \delta v_{\phi }^{2}. \end{aligned} $$(14)

Here, δvr and δvϕ are deviations from Keplerian orbital velocity in the radial and azimuthal components, respectively12.

Magnetic fields also induce transport through the Maxwell stress tensor:

T M R ϕ = v A , R v A , ϕ , Mathematical equation: $$ \begin{aligned}&T^{R\phi }_{\mathrm{M} } = v_{\mathrm{A,R} }v_{\rm A,\phi }, \end{aligned} $$(15)

T M ϕ ϕ = v A , ϕ 2 , Mathematical equation: $$ \begin{aligned}&T^{\phi \phi }_{\mathrm{M} } = v_{\rm A,\phi }^{2}, \end{aligned} $$(16)

where v A = B / 4 π ρ Mathematical equation: $ \boldsymbol{v}_{\mathrm{A}} = \boldsymbol{B}/\sqrt{4\pi\rho} $ is the Alfvén speed. These components yield a Maxwell torque and corresponding internal mass flux ℱM(λ). Finally, we measured the gravitational stress tensor-induced mass flux (ℱ𝒢) using

T G R ϕ = v g , R v g , ϕ , Mathematical equation: $$ \begin{aligned}&T^{R\phi }_{{\mathcal{G} }} = v_{\mathrm{g,R} }v_{\mathrm{g} ,\phi }, \end{aligned} $$(17)

T G ϕ ϕ = v g , ϕ 2 , Mathematical equation: $$ \begin{aligned}&T^{\phi \phi }_{{\mathcal{G} }} = v_{\mathrm{g} ,\phi }^{2}, \end{aligned} $$(18)

where v g = − g / 4 π G ρ Mathematical equation: $ {\boldsymbol{v}}_{\boldsymbol{g}} = -\boldsymbol{g}/\sqrt{4\pi G\rho} $ and g is the gravitational acceleration. Note here that positive values of Ṁ, ℱRey, ℱM, and ℱ𝒢 correspond to decretion, i.e., outward mass transport in the λ direction.

As accretion is highly time-dependent and non-monotonic, measuring mass transport rates from individual snapshots provides limited insight into the cumulative contributions of different terms over the simulation. We therefore performed a statistical analysis of the orbital mass transport. Fixed logarithmic bins in λ are defined, within which we computed Ṁ(λ), ℱRey, ℱM, and ℱ𝒢 for each snapshot using only cells belonging to the disk. All quantities were vertically averaged with density weighting, thus biasing the results with the dense midplane. By this stage of the simulation, the maximum refinement level is reached and time between snapshots becomes roughly constant. The statistical distributions are shown in Fig. 12 (for Ṁ(λ), ℱRey, and ℱM) and Fig. 13 (for ℱ𝒢), where we plot the zeroth (Q0, minimum, bottom gray curves), first (Q1, 25th percentile, blue curves), second (Q2, median, red curves), third (Q3, 75th percentile, black curves), and fourth (Q4, maximum, top gray curves) quantiles for each λ bin, for run R1 (top rows) and R2 (bottom rows). The innermost regions strongly influenced by the sink particle are excluded from the plot.

Thumbnail: Fig. 12. Refer to the following caption and surrounding text. Fig. 12.

Internal kinematics of the disk. Shown here are the orbital mass accretion rate (a and d), turbulent stress tensor induced accretion rate (b and e), and Maxwell stress tensor induced accretion rate (c and f) for runs R1 (first row) and R2 (second row), computed as described in Sect. 5.2. These are computed in fixed logarithmic bins in λ throughout all simulation snapshots using only cells belonging to the disk, and we display the resulting first (blue), second (red), and third (black) quantiles, which respectively represent the 25th, 50th, and 75th percentiles. The gray curves in each plot represent minimal and maximal values.

Thumbnail: Fig. 13. Refer to the following caption and surrounding text. Fig. 13.

Same as Fig. 12 but this time showing the gravitational stress tensor-induced accretion rate (ℱ𝒢).

The orbital mass transport profiles in Figs. 12a and 12d show median Ṁ values of ∼10−5 M⊙ yr−1 for both runs, indicating that material mainly spreads outward in the midplane. Beyond λ ≈ 14 AU for R1 and λ ≈ 18 AU for R2, radial motions are negligible, as the magnetic field truncates the disk (Lee et al. 2021). However, Q1 curves (blue) reveal intermittent inward transport, with a larger spread in R2 due to its higher disk mass. The turbulent stress tensor (Figs. 12b and 12e) drives outward transport in the outer disk (λ ≳ 8 AU for R1; λ ≳ 6 AU for R2) and inward transport in the inner regions, with median rates of ∼10−5 M⊙ yr−1 for R1 and ∼10−4 M⊙ yr−1 for R2. By contrast, the Maxwell stress (Figs. 12c and 12f) consistently drives inward transport, but at much lower amplitudes (median rates of ∼10−6 M⊙ yr−1 for R1, ∼10−5 M⊙ yr−1 for R2), rendering it largely negligible compared to turbulent stresses. ℱ𝒢 displays more complex behavior in run R1 (Fig. 13a), with median accretion rates oscillating between positive and negative values, albeit at negligible amplitudes when compared to turbulent and Maxwell stress. The more dynamical run R2 has ℱ𝒢 display negative median rates (Fig. 13b), indicating inward transport, albeit once again at negligible rates (∼10−5 M⊙ yr−1) when compared to turbulent stresses. The disks are gravitationally stable (Appendix D) as they are hot, which explains the negligible gravitational stress. Notably, in the inner disk, Ṁ(λ) remains positive even when ℱRey is negative, likely due to sink particle effects and limited numerical resolution.

Overall, Fig. 12 indicates that midplane transport is dominated by turbulent stress, leading predominantly to outward spreading. For a typical disk with scale height ≈2.5 AU, sound speed ≈105 cm s−1, and column density ≈102 g cm−2, the measured median outward mass accretion rate of ∼10−5 M⊙ yr−1 corresponds to an effective Shakura & Sunyaev (1973) viscosity of αsh ≈ 0.09 (with α sh = M ˙ / 3 π c s h Σ Mathematical equation: $ \alpha_{\mathrm{sh}} = \dot{M}/3\pi c_{\mathrm{s}} h \Sigma $). We stress, however, that an α-disk prescription is not applicable in this context, as the turbulent transport exhibits strong temporal and spatial variability.

A key question remaining is the origin of this turbulent stress. To test whether accretion drives turbulence, we measured the amount of vertically landed mass in 1 kyr bins, which is long enough for several orbits in the outer disk and dozens in the inner disk. In Fig. 14a, the vertically accreted mass declines from 10−2 to 10−3 M⊙ in R1, while R2 shows a more variable but nevertheless decreasing trend. Figure 14b applies the same statistical method as Fig. 12, but now takes the maximum absolute median accretion rate across orbits in each temporal bin and plots it against the corresponding landed mass. The strong correlation confirms a link between vertical accretion and internal transport through turbulent stress.

Thumbnail: Fig. 14. Refer to the following caption and surrounding text. Fig. 14.

Correlation between vertical accretion and turbulent stress in runs R1 (blue) and R2 (red). Panel (a): Mass accreted vertically onto the disk, shown in 1 kyr bins. Panel (b): Maximum median accretion rate within each bin, plotted against the corresponding accreted mass.

6. Discussions

6.1. Asymmetric feeding of the disk

As in Commerçon et al. (2024), we have shown that the interplay between magnetic fields and turbulence produces anisotropic inflows in the form of streamers onto the disk. These streamers excite and sustain eccentric motions within the disk by continuously supplying material with an AMD. This provides a natural explanation for eccentric protoplanetary disks without requiring extended spiral waves or the presence of planets or stellar companions. Robust models like those presented here can be used to account for the complex kinematic features observed in disks (e.g., Zamponi et al. 2021; Thieme et al. 2023; van’t Hoff et al. 2023; Winter et al. 2025).

Perhaps the most novel result of this study is that asymmetric accretion also enables material to land vertically onto the disk, thereby driving turbulence within it. This turbulence leads to efficient angular-momentum transport inside the disk, which dominates over contributions from the Maxwell and gravitational stress. Thus, while magnetic fields govern angular-momentum transport outside the disk, the asymmetric feeding they facilitate becomes its primary driver within it. Most notably, this process yields highly efficient transport, equivalent to an effective viscosity parameter of α sh = M ˙ / 3 π c s h Σ ≈ 0.1 Mathematical equation: $ \alpha_{\mathrm{sh}} = \dot{M}/3\pi c_{\mathrm{s}} h \Sigma \approx 0.1 $, without the need to invoke magnetic instabilities such as the magneto-rotational instability (Balbus & Hawley 1991).

6.2. A link with cosmochemistry

The efficient angular momentum transport driven by turbulence indicates that disk material predominantly spreads outward, with only occasional episodes of inward motion (Fig. 12).

These results are relevant to cosmochemistry, where an isotopic dichotomy is observed in Solar System meteorites (Nanne et al. 2019; Morbidelli et al. 2024). Refractory inclusions (CAIs and amoeboid olivine aggregates), formed at high temperatures near the Sun, are found in the colder outer disk, a distribution generally interpreted as evidence of outward transport. In addition, carbonaceous chondrites (CCs) and non-carbonaceous chondrites (NCs) define two distinct isotopic reservoirs: NCs, depleted in supernova-derived nuclides, dominate the inner Solar System, whereas CCs are linked to the outer regions. Explaining this dichotomy requires models that couple infall history with disk dynamics.

One-dimensional viscous models (e.g., Pignatale et al. 2018; Morbidelli et al. 2022; Marschall & Morbidelli 2023) can reproduce outward transport of high-temperature condensates if the centrifugal radius remains small (< 1 AU) and if the effective viscosity parameter αsh is ∼0.1 throughout disk formation. However, such models process all material through high temperatures, which is inconsistent with the presence of volatiles (e.g., hydrogen bearing species) showing non-solar isotopic compositions in the outer disk.

Our simulations show that outward transport of solids can occur through disk spreading, reaching median decretion rates of ∼[10−5 − 10−4] M⊙ yr−1 thanks to Reynolds stresses (with an equivalent αsh ≈ 0.1). Because the sub-AU dynamics are unresolved, the spreading material does not attain CAI-forming temperatures. Redistribution by outflows is also possible (Bhandare et al. 2024), though not modeled here as we did not account for jets. Studies that resolve the collapse of the first Larson core find a rapidly expanding disk around the protostar (Vaytet et al. 2018; Machida & Basu 2019; Ahmad et al. 2024, 2025; Bhandare et al. 2024), whose radial spread may be sustained by the vertical accretion-driven turbulence reported here. Said vertical accretion initially delivers material across all radii during the first ∼60 kyr (Fig. 6), with a time-variable median landing radius (Fig. 15). Although transient surges in accretion toward the inner disk occur (Fig. 15), our models do not suggest a sustained inner-disk bias during the Class 0 phase. The vertical accretion patterns reported here may explain the presence of non-solar isotopic material in CC regions but not the predominantly NC composition of the inner disk. The isotopic dichotomy may therefore reflect (a) a time-dependent change in the isotopic composition of infalling material (Morbidelli et al. 2024), and (b) a different accretion mechanism with a sustained inner-disk bias in vertical accretion at later times13. Pressure maxima or proto-planets then preserve this distinction by stopping radial mixing.

Thumbnail: Fig. 15. Refer to the following caption and surrounding text. Fig. 15.

Temporal evolution of the fraction of mass landing at each radius, obtained by normalizing the incoming mass flux from Eq. (4) for each snapshot, for run R1 (a) and run R2 (b). The red curve indicates the median landing radius, where half of the mass lands at larger radii and half at smaller radii. The black curve indicates the location of the snow line. An animated movie, which combines panel (b) with a 3D volume render of the disk of run R2, is available at the CDS.

In summary, our results support a hybrid scenario in which the isotopic dichotomy of Solar System meteorites originates from time-dependent infall composition combined with the later emergence of an inner-disk bias in vertical accretion. This dichotomy is subsequently locked in place by disk structures such as pressure maxima or the birth of a proto-Jupiter.

6.3. Snow line

The location of the water snow line in disks is of interest for planetesimal formation (Stevenson & Lunine 1988; Armitage et al. 2016; Schoonenberg & Ormel 2017; Morbidelli et al. 2022). In 1D models, disks are typically assumed to undergo a gradual buildup over several 105 yr, with the water snow line drifting outward as the disk heats, followed by an inward contraction as accretion subsides (e.g., Pignatale et al. 2018; Drążkowska & Dullemond 2018; Lichtenberg et al. 2021). However, the chaotic accretion seen in our models suggests that the snow line likely evolves in a far less smooth and predictable manner.

We located the water snow line by calculating the water saturation pressure along the disk midplane using the Clausius–Clapeyron relation:

P sat ( T ) = P sat , 0 e − T a / T , Mathematical equation: $$ \begin{aligned} P_{\mathrm{sat} }(T) = P_{\mathrm{sat} ,0}\,e^{-T_{\mathrm{a} }/T}, \end{aligned} $$(19)

where Psat, 0 = 1.14 × 1013 g cm−1 s−2 and Ta = 6062 K (Lichtenegger & Komle 1991). The snow line is defined where the water partial pressure, XH2OP, equals the saturation pressure. Assuming an H2O abundance of XH2O = 5 × 10−4 (Harwit et al. 1998), we identified this location in the disk midplane and plot its evolution over time in Fig. 15 for both runs (black lines). The average midplane density and temperature are shown in Appendix E.

The two runs show distinct snow-line behaviors. In run R1, the gradual decrease in disk mass and slow outward spreading of its radius (Figs. 2a, b) cause the snow line to drift inward from ≈9 to ≈6 AU as the disk cools. In contrast, run R2 exhibits rapid variations in disk mass and radius (Figs. 2a, b), leading to a highly dynamic snow line that follows the disk’s mass and radius evolution. Accretion-driven increases in disk density and temperature push the snow line outward, reaching ≈13 AU, which reflects the hotter disk conditions in R2.

6.4. A dust-enriched disk

Having accounted for dust dynamics during the simulation, we could measure the metallicity of the disk and its temporal evolution. We display in Fig. 16 the dust-to-gas ratio of individual dust species within the disk as a function of time. The total dust-to-gas ratio ϵ (dash-dotted line) begins at 1% in both runs, and shows a similar increase over time, with ϵ reaching ≈1.2% by the final snapshot of run R1. However, despite the enrichment, the dust size is too small for appreciable grain growth (Krapp et al. 2019; Lovascio et al. 2022).

Thumbnail: Fig. 16. Refer to the following caption and surrounding text. Fig. 16.

Disk metallicity of runs R1 (blue) and R2 (red) as a function of time, with t = 0 marking the epoch of sink formation. The dash-dotted line indicates the total dust–to–gas ratio of the disk (ϵ), while the colored solid lines show the ratios of individual dust species, with darker shades corresponding to larger grain sizes.

Once in the disk, all dust species, by virtue of their small size and the high gas density (and hence small Stokes numbers), are tightly coupled to the gas and perfectly follow its kinematics. A more detailed analysis of dust kinematics in the disk would require that we describe its size evolution during the collapse through coagulation and fragmentation (Lombart et al. 2024), as well as using a more robust method to treat the dust dynamics in the poorly coupled regime (Verrier et al. 2025). Our results show that in the absence of grain growth, the gas dynamics dominate the distribution of dust and set the baseline metallicity of the disk.

6.5. Limitations

Our results must be considered in light of their limitations, whether numerical or due to missing physical processes. While the qualitative picture from our simulations is likely robust, the quantitative details should be interpreted with caution.

The resolution, though sufficient for the collapsing cloud and outer disk, cannot capture the inner disk with adequate fidelity, as too few cells resolve the innermost orbits. The sink particle replacing the sub-AU region further affects the global disk properties, most notably the mass (Machida et al. 2014; Vorobyov et al. 2019; Hennebelle et al. 2020b). We also did not account for high-velocity jets, which can remove substantial amounts of angular momentum (Machida & Basu 2019), and the Hall effect, a nonideal MHD process that can either enhance or suppress magnetic braking (Krasnopolsky et al. 2011; Tsukamoto et al. 2015; Wurster et al. 2016; Marchand et al. 2019). Furthermore, studies investigating the evolution of the grain distribution throughout the collapse shed doubt on the standard Mathis-Rumple-Nordsiek (Mathis et al. 1977) size distribution (Kawasaki & Machida 2023; Vallucci-Goy et al. 2024; Lombart et al. 2024).

Taken together, these caveats suggest that while our simulations provide valuable insights into disk evolution, they should be viewed as a step toward a more complete physical picture rather than a definitive model.

7. Conclusion

In this study, we modeled the self-consistent formation and evolution of protoplanetary disks from the gravitational collapse of 1 M⊙ and 3 M⊙ cloud cores, using radiative MHD simulations that include nonideal magnetic effects (ambipolar diffusion). We followed the evolution of the disks for ≈58.5 kyr and ≈27 kyr post-sink formation for both simulations. Our analysis shows how anisotropic accretion governs both the global properties and internal kinematics of nascent disks, and we accounted for eccentric fluid motions to avoid biasing our measurements. Our results can be summarized as follows:

  • (i)

    Magnetic fields, in tandem with turbulence during the collapse, drive asymmetric accretion via mechanisms such as the interchange instability. This process channels material onto the disk primarily through its vertical surfaces, yet the highest mass flux per unit area is from radial infall. This creates a complex accretion geometry that imprints distinct thermodynamic histories on the accreted gas. Gas tracer particles show that gas accreted through polar regions undergoes repeated heating and cooling cycles, while equatorial material follows a more classical adiabatic heating path prior to disk accretion.

  • (ii)

    Nascent disks exhibit significant and persistent eccentricities, e ∼ 0.1. This eccentricity is not transient but continuously generated and sustained by the anisotropic accretion. This supplies material with a net AMD, which sustains eccentricity despite numerical and physical damping.

  • (iii)

    Vertical accretion drives vigorous turbulent motion throughout the disk, which dominates angular momentum transport. The associated Reynolds stress (with an “effective” α ≈ 0.1) far exceeds the contribution from gravitational and Maxwell stresses and causes a rapid spreading of the disk with median mass decretion rates of ∼[10−5 − 10−4] M⊙ yr−1.

  • (iv)

    We place our results in the broader context of the observed isotopic dichotomy of Solar System meteorites: The vertical infall distributes material across a broad range of radii during the Class 0 phase, and the Reynolds stress it generates enables the outward transport of high-temperature condensates (CAIs and AOAs) without requiring all material to pass through the innermost regions of the disk. Although this can account for the presence of material with non-solar isotopic composition in regions inhabited by CC chondrites, the absence of a sustained inner-disk bias at the early times investigated here suggests that an alternative accretion regime is needed to explain the predominance of NC-like material in the inner disk.

Our results establish a direct link between magnetically driven accretion during collapse and the kinematic properties of protoplanetary disks. The resulting eccentricity and efficient turbulence driving have strong implications for the early stages of planet formation, as they set the orbital conditions for planetesimal growth, as well as for our interpretation of the Solar System’s cosmochemical record.

Data availability

Movie associated to Fig. 15 is available at https://www.aanda.org

Acknowledgments

We thank Matthias González for insightful comments during the writing of this paper. Simulation R1 was carried out by FL, and R2 by BC. This work has received funding from the French Agence Nationale de la Recherche (ANR) through the projects DISKBUILD (ANR-20-CE49-0006), and PROMETHEE (ANR-22-CE31-0020). We gratefully acknowledge the support from the CBPsmn (PSMN, Pôle Scientifique de Modélisation Numérique, Quemener & Corvellec 2013) at ENS de Lyon for providing computing resources to carry-out and analyze our simulations. EL was supported by ERC No. 864965 (PODCAST) under the EU’s Horizon 2020 programme. A.M. acknowledges the support from ERC grant 101019380 (HolyEarth).

References

  1. Ahmad, A., González, M., Hennebelle, P., & Commerçon, B. 2024, A&A, 687, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Ahmad, A., González, M., Hennebelle, P., Lebreuilly, U., & Commerçon, B. 2025, A&A, 696, A238 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Armitage, P. J., Eisner, J. A., & Simon, J. B. 2016, ApJ, 828, L2 [Google Scholar]
  5. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
  6. Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650 [Google Scholar]
  7. Bergin, E. A., & Williams, J. P. 2017, Astrophys. Space Sci. Lib., 445, 1 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bhandare, A., Commerçon, B., Laibe, G., et al. 2024, A&A, 687, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015 [Google Scholar]
  10. Bradski, G. 2000, Dr. Dobb’s J. Softw. Tools, 120, 122 [Google Scholar]
  11. Cacciapuoti, L., Macias, E., Gupta, A., et al. 2024, A&A, 682, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Calcino, J., Price, D. J., Hilder, T., et al. 2025, MNRAS, 537, 2695 [Google Scholar]
  13. Charnoz, S., Aléon, J., Chaussidon, M., Sossi, P., & Marrocchi, Y. 2024, LPI Contrib., 86, 6065 [Google Scholar]
  14. Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Commerçon, B., Lebreuilly, U., Price, D. J., et al. 2023, A&A, 671, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Commerçon, B., Lovascio, F., Lynch, E., & Ragusa, E. 2024, A&A, 689, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Drążkowska, J., & Dullemond, C. P. 2018, A&A, 614, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Fitzgibbon, A., Pilu, M., & Fisher, R. 1996, IEEE Trans. Pattern Anal. Mach. Intell., 21, 21 [Google Scholar]
  20. Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Goodchild, S., & Ogilvie, G. 2006, MNRAS, 368, 1123 [NASA ADS] [CrossRef] [Google Scholar]
  22. Greaves, J. S., & Rice, W. K. M. 2010, MNRAS, 407, 1981 [CrossRef] [Google Scholar]
  23. Harwit, M., Neufeld, D. A., Melnick, G. J., & Kaufman, M. J. 1998, ApJ, 497, L105 [Google Scholar]
  24. He, C.-C., & Ricotti, M. 2023, MNRAS, 522, 5374 [Google Scholar]
  25. He, C.-C., & Ricotti, M. 2025, MNRAS, 540, 175 [Google Scholar]
  26. Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016a, ApJ, 830, L8 [Google Scholar]
  27. Hennebelle, P., Lesur, G., & Fromang, S. 2016b, A&A, 590, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Hennebelle, P., Lesur, G., & Fromang, S. 2017, A&A, 599, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020a, A&A, 635, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Hennebelle, P., Commerçon, B., Lee, Y.-N., & Chabrier, G. 2020b, ApJ, 904, 194 [NASA ADS] [CrossRef] [Google Scholar]
  31. Kawasaki, Y., & Machida, M. N. 2023, MNRAS, 522, 3679 [NASA ADS] [CrossRef] [Google Scholar]
  32. Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17 [NASA ADS] [CrossRef] [Google Scholar]
  33. Könyves, V., André, P., Men’shchikov, A., et al. 2015, A&A, 584, A91 [Google Scholar]
  34. Krapp, L., Benítez-Llambay, P., Gressel, O., & Pessah, M. E. 2019, ApJ, 878, L30 [Google Scholar]
  35. Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2011, ApJ, 733, 54 [Google Scholar]
  36. Kuffmeier, M. 2024, Front. Astron. Space Sci., 11, 1403075 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kuznetsova, A., Bae, J., Hartmann, L., & Mac Low, M.-M. 2022, ApJ, 928, 92 [NASA ADS] [CrossRef] [Google Scholar]
  38. Lebreuilly, U., Commerçon, B., & Laibe, G. 2019, A&A, 626, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2024, A&A, 682, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Lee, Y.-N., Charnoz, S., & Hennebelle, P. 2021, A&A, 648, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Lee, Y.-N., Ray, B., Marchand, P., & Hennebelle, P. 2024, ApJ, 961, L28 [NASA ADS] [CrossRef] [Google Scholar]
  42. Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Lesur, G., Flock, M., Ercolano, B., et al. 2023, ASP Conf. Ser., 534, 465 [NASA ADS] [Google Scholar]
  44. Lichtenberg, T., Drążkowska, J., Schönbächler, M., Golabek, G. J., & Hands, T. O. 2021, Science, 371, 365 [NASA ADS] [CrossRef] [Google Scholar]
  45. Lichtenegger, H. I. M., & Komle, N. I. 1991, Icarus, 90, 319 [NASA ADS] [CrossRef] [Google Scholar]
  46. Lombart, M., Bréhier, C.-E., Hutchison, M., & Lee, Y.-N. 2024, MNRAS, 533, 4410 [CrossRef] [Google Scholar]
  47. Lovascio, F., Paardekooper, S.-J., & McNally, C. 2022, MNRAS, 516, 1635 [Google Scholar]
  48. Machida, M. N., & Basu, S. 2019, ApJ, 876, 149 [CrossRef] [Google Scholar]
  49. Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2014, MNRAS, 438, 2278 [NASA ADS] [CrossRef] [Google Scholar]
  50. Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Marchand, P., Tomida, K., Commerçon, B., & Chabrier, G. 2019, A&A, 631, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Marschall, R., & Morbidelli, A. 2023, A&A, 677, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24 [Google Scholar]
  54. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  55. Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Mayer, A. C., Zier, O., Naab, T., et al. 2025a, MNRAS, 537, 379 [Google Scholar]
  57. Mayer, A. C., Naab, T., Caselli, P., et al. 2025b, MNRAS, 543, 3321 [Google Scholar]
  58. Mignon-Risse, R., González, M., & Commerçon, B. 2021a, A&A, 656, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2021b, A&A, 652, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Morbidelli, A., Baillié, K., Batygin, K., et al. 2022, Nat. Astron., 6, 72 [Google Scholar]
  61. Morbidelli, A., Marrocchi, Y., Ali Ahmad, A., et al. 2024, A&A, 691, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130 [NASA ADS] [CrossRef] [Google Scholar]
  63. Mulders, G. D., Pascucci, I., Apai, D., & Ciesla, F. J. 2018, AJ, 156, 24 [Google Scholar]
  64. Najita, J. R., & Kenyon, S. J. 2014, MNRAS, 445, 3315 [NASA ADS] [CrossRef] [Google Scholar]
  65. Nanne, J. A. M., Nimmo, F., Cuzzi, J. N., & Kleine, T. 2019, Earth Planet. Sci. Lett., 511, 44 [CrossRef] [Google Scholar]
  66. Ogilvie, G. I. 2001, MNRAS, 325, 231 [NASA ADS] [CrossRef] [Google Scholar]
  67. Ogilvie, G. I., & Barker, A. J. 2014, MNRAS, 445, 2621 [NASA ADS] [CrossRef] [Google Scholar]
  68. Padoan, P., Pan, L., Pelkonen, V.-M., Haugbølle, T., & Nordlund, A. 2025, Nat. Astron., 9, 9 [Google Scholar]
  69. Pignatale, F. C., Charnoz, S., Chaussidon, M., & Jacquet, E. 2018, ApJ, 867, L23 [Google Scholar]
  70. Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158 [NASA ADS] [CrossRef] [Google Scholar]
  71. Quemener, E., & Corvellec, M. 2013, Linux J., 2013, 2013 [Google Scholar]
  72. Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380 [Google Scholar]
  73. Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188 [Google Scholar]
  74. Savvidou, S., & Bitsch, B. 2025, A&A, 693, A302 [NASA ADS] [Google Scholar]
  75. Schoonenberg, D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  77. Spruit, H. C., Stehle, R., & Papaloizou, J. C. B. 1995, MNRAS, 275, 1223 [NASA ADS] [CrossRef] [Google Scholar]
  78. Stevenson, D. J., & Lunine, J. I. 1988, Icarus, 75, 146 [NASA ADS] [CrossRef] [Google Scholar]
  79. Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
  80. Teyssier, R., & Commerçon, B. 2019, Front. Astron. Space Sci., 6, 6 [NASA ADS] [CrossRef] [Google Scholar]
  81. Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44 [NASA ADS] [CrossRef] [Google Scholar]
  82. Thieme, T. J., Lai, S.-P., Ohashi, N., et al. 2023, ApJ, 958, 60 [NASA ADS] [CrossRef] [Google Scholar]
  83. Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130 [CrossRef] [Google Scholar]
  84. Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
  85. Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, MNRAS, 452, 278 [NASA ADS] [CrossRef] [Google Scholar]
  86. Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, ASP Conf. Ser., 534, 317 [NASA ADS] [Google Scholar]
  87. Tu, Y., Li, Z.-Y., Lam, K. H., Tomida, K., & Hsu, C.-Y. 2024, MNRAS, 527, 10131 [Google Scholar]
  88. Tung, N.-D., Testi, L., Lebreuilly, U., et al. 2024, A&A, 684, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Vallucci-Goy, V., Lebreuilly, U., & Hennebelle, P. 2024, A&A, 690, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. van’t Hoff, M. L. R., Tobin, J. J., Li, Z. Y., et al. 2023, ApJ, 951, 10 [CrossRef] [Google Scholar]
  91. Vaytet, N., Commerçon, B., Masson, J., González, M., & Chabrier, G. 2018, A&A, 615, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Verrier, G., Lebreuilly, U., & Hennebelle, P. 2025, A&A, 701, A174 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Weiss, B. P., Bai, X. N., & Fu, R. R. 2021, Sci. Adv., 7, eaba5967 [NASA ADS] [CrossRef] [Google Scholar]
  95. Williams, J. P. 2012, Meteor. Planet. Sci., 47, 1915 [NASA ADS] [CrossRef] [Google Scholar]
  96. Winter, A. J., Benisty, M., Izquierdo, A. F., et al. 2025, ApJ, 990, L10 [Google Scholar]
  97. Wurster, J., & Lewis, B. T. 2020a, MNRAS, 495, 3795 [CrossRef] [Google Scholar]
  98. Wurster, J., & Lewis, B. T. 2020b, MNRAS, 495, 3807 [CrossRef] [Google Scholar]
  99. Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037 [Google Scholar]
  100. Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 481, 2450 [NASA ADS] [CrossRef] [Google Scholar]
  101. Xu, W., & Kunz, M. W. 2021, MNRAS, 508, 2142 [NASA ADS] [CrossRef] [Google Scholar]
  102. Yang, T. Q., & Federrath, C. 2025, MNRAS, 541, 1969 [Google Scholar]
  103. Zamponi, J., Maureira, M. J., Zhao, B., et al. 2021, MNRAS, 508, 2583 [NASA ADS] [CrossRef] [Google Scholar]
1

Mass estimates remain uncertain (see Bergin & Williams 2017; Tung et al. 2024; Savvidou & Bitsch 2025).

2

This instability generates a turbulent-like velocity field, which, while different from models that include an inherited Kolmogorov cascade from larger scales (see Kuffmeier 2024), leads to significant anisotropic accretion.

3

Done using OpenCV’s (Bradski 2000) implementation of Fitzgibbon et al. (1996)’s algorithm.

4

These are different from the late Class 2 streamers discussed in the literature (e.g., Pineda et al. 2020; Cacciapuoti et al. 2024), which originate from larger scales than that of the dense core.

5

See Appendix A for the disk definition.

6

Analyzed in Sect. 5.2.

7

See Appendix C for details.

8

We note an apparent asymmetry between top and bottom inflow, which is likely circumstantial given the identical initial conditions but differing cloud masses.

9

This method also drives them toward regions of converging flows.

10

See Sect. 5.2.

11

See Appendix B of Ogilvie & Barker (2014) for the mathematical expressions of Rλ and Rϕ.

12

They are contra-variant components.

13

The ambient magnetic field outside the disk will weaken as the envelope dissipates, reducing the strength of magnetic streams feeding the disk and potentially favoring the inner-most regions of the disk (Lee et al. 2021).

Appendix A: Disk definition

Following Commerçon et al. (2024), we defined the disk using a local angular momentum criterion: all cells with angular momentum exceeding 5 × 10−3 |L|max (where |L|max is the maximum local angular momentum in the simulation domain) are identified as disk material. An alternative criterion from Lee et al. (2021) selects cells above the minimum disk density. We compare these two definitions (respectively C1 and C2) in Fig. A.1, which shows that while disk mass and radius agree well, C2 yields consistently higher apparent eccentricities as it better captures the outermost orbit. For simplicity, we adopted C1.

Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Comparison of criteria C1 (blue) and C2 (orange) in run R2 (see text). Shown here are disk mass (a), 90% mass radius (b), and apparent eccentricity (c) as a function of time.

Appendix B: Counter-rotating streamer

In Fig. 2b, the disk radius in run R2 (red curve) shows a dip around t ≈ 15 kyr, caused by accretion from a counter-rotating streamer. Figure B.1 displays the z-component of the angular momentum,

J = r × m v Mathematical equation: $$ \begin{aligned} \boldsymbol{J} = \boldsymbol{r} \times m \boldsymbol{v} \end{aligned} $$(B.1)

(where m is the mass within a cell), at t ∼ 13 kyr, highlighting the streamer with opposite angular momentum. The influx of such material forces the disk to contract.

Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Slice in the z-direction of the z-component of the angular momentum in run R2 at t ∼ 13 kyr, showing a counter-rotating streamer delivering material to the disk. The white curves are velocity vector field streamlines.

Appendix C: Edge-on view of accretion

We have shown in Sects. 4.1 and 4.2 the origin of the accreted material that reaches the disk. This appendix further illustrates those results by presenting velocity vector field streamlines.

In Fig. C.1 we show the azimuthally averaged velocity vector field streamlines at the final simulation snapshots for runs R1 (top) and R2 (bottom). These highlight the outflow, the infall beneath the outflow, and the magnetically regulated infall. The extents of the spherical and cylindrical surfaces used in Sects. 4.1 and 4.2 are indicated by dash-dotted and dashed green lines, respectively.

Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

Azimuthally averaged velocity vector field streamlines of runs R1 (top) and R2 (bottom) at our final simulation snapshots. The dashed and dash-dotted green lines in each plot display the extents of the fixed surfaces used to quantify accretion in Sects. 4.1 and 4.2.

Appendix D: Gravitational stability of the disk

We analyzed the gravitational stability of our disks using the classical Toomre parameter (Toomre 1964), which measures the ratio of thermal and centrifugal support to self-gravity. A value of Q > 1 indicates that the disk is stable against axisymmetric gravitational perturbations.

Figure D.1 shows the radially averaged Toomre profiles for both runs. At all times in both simulations, we find Q ≫ 1, indicating strong stability throughout the simulation’s duration, particularly at later times. The more massive disk in Run 2 (R2) exhibits systematically lower values of Q than the disk in R1, yet remains well above the instability threshold.

Thumbnail: Fig. D.1. Refer to the following caption and surrounding text. Fig. D.1.

Averaged values of the real part of Toomre Q for runs R1 (a) and R2 (b) as a function of radius. Each curve corresponds to a different time, with t = 0 denoting the moment of sink formation. Only cells belonging to the disk were used when computing these curves.

Appendix E: Density and temperature profiles

The thermal history of the gas in the midplane has consequences on the mineralogy of the first condensates (Charnoz et al. 2024). Figure E.1 shows the mass-weighted density and temperature profiles for the disk midplanes in runs R1 and R2. For R1, both density and temperature decrease gradually over time. In contrast, R2 exhibits a more dynamic and chaotic history: episodes of disk accretion cause temporary increases in density and temperature, which are followed by cooling during quieter periods. Averages are computed using only cells belonging to the disk.

Thumbnail: Fig. E.1. Refer to the following caption and surrounding text. Fig. E.1.

Mass-weighted averages of the disk density (a and c) and temperature (b and d) for runs R1 (a and b) and R2 (c and d). Each curve denotes a different time, with t = 0 corresponding to the moment of sink formation. The vertical dotted gray lines correspond to the accretion radius of the sink particle.

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Column density images in the z (panels a, b, d, and e) and x (c and f) directions for runs R1 (first row) and R2 (second row), at both large scales (a and d) and small scales (b, c, e, and f). The lime colored contours in panels (b) and (e) are a least-square fit of an ellipse enveloping the disk’s surface. The associated eccentricity is shown in the top-right corner of panels (b) and (e).
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Global evolution of the protoplanetary disk as a function of time, where t = 0 represents the epoch of sink (i.e., protostellar) formation. The quantities shown are the protostellar mass (solid lines in panel a), disk mass (dotted lines in panel a), disk radius (solid lines in panel b), disk semimajor axis (dotted lines in panel b), and apparent disk eccentricity (c). The disk semimajor axis and apparent eccentricity are inferred from an elliptical fit. The blue (red) curve corresponds to the 1 M⊙ (3 M⊙) run R1 (R2). The shaded regions in panel (d) represent temporal fluctuations in the measurement, and the solid line is an average value.
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Hammer projections displaying the time and surface integrated radial mass flux on a sphere of Rshell = 30 AU for runs R1 (left) and R2 (right), displaying both the inflow (top colorbar) and outflow (bottom colorbar) of material throughout the simulation’s duration post-sink formation (tf), computed by integrating Eq. (1) in time. The lime colored contour delimitates the transition from positive to negative values.
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Total inflow (solid lines) and outflow (dotted lines) of mass (panel a) and specific angular momentum (panel b) through a sphere of radius Rshell = 30 AU around the sink particle as a function of time, where t = 0 corresponds to the epoch of sink formation, for runs R1 (blue) and R2 (red). This figure is complementary to Fig. 3.
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Total inflow of material through a cylinder of radius Rshell = 30 AU and height h = 16 AU throughout the simulation’s duration post-sink formation for runs R1 (≈58 kyr, first row) and R2 (≈27 kyr, second row), computed by integrating Eq. (3) (first column) and Eq. (4) (second and third columns) in time. The red circles in the second and third columns correspond to the minimum and maximum semimajor axis of the disk throughout each simulation’s duration, and thus represent the locations where material has landed in the disk.
In the text
Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Vertical mass accretion onto the disk for runs R1 (blue) and R2 (red), obtained by integrating Eq. (4) over ϕ and time. Dash-dotted and dotted lines show accretion through the top (z = 8 AU) and bottom (z = −8 AU) surfaces, respectively; the solid lines are their sum. Solid vertical lines mark the disk’s minimum and maximum radii during the simulation. The dashed lines represent the cumulative fraction of landed mass (normalized to unity), plotted against the right-hand axis. This figure is complementary to Fig. 5.
In the text
Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Comparison of the total (panel a) and surface-specific (panel b) mass accreted through a cylinder of radius Rcyl = 30 AU (solid lines) and the top and bottom disks of the cylinder at, respectively, z = 8 AU and z = −8 AU (dash-dotted lines). These were computed by integrating Eqs. (3) and (4) in space and time.
In the text
Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

Trajectories of the 5645 tracer particles that settle into the disk during run R2 in the x − y (top) and x − z (bottom) planes. The color scale reflects the cumulative distance traveled by each particle, with darker shades marking the most recent segments of each path.
In the text
Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

Temperature–density histories of randomly selected tracer particles that settle onto the disk from the polar regions (panel a) and from near the equatorial plane (panel b). Each curve corresponds to an individual particle, represented by a different color. The gray curve in panel (b) represents T ∝ ρ5/3 − 1.
In the text
Thumbnail: Fig. 10. Refer to the following caption and surrounding text. Fig. 10.

Distribution of eccentricity, computed from the disk kinematics (Eqs. 5 and 6), for run R1 (a) and run R2 (b). Each curve represents a different time, with t = 0 marking the epoch of sink formation. The dashed line indicates the largest radius enclosing 90% of the disk mass over the entire simulation. The sink accretion radius is indicated by the dotted line.
In the text
Thumbnail: Fig. 11. Refer to the following caption and surrounding text. Fig. 11.

Specific AMD within the disk as a function of time since sink formation for runs R1 (blue) and R2 (red).
In the text
Thumbnail: Fig. 12. Refer to the following caption and surrounding text. Fig. 12.

Internal kinematics of the disk. Shown here are the orbital mass accretion rate (a and d), turbulent stress tensor induced accretion rate (b and e), and Maxwell stress tensor induced accretion rate (c and f) for runs R1 (first row) and R2 (second row), computed as described in Sect. 5.2. These are computed in fixed logarithmic bins in λ throughout all simulation snapshots using only cells belonging to the disk, and we display the resulting first (blue), second (red), and third (black) quantiles, which respectively represent the 25th, 50th, and 75th percentiles. The gray curves in each plot represent minimal and maximal values.
In the text
Thumbnail: Fig. 13. Refer to the following caption and surrounding text. Fig. 13.

Same as Fig. 12 but this time showing the gravitational stress tensor-induced accretion rate (ℱ𝒢).
In the text
Thumbnail: Fig. 14. Refer to the following caption and surrounding text. Fig. 14.

Correlation between vertical accretion and turbulent stress in runs R1 (blue) and R2 (red). Panel (a): Mass accreted vertically onto the disk, shown in 1 kyr bins. Panel (b): Maximum median accretion rate within each bin, plotted against the corresponding accreted mass.
In the text
Thumbnail: Fig. 15. Refer to the following caption and surrounding text. Fig. 15.

Temporal evolution of the fraction of mass landing at each radius, obtained by normalizing the incoming mass flux from Eq. (4) for each snapshot, for run R1 (a) and run R2 (b). The red curve indicates the median landing radius, where half of the mass lands at larger radii and half at smaller radii. The black curve indicates the location of the snow line. An animated movie, which combines panel (b) with a 3D volume render of the disk of run R2, is available at the CDS.
In the text
Thumbnail: Fig. 16. Refer to the following caption and surrounding text. Fig. 16.

Disk metallicity of runs R1 (blue) and R2 (red) as a function of time, with t = 0 marking the epoch of sink formation. The dash-dotted line indicates the total dust–to–gas ratio of the disk (ϵ), while the colored solid lines show the ratios of individual dust species, with darker shades corresponding to larger grain sizes.
In the text
Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Comparison of criteria C1 (blue) and C2 (orange) in run R2 (see text). Shown here are disk mass (a), 90% mass radius (b), and apparent eccentricity (c) as a function of time.
In the text
Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Slice in the z-direction of the z-component of the angular momentum in run R2 at t ∼ 13 kyr, showing a counter-rotating streamer delivering material to the disk. The white curves are velocity vector field streamlines.
In the text
Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

Azimuthally averaged velocity vector field streamlines of runs R1 (top) and R2 (bottom) at our final simulation snapshots. The dashed and dash-dotted green lines in each plot display the extents of the fixed surfaces used to quantify accretion in Sects. 4.1 and 4.2.
In the text
Thumbnail: Fig. D.1. Refer to the following caption and surrounding text. Fig. D.1.

Averaged values of the real part of Toomre Q for runs R1 (a) and R2 (b) as a function of radius. Each curve corresponds to a different time, with t = 0 denoting the moment of sink formation. Only cells belonging to the disk were used when computing these curves.
In the text
Thumbnail: Fig. E.1. Refer to the following caption and surrounding text. Fig. E.1.

Mass-weighted averages of the disk density (a and c) and temperature (b and d) for runs R1 (a and b) and R2 (c and d). Each curve denotes a different time, with t = 0 corresponding to the moment of sink formation. The vertical dotted gray lines correspond to the accretion radius of the sink particle.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.