Collapsing molecular clouds with tracer particles: Part II, Collapse Histories

In order to develop a complete theory of star formation, one essentially needs to know two things: what collapses, and how long it takes. This is the second paper in a series, where we query how long a parcel of gas takes to collapse and the process it undergoes. We embed pseudo-Lagrangian tracer particles in simulations of collapsing molecular clouds, identify the particles that end in dense knots, and then examine the collapse history of the gas. We find a nearly universal behavior of cruise-then-collapse, wherein a core stays at intermediate densities for a significant fraction of its life before finally collapsing. We identify time immediately before each core collapses, $t_{\rm{sing}}$, and examine how it transitions to high density. We find that the time to collapse is uniformly distributed between $0.25 t_{\rm{ff}}$ and the end of the simulation at $\sim 1 t_{\rm{ff}}$, and that the duration of collapse is universally short, $\Delta t \sim 0.1 t_{\rm{ff}}$, where $t_{\rm{ff}}$ is the free-fall time at the mean density. We describe the collapse in three stages; collection, hardening, and singularity. Collection sweeps low density gas into moderate density. Hardening brings kinetic and gravitational energies into quasi-equipartition. Singularity is the free-fall collapse, forming an envelope in rough energy balance and central over density in $\sim 0.1 t_{\rm{ff}}$.


INTRODUCTION
Stars form from clouds of hydrogen that are cold, turbulent, and magnetized.Temperatures tend to be around 10 K, densities around 100 cm −3 , turbulent velocities a few times the sound speed, and magnetic fields are around a few microgauss (Solomon et al. 1987;Crutcher 2012;Miville-Deschênes et al. 2017).Prestellar cores are denser and more magnetized (Crutcher 2012) objects that are often associated with embedded infrared sources.It appears that a population of these cores are coherent, meaning sub-or trans-sonic (Goodman et al. 1998;Singh et al. 2021;Choudhury et al. 2021).The mass distributions of cores (Ladjelate et al. 2020) and the initial mass function of stars (Chabrier 2003) are similar in shape but offset by a factor of  = 3 − 4.This indicates that cores are related to the formation of stars.One interpretation is that a fraction of each core, 1/, given to each star (Alves et al. 2007), though this may be coincidental (Clark et al. 2007).
The simplest model for collapsing gas is pressure free collapse in one dimension.The density behavior can be well approximated as where  ff = (3/32  0 ) 1/2 and  = 1.86 (Carroll & Ostlie 2006;Girichidis et al. 2014).This model neglects thermal pressure, magnetic fields, kinetic energy, and chaotic motions in three dimensions, and it produces infinite densities, which is not physical.Nevertheless is will be useful in describing the collapse.
★ Contact e-mail: dccollins@fsu.edu To improve on this model, Larson (1969) and Penston (1969) simulated a collapsing sphere as a radiating hydrodynamic system in one dimension.They find that the collapse halts when the first hydrostatic core is formed, which happens when the density and opacity are high enough to trap photons in the core, causing a dramatic increase in the internal pressure.These objects are typically 1-20 AU, 0.05 ⊙ , and have number densities of  = 10 11 cm −3 (Young 2023).Larson (1969) also showed that the collapse tends to leave an envelope that has a profile of  ∝  −2 .
The behavior in one dimension is dependent on its initial conditions.Shu (1977) showed that a sphere will collapse in a self similar manner, and an initial density profile  ∝  −2 will collapse from the inside out.Whitworth & Summers (1985) generalized this self similar behavior, and showed that inside out collapse was a result of the initial conditions.These models still neglected the turbulence that has come to dominate our understanding of star formation.Early models had magnetic fields regulating the collapse by slowing the gravitational contraction (e.g.Shu et al. 1987;Mouschovias 1987), but these models were largely ruled out by observations (Crutcher et al. 2009).
The importance of turbulence was first indicated by Larson (1981).Here it was found that velocity and size of molecular clouds are correlated where   = 0.38 and   = −1.1 were found in that work.This has been interpreted to be the result of turbulence in the interstellar medium (eg Mac Low & Klessen 2004;Elmegreen & Scalo 2004).
Another interpretation of the observations puts gravitational collapse (Ballesteros-Paredes et al. 2011) as the primary source of the relationship.Further, the consistency of these relations in the Galaxy have been called into question; Heyer et al. (2009) show that the relation depends on the clouds' column densities.Nevertheless, these relations indicate complex physics at the cloud level that will dictate the cloud's collapse (???).
Many modern turbulent fragmentation models are based on Padoan (1995), wherein some subset of the lognormal density PDF is turned into stars (Padoan & Nordlund 2002;Krumholz & McKee 2005;Hennebelle & Chabrier 2011;Padoan & Nordlund 2011), though the assumptions behind these models are different from one another.More recent models include the inertial flow model (Pelkonen et al. 2021), where large scale flow drives core formation, and the three phase model of Offner et al. (2022), wherein turbulence gives way to coherent structures.Our simulation set up and outcomes are similar to these, with a few notable differences due to the novel way in which we examine the collapse.This is the second paper in a series wherein we inject pseudo-Lagrangian tracer particles in a simulation of turbulent, collapsing molecular clouds.This is an extension of Collins et al. (2023, hereafter Paper I).We allow the cloud to collapse to form several hundred "stars", identify the particles in those "stars", and then examine the cells that those particles occupy during the collapse.At all times prior to the end of the simulation, we refer to zones that contain core particles as preimage gas.In Paper I, we examine the preimage gas at  = 0, the beginning of the simulation.We find that the density PDF is fully covered by the PDF of the preimage gas, but suppressed by a factor of roughly (/ max ) 1/2 : all gas participates in the collapse, dense gas preferentially.We examine length scales of both density and velocity of preimage gas at t=0.We find that the length scale of the preimage gas is several times larger than the density auto correlation length, as many density fluctuations feed a single star.We find that the length scale of the preimage gas is comparable to the velocity correlation length, indicating one or a few large scale flows dictate the collapse.We find that gas is distributed in a fractal manner, with the Minkowski dimension being distributed between 0.25 and 2, peaking at a dimension of 1.6 (see Section 3.3.4for definitions).We also find a high degree of overlap of the preimage gas between different cores.Cores that are close to one another at the end of the simulation, e.g.binaries, come from gas that seems to be mixed in the initial cloud.Many density fluctuations feed each core, and many cores come from each collection of density fluctuations.
In this work, we examine the time history of collapsing gas.We use the same simulations described in Collins et al. (2023), but now follow them beyond  = 0.Here we focus on the appearance of self similar collapse from turbulent initial conditions.We limit our study to cores that collapse in isolation.In future installations of the series we will discuss magnetic fields, binaries, and clusters.
The paper is organized as follows.In Section 2 we discuss the code, simulations, and particle identification.Results are presented in Section 3. In Section 3.1 we discuss the collapse time,  sing , its distribution, and the unexpected lack of correlation between  sing and mean quantities of the preimage gas.Section 3.2 follows in detail the collapse of one fiducial object, and we discuss the three stages of collapse (collection, hardening, and singularity).We then show the behavior of the ensemble of isolated objects in Section 3.3.We summarize and discuss in Section 4.

