# Mixing is easy: New insights for cosmochemical evolution from pre-stellar core collapse

Asmita Bhandare, Benoît Commerçon, Guillaume Laibe, Mario Flock, Rolf Kuiper, Thomas Henning, Andrea Mignone, Gabriel-Dominique Marleau

► To cite this version:

Asmita Bhandare, Benoît Commerçon, Guillaume Laibe, Mario Flock, Rolf Kuiper, et al. Mixing is easy: New insights for cosmochemical evolution from pre-stellar core collapse. Astronomy & Astrophysics - A&A, 2024, 687, 10.1051/0004-6361/202449594. insu-04838864

HAL Id: insu-04838864

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

Submitted on 15 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mixing is easy: New insights for cosmochemical evolution from pre-stellar core collapse

Asmita Bhandare<fnref n="1" />,<fnref n="2" />, Benoît Commerçon<fnref n="1" />, Guillaume Laibe<fnref n="1" />, Mario Flock<fnref n="3" />, Rolf Kuiper<fnref n="4" />, Thomas Henning<fnref n="3" />, Andrea Mignone<fnref n="5" />, and Gabriel-Dominique Marleau<fnref n="4" />,6,3,7

Received 13 February 2024 / Accepted 13 April 2024

## ABSTRACT

**Context.** Signposts of early planet formation are ubiquitous in substructured young discs. Dense, hot, and high-pressure regions that formed during the gravitational collapse process, integral to star formation, facilitate dynamical mixing of dust within the protostellar disc. This provides an incentive to constrain the role of gas and dust interaction and resolve potential zones of dust concentration during star and disc formation stages.

**Aims.** We explore whether the thermal and dynamical conditions that developed during protostellar disc formation can generate gas flows that efficiently mix and transport the well-coupled gas and dust components.

**Methods.** We simulated the collapse of dusty molecular cloud cores with the hydrodynamics code PLUTO augmented with radiation transport and self-gravity. We used a two-dimensional axisymmetric geometry and followed the azimuthal component of the velocity. The dust was treated as Lagrangian particles that are subject to drag from the gas, whose motion is computed on a Eulerian grid. We considered 1, 10, and 100 µm-sized neutral, spherical dust grains. Importantly, the equation of state accurately includes molecular hydrogen dissociation. We focus on molecular cloud core masses of 1 and 3 $M_\odot$ and explore the effects of different initial rotation rates and cloud core sizes.

**Results.** Our study underlines mechanisms for the early transport of dust from the inner hot disc regions via the occurrence of two transient gas motions, namely meridional flow and outflow. The vortical flow fosters dynamical mixing and retention of dust, while the thermal pressure driven outflow replenishes dust in the outer disc. Notably, these phenomena occur regardless of the initial cloud core mass, size, and rotation rate.

**Conclusions.** Young dynamical precursors to planet-forming discs exhibit regions with complex hydrodynamical gas features and high-temperature structures. These can play a crucial role in concentrating dust for subsequent growth into protoplanets. Dust transport, especially, from sub-au scales surrounding the protostar to the outer relatively cooler parts, offers an efficient pathway for thermal reprocessing during pre-stellar core collapse.

**Key words.** equation of state – gravitation – hydrodynamics – radiative transfer – methods: numerical – stars: formation

## 1. Introduction

The advent of unprecedented circumstellar disc surveys with the Atacama Large Millimetre/submillimetre Array (ALMA) and with facilities like the Spectro-Polarimetric High-contrast Exo-planet REsearch on the Very Large Telescope (SPHERE/VLT, e.g. Garufi et al. 2022) as well as observations with the NOrthern Extended Millimetre Array (NOEMA) highlighting environmen- tal effects (e.g. Pineda et al. 2020) are starting to shift the focus from ‘how’ planets form to ‘when’ they form. Over the last decade, with dedicated disc surveys such as Disc Substruc- tures at High Angular Resolution Project (DSHARP; Andrews et al. 2018) and Molecules with ALMA at Planet-forming Scales (MAPS; Öberg et al. 2021), there has been growing obser- vational evidence for substructures such as rings, gaps, and spiral arms via dust continuum emission from Class II discs (see also Pérez et al. 2016; Paneque-Carreño et al. 2021). More

recently, the Early Planet Formation in Embedded Discs sur- vey (eDisc; Ohashi et al. 2023) has revealed substructure even in young Class 0 and I discs (see also ALMA Partnership 2015; Sheehan & Eisner 2018; Tobin et al. 2020; Podio et al. 2020; Segura-Cox et al. 2020; Sheehan et al. 2020). These data suggest, among other things, ongoing planet formation or the presence of already formed young protoplanets carving gaps in these discs (Keppler et al. 2018; Müller et al. 2018; Gratton et al. 2019; Haffert et al. 2019; Fedele et al. 2023; Hammond et al. 2023; Pinte et al. 2023). The early forma- tion of Jupiter in our Solar System is also widely proposed as an explanation for the isotopic dichotomy and radial separa- tion between inner non-carbonaceous chondrites (NC) and outer carbonaceous chondrites (CC) (e.g. Kruijer et al. 2020).

Magnetised, cold, dense, gaseous, and dusty molecular cloud cores provide the birth environment for stars, discs, and planets.

<footnote n="1">¹ ENS de Lyon, CRAL UMR5574, Université Claude Bernard Lyon 1, CNRS, 69007 Lyon, France</footnote>

<footnote n="2">² Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote>e-mail: asmita.bhandare@lmu.de</footnote>

<footnote n="3">³ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="4">⁴ Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1, 47057 Duisburg, Germany</footnote>

<footnote n="5">⁵ Dipartimento di Fisica, Università di Torino via Pietro Giuria 1, 10125 Torino, Italy</footnote>

<footnote n="6">⁶ Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany</footnote>

<footnote n="7">⁷ Physikalisches Institut, Universität Bern, Gesellschaftsstr. 6, 3012 Bern, Switzerland</footnote>
A protostar and its surrounding protostellar disc are formed con-
currently as an outcome of the process of gravitational collapse
within these molecular cloud cores with either a non-zero net
angular momentum (e.g. Chen et al. 2023) or with turbulent
perturbations cascading from larger scales (Pineda et al. 2023).
Early collapse phases lead to high density, pressure, and temper-
ature regions within young protostellar discs (Zhao et al. 2020;
Tsukamoto et al. 2023b). These localised zones provide a unique
setting for spatial mixing and thermal reprocessing of dust early
on. Protostellar and protoplanetary discs are extremely dynam-
ical throughout their evolution due to a multitude of physical
and chemical processes (Aikawa et al. 2022; Lesur et al. 2023;
Bae et al. 2023). The inclusion of dust and gas interactions,
while treating the microphysics at ‘all’ scales, will play a major
role in understanding the earliest stages of star and disc for-
mation and their subsequent evolution. This will also allow the
exploration of the high-temperature sub-au disc region, which
is known to be an ideal location to form refractory minerals
like calcium-aluminium-rich inclusions (CAIs) and amoeboid
olivine aggregates (AOAs). Surprisingly, CAIs and AOAs are the
oldest ingredients mostly found in carbonaceous chondrites that
were accreted in the outer Solar System, implying large-scale
transport after their formation close to the proto-Sun (Grossman
et al. 1988; Krot et al. 2004; MacPherson et al. 2012; Ruzicka
et al. 2012; Desch et al. 2018; Pignatale et al. 2019; Marrocchi
et al. 2019; Hellmann et al. 2023; Piralla et al. 2023).
There remain several open questions for a comprehensive
characterisation of the multi-scale gravitational collapse within
molecular cloud cores: for example, values of the initial mag-
netic field strength and orientation, angular momentum, turbu-
lence, and dust size distribution in the cloud cores (i.e. pre-stellar
cores). These unknown initial conditions introduce various
caveats in our current understanding of the observed star-disc-
planet systems. Filling in the gaps in our current knowledge of
formation and evolution of planet-forming discs will also help
decipher cosmochemical signatures found in different parts of
our Solar System. Zooming in on the smallest sub-au scales to
unravel the effects of various complex physical processes, such
as hydrodynamics; radiative transfer; phase transition, in partic-
ular molecular hydrogen (H₂) dissociation and ionisation; tur-
bulence; chemistry; magnetic fields; and gas-dust (de)coupling,
remains a major challenge for both observations and theory
(Benisty et al. 2023; Tsukamoto et al. 2023b). Despite the recent
wealth of high-resolution observations of young discs, high opti-
cal depths in the surrounding envelope pose several difficulties
to resolve these embedded systems. On the other hand, newly
developed theoretical and numerical models help elucidate and
quantify these nascent phases of star and disc formation.

Since the pioneering work by Larson (1969), there has
been significant progress with numerical studies to follow the
collapse using both grid-based simulations (e.g. Masunaga &
Inutsuka 2000; Tomida et al. 2013; Vaytet & Haugbølle 2017;
Vaytet et al. 2018; Bhandare et al. 2018, 2020; Ahmad et al.
2023) and smoothed particle hydrodynamics (SPH) methods
(e.g. Whitehouse & Bate 2006; Stamatellos et al. 2007; Bate
et al. 2014; Wurster et al. 2018; Tsukamoto et al. 2020) as also
described in the review by Teyssier & Commerçon (2019). How-
ever, despite the vast amount of theoretical work, several aspects
of disc formation and evolution are far from being fully under-
stood, mostly due to numerical and computational constraints in
resolving the disc structure in the protostellar neighbourhood.
These depend strongly on the initial physical conditions and
numerical set-up of simulations (e.g. Machida et al. 2014). There
is no consensus on the formation time of the disc with respect to

Since the pioneering work by Larson (1969), there has
been significant progress with numerical studies to follow the
collapse using both grid-based simulations (e.g. Masunaga &
Inutsuka 2000; Tomida et al. 2013; Vaytet & Haugbølle 2017;
Vaytet et al. 2018; Bhandare et al. 2018, 2020; Ahmad et al.
2023) and smoothed particle hydrodynamics (SPH) methods
(e.g. Whitehouse & Bate 2006; Stamatellos et al. 2007; Bate
et al. 2014; Wurster et al. 2018; Tsukamoto et al. 2020) as also
described in the review by Teyssier & Commerçon (2019). How-
ever, despite the vast amount of theoretical work, several aspects
of disc formation and evolution are far from being fully under-
stood, mostly due to numerical and computational constraints in
resolving the disc structure in the protostellar neighbourhood.
These depend strongly on the initial physical conditions and
numerical set-up of simulations (e.g. Machida et al. 2014). There
is no consensus on the formation time of the disc with respect to

the first and second core formation. The relative importance of
the presence of two discs around these two cores, which could
eventually merge, is yet to be quantitatively determined. Fur-
thermore, dust grains constitute the main ingredients to form
planetary building blocks. Dust properties can also have signif-
icant cumulative dynamical effects on the thermal budget and
kinematic structure of the disc. Even so, present-day numeri-
cal codes using sophisticated techniques rarely account for the
role of drag forces and dust decoupling during protostar (sec-
ond core) and protostellar disc formation stages. The dust size
distribution and concentration in very young protostellar discs
and its influence on the radial drift remain poorly constrained.
These missing pieces are roadblocks on our way to solving the
puzzle of early planet formation, which is strongly affected by
the properties of the disc, the formation of substructures, and the
transport of material within the disc.

