# Formation of low-mass protostars and their circumstellar disks

A. Ahmad, M. González, P. Hennebelle, B. Commerçon

► To cite this version:

A. Ahmad, M. González, P. Hennebelle, B. Commerçon. Formation of low-mass protostars and their circumstellar disks. Astronomy & Astrophysics - A&A, 2024, 687, pp.A90. 10.1051/0004-6361/202449340. cea-04629268

HAL Id: cea-04629268

https://cea.hal.science/cea-04629268v1

Submitted on 28 Jun 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Formation of low-mass protostars and their circumstellar disks*

A. Ahmad<fnref n="1" />, M. González<fnref n="1" />, P. Hennebelle<fnref n="2" />, and B. Commerçon<fnref n="3" />

Received 25 January 2024 / Accepted 21 April 2024

## ABSTRACT

*Context.* Understanding circumstellar disks is of prime importance in astrophysics; however, their birth process remains poorly constrained due to observational and numerical challenges. Recent numerical works have shown that the small-scale physics, often wrapped into a sub-grid model, play a crucial role in disk formation and evolution. This calls for a combined approach in which both the protostar and circumstellar disk are studied in concert.

*Aims.* We aim to elucidate the small-scale physics and constrain sub-grid parameters commonly chosen in the literature by resolving the star-disk interaction.

*Methods.* We carried out a set of very high resolution 3D radiative-hydrodynamics simulations that self-consistently describe the collapse of a turbulent, dense molecular cloud core to stellar densities. We studied the birth of the protostar, the circumstellar disk, and its early evolution (<6 yr after protostellar formation).

*Results.* Following the second gravitational collapse, the nascent protostar quickly reaches breakup velocity and sheds its surface material, thus forming a hot (~$10^3$ K), dense, and highly flared circumstellar disk. The protostar is embedded within the disk such that material can flow without crossing any shock fronts. The circumstellar disk mass quickly exceeds that of the protostar, and its kinematics are dominated by self-gravity. Accretion onto the disk is highly anisotropic, and accretion onto the protostar mainly occurs through material that slides on the disk surface. The polar mass flux is negligible in comparison. The radiative behavior also displays a strong anisotropy, as the polar accretion shock was shown to be supercritical, whereas its equatorial counterpart is subcritical. We also find a remarkable convergence of our results with respect to initial conditions.

*Conclusions.* These results reveal the structure and kinematics in the smallest spatial scales relevant to protostellar and circumstellar disk evolution. They can be used to describe accretion onto regions commonly described by sub-grid models in simulations studying larger-scale physics.

**Key words.** stars: early-type – stars: evolution – stars: formation – stars: low-mass – stars: pre-main sequence – stars: protostars

## 1. Introduction

Circumstellar disks form as a result of the conservation of angular momentum during the collapse of gravitationally unstable pre-stellar cloud cores. Understanding the formation of these disks and their subsequent evolution is of fundamental importance in astrophysics, as they are the birthplace of planets. This task is, however, heavily impeded by numerous challenges in observing star-forming regions, as the dense infalling envelope obscures the nascent disk during the class 0 phase. As such, most observational constraints come from more evolved class I and II disks. From a theoretical standpoint, the prohibitive time-stepping constraints in numerical simulations has made it nearly impossible to self-consistently describe the evolution of a newly formed circumstellar disk over a sufficiently large timescale in order to compare it with observations. To circumvent these constraints, theorists have abandoned the description of the innermost regions (<1 AU) and instead use sink particles (Bate et al. 1995; Bleuler & Teyssier 2014) onto which sub-grid physics are encoded. These particles interact with the surrounding gas through self-gravity and accretion as well as through radiative and mechanical feedback effects such as out-

flows and stellar winds. The parameters of these particles, such as their effective radius and accretion thresholds, have largely been chosen on the grounds of educated guesswork. Although necessary to study the global evolution of the disk, reducing the inner regions (which contain the protostar) to a sub-grid model can produce nonphysical results, especially when much of the actual sub-grid physics that have been encoded remain poorly constrained.

In this respect, Machida et al. (2014) investigated the effects of sink parameters on the formation of circumstellar disks. They found that the choice of the sink radius and its accretion threshold can, in conjunction with the physical model employed, dictate the formation and evolution of a circumstellar disk. Vorobyov et al. (2019) led a similar study but focused on the mass transport rate from the sink cell to the protostellar surface. They found that simulations with a slower mass transport rate would form more massive disks and that the accretion rate onto the protostar displayed more episodic behavior.

Finally, Hennebelle et al. (2020a) studied the influence of the sink accretion threshold on the global evolution of the disk. They found that while the mass contained within the sink is insensitive to this parameter, the disk radius and mass exhibit a strong sensitivity to it. Indeed, they found that the disk mass increases significantly at higher accretion thresholds.

<footnote n="1">¹ Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: adnan.ali.ahmad1998@gmail.com</footnote>

<footnote n="2">² Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Univ. Lyon, Ens de Lyon, Univ. Lyon 1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69007 Lyon, France</footnote>

<footnote>* Movies are available at https://www.aanda.org</footnote>
From these empirical studies, it has become clear that a
deeper understanding of the disk’s inner boundary is a nec-
essary endeavor to pursue in order to understand the global
disk evolution. Of course, this is not the only field of applica-
tion of sink particles. Much larger scale simulations that seek
to provide a protoplanetary disk or stellar population synthe-
sis by modeling the collapse of an entire molecular cloud,
such as those of Bate (2012, 2018), Hennebelle et al. (2020b),
Lebreuilly et al. (2021, 2024a,b), Grudić et al. (2022), also rely
on sink particles for their inclusion of small-scale physics. This
also applies to simulations studying the high-mass regime (e.g.,
Krumholz et al. 2009; Kuiper et al. 2010; Mignon-Risse et al.
2021a,b, 2023; Commerçon et al. 2022; Oliva & Kuiper 2023).
However, understanding the inner boundaries of circumstellar
disks requires a self-consistent description of the inner regions
in which one must model the formation of the protostar fol-
lowing a second gravitational collapse triggered by the disso-
ciation of H₂ molecules and its subsequent interaction with
the disk. Although Larson (1969) had pioneered second col-
lapse calculations in spherical symmetry, the field has since
developed ever more robust codes to tackle the problem in
the three dimensions necessary to describe the formation of a
circumstellar disk while including ever more complex physics
such as radiative transfer (e.g., Whitehouse & Bate 2006; Bate
2010, 2011; Ahmad et al. 2023) and magnetic fields under the
ideal and non-ideal approximations (Machida et al. 2006, 2007,
2008, 2011; Machida & Matsumoto 2011; Tomida et al. 2013,
2015; Bate et al. 2014; Tsukamoto et al. 2015; Wurster et al.
2018, 2021, 2022; Vaytet et al. 2018; Machida & Basu 2019;
Wurster & Lewis 2020). Although the latest studies struggle to
integrate across large timescales due to stringent time-stepping
constraints, an important result they’ve shown is that the higher
density gas ($\rho > 10^{-10} \text{ g cm}^{-3}$) is poorly magnetized due to
magnetic resistivities, thus placing the magnetic pressure orders
of magnitude below the thermal pressure. As such, in addi-
tion to greatly alleviating numerical constraints, omitting mag-
netic fields allows one to describe the inner sub-AU region
with reasonably high fidelity prior to the birth of a stellar
magnetic field through a dynamo process. Additionally, these
studies have not studied in depth the interaction between the
nascent protostar and its surrounding disk, mostly due to the
prohibitive time-stepping. Nevertheless, they have offered valu-
able insight regarding the system’s structure following the sec-
ond collapse phase. Indeed, they seem to indicate that the inner
regions are characterized by a density plateau in the innermost
region ($<10^{-2} \text{ AU}$), which then transitions toward a power-law
distribution (Saigo et al. 2008; Machida & Matsumoto 2011;
Tsukamoto et al. 2015; Vaytet et al. 2018).

In this paper, we investigate the inner boundaries of newly
formed circumstellar disks using high resolution 3D radiation-
hydrodynamics (RHD) calculations of the collapse of a dense
molecular cloud core to protostellar densities under the gray flux
limited diffusion approximation (FLD). To this end, we modeled
the collapse of the pre-stellar core, the formation of a first hydro-
static Larson core, the second gravitational collapse triggered by
the dissociation of H₂ molecules, and the subsequent early evo-
lution of the inner regions. Particular attention was given to the
interaction between the nascent protostar and the disk and how
such a process evolves over time.

## 2. Model

We carried out our simulations using the RAMSES (Teyssier
2002) adaptive mesh refinement (AMR) code and the same setup

