© The Authors 2025

1 Introduction

Nonideal magnetohydrodynamic (MHD) effects, and particularly ambipolar diffusion, play a crucial role in star formation (e.g., Mouschovias & Ciolek 1999). They allow the collapse of cores with a magnetic field initially strong enough to prevent collapse by redistributing the magnetic flux in a cloud (e.g., Basu et al. 2009; Kunz & Mouschovias 2010; Tassis et al. 2012b; Tritsis et al. 2022). Such a redistribution of magnetic flux is necessary to explain why stars do not inherit all the magnetic flux of their parent clouds, thereby providing a solution to the “magnetic-flux problem” in star formation (e.g., Tsukamoto et al. 2015). Finally, a reduced coupling between the magnetic field and the gas is required to mitigate the effect of magnetic breaking, thus allowing rotationally supported disks to form (see e.g., Wurster & Lewis 2020 and references therein).

Numerically, there are two major complications when modeling nonideal MHD effects. The first challenge is that, in explicit codes, the maximum allowed time step for stability, determined by the Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1928), needs to be reduced by many orders of magnitude in comparison to the ideal MHD regime in order to account for the diffusion processes introduced (Mac Low et al. 1995; Desch & Mouschovias 2001). This can lead to prohibitively small time steps and thus very computationally expensive simulations. The second major issue stems from the need to include a chemical network in order to calculate the abundances of ions from which the resistivities can then be self-consistently calculated. This often involves solving hundreds to thousands of chemical reactions in every cell in the simulation box and for every time step. These two issues become especially problematic for highresolution 3D simulations, as they hinder efforts to study such effects in 3D and fully explore the parameter-space.

Several attempts have been made over the past couple of decades to mitigate these issues. Li et al. (2006) proposed the so-called heavy-ion approximation, in which the time step is increased while maintaining the ambipolar diffusion timescale. However, this approach is not physically motivated as the mass of the ions is artificially increased, and its applicability to turbulent clouds is problematic. In other studies (e.g., Abe et al. 2024, Tsukamoto et al. 2022), the calculation of the resistivities is simplified by using simple power-law scaling with the density for either the ion density or the resistivities. However, the extent to which the scaling of the resistivities with the gas density can be accurately approximated by a single power law throughout the evolution of a cloud remains ambiguous.

Here, we aim to solve the second major issue when performing nonideal MHD simulations, that is, the need to include a chemical network. While we also use a power-law for approximating the ionization fraction in the cloud as a function of the H₂ density, in our implementation the power law evolves alongside the collapsing cloud. This approximation and its evolution depend on the physical conditions of any given model. Therefore, an interpolation function over the physical parameters inside the cloud (temperature, cosmic-ray ionization rate, visual extinction, etc.) is used to generalize the method for a large parameter space. With this new method, the resistivities of a cloud can be accurately calculated for various physical conditions without the need to calculate any chemical reactions.

We note here that the contribution of grains in the resistivities is not taken into account in this new approach. Consequently, our approximations are valid up to a density of ~10⁶ cm⁻³ and our method is specifically targeted at parsec-scale star formation simulations that focus on the fragmentation and core formation phase (Clark & Glover 2014; Inoue et al. 2018; Khullar et al. 2019; Appel et al. 2023). Such simulations typically incorporate sink particles at densities of 10⁷ cm⁻³, effectively bypassing the need for a detailed treatment of nonideal MHD effects beyond these densities. Therefore, our method provides an efficient and practical way to include nonideal MHD effects for the majority of the time evolution of such numerical studies.

This paper is organized as follows: in Sect. 2, we derive approximate expressions for the perpendicular, parallel and Hall resistivities. In Sect. 3, we describe the setup of 2D and 3D nonideal MHD (chemo-)dynamical simulations of collapsing molecular clouds and cores used to develop our method and benchmark our results. In Sect. 4, we present our findings, and in Sect. 6, we compare our results with other approximate methods for calculating the resistivities and conclude. A description of the numerical details of our implementation and a general outline of the code we used can be found in Appendices A and B.

2 Basic equations

The expressions for computing the parallel, perpendicular and Hall resistivities that appear in the generalized Ohm’s law (Kunz & Mouschovias 2009) are $$\begin{array}{l}{\eta }_{\Vert }=\frac{1}{{\sigma }_{\Vert }},{\eta }_{\perp }=\frac{{\sigma }_{\perp }}{{\sigma }_{\perp }^2+{\sigma }_H^2},\\ {\eta }_H=\frac{{\sigma }_H}{{\sigma }_{\perp }^2+{\sigma }_H^2}.\end{array}$$(1)

In Eq. (1), σ_‖, σ_⊥, and σ_H are, respectively, the parallel, perpendicular, and Hall conductivities. The conductivities are, in turn, defined (Parks 1991) as $$\begin{array}{l}{\sigma }_{\parallel }={\displaystyle \sum _s{\sigma }_s,}{\sigma }_{\perp }={\displaystyle \sum _s\frac{{\sigma }_s}{1+{(w_s{\text{τ}}_{sn})}^2}},\\ {\sigma }_H=-{\displaystyle \sum _s\frac{{\sigma }_s{w}_s{\text{τ}}_{sn}}{1+{(w_s{\text{τ}}_{sn})}^2}},\end{array}$$(2)