Historically, the effects of dust dynamics have been
neglected in star and disc formation studies based on the assump-
tion of the diffuse interstellar medium (ISM) dust-to-gas mass
ratio being one per cent. Most methods consider the dust and
gas to be strongly coupled, and use a maximum dust size of up
to 0.25 µm. The dust size distribution (number of dust grains
per grain radius *a*) is assumed to be a power law, proportional
to *a*⁻³.⁵ from Mathis, Rumpl & Nordsieck (hereafter MRN.
Mathis et al. 1977; Weingartner & Draine 2001). This has been
recently observed in the Orion A molecular cloud (Uehara et al.
2021). However, discoveries of mid-infrared light scattered from
dense core regions (Steinacker et al. 2010, 2015; Pagani et al.
2010; McClure et al. 2023) and evidence of lower dust opac-
ity spectral index in dense cores compared to the ISM (Henning
et al. 1995; Miettinen et al. 2012; Zhang et al. 2021) and in the
envelopes of young stellar objects (Kwon et al. 2009; Miotello
et al. 2014; Galametz et al. 2019) have shed light on the pres-
ence of larger dust grains of up to 1 µm in pre-stellar cores.
Dust polarisation detection due to self-scattering processes at
(sub)millimetre wavelengths have also provided indirect evi-
dence for up to 100 µm grains in young discs (Kataoka et al.
2015; Valdivia et al. 2019). Additionally, a lower spectral slope
in the spectral energy distribution derived from longer wave-
length observations of spatially resolved young discs has been
linked to the presence of millimetre-sized dust in the disc mid-
plane (e.g. Carrasco-González et al. 2016; Liu et al. 2017).
Details on different observational techniques for estimating dust
sizes in discs can be found in Natta et al. (2007), Testi et al.
(2014), and Birnstiel (2023). Modelling transitions of CO emis-
sion lines within the CO snowline suggest local enhancements
of dust-to-gas mass ratios, especially in the disc midplane (e.g.
Boneberg et al. 2016). Additionally, there has been a growing
indirect evidence for up to millimetre-sized dust lifted by out-
flows, winds, and jets (e.g. Bans & Königt 2012; Ellerbroek et al.
2014; Duchêne et al. 2024; Cacciapuoti et al. 2024).

The assumption of an MRN size distribution has also been
recently challenged by numerical studies that account for inter-
actions between gas and dust components through the drag force
in the Epstein regime where the particle size is much smaller
than the mean free path of the gas (Epstein 1924). To date,
this has been achieved either by treating the neutral dust as a
pressureless fluid (e.g. Lebreuilly et al. 2020) or as Lagrangian
particles (e.g. Bate & Lorén-Aguilar 2017; Cridland et al. 2022;
Koga et al. 2022; Koga & Machida 2023). The current frame-
work used in most simulations excludes an accurate treatment of
the thermodynamics, which should be modelled by incorporat-
ing radiation transport and hydrogen dissociation and ionisation
effects, especially important during the second collapse stage
(Teyssier & Commerçon 2019). Most methods evolve the gas pressure using a barotropic equation of state that does not account for shock heating and optical depth effects such as disc cooling in the vertical direction. Comparisons between the influence of a radiative transfer treatment versus a barotropic approximation on the temporal evolution can be found in Commerçon et al. (2010), Tomida et al. (2010, 2013), and Bate (2011).

Within these numerical restrictions, to date only a handful of research groups have tackled this challenging problem down to protostellar disc scales no smaller than 1 au (Bate & Lorén-Aguilar 2017; Vorobyov et al. 2019; Lebreuilly et al. 2019, 2020; Tsukamoto et al. 2021a,b; Koga et al. 2022; Koga & Machida 2023). These studies have used a mono- or multi-fluid approach or a Lagrangian particle treatment to capture the momentum exchange between gas and dust. Hopkins & Lee (2016), Tricco et al. (2017), and Commerçon et al. (2023) have pointed out the possible decoupling of micrometre dust grains in molecular clouds, whereas Bate & Lorén-Aguilar (2017), Vorobyov et al. (2019), Lebreuilly et al. (2019, 2020), and Tsukamoto et al. (2021a,b) have shown the decoupling of ≥100 µm grains in collapsing pre-stellar cores resulting in enhancements of the dust-to-gas mass ratio within the young protostellar disc. Numerical studies of disc evolution using a dust-to-gas mass ratio higher than one per cent have indicated a pile-up of dust, which is favourable to enhancing dust growth (Laibe 2014; Gonzalez et al. 2017; Vericel & Gonzalez 2020). A recent three-dimensional magneto-hydrodynamics study accounting for dust growth by Tsukamoto et al. (2021b) indicated that outflows on >1 au scales provide a passage to redistribute and recycle large dust grains (≥10 µm) within the disc. The distribution of dust sizes can influence the transfer of angular momentum by altering the coupling between gas and magnetic fields. This influence is attributed to the impact of dust size on resistivity (Kawasaki et al. 2022) and opacity (Ormel et al. 2009, 2011). Various studies that include detailed time-dependent chemical networks also underline the importance of dust grains as efficient charge carriers (Marchand et al. 2016; Zhao et al. 2016, 2018; Dzyurkevich et al. 2017; Tsukamoto et al. 2020). Effects of radiation transport and an accurate gas equation of state to treat the influence of H₂ dissociation and ionisation need to be investigated in combination with resistivity and angular momentum transport induced by turbulence and shear forces, while tracing the evolution of gas and dust separately.

As a first step, we account for the combined effects of hydrodynamics, self-gravity, and radiation on the gas and the influence of gas drag force on the dust. We monitor the behaviour of 1, 10, and 100 µm-sized neutral, spherical dust particles throughout the collapse. The numerical framework and initial set-up is detailed in Sect. 2. In Sects. 3 and 4, we present the results from our study aimed at determining the interplay and impact of gas and dust dynamics at pre-stellar core collapse scales, especially within the inner sub-au region. Future improvements for collapse models are mentioned in Sect. 5. Lastly, we summarise the objectives and outcomes of this project in Sect. 6.

## 2. Numerics and initial set-up

In this section we present details of the hybrid scheme where the gas was treated as a Eulerian fluid and the dust was considered as an ensemble of Lagrangian particles. The fluid-dust method accounts for aerodynamic drag forces and ensures momentum conservation.

### 2.1. Gas treatment

Two-dimensional radiation hydrodynamic (RHD) simulations were performed using a modified version of the PLUTO code (version 4.4: Mignone et al. 2007, 2012). The equations of hydrodynamics are detailed in Bhandare et al. (2018). The HAUMEA module was used for the self-gravity Poisson solver (Kuiper et al. 2010b, 2011). A variable gas equation of state from D'Angelo & Bodenheimer (2013) accounts for effects of H₂ dissociation, ionisation of atomic hydrogen and helium, and molecular vibrations and rotation (Vaidya et al. 2015).

Tabulated dust opacities from Ossenkopf & Henning (1994) and gas opacities from Malygin et al. (2014) are included. We apply a grey, frequency-averaged flux-limited diffusion (FLD) approximation for the radiation transport using the module MAKEMAKE (Kuiper et al. 2010a). This includes updates to consider a time-dependent evolution of the dust dominating at low temperatures as previously described in Kuiper et al. (2020). Gas thermodynamics is treated under the approximation of local thermodynamic equilibrium and a two-temperature approach for the gas and radiation. The two-temperature model uses a linearisation method in which the radiation and medium temperatures evolve as two different quantities (Commerçon et al. 2011). Thermal coupling between gas and dust is assumed to be perfect so that they have the same temperature (Galli et al. 2002). Note that in the MAKEMAKE module, the dust-to-gas mass ratio is fixed to one per cent meaning that the gas and dust are also dynamically well-coupled. This is seen to be a reasonable assumption since we do not find a significant variation in the dust-to-gas mass ratio (see Sect. 4.2).

The mass, momentum, and energy transport equations are solved using second-order Godunov-type schemes for the fluid wherein a shock-capturing Riemann solver is used in a conservative finite volume approach and integrated with a second-order Runge-Kutta method. We used the Harten-Lax-van Leer approximate Riemann solver that restores the middle contact discontinuity (HLLC) and a monotonised central difference flux limiter using piecewise linear interpolation for the flux computation. The FLD and Poisson equations are treated implicitly with a standard generalised minimal residual solver with approximations to the error from previous restart cycles. This is used from the Portable, Extensible Toolkit for Scientific Computation open-source solver library (PETSc; Balay et al. 1997, 2019a,b).

### 2.2. Dust particle treatment

A particle treatment of dust grains proves to be beneficial to trace their dynamical and thermal history in comparison with the Eulerian method. We treated the dust dynamics only in the Epstein regime and did not account for the back-reaction of the dust on the gas nor the interaction between dust grains. Although the back-reaction can modify the gas evolution by decreasing the inward gas velocity compared to pure viscous effects (Dipierro et al. 2018), these feedback effects are found to be negligible in a low-Stokes-number regime covered during the collapse (Lebreuilly et al. 2020, Appendix B).

For the collapse set-up herein, we have modified the dust module to account for the gas self-gravity on the dust. The dust grains are additionally affected by the gas via the drag force defined as

$$F_d = -\frac{m_{dg}}{t_s} \Delta v, \qquad (1)$$
where $m_{dg}$ is the mass of the dust grain and $\Delta\mathbf{v} = \mathbf{v}_{dg} - \mathbf{v}_{g}$ is the differential velocity between the dust grain and gas.

The dust stopping time $t_s$ is defined as

$$t_s = \sqrt{\frac{\pi\Gamma_1}{8} \frac{\rho_{dg}}{\rho_g} \frac{s_{dg}}{c_s}}, \qquad (2)$$

where $\Gamma_1$ is the first adiabatic index, $c_s$ is the sound speed, $\rho_{dg}$ is the intrinsic dust grain density, $s_{dg}$ is the dust grain size, and $\rho_g$ is the gas density.

The time integration of the dust particles equations of motion is performed using an exponential mid-point method. A logarithmic interpolation in the radial direction and a cloud-in-cell weighting scheme in the $\theta$ (polar) direction is applied for depositing particle quantities such as velocities in the grid cells and interpolating fluid properties at the particle positions as detailed in Mignone et al. (2018, 2019).

## 2.3. Initial fluid set-up

The physical and numerical set-up of the collapse is similar to the one previously described in Bhandare et al. (2018, 2020). The initial density distribution has a stable Bonnor–Ebert sphere like profile with the outer $\rho_o$ and central $\rho_c$ density determined as

$$\rho_o = \left( \frac{1.18 c_{s0}^3}{M_0 G^{3/2}} \right)^2, \qquad (3)$$

$$\rho_c = 14.1 \rho_o,$$

where $G$ is the gravitational constant and $M_0$ is the initial pre-stellar core mass (Ebert 1955; Bonnor 1956). The initial uniform sound speed $c_{s0}$ is estimated using

$$c_{s0}^2 = \frac{GM_0}{\ln(14.1) R_{\text{out}}}, \qquad (4)$$

where $R_{\text{out}}$ is the pre-stellar core radius. The density contrast between the inner and outer edge corresponds to a dimensionless radius of $\xi = 6.45$ where $\xi$ is defined as

$$\xi = \sqrt{\frac{4\pi G \rho_o}{c_{s0}^2}} R_{\text{out}}. \qquad (5)$$

The initial thermal pressure was computed for a fixed temperature $T_0$ of 10 K. However, the temperature $T_{\text{BE}}$ for the stable Bonnor–Ebert sphere is slightly higher than 10 K and is computed as

$$T_{\text{BE}} = \frac{\mu c_{s0}^2}{\gamma \mathfrak{R}}, \qquad (6)$$

