© The Authors 2025

1 Introduction

For the better part of half a century, the role of ambipolar diffusion in the fragmentation of molecular clouds has been one of the most hotly debated topics in the field of star formation (see e.g., Mouschovias et al. 2006; Crutcher et al. 2010; Mouschovias & Tassis 2010; Crutcher 2012; Tritsis et al. 2015; Priestley et al. 2019; Krumholz & Federrath 2019; Maury et al. 2022 and references therein). This lack of consensus can be partially attributed to substantial challenges that render the exploration of such effects, both theoretically and observationally, exceedingly difficult.

Thus far, only tentative indications of the neutral-ion drift velocity have been found in observations (Li & Houde 2008; Hezareh et al. 2010; Tang et al. 2018). These studies concluded that neutral species have a higher velocity dispersion than charged ones. On the other hand, Pineda et al. (2021, 2024) reached the opposite conclusion, and reported that the charged species exhibit a higher velocity dispersion than the neutrals. Additionally, they found no evidence that the power spectrum of ions has more power at smaller scales. Three key factors hinder our ability to observationally study ambipolar diffusion. Firstly, the neutral-ion drift velocity typically found in numerical simulations (Desch & Mouschovias 2001; Tritsis et al. 2022) is a factor of two to three smaller than the typical spectral resolution of most spectroscopic surveys. Secondly, the observed charged and neutral species used to study ambipolar diffusion need to be co-spatial in order for the difference in the spectral lines to be attributed to the neutral-ion drift velocity (Tassis et al. 2012b). Finally, the optical depth of the two species needs to be comparable so that the derived neutralion drift velocity is not biased by opacity-broadening and other radiative-transfer effects. Tritsis et al. (2023) extensively studied all these effects and, by employing 2D axisymmetric nonideal magnetohydrodynamic (MHD) chemo-dynamical simulations and nonlocal thermodynamic equilibrium (non-LTE) radiativetransfer simulations, presented a roadmap for future observations of ambipolar diffusion.

From a theoretical perspective, it is certain that nonideal MHD effects have to set in at some stage during the collapse of molecular clouds. This is essential to explain the considerably weaker magnetic field observed in young stars compared to the magnetic field that would emerge from a prestellar core collapsing under flux-freezing (e.g., Tsukamoto et al. 2015; Vaytet et al. 2018). Additionally, the onset of nonideal MHD effects is inevitable in the densest parts of molecular clouds where the degree of ionization is on the order of χ_i = 10^–6–10^–9 (Williams et al. 1998; Bergin et al. 1999; Caselli et al. 2002; Goicoechea et al. 2009).

To date, several 3D simulations have incorporated nonideal MHD effects, often, however, employing varying levels of simplifications in their treatment of these processes (Tsukamoto et al. 2015; Masson et al. 2016; Marchand et al. 2016; Dzyurkevich et al. 2017; Vaytet et al. 2018; Wurster et al. 2018; Hennebelle et al. 2020; Zhao et al. 2021; Machida & Basu 2024; Lebreuilly et al. 2024). Most of these calculations focus on the role of nonideal MHD effects in protostellar disks, while only a few are targeted in prestellar cores (e.g., Yin et al. 2021). Additionally, other simulations that focus on earlier stages during the star formation process are usually limited in terms of dimensionality (Basu et al. 2009; Kunz & Mouschovias 2010; Tassis et al. 2012a,b; Tritsis et al. 2022; Tritsis et al. 2023).

Incorporating nonideal MHD effects into our calculations is critical to study variations (both in time and space) of the mass-to-magnetic flux ratio, one of the key parameters in star-formation theories (Mouschovias & Spitzer 1976). Observationally, many studies have attempted to explore such variations¹ (e.g., Falgarone et al. 2008; Crutcher et al. 2009; Koch et al. 2012; Crutcher 2012; Hwang et al. 2021; Palau et al. 2021; Myers & Basu 2021; Ching et al. 2022; Yen et al. 2023). While all these studies have provided valuable insights, only those utilizing Zeeman observations can measure the relevant component of the magnetic field, needed to calculate the flux. The magnetic flux is defined as the dot product between the magnetic field and the normal vector to a given surface. Polarization observations trace the plane-of-sky (POS) component of the magnetic field, while the surface element, ds, is along the line of sight (LOS). Since the two vectors are orthogonal, the magnetic flux measured by polarization observations is, by definition, zero. In other words, there is no flux of the POS component of the magnetic field through the POS. We note that this is not merely a mathematical “formality” but a fundamental constraint that cannot be ignored. However, even for Zeeman observations, the question remains as to how well they can provide reliable information on the mass-to-flux ratio given 3D effects (e.g., Mouschovias & Tassis 2009).

Here, we present results from a nonideal MHD simulation of a collapsing turbulent molecular cloud in which the resistivities are self-consistently calculated from all charge carriers, the abundances of which are computed by a vast nonequilibrium chemical network. More than four million CPU hours and ~2.5 calendar years were required for this simulation to reach a maximum number density of 10⁵ cm^–3. We used this simulation to study the structure of the neutral-ion drift velocity and that of the mass-to-flux ratio.

Our paper is organized as follows. In Sect. 2 we outline the numerical methods we used to perform the nonideal MHD simulation and the physical setup of the collapsing cloud. In Sect. 3.1 we present first results from our simulation on the neutral-ion drift velocity for two different times during the evolution of the cloud. In Sect. 3.2 we discuss our results on the true “3D” mass-to-flux and demonstrate that the “observed” mass-to-flux ratio that would be inferred from Zeeman observations shows a poor correlation with the true values. Finally, in Sect. 4, we provide a summary of our findings and conclude.

2 Numerical methods and setup

To perform our simulation, we largely followed the numerical methods outlined in Tritsis et al. (2022); Tritsis et al. (2023, 2025a,b) with a few key improvements both in relation to the calculation of the resistivities and in relation to the chemical modeling.

2.1 Nonideal MHD

Under nonideal MHD conditions, the magnetic induction equation is $$\frac{\partial \mathit{B}}{\partial t}=\mathbf{\nabla }\times \left({\upsilon }_n\times \mathit{B}\right)-c\mathbf{\nabla }\times \left({\eta }_{\perp }{j}_{\perp }+{\eta }_{\parallel }{j}_{\parallel }+{\eta }_{H}j\times b\right),$$(1)

where B, b, and u_n denote the magnetic field, the unit vector of the magnetic field, and the velocity of the neutrals, respectively. In Eq. (1), j is the current, and j_∥ and j_┴ are the components of the current parallel and perpendicular to the magnetic field. In turn, η_┴, η_∥, and η_Η are the perpendicular, parallel, and Hall resistivities, defined in terms of the respective conductivities, σ_┴, σ_∥, and σ_Η (Parks 1991), as $${\eta }_{\parallel }=\frac{1}{{\sigma }_{\parallel }},\text{}{\eta }_{\perp }=\frac{{\sigma }_{{}_{\perp }}}{{\sigma }_{{}_{\perp }}^2+{\sigma }_H^2},{\eta }_H=\frac{{\sigma }_H}{{\sigma }_{{}_{\perp }}^2+{\sigma }_H^2}.$$(2)

The conductivities are defined as $$\begin{array}{l}{\sigma }_{\parallel }={\displaystyle \sum _s{\sigma }_s,}{\sigma }_{\perp }={\displaystyle \sum _s\frac{{\sigma }_s}{1+{\left({\omega }_s{\text{τ}}_{sn}\right)}^2}},\\ {\sigma }_H=-{\displaystyle \sum _s\frac{{\sigma }_s{\omega }_s{\text{τ}}_{sn}}{1+{\left({\omega }_s{\text{τ}}_{sn}\right)}^2},}\end{array}$$(3)