where ω_s = e_sB/m_sc is the gyrofrequency of species “s”; m_s and e_s denote the mass and charge of species s, respectively; B is the strength of the magnetic field, and c is the speed of light. In Eq. (2), σ_s is the conductivity of individual species s defined as $${\sigma }_s=\frac{n_s{e}_s^2{\text{τ}}_{sn}}{m_s},$$(3)

where n_s is the number density of species s and τ_sn is the mean collision time between species s and neutral particles given by $${\text{τ}}_{sn}=\frac{1}{{\alpha }_s{}_{\text{He}}}\frac{m_s+m_{{\text{H}}_2}}{{\rho }_{{\text{H}}_2}}\frac{1}{{\langle \sigma w\rangle }_{s{\text{H}}_2}}.$$(4)

In Eq. (4), $m_{{\text{H}}_2}$ is the mass of H₂ particles (the dominant collisional partner of charged species), α_sHe is a factor to account for the slowing-down time of species s due to the presence and collisions with He, ${\rho }_{{\text{H}}_2}$ is the mass density of the cloud, and ${\langle \sigma w\rangle }_{s{\text{H}}_2}$ is the momentum transfer rate. We assumed a constant momentum transfer rate for all ions of ${\langle \sigma w\rangle }_{i{\text{H}}_2}=2\times {10}^{-9}{\text{cm}}^3{\text{s}}^{-1}$ and a constant factor to account for the presence of He of α_iHe = 1.14. For electrons we used ${\langle \sigma w\rangle }_{e{\text{H}}_2}=1.3\times {10}^{-9}{\text{cm}}^3{\text{s}}^{-1}$ and α_eHe = 1.16 (see Pinto & Galli 2008).

We note here that the mean collisional time between different ions and H₂ may be different by up to a factor of seven (see Table 1 from Pinto et al. 2008). However, assuming a constant momentum transfer rate for all ions is still a very well justified approximation since for temperatures relevant to molecular clouds, ${\langle \sigma w\rangle }_{i{\text{H}}_2}$ can only vary by up to a factor of two between different ions (see Table 2 from Pinto et al. 2008 and Figs. 3, 5 and 7 from Pinto & Galli 2008).

To calculate the perpendicular, parallel, and Hall resistivities in star formation simulations, the individual number densities each of charged species s, ns , needs to be known throughout the modeled molecular cloud and throughout the evolution of the simulation. In the following sections, we make certain simplifying assumptions to reduce the complexity of the above expressions and calculate the conductivities and the resistivities of a cloud using only a few key variables.

2.1 Parallel conductivity

From Eqs. (2) and (3), the parallel conductivity can be written as $${\sigma }_{\Vert }=\frac{e_s^2}{{\rho }_{{\text{H}}_2}}{\displaystyle \sum _s\frac{1}{{\alpha }_{s\text{He}}{\langle \sigma w\rangle }_{s{\text{H}}_2}}}\frac{m_s+m_{{\text{H}}_2}}{m_s}{n}_s.$$(5)

Due to their low mass, the contribution of electrons to the parallel conductivity is vastly more important than that of the ions (see Figs. 6–9 in Tritsis et al. 2022). Therefore, the parallel conductivity can be very well approximated as $${\sigma }_{\Vert }\approx \frac{e_s^2}{{\rho }_{{\text{H}}_2}}\frac{1}{{\alpha }_{e\text{He}}{\langle \sigma w\rangle }_{e{\text{H}}_2}}\frac{m_e+m_{{\text{H}}_2}}{m_e}{n}_e\approx 5.62\times {10}^{-7}\frac{n_{i,tot}}{{\rho }_{{\text{H}}_2}},$$(6)

where n_i,tot is the total number density of all ions, and due to charge neutrality, it is equal to the number density of the electrons.

2.2 Perpendicular conductivity

To calculate the perpendicular conductivity, we first made the assumption that ω_sτ_sn ≫ 1. This is an excellent approximation for number densities $n_{{\text{H}}_2}\le {10}^9 {\text{cm}}^{-3}$ where both ions and electrons remain very well attached to the magnetic field (see, e.g., Fig. 5 from Tritsis et al. 2022 and Fig. 6 from Tassis & Mouschovias 2007). With that approximation and using Eq. (2), the perpendicular conductivity can be written as $${\sigma }_{\perp }\approx {\displaystyle \sum _s\frac{{\sigma }_s}{{\left({\omega }_s{\text{τ}}_{sn}\right)}^2}}\approx \frac{c^2{\rho }_{{\text{H}}_2}}{B^2}{\displaystyle \sum _s{\alpha }_{s\text{He}}}{\langle \sigma w\rangle }_{s{\text{H}}_2}\frac{m_s}{m_s+m_{{\text{H}}_2}}{n}_s.$$(7)

In contrast to the parallel conductivity, the electrons contribute very little compared to the ions (see again Figs. 6–9 of Tritsis et al. 2022). As such, the perpendicular conductivity can be approximated as $${\sigma }_{\perp }\approx \frac{1}{C_{\perp }}\frac{n_{i,tot}}{v_A^2},$$(8)

where $v_A=B/\sqrt{4\pi \rho }$ is the Alfvén velocity and C_⊥ is a quantity that depends on the chemical properties of a specific cloud. As such, C_⊥ essentially encapsulates a weighted average of the contributions of different ions to the perpendicular resistivity, mathematically expressed as $$\frac{1}{C_{\perp }}\approx \frac{c^2}{4\pi }{\alpha }_{i\text{He}}{\langle \sigma w\rangle }_{i{\text{H}}_2}{\displaystyle \sum _{s\ne e}\frac{m_s}{m_s+m_{{\text{H}}_2}}}\frac{n_s}{n_{i,tot}}.$$(9)