METHOD: CODE AND CORES
This work uses the simulations and methods described in detail in Collins et al. (2023, Paper I).We will briefly recap the code, simulations, and particle methods here.

Code
We use the open source adaptive mesh refinement code Enzo (Bryan et al. 2014) with the constrained transport magnetohydrodynamics module (Collins et al. 2010).This code adaptively adds resolution as the collapse unfolds, when the local Jeans length  J =  2  /( ) 1/2 exceeds 16 zones.The main solver is a higher order Godunov method, using the linear method of Li et al. (2008) for the differential equations, the HLLD method of Mignone (2007), and the electric field defined by Gardiner & Stone (2005) to maintain the divergence of the magnetic field.The code adaptively adds zones using the methods of Berger & Colella (1989) and Balsara (2001).It has been demonstrated to preserve the divergence of the magnetic field to machine precision (Collins et al. 2010), and has been used in a wide array of applications from star formation (Collins et al. 2011) to galaxy clusters (Xu et al. 2012;Skillman et al. 2013).We additionally include tracer particles that follow the flow by way of a piecewise linear interpolation of the velocity field (see Bryan et al. 2014;Collins et al. 2023, for a description).Details of particle analysis can be found in Section 2.4.

Simulations
Our simulations begin with periodic turbulent boxes that have been stirred in a manner described by Mac Low (1999).Here kinetic energy is added to the large scales of the box, and the nonlinear dynamics carries that energy to small scales.Driving continues until a steady state is reached, at which time gravity, refinement, and tracer particles are added.We will describe the simulations and then the physical dimensions of the simulations in the next two sections.

Simulation Layout
Ideal magnetohydrodynamics is scale free, and can be characterized by dimensionless parameters.To characterize our three simulations, we employ the Mach number M, virial paramter  vir , and the ratio of thermal to magnetic energy,  0 : Here  rms is the r.m.s.velocity,   is the sound speed,  0 is the mean density,  0 is the box size, and  0 is the mean field, taken to be uniform in the  ˆdirection.Simulations continue for slightly less than one free-fall time of the mean density,  ff = (3/32  0 ) 1/2 .The simulations begin with uniform density and magnetic field and are driven until steady state is reached.This was done for three magnetic fields at a resolution of 1024 3 .Once steady state was reached, we downsampled to 128 3 , turned on gravity and refinement, and added the tracer particles (one per zone).Refinement was triggered whenever the local Jeans length became smaller than 16 zones,   =   /( ) 1/2 < 16Δ.This continued for 4 levels for a maximum effective resolution of 2048 3 .

Physical Parameters
While we can formally select the physical units to be what we like, subject to the above dimensionless ratios, it is useful to select a fiducial scale.To do so we assume that our simulations are consistent with Larson's relations (Larson 1981), see Equations 2 and 3. We must also adopt a sound speed.We assume  = 12K, which gives a sound speed of  s = 0.2 km s −1 .
Thus the physical scale for velocity, length, density, mass, global free fall time, and smallest cell size are found to be:

Pseudo-cores and Sink Particles
The first hydrostatic core forms when the opacity of the gas grows large enough to stiffen the equation of state from isothermal (Larson 1969).These objects are only ∼ 20 AU across, and as such not resolvable in our simulations.An approximation must be made.Most often sink particles are employed, which remove some mass from the grid and replace it with a particle that has zero volume (Federrath et al. 2010;Teyssier & Commerçon 2019).This is the industry standard for dealing with the large densities.
We have chosen to not use sink particles for two reasons.First, the creation of the sink disrupts the very process we are trying to measure, since the algorithm will remove mass, kinetic, and thermal energy from the collapsing object.Second, there are presently no sink algorithms that properly treat magnetic fields, each opting to ignore (Teyssier & Commerçon 2019).This violates Alfven's theorem that we work so hard to enforce, and the consequences on the collapsing dynamics have yet to be studied.
Instead we allow the grid pressure to supply the additional support that would be provided by thermal pressure.This allows us to examine the onset of self similar collapse and examine the precipitating prestellar core without the interference of other algorithms.This is philosophically similar to the idea of implicit large eddy simulations (Adams & Hickel 2009).In an such a simulation, the grid viscosity is used in place of an explicit sub grid model of the unresolved scales.The basic idea is that the nature of the turbulent cascade does not depend on the nature of the dissipation, only its existence.In a similar manner, the collapse of a collapsing molecular cloud will not depend on the nature of the approximation we make at the end, within reason.We must then use care in dealing with the numerical artifacts that are leftover, regardless of their nature.
At the end of the collapse, we form small (a few zones across, ∼ 1000 AU) high density ( ∼ 10 5  0 ∼ 10 8 cm −3 ) condensations that we refer to as pseudo-cores.They're several hundred times larger than the first core they aim to replicate and thus not numerically resolved.Therefore we refrain from quantitative analysis of the population of objects, and we cease detailed analysis of each after they have formed.For this reason, among others, we restrict our analysis to isolated cores that do not interact, as it is clear when analysis should be stopped for each core.
To identify pseudo-cores, we use a dendrogram type technique (Rosolowsky et al. 2008;Turk et al. 2011) to find the location of dense peaks.We keep only the peaks above 10 4  0 , as most dense peaks have  peak ∼ 10 5−7  0 , while turbulence only makes peaks with  < 10 3  0 .There are few fluctuations with  ∼ 10 4  0 , these are either turbulent fluctuations or too early in their collapse to learn anything from.For each peak, we take the particles closest to the densest zone.For practicality, we take all the zones such that  >  3/4 peak around each peak, but the definition of which zones to take does not matter substantially as the particles cluster in the densest few zones.We refer to the pseudo core and the particles within collectively as a core, as we are interested in the prestellar core around each pseudo-core, and in this analysis each one is a unique path in space.We refer to zones that are flagged as containing core particles at earlier times as preimage zones.
In our analysis we focus on two things; the process that delivers mass to the pseudo-core, and the envelope around the core.The ∼ 0.1pc envelope around the pseudo-core would potentially be considered a prestellar core, but we do not tie ourselves to these observationally defined objects.
It is worth emphasizing that a pseudo-core is not a resolved object, we introduce them here only to isolate the unresolved endpoint of the collapse from the collapse itself and the remaining envelope in our presentation of the results.

Particles
Tracer particles are evolved using a basic kick-drift approach, where the particle velocity v p is interpolated from the grid velocity in a linear manner, and the position is then updated as x +1 = x  + Δ v p .Particles are initially deposited uniformly throughout the box, one per zone.
Particles have a density (defined by the inter particle spacing) and velocity (interpolated from the grid) that are not of interest.We only use the particles to determine which gas zones to analyze.One of two methods is used: either the zones containing the particles (i.e.preimage gas), or a sphere defined by the centroid and extent of the particles is used.The sphere has a minimum radius of one root grid zone, as the particles eventually occupy only a few fine-grid zones.
Particles, unlike gas, do not feel each others' presence or pressure of any kind.One result of this is over clustering of the particles relative to the gas, which results in a much wider distribution of particle densities (Genel et al. 2013).The reason for this is that two particles, once within the same zone in a converging flow, will be driven to zero separation in a short time (see Collins et al. 2023).Particle velocity is linearly interpolated from the grid velocity, and as such it is unlikely to separate the particles as the nonlinear turbulent field would be were the subgrid velocity resolvable.This results in a loss of the number of degrees of freedom the particle flock has as the collapse proceeds.For our simulations, this is a useful feature rather than a problem.It means that the particles do not inject any noise of their own, so we can reliably trace them back to their origin.We take care to not perform any measurements that rely on the information in the particles after they collapse and all occupy the same few zones.
The net effect of the over clustering is in fact a useful feature of the particles for this study.Once the gas has gotten dense and many particles occupy a single zone, their individual contributions no longer have useful information.By sticking together they effectively remove themselves from consideration as individual pieces of information, eliminating confusion they may cause.
On remedy for the over clustering that has been suggested is to add some effective kinetic temperature to the particles, which would allow the particles to entrain with the turbulence more effectively (Genel et al. 2013).For this study, adding noise to the particle trajectories would destroy the causal connection between the particle at the beginning and the end of the simulation, and is not desirable.