where ω_s and σ_s are, respectively, the gyrofrequency and conductivity of species “s”, defined as $${\omega }_s=\frac{e_{s}B}{m_{s}c},$$(4) $${\sigma }_s=\frac{n_s{e}_s^2{\text{τ}}_{sn}}{m_s},$$(5)

where n_s and m_s are the number density and mass of charged species s, c is the speed of light, and e_s is the charge of each species. We note that we only considered singly ionized species. However, the symbol e_s carries the algebraic sign of the charge of each species considered. In Eqs. (3) and (5), τ_sn is the mean collisionâ time between species s and neutral particles, defined as (Mouschovias 1996) $${\text{τ}}_{sn}=\frac{1}{{\alpha }_{s\text{He}}}\frac{m_s+m_{{\text{H}}_2}}{{\rho }_{{\text{H}}_2}}\frac{1}{{\overline{\sigma w}}_{s{\text{H}}_2}},$$(6)

where $m_{{\text{H}}_2}$ and ${\rho }_{{\text{H}}_2}$ are the mass and mass density of H₂, ${\overline{\sigma w}}_{s{\text{H}}_2}$ is the mean collisional rate of species s with H₂, and α_sHe is an additional factor to account for the deceleration of species s resulting from collisions with He. For ions the mean collisional rate can be obtained from the Langevin approximation. In contrast to Tritsis et al. (2022) and Tritsis et al. (2023), we did not use constant values for ${\overline{\sigma w}}_{s{\text{H}}_2}$ and α_sHe for all ions. Instead, the mean collision rate was individually determined for each ion as (Osterbrock 1961; Pinto & Galli 2008a) $${\overline{\sigma w}}_{s{\text{H}}_2}=2.21\pi {\left(\frac{e_s^2{p}_{{\text{H}}_2}}{\mu }\right)}^{1/2},$$(7)

where μ is the reduced mass of species s and H₂, and $p_{{\text{H}}_2}$ is the polarizability of H₂. For the factor, α_sHe, we adopted (Pinto & Galli 2008a) $${\alpha }_{s\text{He}}=1+{\left[\frac{\left(m_s+m_{{\text{H}}_2}\right)m_{\text{He}}{p}_{\text{He}}}{\left(m_s+m_{\text{He}}\right)m_{{\text{H}}_2}{p}_{{\text{H}}_2}}\right]}^{1/2}\left(\frac{n_{\text{He}}}{n_{{\text{H}}_2}}\right),$$(8)

where p_He and m_He are the polarizability and mass of He, and $n_{{\text{H}}_2}$ and n_He are the number densities of H₂ and He. Here, we adopted $p_{{\text{H}}_2}=8.04\times {10}^{-25}$ cm³ and p_He = 2.07 × 10^–25 cm³ (Osterbrock 1961). We note that the maximum variation in collisional rate coefficients among different ion species can reach up to ~67% (see, e.g., Figs. 5 and 7 from Pinto & Galli 2008a). In turn, this variation (along with variations in α_sHe) translates into differences in the perpendicular resistivity of ≲ 15%² compared to, say, adopting a constant value for all species. While this is not a drastic deviation at any given time, it is a systematic one that can accumulate over the course of tens of thousands of integration steps. Our initial elemental abundances are listed in Table 1 of Tritsis et al. (2022). For electrons, the Langevin approximation is not applicable. Mouschovias (1996) extrapolated the results by Motte & Massey (1971) to temperatures relevant to molecular clouds, and found a value for the mean collisional rate between electrons and H₂ of ${\overline{\sigma w}}_{e{H}_2}=1.3\times {10}^{-9}$ cm³ s^–1. For the factor α_eHe, we adopted a value of 1.16. Finally, for grains we have that ${\alpha }_{g^{-}He}=1.28$ (Tassis & Mouschovias 2005) and ${\overline{\sigma w}}_{g^{-}{H}_2}$ is given by (Epstein 1924; Ciolek & Mouschovias 1993; Pinto & Galli 2008a) $${\overline{\sigma w}}_{g^{-}{H}_2}=\frac{4}{3}\pi r_g^2\delta {\left(8k_B\text{T}/\pi m_{H_2}\right)}^{1/2},$$(9)

where k_B is the Boltzmann constant, T is the temperature, r_g is the grain radius, and δ is a dimensionless parameter determined from laboratory experiments to be equal to (1.3) (Liu et al. 2003). We adopted an MRN grain size distribution Mathis et al. (1977) following the methods outlined in Kunz & Mouschovias (2009) and considering thirty bins of grain radii with $r_g^{min}=0.0181$ μm and $r_g^{max}=0.9049$ μm. Additionally, we adopted a constant grain abundance of $n_{g^{-}}/n_{{\text{H}}_2}={10}^{-12}$ (Tassis et al. 2012a) and a mass density of ${\rho }_{{\text{g}}^{-}}=2.3$ g cm^–3. A growing body of theoretical and observational studies suggests that grain growth could begin prior to the protostellar phase (e.g., Ormel et al. 2009; Steinacker et al. 2010; Pagani et al. 2010; Miotello et al. 2014; Galametz et al. 2019; Valdivia et al. 2019; Guillet et al. 2020; Bate 2022). At the same time, during the early phases of molecular cloud evolution, grains can contribute up to ~10% to the perpendicular conductivity, depending on the underlying assumptions adopted. While considering grain coagulation and/or a different grain size distribution could moderately influence our results, adopting an MRN grain-size distribution remains a reasonable assumption for the range of densities considered here.

We note here that the Langevin approximation is valid as long as the neutral-ion drift velocity (v_dr) is less than the speed of sound (Mouschovias 1996). However, if magnetic-field gradients become steep enough, this assumption may be invalidated (Hopkins et al. 2024). In our implementation we accounted for such discrepancies using the following approach. For every grid point we computed the resistivities and drift velocities based on the Langevin rates for the mean collisional rates. If the drift-velocity in the grid point was supersonic, we computed new mean collisional rates taking into account the dependence of the mean collisional rate on the drift velocity itself (see e.g., Figs. 3, 5, and 20 from Pinto & Galli 2008a). To calculate the new collisional rates for the ions and the grains we used the relations provided by Pinto & Galli (2008a,b), whereas for electrons we used Eq. (14) from Draine et al. (1983; see also Phelps 1979), adjusted such that, when the electron-neutral drift velocity is zero, the value of the mean collisional rate was ${\overline{\sigma w}}_{e{H}_2}=1.3\times {10}^{-9}$ cm³ s^–1, as was stated above. With these new collisional rates, we calculated new values for the resistivities. Finally, the process continued iteratively until the resistivity values converged to within a 1% tolerance. Given the weak dependence of the mean collisional rates on the drift velocities and the fact that the drift velocities in our simulation only become marginally supersonic, only very minor adjustments on the final values of the resistivities were observed.

2.2 Chemical network

The inclusion of a chemical network is imperative to accurately calculate the resistivities. Specifically, the chemical network was used to compute the number densities of species s (i.e., n_s that first appeared in Eq. (5)) in every point in space and time. We used a nonequilibrium gas-grain chemical network to model the evolution of 115 species described in Tritsis et al. (2025a; 2025b; see also, Tritsis et al. 2016; Tritsis et al. 2022; Tritsis et al. 2023). We adopted reaction rates from the fifth release of the UMIST database (McElroy et al. 2013). The network consisted of ~ 1650 chemical reactions and included a range of processes from freeze-out to cosmic-ray and photo-related processes. The visual extinction that is used to determine the reaction rates of photo-related processes was calculated “on the fly” for every grid point based on the six-ray approximation, following a similar approach as Glover & Clark (2012). For the cosmic-ray ionization rate we adopted the “standard” value of ζ₀ = 1.3 × 10^–17 s^–1 (Caselli et al. 1998). Given its role in setting the ionization fraction, the cosmic-ray ionization rate is a critical free parameter in nonideal MHD simulations. Even modest variations, as small as a factor of two, can significantly impact the dynamical evolution of a cloud, as well as the magnitude and spatial structure of the drift velocity (see e.g., Figs. 1 and 3 from Tritsis et al. 2023). Further observational efforts, such as those by Sabatini et al. (2023), are therefore crucial for placing tighter constraints on its value.