employed in Ahmad et al. (2023) but with a notable difference
being the presence of angular momentum in the system through
the inclusion of an initially turbulent velocity vector field param-
eterized by a turbulent mach number $M_a$. Commerçon et al.
(2011, 2014) and González et al. (2015) implemented FLD in
the code. We used the equation of state of Saumon et al. (1995)
and the opacity tables of Semenov et al. (2003), Ferguson et al.
(2005), and Badnell et al. (2005), which were pieced together
by Vaytet et al. (2013). The initial conditions consist of a uni-
form density sphere of mass $M_0 = 1 M_\odot$, initial temperature
$T_0 = 10 \text{ K}$, and a radius of $R_0 = 2.465 \times 10^3 \text{ AU}$, thus yield-
ing a ratio of thermal to gravitational energy of $\alpha = 0.25$. We
present four runs in the main body of this paper where the ini-
tial amount of turbulence varies: $M_a = 0.2$ for run G1, 0.4 for
G2, 0.8 for G3, and 1 for G4. Although $M_a$ varies, the turbulent
seed does not. This means that run G4 has five times as much
initial angular momentum as G1. For comparative purposes, we
also ran two additional simulations, labeled G5 and G6. The
simulation of G5 possesses the same parameters as run G2 but
also includes solid body rotation in the initial cloud core, whose
ratio of rotational to gravitational energy is $\beta_{rot} = 10^{-2}$. Run G6
also contains the same parameters as G2, but it has a higher $\alpha$
of 0.5.

We used the same refinement strategy as in Ahmad et al. (2023); however, since angular momentum is present in the system, the simulations do not require a resolution as stringent as their spherically symmetrical counterpart, as the properties of the protostar, such as its central density and radius, are easier to resolve. Thus, we lowered the maximum refinement level $l_{max}$ to 25 in order to alleviate our time-stepping constraints, thus yielding a spatial resolution of $\Delta x = 2.9 \times 10^{-4} \text{ AU}$ at the finest level (the coarse grid is 64³ cells, with $l_{min} = 6$). Nevertheless, the stringent refinement criterion yielded some of the best-resolved disks in the literature. Indeed, the circumstellar disks in our simulations have $[95-2.7 \times 10^3]$ cells per Jeans length, with $\sim 10^7$ cells within their volume.

## 3. Results

Our simulations cover the initial isothermal contraction of the
cloud, the birth of the first Larson core, the second gravitational
collapse, and the subsequent evolution of the star-disk system.
The physics at large scales (i.e., from the cloud core to the first
Larson core) have been thoroughly discussed in the literature,
and as such, they are only briefly covered in our paper. Below, we
focus our study on the behavior of the system in the innermost
regions that contain the protostar.

### 3.1. *The dynamical range*

Herein, we illustrate the full dynamical range covered by all of
our simulations using run G4 as an example. To this end, we
display in Fig. 1 column densities at various scales (panels a–
d), and slices through the center of the domain displaying den-
sity (panels e–f), radiative flux (panels g–h), and temperature
(panels i–j). Panel a displays the column density at the scale
of the dense molecular cloud core. Here, a filamentary struc-
ture of size $\sim 10^3 \text{ AU}$ formed by gravo-turbulence can be seen
(Tsukamoto & Machida 2013), and within this structure a first
Larson core is born. In this run, the first core lifetime is sig-
nificantly extended thanks to ample angular momentum, which
reduces the mass accretion rate onto it. As such, a disk was able
<figure><img src="image_2.png" /><figcaption>Fig. 1. Visualization of the entire dynamical range covered in our simulations. The figure displays data taken from the final snapshot of run G4. Panels a–d show the projected gas column density along the y direction across multiple scales, ranging from ≈9.35 × 10³ AU in panel a to 3 AU in panel d. Panels e–j are slices through the center of the domain along the y direction for panels e, g, and i, and along the z direction for panels f, h, and j. These display density (panels e and f), radiative flux (panels g and h), and temperature (panels i and j). The color bars in panels e and f are centered on the density of the inner disk’s shock front (≈1.5 × 10⁻⁹ g cm⁻³ at this snapshot).</figcaption></figure>

to form around it<fnref n="1" />, as seen in panel c. Within this disk, the second collapse takes place and gives birth to the protostar and the circumstellar disk, as seen in the lower panels.

## 3.2. Protostellar breakup

We begin by describing the structure of the system shortly after the onset of the second gravitational collapse using Fig. 2. We

use data from run G6 as an example, although the evolutionary sequence displayed here applies to all other runs as well. The top row of this figure (panels a–d) displays the local radiative flux, an excellent tracer of shock fronts. It is computed as

$$F_{\text{rad}} = - \frac{c\lambda\nabla E_r}{\rho\kappa_R}, \qquad (1)$$

where c is the speed of light, λ the Minerbo (1978) flux limiter, E_r the radiative energy, ρ the gas density, and κ_R the Rosse-land mean opacity. At a shock front, kinetic energy is converted into radiation and as such it is accompanied by an increase in

<footnote n="1">¹ In the case of runs G1 and G2, no disk was formed at these scales prior to the onset of the second collapse. This is discussed further in Sect. 3.5.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Demonstration of the breakup of the protostar due to excess angular momentum. The data is taken from run G6. Each column represents a different time, where t = 0 marks the birth of the protostar. Panels a–d are slices along the z direction through the center of the domain that display the local radiative flux emanating from the cells (see Eq. (1)). The scale bar in panel a applies to panels b–d as well. The mass of the protostar at each snapshot is displayed in the top-right corners of panels a–d. Panels e–h display 2D histograms binning all the cells in our computational domain, which show the distribution of azimuthal velocities divided by $\sqrt{|g_r|}r$ with respect to radius. The color code in the histograms represents the co-latitude $\theta$ divided by $\pi$, where $\theta/\pi = 0.5$ corresponds to the equator and $\theta/\pi = 1$ (respectively 0) corresponds to the south (respectively north) pole. The dotted black lines in panels e–h display $v_\phi/v_{crit} = 1$ (see Eq. (2)). An animated version of this plot is available online.</figcaption></figure>

radiative flux. Hence, this quantity prominently displays accretion shocks and spiral waves.

The bottom row of the figure (panels e–h) displays the azimuthal velocity<fnref n="2" /> distribution of all cells in our computational domain with respect to radius, which we have divided by the critical velocity beyond which the centrifugal force exceeds the radial component of the gravitational force:

$$v_{\text{crit}} = \sqrt{|g_r|}r, \qquad (2)$$

where

$$g_r = -\frac{\partial\phi}{\partial r}. \qquad (3)$$

Here, $\phi$ is the gravitational potential obtained through the Poisson equation.

Once the gas has completed its dissociation of H₂ molecules and ample thermal pressure support is gathered, the second Larson core (i.e., the protostar) is formed. At birth (first column), the protostar is a thermally supported spherical object, and its azimuthal velocities are well below $v_{crit}$. A mere eight and a half days later (second column), the protostar has nearly doubled in mass, and an equatorial bulge is now visible in panel b. This is due to the fact that as the protostar accretes, it is also accumulating angular momentum. Nevertheless, it is still rotating below

$v_{crit}$. A month later (third column), the outer shells of the protostar finally exceed $v_{crit}$, after which material spreads outward and transitions to a differential rotation profile in which the centrifugal force is now the main counterbalance to gravity. The sustained accretion ensures a constant flow of material within the protostar exceeding $v_{crit}$.

### 3.3. An embedded protostar

The expulsion of material by the protostellar surface will naturally lead to the formation of a circumstellar disk. Here, we study how such a disk grows and evolves while analyzing the accretion mechanism onto the protostar and the disk. To this end, we studied Fig. 3, which displays density slices through the center of the computational domain for run G1 at different times. The velocity streamlines are color coded with the local radial mass-flux $-\rho v_r$, where red colors denote inward transport of material, blue colors denote outward transport, and white signifies very little transport. The evolutionary sequence displayed here applies to all other runs.

Panels a and f display the system once temperatures exceed 2000 K and the dissociation of H₂ is triggered, where the gas spirals inward almost isotropically. The central region accumulates material quickly, and once ample thermal pressure support is gathered, the protostar forms. We display in panels b and g the structure of the system once the protostar reaches breakup velocity. As the protostar’s surface begins expelling

<footnote n="2">² $v_\phi$ was computed along the angular momentum vector of the gas within 0.5 AU.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Top-down (top row, panels a–e) and edge-on (bottom row, panels f–j) slices through the center of the domain of run G1 displaying density and velocity streamlines. The color coding in the velocity streamlines displays the local radial mass flux $-\rho v_r$. Each column displays a different epoch, where $t = 0$ (panels b and g) corresponds to the moment of protostellar breakup. The scale bar in panel a applies to all other panels. An animated version of this plot is available online.</figcaption></figure>

material, a disk immediately forms afterward, and as time progresses, an increasing amount of material collides with the disk instead of the protostar, thus causing the former to grow significantly over time. Hereafter, we refer to this newly formed circumstellar disk as the inner disk (in accordance with the terminology of Machida & Matsumoto 2011). In panel d, we observed the development of spiral waves, which seem to have subsided into near-circular waves in panel e<fnref n="3" />. We note that during this phase, the accretion timescale of the disk ($M_d/\dot{M}_d \sim 10^{-2} M_\odot/10^{-2} M_\odot \text{yr}^{-1} = 1 \text{ yr}$) is shorter than its dynamical timescale ($2\pi R_d/v_\phi \sim 2\pi \times 1 \text{ AU}/3 \text{ km s}^{-1} \approx 10 \text{ yr}$). This means that any angular momentum redistribution process within the inner disk occurs on longer timescales than accretion. Thus, accretion is the dominant process behind the expansion of the disk.