where the mean molecular weight $\mu = 2.353 \text{ g mol}^{-1}$, the adiabatic index $\gamma = 5/3$, and $\mathfrak{R}$ is the universal gas constant. Thus the reduced pressure support of the stable Bonnor–Ebert sphere initiates the collapse. The radiation temperature was initially in equilibrium with the gas temperature.

The outer radius was fixed to 3000 au for the fiducial case of an initial $M_\odot$ and also for an initial $3 M_\odot$ pre-stellar core. The initial ratio of thermal-to-gravitational energy was fixed to 0.29 and 0.098 for the 1 and 3 $M_\odot$ cases, respectively and computed as

$$\frac{E_{\text{th}}}{E_{\text{grav}}} = \frac{5R_{\text{out}} k_{\text{B}} T_0}{2GM_0 \mu m_{\text{H}}}, \qquad (7)$$

<table><caption>**Table 1.** Initial pre-stellar core properties.</caption><thead><tr><th><i>M</i><sub>0</sub> (<i>M</i><sub>☉</sub>)</th><th><i>R</i><sub>out</sub> (au)</th><th><i>E</i><sub>rot</sub>/<i>E</i><sub>grav</sub></th><th>Dust size (μm)</th><th>Discussion</th></tr></thead><tbody><tr><td>1</td><td>3000</td><td>0.007</td><td>1</td><td>Sect. 3.1, A.1</td></tr><tr><td>1</td><td>3000</td><td>0.01</td><td>1</td><td>Sects. A.2, A.1</td></tr><tr><td>1</td><td>5000</td><td>0.007</td><td>1</td><td>Sects. A.3, A.1</td></tr><tr><td>3</td><td>3000</td><td>0.007</td><td>1</td><td>Sect. 3.2</td></tr><tr><td>1</td><td>3000</td><td>0.007</td><td>1, 10, 100</td><td>Sect. 4</td></tr></tbody></table>

**Notes.** Listed above are the pre-stellar core properties for runs with different initial pre-stellar core mass $M_0$ ($M_\odot$), outer radius $R_{\text{out}}$ (au), rotational-to-gravitational energy ratio $E_{\text{rot}}/E_{\text{grav}}$, fixed dust size (μm), and the corresponding discussion section.

where $k_B$ is Boltzmann constant and $m_H$ is the hydrogen mass. We compare simulations with two different values for the initial solid-body rotation rate of $\Omega_0 = 1.77 \times 10^{-13} \text{ rad s}^{-1}$ and $\Omega_0 = 2.099 \times 10^{-13} \text{ rad s}^{-1}$. The ratio of rotational-to-gravitational energy is approximated as

$$\frac{E_{\text{rot}}}{E_{\text{grav}}} = \frac{R_{\text{out}}^3 \Omega_0^2}{3GM_0} \qquad (8)$$

and for the two rotation rates is 0.007 and 0.01, respectively. These are below the disc fragmentation limit of $E_{\text{rot}}/E_{\text{grav}} > 0.01$ (Matsumoto & Hanawa 2003; Bate 2011). The fiducial run with the slower rotation rate is compared to a case with a larger pre-stellar core with an outer radius of 5000 au detailed in Appendix A.1. The varied initial pre-stellar core properties discussed here are listed in Table 1.

Mechanisms such as gravitational torques exerted from the spiral arms or bar-like structures, magneto-rotational instability, hydrodynamically driven turbulence, outflows, disc winds, and jets are responsible for the angular momentum transport during the collapse and protostellar disc formation. A physical shear viscosity was implemented to mimic some of these effects of angular momentum transport (Kuiper et al. 2010b, 2011). The shear viscosity is described using the so-called $\alpha$-parametrisation from Shakura & Sunyaev (1973) and is given by

$$\nu = \alpha c_s H \rho \approx \alpha \Omega_{\text{K}}(r) R^2 \left(\frac{H}{R}\right)^2 \rho. \qquad (9)$$

The local sound speed $c_s \approx H \Omega_{\text{K}}(r)$ with the Keplerian angular velocity $\Omega_{\text{K}}(r) \approx \sqrt{GM(r)/r^3}$. The local pressure scale height $H = (H/R)R$ where the cylindrical radius $R = r \sin(\theta)$. The dimensionless parameters $(H/R) = 0.05$ and $\alpha = 1.0$ were fixed in space and time for calculating the viscosity and to mimic the net outward angular momentum transport within self-gravitating discs. This shear viscosity is similar to the $\beta$-viscosity prescription for self-gravitating discs by Duschl et al. (2000) with a $\beta$-parameter of $\beta = 2.5 \times 10^{-3}$ for the fiducial case at time = 0. This temperature-independent $\beta$-viscosity prescription was also previously tested in the low-mass core collapse studies by Schönke & Tscharnuter (2011) and in the high-mass core collapse studies by Kuiper & Yorke (2013); Kuiper et al. (2015, 2016), and Kuiper & Hosokawa (2018).

## 2.4. Initial particle set-up

We performed simulations using a single fixed size of 1 µm dust and also a run with three fixed sizes of 1, 10, and 100 µm
initialised at the same location on the grid. We considered spher-
ical, neutral dust with an intrinsic density ρ<sub>dg</sub> of 3 g cm<sup>-3</sup> to
mimic a combination of carbonaceous and silicate grains as pre-
viously used in Tricco et al. (2017) and Lebreuilly et al. (2019).
The initial dust density was fixed to one per cent of the gas den-
sity using the conservative estimate of the averaged dust-to-gas
mass ratio of the ISM. This results in a particle mass distribu-
tion that depends on the grid cell density and volume, which is
reasonable when neglecting the dynamical back-reaction of the
dust onto the gas. In this case, the mass of the particle does not
play a role for the dynamics. In order to ensure a uniform par-
ticle sampling throughout the pre-stellar core, especially in the
innermost regions, we used 200 particles per cell for all the sim-
ulations with the fixed 1 µm dust and outer radius of 3000 au.
For the multi-size simulation, we assigned 50 particles per cell
per size for the three different sizes, that is, a total of 150 parti-
cles per cell. The particles are randomly distributed between the
cell interfaces and assigned a fraction of the grid cell mass. This
leads to a mass variation between the particles also within a grid
cell. The local and global mass conservation is ensured. Obser-
vational evidence suggests an MRN-like size distribution with a
larger population of smaller sized dust within a pre-stellar core.
In the simulations presented here, the choice of the size distribu-
tion does not affect the results since back-reaction is neglected.
As such, we choose an equal number of particles for each size in
order to compare their relative behaviour during the collapse.

2.5. Grid

We used a two-dimensional (2D) axially and midplane symmet-
ric Eulerian grid. The radial grid extends from 10⁻² au to 3000 au
and consists of 241 logarithmically spaced cells. This resolves
better the central parts of the computational domain. We used
30 uniformly spaced cells in the polar direction stretching from
the pole (θ = 0°) to the midplane (θ = 90°). The number of grid
cells is selected to have the same spatial extent in the radial and
polar direction. The smallest cell size is Δx<sub>min</sub> = Δr = rΔθ =
5.37 × 10⁻⁴ au. This resulted in 84 cells per Jeans length at the
highest central density achieved in the course of the fiducial sim-
ulation. Otherwise, we used 21–573 cells per Jeans length over
the computational domain estimated using the midplane density
and sound speed at the final simulation time snapshot.

## 2.6. Boundary conditions

At the inner radial edge, we adopted reflective boundary con-
ditions for the hydrodynamics, which sets the density, thermal
pressure, and radial velocity. A zero gradient condition was
adopted for the polar and azimuthal velocity components as well
as for the radiation energy, which implies that no radiative flux
can cross the inner boundary interface. At the outer radial edge, a
Dirichlet boundary condition was applied for the radiation tem-
perature with a constant boundary value of 10 K. Additionally,
we imposed an outflow–no-inflow condition for the hydrody-
namics, which includes a zero-gradient boundary condition for
the thermal pressure, the polar, and the azimuthal velocity com-
ponents. Axisymmetric boundaries were used at the pole and
mirror-symmetric boundaries at the equator.

**3. Results I: Features formed during the collapse**

Studies over the last ~55 yr have established that low-mass stars
form via a two-step process involving the formation of first and
second hydrostatic cores (see summary in Bhandare et al. 2018,
2020). Angular momentum conservation in these collapsing

pre-stellar cores initiates the formation of protostellar discs as
early as the first core stage. This section describes transient gas
features, namely meridional flows and outflow, formed during
the collapse process for cases with different initial conditions.

3.1. Fiducial case: 1 $M_{\odot}$ pre-stellar core collapse

The 2D RHD simulation for a 1 $M_{\odot}$ pre-stellar core collapse shows an oblate first hydrostatic core transitioning into a rotationally supported disc over a period of 1218 yr even ‘before’ the second collapse, which is the protostellar formation stage. The presence of a disc before the second collapse phase is also noted by previous hydrodynamic studies (e.g. Bate 1998, 2010, 2011; Saigo et al. 2008; Tscharnuter et al. 2009; Machida et al. 2010; Machida & Matsumoto 2011).

We extract the disc based on a combination of definitions
from Joos et al. (2012) and Koga et al. (2022). The disc is defined
as the region which satisfies all the following conditions:

- Azimuthal (rotational) velocity must be larger than twice the radial velocity $v_{\phi} > 2|v_r|$.

- Rotational support must be higher than the thermal support ($\rho v_{\phi}^2/2 > 2 P_{\text{th}}$ and $v_{\phi} > 0.6 v_{\text{K}}$, where $v_{\text{K}}$ is the Keplerian velocity).

- Gas density must be above the minimum limit set by the mid-plane density immediately outside of the first (second) shock location for the outer (inner) disc.

Figure 1 is a zoom-in within 15 au of a 3000 au pre-stellar core. The initial temperature was fixed to 10 K and the initial rotation rate to $\Omega_0 = 1.77 \times 10^{-13}$ rad s$^{-1}$ (or $E_{\text{rot}}/E_{\text{grav}} = 0.007$). Comparisons to a faster rotation case ($E_{\text{rot}}/E_{\text{grav}} = 0.01$) are detailed in the Appendix A.2 and to a larger initial cloud core (5000 au) are discussed in Appendix A.3. The dust size was fixed to a constant value of 1 µm throughout the evolution. The four panels show various dust and gas properties at 1221 yr from the start of the first core formation as it transitions into a rotationally supported disc, which is still embedded in its surrounding infalling envelope. Panel a displays the Stokes number at the dust particle location defined using the stopping time and orbital frequency as $St = t_s \times \Omega_{\text{orbital}}$, where $\Omega_{\text{orbital}} = v_\phi/r \sin(\theta)$ and $t_s$ is defined above in Eq. (2). A Stokes number much below unity indicates a strong dust–gas coupling. Gas temperature is shown in panel b with contours at a few specific temperatures for silicate sublimation (1400 K; Isella & Natta 2005; Izidoro et al. 2022), ice line (150 K), and ammonia snowline (80 K; Pinilla et al. 2017). Panel c shows the non-homologous nature of the collapse with gas density increasing towards the centre. Gas radial velocity highlighting both inward (blue) and outward (red) motion is shown in panel d. Gas falls in faster with a higher radial velocity onto the central part via the polar regions compared to the midplane. Streamlines in all four panels help visualise material falling onto the disc and mixing within, especially the meridional flow close to the first accretion shock at around 10 au. The turbulent feature moves radially outwards during the 1221 yr of evolution, while trapping some dust and enabling an outward transport. Mixing and outward gas transport via circulation currents within early disc formation stages was also previously reported in Tscharnuter et al. (2009).