2.3 Numerical setup

To perform our simulation of a collapsing turbulent molecular cloud, we made use of the adaptive mesh refinement (AMR) code FLASH (Fryxell et al. 2000; Dubey et al. 2008). We used an identical numerical and physical setup as in Tritsis et al. (2025a,b) with the only exception being the inclusion of nonideal MHD effects and a different initial condition for the value of the magnetic field. Below, we provide only a brief description of our setup.

Our initial values for the number density and magnetic-field strength were set equal to 500 cm^–3 and 10 μG, respectively. Given our initial strength for the magnetic field, which was initially oriented along the ɀ axis, the value of mass-to-flux ratio (in units of the critical value for collapse Mouschovias & Spitzer 1976), was ~1.7. We assumed isothermality, with the temperature of the cloud being set equal to T = 10 K. The dimensions of our simulated cloud were equal to 2 pc, we used open boundary conditions, and the initial size of our grid was 64³ grid points.

However, we also included two levels of AMR and, as a result, the resolution of our smallest cell was ~8 × 10^–3 pc. The total mass in our simulation box was 240 M_⊙. Turbulence was also initiated as described in Tritsis et al. (2025a), based on the publicly available code TurbGen (Federrath et al. 2022; see also Federrath et al. 2010). Given our initial values for the temperature and strength of the magnetic field, the initial sonic and Alfvén Mach numbers were ℳ_s = σ_υ/c_s ~ 3 and ℳ_A = σ_υ/υ_A ~ 0.85, respectively. We emphasize that by adopting the initial conditions described above we do not imply that these necessarily represent the typical initial conditions of real molecular clouds. Instead, this simulation is part of a broader suite designed to systematically explore the impact of different parameters on the dynamics and chemistry.

3 Results

3.1 The spatial structure and time evolution of the neutral-ion drift velocity

In Fig. 1 we present H₂ number density slices through the computational cube perpendicular to the y (upper row), x (middle row), and z (bottom row) directions, respectively. In the left column we present our results for one free-fall time (hereinafter t_ff), and in the right column we show our results at the end of our simulation, corresponding to 1.44 t_ff. Slices are taken at the location of the maximum density. The overplotted cyan lines show the morphology of the magnetic field for each slice and the blue contours mark the 30, 50, 70, and 90% of the maximum density in this specific slice (dotted, dashed, dash-dotted, and solid lines, respectively). The velocity field of the neutrals is shown with the red-yellow arrows and respective color maps, and the neutral-ion drift velocity is shown with the cyan-fuchsia arrows.

From Fig. 1 several notable features are visible. Firstly, the neutral-ion drift velocity peaks in the surroundings of the densest regions of the cloud rather than precisely at the location of the maximum density. This trend stands for both times during the evolution of the cloud examined here. For purely gravitationally-driven ambipolar diffusion, the neutral-ion drift velocity should scale as υ_dr ∝ –τ_ni┴g_┴, where g_┴ is the component of the gravitational field perpendicular to magnetic field lines (Mouschovias 1987). Given that g_┴ vanishes at the center of the cloud, so should the neutral-ion drift velocity. Additionally, to first order, the perpendicular resistivity scales as ${\eta }_{\perp }\text{}\propto {\upsilon }_a^2/n_{i,tot}$, where υ_a is the Alfvén velocity and n_i,tot is the total number density of all ions (see Eq. (11) from Agianoglou et al. 2025). At a time of 1.44 t_ff, both the total ion density and the Alfvén velocity decrease by a similar factor (~5–7) as the density increases from few times 10⁴ cm^–3 to 10⁵ cm^–3. However, since the perpendicular resistivity depends on the square of the Alfvén velocity, the reduction in υ_a has a stronger impact on η_┴ than the change in total ion density. Additionally, the “midplane” (e.g., ɀ ≈ 0) of the cloud, where the density reaches its maximum value, is a zero-crossing point for the perpendicular components of the current. Since υ_dr ∝ η_┴ J_┴ (see Eq. (A.1)), the drift velocity is suppressed in the densest region and peaks in the surroundings. Finally, the low-density regions of the cloud (~10²–10³ cm^–3 ) remain close to the ideal MHD regime, as expected.

The neutral-ion drift velocity is also found to be primarily perpendicular to the velocity of the neutrals. This can be understood as follows. Along magnetic field lines the neutral-ion drift velocity is nearly zero since both the ions and the neutrals are subject to the same forces. Therefore, the only significant velocity separation occurs for v_dr,⊥. Given that, for the most part, the velocity of the neutrals follows magnetic-field lines, the drift velocity is perpendicular to both the velocity of the neutrals and the magnetic field. This is particularly evident in regions where the magnetic field has the most pronounced hourglass morphology.

From the left column of Fig. 1 it can also be seen that the neutral-ion drift velocity is in very good agreement with the expectations for the gravitationally driven ambipolar diffusion from 2D axisymmetric simulations. That is, the vectors of the neutral-ion drift velocity point inward, toward the center of the cloud with this being the case both above and below the mid-plane of the cloud. This is also evident in the ɀ slice through our computational cube (lower left panel in Fig. 1) where, again, the majority of the vectors of the neutral-ion drift velocity point toward the densest regions of our simulated cloud. However, even at this stage during the evolution of the cloud, a few vectors are more randomly oriented compared to the case of a simple isolated collapsing core with an hourglass morphology. Consequently, some kinetic (or “gravo-kinetic”) ambipolar diffusion is also present, although this does not seem to be the dominant effect at this stage. A radial profile of the magnitude of the drift velocity perpendicular to the mean direction of the magnetic field (not shown here) also looks remarkably similar to 2D axisymmetric nonideal MHD simulations.

Despite this very good agreement with gravitationally driven ambipolar diffusion at early stages during the evolution of the cloud, the situation is drastically different at later times. Especially in the slice perpendicular to the ɀ axis and for t = 1.44 t_ff (bottom right panel in Fig. 1), the neutral-ion drift velocity exhibits a very complicated structure. Of particular interest are the two elongated structures at x ≈ –0.2 pc (with y ≈ –1–0 pc), and x ≈ –0.1 pc (with y ≈ 0.2–0.75 pc) created by solenoidal motions where the magnetic field also exhibits a highly pinched morphology. In both these regions, the neutral-ion drift velocity is oriented perpendicular to the long axis of these elongated structures, suggesting that the neutrals are dispersing outward from the spine of these filamentary structures compared to the ions. However, the most striking feature from the right column of Fig. 1 is the “messy” structure of the neutral-ion drift velocity in the x and y slices. That is, many vectors point outward instead of inward, toward the densest parts of the cloud. Additionally, the direction of the neutral-ion drift velocity above and below the miplane of the cloud, is in “antiphase”.

To further emphasize this point, we show in Fig. 2 1D profiles of the x and y components of the neutral-ion drift velocity (left and right columns, respectively) for t = 1. t_ff and t = 1.44 t_ff (upper and lower rows, respectively). With the thin solid lines we show profiles on the midplane of the cloud and with the thick dash-dotted and dotted lines we show profiles above and below the midplane, respectively. All 1D profiles were obtained by averaging over a 0.05 pc slice in the z direction. The profiles above and below the midplane were taken at a vertical distance of 0.1 pc from the midplane. As is evident from the bottom row of Fig. 2, at a time of t = 1.44 t_ff, the 1D profiles of the drift velocity above and below the midplane are approximately in antiphase. We have confirmed that this is the case in other cuts through our computational cube.