Mode Separation
We begin our analysis by classifying the collapse mode of each core; single, binary, or cluster.This allows us to avoid the effects of unresolved objects and tidal interactions.
For each core, we measure the distance, , between its centroid and the centroid of all other cores for all time steps.The number of cores within  = 0.05 0 ,   , indicates the mode of collapse; single (  = 1), binary (  = 2) or cluster (  > 2).The value of  was determined by examining the path of the centroid for each core and its neighbors.The results presented here are insensitive to the exact value of .This simple grading scheme was verified by visual inspection of movies of the collapsing cores.While more sophisticated methods for identifying neighbors can be concocted, this simple scheme is sufficient for our purposes here.Table 1 shows the number of cores in each mode for each simulation, though these demographics are as approximate as our separation scheme.

RESULTS
Now we describe the behavior of the collapse.In Section 3.1, we define and measure  sing and  sung , which indicate the beginning and end of the free-fall collapse.We also measure the distribution of collapse times, showing a universally short collapse, faster than free fall of the mean cloud.We focus our attention on a careful examination of one representative core in Section 3.2.We then show this same behavior for all single cores in Section 3.3.

The singularity time, 𝑡 sing
Figure 1 shows density vs. time for two cores.Core 214 is quite slow, not collapsing until near the end of the simulation, while core 74 begins to get dense much earlier, with a few particles reaching high density early, and the remaining particles accreting over the  remainder of the simulation.It is useful to identify the moment immediately before the collapse, which we denote  sing , and the moment when the collapse ends, which we denote  sung .These are found by examining when the time derivative of the maximum of density over the preimage gas;  sing is defined as the time when the derivative gets above a threshold value, and  sung is defined as the time when the density derivative is zero again.That is, where the threshold value is found to be 10 5 in code units, or 7 × 10 7 cm −3 Myr −1 using the scaling defined earlier.These are indicated by vertical lines in Figure 1.
The cumulative distribution of  sing for all cores can be seen in the left column of Figure 2.Each row shows a different simulation.Red lines show single cores, green shows binaries, and blue shows clusters.For each population, the cumulative distribution indicates that the distribution of  sing is relatively uniform between 0.25 ff and 1 ff when the simulation ends.Some fluctuations exist, but no consistent trends are noticeable.
The right column of Figure 2 shows the distribution of the duration of the initial singularity, Δ sing =  sung − sing .Note that the horizontal axes between the left and right columns are not the same.The collapse shows a nearly universal short duration.
We can explain the typical behavior of Δ sing by noting the almost free-fall behavior at the end of the collapse.The typical density of the preimage gas at  sing is 500 0 , and the grid pressure sets in at roughly 10 5  0 .Free fall is self similar, so the timescale to go between two densities is proportional to the square root of the densities.This gives a collapse time of ∼ 0.1 ff to go from the density at  sing tolthat at  sung .This does not, though, imply a critical threshold for collapse; by the time the gas reaches the density it has at  sing the collapse has already been triggered.

What sets 𝑡 sing ?
It would be extremely useful to understand what made some cores collapse faster than others.Unfortunately, after examining many possible correlations, we do not find any that show strong correlations with the collapse time.For each simulation, we examine Pearson  for  sing vs. several properties of the cores.Pearson  measures the degree of linear correlation between two quantities, with -1 being perfectly anti-correlated.
is defined as Here the average ⟨⟩  is taken over all cores, and  is the average or total of the quantity over the preimage gas for each core.For the extrinsic quantities,  represents the total over preimage, which includes the number of particles,  particles , and the volume of the convex hull bounding the preimage,  hull .For other quantities,  is the mass-weighted average of the quantity over the preimage gas for core .
We examine the correlation coefficient for the following quantities: the number of particles  particles ; the volume of the convex hull bounding the particles  hull ; the mean density, ; the mean gravitational binding energy,  G ; the mean kinetic energy,   ; the r.m.s.velocity,  rms ; the mean radial velocity,  R ; the mean tangential velocity  T ; and the free-fall-time at the mean density of the core, Table 2 shows  for these quantities for each of the three simulations.Only the single cores are included in this table.There is only one correlation that is as high (in magnitude) as 0.5,  for 3, The number of particles has a mild trend for anti-correlation (-0.41, -0.37, and -0.24) for the three simulations, respectively.The local free fall time shows a mild positive correlation, with coefficients of (0.37, 0.38, and 0.48), respectively.The velocity and energies are all weakly correlated, with no consistent trend in the simulations.
While the trends are in the direction one would expect (e.g. sing is positively correlated with  ff , and the velocity dispersion is negatively correlated) the scatter is too large to be predictive.This is not to say that  sing is not predictable, but the obvious things to check are individually not able to predict  sing .

Collapse Paradigm
Here we examine one core in detail.To further understand the processes at work, we would like to examine   () for every particle, , in every core, as well as velocity, kinetic, and gravitational energies.This is several thousand plots.We summarize these plots by way of a case study (this section), followed by an ensemble of cores (next section).