Figure 2 is a further zoom-in showing different properties within 1.5 au of the 3000 au pre-stellar core at the same time snapshot as Fig. 1. The four panels show the same dust and gas properties at 3 yr after the pseudo-disc<fnref n="1" /> begins to form around the second core and is embedded within the first disc.

<footnote n="1">¹ A pseudo-disc is a dense, non-rotationally supported region.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** 2D view of the first hydrostatic core evolved into a rotationally supported disc at 1221 yr after its formation as a result of the fiducial collapse of a $1\ M_{\odot}$ pre-stellar core with an initial rotation rate of $\Omega_0 = 1.77 \times 10^{-13} \text{ rad s}^{-1}$ (same time snapshot as Fig. 2). The dust size is fixed to a constant value of 1 µm. The four panels show the (a) Stokes number, (b) gas temperature, (c) gas density, and (d) radial gas velocity within the inner 15 au of the 3000 au collapsing pre-stellar core. The gas velocity streamlines indicate the material falling onto the disc and the mixing within. The meridional flow exhibited in the plot between 5–10 au forms in the inner regions and travels outwards. This motion results in an outward transport and retains dust in the outer disc.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** 2D view of a pseudo-disc around the second core at 3 yr after its formation resulting from the same initial conditions and at the same time snapshot as Fig. 1. The vortical feature forming at sub-au scales in the vicinity of the unresolved protostar circulates dust in the high-temperature inner regions.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Polar-angle averaged radial profile of the entropy at different time snapshots during the 1 $M_\odot$ pre-stellar core collapse. A negative radial gradient generates turbulent flows within the two discs.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Thermal pressure $P_\text{th}$ in colour and isodensity contours within the first core disc at 1221 yr. The infalling envelope has been masked to highlight the non-zero baroclinicity within the disc. The gas velocity streamlines in the background indicate the material falling onto the disc and the circulation within.</figcaption></figure>

Similar examples of nested discs can be found in Saigo et al. (2008) and Machida & Matsumoto (2011). The second shock location coincides with the silicate sublimation temperature of 1400 K. A second vortical feature seen in the figure, although lasting only for a few years, helps circulate and retain dust in the high-temperature sub-au region. Further discussion on the cause and effects of both the turbulent flows can be found next in Sects. 3.1.1 and 4.1.

### 3.1.1. Formation of dust traps in turbulent gas flows

Our results show the formation of two turbulent gas flows during the 1 $M_\odot$ collapse. A meridional flow starts building up within the first core and moves radially outwards as this first core transitions into a rotationally supported disc. Another vortical feature forms within the pseudo-disc around the protostar. Figure 3 shows the radial profile of the polar-angle averaged entropy

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Baroclinicity ($s^{-2}$) at 1221 yr after the transition of the first core into a rotationally supported disc resulting from a 1 $M_\odot$ pre-stellar core collapse. The infalling envelope has been masked to highlight the baroclinic contribution in generating vortical flows within the disc. The gas velocity streamlines indicate the material falling onto the disc and the circulation within.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Solberg–Høiland unstable region within the first core disc at 1221 yr found mainly around the meridional feature.</figcaption></figure>

over different time snapshots. The first core begins to form at 18 000 yr and the pseudo-disc at 19 218 yr after the onset of the collapse. At both these times, a negative radial entropy gradient triggers a turbulent flow. In case of the first core disc, this turbulent gas flow eventually generates a baroclinically amplified vortical feature. Formation of such meridional gas flows due to a mismatch in the planes of constant density and pressure (i.e. non-zero baroclinicity) were previously reported in Klahr & Henning (1997); Klahr & Bodenheimer (2003); Klahr (2004); Petersen et al. (2007), and Lesur & Papaloizou (2010). This is depicted in Fig. 4 as a slight misalignment between isopressure surfaces and isodensity contours. Figure 5 shows the baroclinic term from the vorticity equation within the first core disc at 1221 yr. The infalling envelope surrounding the rotationally supported disc is masked to highlight the disc vortical feature. A flip in the baroclinic term occurs gradually during the radial outward motion of this meridional flow. This change of sign is visible at the
locations where there is a deviation from the inward gas flow. Additionally, in Fig. 6 we explore the origin of this baroclinicity. For that purpose, we compute the Solberg–Høiland criterion for dynamic stability of an incompressible rotating flow with respect to small isentropic disturbances (Tassoul 2000; Rüdiger et al. 2002). In cylindrical coordinates, this condition requires

$$ \frac{\partial P_{\text{th}}}{\partial z} \left( \frac{\partial j^2}{\partial R} \frac{\partial S}{\partial z} - \frac{\partial j^2}{\partial z} \frac{\partial S}{\partial R} \right) < 0, \qquad (10) $$

where $P_{\text{th}}$ is the thermal gas pressure and $j = R^2\Omega$ is the angular momentum per unit mass. The entropy $S$ is computed<fnref n="2" /> using the Sackur–Tetrode equation for a perfect gas consistent with the D’Angelo & Bodenheimer (2013) equation of state used in the simulation (see Vaidya et al. 2015 for the details of the implementation in PLUTO). It accounts for molecular, atomic, and ionised hydrogen as well as the contribution from the electrons and thus represents a straightforward extension of the expressions in Berardo et al. (2017, Appendix A). We only highlight the unstable regions found within the rotationally supported first core disc and mask out the stable envelope. On the other hand, the turbulent feature within the pseudo-disc spreads horizontally across the pseudo-disc within its short 3-yr lifetime, only to be destroyed a year later due to the launch of an outflow discussed next in Sect. 3.1.2. The short-lived nature of the pseudo-disc makes it difficult to follow the evolution and development of a baroclinic and Solberg–Høiland unstable region.

Deviation from a radial inward motion creates transient dust pockets that trap particles and circulate or sweep them outwards in both the discs before they can either settle at the midplane or drift towards the protostar. Furthermore, the longer-lived dust pocket close to the first core shock can act as a viable location for dust growth and eventually promote protoplanet formation at earlier times (see also Koga & Machida 2023). Dust traps in gas meridional flows have been previously reported in simulations resulting from strong poloidal magnetic fields in protoplanetary discs (Hu et al. 2022), inclined isodensity and isotemperature (or isopressure) surfaces (i.e. non-zero baroclinicity) in planetary envelopes (Schulik et al. 2019, 2020), and planetary spiral wakes in circumplanetary discs (Szulágyi et al. 2022).

### 3.1.2. Long-term evolution: tracing an outflow

As the collapsing core evolves further, we find a thermal pressure driven outflow launched from close to the protostar at 19 222 yr after the onset of the 1 $M_{\odot}$ pre-stellar core collapse, which is 4 yr after the second collapse and pseudo-disc formation. We follow its evolution for ≈55 yr limited by computational time expenses. An outflow is defined as the region where the radial velocity component becomes stronger than the sound speed ($v_r > c_s$) as also used in Koga et al. (2022). The polar-angle averaged radial velocity of the outflow reaches a maximum of ≈9 km s⁻¹ and decreases over its lifetime. This transient nature of the outflow is shown in Fig. 7 over a few different time snapshots. Figure 8 shows the radial Mach number during the first and last snapshots highlighting supersonic regions within the outflow. Regions with a stronger thermal pressure support defined when $\rho v_{\phi}^2/2 < 2 P_{\text{th}}$ within the outflow can be seen in Fig. 9 at different time snapshots from the onset until 55 yr after the outflow is launched. Thermal pressure building up in the pseudo-disc triggers this

<figure><img src="image_8.png" /><figcaption>Fig. 7. Polar-angle averaged radial velocity at different time snapshots showing the transient nature of the outflow during the 1 $M_{\odot}$ pre-stellar core collapse.</figcaption></figure>

outflow as seen in the topmost panel. This high-pressure, supersonic outflow destroys the vortical features in both the discs, while pushing gas and dust upwards from the midplane and radially outwards. As a result of the material being swept outwards we also find an eventual merger of the first core disc and pseudo-disc around the second core (see also Bate 1998; Machida & Matsumoto 2011). During this period, thermal pressure becomes more dominant in most parts of the initially rotationally supported first core disc. The vertical white dashed line in Fig. 9 marks the gradually increasing radius of the first core disc or merged disc defined according to the criteria listed in Sect. 3.1.

Thermal pressure driven outflows have also been reported in previous radiation hydrodynamic simulations (Schönke & Tscharnuter 2011; Bate 2011; Rodenkirch & Dullemond 2022). The closest comparison to our study are results from three-dimensional SPH simulations by Bate (2010, 2011) and 2D grid-based hydrodynamic simulations by Schönke & Tscharnuter (2011), including FLD and a realistic gas equation of state to accurately treat the effects of $H_2$ dissociation and ionisation. We focus on a qualitative comparison to these works and do not pursue a quantitative comparison due to differences in initial pre-stellar core properties, which have a strong effect on the protostellar and disc properties. Bate (2010, 2011) investigated a 1 $M_{\odot}$ pre-stellar core collapse but with different initial size (≈4679 au), rotation rates, and a uniform initial density profile. These studies report transient bipolar outflows launched after the second core formation with velocities reaching 5–10 km s⁻¹ depending on the initial rotation rate. These outflows spreading to distances of about 30–60 au are stated to be an outcome of the release of accretion energy from the protostar. It also results in an eventual merger of the inner and outer discs and a shockwave propagating radially outwards. Bate (2011) suggests that transient hydrodynamic outflows may reoccur if the accretion rate onto the stellar core increases and that this could be a source of episodic accretion and outburst events. Similarly, Schönke & Tscharnuter (2011) studied a 1 $M_{\odot}$ pre-stellar core collapse with an outer radius of 8700 au and different values of artificial viscosity to drive the angular momentum transport within the

<footnote n="2">² https://github.com/gabrielastro/idealgasentropy</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Radial Mach number ($v_r/c_s$) at first and last time snapshots of outflow formation during the 1 $M_\odot$ pre-stellar core collapse. The gas velocity is indicated via the streamlines. The colour gradients in red and blue indicate the supersonic outflow and inflow regions, respectively. The colour gradient in green shows the subsonic regions mainly within the disc.</figcaption></figure>

collapsing pre-stellar core. They also find an accretion energy driven outflow launched close to the protostar reaching distances of about 500 au. The material blown-out by the outflow is seen to be re-accreted onto the disc over an evolution period of 11 590 yr after the onset of second core formation. This long-term collapse calculation for a total of 56 kyr is achieved by restarting the simulation at 70 yr after the protostar formation with a central sink neglecting the inner 0.7 au. Our findings of a short-lived outflow launched even in hydrodynamic simulations are in agreement with both these studies.

## 3.2. Effect of initial mass: 3 $M_\odot$ pre-stellar core collapse