To identify the physical mechanism behind this approximate anticorrelation, we decomposed the drift velocity into two components: one driven by the magnetic tension force per unit volume and one driven by the magnetic pressure force per unit volume (see Appendix A). We note that this separation does not imply that one component is exclusively associated to gravita-tionally driven ambipolar diffusion and the other with kinetic ambipolar diffusion, as both gravity and turbulence can bend and/or compress the magnetic field. In Fig. 3, we show these two components of the drift velocity for a time of 1.44 t_ff. In the upper row we show the x and y components of the drift velocity (left and right panels, respectively) driven by the magnetic tension force per unit volume. In the lower row we show the corresponding components of neutral-ion drift velocity associated to the magnetic pressure force per unit volume. The linestyles have the same meaning as in Fig. 2. From Fig. 3 we can draw two conclusions. Firstly, the components of the neutral-ion drift velocity associated with the magnetic tension force per unit volume are notably higher (i.e., by a fact of ≳2). Secondly, the approximate anticorrelation of the drift velocity above and below the midplane of the cloud is primarily driven by the magnetic tension force per unit volume.

Upon inspecting the 3D structure of the magnetic field in our simulation at a time of 1.44 t_ff, we identified that a significant portion of magnetic field lines exhibit helical “loops” near the midplane of the cloud. These helical loops are in addition to the hourglass morphology which is clearly evident in the upper right, and middle right panels of Fig. 1. In turn, such helical structures arise in regions of high vorticity (see also the velocity of the neutrals in the bottom right panel of Fig. 1). In such magnetic field morphologies, the vector of the magnetic tension force per unit volume points toward the inner regions of the loop, both above and below the midplane. If we ignore any “contributions” by the magnetic pressure force per unit volume to the neutral-ion drift velocity, the latter must have an opposite sign to the magnetic tension force per unit volume. Therefore, the neutral-ion drift velocity should point in one direction above the midplane and in the opposite direction below the midplane. This leads to the anticorrelation of the 1D profiles observed in Figs. 1–3. Additionally, such helical loops lead to a situation where a significant portion of vectors of the neutral-ion drift velocity do not point toward the inner regions of the cloud.

This physical picture is schematically explained on the right side of Fig. 4. With the green arrows we show the direction of the magnetic tension force per unit volume and with the orange vectors we show the direction of the neutral-ion drift velocity. The morphology of the magnetic field lines is depicted with the color-coded tubes. To avoid cluttering, we have not indicated the direction of the magnetic tension force at the region of the helical loop closer to the center of the cloud. The magnetic tension force in this location would point away from the center of the cloud, leading to a drift velocity pointing toward the cloud’s center. Hence, some vectors of the neutral-ion drift velocity should still point inward. On the left side of Fig. 4, we show the structure of the neutral-ion drift velocity (and that of the magnetic tension force per unit volume) expected from a simple hourglass morphology in the magnetic field (consistent with 1 t_ff). We emphasize here again that we have ignored any contributions from the magnetic pressure force per unit volume which are, in reality, non-negligible, and can further complicate the picture.

Fig. 1 Density slices (black-white color map) through our computational cube at the location of the maximum density at one free-fall time (left column) and at the end of our simulation (right column). From top to bottom we show slices through the y, x, and z directions. Cyan streamlines show the morphology of the magnetic field in each slice and blue contours mark the density with dotted, dashed, dash-dotted and solid lines marking, respectively, the 30, 50, 70, and 90% of the maximum density. The red-yellow color map in each panel and respective arrows show the velocity of the neutrals in the simulation and the cyan-fuchsia color map and respective arrows show the neutral-ion drift velocity.

Fig. 2 One-dimensional profiles of the x and y components of the neutral-ion drift (left and right panels, respectively) for two different times during the evolution of the cloud. In the upper row, we present our results for one free-fall time and in the bottom row we show our results at the end of the simulation. With the thin solid lines, we show the 1D profiles on the midplane of the cloud (i.e., z ≈ 0), and with the dash dotted and dotted lines we show the profiles above and below the midplane. Notice the approximate anticorrelation of the profiles above and below the midplane at the end of simulation.

Fig. 3 One-dimensional profiles of the x and y components of the neutral-ion drift (left and right panels, respectively) at the end of the simulation. In the upper row, we show the components of the neutralion drift velocity driven by the magnetic tension force per unit volume, and in the lower row we show the components driven by the magnetic pressure force per unit volume. Note that, to zeroth order, the features observed in the lower row of Fig. 2 can be explained by the component of the neutral-ion drift velocity driven by the magnetic tension force per unit volume.

Fig. 4 Illustration of the physical picture that explains (to zeroth order) the features seen in the drift velocity in Fig. 1 at two different times during the evolution of the cloud. Color-coded tubes represent magnetic field lines, green arrows show the direction of the magnetic tension force per unit volume, and orange arrows show the direction of the neutralion drift velocity. On the left side of the plot we depict the situation at one free-fall time, when magnetic-field lines have approximately an hourglass morphology. The right side of the plot shows the evolved magnetic-field morphology (at 1.44 tff), which explains the anticorre-lation seen in the neutral-ion drift velocity with respect to the midplane of the cloud.

Fig. 5 Evolution of the mass-to-flux ratio as a function of time in a slice perpendicular to the z axis, centered on the location of the maximum density at a time of 1.44 tff. The mass-to-flux ratio was computed for all cells within a 0.2×0.2 pc2 area. With the black points we show the 50th percentile from all pixels in the region and the error bars represent the 16th and 84th percentiles, respectively. With the red points we have zoomed further in to the location where the density eventually increases. The area of this region is 0.01 pc2. The dashed red line marks the initial value of the mass-to-flux ratio in the simulation.

3.2 The mass-to-flux ratio

3.2.1 Average properties of the mass-to-flux ratio in time and space

The neutral-ion drift velocity is of paramount importance in the evolution of the cloud, as it governs the spatial variations in the mass-to-flux ratio. Specifically, the mass-to-flux ratio should increase in the same direction as the direction of the neutral-ion drift velocity. The particularly messy morphology of the neutral-ion drift velocity at later times, and the fact that a significant fraction of the vectors point outward from the dense region raise an important and provocative question, of whether the chaotic structure of the neutral-ion drift velocity could lead to a reduction in the mass-to-flux ratio in the densest region.

To answer this question, we developed a new method to measure the true 3D mass-to-flux ratio. By true mass-to-flux ratio, we refer to what is also known in the literature as the differential mass-to-flux ratio (e.g., Mouschovias 1991). Our approach is described in Appendix B.1. In Fig. 5 we show the time evolution of the mass-to-flux ratio in a region centered on the location of the maximum density at a time of 1.44 t_ff, and spanning 0.2 pc in the x and y directions. We then follow the evolution of the mass-to-flux ratio in this region through time. The black points represent the 50th percentile from all pixels in the region and the error bars represent the 16th and 84th percentiles, respectively. With the horizontal red line we depict our initial condition for the mass-to-flux ratio. As is evident from Fig. 5, the mass-to-flux ratio in the region monotonically increases, with only a very tentative indication for a slight reduction of the median value from all pixels at very late times. To examine to what extent this slight reduction depends on our choice of the region and averaging effects, we repeat the procedure outlined above, this time selecting a zoomed-in region that spans 0.1 pc in the x and y directions, and it is again centered on the location of the maximum density at a time of 1.44 t_ff. We present our results with the red points in Fig. 5. In this zoomed-in region the monotonie increase in the mass-to-flux ratio as a function of time is even clearer.