Three stages
We find that each core goes through three stages of collapse to one degree or another.We refer to these as collection, hardening, and singularity.The collection stage sweeps low density gas together by converging supersonic flows; kinetic energy is typically larger than gravitational, but not universally.This stage does not appear to happen for every core.During the hardening stage, the density becomes moderately high (∼ 10 0 ), the velocity decays to trans-to slightly super-sonic, and gravitational and kinetic energies equilibrate.During the singularity, density rises extremely fast and the majority of the mass is delivered to the core in ∼ 0.1 ff ; gravitational energy dominates.By the end of the singularity, the envelope has a density profile proportional to  −2 , but not generally spherical, and the ratio of kinetic to gravitational energy follows a remarkably consistent pattern dropping from  K / G ∼ 4 at the inner radius to  K / G = 0.5 − 2 at outer radii.

Case Study
In Figure 3 we show the collapse history of one representative core (core 114 from 2) outlining the three stages of collapse.This plot summarizes how density, velocity, gravitational energy, and kinetic energy change in concert with one another.The top row shows the preimage density, one line for each particle (grey lines, right axis) and the average velocity (r.m.s in black, radial velocity in red, and tangential velocity in teal, left axis).Also indicated are six snapshots in time (vertical bars).The second row shows the cumulative density distribution for preimage density at each of those snapshots.The third row shows radially averaged kinetic and gravitational energy for a sphere surrounding the particles at each snapshot.The first two rows utilize only the preimage gas, while the third row measures  G and  K on spheres centered on the centroid of the particles.We will elaborate on each of these in turn.
Six snapshots in time, indicated by the vertical bars in the top row of Figure 3, highlight the stages of the collapse; the initial conditions (navy), the collection stage (light blue), the hardening stage (teal and light green), the beginning of singularity ( sing , dark green), and the end of singularity ( sung , red).We do not presently have a quantitative boundary between collection and hardening, as no consistent trends are obvious (see Section 3.3).
The top row in Figure 3 also shows the density and velocity behavior of this core.Density begins low, sampling the turbulent density field, as discussed in Paper I. The total velocity is initially quite large,    = 6  in this core (not universal).During the collection stage (blue lines), the large scale flow collects the material into a smaller, higher density patch.During the hardening stage (teal, light green lines), the density distribution becomes more narrow and the velocity decays.For this core, the radial velocity becomes transsonic by  = 0.3 ff , at which time the tangential velocity increases over the radial.The r.m.s velocity is subsonic at the start of singularity, but this is also not a universal feature among all cores.After the end of the singularity, the core accrets the remaining particles.As we discuss in Section 3.2.3,this stage does not accrete much mass, the vast majority of the mass is assembled during the singularity.
The second row in Figure 3 shows the cumulative preimage density distribution at each of the six snapshots.At the first time (navy line), the density is roughly lognormal, as it samples the turbulent initial conditions.This core begins in a large converging flow which sweeps most (80%) of the particles to moderate density ( = 10 − 100 0 ), while some remains lower ( = 10 0 ).During the hardening stage, the moderate-density gas gets moderately more dense, and the low density gas is compressed, but the peak density does not increase dramatically.During the singularity (between the dark green and red lines), the peak density increases rapidly.Roughly 20% of the particles are above a density of 10 7 cm −3 at the end of the collapse.
The third row in Figure 3 shows the gravitational and kinetic energies of the gas on spheres surrounding the particles.The dotted line shows kinetic energy, averaged over the sphere of volume  (), and the solid line shows the magnitude of gravitational energy: where  () is the volume of a sphere with radius .Initially the gravitational energy is flat and kinetic is several times larger.By the end of the hardening stage (light green line) the profiles of  G and  K have been driven toward each other.Singularity drives an enormous increase in both kinetic and gravitational energies, and we're left with a relaxed envelope with a rapidly rotating pseudo-core in the center.

Mass Delivery
We now examine the nature of the gas mass during the collapse.Figure 3 only shows density of the preimage gas, which covers a sparse subset of the surrounding material.It is valuable to also make a complete census of the gas around the collapsing object.
Figure 4 shows mass vs. time and radius.At each time, a sphere of gas is extracted that is centered on the centroid of the particles and has a radius equal to the most distant particle.The particles all collapse to a few zones at the end, so the analysis sphere has a minimum radius of one root-grid-zone.The vertical lines show  sing and  sung .
The top panel shows interior mass vs time, The edge of a clump is poorly defined, so to avoid murky clump edge definitions we adopt a fiducial radius of 1000AU and fiducial mass,  1000 , the mass within that radius at the final snapshot.The colorbar is arranged so that the white line follows a mass  1000 .It should be noted that this is not, strictly speaking, a Lagrangian surface, since there is substantial flux both in and out of this surface.It does show that the mass is assembled at the center of the core rapidly during the singularity, and then remains relatively constant.The middle panel shows the normalized relative accretion rate: This is computed by way of finite differences of the binned  < in the top plot.We see a large positive accretion rate during the collection phase ( < 0.2 ff ), followed by a high degree of volatility, fluctuating around zero.Beginning at  sing and ending at  sung ,  grows to a factor of several as most of a solar mass is delivered to the center in 0.1  ff .Then the high frequency rotational motions of the pseudo-core make the mass fluctuate in its location, but not gain much overall mass.The lack of growth in mass is due to the fact that the envelope is virialized, with large kinetic energy.The third panel of Figure 4 shows the ratio of kinetic to gravitational energy,  K ()/| G ()|, both defined in Equation 17.The colorbar is tuned so that  K ()/| G ()| = 1 is white; gravitationally dominated gas is in blue, and kinetic energy dominated gas is in red.It is seen that at the beginning of the simulation, kinetic energy dominates the bulk of the flow.As time progresses during the hardening phase, the ratio grows to roughly unity.During the singularity, gravity dominates the entire system, and immediately after  sung we find a pattern slightly dominated by kinetic energy ( K /| G | = 2), and transitions to  K /| G | = 0.5 in the outer regions.
It should be noted that these dynamics are visible to us because we did not use sink particles.Sink particles (Teyssier & Commerçon 2019) are a useful tool to follow self-gravitating simulations in time, but they are a subgrid model selected by the simulator, and the accretion rate they give is, while generally reasonably constructed, dramatically influenced by the parameters of the sink.Here we allow the dynamics to unfold self-consistently, at the cost of shortened simulation time.This allows us to measure how singular collapse from turbulent initial conditions unfolds.