In this section, we describe results from a comparison simulation starting with similar initial conditions to those described

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Rotational support ($\rho v_\phi^2/2 > 2 P_{\text{th}}$) at different time snapshots highlighting the merger of the first core disc and pseudo-disc around second core formed during the 1 $M_\odot$ pre-stellar core collapse. The topmost plot shows the presence of both discs, which eventually merge after the outflow is launched. The white dashed vertical line marks the radius of the rotationally supported first core disc or merged disc using the definition detailed in Sect. 3.1. The gas velocity streamlines show the infalling envelope, mixing within the discs, and the outflowing gas.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** 2D view of a first hydrostatic core evolved into a rotationally supported disc at 533 yr after its formation as a result of the collapse of a $3 M_{\odot}$ pre-stellar core with an initial rotation rate of $\Omega_0 = 3.042 \times 10^{-13}$ rad s$^{-1}$ (same time snapshot as Fig. 11). The dust size is fixed to a constant value of 1 µm. The four panels show the (a) Stokes number, (b) gas temperature, (c) gas density, and (d) radial gas velocity within the inner 15 au of the 3000 au collapsing pre-stellar core. The gas velocity streamlines indicate the material falling onto the disc.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** 2D view of a pseudo-disc around the second core at 2 yr after its formation resulting from the same initial conditions and at the same time snapshot as Fig. 10. The vortical feature forming at sub-au scales in the vicinity of the unresolved protostar circulates dust in the high-temperature inner regions.</figcaption></figure>

so far but with a more massive initial pre-stellar core of $3 M_{\odot}$. The rotation rate fixed to $3.042 \times 10^{-13}$ rad s$^{-1}$ gives the same $E_{\text{rot}}/E_{\text{grav}} = 0.007$ as for the $1 M_{\odot}$ case. 1 µm dust is used throughout the evolution. The $3 M_{\odot}$ pre-stellar core collapse proceeds alike the $1 M_{\odot}$ but at a higher rate or a lower free-fall time.

The first hydrostatic core transitions into a rotationally supported disc and a snapshot at 533 yr after its onset is depicted in Fig. 10. A further zoom-in within 1.5 au in Fig. 11 magnifies the pseudo-disc around the second core at the same time snapshot and 2 yr into its evolution. The four panels in both these plots show the
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Polar-angle averaged radial profile of the entropy at different time snapshots during the 3 $M_\odot$ pre-stellar core collapse. A negative radial gradient generates turbulent flows within the two discs.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Polar-angle averaged radial velocity at different time snapshots showing the transient nature of the outflow during the 3 $M_\odot$ pre-stellar core collapse.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Radial Mach number ($v_r/c_s$) at first and last time snapshots of outflow formation during the 3 $M_\odot$ pre-stellar core collapse. The gas velocity is indicated via the streamlines. The colour gradients in red and blue indicate the supersonic outflow and inflow regions, respectively. The colour gradient in green shows the subsonic regions mainly within the disc.</figcaption></figure>

same gas and dust properties as in Figs. 1 and 2. The structure and properties of both the discs are different in comparison to the 1 $M_\odot$ collapse. Dust with comparatively lower Stokes number is found at the outer edge of the first core disc whereas distribution in the pseudo-disc is similar to the 1 $M_\odot$ case. The temperature within the first core disc and pseudo-disc is much higher as well. For example, note the shift in the water iceline that sits further outside the first core shock and the silicate sublimation line, which has also moved to a larger radius. Both the discs are as dense as in the 1 $M_\odot$ collapse. The first core disc has a comparatively larger vertical spread and a slightly smaller radial extent. There is only a slight deviation from the radial inward motion of the gas and no presence of a meridional flow close to the first core shock. However, a turbulent flow is present in

the pseudo-disc circulating material within the high-temperature interiors and is described in the next Sect. 3.2.1.

### 3.2.1. Formation of dust traps in turbulent gas flows

In the case of a 3 $M_\odot$ pre-stellar core collapse we find the presence of turbulent gas flows within the pseudo-disc surrounding the second core and as the pseudo-disc mergers with the outer first core disc. Figure 12 shows the polar-angle averaged radial entropy over different time snapshots. The first core begins to form at 8700 yr after the collapse is initiated. This is followed by the pseudo-disc building up around the second core at 9231 yr, which evolves for a few years. The two nested discs merge by 9252 yr. A negative entropy gradient generates turbulent gas flows within the two discs. In comparison with the 1 $M_\odot$ case, we do not find the development of a meridional flow at the outer edge of the first core disc. However, within
the inner pseudo-disc, a negative entropy gradient leads to vortical features. During the merger of the two discs this feature is destroyed only to be followed by another vortical feature that develops in the upper layers. This then moves along the merged disc surface towards the midplane before being eventually destroyed by the outflow discussed next in Sect. 3.2.2. These transient features trap, circulate, and sweep dust outwards from the high-temperature interiors. We do not find a change in the baroclinicity at the locations of these vortical features.

### 3.2.2. Long-term evolution: tracing an outflow

In comparison with the fiducial case, a thermal pressure driven outflow is launched after the two discs merge. The evolution of the outflow is followed for 61 yr, reaches a maximum of ≈6 km s⁻¹, and its velocity is seen to decrease over this time frame. The polar-angle averaged radial velocity profiles over a few different time snapshots are shown in Fig. 13. The supersonic outflow region at its start and final simulation snapshot is shown in Fig. 14. Thermal pressure builds up within the pseudo-disc at sub-au scales and the merged disc, which eventually triggers an outflow from the innermost regions. As in the fiducial case, the rotational and thermal support are compared in Fig. 15 at different time snapshots. The vertical dashed white line marks the disc radius following the definition listed in Sect. 3.1. The outflow launched at 9266 yr eventually destroys the vortical flow. During the evolution of the outflow, high thermal pressure also reduces the rotational support within the merged disc. Circulation currents also occur within the outer wings of outflow as seen in the bottom-most panel. This replenishes the outer disc region with material that has been lifted up by the outflow. Once this transient outflow is quenched, the disc is expected to transition into a rotationally supported structure.

## 4. Results II: Dust dynamics

This section focuses on the behaviour of dust particles within transient meridional flows and outflow described in the previous section. We performed a simulation with the exact same initial conditions as for the fiducial 1 $M_\odot$ pre-stellar core collapse case discussed in Sect. 3.1 but with three different dust sizes of 1, 10, and 100 µm. We assign 50 particles per size within every grid cell to compare the behaviour of different sized dust. The three fixed dust sizes share the exact same initial location.

### 4.1. Mixing and transport within high-pressure and high-temperature regions

Dust settling ($t_{dust, sett}$) and eddy turnover times ($t_{eddy}$) within the first core disc are compared in Fig. 16 and computed as

$$t_{\text{dust, sett}} = \frac{1 + \text{St}^2}{\text{St} \Omega_K}, \qquad (11)$$

$$t_{\text{eddy}} = \frac{1}{\Omega_K}, \qquad (12)$$

where Stokes number St is previously defined in Sect. 3.1 and $\Omega_K = \sqrt{GM(r)/(r \sin \theta)^3}$ (Youdin & Lithwick 2007). The infalling gas and dust from the surrounding envelope have been masked out in Fig. 16. The three plots show the behaviour of three dust sizes. As expected 100 µm dust settles faster compared to the 10 µm and 1 µm dust. However, owing to the high

<figure><img src="image_16.png" /><figcaption>Fig. 15. Rotational support ($\rho v_\phi^2/2 > 2 P_{th}$) at different time snapshots highlighting the merger of the first core disc and pseudo-disc around second core formed during the 3 $M_\odot$ pre-stellar core collapse. The top-most plot shows the presence of both discs, which eventually merge before the outflow is launched (see third plot). The white vertical line marks the radius of the rotationally supported first core disc or merged disc using the definition detailed in Sect. 3.1. The gas velocity stream-lines show the infalling envelope, mixing within the discs, and the outflowing gas.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

**Fig. 16.** Ratio of dust settling to eddy turnover time only within the first core rotationally supported disc formed during the fiducial $1 M_\odot$ pre-stellar core collapse. The three plots highlight the behaviour of different dust sizes. Also overplotted is a single particle trajectory as a black dashed line for three different sizes of dust particles starting at the same initial grid cell location. The gas velocity streamlines indicate the material falling onto the disc and the deviation of gas and dust flow from a standard radial inward motion.

gas density within the disc, all three sizes have a Stokes number lower than $10^{-2}$ and can be trapped within the meridional flow instead of moving radially inwards. As an example, trajectories for the spatial motion of different sized dust initially located in the same grid cell are marked in each plot using a black dashed line.

Regions with a stronger rotational support defined when $\rho v_\phi^2/2 > 2 P_{\text{th}}$ within the disc can be seen in Fig. 17 at different time snapshots from the start of the outflow until 55 yr after it is launched. Given the exact same initial conditions as used for the fiducial case, gas shows the same behaviour as in Fig. 9. A selected sample of 1 µm (green), 10 µm (orange), and 100 µm (blue) dust located within the very young disc are overplotted. Dust grains are well-coupled to the dense gas within both the discs irrespective of their size and hence follow the gas motion instead of gradually moving towards or accreting onto the central protostar. The innermost disc regions experience an influx of dust either directly from regions closer to the pole or via a sliding

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

**Fig. 17.** Rotational support ($\rho v_\phi^2/2 > 2 P_{\text{th}}$) at different time snapshots during the $1 M_\odot$ pre-stellar core collapse. The topmost plot shows the presence of both discs, which eventually merge after the outflow is launched. The white vertical line marks the radius of the rotationally supported first core disc or merged disc using the definition detailed in Sect. 3.1. A thermal pressure driven outflow is also seen to transport dust from the inner regions. A selected sample of 1 µm (green), 10 µm (orange), and 100 µm (blue) dust from within the young protostellar disc that move radially outwards and/or above the midplane are overplotted. The gas velocity streamlines indicate the material falling onto the disc, mixing within the two discs, and outflowing. The size-independent well-coupled dust follows the gas motion in all the different regions.
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Tracks for a selected sample of 1 µm (top), 10 µm (middle), and 100 µm (bottom) dust that end up within the two discs during the 1 $M_\odot$ pre-stellar core collapse. The colourmap shows the Stokes number, which changes as the fixed-size dust moves through different density regimes within the envelope, discs, and outflow. The vertical lines indicate the formation times of the two discs, and the outflow, and merger of the discs from the onset of the first core formation ($t=0$). The horizontal dashed black line marks the first core disc or merged disc radius defined using the conditions listed in Sect. 3.1. A Stokes number much smaller than one throughout the evolution in all three plots indicates the size-independent coupling to the gas.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Tracks for a selected sample of 1 µm dust that end up within the two discs during the 1 $M_\odot$ pre-stellar core collapse. The colourmap is split below and above silicate sublimation at 1400 K and indicates the gas temperature at dust location. The vertical lines indicate the formation times of the two discs, and the outflow, and merger of the discs from the onset of the first core formation ($t=0$). The horizontal dashed black line marks the first core disc or merged disc radius defined using the conditions listed in Sect. 3.1. The temporal evolution of the dust provides a clear indication of thermal reprocessing within the young protostellar disc. Tracks for 10 and 100 µm can be found in Appendix B for comparison.</figcaption></figure>

effect over the disc surface. Once the outflow is launched dust particles from the innermost region are lifted upwards and/or swept radially outwards. The outflow velocity decreases over time (see Fig. 7) and will eventually be quenched, returning the entrained dust back in the disc and potentially at larger radii. Similar to our findings, previous works have also suggested (Bate 2010) and shown protostellar disc outflows (Wong et al. 2016; Tsukamoto et al. 2021b; Koga & Machida 2023) and magneto-centrifugal disc winds (Giacalone et al. 2019) to be carriers of $\ge 10$ µm dust.

The Lagrangian particle treatment allows us to trace the thermal history of different dust grains. Figures 18 and 19 display tracks for a selected sample of dust that end up within the two discs during the collapse and their dynamical changes due to the turbulent flows and outflow. The time $t=0$ in these figures marks the beginning of the first hydrostatic core formation. Different vertical dashed lines mark the times of onset of the two discs, outflow, and merger of the discs.