We now proceed to examine the spatial features of the mass-to-flux ratio at a time of 1.44 t_ff. To this end, we select a larger region spanning, 0.6 pc in the x and y directions and is again centered on the location of the maximum density. We then create radial profiles of the true mass-to-flux ratio and the density by considering concentric radial shells. We present our results on the true mass-to-flux ratio with the blue points in the upper panel of Fig. 6. The points represent the 50th percentile from all pixels in the shell and the error bars represent the 16th and 84th percentiles, respectively. Our results on the radial profile of the density are shown in the bottom panel of Fig. 6 (radial profiles of the mass-to-flux ratio and the density at a time of 1 t_ff are presented in Appendix C). As is evident, the average mass-to-flux ratio is in good agreement with the 1D density profile in that it peaks in the same location as the density and gradually declining outward, until it reaches a value close to the initial condition. Therefore, based on Figs. 5 and 6 we conclude that, at least in terms of its average properties, the complex structure of the neutral-ion drift velocity cannot lead to a reduction in the mass-to-flux ratio.

Despite this conclusion, we note two important points. Firstly, there seem to be additional small features in the radial profile of the mass-to-flux ratio (i.e., at r ≈ 0.1 pc) which are not present in the radial profile of the density. Additionally, other features (i.e., at r ≈ 0.2 pc) seem to have an offset with respect to the density profile. However, given our systematic uncertainties on the individual pixel level (see the discussion in Appendix B.2), it remains uncertain whether these minor discrepancies hold any significance. Secondly, and most importantly, while the average 1D profiles of the mass-to-flux ratio and density align well, their peak locations in the full 2D maps do not coincide. We will further explore this issue in a future study.

To understand why the mass-to-flux ratio does not decrease at late stages, we quantified the direction of the drift velocity in the cloud. First, we considered slices perpendicular to the z axis in the range | z |≤ 0.25 pc. In each slice, we definee a reference point (x_m, y_n), corresponding to the location of maximum density in the midplane slice. For every pixel within a given slice, we constructed a vector pointing toward (x_m, y_n) and computed the angle ${\theta }_{\measuredangle }$ between this vector and the local drift velocity in the x – y plane. For simplicity, we set the z component of the drift velocity to zero, as vertical motions will ultimately be regulated by gravity. In other words, the important physical process here is the direction in which the neutrals will cross magnetic field lines in different planes perpendicular to the mean component of the field.

We present our results on ${\theta }_{\measuredangle }$ in Fig. 7. With the black and red lines we show the probability density functions of the ${\theta }_{\measuredangle }$ for a time of 1 t_ff and 1.44 t_ff, respectively. ${\theta }_{\measuredangle }=0°$ signifies that the drift velocity points toward the dense region, and ${\theta }_{\measuredangle }=180°$ means that the drift velocity points away from the cloud’s center. As is evident from Fig. 7, at 1 t_ff the vast majority of the vectors of the drift velocity point toward the center of the cloud, as expected from the simple physical picture shown on the left side of Fig. 4. At later times (red line in Fig. 7), a significant fraction of vectors point outward from the cloud’s center. Yet, the majority of the vectors still point inward or at ${\theta }_{\measuredangle }\approx 90°$. This can again be explained based on the simple physical picture shown on the right side of Fig. 4. We note that ${\theta }_{\measuredangle }\approx 90°$ should not significantly affect the mass-to-flux ratio. A combination of drift velocities oriented at ${\theta }_{\measuredangle }\approx 90°$ across different locations would result in a “flow” pattern resembling a circular vortex causing neutrals to “circulate” around the central region (relative to the ions), rather than escape it. Our results do not qualitatively change if we consider fewer slices and/or if we consider only the drift velocities with magnitude above some threshold (e.g., ∥υ_dr∥ ≥ 0.01 km s^–1).

Fig. 6 Average profiles of the mass-to-flux ratio (upper panel), and the number density (lower panel) of the cloud in concentric radial shells at a time of 1.44 tff The averages were created by considering a slice perpendicular to the z axis centered at the location of the maximum density. With the blue points in the upper panel we show the true mass-to-flux ratio, and with the magenta points we show the observed mass-to-flux ratio. Even though the drift velocity exhibits complex features, there is still a good correlation in terms of averages between the true mass-to-flux ratio and the number density of the cloud. On the other hand, the radial profile of the observed mass-to-flux ratio does not correspond well to either the true mass-to-flux ratio or the number density of the cloud.

3.2.2 The observed mass-to-flux ratio

Examining all the complications that enter in observational measurements of the mass-to-flux ratio through Zeeman observations requires non-LTE spectropolarimetric radiative transfer calculations of species such as CN and/or OH. Even though both these species are present in our chemical network, we avoid such complications and consider the following idealized scenario. Firstly, we assume that the cloud is viewed along the mean direction of the magnetic field (i.e., along the z axis of our simulation). Secondly, we assume that our “Zeeman observations” can perfectly probe the LOS component of the magnetic field (i.e., the z component) at the midplane of the cloud, where the density of the cloud peaks. With these assumptions we essentially consider a hypothetical molecule that is perfectly concentrated at the midplane of the cloud, and we ignore any complications arising from chemical and/or radiative transfer effects. We then measure the observed mass-to-flux ratio as described in Appendix B.3. As such, our results should be considered the best-case scenario regarding how well we can probe the mass-to-flux ratio with Zeeman observations.

With the magenta points in the upper panel of Fig. 6, we show the radial profile of the observed mass-to-flux ratio averaged in concentric radial shells. The points and error bars have the same meaning as in the case of the true mass-to-flux ratio. As is evident, the observed mass-to-flux ratio does not accurately reflect the true values, overestimating them in some regions and underestimating them in others. Most notably, the observed mass-to-flux ratio does not monotonically decrease outward. That could lead into a scenario where, if the beam of a telescope was placed at r ≈ 0.24 pc, and at r ≈ 0.17 pc, the mass-to-flux ratio would appear to be decreasing toward the center of the core (Crutcher et al. 2009). These discrepancies arise because of the fact that, observationally, we cannot follow the true 3D structure of magnetic flux tubes or, in other words, we cannot perform a line integral along a flux tube. Therefore, this is a geometrical effect, unrelated to nonideal MHD effects and/or turbulence, which needs to be considered when interpreting observations. We have also confirmed that using a thicker slice (i.e., averaging the magnetic field over multiple layers close to the midplane of the cloud, instead of using a single plane) does not qualitatively alter our results and the correlation of the observed mass-to-flux ratio with the true values does not improve. Finally, we note that the same discrepancies between the observed mass-to-flux ratio and the true one can also be realized in HI narrow self-absorption Zeeman observations (Ching et al. 2022) for the same reason as outlined above.

Fig. 7 Probability density functions of the angle between the drift velocity and vectors pointing toward the center of the cloud for the two different times examined. ${\theta }_{\measuredangle }=0°$ means that the neutral-ion drift points toward the center of the cloud and ${\theta }_{\measuredangle }=180°$ means that the drift velocity points away from the dense region.

4 Summary and conclusions

We performed a nonideal MHD chemo-dynamical simulation of a supercritical transAlfvénic supersonic collapsing molecular cloud. The resistivities were self-consistently calculated from the abundances of the charge carriers in our chemical network. We studied the structure of the neutral-ion drift velocity in two different times during the evolution of the cloud. We developed a novel approach for measuring the true mass-to-flux ratio in MHD simulations. Using this method, we examined the time evolution of the mass-to-flux ratio, and computed radially averaged profiles from the cloud’s center to assess global trends.

At early times, the structure of the neutral-ion drift velocity is in very good agreement to 2D axisymmetric nonideal MHD models that have a simple hourglass magnetic field morphology. That is, the vast majority of the vectors of the neutral-ion drift velocity point toward the center of the cloud, with this being the case both above and below the cloud’s midplane. However, as the cloud evolves, the structure of the drift velocity becomes significantly more complicated. Specifically, we found that the drift velocity above and below the cloud’s midplane is in antiphase. Additionally, we analyzed the angle between the drift velocity and the direction toward the center of the cloud. We found that a substantial fraction of the drift velocity vectors point away from the center of the cloud, albeit the number of vectors pointing inward is still significantly higher. Finally, we found that a another significant fraction of neutrals cross magnetic field lines at angles of ${\theta }_{\measuredangle }\approx 90°$ relative to the direction toward the cloud’s center. We show that this behavior of the drift velocity at later times is caused by the magnetic tension force per unit volume, and the presence of magnetic helical loops, which are formed in high vorticity regions. Based on these results, we conclude that the presence of turbulence can significantly affect the structure of the neutral-ion drift velocity, further challenging the interpretation of observations aiming to quantify ambipolar diffusion.