Since the quantity $m_s/\left(m_s+m_{{\text{H}}_2}\right)$ in the sum does not significantly change between the most dominant ion species (0.86 for C⁺, 0.6 for ${\text{H}}_3^{+}$) we can approximate C_⊥ as being dependent only on the physical properties of the cloud and constant throughout the cloud’s evolution.

2.3 Hall conductivity

Similar to the perpendicular conductivity, to derive a simplified expression for the Hall conductivity, we also assumed that ω_sτsn ≫ 1. With that approximation and taking into account charge neutrality, the expression for the Hall conductivity from Eq. (2) reads $${\sigma }_H\approx -{\displaystyle \sum _s\frac{{\sigma }_s}{{\omega }_s{\text{τ}}_{sn}}}\approx -\frac{c}{B}\left(e_e{n}_e+e_i{n}_{i,tot}\right)=0.$$(10)

This is to be expected since |σ_H| ≪ |σ_⊥|, as also supported by numerical calculations for central number densities below n_c ≤ 10⁶–10⁷ cm⁻³ (Nakano et al. 2002; Kunz & Mouschovias 2009; Marchand et al. 2016; Zhao et al. 2018). Therefore, for the range of densities our method considers, the Hall conductivity can be ignored with minimal errors to the perpendicular resistivity (assuming a strict case where the Hall conductivity becomes equal to the perpendicular at a density of 10⁶ cm⁻³, the error produced would be 25% for the last half order of magnitude where we consider our method to be applicable).

Putting the above approximations in Eq. (1), the expressions for the resistivities of the cloud can be written as follows: $$\begin{array}{l}{\eta }_{\Vert }=\frac{1}{{\sigma }_{\Vert }}\approx C_{\parallel }\frac{{\rho }_{{\text{H}}_2}}{n_{i,tot}},{\eta }_{\perp }\approx \frac{1}{{\sigma }_{\perp }}\approx C_{\perp }\frac{v_A^2}{n_{i,tot}},\\ {\eta }_H\approx 0\end{array}$$(11)

where C_‖ = 1.78 × 10⁶g⁻¹s and C_⊥ depends on the specific physical conditions of a simulated cloud. From Eq. (11), it is clear that to calculate the resistivities, we only need to know two quantities: the total number density of all ions n_i,tot and the constant C_⊥. The total ion density is calculated from the total ion abundance χ^i,tot, which is approximated as a function of $n_{{\text{H}}_2}$ by a power law that evolves alongside the cloud depending on its physical conditions. Both the power law and the constant C_⊥ are calculated by interpolations using results from simulations with a detailed chemical network. More details about these calculations are discussed in Appendices A and B.

3 Numerical setup

We followed Tritsis et al. (2022) and Tritsis et al. (2023) to perform a series of 2D nonideal MHD simulations with nonequilibrium chemistry in cylindrical geometry with a radius and half-height of 0.75 pc. Here, we briefly summarize the initial and boundary conditions and refer the reader to Tritsis et al. (2022) for more details. For all models, we used an isothermal equation of state, the magnetic field was initially uniform and aligned with the axis of symmetry (z-axis), and all velocity components were initially set to zero. For our chemodynamical simulations, we used the chemical network without Deuterium chemistry, which was first presented in Tritsis et al. (2022). All of our initial 2D simulations used to develop our interpolation function were performed using an adaptive mesh grid with six levels of refinement and an initial grid size of 32×64 cells in the r and z directions, respectively.

To generalize our method to a large parameter space, we simulated various clouds where we altered the cosmic-ray ionization rate, visual extinction, temperature, and initial H₂ density. In Table 1, we list the physical conditions used in each of these simulated clouds. The mass-to-flux ratio in all of these models was equal to 0.5 in units of the critical value for collapse (Mouschovias & Spitzer 1976). In our fiducial model, the cosmic-ray ionization rate is equal to the standard value of ζ₀ = 1.3 × 10⁻¹⁷ s⁻¹ (Caselli et al. 1998), the visual extinction is equal to 10 mag, the temperature is equal to 10 K, and the initial number density is equal to 300 cm⁻³ .

To test the validity of our approximations and the accuracy of the interpolation function developed for an even wider range of physical conditions, we performed additional simulations where we further altered the physical and initial conditions. However, results from these models never enter the interpolation function. These simulations are listed in the bottom half of Table 1 and also include models with different initial mass-to-flux ratios. In the interest of computational efficiency, these models were performed considering only one level of refinement, as opposed to the six levels used for the models that enter in the interpolation function. For each model listed in Table 1, the resistivities were calculated by using a non-equilibrium chemical network and via our approximation. Therefore, a total of 30 2D simulations were performed in this study.

Finally, we performed a pair of 3D simulations to test the validity of our method in a 3D simulation as well. A thorough description of these simulations will be provided in a future publication. Consequently, in this work we focus on the key differences in comparison to the 2D simulations. The cosmic-ray ionization and visual extinction are not constant, but are calculated for every time step and at each grid point inside the computational domain based on a six-ray approximation. Additionally, we included an initial turbulent power spectrum, with ℳ = σ_υ/c_s = 3 and ℳ_A ≈ 0.45. Finally, the value of the mass- to-flux ratio (in units of the critical value for collapse) was set equal to 0.85.

4 Results

4.1 Time evolution-interpolation models