Figure 18 highlights the change in the Stokes number as the dust grains move through the envelope and mix within the disc and outflow regions with different densities. The Stokes number is found to be well below unity throughout the evolution, even for the largest 100 µm grains, which indicates a size independent strong coupling to the gas.

During post-processing we locate particles closest to a grid cell centre and tag them with gas properties such as temperature from the corresponding cell. In Fig. 19 and Appendix B, gas temperature at the dust location is shown in colour with two different colourmaps differentiating between the temperature below and above silicate sublimation at 1400 K. A thermal pressure driven outflow and outward gas motion during the merger of the nested discs transports high-temperature sublimated dust from the innermost sub-au parts of the disc to outer regions, while cooling it down to below 1400 K. The hotter innermost regions of the disc provide ideal conditions to condensate CAIs and AOAs. Subsequent mixing and outward transport of dust via vortical gas flows and outflow could act as potential mechanisms that enable the presence of CAIs and AOAs embedded within carbonaceous chondrites found in the outer Solar System (Morfill et al. 1978;

Tscharnuter et al. 2009; Haugbølle et al. 2019; Tsukamoto et al. 2021b; Jongejan et al. 2023).

## 4.2. Dust mass to gas mass variations during collapse

In Fig. 20, we plot the ratio of the dust mass and enclosed gas mass summed within the nested and merged discs, infalling envelope, and the outflow, normalised to its initial value. We use global values within different regions since we observe that the local dust density is highly sensitive to the particle sampling and grid resolution for a Lagrangian treatment (Commerçon et al. 2023). Time $t=0$ in these plots indicates the beginning of transition of the first core into a rotationally supported disc. As the collapse forms a protostar-disc-outflow system, we find an increase in the normalised dust mass to gas mass ratio within the disc until the outflow is launched. As more of the 100 µm (blue) dust settles in the disc, an increase in the ratio is slightly higher compared to the 10 µm (orange) and 1 µm (green) dust (see also Fig. 21). During the same period, dust is depleted in the envelope and the ratio slowly decreases over time until the outflow is triggered. As the outflow evolves while lifting dust from the disc there is a slight enhancement of the ratio in the outflow and a corresponding depletion in the disc. As mentioned earlier in Sect. 3.1.2, the outflow velocity decreases over time marking its transient nature, which enables a replenishment of dust back into the disc. This is seen as a decrease in the normalised dust mass to gas mass ratio within the outflow and a subsequent increase in the disc. The outflow shock front also prevents the infalling dust to be transported in the disc and creates a short-term increase in the normalised dust mass to gas mass ratio in the surrounding envelope.

## 4.3. Comparisons with previous work

Recent work by Koga et al. (2022) and Koga & Machida (2023) made use of a fluid method for the gas evolution and a Lagrangian particle treatment for the dust (without dust feedback) on a three-dimensional nested grid. This study treats the effects of non-ideal magneto-hydrodynamics via the Ohmic dissipation term and uses a barotropic equation of state. They
<figure><img src="image_21.png" /></figure>

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

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

**Fig. 20.** Time evolution of the dust mass to gas mass ratio normalised to its initial value for the three dust sizes within the two nested and merged discs (top), infalling envelope (middle), and outflow (bottom) formed during the collapse of a 1 $M_{\odot}$ pre-stellar core. Time $t = 0$ indicates the beginning of transition of the first hydrostatic core into a rotationally supported disc. The vertical dashed magenta line marks the launch of the outflow.

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

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

**Fig. 21.** Fraction of dust particles within the disc and outflow over time with respect to its reservoir. Top: time evolution of the fraction of dust particles within the disc with respect to the total number of particles used in the simulation. Bottom: time evolution of the fraction of dust particles entrained in the outflow with respect to those present in the disc just before the outflow is launched. Time $t = 0$ indicates the beginning of transition of the first hydrostatic core into a rotationally supported disc formed during the collapse of a 1 $M_{\odot}$ pre-stellar core. The vertical dashed magenta line marks the launch of the outflow.

follow the evolution of the protostar–disc–outflow system for 85 000 yr and track the motion of neutral dust particles (0.01–1000 µm) throughout the collapse spanning the 6130 au-wide envelope down to 1 au. During a 1.25 $M_{\odot}$ pre-stellar core collapse they also find $\le 10$ µm dust to be very well coupled to the gas, in agreement with Bate & Lorén-Aguilar (2017) and Lebreuilly et al. (2020). Larger $\ge 100$ µm dust concentrates comparatively faster in the central highest density region and within the rotationally supported ~20 au disc due to a high Stokes number and hence weak coupling in the low-density envelope. Additionally, dust up to 100 µm is entrained in the outflow driven from the disc surface reaching ~1000 au. They report a size-independent dust motion within the young dense protostellar disc controlled by the gas flow, similar to our findings. In their three-dimensional simulation, a radial inward drift outside of the spiral structure found in the outer disc is prevented as this region receives angular momentum from the inner part via the gravitational instability. On the other hand, in our
two-dimensional simulation, we find dust concentration in the outer disc via meridional gas flows.

In Koga et al. (2022) and Koga & Machida (2023), dust-to-gas mass ratio normalised to the initial value is computed by introducing gas tracer particles and mass weighting. In order to overcome the limitations of the Lagrangian particle treatment, the dust-to-gas mass ratio is summed up within specific regions (such as disc, outflow, envelope) and, similarly to our treatment, is not evaluated at each location. In agreement with the results presented here, a variation within 10% from the initial value is seen in the dust-to-gas mass ratio for ≤10 µm grains. On the other hand, 100 µm grains show a larger increase in the dust-to-gas mass ratio within the disc and outflow and a stronger depletion in the envelope. In both studies using a Lagrangian particle treatment, dust-to-gas mass ratio variation in different regions is a few times lower than that reported in Lebreuilly et al. (2020) who treat the dust using a fluid method.

## 5. Future additions to the collapse models

The gravitational collapse scenario of a pre-stellar core is outlined through the formation stages of the first core transitioning into a disc, protostar (or second core), a second pseudo-disc, and an outflow driven by thermal pressure. Here we focus on the effects of self-gravity, radiation transport, an accurate gas equation of state to model H₂ dissociation, and solid-body rotation during these different evolutionary phases. The influence of the fluid gas drag on different sized dust particles is also detailed. Future improvements to our current set-up will be able to provide further understanding of complex physical processes affecting gas and dust dynamics within young protostellar discs.

The importance of non-ideal magneto-hydrodynamics effects (Lesur et al. 2023; Tsukamoto et al. 2023b) on gas and dust interaction in two- and three-dimensional simulations, especially on formation of vortical gas features and outflows, will be studied in our follow-up work. Modelling long-term disc evolution (Vorobyov et al. 2018) through the episodic accretion phase is essential to study the origin of FU Orionis and EX Orionis outbursts that have been widely observed in low-mass star-disc systems (e.g. Cieza et al. 2018) and their impact on dust dynamics. Additionally, the consequences of aerodynamic drag forces including dust back-reaction (Hennebelle & Lebreuilly 2023) and the Lorentz force for charged dust on formation of substrates and instabilities through disc evolution stages and for a dust size distribution need to be investigated.

A further step to interpret the influence of dust evolution (Birnstiel et al. 2012; Testi et al. 2014; Birnstiel 2023) will be to combine ongoing independent efforts to include charged dust species (e.g. Tsukamoto et al. 2021a), dust growth (e.g. Tsukamoto et al. 2021b, 2023a; Beitia-Antero & Gómez de Castro 2021; Tu et al. 2022; Vorobyov et al. 2024), coagulation (e.g. Guillet et al. 2020; Tominaga et al. 2021; Marchand et al. 2022; Bate 2022), fragmentation (e.g. Kawasaki et al. 2022; Kawasaki & Machida 2023), evaporation (e.g. Grassi et al. 2017), and also the porosity (e.g. Ormel et al. 2009, 2011) of the dust in the current collapse models. These effects can significantly alter the coupling between magnetic fields and gas flow via changes in opacity, resistivity, and ionisation and in turn impact disc formation, disc evolution, and dust dynamics.

## 6. Summary

A self-consistent numerical analysis of the role of gas and dust dynamics in the context of gravitational pre-stellar core collapse

has received attention only over the last few years and is slowly gaining momentum. In this study we shed some light on the effect of gas drag on dust transport during protostar and protostellar disc formation stages. We explored different phases of a collapsing pre-stellar core with a microphysical lens accounting for interactions between fluid gas and Lagrangian dust. We performed two-dimensional radiation hydrodynamics simulations for different initial pre-stellar core masses, sizes, and rotation rates. We followed the evolution for a few tens of years after the protostellar formation, while resolving the innermost sub-au regions (i.e. no sink). Spatial and temporal evolution of neutral, spherical dust with three fixed sizes of 1, 10, and 100 µm was traced through the first core, disc, protostar, and outflow formation stages. Our main findings are summarised below:

- Pressure and temperature within young protostellar discs are much higher than the values assumed for an isolated cold nebula, thereby suggesting revisions to cosmological findings.

- Qualitatively, for all initial pre-stellar core configurations discussed here, disc formation starts even 'before' the formation of the protostar. Following the transition of the first core into a rotationally supported disc, another pseudo-disc builds up around the forming protostar. The nested inner and outer discs eventually merge into a single disc. During further evolution, this protostar-disc system is also disrupted by a transient outflow, launched in the vicinity of the protostar and driven by thermal pressure.

- The order of the disc merger and outflow launching depends on the initial pre-stellar core mass. The evolutionary pathway for a $1 M_\odot$ collapse unfolds as follows: First core; first core disc; protostar and second inner pseudo-disc; outflow; inner and outer disc merger. Using the same initial pre-stellar core properties but with a higher mass of 3 $M_\odot$ the two discs merge before the onset of the outflow. The sequence then becomes first core; first core disc; protostar and second inner pseudo-disc; inner and outer disc merger; outflow.

- The size, mass, and lifetime of the young protostellar disc have a strong dependence on the initial pre-stellar core size, mass, and ratio of rotational-to-gravitational energy.

- The early stages of disc formation and evolution are extremely dynamic during which the dust mostly follows the gas flow. The low Stokes number (far below unity) due to a high gas density within the disc keeps even large dust (up to 100 µm) very well coupled to the gas. As the pre-stellar core evolution proceeds through the first and second collapse stages, we find the presence of two transient gas features, namely meridional gas flows and outflow, that enable mixing and transport of dust and gas within the nested discs.

- Meridional gas flows: In the case of a $1 M_\odot$ collapse a meridional flow is generated at the outer edge of the first core disc and another vortical feature within the inner pseudo-disc. On the other hand, in the $3 M_\odot$ case, this transient dust trap only appears in the inner pseudo-disc. These turbulent flows are a consequence of a negative radial entropy gradient. The meridional gas flow at the outer edge of the first core disc is further enhanced due to a non-zero baroclinicity that develops over its evolution. Such meridional gas flows circulating material for a few hundred years in the outer disc and vortical features in the inner sub-au regions lasting for a few years can act as ideal locations to concentrate and mix well-coupled dust. This can provide building material for early dust growth at varied locations within the protostellar disc leading to a heterogeneous planet population.
- Outflow: a thermal pressure driven outflow is launched a few years after the protostellar birth. This outflow is seen to transport up to 100 µm dust upwards to distances of few tens of au and radially outwards until the outer edge of the first core disc. The entrained dust can be resupplied back into the very young disc due to the short lifetime of the outflow. This provides a possible mechanism to transport dust from the hotter (≥1400 K) innermost sub-au regions to the colder (≥150 K) outer parts of the disc. An analogue of this remixing could be a solution to condensate, transport, and store CAIs and AOAs mostly found in the outer Solar System.