Despite these complexities, we demonstrate that, when averaged over a central region, the true mass-to-flux ratio still follows the expected monotonic increase over time, and a radial decrease toward the edges of the cloud. This is substantial evidence that the mass accumulation (relative to the magnetic flux) in the central parts of the cloud is not disrupted by the turbulence-induced changes in the structure of the drift velocity, and that the predictions of the ambipolar diffusion theory of star formation (Mouschovias & Ciolek 1999) are still relevant.

Finally, we show that the observed mass-to-flux ratio, even in an idealized scenario and even upon averaging over a region, fails to accurately capture the spatial variations in the true mass-to-flux ratio and shows poor correlation with the underlying density structure of the cloud. This discrepancy highlights the challenges in interpreting observational measurements and the necessity to account for 3D effects in Zeeman measurements of the mass-to-flux ratio.

Acknowledgements

We thank K. Tassis, T. Mouschovias, S. Basu, C. Federrath, J. Schober, G. V. Panopoulou, R. Skalidis, and V. Andrianatou for stimulating discussions and comments. We also thank the referee for a useful report that helped us improve the manuscript. A.T. acknowledges support by the Ambizione grant no. PZ00P2_202199 of the Swiss National Science Foundation (SNSF). The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. This research was enabled in part by support provided by SHARCNET (Shared Hierarchical Academic Research Computing Network) and Compute/Calcul Canada and the Digital Research Alliance of Canada. We also acknowledge use of the following software: MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020) and the YT analysis toolkit (Turk et al. 2011).

Appendix A The neutral-ion drift velocity

By combining the momentum equation for charged species with the generalized Ohm’s law (see Tritsis et al. 2023 for a detailed derivation), it can be shown that $$\left[\begin{array}{ccc}\frac{m_{\text{s}}c}{{\text{τ}}_{\text{sn}}{e}_s}& B_z& -B_y\\ -B_z& \frac{m_{\text{s}}c}{{\text{τ}}_{\text{sn}}{e}_s}& B_x\\ B_y& -B_x& \frac{m_{\text{s}}c}{{\text{τ}}_{\text{sn}}{e}_s}\end{array}\right]\left[\begin{array}{c}{\upsilon }_{\text{dr},x}\\ {\upsilon }_{\text{dr},y}\\ {\upsilon }_{\text{dr},z}\end{array}\right]=c\left[\begin{array}{c}{\eta }_{\perp }{j}_{\perp ,x}+{\eta }_{\parallel }{j}_{\parallel ,x}+{\eta }_{\text{H}}{\left(j+b\right)}_x\\ {\eta }_{\perp }{j}_{\perp ,y}+{\eta }_{\parallel }{j}_{\parallel ,y}+{\eta }_{\text{H}}{\left(j+b\right)}_y\\ {\eta }_{\perp }{j}_{\perp ,z}+{\eta }_{\parallel }{j}_{\parallel ,z}+{\eta }_{\text{H}}{\left(j+b\right)}_z\end{array}\right],$$(A.1)

where all of the symbols have the same meaning as in Sect. 2.1. Additionally, by neglecting the terms containing the parallel and Hall resistivities, as well as other subdominant terms it follows that $$\begin{array}{l}{\upsilon }_{dr,x}=-\frac{C{\eta }_{\perp }}{B^2}{\left(\mathit{j}+\mathit{B}\right)}_x,{\upsilon }_{dr,y}=-\frac{C{\eta }_{\perp }}{B^2}{\left(\mathit{j}+\mathit{B}\right)}_y,\\ {\upsilon }_{dr,z}=-\frac{C{\eta }_{\perp }}{B^2}{\left(\mathit{j}+\mathit{B}\right)}_z,\end{array}$$(A.2)

which can, in turn, be written as $$\begin{array}{c}{\upsilon }_{dr,x}=-\frac{c^2{\eta }_{\perp }}{4\pi B^2}{\left(\left(\mathit{B}\cdot \nabla \right)\mathit{B}\right)}_x+\frac{c^2{\eta }_{\perp }}{8\pi B^2}{\left(\nabla B^2\right)}_x,\\ {\upsilon }_{dr,y}=-\frac{c^2{\eta }_{\perp }}{4\pi B^2}{\left(\left(\mathit{B}\cdot \nabla \right)\mathit{B}\right)}_y+\frac{c^2{\eta }_{\perp }}{8\pi B^2}{\left(\nabla B^2\right)}_y,\\ {\upsilon }_{dr,z}=-\frac{c^2{\eta }_{\perp }}{4\pi B^2}{\left(\left(\mathit{B}\cdot \nabla \right)\mathit{B}\right)}_z+\frac{c^2{\eta }_{\perp }}{8\pi B^2}{\left(\nabla B^2\right)}_z.\end{array}$$(A.3)

Therefore, the neutral-ion drift velocity can be decomposed into two components; one associated to the magnetic tension force per unit volume and the other to the magnetic pressure force per unit volume.

Appendix B Measuring the mass-to-flux ratio

B.1 Calculating the true mass-to-flux ratio

The basic procedure to calculate the true mass-flux-ratio is schematically described in Fig. B.1. We begin by defining a slice perpendicular to the z axis of our simulation at z = z_k. Each grid cell (x_i, y_j) within this slice defines an initial position used to trace magnetic field lines through the computational domain. Integration is performed using a fourth-order accurate Runge-Kutta method (Fehlberg 1969). Since the magnetic-field-line tracing involves arbitrary positions (x′, y′, z′) that may fall in between the cell centers, we use trilinear interpolation to determine the magnetic field components at these locations while preserving the divergence-free condition³. The integration process continues until a traced point is beyond the simulation domain. The same procedure is repeated in the opposite direction (i.e., for z < z_k) by reversing the field direction.

Once we have the points r_n that trace a magnetic field line, we compute the orthogonal unit vectors ê_1,i, ê_2,i, and ê_3,i, at a position “i” along the magnetic field line as $${\widehat{\mathbf{e}}}_{1,i}=\frac{\mathbf{B}\left(r_i\right)}{\Vert \mathbf{B}\left(r_i\right)\Vert },{\widehat{\mathbf{e}}}_{2,i}=\frac{\left(\begin{array}{c}-e_{1y,i}\\ e_{1x,i}\\ 0\end{array}\right)}{\sqrt{e_{1x,i}^2+e_{1y,i}^2}},{\widehat{\mathbf{e}}}_{3,i}={\widehat{\mathbf{e}}}_{1,i}\times {\widehat{\mathbf{e}}}_{2.i}.$$(B.1)

We then extend outward from the field line in the directions given by the vectors ±ê_2,i, ±ê_3,i and the diagonal directions $\pm \left({\widehat{e}}_{2,i}+{\widehat{e}}_{3,i}\right)/\sqrt{2}$ and $\pm \left({\widehat{e}}_{2,i}-{\widehat{e}}_{3,i}\right)/\sqrt{2}$. The same process is repeated for the next location “i + 1” along the magnetic field line. To scale these displacements appropriately, we introduce the scaling factors f_i and f_i+1, which depend on the local magnetic field strength. These factors ensure that the cross-sectional shape adapts to variations in the field strength along the field line (see Fig. B.1), thereby also ensuring that the magnetic flux is constant regardless of the location along the field line. The maximum allowed displacement is one cell.

