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

,


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 Exoplanet 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 environmental 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 Substructures 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 observational evidence for substructures such as rings, gaps, and spiral arms via dust continuum emission from Class II discs (see also Pérez et al. concurrently 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 temperature 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 dynamical 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 formation 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 calciumaluminium-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 magnetic field strength and orientation, angular momentum, turbulence, and dust size distribution in the cloud cores (i.e.prestellar cores).These unknown initial conditions introduce various caveats in our current understanding of the observed stardisc-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 particular molecular hydrogen (H 2 ) dissociation and ionisation; turbulence; 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 optical 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. 2018Bhandare et al. , 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).However, despite the vast amount of theoretical work, several aspects of disc formation and evolution are far from being fully understood, 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.Furthermore, dust grains constitute the main ingredients to form planetary building blocks.Dust properties can also have significant cumulative dynamical effects on the thermal budget and kinematic structure of the disc.Even so, present-day numerical codes using sophisticated techniques rarely account for the role of drag forces and dust decoupling during protostar (second 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 assumption 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 micron.The dust size distribution (number of dust grains per grain radius a) is assumed to be a power law, proportional to a −3.5 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;Pagani et al. 2010;Steinacker et al. 2015;McClure et al. 2023) and evidence of lower dust opacity 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 presence of larger dust grains of up to 1 micron in pre-stellar cores.Dust polarisation detection due to selfscattering processes at (sub)millimetre wavelengths have also provided indirect evidence for up to 100 micron 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 wavelength observations of spatially resolved young discs has been linked to the presence of millimetre-sized dust in the disc midplane (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), andBirnstiel (2023).Modelling transitions of CO emission 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 outflows, winds, and jets (e.g.Bans & Königl 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 interactions 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 framework used in most simulations excludes an accurate treatment of the thermodynamics, which should be modelled by incorporating radiation transport and hydrogen dissociation and ionisation effects, especially important during the second collapse stage (   & 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), Bate (2011), andTomida et al. (2013).
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. 2019Lebreuilly et al. , 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), andCommerç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. (2019Lebreuilly et al. ( , 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 threedimensional 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(Ormel et al. , 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;Dzyurkevich et al. 2017;Zhao et al. 2018;Tsukamoto et al. 2020).Effects of radiation transport and an accurate gas equation of state to treat the influence of H 2 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 micron-sized neutral, spherical dust particles throughout the collapse.The numerical framework and initial set-up is detailed in Sect. 2. In Sect. 3 and Sect. 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.

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.

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;Mignone et al. 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(Kuiper et al. , 2011)).A variable gas equation of state from D'Angelo & Bodenheimer (2013) accounts for effects of H 2 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. 1997Balay et al. , 2019a,b),b).

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 where m dg is the mass of the dust grain and ∆v = v dg − v g is the differential velocity between the dust grain and gas.
Article number, page 3 of 24 A&A proofs: manuscript no.Papercorr The dust stopping time t s is defined as where Γ 1 is the first adiabatic index, c s is the sound speed, ρ dg is the intrinsic dust grain density, s dg is the dust grain size, and ρ 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 θ (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. (2018Mignone et al. ( , 2019)).

Initial fluid set-up
The physical and numerical set-up of the collapse is similar to the one previously described in Bhandare et al. (2018Bhandare et al. ( , 2020)).The initial density distribution has a stable Bonnor-Ebert sphere like profile with the outer ρ o and central ρ c density determined as where G is the gravitational constant and M 0 is the initial prestellar core mass (Ebert 1955;Bonnor 1956).The initial uniform sound speed c s0 is estimated using where R out is the pre-stellar core radius.The density contrast between the inner and outer edge corresponds to a dimensionless radius of ξ = 6.45where ξ is defined as The initial thermal pressure was computed for a fixed temperature T 0 of 10 K.However, the temperature T BE for the stable Bonnor-Ebert sphere is slightly higher than 10 K and is computed as where the mean molecular weight µ = 2.353 g/mol, the adiabatic index γ = 5/3, and ℜ 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 1 M ⊙ and also for an initial 3 M ⊙ 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 ⊙ cases, respectively and computed as where k B is Boltzmann constant and m H is the hydrogen mass.We compare simulations with two different values for Notes.Listed above are the pre-stellar core properties for runs with different initial pre-stellar core mass M 0 (M ⊙ ), outer radius R out (au), rotational-to-gravitational energy ratio E rot /E grav , fixed dust size (µm), and the corresponding discussion section.
the initial solid-body rotation rate of Ω 0 = 1.77 × 10 −13 rad s −1 and Ω 0 = 2.099 × 10 −13 rad s −1 .The ratio of rotational-togravitational energy is approximated as and for the two rotation rates is 0.007 and 0.01, respectively.These are below the disc fragmentation limit of E rot /E grav > 0.01 (Matsumoto & Hanawa 2003;Bate 2011).The fiducial run with the slower rotation rate is compared to a case with a larger prestellar 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(Kuiper et al. , 2011)).The shear viscosity is described using the so-called αparametrisation from Shakura & Sunyaev (1973) and is given by The local sound speed c s ≈ H Ω K (r) with the Keplerian angular velocity Ω K (r) ≈ GM(r)/r 3 .The local pressure scale height H = (H/R)R where the cylindrical radius R = r sin(θ).The dimensionless parameters (H/R) = 0.05 and α = 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 β-viscosity prescription for self-gravitating discs by Duschl et al. (2000) with a βparameter of β = 2.5×10 −3 for the fiducial case at time = 0.This temperature-independent β-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. (2015Kuiper et al. ( , 2016)), and Kuiper & Hosokawa (2018).

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 spherical, neutral dust with an intrinsic density ρ dg of 3 g cm −3 to Article number, page 4 of 24 mimic a combination of carbonaceous and silicate grains as previously used in Tricco et al. (2017) and Lebreuilly et al. (2019).The initial dust density was fixed to one per cent of the gas density using the conservative estimate of the averaged dust-to-gas mass ratio of the ISM.This results in a particle mass distribution 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 particle sampling throughout the pre-stellar core, especially in the innermost regions, we used 200 particles per cell for all the simulations 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 particles 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.Observational 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 distribution 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.

Grid
We used a two-dimensional (2D) axially and midplane symmetric Eulerian grid.The radial grid extends from 10 −2 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 min = ∆r = r∆θ = 5.37 × 10 −4 au.This resulted in 84 cells per Jeans length at the highest central density achieved in the course of the fiducial simulation.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.

Boundary conditions
At the inner radial edge, we adopted reflective boundary conditions 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 temperature with a constant boundary value of 10 K. Additionally, we imposed an outflow-no-inflow condition for the hydrodynamics, which includes a zero-gradient boundary condition for the thermal pressure, the polar, and the azimuthal velocity components.Axisymmetric boundaries were used at the pole and mirror-symmetric boundaries at the equator.

Results I: Features formed during the collapse
Studies over the last ∼55 years 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. 2018Bhandare et al. , 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 ⊙ pre-stellar core collapse The 2D RHD simulation for a 1 M ⊙ pre-stellar core collapse shows an oblate first hydrostatic core transitioning into a rotationally supported disc over a period of 1218 years 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.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 ϕ > 2 |v r |. -Rotational support must be higher than the thermal support (ρv 2 ϕ /2) > 2 P th and v ϕ > 0.6 v K , where v K is the Keplerian velocity.
-Gas density must be above the minimum limit set by the midplane 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 Ω 0 = 1.77 × 10 −13 rad s −1 (or E rot /E grav = 0.007).Comparisons to a faster rotation case (E rot /E 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 years 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 × Ω orbital , where Ω orbital = v ϕ /r sin(θ) 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 years 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 Article number, page 5 of 24  properties at three years after the pseudo-disc 1 begins to form around the second core and is embedded within the first disc.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 Sect.3.1.1and Sect.4.1.

Formation of dust traps in turbulent gas flows
Our results show the formation of two turbulent gas flows during the 1 M ⊙ 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 over different time snapshots.The first core begins to form at 18000 years and the pseudo-disc at 19218 years 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 years.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 ro-Article number, page 7 of 24 tating flow with respect to small isentropic disturbances (Tassoul 2000;Rüdiger et al. 2002).In cylindrical coordinates, this condition requires where P th is the thermal gas pressure and j = R2 Ω is the angular momentum per unit mass.The entropy S is computed 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 pseudodisc spreads horizontally across the pseudo-disc within its short three-year lifetime, only to be destroyed a year later due to the launch of an outflow discussed next in Sect.3.1.2.The shortlived 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(Schulik et al. , 2020)), and planetary spiral wakes in circumplanetary discs (Szulágyi et al. 2022).

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 19222 years after the onset of the 1 M ⊙ pre-stellar core collapse, which is four years after the second collapse and pseudo-disc formation.We follow its evolution for ≈ 55 years 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 polarangle averaged radial velocity of the outflow reaches a maximum of ≈ 9 km s −1 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 ρv 2 ϕ /2 < 2 P th within the outflow can be seen in Fig. 9 at different time snapshots from the onset until 55 years after the outflow is launched.Thermal pressure building up in the pseudo-disc triggers this 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 threedimensional SPH simulations by Bate (2010Bate ( , 2011) ) and 2D gridbased 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 (2010Bate ( , 2011) ) investigated a 1 M ⊙ 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 −1 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 ⊙ 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 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 11590 years 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 years 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.

Effect of initial mass: 3 M ⊙ pre-stellar core collapse
In this section we describe results from a comparison simulation starting with similar initial conditions to those described so far but with a more massive initial pre-stellar core of 3 M ⊙ .The rotation rate fixed to 3.042 ×10 −13 rad s −1 gives the same ϕ /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 1 M ⊙ 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.
Article number, page 9 of 24 Fig. 10.2D view of a first hydrostatic core evolved into a rotationally supported disc at 533 years after its formation as a result of the collapse of a 3 M ⊙ pre-stellar core with an initial rotation rate of Ω 0 = 3.042 × 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.E rot /E grav = 0.007 as for the 1 M ⊙ case. 1 µm dust is used throughout the evolution.The 3 M ⊙ pre-stellar core collapse proceeds alike the 1 M ⊙ 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 years 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 two years into its evolution.The four panels in both these plots show the same gas and dust properties as in Fig. 1 and Fig. 2. The structure and properties of both the discs are    different in comparison to the 1 M ⊙ 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 ⊙ 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 ⊙ 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.

Formation of dust traps in turbulent gas flows
In the case of a 3 M ⊙ 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 years after the collapse is initiated.This is followed by the pseudo-disc building up around the second core at 9231 years, which evolves for a few years.The two nested discs merge by 9252 years.A negative entropy gradient generates turbulent gas flows within the two discs.In comparison with the 1 M ⊙ 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 Article number, page 11 of 24 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.

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 years, reaches a maximum of ≈ 6 km s −1 , 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 pseudodisc 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 years 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.

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

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 where Stokes number St is previously defined in Sect.3.1 and Ω K = GM(r)/(r sin θ) 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 gas density within the disc, all three sizes have a Stokes number Fig. 15.Rotational support (ρv 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 ⊙ pre-stellar core collapse.The topmost 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 streamlines show the infalling envelope, mixing within the discs, and the outflowing gas.
Article number, page 12 of 24 Fig. 16.Ratio of dust settling to eddy turnover time only within the first core rotationally supported disc formed during the fiducial 1 M ⊙ 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.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 ρv 2 ϕ /2 > 2 P th within the disc can be seen in Fig. 17 at different time snapshots from the start of the outflow until 55 years 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 effect over the disc surface.Once the outflow is launched dust particles from the innermost region are ϕ /2 > 2 P th ) at different time snapshots during the 1 M ⊙ 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 wellcoupled dust follows the gas motion in all the different regions.
Article number, page 13 of 24 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 ⊙ 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.Article number, page 14 of 24 Fig. 19.Tracks for a selected sample of 1 µm dust that end up within the two discs during the 1 M ⊙ 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.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 ≥ 10 µm dust.
The Lagrangian particle treatment allows us to trace the thermal history of different dust grains.Figure 18 and Fig. 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).

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.

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 follow the evolution of the protostar-disc-outflow system for 85000 years 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 ⊙ pre-stellar core collapse they also find ≤ 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 ≥ 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 mo-Article number, page 16 of 24 mentum 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-togas 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-togas 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.

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 2 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 substructures and instabilities through disc evolution stages and for a dust size distribution need to be investigated.

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 cosmochemical 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 ⊙ 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 ⊙ 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 ⊙ 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 ⊙ 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 Article number, page 17 of 24 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.