Ensemble properties
The previous section examined one core in detail.Now we examine a larger sample of cores to see how many of these traits are universal, and how many have some variation.We restrict our examination here to single cores, from the first simulation 1.Binary cores, and cores from the other simulations, follow similar trends.A fraction of cluster cores do as well but many cluster cores have much more complex velocity signatures and are likely not actually relaxed at the end of the simulation.We revisit binaries and clusters in a future work.

Mean density and velocity of preimage gas
Figure 5 shows mean of preimage gas quantities for every core in the sample.Each core has been scaled in time to its own  sing .The top row shows mean density, , which shows a universal behavior of slow increase followed by a transition to a free-fall solution (cruise-thencollapse).The second row shows the r.m.s.velocity of the preimage gas,  rms .The flow begins highly supersonic in total.Velocity decays by  sing to a Mach number of a few.The third row shows the mean radial velocity, | R |.Interestingly, there is not a universal trend of large converging flow.In fact, some cores have very small radial velocities initially.To the extent that the collection stage is real, this can be seen in this plot as large radial velocities that decay to transsonic to slightly supersonic.Finally, note that this is the absolute value of the radial velocity, but all mean velocities are negative.By  sing , the radial velocity is clustered around transsonic.Some are subsonic, as expected by others, while some are slightly supersonic.The fourth row shows the tangential velocity, | T |, which is universally supersonic initially, and decays to transsonic before  sing .

Mass Delivery
The fast mass delivery for each core can be seen in Figure 6.Here we define  1000 as the radius that contains a mass of  1000 , which is defined as the mass within 1000 AU at the end of the simulation (see Section 3.2.3.)This can be seen as the white line in the top panel of Figure 4 for one core.In Figure 6, we show  1000 for each core vs. time, scaled to the cores own  sing .Most cores deliver the bulk of the mass of the final object beginning at  sing and taking about 0.1  ff to get there.A few behave slightly differently, as these are ultimately chaotic dynamics.Note that this plot measures all gas, not just preimage gas.

Scaled Free Fall
We find that density is reasonably described by scaled free fall.Given an isolated pressure free sphere of gas collapsing under its own gravity, it's radius can be computed as a function of time.This yields a transcendental equation for radius that can be approximated by where  = 1.8614 (Girichidis et al. 2014).We find that the density of a collapsing object is reasonably approximated by scaled free fall, This is faster than free fall, as  sung <  ff for all cores in our study.
Figure 7 shows scaled free fall and the renormalized density for each core.The top panel shows the maximum of the density for each core, normalized to its value at the start, and the time is stretched to the end of the singularity,  sung : This is plotted for each core, along with Equation 22. Scaled free fall forms a lower bound for the suite of curves; clearly free fall is not the only physics at work, but it serves as the stage on which the dynamics unfold.
The lower panel of Figure 7 shows a similar plot, but with the maximum over density replaced with the mean of the preimage density.The mean densities do not decrease past  = 0 the way the maximum density does, and the curves are all somewhat denser that what free fall at this rate would collapse.Again, free fall collapse is not the only physics at work, but is a useful framework.
If a useful prediction of  sung could be found, this could be used in a predictive theory of star formation.