Fig. B.1 Schematic representation of the process followed to measure the mass-to-flux ratio. First, we follow magnetic field lines in our simulation box (black line). Using some basis vectors in each set of positions “i” and “i + 1” along a magnetic field line, we define convex regions as shown in the inset panel. The volume of these regions is calculated using Delaunay triangulation and the density is interpolated from data in the simulation. The surface “dsi,” used to compute the magnetic flux, is also calculated using triangulation (see text for more details).

Through this process we have a set of 16 vertices that define a convex shape with volume dV_i. The volume dV_i of these shapes is calculated using Delaunay triangulation to break the region into tetrahedra, with the volume of each tetrahedron being $$V_{tetrahedron}=\frac{1}{6}\left|\det \left(a-d,b-d,c-d\right)\right|,$$(B.2)

where a, b, c, and d are the vertices of each individual tetrahedron. The mass density inside each convex shape ${\rho }_{x\prime ,y\prime ,z\prime }$ is calculated using trilinear interpolation. Having computed the mass density and the volume of the shape in question, we then simply compute dm_i. The total mass loading of the flux tube is calculated by repeating the process for all positions along the magnetic field line and adding all the infinitesimal mass elements, as shown in the inset panel of Fig. B.1. The area ds_i, used in the calculation of the magnetic flux, is again calculated using triangulation. The magnetic flux is then computed as $${\Phi }_B=B\left(r_i\right)\cdot \widehat{B}\left(r_i\right)d{s}_i,$$(B.3)

and, as mentioned, is constant regardless of the location at which is measured along an individual magnetic field line.

B.2 Benchmarking the approach for measuring the mass-to-flux ratio

The process introduced in Appendix B.1 to measure the true mass-to-flux ratio may be subject to several sources of errors.

The most significant of such errors stem from the accuracy with which we can trace magnetic field lines. This becomes increasingly challenging as the cloud is collapsing and the magnetic field morphology gets more complex. Another source of error is related to the accuracy with which we can calculate the volume of the individual volume elements. Interpolation in both the magnetic field and the density, can also introduce small errors in the calculation of the individual mass elements, which can however accumulate, especially in regions of the cloud with strong density and/or magnetic-field gradients. Finally, exporting the simulation data into a uniform grid also introduces some errors.

To verify and benchmark our approach for calculating the mass-to-flux ratio, and access the significance of these errors, we analyze the ideal MHD simulation described in Tritsis et al. (2025a,b). In contrast to the simulation described in the present paper, that simulation was supercritical (with a mass-to-flux ratio of 2.28), but was performed using the same base resolution and the same levels of refinement. Given that no nonideal MHD effects were introduced in the simulation performed by Tritsis et al. (2025a,b) we expect to recover the initial mass-to-flux ratio of 2.28.

We begin by defining a slice perpendicular to the z axis that passes through the location of the maximum density at a time of ~1.75 Myrs, when the central number density is 10⁵ cm^–3. This slice spans 0.2 pc in both the x and y directions and its centered in the location of the maximum density. We compute the mass-to-flux ratio in this region, both at a time of ~1.75 Myrs and also follow the evolution of the mass-to-flux ratio in the same region through time. The upper panel of Fig. B.2 shows this time evolution, where the points represent the 50th percentile from all pixels in the region, and the error bars correspond to the 16th and 84th percentiles. The dashed red line in Fig. B.2 marks the initial condition of Μ/Φ_Β = 2.28.

As is evident from the upper panel of Fig. B.2, the measured mass-to-flux ratio is slightly below the initial value of Μ/Φ_Β = 2.28 for t ≲ 1 Myr (corresponding to 0.7 t_ff). This discrepancy is not necessarily a systematic bias as, due to our open boundary conditions and the initial turbulent conditions, some mass can leave the computational domain. At later times, however, the error bars indicate that some cells exhibit larger mass-to-flux ratios than the initial condition. This is due to two factors. The first is some unavoidable numerical diffusion in the simulation itself. The second is that we are more prone to errors as the cloud is collapsing. At late times, as the cloud is collapsing along magnetic field lines, most of the mass is concentrated within a thin slice which is approximately perpendicular to the z direction. As such, interpolation errors in only a small number of dm_i elements are sufficient to overestimate or underestimate the mass-to-flux ratio. Despite these challenges, the measured values for the true mass-to-flux ratio remain very close to the initial value of 2.28 throughout the time evolution of the simulation.

In the lower panel of Fig. B.2 we present radially averaged profiles in concentric rings centered again around the location of the maximum density. The error bars correspond to the 16th and 84th percentiles from all pixels in each concentric ring. Once more, the values measured are very close to the initial condition, albeit somewhat lower due to the same reasons explained above. Notably, there is very little variation as a function of radius, suggesting that no significant spatial biases are present in our mass-to-flux ratio measurements. Based on these tests, we conclude that the systematic errors on the measured mass-to-flux ratio for individual pixels, excluding changes due to nonideal MHD effects, are on the order of 0.2–0.3.

Fig. B.2 Time evolution (upper panel) and radial profiles (lower panel) of the mass-to-flux ratio in an ideal MHD simulation. The dashed red line in both panels shows the initial value of the mass-to-flux ratio in the simulation. While there are some deviations from the original value (see Sect. B.2 for further explanation), the initial mass-to-flux ratio is recovered with notable accuracy.

B.3 Calculating the observed mass-to-flux ratio

In contrast to the process described in Appendix B.1, the “observed” magnetic flux, when the cloud is viewed along the z axis, is calculated as $${\text{Φ}}_{B,obs}=B_z\left(z_k\right)d{x}^2,$$(B.4)

where dx is the resolution of each cell. The observed mass that is (erroneously) assigned to that “flux tube” is calculated as $${\text{M}}_{{\text{ΦB}}_{},\text{obs}}={\displaystyle \sum _{{\text{zk}}_{}=0}^{\text{N}}{\rho }_{{\text{xi}}_{},{\text{yj}}_{},{\text{zk}}_{}}{\text{dx}}^3},$$(B.5)

where Ν is the number of cells in the z direction.

Appendix C Density and mass-to-flux ratio radial profiles at 1 t_ff

Here, we present radial profiles of the cloud’s number density, along with the true and observed mass-to-flux ratios, at a time of 1 t_ff. However, a radial profile of the number density created in the same manner as for 1.44 t_ff is ill defined for two reasons. Firstly, the “denser” structure (the maximum number density at this stage is ~ 7.5 × 10³ cm^–3) is significantly tilted with respect to the principal axes of the grid (see the blue contours in the upper left and middle left panels of Fig. 1). As such, considering a radial profile on the x – y plane would misrepresent the cloud’s number density. Secondly, the cloud has not yet fully thermally relaxed along magnetic field lines to form a flattened configuration with approximately uniform density along its minor axis. Therefore, a one-to-one correspondence between the radial profiles of the true mass-to-flux ratio and the number density should also not be expected since the former is more affected (compared to later times) from the gas distribution further away from the midplane.

To provide a meaningful comparison with the results presented at later times, we constructed the radial profiles using a modified procedure. We first selected a subregion of the cloud, spanning 0.6 pc in each of the x, y, and z directions (i.e., similar to the subregion used in Sect. 3.2), but now extended across all three dimensions. We then performed a weighted principal component analysis (PCA) to the density distribution in order to determine the orientation of the “denser” region within the cloud. In this context, by denser gas we define regions with density exceeding the 84th percentile of the density distribution within the selected subregion. As weights, we used the actual density values. We then rotated the full 3D density structure such that the minor axis of the denser region at this stage aligns with the z direction. This allowed us to define a quasi-planar frame in which we averaged the density (within 0.1 pc) to mitigate density variations along the minor axis. To avoid contamination from interpolation artifacts near the edges of the rotated cube, we masked regions affected by zero-padding by assigning undefined values, which were then excluded from the radial averaging. The radial density profile was finally computed in concentric annuli centered on the peak of the averaged density. The same PCA-based rotation was applied to the 2D maps of the true and observed mass-to-flux ratios prior to creating the radial profiles so that there is one-to-one spatial correspondence with the density.