In Fig. 1, we show the evolution of the central H₂ number density of the modeled clouds as a function of time up to a central number density of 10⁶ cm⁻³. Different lines (described in the legend of the figure) show simulations with different physical conditions (see also Table 1). With black lines, we show our results from the nonideal MHD simulations with non-equilibrium chemistry, thereby representing the “ground truth.” With the red lines, we plot our results using our newly developed approximation for calculating the resistivities. As can be seen from Fig. 1, the simulations where the resistivities are calculated using our approximation are in excellent agreement with the simulations where the resistivities are calculated from a non-equilibrium chemical network for all physical parameters shown. The agreement, both in regard to the time required for the central H₂ number density to reach 10⁶ cm⁻³ and in terms of the overall evolution of the central density throughout time, shows that our approximation can correctly reproduce nonideal MHD effects. In other words, the calculation of the resistivities using our approximation leads to an accurate redistribution of magnetic flux throughout the evolution of the cloud.

Fig. 1 Comparisons of the time evolution of the central H2 number density for models used in the interpolation. The black lines represent simulations using a full chemical network, and the red lines represent simulations using our approximate method. Different line styles are used to denote models with different physical parameters (see legend). The green line represents a model with the same physical parameters as the fiducial model, which also includes the effects of grains. Different physical conditions have a more drastic effect on the evolution of the cloud than the inclusion of grains.

Fig. 2 Same as Fig. 1, but for models where no prior knowledge was present in the interpolating function. The blue lines represent simulations using a full chemical network, and the magenta lines represent simulations using our approximate method.

4.2 Time evolution-generalized models

Figure 1 shows that our method can reproduce results produced by simulations using non-equilibrium chemistry. However, results from these simulations were already present in our interpolating function. To test the applicability of our method to a large parameter space, in Fig. 2 we show a similar comparison of the evolution of the cloud’s central H₂ density as a function of time for simulations with different parameters tested “blindly” (listed in the second half of Table 1). That is, no prior knowledge from these chemodynamical simulations was present in our interpolating function and no further optimization was performed in the interpolating function. Here, we have once again changed various physical parameters along with the mass-to-flux ratio. Despite the fact that the mass-to-flux ratio is not considered as a parameter in the interpolation, our approximation still produces very accurate results. Even in the case of the model with random initial conditions, where every parameter is changed simultaneously, our method can still accurately reproduce the results of a simulation with a full chemical network.

Fig. 3 Magnetic field at the center of the cloud as a function of the central H2 number density for the fiducial and random ICs models. The black and blue lines represent the fiducial and random ICs simulations, respectively, using a full chemical network. The red and purple lines represent the fiducial and random ICs simulations, respectively, using the method described in this paper. One can see that the function approaches a power law B ∝ ρk toward higher densities.

4.3 B ∝ ρ^κ relation

Another test that we conducted to ensure the accuracy of our model was to compare the evolution of the magnetic field at the center of the domain for simulations of various models using non-equilibrium chemistry and the approximate method developed here. In ideal MHD models, the magnetic field evolves with the square root of the density (B ∼ ρ^κ, κ = 0.5; e.g., Tritsis et al. 2015). In nonideal MHD models, the factor κ tends to be slightly lower, depending on how efficiently the magnetic flux can be redistributed. As shown in Fig. 3, the B–ρ relation between simulations using a full chemical network and our model is in very good agreement, for both the fiducial model and the model with random initial conditions.

4.4 Spatial comparison

Going beyond the central point of the cloud, the approximation we developed also needs to preserve the morphology and spatial distribution of the density of a cloud in order for it to be generally applicable to star formation simulations. To test this, as shown in Fig. 4, we made a comparison between a simulation with full non-equilibrium chemical modeling (top row) and a simulation performed using our approximation (bottom row) for our fiducial model. From left to right, the columns show the density distribution in the simulation when the central number density of the cloud is 10⁴, 10⁵ , and 10⁶ cm⁻³ , respectively. The number density of the cloud is shown with the orange color map, and the black lines show magnetic field lines, while the velocity field in the cloud is shown with the blue color map corresponding to the small arrows.

From Fig. 4, it is evident that the simulation where the resistivities are calculated using our approximation is in excellent agreement with the simulation where the resistivities are calculated exactly using the full non-equilibrium chemical network. Specifically, the simulation performed using our approximation yields very accurate results in terms of the morphology and size of the core as well as its kinematic properties.

4.5 Comparisons of 3D simulations

To test how well the approximate method we developed generalizes to 3D simulations, in Fig. 5 we compare the H₂ number density in the central x–y slice of two 3D simulations. In the upper panel of Fig. 5, we present our results from the simulation where the resistivities are calculated using full non-equilibrium chemistry, and in the bottom panel we show the same slice using the approximation described herewith. The results from the chemodynamical simulation presented in the upper panel of Fig. 5 also include the effect of grains in the calculation of the resistivities, which we do not take into account in our approximation. Further discussion about this decision can be found in Sect. 5.2.

As is evident from Fig. 5, our method approximates the behavior of ambipolar diffusion very well, even in a 3D simulation. At a time of 2.5 Myr, the spatial comparison of the two simulations is in very good agreement. The main differences are that in the simulation where we used our approximation, the maximum H₂ density is slightly lower, and the underdense region that can be seen in the lower-right part of both images is more prominent. In general however, even with the presence of turbulence and without accounting for the presence of grains, our approximation can accurately reproduce the results from a “full” nonideal MHD chemodynamical simulation, thus preserving the morphology and fragmentation properties of a cloud.