We now turn to studying the accretion process with the aid of the streamlines in Figs. 3 and 4, which displays the radial mass flux in slices through the center of the domain. In addition, we display unbroken velocity vector field streamlines in Fig. 5 at a curated moment. Although the polar regions initially bring a large amount of material to the central protostar, the polar reservoir of gas is quickly depleted and by $t \approx 268$ days (fourth columns and onward of Figs. 3 and 4), very little mass is accreted through the poles. Indeed, most of the material landing at the protostellar surface is sliding on the inner disk’s surface, as its velocity component normal to the disk surface is not strong enough to break through the shock front. We note, however, that some material landing on the disk surface can sometimes break through the shock front and is then transported into the inner disk, as can be seen in panel b of Fig. 5. The previously men-

tioned spiral waves can be seen transporting material radially in panels d and e. This is more apparent in Fig. 6, which displays the radial mass flux averaged in radial bins and measured separately for both the upper layers of the disk and its main body for all our runs. One can see that in all runs and across all radii, the upper layers of the disk have a strictly positive radial mass flux, whereas the main body shows alternating inflows and outflows of material.

What is most visually striking from the edge-on views is the vertical extent of the disk: it is substantially flared, giving it the shape of a torus. This is more so apparent in the 3D rendering of the system displayed in Fig. 7. The inner disk’s surface (rendered in blue) completely engulfs the protostar (rendered in green). Figure 7 also displays the 3D velocity streamlines, which show that the material accreted through the poles carries with it angular momentum as the gas is spiraling inward. The cross-sectional slices in the figure display the radiative flux, which reveals shock fronts and spiral waves. Interestingly, there does not seem to be any shock fronts separating the protostar from the inner disk (barring the spiral waves). This means that the accretion shock envelopes both the protostar and the inner disk, and the two act as a continuous fluid system. As such, differentiating the protostar and the inner disk is rather difficult, but the rotational profiles seem to indicate that the protostar is in solid body rotation and the inner disk exhibits differential rotation<fnref n="4" />. All other runs have displayed an identical structure of the inner disk.

<footnote n="3">³ The gravitational stability of the inner disk is discussed in Sect. 3.6.</footnote>