Our results showcase fresh insights from pre-stellar core collapse for gas evolution and its influence on dust dynamics during the birth of a protostar and protostellar disc that could have significant implications for thermal reprocessing, cosmochemical evolution, and early protoplanet formation.

**Acknowledgements.** The authors thank the referee for their positive and constructive feedback. AB would like to thank Manuel Riener, Hans Baehr, Francesco Lovascio, David Melon Fuksman, Hubert Klahr, Yusuke Tsukamoto, Kengo Tomida, Yves Marrocchi, Enrico Ragusa, and Benedetta Veronesi for helpful comments and discussions during this project. We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. This project was partly supported by the IDEX-Lyon Turbulent project (contract ANR-16-IDEX-0005). AB and GL acknowledge funding from the ERC CoG project PODCAST No. 864965. AB also acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311. BC acknowledges funding from the French Agency Nationale de la Recherche (ANR) through the projects DISKBUILD (ANR-20-CE49-0006) and PROMETHEE (ANR-22-CE31-0020). GDM acknowledges the support of the DFG priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (MA 9185/1) and from the Swiss National Science Foundation under grant 200021_204847 "PlanetsInTime". Parts of this work have been carried out within the framework of the NCCR Planets supported by the Swiss National Science Foundation. Data availability: https://asmitabhandare.github.io/#tara. Software: Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), seaborn (Waskom 2021), and Jupyter Notebook (Kluyver et al. 2016).

## References

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

Aikawa, Y., Okuzumi, S., & Pontoppidan, K. 2022, arXiv e-prints, [arXiv:2212.14529]

ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41

Bae, J., Isella, A., Zhu, Z., et al. 2023, in Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, Astronomical Society of the Pacific Conference Series, 534, 423

Balay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F. 1997, in Modern Software Tools in Scientific Computing, eds. E. Arge, A. M. Bruaset, & H. P. Langtangen (Birkhäuser Press), 163

Balay, S., Abhyankar, S., Adams, M. F., et al. 2019a, PETSc Users Manual Tech. Rep. ANL-95/11 – Revision 3.11, Argonne National Laboratory

Balay, S., Abhyankar, S., Adams, M. F., et al. 2019b, PETSc Web page, https://www.mcs.anl.gov/petsc

Bans, A., & Königl, A. 2012, ApJ, 758, 100

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

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

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

Bate, M. R. 2022, MNRAS, 514, 2145

Bate, M. R., & Lorén-Aguilar, P. 2017, MNRAS, 465, 1089

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

Beitia-Antero, L., & Gómez de Castro, A. I. 2021, ApJ, 909, 206

Benisty, M., Dominik, C., Follette, K., et al. 2023, in Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, Astronomical Society of the Pacific Conference Series, 534, 605

Berardo, D., Cumming, A., & Marleau, G.-D. 2017, ApJ, 834, 149

Bhandare, A., Kuiper, R., Henning, T., et al. 2018, A&A, 618, A95

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

Birnstiel, T. 2023, ARA&A, submitted [arXiv:2312.13287]

Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

Boneberg, D. M., Panić, O., Haworth, T. J., Clarke, C. J., & Min, M. 2016, MNRAS, 461, 385

Bonnor, W. B. 1956, MNRAS, 116, 351

Cacciapuoti, L., Testi, L., Podio, L., et al. 2024, ApJ, 961, 90

Carrasco-González, C., Henning, T., Chandler, C. J., et al. 2016, ApJ, 821, L16

Chen, C.-Y., Friesen, R., Li, J., et al. 2023, MNRAS, submitted, [arXiv:2312.12504]

Cieza, L. A., Ruíz-Rodríguez, D., Perez, S., et al. 2018, MNRAS, 474, 4347

Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3

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

Commerçon, B., Lebreuilly, U., Price, D. J., et al. 2023, A&A, 671, A128

Cridland, A. J., Rosotti, G. P., Tabone, B., et al. 2022, A&A, 662, A90

D'Angelo, G., & Bodenheimer, P. 2013, ApJ, 778, 77

Desch, S. J., Kalyaan, A., & O'D. Alexander, C. M. 2018, ApJS, 238, 11

Dipierro, G., Laibe, G., Alexander, R., & Hutchison, M. 2018, MNRAS, 479, 4187

Duchêne, G., Ménard, F., Stapelfeldt, K. R., et al. 2024, AJ, 167, 77

Duschl, W. J., Strittmatter, P. A., & Biermann, P. L. 2000, A&A, 357, 1123

Dzyurkevich, N., Commerçon, B., Lesaffre, P., & Semenov, D. 2017, A&A, 603, A105

Ebert, R. 1955, ZAp, 37, 217

Ellerbroek, L. E., Podio, L., Dougados, C., et al. 2014, A&A, 563, A87

Epstein, P. S. 1924, Phys. Rev., 23, 710

Fedele, D., Bollati, F., & Lodato, G. 2023, A&A, 672, A125

Galametz, M., Maury, A. J., Valdivia, V., et al. 2019, A&A, 632, A5

Galli, D., Walmsley, M., & Gonçalves, J. 2002, A&A, 394, 275

Garufi, A., Dominik, C., Ginski, C., et al. 2022, A&A, 658, A137

Giacalone, S., Teitler, S., Königl, A., Krijt, S., & Ciesla, F. J. 2019, ApJ, 882, 33

Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

Grassi, T., Bovino, S., Haugbølle, T., & Schleicher, D. R. G. 2017, MNRAS, 466, 1259

Gratton, R., Ligi, R., Sissa, E., et al. 2019, A&A, 623, A140

Grossman, J. N., Rubin, A. E., Nagahara, H., & King, E. A. 1988, in Meteorites and the Early Solar System, eds. J. F. Kerridge, & M. S. Matthews (Tucson, AZ: University of Arizona Press), 619

Guillet, V., Hennebelle, P., Pineau des Forêts, G., et al. 2020, A&A, 643, A17

Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749

Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Haugbølle, T., Weber, P., Wielandt, D. P., et al. 2019, AJ, 158, 55

Hellmann, J. L., Schneider, J. M., Wölfer, E., et al. 2023, ApJ, 946, L34

Hennebelle, P., & Lebreuilly, U. 2023, A&A, 674, A149

Henning, T., Michel, B., & Stognienko, R. 1995, Planet. Space Sci., 43, 1333

Hopkins, P. F., & Lee, H. 2016, MNRAS, 456, 4174

Hu, X., Li, Z.-Y., Zhu, Z., & Yang, C.-C. 2022, MNRAS, 516, 2006

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Isella, A., & Natta, A. 2005, A&A, 438, 899

Izidoro, A., Dasgupta, R., Raymond, S. N., et al. 2022, Nat. Astron., 6, 357

Jongejan, S., Dominik, C., & Dullemond, C. P. 2023, A&A, 679, A45

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

Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78

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

Kawasaki, Y., Koga, S., & Machida, M. N. 2022, MNRAS, 515, 2072

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Klahr, H. 2004, ApJ, 606, 1070

Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869

Klahr, H. H., & Henning, T. 1997, Icarus, 128, 213

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides & B. Schmidt (IOS Press), 87

Koga, S., & Machida, M. N. 2023, MNRAS, 519, 3595

Koga, S., Kawasaki, Y., & Machida, M. N. 2022, MNRAS, 515, 6073

Krot, A. N., Petaev, M. I., Russell, S. S., et al. 2004, Chem. Erde/Geochemistry, 64, 185

Kruijer, T. S., Kleine, T., & Borg, L. E. 2020, Nat. Astron., 4, 32

Kuiper, R., & Hosokawa, T. 2018, A&A, 616, A101

Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61

Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010a, A&A, 511, A81

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

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20

Kuiper, R., Yorke, H. W., & Turner, N. J. 2015, ApJ, 800, 86

Kuiper, R., Turner, N. J., & Yorke, H. W. 2016, ApJ, 832, 40

Kuiper, R., Yorke, H. W., & Mignone, A. 2020, ApJS, 250, 13
Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.-F., & Kemball, A. J. 2009, ApJ, 696, 841

Laibe, G. 2014, MNRAS, 437, 3037

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

Lebreuilly, U., Commerçon, B., & Laibe, G. 2019, A&A, 626, A96

Lebreuilly, U., Commerçon, B., & Laibe, G. 2020, A&A, 641, A112

Lesur, G., & Papalozou, J. C. B. 2010, A&A, 513, A60

Lesur, G., Flock, M., Ercolano, B., et al. 2023, in Astronomical Society of the Pacific Conference Series, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, Astronomical Society of the Pacific Conference Series, 534, 465

Liu, Y., Henning, T., Carrasco-González, C., et al. 2017, A&A, 607, A74

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

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

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

MacPherson, G., Kita, N., Ushikubo, T., Bullock, E., & Davis, A. 2012, Earth Planet. Sci. Let., 331, 43

Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, A91

Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18

Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2022, A&A, 666, A27

Marrocchi, Y., Villeneuve, J., Jacquet, E., Piralla, M., & Chaussidon, M. 2019, PNAS, 116, 23461

Masunaga, H., & Inutsuka, S.-i. 2000, ApJ, 531, 350

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

Matsumoto, T., & Hanawa, T. 2003, ApJ, 595, 913

McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, Nat. Astron., 7, 431

Miettinen, O., Harju, J., Haikala, L. K., & Juvela, M. 2012, A&A, 538, A137

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7

Mignone, A., Bodo, G., Vaidya, B., & Mattia, G. 2018, ApJ, 859, 13

Mignone, A., Flock, M., & Vaidya, B. 2019, ApJS, 244, 38

Miotello, A., Testi, L., Lodato, G., et al. 2014, A&A, 567, A32

Morfill, G., Roeser, S., Voelk, H., & Tscharnuter, W. 1978, Moon and Planets, 19, 211

Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, A2

Natta, A., Testi, L., Calvet, N., et al. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 767

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Ohashi, N., Tobin, J. J., Jørgensen, J. K., et al. 2023, ApJ, 951, 8

Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845

Ormel, C. W., Min, M., Tielens, A. G. G. M., Dominik, C., & Paszun, D. 2011, A&A, 532, A43

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Pagani, L., Steinacker, J., Bacmann, A., Stutz, A., & Henning, T. 2010, Science, 329, 1622

Paneque-Carreño, T., Pérez, L. M., Benisty, M., et al. 2021, ApJ, 914, 88

Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519

Petersen, M. R., Julien, K., & Stewart, G. R. 2007, ApJ, 658, 1236

Pignatale, F. C., Jacquet, E., Chaussidon, M., & Charnoz, S. 2019, ApJ, 884, 31

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158

Pineda, J. E., Arzoumanian, D., Andre, P., et al. 2023, in Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, Astronomical Society of the Pacific Conference Series, 534, 233

Pinilla, P., Pohl, A., Stammler, S. M., & Birnstiel, T. 2017, ApJ, 845, 68

Pinte, C., Hammond, I., Price, D. J., et al. 2023, MNRAS, 526, L41

Piralla, M., Villeneuve, J., Schnuriger, N., Bekaert, D. V., & Marrocchi, Y. 2023, Icarus, 394, 115427

Podio, L., Garufi, A., Codella, C., et al. 2020, A&A, 642, L7

Rodenkirch, P. J., & Dullemond, C. P. 2022, A&A, 659, A42

Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781