5 Discussion

5.1 Comparison with other methods

To demonstrate the accuracy of our implementation in comparison to other approximate methods proposed in the literature, we once again show in Fig. 6 the evolution of the central H₂ number density as a function of time for our fiducial model. The black line shows our results from the full chemodynamical simulation, thus representing the “ground truth,” and the red line shows our results from a simulation where the resistivities are calculated using the interpolation function we developed. With the blue lines, we show the results from other approximate expressions for calculating the resistivities.

A common expression used in the literature (e.g., Abe et al. 2024) for the calculation of the ambipolar-diffusion resistivity is that introduced by Shu (1992). In this approximation, the resistivity is calculated as ${\eta }_{AD}=\frac{B^2}{4\pi {\gamma }_{in}{\rho }_n{\rho }_i}$, where the density of ions is calculated as ${\rho }_i=3\times {10}^{-16}{\rho }_{{\text{H}}_2}^{1/2}$ (Elmegreen 1979). The coupling constant between the ions and the neutrals was taken to be γ_in = 3.5 × 10¹³ cm³ g⁻¹ s⁻¹. This relation is thought to be valid up to neutral densities of ~10⁸ cm⁻³ (Shu et al. 1987), that is in approximately the same density range considered here.

Tassis et al. (2012b, blue dash-dotted line in Fig. 6) used results from models of collapsing cores from Tassis et al. (2012a) to arrive at an approximation for the ionization fraction of the form ${\chi }_i=3\times {10}^{-8}{n}_{{\text{H}}_2}^{-0.6}$ (see Fig. 5 in Tassis et al. 2012b). The models used simulated cores of radius 0.4 pc, an initial number density of 1000 cm⁻³, temperatures of T = 7–15 K, mass-to- flux ratios of M/Φ = 0.7–1.3 (M/Φ)₀, and cosmic-ray ionization rates of ζ/ζ₀ = 0.25–4. Therefore, the phase of prestellar-core evolution considered in Tassis et al. (2012b) (and the resulting expression for the ionization fraction) can be directly compared to the setup considered in our study.

Tielens (2005, blue densely dotted line in Fig. 6) considered dense cores $\left(n_{{\text{H}}_2}={10}^3-{10}^8 {\text{cm}}^{-3}\right)$ with A_υ >15, temperatures of T ≈ 10 K, and a cosmic-ray ionization rate of ζ/ζ₀ ≈ 2.3. The ionization fraction was approximated as ${\chi }_i\approx \sqrt{\frac{\zeta }{3\times {10}^{-7}{n}_{{\text{H}}_2}}}$ (see Fig. 10.2 in Tielens 2005). Even though the study by Tielens (2005) considered somewhat more evolved cores compared to our fiducial case to arrive at an expression for the ionization fraction, their resulting functional form leads to higher values for the ionization fraction, and therefore a longer timescale is required for the cloud to collapse.

Tsukamoto et al. (2022, loosely dotted line in Fig. 6) presented three simple power law approximations of η_AD for three different H₂ density ranges. For our comparison to their expression, we only used the derived power law corresponding to an H₂ number density lower than 10⁷ cm⁻³ $\left({\eta }_{AD}\approx 2\times {10}^{18}\sqrt{\frac{{10}^{-16}}{{\rho }_{{\text{H}}_2}}}\right)$.

All the approximations we compare our method to are thus aimed at molecular clouds with physical conditions similar to the ones considered in this paper. For the ones that calculate only the ionization fraction, we calculated the number density of ions and inserted it into Eq. (11). As is evident, most methods fall short in terms of reproducing the dynamical evolution of a cloud by several millions of years, while other methods do not even lead to the collapse of the cloud.

Fig. 4 Comparison of the spatial distribution of the density between a simulation where the resistivities are calculated using our non-equilibrium chemical network (top row) and a simulation where the resistivities are computed using our approximation (bottom row) for our fiducial model. We compare our results for three different central number densities of 104 (left column), 105 (middle column), and 106 cm−3 (right column). The black arrows represent the magnetic field, and the blue arrows represent the velocity of the plasma.

5.2 Effects of grains

In Fig. 1, we include a model that has the same physical parameters as the fiducial model and includes the effects of grains in the calculation of the resistivities. As is evident from Fig. 1, the effect of grains in the overall time evolution of the cloud is dwarfed in comparison to varying the physical conditions of the cloud.

Previous studies have shown that grains play a significant role in the ionization of a cloud and the resulting resistivities. However, the role of grains most likely becomes dominant at densities higher than the ones considered here. Marchand et al. (2023) described the effects of coagulation on the ionization fraction and the resistivities (see their Fig. 7) and demonstrated that these effects become significant for H₂ number densities larger than 10⁶ cm⁻³. Zhao et al. (2018) calculated the resistivities for different grain size distributions. The variation of the ambipolar diffusion resistivity between those different distributions was shown to become significant for H₂ number densities larger than 10⁵ cm⁻³. The perpendicular conductivity was shown to be larger than the Hall conductivity for densities smaller than 10⁶ cm⁻³ regardless of distribution (see their Figs. 6, 7, and 8).