<footnote n="4">⁴ See Appendix A for an overview of how each object was defined.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Same as Fig. 3 but displaying only the local radial mass flux $-\rho v_r$ in order to better demonstrate how accretion occurs in the star and inner disk system. An animated version of this plot is available [online](#).</figcaption></figure>

### 3.4. A convergent structure

We display in Figs. 8 and 9 the radiative flux in slices through the center of the domain for runs G1–G6 at a moment in time where the inner disk has reached a radius of ≈0.5 AU. We note the almost identical structure of the inner disk in all runs: it is toroidal and highly flared. We also note that run G3 (panels c and q) displays a strong eccentricity, as the outer disk that formed around the first Larson core in this run was already highly eccentric prior to the second collapse.

In the top-down slices (panels a–d), we notice ripples in radiative flux. These are spiral waves which have relaxed into nearly circular perturbations. We note that runs G1 and G2 have stronger spirals waves (i.e., the radiative flux emanating from them is stronger) due to their higher mass. (For a more in-depth analysis of these spiral waves, see Sect. 3.6).

An interesting observation from panels e–h of Fig. 8 is the prominence of the radiative flux along the poles. Indeed, the polar region is much less dense than the equator, which causes the radiative flux to escape much more easily along this direction. This also causes the gas to heat up in the polar direction<fnref n="5" />.

### 3.5. Evolution

Having ascertained the structure of the system in the innermost regions following the second gravitational collapse, we follow the temporal evolution of the inner disk with the aid of Fig. 10. The properties of the inner disk once it has reached a radius of 1 AU are of particular interest, as that is the most commonly chosen sink radius.

First, we point out the different evolutionary history of each simulation in panel a. This figure displays the maximum density in our computational domain as a function of time since the birth of the first core. The moment in which each curve exhibits a sharp rise in central density corresponds to the onset of the sec-

ond gravitational collapse. Here, we observed clearly that simulations with higher initial amounts of angular momentum have a delayed onset of second collapse, as the additional centrifugal support significantly extends the first core’s lifetime by reducing its mass accretion rate. In runs G3, G4, G5, and G6, the first core lifetime is long enough for it to have a disk built around it such that the inner disk forms within a disk (Machida & Matsumoto 2011). Despite the differing evolutionary histories, the resulting properties of the inner disk, displayed in panels b, c, e, and f, show remarkable similarity. Indeed the temporal evolution of the inner disk equatorial radius $R_d$ (panel f), shows very little spread. The specific angular momentum of the inner disks, displayed in panel e, also exhibit striking similarity. This shows that the amount of specific angular momentum in the inner disk is independent of the initial amount of angular momentum in the parent cloud core, a result that is in agreement with Wurster & Lewis (2020)'s non ideal MHD and hydro simulations. Furthermore, the entirety of the angular momentum budget of the first core is found within the inner disk and protostar after it is accreted.

The curves in panel c display the inner disk’s mass ($M_d$), which exhibit the same evolutionary trend and $M_d$(1 AU) ∈ [1.634 × 10⁻²; 2.755 × 10⁻²] M☉. The mass of the protostar (panel d), although very similar in runs G1–G4, has runs G5–G6 as outliers since $M_*$ is about ≳40% larger in these runs. Interestingly, $M_*$ also seems to be decreasing in most runs, meaning that the protostar is shedding its mass to the disk due to excess angular momentum. The notable exception is run G1, in which the protostar’s mass is increasing due to strong gravitational torques.

Finally, we turn our attention to panel b, which displays the density of the inner disk’s equatorial shock front ($\rho_s$). We measure this quantity along the equator since that is the region where most of the incoming mass flux lands on the inner disk (as shown in Fig. 3). As such, this quantity is an equivalent to the accretion threshold used in sink particles ($n_{acc}$). We report $\rho_s$(1 AU) ∈ [5.35 × 10⁻¹⁰; 2 × 10⁻⁹] g cm⁻³. The most commonly adopted accretion threshold in the literature is 1.66 × 10⁻¹¹ g cm⁻³ ($n_{acc}$ = 10¹³ cm⁻³). We thus suggest higher values of $n_{acc}$ when possible

<footnote n="5">⁵ The radiative behavior of the system is studied in Sect. 5.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Unbroken velocity vector field streamlines of run G1 at $t \approx 322$ days after the birth of the inner disk illustrating the dynamics of accretion onto the disk and protostar in a top-down (panel a) and edge-on (panel b) view. The scale bar in panel b applies to panel a as well. An animated version of these streamlines, made using windmap (https://github.com/rougier/windmap) is available online.</figcaption></figure>

for studies employing sink particles whose radius is 1 AU. However, we acknowledge that this can significantly increase the numerical cost of simulations, and thus might be too constraining for certain studies, particularly those that aim for very long temporal evolution.

The reason for such a convergence in the inner disk properties is the first core itself. As discussed in Ahmad et al. (2023), this hydrostatic object halts any inward accretion until temperatures can exceed the H₂ dissociation temperature of 2000 K, after which the second collapse ensues. As such, the mass accretion rate asymptotically reaches $c_s^3/G \sim 10^{-2} M_\odot \text{yr}^{-1}$, Larson 1969; Penston 1969) independently of initial conditions, provided that a first core forms. The small spread we observed in Fig. 10 is the result of our turbulent initial conditions, as no discernible trend can be inferred from their differences. We expect the large-scale initial conditions to play a more significant role later on when the entirety of the first Larson core is accreted and that mass accretion rates onto the innermost regions are dictated by transport processes within the disk and by the infall of material onto said disk.

We note that in the case where an outer disk exists prior to the onset of a second collapse, the inner disk will simply merge with it (Machida & Matsumoto 2011). This was the case for runs

<figure><img src="image_7.png" /><figcaption>Fig. 6. Average radial mass flux measured in radial bins for runs G1 (blue), G2 (orange), G3 (green), G4 (red), G5 (black), and G6 (magenta) at a moment in time where the inner disk has reached $\approx 1$ AU in radius. Only cells belonging to the inner disk were considered (see Appendix A for information on how the inner disk was defined). The solid lines represent measurements made on the upper layers of the disk, whereas the dashed lines represent measurements made in its main body.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Three-dimensional view of the inner disk and protostar at a moment in time when the former has reached $\approx 0.5$ AU in radius. The blue structure is the surface of the inner disk. The inner $r < 0.1$ AU region has been cut out in order to reveal the flow onto the protostar (rendered in green). The white curves are velocity vector field streamlines launched along the poles to reveal polar accretion. The bottom, left, and right panels are cross sections through the center of the domain displaying the radiative flux. The visualized volume is $1.6 \times 1.6 \times 1.6$ AU$^3$. An animated version of this plot is available online.</figcaption></figure>

G3, G4 and G5. The first core lifetime in runs G1 and G2 was not long enough for it to gather enough material around it to form a disk.

### 3.6. Gravitational stability of the inner disk

The existence of such a massive circumstellar disk naturally begs the question of whether it will undergo fragmentation or not. In order to determine the gravitational stability of the inner disk, we use the classical Toomre $Q$ parameter (Toomre 1964):

$$Q = \frac{\omega c_s}{\pi G \Sigma}, \qquad (4)$$
<figure><img src="image_9.png" /><figcaption>Fig. 8. Slices through the center of the domain for runs G1, G2, G3, and G4 (respectively first, second, third, and fourth columns) showing the local radiative flux in a top-down (top row, panels a–d) and edge-on (bottom row, panels e–h) view. The slices illustrate the structure of the accretion shock as well as the presence of spirals in the inner disk (see Eq. (1)). The scale bar in panels d and h apply to all other panels as well. The slices are shown at a moment in time where the inner disks have reached a radius of ≈0.5 AU. The mass of each inner disk is displayed in the top-right corners of panels a–d.</figcaption></figure>

where $c_s$ is the gas sound speed, $\Sigma$ its surface density, and $\omega$ its epicyclic frequency, defined as

$$ \omega^2 = \frac{1}{r^3} \frac{\partial(r^4 \Omega^2)}{\partial r}, \qquad (5) $$

where $\Omega$ is the angular velocity of the gas. This parameter represents the ratio of the outward pointing forces on the gas, namely the centrifugal and pressure gradient forces, to the inward pointing gravitational force. If $Q < 1$, the disk is unstable and a fragmentation is likely. We measure the sound speed and the angular velocity by averaging along the vertical axis of the inner disk.

We display the real part of $Q$ in top-down slices through the center of the domain for run G1 in Fig. 11 at curated moments. In panel a, we display the circumstellar disk at its birth, just after the protostar exceeded breakup velocity and began shedding its mass. This results in a ring of gas surrounding the protostar, whose $Q$ value is above unity. In panel b, we display the system just prior to the onset of the first large amplitude spiral wave. Here, the $Q$ parameter remains above unity in the innermost regions of the disk; however, an hourglass-shaped region has $Q < 1$ at slightly larger radii. The ratio of inner disk mass to protostellar mass has also increased by a factor ≈7. In panel c, a coherent two-armed spiral wave is launched from the center, and it grows radially as it is sheared apart by differential rotation when propagating outward. Finally, panels d and e show that these spiral waves relax into nearly circular ripples as a result of the increase in temperature. The $Q$ values in panels b–e hovers around unity throughout the disk, meaning that the disk remains marginally stable against gravitational instabilities despite its high mass relative to that of the protostar. This is due to its very high temperature.

A recent numerical study by Brucy & Hennebelle (2021) has shown that the fragmentation barrier of disks is quite blurry; it is better described by a probabilistic approach and said probabili-

<figure><img src="image_10.png" /><figcaption>Fig. 9. Same as Fig. 8 but for runs G5 and G6 (respectively first and second columns).</figcaption></figure>

ties strongly depend on how efficiently the disk is able to cool. As such, the fragmentation of the disk is also set by the cooling criterion. In our case, the inner disk is still strongly accreting, and a considerable amount of accretion energy is absorbed by the disk, particularly along the optically thick equator (this is shown in Sect. 5). This ensures that the inner disk remains hot and thus mostly stable against fragmentation despite its very high mass relative to the protostar. The primary regulator of the inner disk
<figure><img src="image_11.png" /><figcaption>Fig. 10. Temporal analysis of the inner disk in all our runs. Blue is for G1, orange for G2, green for G3, red for G4, black for G5, and magenta for G6. Panel a provides the context in which the inner disk is born by displaying the central density's evolution since the formation of the first Larson core (defined as the moment where $\rho_c > 10^{-10}$ g cm$^{-3}$). Panels b–e respectively display as a function of the inner disk's equatorial radius $R_d$ (analogous to time) the density of the inner disk's equatorial shock front, the mass of the inner disk, the protostellar mass, and the specific angular momentum of the inner disk. Panel f displays the evolution of the inner disk radius with respect to time, where $t = 0$ marks the moment of birth of the inner disk. The gray dotted line in panels b–e marks $R_d = 1$ AU. (For information on how the inner disk was defined, see Appendix A).</figcaption></figure>

temperature during this phase is the endothermic dissociation of H$_2$, which places $\gamma_{eff}$ at $\approx 1.1$.

## 4. Gas structure and kinematics

Herein, we provide a more quantitative analysis of the structure and kinematics of the inner disk with the aid of Fig. 12. Panels a–e of this figure display various quantities azimuthally averaged in radial bins in the equatorial region, where the equator was defined as the region in which $\theta/\pi \in [0.45; 0.55]$, where $\theta$ is the co-latitude. The curves are shown at a moment in time in which the inner disk has reached a radius of $\approx 1$ AU.

Panel a displays the equatorial density curve, which exhibits a plateau in the innermost regions ($r < 10^{-2}$ AU) and a power law tail$^6$. This behavior has been described by all previous 3D studies in the literature that employ either pure hydro or non-ideal MHD. This density structure, as well as the radial velocity curves displayed in panel b, suggest that no discontinuity in the flow separates the protostar from the inner disk. The only way for us to differentiate the two is by studying the rotational behavior of the system: we observed an object in solid body rotation in the inner-most regions ($j \propto r^2$, panel d) and a transition to a differential rotation profile (panel c).

### 4.1. Deviations from Keplerian rotation

We observed a significant deviation from Keplerian rotation ($v_\phi \propto r^{-0.5}$) in the inner disk. Indeed, it seems that $v_\phi \propto r^{-0.3}$.

In order to explain such a difference, one must start by analyzing the balance equation between the centrifugal, pressure, and gravitational forces (Pringle 1981; Lodato 2007):

$$ \frac{v_{\phi}^{2}}{r} = \frac{1}{\rho} \frac{\partial P}{\partial r} - g_{r}, \qquad (6) $$

where $P$ is the thermal pressure. By assuming a radial density profile where $\rho \propto r^{-\beta}$ and radial isothermality ($\partial c_s^2 / \partial r \approx 0$), we may write

$$ v_{\phi}^{2} \approx -g_{r}r - \beta c_{s}^{2}, \qquad (7) $$

where $c_s$ is the gas sound speed. We may approximate $g_r$ from the column density profile of the inner disk:

$$ g_r \approx -\frac{GM_{\text{enc}}(r)}{r^2}, \qquad (8) $$

where $G$ is the gravitational constant and

$$ M_{\text{enc}}(r) = M_{*} + 2\pi \int_{R_{*}}^{r} \Sigma(r')r'dr'. \qquad (9) $$

Here, $\Sigma$ is the disk's surface density. We note that this is merely an analytical estimate of $g_r$, which is different from the true potential computed in Eq. (3). Now let us assume a column density profile of $\Sigma \propto r^{-\xi}$. This means that the $g_r$ term in Eq. (7) scales with $r^{-\xi+1}$, whereas the $\beta c_s^2$ term scales with $r^{-\beta(\gamma-1)}$ ($T \propto \rho^{\gamma-1}$). From Fig. 12, we can write $\xi \approx 3/2$ and $\beta \approx 3$. Additionally, $\gamma \approx 1.1$ in the inner disk (panel e of Fig. 12). Thus, $g_r \propto r^{-0.5}$ and $\beta c_s^2 \propto r^{-0.3}$. This means that we may expect stronger deviations from Keplerian velocity at larger radii in the inner disk. At $R_d = 1$ AU, we have $M_d \sim 10^{-2} M_\odot$

<footnote>$^6$ We note that the disk in which the protostar was born in run G6 had fragmented prior to second collapse, and three first Larson cores exist within it. This causes the density spikes seen at larger radii.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Top-down slices through the center of the domain of run G1 displaying the real part of the Toomre Q parameter (see Eq. (4)) at different times. Here, $t = 0$ (panel a) corresponds to the moment of birth of the inner disk. Only cells belonging to the inner disk were used in the calculation of Q. The ratio of the inner disk mass to protostellar mass is displayed on the bottom left of each panel. The scale bar in panel a applies to all other panels. An animated version of this plot is available online.</figcaption></figure>

and $c_s \sim 1 \text{ km s}^{-1}$. As such, $\sqrt{1 + \beta c_s^2/g_t r} \approx 0.8$, and we thus expect the deviation from Keplerian rotation to be of the order of $(v_\phi - v_K)/v_K \approx -20\%$ (where $v_K = \sqrt{GM_{enc}/r}$). Panel g of Fig. 12 shows us that these are reasonable approximations.

Note that what we call “Keplerian” rotation in the present paper is different from its commonly adopted meaning in the literature. Indeed, the literature defines the Keplerian velocity as $v_{K,lit} = \sqrt{GM_{enc}/r}$; however, as the disk’s mass exceeds that of the protostar by up to a factor of seven, the measured deviations from $v_{K,lit}$ would be ~100%.

## 4.2. Analytical description of the inner disk structure

Herein, we set out to provide an analytical description of the structure of the inner disk, which is rather exotic by virtue of its very high mass relative to that of the protostar, as well as its very high temperature. This is not meant to be a full analytical development, but rather, one that seeks to deepen our understanding of the structure witnessed in our simulations. Namely, we seek to understand why $\Sigma \propto r^{-3/2}$, and in turn provide an analytical prediction of $n_{acc}$. To do so, we use as our starting point the important result shown in Sect. 3.6, namely that $Q \approx 1$ throughout most of the inner disk. This allowed us to link the sound speed $c_s$ and the disk’s column density profile $\Sigma$ through

$$c_s \approx \frac{\pi G \Sigma}{\Omega}. \qquad (10)$$

We can write the following power-law descriptions of these two quantities

$$\Sigma = \Sigma_0 \left(\frac{r}{R_*}\right)^{-\xi}, \qquad (11)$$

$$c_s = c_{s0} \left(\frac{r}{R_*}\right)^{-\beta(\gamma-1)/2}, \qquad (12)$$

where we have made use of the fact that $\rho \propto r^{-\beta}$. We note that $\Omega$ is constrained by $\Sigma$ through the equation

$$\Omega \approx \sqrt{\frac{G}{r^3}} \left( M_* + 2\pi \int_{R_*}^{r} \Sigma r' dr' \right). \qquad (13)$$

This allowed us to write

$$\Sigma_0 \approx \frac{c_{s0}}{\pi G} \sqrt{\frac{GM_*}{R_*^3}}. \qquad (14)$$

To describe $\Sigma$, one must thus first obtain the boundary values at $r = R_*$ and $r = 1$ AU. Since $\Sigma$ and $c_s$ are linked through Eq. (10), this is equivalent to finding the boundary values of $c_s$. Our simulations indicate $\beta \approx 3$, and $\gamma \approx 1.1$ in the inner disk. Hence, $c_s \propto r^{-0.15}$. This allowed us to write

$$\log \left( \frac{c_{s0}}{c_s} \right) = \frac{\beta(\gamma - 1)}{2} \log \left( \frac{r}{R_*} \right) = 0.15 \times \log \left( \frac{r}{R_*} \right). \qquad (15)$$

Since $R_* \sim 10^{-2}$ AU,

$$c_{s0} = 10^{0.3} \times c_s(1 \text{ AU}) \approx 2 \times c_s(1 \text{ AU}). \qquad (16)$$

This means that when traveling from the protostellar surface to $r = 1$ AU, $c_s$ reduces by a factor of approximately two. The temperature of the inner disk is $\sim 10^3$ K, hence $c_s(1 \text{ AU}) \sim 1 \text{ km s}^{-1}$ and $c_{s0} \sim 2 \text{ km s}^{-1}$. By using Eq. (10), we can link this boundary condition to $\xi$:

$$\left(\frac{c_s}{c_{s0}}\right)^2 = \left(\frac{r}{R_*}\right)^{3-2\xi} \frac{M_*}{M_{enc}}, \qquad (17)$$

$$\Rightarrow \left(\frac{c_s}{c_{s0}}\right)^2 = \left(\frac{r}{R_*}\right)^{3-2\xi} \frac{M_*}{M_* + 2\pi \frac{\Sigma_0}{R_*^{\xi}(-\xi+2)}(r^{-\xi+2} - R_*^{-\xi+2})}. \qquad (18)$$

Numerically solving Eq. (18) for $\xi$ such that $(c_s/c_{s0})^2 \approx 1/4$ yields $\xi \approx 1.38$ (with $M_* \sim 10^{-3} M_\odot$), which is close to the value witnessed in most of our runs ($\approx 3/2$, panel h of Fig. 12).

Once we had a description of both $\Sigma$ and $c_s$, we could describe the density profile with the aim of predicting $n_{acc}$ as

$$\rho = \frac{\Sigma}{2h}, \qquad (19)$$

where $h$ is the disk scale height:

$$h = \frac{c_s}{\Omega}. \qquad (20)$$

Thus, this analytical framework provides $\rho(1 \text{ AU}) \approx 5.37 \times 10^{-10} \text{ g cm}^{-3}$ ($n_{acc} \approx 3.23 \times 10^{14} \text{ cm}^{-3}$, when considering a mean molecular weight of one). This result is within the bounds provided by the simulations ($[5.35 \times 10^{-10}; 2 \times 10^{-9}] \text{ g cm}^{-3}$), and thus we consider it to be satisfactory.
<figure><img src="image_13.png" /><figcaption>Fig. 12. Studying the structure and kinematics of the gas in our simulations. Panels a–e display a set of azimuthally averaged quantities along the equatorial regions with respect to radius, where the equator was defined as the region where $\theta/\pi \in [0.45; 0.55]$ at a moment in time where the inner disk has reached ≈1 AU in radius. Panel a displays the density, panel b shows the radial velocity, panel c presents the azimuthal velocity, panel d shows the specific angular momentum, and panel e presents the temperature. Panel f displays the enclosed angular momentum profile. Panel g displays the deviation from Keplerian rotation along the disk midplane. Panel h displays the column density profile of the inner disk, where non-disk cells were masked.</figcaption></figure>

### 4.3. Radial transport within the inner disk

Once we had an analytical framework with which we could describe the inner disk, it was of interest to quantify the transport of material within it. More specifically, we wished to describe the transport processes using a simple alpha-disk model (Shakura & Sunyaev 1973) and compare it with the measured values within our simulations. A common approach in this regard is to measure the stress tensors induced by turbulent fluctuations and the self-gravity of the disk, respectively $\alpha_R$ and $\alpha_{grav}$ (e.g., Lodato & Rice 2004; Brucy & Hennebelle 2021; Lee et al. 2021). However, in our case, we found these measurements difficult to interpret, as the disk is awash with eccentric motions (Lovascio et al., in prep.), turbulent eddies, spiral waves, and erratically infalling material (see for example Fig. 5), which resulted in unreliable values of $\alpha_R$ and $\alpha_{grav}$. Instead, we opted for a simpler approach in which we attempted to fit an approximate analytical description of the transport within the inner disk with the measured values of our simulations. The mass accretion rate at any given radius in our simulations can be obtained using

$$ \dot{M}_s(r) = -2\pi r \Sigma \langle v_r \rangle, \quad (21) $$

where $\langle v_r \rangle$ is the vertically averaged radial velocity. We note that the majority of the inward transport of material occurs in the upper layers of the disk, as seen in Fig. 6. Along the midplane, material tends to spread outward, causing negative mass accretion rates. As such, when performing the vertical average of $v_r$, we weighed $v_r$ in two different ways:

$$ \langle v_r \rangle_z = \frac{1}{z_{\max} - z_{\min}} \int_{z_{\min}}^{z_{\max}} v_r dz, \quad (22) $$

$$ \langle v_r \rangle_{\rho} = \frac{1}{\Sigma} \int_{z_{\min}}^{z_{\max}} \rho v_r dz, \quad (23) $$

where $z_{\min}$ and $z_{\max}$ are respectively the minimal and maximal heights of the disk at a given radius and azimuth. $\langle v_r \rangle_z$ will thus be biased by the infall in the upper layers of the disk, whereas $\langle v_r \rangle_{\rho}$ will be biased by the dense midplane.

The analytical estimate of the mass accretion rate is (Shakura & Sunyaev 1973; Pringle 1981)

$$ \dot{M}(r) = 3\pi v \Sigma, \quad (24) $$

where $v$ is the effective viscosity of the inner disk

$$ v = \frac{\alpha_{\text{sh}} c_s^2}{\Omega}. \quad (25) $$

Here, $\alpha_{\text{sh}}$ is the Shakura & Sunyaev (1973) alpha which dictates the vigour of radial mass transport. Using the equations developed in the previous section, we may write Eq. (24) as

$$ \dot{M}(r) = 3\alpha_{\text{sh}} \frac{c_{s0}^3}{G} \left(\frac{r}{R_*}\right)^{-\beta(\gamma-1)-\xi+3/2} \sqrt{\frac{M_*}{M_{\text{enc}}}}, \quad (26) $$

$$ \Rightarrow \dot{M}(r) = 3\alpha_{\text{sh}} \frac{c_{s0}^3}{G} \left(\frac{r}{R_*}\right)^{-\beta(\gamma-1)-\xi+3/2} \times \sqrt{\frac{M_*}{M_* + 2\pi \frac{\Sigma_0}{R_*^\xi(-\xi+2)} (r^{-\xi+2} - R_*^{-\xi+2})}} \quad (27) $$

We may obtain an estimate of the $\alpha_{\text{sh}}$ parameter in Eq. (26) with the aid of Fig. 13, which displays the mass accretion rate in our simulations computed using Eq. (21) at a moment in time when the inner disk radius has reached ≈ 1 AU. We measure $\dot{M}_s$ using both $\langle v_r \rangle_z$ (panel a) and $\langle v_r \rangle_{\rho}$ (panel b). Using average stellar parameters reported in the figure caption, we over-plot Eq. (26) (solid gray lines). The figure shows that the vigour of inward transport varies across all radii, and in the case of runs G2-6, the
<figure><img src="image_14.png" /><figcaption>Fig. 13. Radial mass accretion rate within the inner disk for all runs (solid colored lines) as a function of cylindrical radius and computed using Eq. (21) at a moment in time when R<sub>d</sub> has reached ≈1 AU. Only cells belonging to the inner disk were used in the computation of M&#x0307;<sub>s</sub>. Panel a displays the measurements of M&#x0307;<sub>s</sub> made using 〈v<sub>r</sub>〉<sub>z</sub> (Eq. (22)), whereas panel b displays the same measurements made using 〈v<sub>r</sub>〉<sub>ρ</sub> (Eq. (23)). The solid gray lines are analytical estimates of radial transport computed using Eq. (26), with α<sub>sh</sub> = 10<sup>-1</sup> (top), 5 × 10<sup>-2</sup> (middle), and 10<sup>-2</sup> (bottom). The other parameters for the solid gray lines are M<sub>*</sub> ≈ 7 × 10<sup>-3</sup> M<sub>☉</sub>; R<sub>*</sub> ≈ 3 × 10<sup>-2</sup> AU; c<sub>s0</sub> ≈ 5 km s<sup>-1</sup>; β ≈ 3; ξ ≈ 3/2, and γ ≈ 1.1.</figcaption></figure>

outer layers ($r > 0.2$ AU) of the inner disk have negative mass accretion rates, meaning that material is mostly spreading outward. In panel b, we saw that the transport of material in the main body of the inner disk fluctuates wildly. Furthermore, any inward transport in this region consistently has lower mass accretion rates than in the upper layers of the disk (panel a). Indeed, the upper layers of the disk have a mass accretion rate that can be approximated by Eq. (26) with $\alpha_{sh} \sim 5 \times 10^{-2}$.

We note, however, that M&#x0307;<sub>s</sub> varies not only in space but also in time, and so the value of α<sub>sh</sub> reported here is not valid throughout the entirety of the class 0 phase. Indeed, once most of the remnants of the first Larson core are accreted, and that the mass accretion rate onto the star-disk system significantly reduces, this will inevitably cause a decrease in α<sub>sh</sub>, as M&#x0307;<sub>s</sub> is mainly dominated by the infall and will thus drop by several orders of magnitude. We would thus have a disk whose angular momentum transport is very weak, with α<sub>sh</sub> ~ [10<sup>-5</sup>; 10<sup>-3</sup>], but whose mass is significantly higher than that of the protostar. In addition, panel d of Fig. 10 shows that despite the inward transport of material within the inner disk, the protostellar mass is decreasing in most runs as a result of excess angular momentum. This allowed us to conclude that although the inner disk can be described by the physics of alpha disks, such a description is a first-order approximation whose results are not entirely reliable.

## 5. Radiative behavior

Although the protostar has not reached the temperatures required for fusion yet (>10<sup>6</sup> K), it will dominate the radiative output of

the system by virtue of its high temperature and its accretion luminosity. A quantitative analysis of said radiative behavior is seldom provided in the literature, and it is the purpose of this section.

To this end, we first begin by providing a qualitative overview of the radiative behavior of the system at medium (i.e., 10<sup>2</sup> AU) scales. As discussed previously, the structure of the accretion shock is rather complex, as it envelopes both the protostar and the inner disk. Nevertheless, the polar accretion shock, as seen in the bottom rows of Figs. 8 and 9, produces the majority of the radiative flux. In addition, the density cavity along the poles allows said radiation to escape much more easily than along the optically thick equator. This is reflected in Fig. 14, which displays the radiative flux emanating from each cell in edge-on slices across the center of our domain at curated moments for run G5. The protostar is born embedded within a disk that formed around the first Larson core (panel a). As time progresses, the radiation produced at the protostellar accretion shock escapes and brightens the polar regions significantly (panels b–d). However, the equatorial regions remain dark as the radiative flux struggles to pierce through the highly opaque disk. As such, the disk is almost unaffected by the protostar’s radiation.

We now turn to providing a quantitative analysis of the radiative behavior of the protostar and inner disk. To do so, as in Vaytet et al. (2018), Bhandare et al. (2020), and Ahmad et al. (2023), we compare two quantities: the radiative flux just upstream of the accretion shock ($F_{rad}$, see Eq. (1)), and the incoming accretion energy flux. However, since the shock front’s structure is complex, we measure these quantities in only four directions: north-south along the poles, and east-west along the equator (see panel a of Fig. 15):

$$F_{\text{acc,pol}} \approx -\rho v_r \frac{GM_{\text{enc}}(R_*)}{R_*} - E v_r - P v_r - \frac{\rho v_r^3}{2}, \quad (28)$$

$$F_{\text{acc,eq}} \approx -\rho v_r \frac{GM_{\text{enc}}(R_d)}{R_d} - E v_r - P v_r - \frac{\rho v_r^3}{2}, \quad (29)$$

where $R_*$ (resp. $R_d$) is the protostellar (resp. inner disk) radius, E the gas internal energy, P its thermal pressure, and

$$M_{\text{enc}}(r) = 4\pi \int_{0}^{r} \rho r'^{2} dr'. \quad (30)$$

This allowed us to define the radiative efficiency as

$$f_{\text{acc,pol}} \approx \frac{F_{\text{rad,pol}}}{F_{\text{acc,pol}}}, \quad (31)$$

$$f_{\text{acc,eq}} \approx \frac{F_{\text{rad,eq}}}{F_{\text{acc,eq}}}. \quad (32)$$

These are approximate measurements of the radiative efficiency of the accretion shock because the radiative flux also contains the cooling flux emanating from the protostellar interior, although the later remains smaller than that produced at the accretion shock because of the low temperature of the protostar (~10<sup>4</sup> K) prior to deuterium burning.

We display in panels b and c of Fig. 15 the resulting measurements of $f_{acc,pol}$ and $f_{acc,eq}$ obtained through ray-tracing. A schematic drawing displaying the location at which each quantity was measured is presented in panel a. In contrast to spherically symmetrical calculations, the radiative efficiency of the accretion shock displays a strong anisotropy: the polar accretion
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Edge-on slices through the center of the domain displaying the local radiative flux at different times (see Eq. (1)). Here, $t = 0$ (panel a) corresponds to the moment of protostellar birth. The data is taken from run G5. The scale bar in panel a applies to all other panels. An animated version of this plot is available [online](#).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Quantitative analysis of the radiative behavior of the accretion shock. Panel a is a schematic representation of the star and inner disk system, where the solid line is the location of the accretion shock that envelopes both the protostar and the disk. The dashed lines represent four rays launched from the center of the system along which the location of the accretion shock is found and its radiative efficiency measured just upstream from it (numbered red dots). The resulting measurements are presented in panels b and c. North and South (respectively solid and dashed lines in panel b and East and West (respectively solid and dashed lines in panel c) are shown for all simulations. These display the polar and equatorial radiative efficiencies (respectively $f_{acc,pol}$ and $f_{acc,eq}$; see Eqs. (31) and (32)).</figcaption></figure>

shock (panel b) is much more efficient than its equatorial counterpart (panel c). Indeed, the polar shock front reaches supercriticality ($f_{acc,pol} = 1$) in less than 2.5 years for most runs, whereas $f_{acc,eq} < 1$ throughout all simulations. This is due to the polar density cavity that allows the accretion shock to shine into an optically thin medium, whereas the equatorial accretion shock remains optically thick due to the remnants of the first Larson core and the presence of an extended outer disk around it.

The drop in $f_{acc,eq}$ seen in most runs (e.g., at $t \approx 1.2$ yr for run G1) is due to the equatorial shock front expanding out of the opacity gap, a region where temperatures are high enough to sublimate dust (see Fig. 1 of Ahmad et al. 2023). This causes the shock front to shine into a region of higher opacity, which in turn reduces its luminosity. Eccentric inner disks can also cause an east-west anisotropy, as is most apparent in run G3. In contrast,

we witness very minor anisotropies when comparing north-south radiative efficiencies.

The low radiative efficiency of the equatorial shock front has an important consequence on the structure of the inner disk: as it accretes, the majority of the accretion energy is dumped into its thermal budget, thus causing it to maintain a high temperature and to swell along the vertical extent.

# 6. Discussions

## 6.1. The star-disk structure

In the common star formation paradigm, the protostar is often seen as a standalone object wholly separate from the circumstellar disk due to an accretion shock separating the two. Our results
show that such a boundary does not exist: the transition from the protostar to the inner disk is a smooth one, and the two act as a continuous fluid system. This means that until either torque mechanisms transport sufficient material inward, or a protostellar magnetosphere arises, the protostar and the inner disk will continue to behave as a continuous fluid system.

Unfortunately, obtaining observational constraints on the structure of the system in such a young and deeply embedded object is very difficult, but two recent studies by Laos et al. (2021) and Le Gouellec et al. (2024) seem to be offering hints on the accretion mechanism during the class 0 phase. Whereas Laos et al. (2021) argues for magnetospheric accretion similarly to class I objects by basing their arguments on the similar shapes of the Bry line profiles, Le Gouellec et al. (2024) additionally analyzed the velocity profiles of these lines which exhibited clear differences. Based on this result, they concluded that the accretion mechanism in class 0 sources must be of a different nature. Additionally, they report very vigorous accretion unto the central regions. These observations are an indirect probe of the structure of the system following the second collapse. If a stellar magnetosphere cannot be inferred from them, then we argue that a structure similar to that described in this paper may be present.

## 6.2. Toward an eventual fragmentation

Although our analysis of the gravitational stability of the inner disk indicates that it is currently marginally stable against gravitational instabilities, it is also in a transient state in which it maintains an incredibly high temperature through accretion. This is by virtue of the low radiative efficiency of the inner disk's shock front. The protostar, as shown in Sect. 5, contributes almost nothing to heating the inner disk because it is currently embedded within it and the majority of its radiative flux escapes along the poles. One can imagine that such a disk will inevitably cool down once accretion subsides, and thus would be much more prone to fragmentation. The high mass budget of the inner disk also indicates that such fragmentations can occur multiple times, and thus lead to the formation of multiple star systems, or perhaps of Jupiter-like planets in close proximity to the central star.

## 6.3. Caveats

Although the radiative-hydro approximation is valid for the timescales described here, this obviously cannot be the case on longer timescales, as the stellar magnetic field will undoubtedly increase once a dynamo process begins. Furthermore, the hot, dense, and highly turbulent inner disk may also generate a dynamo process (Balbus & Hawley 1991; Wardle 1999; Lesur et al. 2014; Riols & Latter 2019; Deng et al. 2020). The onset of strong magnetic fields will bring about outflows and jets, as well as induce strong magnetic torques that can transport material inward and also reduce the rotation rate of the protostar to observed values in more evolved systems (~10% of breakup velocity, Hartmann & Stauffer 1989; Herbst et al. 2007). These will likely heavily influence the evolution of the system in the innermost regions. Nevertheless, the simulation presented in Machida & Basu (2019) was integrated for 2000 years after protostellar birth with non-ideal MHD, and their results seem to indicate that $\rho_s(1\,\text{AU}) \sim 10^{-9}\,\text{g cm}^{-3}$ (see their Fig. 13). Their protostar also seems to progressively separate from the inner disk (see panels m–p of their Fig. 5 and $v_r$ in their Fig. 11), likely

as a result of magnetic torques. Although state of the art non-ideal MHD simulations indicate that thermal pressure support far outweighs its magnetic counterpart in circumstellar disks (e.g., Machida et al. 2010; Vaytet et al. 2018; Machida & Basu 2019; Lee et al. 2021), magnetic torques likely outweigh turbulent viscosity or gravitational torques (as quantified by Machida & Basu 2019).

In addition, although we notice relatively little spread in our results, more exhaustive constraints on the inner boundaries of circumstellar disks can be obtained by exploring a more varied range of $M_0$, $\alpha$, and $\beta_{\text{rot}}$. This was not done due to computational constraints, as the simulations are expensive to run at such high resolution (each run used at least 100 000 CPU hours in total).

## 6.4. Consequences on global disk evolution scenarios

As a result of the high values of $n_{\text{acc}}$ hereby reported, the mass of circumstellar disks in numerical simulations may have been underestimated by about an order of magnitude (Hennebelle et al. 2020a). This appears to be inconsistent with observations of class 0 disks, which report lower disk masses of $\sim 10^{-2} M_{\oplus}$ (Tobin et al. 2020). We note, however, that a recent study by Tung et al. (2024) has found that class 0 disk masses are routinely underestimated with current observational techniques. Nevertheless, the results indicate that the disk is very massive at birth, and naturally one must question how such a disk might evolve over time.

As such, we put forward two speculative evolutionary scenarios that we believe to be possible. The first scenario rests on the previously discussed fragmentation of the disk; as accretion subsides and the disk cools, it will become more prone to gravitational instabilities. Should the inner disk fragment, a significant amount of angular momentum could be extracted from the system.

Another scenario would rest on the strength of magnetic fields in the inner disk. Perhaps the results of current state of the art papers that report high thermal to magnetic pressure ratios are overestimating the strength of magnetic resistivities in the gas whose density exceeds first Larson core densities ($\rho > 10^{-13}\,\text{g cm}^{-3}$), and thus would be underestimating the magnetic field strength within the protostar's and inner disk's precursor. This is in light of new studies (Lebreuilly et al. 2023; Kawasaki & Machida 2023; Tsukamoto et al. 2023; Vallucci-Goy et al., in prep.) that cast doubt on the validity of the MRN (Mathis-Rumple-Nordsiek, Mathis et al. 1977) dust size distribution used in said state of the art papers. In the case where the resistivities inside the first Larson core are overestimated, we would expect stronger coupling within the gas prior to hydrogen ionization. Although the ratio of thermal to magnetic pressure is still expected to greatly exceed unity within the inner disk, this would undoubtedly increase the strength of magnetic torques and cause a more central distribution of material, in which the protostar quickly exceeds the inner disk's mass and separates itself from it. However, simulations employing MRN resistivities create disk radii in broad agreement with observational surveys (Maury et al. 2019; Tobin et al. 2020), in which the magnetic field truncates the disk radius at $\sim 10^1\,\text{AU}$. Perhaps the disk radii can serve as a lower bound on magnetic resistivities for low density gas, whereas the strength of magnetic fields in young stellar objects ($\sim 10^3\,\text{G}$, Durney et al. 1993; Chabrier & Küker 2006) can serve as an upper bound on resistivities within the first Larson core.
## 6.5. Comparison with previous works

We now turn to comparing our results with previous calculations in the literature that have resolved the birth of the protostar and the inner disk structure described in this paper. Due to the very stringent time-stepping, only a handful of such studies exist. To our knowledge, the first report of the existence of a swollen disk-like structure around the second core is that of Bate (1998), who carried out their calculations using the SPH numerical method under the barotropic approximation. Saigo et al. (2008) later confirmed the existence of such a structure with their own barotropic calculations, this time using a nested-grid code. Although these papers do not have enough details presented for us to quantitatively compare our results, the description they provide of the structure of the second core and its surroundings seem to be qualitatively similar to ours. Following these two studies, Machida & Matsumoto (2011) led a detailed study of the inner disk under the barotropic approximation in which they include magnetic fields with ohmic dissipation, where they found qualitatively similar results to the hydro case. Notably, they report an inner disk mass of ~$10^{-2} M_\odot$, and a protostellar mass of ~$10^{-3} M_\odot$, in accordance with our results.

More recent studies include those of Wurster et al. (2018, 2022), Vaytet et al. (2018), Machida & Basu (2019), Wurster & Lewis (2020). The simulations presented in Vaytet et al. (2018) were also run with the RAMSES code while including magnetic fields. They reported the existence of a circumstellar disk around the protostar when including magnetic resistivities. They also reported that they could not resolve the shock front separating the protostar from the inner disk, but as we have demonstrated in this paper, such a shock front does not exist. The mass of the disk reported in their paper is ~$10^{-4} M_\odot$, which is to be expected given their short simulation time following protostellar birth (24 days) that results from stringent time-stepping.

As stated previously, Machida & Basu (2019) have also studied the inner disk, in which they were able to follow the calculations for a period of 2000 years following protostellar birth. Their results seem to indicate similar ratios of protostellar to disk mass than our paper, and the structure of their inner disk also seems to be similar to ours despite the presence of a magnetically launched outflow and a high velocity jet.

Finally, Wurster et al. (2018, 2021, 2022), Wurster & Lewis (2020) ran these calculations using SPH and all non-ideal MHD effects: ohmic dissipation, ambipolar diffusion, and the hall effect. They found the same disk structure following the second collapse phase. In contrast to Machida & Basu (2019), they report no high velocity jets and argue that magnetically launched outflows are circumstantial and depend on the initial turbulent velocity vector field, as well as the non-ideal effects at play. The only way in which we may quantitatively compare our results to theirs is through the protostellar mass, which seems to be ~$10^{-3} M_\odot$. The inner disk mass was not measured in these papers.

In summary, the circumstellar disk structure described in this paper is routinely found in previous papers simulating the second collapse in 3D. Although the quantitative details of the properties of said disk may differ due to different numerical methods and physical setups, it is guaranteed to be recovered in simulations that possess sufficient amounts of angular momentum in the first core. As such, the only simulations that do not recover it are those that make use of the ideal MHD approximation, which extracts too much angular momentum from the system and prevents the protostar from ever reaching breakup velocity.

## 7. Conclusion

We have carried out a set of high resolution 3D RHD simulations that self-consistently model the collapse of a $1 M_\odot$ dense molecular cloud core to stellar densities with the goal of studying the innermost (<1 AU) regions. Our results are summarized as follows:

(i) Following the second gravitational collapse, the protostar is formed through hydrostatic balance. Through accretion, the protostar accumulates angular momentum and reaches breakup velocity, after which it sheds some of its mass to form a hot, dense, and highly turbulent circumstellar disk, which we call the inner disk. The protostar is embedded within this disk, and no shock front separates the two. As accretion continues, the disk completely engulfs the protostar and spreads outward due to a combination of excess angular momentum and accretion. The disk mass exceeds that of the protostar by a factor of approximately seven, which means that the majority of the mass following the second collapse resides in the disk and its self-gravity dominates, with notable contributions from thermal pressure on its dynamics. In the case where an outer disk exists prior to the second collapse, this circumstellar disk forms within it, and the two merge after the inner disk spreads to sufficiently large radii.

(ii) Despite the differing evolutionary histories at larger spatial scales, the star-disk structure formed after the onset of the second collapse is identical, with a small spread caused by the turbulent initial conditions. This is due to the formation of the first Larson core in all our simulations, a hydrostatic object that ensures that the second collapse occurs in approximately similar conditions.

(iii) Accretion onto the protostar mainly occurs through material that slides on the disk's surface, as polar accretion has a low-mass flux in comparison. Along the equator, material spreads outward due to excess angular momentum. Accretion onto the inner disk is highly anisotropic.

(iv) The radiative emissions of the star-disk system are anisotropic. The radiative efficiency of the accretion shock is supercritical along the poles, whereas the inner disk's equatorial shock front is subcritical and on the order of $10^{-3}$.

(v) The density of the inner disk's shock front at 1 AU (the most commonly used sink radius) is in the range of $[5.35 \times 10^{-10}; 2 \times 10^{-9}] \text{ g cm}^{-3}$, which is about an order of magnitude higher than the commonly used sink accretion threshold of $1.66 \times 10^{-11} \text{ g cm}^{-3}$. Thus, we suggest higher accretion thresholds for studies employing sink particles whenever possible. We note, however, that our results correspond to very early times and may not be applicable throughout the entirety of the class 0 phase.

(vi) In order to physically decouple the protostar from its disk and reduce its rotation rate to observed values, torque mechanisms need to transport a sufficient amount of angular momentum outward. Magnetic fields are likely to play this role.

These results reveal the structure and kinematics of the innermost regions of circumstellar disks, which are often omitted from simulations due to computational constraints. Although our results are valid for the timescales described here (<6 years following protostellar birth), we expect magnetic fields to play a more significant role later on, particularly in creating powerful torques that transport material toward the protostar.
**Acknowledgements.** This work has received funding from the French Agence Nationale de la Recherche (ANR) through the projects COSMHIC (ANR-20-CE31-0009), DISKBUILD (ANR-20-CE49-0006), and PROMETHEE (ANR-22-CE31-0020). We have also received funding from the European Research Council synergy grant ECOGAL (Grant: 855130). We thank Ugo Lebreuilly, Anaëlle Maury, and Thierry Foglizzo for insightful discussions during the writing of this paper. The simulations were carried out on the Alfven supercomputing cluster of the Commissariat à l'Énergie Atomique et aux énergies alternatives (CEA). Post-processing and data visualization was done using the open source Osyris (https://github.com/osyris-project/osyris) package. 3D visualizations were done using the open-source PyVista (https://docs.pyvista.org/version/stable/) package (Sullivan & Kaszynski 2019).

## References

Ahmad, A., González, M., Hennebelle, P., & Commerçon, B. 2023, A&A, 680, A23

Badnell, N. R., Bautista, M. A., Butler, K., et al. 2005, MNRAS, 360, 458

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Bate, M. R. 1998, ApJ, 508, 508

Bate, M. R. 2010, MNRAS, 404, L79

Bate, M. R. 2011, MNRAS, 417, 2036

Bate, M. R. 2012, MNRAS, 419, 3115

Bate, M. R. 2018, MNRAS, 475, 5618

Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362

Bate, M. R., Tricco, T. S., & Price, D. J. 2014, MNRAS, 437, 77

Bhandare, A., Kuiper, R., Henning, T., et al. 2020, A&A, 638, A86

Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

Brucy, N., & Hennebelle, P. 2021, MNRAS, 503, 4192

Chabrier, G., & Küker, M. 2006, A&A, 446, 1027

Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35

Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11

Commerçon, B., González, M., Mignon-Risse, R., Hennebelle, P., & Vaytet, N. 2022, A&A, 658, A52

Deng, H., Mayer, L., & Latter, H. 2020, ApJ, 891, 154

Durney, B. R., De Young, D. S., & Roxburgh, I. W. 1993, Sol. Phys., 145, 207

Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585

González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12

Grudić, M. Y., Guszejnov, D., Offner, S. S. R., et al. 2022, MNRAS, 512, 216

Hartmann, L., & Stauffer, J. R. 1989, AJ, 97, 873

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020a, A&A, 635, A67

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Chabrier, G. 2020b, ApJ, 904, 194

Herbst, W., Eislöffel, J., Mundt, R., & Scholz, A. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 297

Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128

Kawasaki, Y., & Machida, M. N. 2023, MNRAS, 522, 3679

Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010, ApJ, 722, 1556

Laos, S., Greene, T. P., Najita, J. R., & Stassun, K. G. 2021, ApJ, 921, 110

Larson, R. B. 1969, MNRAS, 145, 271

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326

Lebreuilly, U., Hennebelle, P., Maury, A., et al. 2024a, A&A, 683, A13

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2024b, A&A, 682, A30

Lee, Y.-N., Charnoz, S., & Hennebelle, P. 2021, A&A, 648, A101

Le Gouellec, V. J. M., Greene, T. P., Hillenbrand, L. A., & Yates, Z. 2024, ApJ, 966, 91

Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

Lodato, G. 2007, Nuovo Cimento Rivista Serie, 30, 293

Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630

Machida, M. N., & Basu, S. 2019, ApJ, 876, 149

Machida, M. N., & Matsumoto, T. 2011, MNRAS, 413, 2767

Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2006, ApJ, 647, L151

Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2007, ApJ, 670, 1198

Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2008, ApJ, 676, 1088

Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2010, ApJ, 724, 1006

Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2011, PASJ, 63, 555

Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2014, MNRAS, 438, 2278

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

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

Mignon-Risse, R., González, M., & Commerçon, B. 2021a, A&A, 656, A85

Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2021b, A&A, 652, A69

Mignon-Risse, R., González, M., & Commerçon, B. 2023, A&A, 673, A134

Minerbo, G. N. 1978, J. Quant. Spectr. Radiat. Transf., 20, 541

Oliva, A., & Kuiper, R. 2023, A&A, 669, A80

Penston, M. V. 1969, MNRAS, 144, 425

Pringle, J. E. 1981, ARA&A, 19, 137

Riols, A., & Latter, H. 2019, MNRAS, 482, 3989

Saigo, K., Tomisaka, K., & Matsumoto, T. 2008, ApJ, 674, 997

Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713

Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Sullivan, C., & Kaszynski, A. 2019, J. Open Source Softw., 4, 1450

Teyssier, R. 2002, A&A, 385, 337

Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130

Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6

Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117

Toomre, A. 1964, ApJ, 139, 1217

Tsukamoto, Y., & Machida, M. N. 2013, MNRAS, 428, 1321

Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., & Inutsuka, S. 2015, MNRAS, 452, 278

Tsukamoto, Y., Machida, M. N., & Inutsuka, S.-I. 2023, PASJ, 75, 835

Tung, N.-D., Testi, L., Lebreuilly, U., et al. 2024, A&A, 684, A36

Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90

Vaytet, N., Commerçon, B., Masson, J., González, M., & Chabrier, G. 2018, A&A, 615, A5

Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627, A154

Wardle, M. 1999, MNRAS, 307, 849

Whitehouse, S. C., & Bate, M. R. 2006, MNRAS, 367, 32

Wurster, J., & Lewis, B. T. 2020, MNRAS, 495, 3807

Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 481, 2450

Wurster, J., Bate, M. R., & Bonnell, I. A. 2021, MNRAS, 507, 2354

Wurster, J., Bate, M. R., Price, D. J., & Bonnell, I. A. 2022, MNRAS, 511, 746
# Appendix A: Defining the protostar and inner disk in our simulations

The definition of the protostar in our simulations is rather semantic, as it is not a separate object from the inner disk. We use the same criterion as Vaytet et al. (2018) for defining the protostar: it is simply the gas whose density is above $10^{-5}$ g cm$^{-3}$. We note that in Vaytet et al. (2018), the authors mention that they do not have enough resolution to resolve the protostellar accretion shock separating it from the circumstellar disk, although as we have shown in Sect. 3.2, such a discontinuity does not exist.

In order to define the inner disk, none of the criteria currently used in the literature were adequate in our case, mostly due to our use of turbulent initial conditions. Additionally, in some simulations, the second collapse occurred within a larger-scale disk,

which further complicated our disk selection criterion. Inspired by Joos et al. (2012), we used the following criterion:

$$ \begin{cases} P > \rho v_r^2, & \text{Thermal support against radial infall} \\ \rho v_\phi^2 > P, & \text{Centrifugal support exceeds pressure} \\ \rho_s < \rho < 10^{-5} \text{ g cm}^{-3}, & \text{Density threshold} \end{cases} $$

where $\rho_s$ is the density of the inner disk's equatorial shock front, obtained at each snapshot through ray-tracing. This criterion ensures that no cells currently undergoing the second gravitational collapse phase are selected and that sufficient angular momentum is present in the gas to qualify as a disk. The third item in this criterion states that $\rho < 10^{-5} \text{ g cm}^{-3}$, but this is again an arbitrary choice to separate the protostar from its circumstellar disk.