Fig. 1 .
Fig.1.2D view of the first hydrostatic core evolved into a rotationally supported disc at 1221 years after its formation as a result of the fiducial collapse of a 1 M ⊙ pre-stellar core with an initial rotation rate of Ω 0 = 1.77 × 10 −13 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.

Fig. 2 .
Fig. 2. 2D view of a pseudo-disc around the second core at three years 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 hightemperature inner regions.

Fig. 3 .
Fig. 3. Polar-angle averaged radial profile of the entropy at different time snapshots during the 1 M ⊙ pre-stellar core collapse.A negative radial gradient generates turbulent flows within the two discs.

Fig. 4 .
Fig. 4. Thermal pressure P th in colour and isodensity contours within the first core disc at 1221 years.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.

Fig. 5 .
Fig. 5. Baroclinicity (s −2 ) at 1221 years after the transition of the first core into a rotationally supported disc resulting from a 1 M ⊙ 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.

Fig. 6 .
Fig. 6.Solberg-Høiland unstable region within the first core disc at 1221 years found mainly around the meridional feature.

Fig. 7 .
Fig. 7. Polar-angle averaged radial velocity at different time snapshots showing the transient nature of the outflow during the 1 M ⊙ pre-stellar core collapse.

Fig. 8 .
Fig. 8. Radial Mach number (v r /c s ) at first and last time snapshots of outflow formation during the 1 M ⊙ 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.

Fig. 9 .
Fig. 9. Rotational support (ρv 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 1 M ⊙ 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.

Fig. 11 .
Fig. 11.2D view of a pseudo-disc around the second core at two years 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.

Fig. 12 .
Fig. 12. Polar-angle averaged radial profile of the entropy at different time snapshots during the 3 M ⊙ pre-stellar core collapse.A negative radial gradient generates turbulent flows within the two discs.

Fig. 13 .
Fig. 13.Polar-angle averaged radial velocity at different time snapshots showing the transient nature of the outflow during the 3 M ⊙ pre-stellar core collapse.

Fig. 14 .
Fig. 14.Radial Mach number (v r /c s ) at first and last time snapshots of outflow formation during the 3 M ⊙ 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.

Fig. 17 .
Fig. 17.Rotational support (ρv 2ϕ /2 > 2 P th ) at different time snapshots during the 1 M ⊙ 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 wellcoupled dust follows the gas motion in all the different regions.

Fig. 20 .Fig. 21 .
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 ⊙ 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.
Article number, page 2 of 24 A.Bhandare et al.:Gas & dust dynamics during pre-stellar core collapse