Fractal dimension
In Paper I, we examined the fractal dimension of the preimage gas at the initial snapshot.We use the Minkowski (box counting) dimension, defined as () is the number of boxes of size  that are needed to cover the preimage gas as  decreases.We find in the first paper that the distribution of dimensions for the preimage gas spreads between  = 0.25 and  = 2, peaked at 1.6.Figure 8 shows  vs. time for each of the single cores from 1, normalized to the collapse time  sung .Here, a grid at the coarsest resolution ( = 1/128 in code units) was used to cover each grid, and  is computed by increasing .As the collapse proceeds,  decreases nearly monotonically for every core.This is not surprising, as gravity tends to reduce the dimensionality of a collapsing object: a triaxial object will first collapse along the shortest axis forming a pancake, and then along the next shortest axis to form a filament, and then along the filament to form a knot (Zel'dovich 1970).This is mirrored in the fractal dimension behavior, as the dimension decreases.

Radial Profiles at 𝑡 sing and 𝑡 sung
Figure 9 shows azimuthally averaged quantities for spheres around every core (grey lines) at four different times (columns).The four columns show  = 0, 0.5 sing ,  sing , and  sung .Note that  sing and  sung are different for every core, see Section 3.1 for definitions and distributions.The rows show azimuthal averages vs. radius for density, radial velocity, tangential velocity, and ratio of energies.Horizontal lines mark relevant quantities; radial velocities of 0 and -1 and energy ratios of 0.5 and 2. The vertical line in the last column shows the approzimate radius of the pseudo-core, inside of which the result is dominated by numerics.
The first row shows azimuthally averaged density, for every core at four frames.Here  () is a sphere centered on the centroid of the particles with radius , and  () is its surface.At  = 0, density is roughly flat, showing a mild tendency for increased density at the center.As discussed in Collins et al. (2023), this is a result of the density correlation length being significantly shorter than the size of the preimage gas.By 0.5 sing , a centrally condensed object is beginning to form, some of them have developed an outer region that follows a  ∼  −2 behavior.By  sing , the central densities have grown somewhat and all objects begin to resemble Bonner-Ebert spheres (Ebert 1955;Bonnor 1956).In the small window between  sing and  sung , free-fall collapse sets in and the remaining envelope is roughly  ∼  −2 as predicted by Penston (1969).The red line in the last panel shows  −2 for comparison.
The second row of Figure 9 shows mean interior radial velocity, This is the average radial motion within a radius , relative to the central velocity for the core, v  , along the radial unit vector for that core,  ˆ.At  = 0, the radial velocity does not have a typical trend, some contracting and some weakly expanding.By the end of the collection phase at 0.5 sing , the radial velocity is negative for all particles.Note that at early phases, the curve of   () does not, in general, go to zero as  → 0. This is due to the fractal and chaotic nature of the preimage gas; the geometric center, the center of mass, and the center of velocity are not at the same point in space.By  sing , the radial velocity tends to be logarithmic in radius, and subsonic.We posit that this phase is what is observed as a coherent core.Between  sing and  sung , the majority of the mass is delivered to the center by negative radial velocity.After the relaxation event, at  sung , the radial velocity behavior is now increasing with radius, and generally  Radial profiles for several quantities (rows) at several frames (columns) for each single core (grey lines).Four times at t=0,  = 0.5 sing ,  sing , and  sung are shown.Note that  sing and  sung are different for each core.The rows show mean density, radial velocity, tangential velocity, and ratio of gravitational to kinetic energy.Notable features are highlighted by red lines, see text for details.supersonic: This reversal happens by way of a wave that travels inwards.The wave follows the knee in the density profile.The red lines in the third and fourth panel show the approximate scaling indicated by Equation 28.
The third row of Figure 9 shows tangential velocity, which is just the full velocity less the radial in the frame of the core.Also plotted in the third panel is a typical turbulent velocity expected from theory; for an r.m.s.Mach number of 9,  turb = 9(/ 0 ) 1/2 (Kritsuk et al. 2007;Collins et al. 2023).This can be seen as the red line in the third panel.Core profiles are bounded below by this turbu-lent velocity.During the hardening, the tangential velocity maintains roughly a profile like  1/2 , decreasing in magnitude.After the singularity, the rotation of the pseudo-core becomes quite large and no longer centered on zero velocity.
In the fourth row of Figure 9 we show the ratio of kinetic to gravitational energy.We define gravitational energy  G and kinetic energy  K as in Equation 17and we plot the ratio of their magnitudes.In the first panel, at  = 0, we find that the ratio tends to decrease with radius (though not universally).This can be understood again as a signature of the turbulent initial conditions.The kinetic energy,  K , increases as roughly , since the total velocity is subject to turbulent statistics, and () ∝  1/2 .On the other hand, we integrate over several density correlation lengths, so that average density is roughly constant and tends to the mean density.Thus their ratio is roughly linear in radius.This is a little counter intuitive, and can be understood because one preimage region covers many density fluctuations, so there is on average no net acceleration vector.During the hardening phase,  G and  K both increase until they are moreor-less equal.Large values of  K / G are seen at the center, while the ratio adjusts to be between 0.5 and 2 by 10,000 AU.During the singularity between  sing and  sung , the material in the pseudo-core begins to rotate, and the ratio of energies in the envelope aligns to be roughly 2 at the center, between 0.5 and 2 for larger radii.

CONCLUSIONS
In this work, we examine how the high vacuum of outer space becomes stars.We did this by embedding a fleet of pseudo-Lagrangian tracer particles in a turbulent, collapsing molecular cloud; identifying dense cores and the particles within; and then examining the density, velocity, and energetics of each particle as they collapse.
We find that collapse happens in three stages; collection, hardening, and singularity.The specifics are unique to each core, and many cores do not experience the collection stage.During collection, low density, high velocity gas is swept into moderate density object by converging flows.Sometimes these are initially dominated by gravity, sometimes by kinetic energy.During the hardening phase, gravitational and kinetic energy both increase to rough equipartition and the velocity decays to near transsonic values.During the singularity, the density gets very high, mass is delivered to the central object.This is a universally fast event, lasting ∼ 0.1 ff .At the end of singularity, a relaxed envelope remains, with kinetic-to-gravitational energy ratio that varies logarithmically from about 2 at the edge of the pseudocore to about 0.5-2 at 10,000 AU.We demonstrate this in detail for one representative object in Section 3.2, and the variation in mean quantities in Section 3.3.
This can be seen as two related transitions to coherence: one in which velocity of the gas that will collapse becomes low; and one in which the kinetic and gravitational energies equilibrate.The first is analogous to those subsonic or trans-sonic cores (Pineda et al. 2010;Singh et al. 2021;Choudhury et al. 2021).The second is the fast equillibration of energies in the envelope.
Our study shows similarities with other recent works.The three phase view of Offner et al. (2022) used machine learning techniques on similar simulations to organize the flow.They find that core experience three phases, turbulence, coherence, and collapse, which are similar in spirit to our collection, hardening, and singularity.In their work, cores move between all phases, while our cores do not.One particular difference is our selection bias.We only examine cores that do not interact with other dense objects.Similarly the inertial flow model (Pelkonen et al. 2021) has high velocity streamers delivering the bulk of the mass to the core, which can be seen during the singularity stage in our work.These works are quite similar in setup and simulation software to our own, differences in conclusions are due mostly to differences in the lenses we use to view our simulations.The physics of self-gravitating turbulence is central to astrophysics, but underpinned by chaos, so many views are necessary to unpack the true dynamics.
This work is a low resolution pilot to gain basic insight into this technique and what it might show us.Simulations with higher resolution and improved sink particles are on the way.Here we have made an attempt to only present results that will stand up to improved resolution.By improving resolution, one increases the power of both density and velocity at the small scales.Small scale density power is not likely to affect our results, as we have shown many density fluctuations are gathered together in the collapse process.The increased velocity power, on the other hand, is known to delay star formation, so our discussion of the distribution of  sing will be affected.It is unlikely that more resolution will introduce any features into this essentially featureless distribution.The most that increased resolution will do is delay the first collapse, shifting the distribution to higher  sing .It is also unlikely to alter the duration of the collapse as this process is dominated by gravity which only depends on the integrated mass.More resolution will certainly allow the collapse to reach higher densities, but again the error in the collapse duration will be small, as the time derivative of density is also singular.The fractal nature of the preimage gas will also be unchanged, as further resolving a fractal yields more fractal, not less, by the nature of a fractal.

Figure 2 .
Figure 2. (Left)The cumulative distribution of the singularity time,  sing .For each mode of collapse, the singularity time is nearly uniformly distributed.(Right) The duration of the singularity, which is distributed between 0.05 ff and 0.15 ff .

Figure 3 .
Figure 3.The collapse of one core and its phases.(Top row) Density (grey lines, right axis) and velocity (radial in red, tangential in teal, and total in black, left axis).Vertical lines denote the times of analysis for the second and third rows.These are the beginning (navy) the middle of the hardening phase (light blue), the end of the hardening phase (light green), the beginning of singularity (dark green) and the "end" of singularity (red).(Middle row) Cumulative density distribution at the 5 indicated snapshots.(Bottom row) Radial profiles of of   and   for each of the 5 snapshots.

Figure 4 .
Figure 4. Extremely rapid assembly of mass.All three plots show azimuthally averaged quantities on a sphere centered on the particles, but taking all zones.The left vertical axis is radius in AU, the horizontal axis is time in units of the free fall time.Vertical black lines show  sing and  sung .(Top) Interior mass,  < .The mass is assembled to the center extremely rapidly between  sing and  sung , with relatively little accretion after  sung .(Center) The relative mass accretion rate.Very little accretion happens after  sing .(Bottom) The ratio of kinetic and gravitational energies.Immediately after  sung , the energies show that the core forms a relaxed envelope.During the singularity, gravitational energy dominates.

Figure 5 .
Figure5.Tracks for all of the single cores from 1, each one scaled to its respective  sing .Each line represents a core.(Top Row) Mean density of preimage gas.A consistent profile can be seen of slow increase followed by free-fall.(Second Row) r.m.s velocity of preimage gas.The total velocity remains quite large.(Third Row) Radial velocity of the preimage.Radial velocity has no preferred initial configuration.Some begin with large convergence, some do not.By  sing , the radial velocity is clustered around transsonic. (Bottom Row) Tangential velocity of the preimage, which is universally supersonic at the beginning, and becomes transsonic by  sing .

Figure 6 .
Figure 6. 1000 is the radius containing a mass of  1000 , which is defined as the total mass within 1000 AU at the end of the simulation.Most cores deliver the bulk of their mass to the center within 0.1  ff of  sing (see Figure2).

Figure 7 .Figure 8 .
Figure 7. (Top) Peak density for each core vs. time normalized to the initial maximum density and stretched in time to  sung .The green line shows the free fall solution, Equation 23. (Bottom) Average density for each core, scaled to the initial average value and stretched in time.

Figure 9 .
Figure9.Radial profiles for several quantities (rows) at several frames (columns) for each single core (grey lines).Four times at t=0,  = 0.5 sing ,  sing , and  sung are shown.Note that  sing and  sung are different for each core.The rows show mean density, radial velocity, tangential velocity, and ratio of gravitational to kinetic energy.Notable features are highlighted by red lines, see text for details.

Table 1 .
Mode of star formation by simulation.

Table 2 .
Pearson R between properties, Q, for each single core and their singularity time,  sing .No substantial correlations are seen.