Nakano et al. (2002) found that grains can be the dominant contributor to the perpendicular resistivity for densities larger than 10⁴ cm⁻³. However, in their implementation, they assumed a grain size distribution skewed toward smaller grains. Therefore, the vast majority of the grains remained well attached to the magnetic field at these densities. Lefèvre et al. (2014) and Steinacker et al. (2015) used observations of the coreshine effect to deduce the grain distribution and found that the maximum grain size in starless cores is a factor of approximately five higher than what is used in Nakano et al. (2002) and much closer to the values considered in Kunz & Mouschovias (2009) and Marchand et al. (2016). The range of densities over which our approximation remains accurate would thus need to be lowered should any future observations provide evidence of a smaller grain size distribution.

Regardless, such issues related to the grain size distribution cannot be accurately modeled in any approximate method aiming to model the resistivities. Instead, such complications require a “full treatment” of the resistivities where, apart from the ionization fraction, the properties of the grains can also be handled as free parameters. Additionally, given that clouds, either super or sub-critical, spend the majority of their lifetime at lower densities (i.e., ≤10⁴ cm⁻³), we do not expect that such effects will have a big impact in the overall evolution of a cloud for the range of densities considered here.

Finally, in Sect. 4.5, we presented a comparison between two 3D simulations, one using detailed chemistry and including the effects of grains and another using our method, which does not take grains into account. The two simulations are in close qualitative and quantitative agreement, with only minor differences in their spatial structure. Due to the reasons described above, we are confident that our method can produce reliable results up to densities of ∼10⁶ cm⁻³ for various physical conditions even without considering the effect of grains.

Fig. 5 H2 density on the central x–y slice for a chemodynamical 3D simulation (top) and a simulation using our method (bottom) at 2.5 Myr simulation time. The two simulations seem to have evolved very similarly, showing only minor differences. The simulation in the bottom panel has collapsed at a slightly slower rate, and the spots where the density is lower are exaggerated compared to the simulation in the top panel.

Fig. 6 Central H2 number density evolution in simulations with identical initial conditions (fiducial model; see Table 1) but using different methods for calculating the resistivities. The black line represents a simulation using a full chemical network, the red line represents a simulation using the method described in this paper, and the blue lines represent methods proposed in the literature.

6 Summary

In this paper, we have derived approximate expressions for the resistivities that only depend on a handful of variables (magnetic field, H₂ density, total ion density) and do not require a chemical network to be calculated. Our approximations are very well physically motivated for number densities $n_{{\text{H}}_2}\lesssim {10}^6 {\text{cm}}^{-3}$. These approximations aim to accurately reproduce the resistivities throughout the entire domain, as opposed to only the “central” densest point of a collapsing cloud.

We carried out a series of simulations using these approximate expressions for the resistivities, and we compared them against full nonideal MHD chemodynamical simulations while simultaneously exploring a large part of the parameter space. We find that our approximation leads to excellent results both in terms of the ambipolar-diffusion timescale and in terms of the morphology and shape of the cores formed. At the same time, our approximation leads to a factor of ∼10² increase in computational speed in 2D simulations and ∼10⁴ for 3D simulations.

In contrast to other approximations aiming to address the time step constraint when performing nonideal MHD simulations (i.e., Li et al. 2006), with the method developed here, our aim was to eliminate the need to include a chemical network.

We therefore trust that the method we developed for calculating the resistivities can be proven very valuable to the computational star formation community focusing on the early stages of molecular-cloud evolution and collapsing prestellar cores.

Data availability

The tabulated data required for integrating our method in hydrodynamical codes, along with a FORTRAN implementation of the interpolating function are publicly available here.

Acknowledgements

We thank the referee, whose detailed and insightful comments helped improve our manuscript. E. A. and K. T. acknowledge support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 7712821. A. Tritsis 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. We also acknowledge use of the following software: MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), NUMBA (Lam et al. 2015), and the YT analysis toolkit (Turk et al. 2011).

Appendix A Approximating the evolution of dominant species

Fig. A.1 shows the relation between the total ion abundance and the H₂ number density for every point in the domain at different times throughout the evolution of the cloud. Comparing the abundance-H₂ number density relation for different times, it is evident that the abundance of the dominant species does not follow a single power law with the H₂ density. Instead, the ${\chi }_{i,tot}-{\rho }_{{\text{H}}_2}$ relation evolves as the simulation progresses and converges toward a single relation as the cloud collapses. That single relation, however, gives values of χ_i,tot that are up to two orders of magnitude lower than the true values at the beginning of the simulation. An evolving approximation of the abundance is therefore needed. A simple yet accurate approximation for the abundance, shown in Fig. A.1, is a power law $${\chi }_{i,tot}=A{\left(\frac{{\rho }_{{\text{H}}_2}}{{\rho }_{int}}\right)}^B,$$(A.1)

where the constants A, B are evolving along with the cloud. The method for approximating the abundance is then the following. We use the geometric mean of the maximum density (ρ_max) and logarithmic mean of the H₂ density within the cloud at any time step $\left({\rho }_{\text{int }}=\sqrt{{\rho }_{\max ,{\text{H}}_2}{10}^{\langle {\log }_{10}{\rho }_{{\text{H}}_2}\rangle }}\right)$ to track the evolution of the cloud. This quantity was chosen over simpler ones (e.g., the maximum density of the cloud at any given time), in order to track the cloud’s evolution more accurately under turbulent conditions, since our choice of ρ_int better reflects the entire distribution of densities over the entire cloud. We then calculate the constants A, B as a function of ρ_int for the chemodynamical simulations of the models listed in the first half of Table 1, which we keep in tabulated form. For a given time step with a given value of ρ_int , we can then calculate the values of A, B at all points in the domain using our multivariate interpolation function that also accounts for the physical conditions (e.g., temperature, visual extinction, etc) of the grid point under consideration. The interpolation function is further explained in Appendix. B.