Ruzicka, A., Floss, C., & Hutson, M. 2012, Geochim. Cosmochim. Acta, 79, 79

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

Schönke, J., & Tscharnuter, W. M. 2011, A&A, 526, A139

Schulik, M., Johansen, A., Bitsch, B., & Lega, E. 2019, A&A, 632, A118

Schulik, M., Johansen, A., Bitsch, B., Lega, E., & Lambrechts, M. 2020, A&A, 642, A187

Segura-Cox, D. M., Schmiedeke, A., Pineda, J. E., et al. 2020, Nature, 586, 228

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

Sheehan, P. D., & Eisner, J. A. 2018, ApJ, 857, 18

Sheehan, P. D., Tobin, J. J., Federman, S., Megeath, S. T., & Looney, L. W. 2020, ApJ, 902, 141

Stamatellos, D., Whitworth, A. P., Bisbas, T., & Goodwin, S. 2007, A&A, 475, 37

Steinacker, J., Pagani, L., Bacmann, A., & Guieu, S. 2010, A&A, 511, A9

Steinacker, J., Andersen, M., Thi, W. F., et al. 2015, A&A, 582, A70

Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1

Tassoul, J.-L. 2000, Stellar Rotation (New York, NY: Cambridge University Press)

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339

Teyssier, R., & Commerçon, B. 2019, Front. Astron. Space Sci., 6, 51

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

Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2010, ApJ, 714, L58

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

Tominaga, R. T., Inutsuka, S.-i., & Kobayashi, H. 2021, ApJ, 923, 34

Tricco, T. S., Price, D. J., & Laibe, G. 2017, MNRAS, 471, L52

Tscharnuter, W. M., Schönke, J., Gail, H. P., Trieloff, M., & Lüttjohann, E. 2009, A&A, 504, 109

Tsukamoto, Y., Machida, M. N., Susa, H., Nomura, H., & Inutsuka, S. 2020, ApJ, 896, 158

Tsukamoto, Y., Machida, M. N., & Inutsuka, S. 2021a, ApJ, 913, 148

Tsukamoto, Y., Machida, M. N., & Inutsuka, S.-i. 2021b, ApJ, 920, L35

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

Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023b, in Astronomical Society of the Pacific Conference Series, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, Astronomical Society of the Pacific Conference Series, 534, 317

Tu, Y., Li, Z.-Y., & Lam, K. H. 2022, MNRAS, 515, 4780

Uehara, H., Dobashi, K., Nishiura, S., Shimoikura, T., & Naoi, T. 2021, ApJ, 915, 74

Vaidya, B., Mignone, A., Bodo, G., & Massaglia, S. 2015, A&A, 580, A110

Valdivia, V., Maury, A., Brauer, R., et al. 2019, MNRAS, 488, 4897

Vaytet, N., & Haugbølle, T. 2017, A&A, 598, A116

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

Vericel, A., & Gonzalez, J.-F. 2020, MNRAS, 492, 210

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Vorobyov, E. I., Akimkin, V., Stoyanovskaya, O., Pavlyuchenkov, Y., & Liu, H. B. 2018, A&A, 614, A98

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

Vorobyov, E., Kulikov, I., Elbakyan, V., McKevitt, J., & Guedel, M. 2024, A&A, 683, A202

Waskom, M. L. 2021, J. Open Source Software, 6, 3021

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

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

Wong, Y. H. V., Hirashita, H., & Li, Z.-Y. 2016, PASJ, 68, 67

Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 475, 1859

Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588

Zhang, C.-P., Launhardt, R., Liu, Y., Tobin, J. J., & Henning, T. 2021, A&A, 646, A18

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050

Zhao, B., Caselli, P., & Li, Z.-Y. 2018, MNRAS, 478, 2723

Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 43
# Appendix A: Dependence on initial pre-stellar core properties

## A.1. Disc properties for the fiducial $1 M_\odot$ case

Protostellar disc properties are compared for cases of a $1 M_\odot$ collapse with different initial ratio of rotational-to-gravitational energy (or rotation rate) and pre-stellar core size until they reach the same central density of $0.01 \text{ g cm}^{-3}$. The evolution is shown from the beginning ($t=0$) of transition of the first hydrostatic core into a rotationally supported disc. For the same initial conditions except the ratio of rotational-to-gravitational energy (compare red and yellow lines), the disc radius, lifetime, and the final enclosed disc mass are larger for a faster rotating pre-stellar core in agreement with results from Machida et al. (2010), Bate (2011), and Machida & Matsumoto (2011).

A larger initial pre-stellar core with the same ratio of rotational-to-gravitational energy (but different rotation rate) also leads to longer lived and larger disc (compare blue and yellow lines). However, a slower collapse for an initially larger pre-stellar core leads to a slower build up of disc mass. These comparisons emphasise the strong dependence of protostellar disc properties on the initial physical conditions.

## A.2. Dependence on initial rotation rate

In order to investigate the dependence on the initial rotation rate we performed a simulation with the same initial conditions as for the fiducial case discussed in Sect. 3.1 but with a faster initial rotation $\Omega_0 = 2.099 \times 10^{-13} \text{ rad s}^{-1}$. The ratio of rotational-to-gravitational energy was fixed to 0.01. Dust size was fixed to 1 µm throughout the evolution. Figure A.2 shows a 2D view of the first hydrostatic core evolved into a rotationally supported disc at 1427 years after the onset of first core formed during a $1 M_\odot$ collapse. Both Fig. 1 and Fig. A.2 are shown at the snapshot when the central density reaches $\approx 0.01 \text{ g cm}^{-3}$. An initially faster rotating pre-stellar core leads to a comparatively flattened and slightly larger disc. The disc radius and mass for the two runs are compared in the previous Appendix A.1. However, different gas and dust properties within the first core disc in both cases are similar, with the meridional flow seen at the outer edge. Figure A.3 is a further zoom-in within 1.5 au at the same time snapshot as Fig. A.2 magnifying the turbulent flow within the pseudo-disc and the onset of a thermal pressure driven outflow. We do not pursue a further evolution of the outflow but expect a similar behaviour as in the fiducial case.

## A.3. Dependence on initial pre-stellar core radius

All the simulations discussed in the main text were performed with an outer pre-stellar core radius of 3000 au. Here we discuss the effects of initiating a $1 M_\odot$ collapse with a 5000 au outer radius but setting all other cloud core properties to be the same as the fiducial case (see Sect. 3.1). The radial grid extends from $10^{-2}$ au to 5000 au and consists of 167 logarithmically spaced cells. We used 20 uniformly spaced cells in the polar direction stretching from the pole ($\theta = 0^\circ$) to the midplane ($\theta = 90^\circ$). In this case, we used 50 dust particles per cell (167000 particles in total) with a fixed size of 1 µm.

Figure A.4 shows different quantities at 1646 years after the onset of the first core. The time snapshot is selected such that the central density is $0.01 \text{ g cm}^{-3}$ to compare with the fiducial case shown in Fig. 1. Similar meridional gas flow acting as a dust trap at the outer disc edge is visualised in the radial velocity

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

**Fig. A.1.** Time evolution of the nested protostellar disc radius (top) and enclosed mass (bottom) formed during the collapse of a $1 M_\odot$ pre-stellar core. Comparisons between two different initial rotation rates of $\alpha = 0.007$ and 0.01 and an initial pre-stellar core radius of 3000 au vs 5000 au are shown. Time $t=0$ indicates the beginning of transition of the first hydrostatic core into a rotationally supported disc. The three cases are compared until central density reaches $0.01 \text{ g cm}^{-3}$ and before the launch of an outflow.

panel. Figure A.1 highlights the dependence of the disc radius, mass, and lifetime on the initial pre-stellar core size. A larger initial pre-stellar core leads to a slower collapse due to a longer free-fall time. The resulting protostellar disc is relatively larger and flattened. Figure A.5 is a further zoom-in within 1.5 au at the same time snapshot as Fig. A.4 magnifying the second pseudo-disc and the thermal pressure driven outflow. The vertical feature within the inner disc is already destroyed by the outflow at this time snapshot. We do not pursue a further evolution of the outflow but expect a similar behaviour as in the fiducial case.
<figure><img src="image_27.png" /><figcaption>**Fig. A.2.** 2D view of a first hydrostatic core evolved into a rotationally supported disc at 1427 years after its formation as a result of the collapse of a 1 $M_\odot$ pre-stellar core with an initial rotation rate of $\Omega_0 = 2.099 \times 10^{-13}$ rad s$^{-1}$ (same time snapshot as Fig. A.3). The dust size is fixed to a constant value of 1 µm. The four panels show the **a)** Stokes number, **b)** gas temperature, **c)** gas density, and **d)** radial gas velocity within the inner 15 au of the 3000 au collapsing pre-stellar core. The gas velocity streamlines indicate the material falling onto the disc and the mixing within the meridional flow located at around 10 au.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. A.3.** 2D view of a pseudo-disc forming around the second core at three years after its formation as a result of the same initial conditions and at the same time snapshot as Fig. A.2. The gas velocity streamlines indicate material falling onto the pseudo-disc, circulation within the pseudo-disc via a turbulent flow, and outflowing gas from the vicinity of the unresolved protostar.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Fig. A.4.** 2D view of a first hydrostatic core evolved into a rotationally supported disc at 1646 years after its formation as a result of the collapse of a $M_\odot$ pre-stellar core with an initial rotation rate of $\Omega_0 = 8.16 \times 10^{-14}$ rad s$^{-1}$ or $E_{\text{rot}}/E_{\text{grav}} = 0.007$ (same time snapshot as Fig. A.5). The dust size is fixed to a constant value of 1 µm. The four panels show the **a)** Stokes number, **b)** gas temperature, **c)** gas density, and **d)** radial gas velocity within the inner 15 au of the 5000 au collapsing pre-stellar core. The gas velocity streamlines indicate the material falling onto the disc and the mixing within the meridional flow located at around 10 au.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. A.5.** 2D view of a pseudo-disc forming around the second core at three years after its onset as a result of the same initial conditions and at the same time snapshot as Fig. A.4. The gas velocity streamlines indicate material falling onto the pseudo-disc and outward transport within the pseudo-disc, due to the outflowing gas from the vicinity of the unresolved protostar. At this stage the vortical flow within the pseudo-disc has already been destroyed by the outflow.</figcaption></figure>
**Appendix B: Dust temporal evolution**

In Figure B.1, we display gas temperature at the location of a selected sample of 10 and 100 µm dust for comparison with the 1 µm case previously also shown in Fig. 19. The strong gas-dust coupling within the dense disc leads to a similar temporal evolution for the three different dust sizes. We find regions (see yellow

part of the tracks) where some dust particles could undergo sublimation/evaporation. This will be investigated in our follow-up work. Note that although dust evaporation is not included in our current treatment for Lagrangian particles, it is considered in the opacity tables used in the radiative transfer module.

<figure><img src="image_31.png" /><figcaption>**Fig. B.1.** Tracks for a selected sample of 1 µm (top), 10 µm (middle), and 100 µm (bottom) dust that end up within the two discs during the 1 $M_\odot$ pre-stellar core collapse. The colourmap is split below and above silicate sublimation at 1400 K and indicates the gas temperature at dust location. The vertical lines indicate the formation times of the two discs, and the outflow, and merger of the discs from the onset of the first core formation ($t=0$). The horizontal dashed black line marks the first core disc or merged disc radius defined using the conditions listed in Sect. 3.1. The temporal evolution of the dust provides a clear indication of thermal reprocessing within the very young disc.</figcaption></figure>