We present our results on the radial profiles of the mass-to-flux ratio and the number density, in the upper and lower panels of Fig. C.l, respectively. Radial profiles were calculated following the same procedure as in Sect. 3.2. That is, the points represent the 50th percentile from all pixels in each shell and the error bars represent the 16th and 84th percentiles, respectively. Fig. C.l reveals two interesting features. First, similarly to later times, the observed mass-to-flux ratio does not correspond exactly to the true one. This indicates that the projection effects in mass-to-flux ratio observational measurements, discussed in Sect. 3.2, are relevant throughout the cloud’s evolution. Second, the radial profile of the true mass-to-flux ratio does not peak at the same location as the density. The difference in the value of the mass-to-flux ratio between the center and its peak is relatively small (~0.1–0.15), but nonetheless noticeable. As was previously mentioned, this could be attributed to the fact that the true mass-to-flux ratio at this stage would be more affected by the density distribution further away from the mid-plane of the cloud. However, the fact that the true mass-to-flux ratio peaks ~0.1 pc away from the peak location of the density could also be a consequence of the structure of the drift velocity at this, and previous, evolutionary stages. To explore this possibility, we compute the radial profile of the mass-to-flux ratio in a 2D axisymmetric supercritical simulation (i.e., model “Μ/Φ2.6_ζ/ζ₀ l_A_υ 10_T10” from Tritsis et al. 2023). The initial mass-to-flux ratio in this simulation (in units of the critical value) was 2.6 and the radial profile was calculated in the midplane of the cloud. The corresponding profile of the number density in the midplane of the cloud is shown with the dash-dotted black line in the lower panel of Fig. C.l. Despite the significantly simpler morphology of this cloud, and the fact that only gravitationally driven ambipolar diffusion was present, the mass-to-flux ratio is still not centrally peaked, although the effect is even smaller than the 3D simulation presented herein. At later times, the mass-to-flux ratio profile in the 2D simulation (not shown here) also peaks at the center of the cloud (i.e., at r = 0). These results suggest that the structure of the drift velocity in supercritical clouds could contribute to spatial offsets between the peak of the mass-to-flux ratio and that of the density.

Fig. C.1 Average profiles of the mass-to-flux ratio (upper panel), and the number density (lower panel) of the cloud in concentric radial shells at a time of 1 tff. With the blue points in the upper panel we show the true mass-to-flux ratio, and with the magenta points we show the observed mass-to-flux ratio. The dash-dotted black lines in the upper and lower panels show, respectively, radial profiles of the true mass-to-flux ratio and the density from a 2D axisymmetric supercritical simulation (Tritsis et al. 2023).

References

All Figures

Fig. 1 Density slices (black-white color map) through our computational cube at the location of the maximum density at one free-fall time (left column) and at the end of our simulation (right column). From top to bottom we show slices through the y, x, and z directions. Cyan streamlines show the morphology of the magnetic field in each slice and blue contours mark the density with dotted, dashed, dash-dotted and solid lines marking, respectively, the 30, 50, 70, and 90% of the maximum density. The red-yellow color map in each panel and respective arrows show the velocity of the neutrals in the simulation and the cyan-fuchsia color map and respective arrows show the neutral-ion drift velocity.

In the text

Fig. 2 One-dimensional profiles of the x and y components of the neutral-ion drift (left and right panels, respectively) for two different times during the evolution of the cloud. In the upper row, we present our results for one free-fall time and in the bottom row we show our results at the end of the simulation. With the thin solid lines, we show the 1D profiles on the midplane of the cloud (i.e., z ≈ 0), and with the dash dotted and dotted lines we show the profiles above and below the midplane. Notice the approximate anticorrelation of the profiles above and below the midplane at the end of simulation.

In the text

Fig. 3 One-dimensional profiles of the x and y components of the neutral-ion drift (left and right panels, respectively) at the end of the simulation. In the upper row, we show the components of the neutralion drift velocity driven by the magnetic tension force per unit volume, and in the lower row we show the components driven by the magnetic pressure force per unit volume. Note that, to zeroth order, the features observed in the lower row of Fig. 2 can be explained by the component of the neutral-ion drift velocity driven by the magnetic tension force per unit volume.

In the text

Fig. 4 Illustration of the physical picture that explains (to zeroth order) the features seen in the drift velocity in Fig. 1 at two different times during the evolution of the cloud. Color-coded tubes represent magnetic field lines, green arrows show the direction of the magnetic tension force per unit volume, and orange arrows show the direction of the neutralion drift velocity. On the left side of the plot we depict the situation at one free-fall time, when magnetic-field lines have approximately an hourglass morphology. The right side of the plot shows the evolved magnetic-field morphology (at 1.44 tff), which explains the anticorre-lation seen in the neutral-ion drift velocity with respect to the midplane of the cloud.

In the text

Fig. 5 Evolution of the mass-to-flux ratio as a function of time in a slice perpendicular to the z axis, centered on the location of the maximum density at a time of 1.44 tff. The mass-to-flux ratio was computed for all cells within a 0.2×0.2 pc2 area. With the black points we show the 50th percentile from all pixels in the region and the error bars represent the 16th and 84th percentiles, respectively. With the red points we have zoomed further in to the location where the density eventually increases. The area of this region is 0.01 pc2. The dashed red line marks the initial value of the mass-to-flux ratio in the simulation.

In the text

Fig. 6 Average profiles of the mass-to-flux ratio (upper panel), and the number density (lower panel) of the cloud in concentric radial shells at a time of 1.44 tff The averages were created by considering a slice perpendicular to the z axis centered at the location of the maximum density. With the blue points in the upper panel we show the true mass-to-flux ratio, and with the magenta points we show the observed mass-to-flux ratio. Even though the drift velocity exhibits complex features, there is still a good correlation in terms of averages between the true mass-to-flux ratio and the number density of the cloud. On the other hand, the radial profile of the observed mass-to-flux ratio does not correspond well to either the true mass-to-flux ratio or the number density of the cloud.

In the text

	Fig. 7 Probability density functions of the angle between the drift velocity and vectors pointing toward the center of the cloud for the two different times examined. ${\theta }_{\measuredangle }=0°$ means that the neutral-ion drift points toward the center of the cloud and ${\theta }_{\measuredangle }=180°$ means that the drift velocity points away from the dense region.
In the text

Fig. B.1 Schematic representation of the process followed to measure the mass-to-flux ratio. First, we follow magnetic field lines in our simulation box (black line). Using some basis vectors in each set of positions “i” and “i + 1” along a magnetic field line, we define convex regions as shown in the inset panel. The volume of these regions is calculated using Delaunay triangulation and the density is interpolated from data in the simulation. The surface “dsi,” used to compute the magnetic flux, is also calculated using triangulation (see text for more details).

In the text

	Fig. B.2 Time evolution (upper panel) and radial profiles (lower panel) of the mass-to-flux ratio in an ideal MHD simulation. The dashed red line in both panels shows the initial value of the mass-to-flux ratio in the simulation. While there are some deviations from the original value (see Sect. B.2 for further explanation), the initial mass-to-flux ratio is recovered with notable accuracy.
In the text

Fig. C.1 Average profiles of the mass-to-flux ratio (upper panel), and the number density (lower panel) of the cloud in concentric radial shells at a time of 1 tff. With the blue points in the upper panel we show the true mass-to-flux ratio, and with the magenta points we show the observed mass-to-flux ratio. The dash-dotted black lines in the upper and lower panels show, respectively, radial profiles of the true mass-to-flux ratio and the density from a 2D axisymmetric supercritical simulation (Tritsis et al. 2023).

In the text