Appendix B General Code Description

The table provided alongside the multivariable interpolation function consists of 17 columns. The first column (IntDens) consists of the geometric mean of the maximum and the logarithmic mean of the H₂ density ρ_int, and the other eight pairs of columns represent the values of A and B corresponding to the values of ρ_int for the models used in the interpolation. The models used are listed in the first half of Table 1 and have varying initial density $n_{{\text{H}}_2,0}$, temperature T, cosmic-ray ray ionization ζ/ζ₀ and visual extinction A_υ. As shown by their names each model focuses on varying one variable compared to the fiducial model. Knowing the values of A, B and ρ_int , the dominant ion density can then be approximated through a power law as in Eq. (A.1).

Given the initial H₂ number density $n_{{\text{H}}_2,0}$, the value of ρ_int at any time step and the values of ζ/ζ₀, A_υ, T of any given point in the computational domain, the constants A, B can be interpolated as follows. Firstly, the density ρ_int is adjusted so that the interpolation can work for different starting densities. The adjusted density is calculated as $${\rho }_{adj}={\rho }_{int}\frac{n_{Fid}/n_{{\text{H}}_2,0}+u^2}{1+u^2},$$(B.1) $$u=\frac{lo{𝑔}_{10}\left({\rho }_{int}/\left(2m_p{n}_{{\text{H}}_2,0}\right)\right)}{3n_{{\text{H}}_2,0}/n_{Fid}}.$$(B.2)

To find the two rows used in the interpolation we need to find the two consecutive rows whose densities bracket the given ρ_int. To that end, we simply do $$\begin{array}{l} foriin [2, size ] : \hfill \\ ifIntDens[i]>{\rho }_{adj } then \hfill \\ N=i\hfill \\ stop. \hfill \end{array}$$(B.3)

The rows used in the interpolation will then be the rows N and N − 1.

Having selected the rows, we then select the columns (i.e., the models) that will be used in the interpolation. The fiducial model is used as one of the models for the interpolation. For ζ/ζ₀, A_υ, T there are two alternative models where each variable differs from the fiducial. For example, for the cosmic-ray ionization rate ζ/ζ₀ we have one model with higher and one model with lower rate than the fiducial model’s value of ζ/ζ₀ = 1. We select which of these two models is more appropriate for the interpolation as follows $$\begin{array}{l}\text{ If }\left(\zeta /{\zeta }_0>=1\right)\text{ then }\\ f_{\zeta }=2\\ A_{\zeta ,N-1}=A_{High\zeta }[N-1]\\ B_{\zeta ,N-1}=B_{High\zeta }[N-1]\\ A_{\zeta ,N}=A_{\text{High}\zeta }[N]\\ B_{\zeta ,N}=B_{High\zeta }[N]\\ C_{\zeta }=C_{\perp ,High}{}_{\zeta }\\ Else\\ f_{\zeta }=0.5\\ A_{\zeta ,N-1}=A_{Low\zeta }[N-1]\\ B_{\zeta ,N-1}=B_{Low\zeta }[N-1]\\ A_{\zeta ,N}=A_{Low\zeta }[N]\\ B_{\zeta ,N}=B_{Low\zeta }[N]\\ C_{\zeta }=C_{\perp ,\text{ Low }\zeta }\\ EndIf\end{array}$$(B.4)

where the subscripts denote the model name and the values of the constant C_⊥ are given in Table B.1. With this process, we have selected four values A_ζ,N−1, B_ζ,N−1, A_ζ,N, B_ζ,N, which correspond to the values of A, B of the appropriate model for the densities bounding ρ_int. We repeat the selection process for the variables A_υ and T. In the models where we varied the temperature we also changed the starting number density. Therefore, the initial number density of the temperature models (n_T) also needs to be defined in the selection process and used in the interpolation. Given that there is only one model focused on varying the starting number density of the cloud $n_{{\text{H}}_2,0}\left(\text{ High }{n}_{{\text{H}}_2}\right)$, no selection process is needed for the starting density.

With the interpolation models chosen, we first calculate the constants A, B corresponding to the two selected rows, N-1 and N for the given physical conditions. The interpolation function for A for the first row (N-1) is of the form $$\begin{array}{l}{A}_{N-1}=\left(A_{Fid,N-1}+\frac{n_{Fid}}{n_{{\text{H}}_2,0}}\left(A_{\zeta ,N-1}-A_{Fid,N-1}\right)\frac{\zeta -1}{f_{\zeta }-1}\right)\\ {\left(\frac{n_{{\text{H}}_2,0}}{n_{Fid}}\right)}^{{\alpha }_{A_{N-1}}}{\left(\frac{T}{T_{Fid}}\right)}^{{\beta }_{A_{N-1}}}{\left(\frac{A_{A_v,N-1}}{A_{Fid,N-1}}\right)}^{{\mathit{exp}}_{A_v}},\end{array}$$(B.5)

where $$\begin{array}{l}{\alpha }_{A_{N-1}}=\frac{lo{g}_{10}\left(A_{n_{{\text{H}}_2,0},N-1}/A_{Fid,N-1}\right)}{lo{g}_{10}\left(f_{n_{{\text{H}}_2,0}}/n_{Fid}\right)}\\ {\beta }_{A_{N-1}}=\frac{lo{g}_{10}\left(A_{T,N-1}/A_{Fid,N-1}\right)-{\alpha }_{A_{N-1}}lo{g}_{10}\left(n_T/n_{Fid}\right)}{lo{g}_{10}\left(f_T/T_{Fid}\right)}\\ ex{p}_{A_v}=\frac{e^{-A_v}-e^{-A_{v_{Fid}}}}{e^{-f_{A_v}}-e^{-A_{v_{Fid}}}}.\end{array}$$(B.6)

Fig. A.1 Evolution of the ion abundance as a function of the H2 number density for the fiducial model at times of, from top right to bottom left, 1 Myr, 3 Myrs, 9 Myrs, and 11.3 Myrs. The black line represents the abundance-density relation at the center of the cloud for the entire duration of the simulation, the red points correspond to the ionization fraction from all cells in the domain at the specific point in time, and the blue dashed lines represent the power law that best fits the abundance-density relation at that point in time.

We repeat that process for the second row (N). To calculate A, we linearly interpolate between A_N−1 and A_N $$\text{A}={\text{A}}_{\text{N}-1}\frac{{\rho }_{\text{int }}-{\rho }_{\text{int,N}-1}}{{\rho }_{\text{int},\text{N}}-{\rho }_{\text{int},\text{N}-1}}+{\text{A}}_{\text{N}}\frac{{\rho }_{\text{int},\text{N}}-{\rho }_{\text{int}}}{{\rho }_{\text{int},\text{N}}-{\rho }_{\text{int},\text{N}-1}}.$$(B.7)

The above process is then repeated for the constant B.

The constant C_⊥ used in calculating the perpendicular resistivity is multi-linearly interpolated from the constants in Table B.1, using the same selection process that we used for A, B. Since C_⊥ is independent of ρ_int, the interpolation is performed using only the initial conditions (ζ/ζ₀, A_υ, T, and $n_{{\text{H}}_2,0}$) as $$\begin{array}{l}{C}_{\perp }=C_{Fid}+\left(C_T-C_{Fid}\right)\frac{T-T_{Fid}}{f_T-T_{Fid}}\\ +\left(C_{A_v}-C_{Fid}\right)\frac{A_v-A_{v_{\text{Fid }}}}{f_{A_v}-A_{v_{\text{Fid }}}}+\left(C_{\zeta }-C_{Fid}\right)\frac{\zeta -{\zeta }_{\text{Fid }}}{f_{\zeta }-{\zeta }_{Fid}}\\ +\left(C_{n_{{\text{H}}_2,0}}-C_{Fid}\right)\frac{n_{{\text{H}}_2,0}-n_{Fid}-\left(n_T-n_{Fid}\right)\left(T-T_{Fid}\right)/\left(f_T-T_{Fid}\right)}{f_{n_{{\text{H}}_2,0}}-n_{Fid}}.\end{array}$$(B.8)

With the values of A, B, and C_⊥ calculated, the ion number density can then be computed as in Appendix A and the resistivities can be determined as in Sect. 2.

References

All Tables

All Figures

Fig. 1 Comparisons of the time evolution of the central H2 number density for models used in the interpolation. The black lines represent simulations using a full chemical network, and the red lines represent simulations using our approximate method. Different line styles are used to denote models with different physical parameters (see legend). The green line represents a model with the same physical parameters as the fiducial model, which also includes the effects of grains. Different physical conditions have a more drastic effect on the evolution of the cloud than the inclusion of grains.

In the text

	Fig. 2 Same as Fig. 1, but for models where no prior knowledge was present in the interpolating function. The blue lines represent simulations using a full chemical network, and the magenta lines represent simulations using our approximate method.
In the text

Fig. 3 Magnetic field at the center of the cloud as a function of the central H2 number density for the fiducial and random ICs models. The black and blue lines represent the fiducial and random ICs simulations, respectively, using a full chemical network. The red and purple lines represent the fiducial and random ICs simulations, respectively, using the method described in this paper. One can see that the function approaches a power law B ∝ ρk toward higher densities.

In the text

Fig. 4 Comparison of the spatial distribution of the density between a simulation where the resistivities are calculated using our non-equilibrium chemical network (top row) and a simulation where the resistivities are computed using our approximation (bottom row) for our fiducial model. We compare our results for three different central number densities of 104 (left column), 105 (middle column), and 106 cm−3 (right column). The black arrows represent the magnetic field, and the blue arrows represent the velocity of the plasma.

In the text

Fig. 5 H2 density on the central x–y slice for a chemodynamical 3D simulation (top) and a simulation using our method (bottom) at 2.5 Myr simulation time. The two simulations seem to have evolved very similarly, showing only minor differences. The simulation in the bottom panel has collapsed at a slightly slower rate, and the spots where the density is lower are exaggerated compared to the simulation in the top panel.

In the text

	Fig. 6 Central H2 number density evolution in simulations with identical initial conditions (fiducial model; see Table 1) but using different methods for calculating the resistivities. The black line represents a simulation using a full chemical network, the red line represents a simulation using the method described in this paper, and the blue lines represent methods proposed in the literature.
In the text

Fig. A.1 Evolution of the ion abundance as a function of the H2 number density for the fiducial model at times of, from top right to bottom left, 1 Myr, 3 Myrs, 9 Myrs, and 11.3 Myrs. The black line represents the abundance-density relation at the center of the cloud for the entire duration of the simulation, the red points correspond to the ionization fraction from all cells in the domain at the specific point in time, and the blue dashed lines represent the power law that best fits the abundance-density relation at that point in time.

In the text