Extreme ionisation rates in star formation

The effect of extreme ionisation rates during the initial collapse of a molecular cloud core

James Wurster, Matthew R. Bate, and Daniel J. Price
School of Physics and Astronomy, University of Exeter, Stocker Rd, Exeter EX4 4QL, UK
Monash Centre for Astrophysics and School of Physics and Astronomy, Monash University, Vic 3800, Australia
Submitted: Revised: Accepted:

What cosmic ray ionisation rate is required such that a non-ideal magnetohydrodynamics (MHD) simulation of a collapsing molecular cloud will follow the same evolutionary path as an ideal MHD simulation or as a purely hydrodynamics simulation? To investigate this question, we perform three-dimensional smoothed particle non-ideal magnetohydrodynamics simulations of the gravitational collapse of rotating, one solar mass, magnetised molecular cloud cores, that include Ohmic resistivity, ambipolar diffusion, and the Hall effect. We assume a uniform grain size of m, and our free parameter is the cosmic ray ionisation rate, . We evolve our models, where possible, until they have produced a first hydrostatic core. Models with s are indistinguishable from ideal MHD models and the evolution of the model with s matches the evolution of the ideal MHD model within one per cent when considering maximum density, magnetic energy, and maximum magnetic field strength as a function of time; these results are independent of . Models with very low ionisation rates ( s) are required to approach hydrodynamical collapse, and even lower ionisation rates may be required for larger . Thus, it is possible to reproduce ideal MHD and purely hydrodynamical collapses using non-ideal MHD given an appropriate cosmic ray ionisation rate. However, realistic cosmic ray ionisation rates approach neither limit, thus non-ideal MHD cannot be neglected in star formation simulations.

magnetic fields — MHD — methods: numerical — stars: formation
pagerange: The effect of extreme ionisation rates during the initial collapse of a molecular cloud coreThe effect of extreme ionisation rates during the initial collapse of a molecular cloud corepubyear: 2018

1 Introduction

Molecular clouds contain magnetic fields (e.g. Crutcher, 1999; Bourke et al., 2001; Heiles & Crutcher, 2005; Troland & Crutcher, 2008) with low ionisation fractions (Mestel & Spitzer, 1956) as low as in dense cores (Nakano & Umebayashi, 1986a; Umebayashi & Nakano, 1990). Prior to star formation, ionisation is mostly driven by cosmic rays interacting with the gas and dust, with contributions from radionuclide decay. After a protostar forms, the protostar itself is thermally ionised and ionises its immediate environment through X-rays. The ionisation rate depends on the source, with typical rates for cosmic rays, radionuclide decay and X-rays given by s (Spitzer & Tomasko, 1968; Umebayashi & Nakano, 1981), s (Umebayashi & Nakano, 2009) and s (e.g. Igea & Glassgold, 1999; Turner & Sano, 2008), respectively, where is the surface density of the gas, and and are the characteristic attenuation depths of cosmic rays and X-rays, respectively.

A completely ionised medium is well represented by ideal magnetohydrodynamics (MHD), while a completely unionised fluid embedded in a magnetic field should be well represented by pure hydrodynamics. In a partially ionised medium, non-ideal MHD is required, where the three non-ideal effects are electron-ion/neutral drift (Ohmic resistivity), ion-electron drift (Hall effect) and ion-neutral drift (ambipolar diffusion). Their relative importance depends, amongst other things, on the gas density, number density of charged species (including grains), gas temperature, and magnetic field strength (e.g. Wardle & Ng, 1999; Nakano et al., 2002; Tassis & Mouschovias, 2007; Wardle, 2007; Pandey & Wardle, 2008; Keith & Wardle, 2014). The Hall effect also depends on the direction of the magnetic field with respect to the rotation axis (e.g. Braiding & Wardle, 2012a, b; Tsukamoto et al., 2015b; Wurster et al., 2016; Tsukamoto et al., 2017).

Many studies have modelled the collapse of a molecular cloud to the first or second Larson core (Larson, 1969) using non-ideal MHD (e.g. Nakano & Umebayashi, 1986b; Fiedler & Mouschovias, 1993; Ciolek & Mouschovias, 1994; Li & Shu, 1996; Mouschovias, 1996; Mouschovias & Ciolek, 1999; Shu et al., 2006; Mellon & Li, 2009; Duffin & Pudritz, 2009; Dapp & Basu, 2010; Machida et al., 2011; Li et al., 2011; Krasnopolsky et al., 2012; Dapp et al., 2012; Tomida et al., 2013; Tomida et al., 2015; Tsukamoto et al., 2015a, b; Wurster et al., 2016; Tsukamoto et al., 2017). For efficiency, these studies typically assumed that the dominant ionisation source at early times was cosmic ray ionisation and that there was no attenuation. Thus, the canonically used cosmic ray ionisation rate is s.

The first three-dimensional models of collapsing magnetised molecular clouds were performed using ideal MHD (e.g. Price & Bate, 2007; Hennebelle & Fromang, 2008; Duffin & Pudritz, 2009; Hennebelle & Ciardi, 2009; Commerçon et al., 2010; Seifried et al., 2011), despite ideal MHD being a poor approximation to observed molecular cloud environment. However, these studies provided useful insight into the behaviour of magnetic fields and provided important benchmarks for future studies. These models were effectively fully ionised, thus what cosmic ray ionisation rate would be required to reproduce these results, assuming non-ideal effects MHD were included?

As a rotating molecular cloud collapses, a dense disc forms (e.g. Larson, 1972; Tscharnuter, 1987). Given a realistic cosmic ray attenuation rate, the centre of the dense disc should be very weakly ionised or completely neutral, forming a magnetic dead zone (Gammie, 1996). A very weakly ionised medium can be self-consistently modelled with non-ideal MHD, however this can be very expensive to run. Thus, at what ionisation rate can a medium be treated as purely hydrodynamical?

The goal of this study is to model the early collapse of a rotating, magnetised molecular cloud core using non-ideal MHD to determine at what cosmic ray ionisation rates (if any) a purely hydrodynamical or an ideal MHD collapse can be recovered. The free parameter is the cosmic ray ionisation rate, , which we held constant throughout each simulation. Due to the computational expense when low ionisation rates are used, we only model the collapse up to the formation the first hydrostatic core, except in our two lowest ionisation rate models, which never evolved out of the isothermal collapse phase. In Wurster et al. (2018), we examined how the collapse to stellar core formation changes if one assumes cosmic ray ionisation rates s.

This paper is organised as follows: In Section 2, we present our numerical methods, and in Section 3 we present our initial conditions and discuss how the initial environment is affected by different cosmic ray ionisation rates. In Section 4 we present and discuss our results, and we conclude in Section 5.

2 Numerical method

2.1 Non-ideal magnetohydrodynamics

We solve the equations of self-gravitating, non-ideal magnetohydrodynamics given by


where is the Lagrangian derivative, is the density, is the velocity, the hydrodynamic pressure, is the magnetic field (which has been normalised such that the Alfvén velocity is defined as in code units; see Price & Monaghan 2004), is the gravitational potential, is the gravitational constant, and is the identity matrix. The equation set is closed by a barotropic equation of state,


where is the initial isothermal sound speed, and the density thresholds are g cm and g cm.

The non-ideal MHD term in (3) is


where the non-ideal coefficients for Ohmic resistivity, the Hall effect, and ambipolar diffusion terms are given in (e.g.) Wardle & Ng (1999); Wardle (2007).

To calculate the number densities of the charged species and thus the non-ideal MHD coefficients, we use Version 1.2.1 of the Nicil library (Wurster, 2016). The maximum temperature reached in this study will be  K, thus cosmic rays will be the only ionisation source, since we intentionally ignore ionisation from radionuclide decay in order to test the effect of low ionisation rates. Cosmic rays can create two species of negatively charged ions: a light ion species based upon hydrogen and helium components and a heavy ion species with the mass of magnesium (e.g. Asplund et al., 2009). We include three species of grains which can absorb free electrons to become negatively charged , or lose electrons through collisions to become positively charged , or remain neutral . The total number density of grains is dependent on the local gas density, viz.,


where is the gas to dust ratio, and are the masses of a neutral particle and dust grain, respectively, and is the gas number density. To conserve gain number density, .

2.2 Smoothed particle non-ideal magnetohydrodynamics

Our calculations are carried out using the 3D smoothed particle magnetohydrodynamics (SPMHD) code Phantom (Price et al., 2017) with the inclusion of self-gravity and non-ideal MHD (Wurster, 2016). The density of each SPH particle is calculated by iteratively solving


using the Newton-Raphson method, where we sum over all neighbours , and are the particle’s mass and smoothing length, respectively, is the smoothing kernel, and is a coefficient required to obtain 58 neighbours when using the adopted cubic spline kernel.

The remainder of the discretised SPMHD equations are readily available in the literature (e.g. see review by Price, 2012), and we use the same form as given in Wurster et al. (2016). We enforce the divergence-free condition on the magnetic field using the constrained hyperbolic/parabolic divergence cleaning algorithm described in Tricco & Price (2012) and Tricco, Price & Bate (2016).

In ideal MHD, artificial resistivity is required for magnetic stability (i.e. the final term in Eqn. 3); as per convention, artificial resistivity is included in all of our simulations both for consistency and for the possibility that physical and artificial resistivity are important in different regions. We use the form given by Price & Monaghan (2004, 2005), however, the signal velocity is instead given by (Price et al., 2017). A comparison of artificial resistivity algorithms presented in Wurster et al. (2017a) showed that the method used here is the least dissipative of all SPMHD algorithms used to date, especially during the collapse to form the first core. Hence, our results are dominated by physical and not artificial resistivity.

2.3 Timestepping

In ideal MHD, the limiting timestep for particle is typically the Courant-Friedrichs-Lewy (CFL) condition,


where is the dimensionless Courant number and is the sound speed. However, non-ideal effects each add a new time-constraint, viz.,


where is a dimensionless coefficient analogous to the Courant number. Test cases with ambipolar diffusion show that the non-ideal MHD timestep can be 40-50 shorter than the CFL timestep (e.g. Mac Low et al., 1995; Wurster et al., 2014), however, in realistic problems, the minimum non-ideal MHD timestep can be several hundred times shorter in quickly evolving, dense regions.

Super-timestepping (Alexiades et al., 1996) is used to relax the conditions given by (12) for the diffusive terms that are parabolic in nature (i.e. Ohmic resistivity and ambipolar diffusion). This involves taking steps of d where and and requiring stability only at the end of steps rather than at the end of every step. The best possible speed-up yields (e.g. Choi et al., 2009),


where we set . For added stability, we first sub-divide by positive integer such that and then take steps per ; this subdivision by is only required in extreme environments where multiple physical processes are simultaneously contributing to a complex evolution, or for very low ionisation rates (i.e. s).

Given the hyperbolic nature of the Hall effect, we are required to solve this timestep explicitly. Methods (e.g. O’Sullivan & Downes, 2007; Meyer et al., 2012) have been proposed to sub-step with the Hall term, but these have yet to be implemented into Phantom.

3 Initial conditions

Our models are similar to those used in our previous studies (Price & Bate, 2007; Bate et al., 2014; Wurster et al., 2016) and consists of a spherical cloud of radius  cm = 0.013 pc and density g cm which is placed inside a low-density box of edge length and a density contrast of 30:1; the cloud and surrounding medium are in pressure equilibrium. This allows the cloud to be modelled self-consistently, and we use quasi-periodic boundary conditions at the edge of the box, in which SPH particles interact magnetohydrodynamically ‘across the box’, but not gravitationally. Our simulations use particles in the sphere which are initialised on a regular close-packed lattice.

The initial cloud has mass  M, rotational velocity rad s, and sound speed  cm s (i.e.  K). We thread the cloud with a uniform magnetic field that is anti-aligned with the axis of rotation, i.e. , which will promote disc formation in the presence of the Hall effect (e.g. Braiding & Wardle, 2012a, b; Tsukamoto et al., 2015b; Wurster et al., 2016); this configuration yields . The magnetic field has an initial strength of G, which corresponds to a normalised mass-to-flux ratio of , where is the initial mass-to-flux ratio and is the critical value where magnetic fields prevent gravitational collapse altogether; is the total mass contained within the cloud, is the magnetic flux threading the surface of the (spherical) cloud at radius assuming a uniform magnetic field of strength , is the gravitational constant and is a parameter numerically determined by Mouschovias & Spitzer (1976). The free-fall time is  yr, which is the characteristic timescale for this study.

The non-ideal MHD models use the default values included in the Nicil library (Wurster, 2016). The dust grains have a radius and bulk density of m and g cm (Pollack et al., 1994), respectively, and the dust-to-gas ratio is . The mass of the neutral particle is based upon the hydrogen and helium abundance, thus m, where m is the mass of a proton. We test 15 cosmic ray ionisation rates in the range s, which are indicated by the vertical lines in Fig. 1. Our models will be named after their ionisation rate such that model has a cosmic ray ionisation rate of s. The ideal MHD and purely hydrodynamical models will be referred to as iMHD and HD, respectively.

Since the time constraint imposed by non-ideal MHD quickly becomes prohibitively expensive for low ionisation rates, our focus is on comparing the very early phases of the collapses, typically just prior to or just after entering the first hydrostatic core phase, which begins at g cm. The maximum density analysed in this study is g cm, however, several of the low ionisation rate models end at g cm.

3.1 Initial behaviour of the charged species and non-ideal MHD coefficients

Fig. 1 shows the species number densities and the non-ideal MHD coefficients, calculated using our initial conditions.

Figure 1: The species number densities (top) and non-ideal MHD coefficients (bottom) calculated using our initial conditions of g cm, G (i.e. ) and  cm s (i.e. 13.5 K). The vertical lines represent the values of the cosmic ray ionisation rate, , that are included in our suite. The grain populations are the dominant species at low ionisation rates, and the positively charged ions and electrons are dominant at high ionisation rates. This turnover is reflected in , where . Increasing  does not lead to a monotonic change in , thus models with different initial values of  will start with different non-ideal effects controlling the evolution.

At constant density, temperature and magnetic field strength, the number densities and non-ideal MHD coefficients are strongly dependent on the cosmic ionisation rate, .

At ionisation rates of s, cosmic rays are unable to ionise ions rapidly enough for the ions and electrons to significantly contribute to the charged species populations; at these rates, the charged species are from grain collisions that transfer electrons to make a positively and negatively charged grain population, with . At ionisation rates of s, the ion and electron populations are several orders of magnitude more populous than the charged grain number densities. However, the grains have a much larger mass than the ions (i.e. m compared to m and m), thus still contribute non-trivially to the value of the non-ideal MHD coefficients even at very high ionisation rates.

For s (recall that the canonical cosmic ionisation rate is s), the number density of ions and electrons is similar, and the ionisation fraction reaches at s; thus, even at high cosmic ray ionisation rates, thermal ionisation or another source is required to fully ionise the medium.

Using our given initial conditions, the initial non-ideal MHD coefficients have six regimes:

  1. : ;

  2. : ;

  3. : ;

  4. : ;

  5. : ;

  6. :                     .

Although ambipolar diffusion is typically the dominant effect, the Hall effect is the dominant term in region (iv), which is the same region where grains transition from higher number densities compared to the ions to lower number densities. At very low ionisation rates, the Ohmic coefficient is approximately constant, since it is dependent on the Ohmic conductivity , which is approximately constant due to the grain number densities. The Hall coefficient rapidly decreases at low ionisation rates due to its dependence on the Hall conductivity , which will rapidly decreases for . The ambipolar coefficient is also dependent on the Hall conductivity, but via the perpendicular conductivity, hence its delayed decrease. At high ionisation rates, all three terms decrease rapidly as the cloud becomes more ionised. Thus, all three terms have less of an effect on the evolution of the cloud in an absolute sense.

These results are qualitatively similar to those that are obtained from using different , and . Therefore, the numerical results will necessarily differ if we change our initial conditions, however, our qualitative results will be independent of them.

3.2 Grain properties

Although a uniform grain size is not realistic, they are common in numerical models. We use the uniform grain size of m to match our previous studies (Wurster et al., 2016, 2017b, 2018) and to agree with the fiducial value suggested by Pollack et al. (1994). Uniform grain sizes of smaller radius were used in Tsukamoto et al. (2015a, b, 2017).

An alternative to the uniform grain size is the Mathis, Rumpl & Nordsieck (1977) (MRN) grain distribution,


where is the number density of the hydrogen nucleus, is the number density of grains with a radius smaller than , and cm (Draine & Lee, 1984). In this distribution, there are more grains with smaller radii, thus the smaller grains will more strongly influence the evolution than the larger grains. Fig. 2 shows the non-ideal MHD coefficients as a function of  using a uniform grain size of and m (top two panels), and using the MRN grain distribution using the ranges suggested in Kunz & Mouschovias (2009) and Wardle & Ng (1999) (bottom two panels, respectively).

Figure 2: The non-ideal MHD coefficients using using our initial conditions as in the bottom panel of Fig. 1. The dotted lines in each panel are calculated using our fiducial uniform grain size of m, and the solid lines are calculated using the grain size/distribution listed in the panel. The MRN distribution is calculated assuming 40 bins of equal width in log-space. The coefficients have a greater dependence on grain properties at lower cosmic ray ionisation rates, where the coefficients can differ by up to 9 dex; at high ionisation rates ( s) the coefficients typically differ by less than 10 per cent.

At high ionisation rates ( s), the coefficients differ by less than 10 per cent, except for m where the 10 per cent agreement is only for s. Thus, for high ionisation rates, we expect the grain properties to play a minimal role in the evolution of the system. As  decreases to realistic rates ( s), the coefficients become more dependent on the grain properties, although the coefficients for our fiducial grain size and the MRN distribution using the Kunz & Mouschovias (2009) range differ by less than a factor of 1.3.

At low ionisation rates ( s), the coefficients can differ by up to nine orders of magnitude. For m, the coefficients are larger than using our fiducial grain size, suggesting that these systems will approach the hydrodynamical limit at higher ionisation rates due to greater non-ideal MHD effects. For m and the MRN distributions, the coefficients are typically lower than for our fiducial grain size suggesting these models will be slightly more ideal than our fiducial models. Given that the MRN distributions are the more realistic models, the results we present in Section 4.4 will be upper limits such that using the MRN distribution would require even lower ionisation rates to approximate the hydrodynamical case than models using our fiducial uniform grain size of m.

4 Results

4.1 Limiting models: Hydrodynamic and Ideal MHD

The iMHD and HD models represent the limiting models such that the models are totally ionised and totally neutral, respectively. Fig. 3 shows their evolution of the total kinetic energy and the maximum density.

Figure 3: The evolution of the total kinetic energy (top) and the maximum density (bottom) for the purely hydrodynamic (HD) and ideal MHD (iMHD) models. The kinetic energy differs by more than 10 per cent for , and the maximum density differs by more than 10 per cent at .

Qualitatively, both models follow similar trends, with a slow evolution until , at which time the collapse occurs very rapidly. The presence of strong magnetic fields delays the collapse of the molecular cloud, such that, by , the maximum densities differs by 10 per cent, and it then takes iMHD 370 yr longer to reach g cm than HD.

The kinetic energy begins to diverge almost immediately with it differing by 10 per cent after , with more kinetic energy in HD than in iMHD at any given time. This is expected since the magnetic field supports the cloud against gravitational collapse, thus the gas collapses slightly slower.

4.2 Transition models: Non-ideal MHD

We define any model that includes non-ideal MHD as a transition model since the ionisation fraction is neither one (iMHD) nor zero (HD). Fig. 4 shows the late evolution () of the maximum density for selected models, and Figs. 5-8 show cross sections of the density, magnetic field strength, radial and azimuthal velocities of selected models at , and g cm.

Figure 4: The evolution of the maximum density for selected models for . The thick green and cyan lines represent models iMHD and HD, respectively, whose full evolution is shown in the bottom panel of Fig. 3. Models with s all lie on top of the iMHD curve, thus only two have been shown for clarity. Model evolves slightly faster than iMHD, while - evolve slower. The low ionisation rate models with s evolve similarly to HD. The maximum densities for HD and - do not coincide with the centre of the core. Not all models have reached g cm due to computational limitations. The horizontal lines match the maximum densities shown in Figs. 5-8.
Figure 5: Density slice through the core of the clouds of selected models at , and g cm. Frame sizes and colour bar range change with each row to better show the structure at each . At g cm, a dense column of gas has formed along the rotation axis in the higher ionisation rate models ( s), an oblate spheroid has formed for mid-range ionisation rates, and the collapse is approximately spherical for low ionisation rates ( s); the maximum density is not in the core for the models with s. At g cm, the scale height of the discs decreases from iMHD to and then increases again.
Figure 6: Magnetic field strength through the core of the clouds of selected models as in Fig. 5. At g cm, the dense column of gas has an enhanced magnetic field strength in the higher ionisation rate models ( s), while the magnetic field strength is approximately uniform at low ionisation rates ( s). By g cm, the models with s have an unstructured magnetic field that is stronger in the midplane.
Figure 7: Radial velocity slices through the core of the clouds of selected models as in Fig. 5. The initial radial infall (i.e. where ) is approximately spherical for the low ionisation rate models ( s) and faster at larger radii from the core than for models with higher rates ( s). At larger , the infall becomes less spherical in all models, and the infall rate is faster in the midplane for the magnetised models while the vertical infall is faster in HD.
Figure 8: Azimuthal velocity slices through the core of the clouds of selected models as in Fig. 5. The rotation speed of the cloud increases as the cosmic ray ionisation rate decreases and as time increases. At g cm, the rotational speed is faster above and below the disc for s; a counter-rotating envelope is beginning to form in and .

The models with s evolve similarly to iMHD and will be discussed in more detail in Section 4.3 below. These clouds are magnetically supported, thus form a dense collimated structure with strong magnetic fields, since initially, the gas is free to collapse along the rotation axis. The models with s follow a similar evolutionary path as HD and will be discussed in Section 4.4. These models follow an initially spherical collapse, however, by g cm, a thick, rotationally supported disc has formed; throughout their evolution, they retain an approximately uniform magnetic field.

Models to do not represent a smooth transition between the evolutionary paths of iMHD and HD as the initial ionisation rate is decreased, with the exception of the magnetic field strength and geometry. Agreeing with intuition, evolves similarly to but slightly faster than iMHD, thus is in the region bracketed by iMHD and HD in Fig. 4. As expected, at our selected  snapshots, its central gas distribution is more diffuse and has a weaker magnetic field strength than the models with s.

The models with evolve slower than iMHD, with having the longest evolutionary time. This is a result of the Hall effect. At these early times, the rotation is beginning to convert the poloidal magnetic field, , into a toroidal magnetic field, . Since the field lines remain closed, there are both components, and the Hall effect enhances one component and decreases the other. At this stage, the initial direction of the magnetic field is relatively unimportant for the characteristics that we investigate, and we find that models initialised with collapse only slightly faster than their counterparts with , but still do not collapse faster than iMHD. By g cm, a weak counter-rotating envelope has formed in and , and we have verified that this does not form for the models.

Fig. 9 shows the maximum magnetic field strength and the magnitude of the maximum toroidal field as a function of maximum density , which we use as a proxy for time; recall that .

Figure 9: The evolution of the maximum magnetic field strength (top) and the magnitude of the maximum toroidal field (bottom) for selected models. The thick cyan line represents the initial magnetic field strength. The legend is split across the panels for clarity. Initially, the entire magnetic field is poloidal , and the evolution of the poloidal and total magnetic field strengths are similar. The maximum magnetic field strength increases for all models, although only increases by a factor of 1.7 for the low ionisation models. For the majority of the models, , however, , , , , and all switch to having a dominant toroidal field at increasing maximum densities.

Although the absolute value of the Hall effect is not the strongest in (i.e. there are larger values of at lower ionisation rates; see Fig. 1), the net magnetic field in this model is stronger than in (i.e. the model with the largest ) since that model has large ambipolar diffusion. Thus, more of the net magnetic field is converted into the toroidal magnetic field in than (or any other model in our suite). In several models, the maximum toroidal field becomes stronger than the maximum poloidal field at a given location, with occurring for , , , , and at increasing maximum densities.

The maximum magnetic field strength continually increases as the cloud collapses for the models with s, although the growth rate is much slower for and . By g cm, , and for , and iMHD, respectively; the low ionisation rate models (i.e. , and ) have magnetic field strengths that asymptote at .

Fig. 10 shows the evolution of the total kinetic and magnetic energies for selected models. These values include the gas in both the cloud and background medium.

Figure 10: Evolution of the total kinetic (top) and magnetic (bottom) energies for selected models. The thick green and cyan lines represent iMHD and HD, respectively; the thick cyan line in the bottom panel represents the initial magnetic energy, . The legend is split across the panels for clarity. When comparing kinetic energy at similar maximum densities, there is typically a progression from the iMHD to HD models as the ionisation rate decreases. For s, the total magnetic energy increases as the molecular cloud collapses, whereas it slightly decreases for s.

Magnetic fields support the molecular cloud against collapse, thus at any given , there is more kinetic energy in HD than in iMHD, with the kinetic energy decreasing from the HD value to the iMHD value as  is increased; this can also be seen in Figs. 7 and 8, which show decreasing radial and azimuthal velocities for increasing . The total magnetic energy increases for models with s, with the magnetic energy growing more slowly for models with lower . The total magnetic energy decreases for the models with s; at the end of the simulation, the models with s have decreased by 2 per cent, while has decreased by only 0.4 per cent.

4.3 Approaching ideal MHD

As ionisation rates increase, the medium becomes more ionised, thus the ionisation fractions begin to approach the ideal MHD limit; however, at temperatures and densities presented here, ionisation fraction remains . The question we ask in this section is how high of an ionisation rate is required to safely approximate the ideal MHD limit?

From the induction equation (Eqn. 3), a one-fluid non-ideal MHD system can be approximated as ideal when . However, we want to know the largest possible that will still result in a system equivalent to iMHD. Since different quantities may show equivalence at different , we will compare several different quantities below, including 2D cross sections, total and maximum values. A non-ideal MHD model must be equivalent in all of our quantities to be deemed equivalent to iMHD.

Upon visual inspection of the cross section plots in Figs. 58, we find that the models with s look similar to iMHD for all properties. This also appears true for the evolution of the maximum density with time (Fig. 4) and for total kinetic and magnetic energy (Fig. 10). To quantify this, Fig. 11 shows the relative difference,


where and are the values for the non-ideal model and iMHD, respectively, of the maximum density, , and Fig. 12 shows the relative difference of magnetic energy and total magnetic field strength of a high ionisation rate model compared to the iMHD model; for comparison, we intentionally include models that clearly do not approach the ideal limit.

Figure 11: The relative difference, as defined in (15), of the maximum density with respect to time for the high ionisation rate models compared to iMHD. The relative differences are less than 10 per cent for , however, even slight differences in evolution times at high densities (i.e. at ) cause large relative differences. The evolution times of the models with s agree within

for the majority of the collapse.

Figure 12: The relative differences of the magnetic energy (top) and the maximum magnetic field strength (bottom) of the high ionisation rate models compared to iMHD. For g cm, the magnetic energy of the high ionisation rate models differ by less than 10 per cent from iMHD, and the maximum magnetic field strength differs by less than 40 per cent.

The maximum density in the models with s differs from iMHD by less than 10 per cent for the entire calculation, with only yielding a difference larger than one per cent near the end.

The relative difference of the magnetic energy and maximum magnetic field strength compared to decreases for increasing ionisation rates; models with higher ionisation rates are more similar to iMHD. The magnetic energy of all models with s differs from iMHD by less than 10 per cent, while the energy in models with s differ by less than one per cent.

Total magnetic energy is a global property of the simulation, including both the quickly collapsing inner region and the slowly evolving background medium. The maximum magnetic field strength, however, is localised in or near the core, and is sensitive to the evolution of the collapse and is obtained from a single particle; thus, we do not expect as low of relative difference as in the energies. Models with s have a maximum magnetic field strengths that differ from iMHD by less than , while models with s differ by less than one per cent. The relative difference between and iMHD increases to almost 10 per cent in the first hydrostatic core; this large difference and the increasingly large difference in RD() suggests that does not approach the ideal MHD limit.

For s, the grain properties have minimal effect on the non-ideal MHD coefficients (see Section 3.2), except when decreasing to smaller grains of uniform size. Although our conclusions will be qualitatively unaffected by switching to larger grains or to an MRN grain distribution, the agreement between (e.g.) and iMHD may not be as robust for the smaller grain size of m.

For the collapse up to g cm, we conclude that non-ideal MHD models with s agree with ideal MHD within one per cent, and models with s agree with ideal MHD within of 0.1 per cent, suggesting that they are essentially indistinguishable from ideal MHD models.

4.4 Approaching pure hydrodynamics

As ionisation rates decrease, the medium becomes more neutral, thus begins to approach the purely hydrodynamic limit. Unlike approaching the ideal limit where all simulations include magnetic fields and non-ideal effects begin to negligibly contribute to their evolution, the hydrodynamic models by definition exclude magnetic fields; thus, for a low-ionisation rate model to approach the HD model, their evolution must essentially ignore the magnetic field.

In order that the one-fluid non-ideal MHD equations reduce to the ideal MHD limit, all that has to happen is for the third term in Eqn. 3 to become negligible. Thus, it is unsurprising that with a high enough ionisation rate, the ideal MHD limit is recovered to a high level of accuracy (see Section 4.3). To recover the hydrodynamical limit, however, all four terms in the induction equation (Eqn. 3), which are calculated using the magnetic field and density of a particle and its neighbours must sum to exactly zero, and the second term in the momentum equation (Eqn. 2) must also be zero. Thus, we do not expect to be able to approximate the hydrodynamic limit using non-ideal MHD to the same level of accuracy that was achieved for the ideal limit.

As can be seen in the cross section plots (Figs. 58), shares more characteristics with HD than iMHD or even ; however, there are still noticeable structure differences. Model HD has a slightly larger scale height, and slightly faster rotational and radial velocities. Based upon visual inspectional alone, we conclude that by decreasing the ionisation rate, the models approach HD, however, is not equivalent to HD. At g cm, , and are indistinguishable from one another, although, neither of the latter two models progressed to g cm for a later comparison.

For our quantitative analysis, Fig. 13 shows  and the relative difference in  of the low ionisation rate models compared HD.

Figure 13: Top: Evolution of the maximum density over a short range of time for our low ionisation rate models, as in Fig. 4. Bottom: The relative difference of the maximum density with respect to time over the entire evolution, as in Fig. 11. The collapse rates are slightly faster for and compared to HD, with the relative difference never surpassing 10 per cent. Model collapses slightly slower than HD, and the relative difference surpasses 10 per cent at .

The evolution of has been included for reference, but it clearly cannot approximate HD. The collapse rate for these low ionisation rate models varies in time with respect to HD, collapsing slightly faster or slower depending on the time; the transitions are marked by the relative difference approaching zero in these plots. Over the entire evolution, collapses slower than HD, with the relative difference surpassing 10 per cent by . Given the slower collapse and the visual differences, is not equivalent to HD. Models and collapse faster than HD but evolve at the same rate as one another, and at any given time, their maximum density differs from HD by less than 10 per cent.

In , and , the maximum magnetic field strength grows by a factor of 1.7 (bottom panel of Fig. 9); for reference, by this density of g cm, the magnetic field strength in iMHD has increased by a factor of 60. Despite the maximum magnetic field strength slightly increasing, the total magnetic energy in these three models decreases to (bottom panel of Fig. 10).

If the grain size was decreased to m or an MRN distribution was used, then the non-ideal MHD coefficients would be lower thus these models would be slightly more ideal (Fig. 2). Therefore, modelling smaller grains would yield low ionisation rate models that are even less similar to the hydrodynamical case than those presented here. If an m grain size was modelled, then the non-ideal MHD coefficients would be larger, suggesting these models would be more similar to HD. However, the larger coefficients results in a smaller timestep (recall Eqn. 12), which makes it prohibitively expensive to model; see further discussion in Section 4.5.

Thus, even with low cosmic ray ionisation rates of s, the relative differences of the energies, magnetic field strengths and the maximum densities remain between one and 10 per cent. However, given the slow growth of the magnetic energy and maximum magnetic field strength, we can conclude that models with s approach the hydrodynamic model. At these low ionisation rate, unless precision results are required, pure hydrodynamics can be used in place of non-ideal MHD.

4.5 Timestepping

Given unlimited resources and time, one should simply use non-ideal MHD with the desired cosmic ray ionisation rate rather than ideal or hydrodynamic approximations. However, including non-ideal MHD adds to the computational expense both by additional calculations and by decreasing the minimum timestep.

Fig. 14 shows the minimum and CFL timesteps at each density.

Figure 14: The evolution of the minimum timestep as calculated from the non-ideal terms (see Eqn. 12); the thick green line is the timestep calculated from the Courant-Friedrichs-Lewy (CFL; see Eqn. 11) condition, and is included for reference. For the models whose timesteps are controlled by Ohmic resistivity or ambipolar diffusion, we plot d, where is given in (13). Initially, requires the smallest timesteps, while uses a timestep times larger due to super-timestepping with the Ohmic time constraint. Models with s have d thus are excluded.

All model with s have d, while has d. Models with are always limited by the Hall timestep, while and are Hall-limited for g cm and g cm, respectively. Model is limited by the super-timestepping algorithm using the Ohmic timestep. Note that the ambipolar or Ohmic timesteps would be the limiting case in most models if super-timestepping were not used.

At g cm, the non-ideal timestep in is limited by the Hall timestep and is 3900 times shorter than the CFL timestep. At g cm, the non-ideal timestep in (which uses the fiducial cosmic ionisation ray rate) is 15 times shorter than the CFL timestep. Given that studies of collapsing molecular clouds need to reach maximum densities at least a few orders of magnitude larger than presented here, the slow-down of a factor of a few for the canonical ionisation rate of s can be tedious, while it can be completely prohibitive for lower ionisation rates.

The effect of non-ideal MHD on performance is very evident in the number of CPU hours used, as shown in Fig. 15.

Figure 15: The cumulative CPU time used for each model. HD, iMHD and models with s run very quickly to g cm. As the ionisation rate decreases, the required resources becomes prohibitively expensive, especially for s. To reach g cm, it takes 2200 times longer than iMHD. All the non-ideal models with s take the same length of time, which is 1.4 times longer than iMHD to reach g cm.

The HD model runs 53 per cent faster than the iMHD model in part due to the reduce number of calculations required and the slightly longer CFL timestep since . The models with s take 1.4 times longer to reach g cm than iMHD; since their non-ideal timesteps are similar to the CFL timestep, the additional time used is mostly taken up by the Nicil library. Given that the models with s well approximate the ideal MHD model (see Section 4.3), an ideal MHD model can instead be run at a speed-up of 1.4.

As the ionisation rate decreases, the cumulative CPU time increases, with being the most expensive simulation followed by and then the remainder of the models in order of increasing ionisation rate. To reach g cm, it takes 2200 times longer than iMHD, thus these models are prohibitively expensive to run to any useful maximum density. However, since the models with s have small relative differences for their total magnetic energies and maximum magnetic field strengths, respectively, it would be reasonable to run purely hydrodynamical models in their place.

5 Summary and conclusion

We have presented a suite of non-ideal magnetohydrodynamics simulations with various cosmic ray ionisation rates  to determine what rate is required to recover a hydrodynamical collapse and an ideal MHD collapse. Our models were initialised as a 1 M, spherically symmetric, rotating molecular cloud core; the cloud was magnetised with a magnetic field initially aligned anti-parallel to the rotation axis, and had an initial strength of G, or . Our models used a uniform grain size of m, but we discussed how different grain models are expected to change the results. In particular, we found that at high cosmic ray ionisation rates the results will be approximately independent of grain properties, but at very low ionisation rates the non-ideal MHD coefficients that we uses tend to be larger than those produced by MRN grain size distributions, suggesting that the ionisation rates required to approximate hydrodynamical evolution may be even lower than the values we find using m. All three non-ideal MHD terms (Ohmic resistivity, Hall effect and ambipolar diffusion) were included. We tested 15 different cosmic ray ionisation rates, which were held constant for the entire simulation. The initial density of the molecular cloud core was g cm, and it was evolved until g cm, depending on the value of .

Our two key results are as follows:

  1. Approaching the ideal MHD limit: We evolved models with high ionisation rates until they entered the first hydrostatic core and reached g cm. Models with s were indistinguishable from the ideal MHD model when considering the evolution of their maximum density, magnetic energy and magnetic field strengths i.e. all properties matched the evolution of the ideal MHD model within 0.1 per cent. The evolution of the model with s matched the evolution of the ideal MHD model within 1 per cent.

  2. Approaching the hydrodynamic limit: The models with s look similar to the hydrodynamical model, but with noticeable differences; in these models, the total magnetic energy grew by a factor of 1.7 and the maximum magnetic field strength decreased by 2 per cent. Our lowest ionisation rate models were not followed beyond the isothermal collapse phase due to the non-ideal MHD constraints on the timestep. Those with s had maximum densities at a given time that agreed with the hydrodynamical model within 10 per cent. Given the reasonable agreement with the hydrodynamical model and the orders of magnitude increase in runtime, we conclude that hydrodynamical models can be used to approximate non-ideal MHD models with s unless precision results are required.

We conclude that it is possible to reproduce ideal MHD and purely hydrodynamical collapses using non-ideal MHD given an appropriate cosmic ray ionisation rate. However, reaching either limit by cosmic ray ionisation alone is unlikely, since molecular clouds in the local neighbourhood have cosmic ray ionisation rates of s (e.g. Padovani et al., 2009; Neufeld & Wolfire, 2017). Even if the comic ray ionisation rate were to temporarily increase (e.g. due to a nearby supernova), the increase will not be sustained enough to shift the evolution into the ideal MHD regime for the lifetime of the initial collapse.

For cosmic ray ionisation rates of s, we can approximate a hydrodynamic collapse, however, it is improbable to reach these low ionisation rates in local diffuse molecular clouds given that this rate is well below the rate expected from radionuclide decay. Studies of the early universe by Susa et al. (2015) found that, even in the absence of cosmic rays and heavy metals, the primordial gas is partially ionised (i.e. contains e, H, Li). Their electron fractions for this case are similar to ours for s, however, their H, Li fractions are higher due to the absence of dust grains. Although our low ionisation rate models are unlikely to be relevant in the local Universe, they may have implications for the early Universe.

While our quantitative conclusions depend on our initial conditions and our chosen maximum densities, we find that the canonical cosmic ray ionisation rate of s approaches neither the ideal MHD limit nor the hydrodynamical limit. For star formation simulations to properly evolve the magnetic field, non-ideal MHD is essential.


JW and MRB acknowledge support from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007- 2013 grant agreement no. 339248). DJP received funding via Australian Research Council grants FT130100034, DP130102078 and DP180104235. This work was supported by resources on the swinSTAR national facility at Swinburne University of Technology. swinSTAR is funded by Swinburne and the Australian Government’s Education Investment Fund.


  • Alexiades et al. (1996) Alexiades V., Amiez G., Gremaud P.-A., 1996, Commun. Numer. Meth. Eng., 12, 31
  • Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  • Bate et al. (2014) Bate M. R., Tricco T. S., Price D. J., 2014, MNRAS, 437, 77
  • Bourke et al. (2001) Bourke T. L., Myers P. C., Robinson G., Hyland A. R., 2001, ApJ, 554, 916
  • Braiding & Wardle (2012a) Braiding C. R., Wardle M., 2012a, MNRAS, 422, 261
  • Braiding & Wardle (2012b) Braiding C. R., Wardle M., 2012b, MNRAS, 427, 3188
  • Choi et al. (2009) Choi E., Kim J., Wiita P. J., 2009, ApJS, 181, 413
  • Ciolek & Mouschovias (1994) Ciolek G. E., Mouschovias T. C., 1994, ApJ, 425, 142
  • Commerçon et al. (2010) Commerçon B., Hennebelle P., Audit E., Chabrier G., Teyssier R., 2010, A&A, 510, L3
  • Crutcher (1999) Crutcher R. M., 1999, ApJ, 520, 706
  • Dapp & Basu (2010) Dapp W. B., Basu S., 2010, A&A, 521, L56
  • Dapp et al. (2012) Dapp W. B., Basu S., Kunz M. W., 2012, A&A, 541, A35
  • Draine & Lee (1984) Draine B. T., Lee H. M., 1984, ApJ, 285, 89
  • Duffin & Pudritz (2009) Duffin D. F., Pudritz R. E., 2009, ApJ, 706, L46
  • Fiedler & Mouschovias (1993) Fiedler R. A., Mouschovias T. C., 1993, ApJ, 415, 680
  • Gammie (1996) Gammie C. F., 1996, ApJ, 457, 355
  • Heiles & Crutcher (2005) Heiles C., Crutcher R., 2005, in Wielebinski R., Beck R., eds, Lecture Notes in Physics, Berlin Springer Verlag Vol. 664, Cosmic Magnetic Fields. p. 137 (arXiv:astro-ph/0501550), doi:10.1007/11369875_7
  • Hennebelle & Ciardi (2009) Hennebelle P., Ciardi A., 2009, A&A, 506, L29
  • Hennebelle & Fromang (2008) Hennebelle P., Fromang S., 2008, A&A, 477, 9
  • Igea & Glassgold (1999) Igea J., Glassgold A. E., 1999, ApJ, 518, 848
  • Keith & Wardle (2014) Keith S. L., Wardle M., 2014, MNRAS, 440, 89
  • Krasnopolsky et al. (2012) Krasnopolsky R., Li Z.-Y., Shang H., Zhao B., 2012, ApJ, 757, 77
  • Kunz & Mouschovias (2009) Kunz M. W., Mouschovias T. C., 2009, ApJ, 693, 1895
  • Larson (1969) Larson R. B., 1969, MNRAS, 145, 271
  • Larson (1972) Larson R. B., 1972, MNRAS, 156, 437
  • Li & Shu (1996) Li Z.-Y., Shu F. H., 1996, ApJ, 472, 211
  • Li et al. (2011) Li Z.-Y., Krasnopolsky R., Shang H., 2011, ApJ, 738, 180
  • Mac Low et al. (1995) Mac Low M.-M., Norman M. L., Konigl A., Wardle M., 1995, ApJ, 442, 726
  • Machida et al. (2011) Machida M. N., Inutsuka S.-I., Matsumoto T., 2011, PASJ, 63, 555
  • Mathis et al. (1977) Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425
  • Mellon & Li (2009) Mellon R. R., Li Z.-Y., 2009, ApJ, 698, 922
  • Mestel & Spitzer (1956) Mestel L., Spitzer Jr. L., 1956, MNRAS, 116, 503
  • Meyer et al. (2012) Meyer C. D., Balsara D. S., Aslam T. D., 2012, MNRAS, 422, 2102
  • Mouschovias (1996) Mouschovias T. C., 1996, in Tsinganos K. C., Ferrari A., eds, NATO Advanced Science Institutes (ASI) Series C Vol. 481, NATO Advanced Science Institutes (ASI) Series C. pp 505–538
  • Mouschovias & Ciolek (1999) Mouschovias T. C., Ciolek G. E., 1999, in Lada C. J., Kylafis N. D., eds, NATO Advanced Science Institutes (ASI) Series C Vol. 540, NATO Advanced Science Institutes (ASI) Series C. p. 305
  • Mouschovias & Spitzer (1976) Mouschovias T. C., Spitzer Jr. L., 1976, ApJ, 210, 326
  • Nakano & Umebayashi (1986a) Nakano T., Umebayashi T., 1986a, MNRAS, 218, 663
  • Nakano & Umebayashi (1986b) Nakano T., Umebayashi T., 1986b, MNRAS, 221, 319
  • Nakano et al. (2002) Nakano T., Nishi R., Umebayashi T., 2002, ApJ, 573, 199
  • Neufeld & Wolfire (2017) Neufeld D. A., Wolfire M. G., 2017, ApJ, 845, 163
  • O’Sullivan & Downes (2007) O’Sullivan S., Downes T. P., 2007, MNRAS, 376, 1648
  • Padovani et al. (2009) Padovani M., Galli D., Glassgold A. E., 2009, A&A, 501, 619
  • Pandey & Wardle (2008) Pandey B. P., Wardle M., 2008, MNRAS, 385, 2269
  • Pollack et al. (1994) Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush T., Fong W., 1994, ApJ, 421, 615
  • Price (2012) Price D. J., 2012, Journal of Computational Physics, 231, 759
  • Price & Bate (2007) Price D. J., Bate M. R., 2007, MNRAS, 377, 77
  • Price & Monaghan (2004) Price D. J., Monaghan J. J., 2004, MNRAS, 348, 123
  • Price & Monaghan (2005) Price D. J., Monaghan J. J., 2005, MNRAS, 364, 384
  • Price et al. (2017) Price D. J., et al., 2017, preprint, (arXiv:1702.03930)
  • Seifried et al. (2011) Seifried D., Banerjee R., Klessen R. S., Duffin D., Pudritz R. E., 2011, MNRAS, 417, 1054
  • Shu et al. (2006) Shu F. H., Galli D., Lizano S., Cai M., 2006, ApJ, 647, 382
  • Spitzer & Tomasko (1968) Spitzer Jr. L., Tomasko M. G., 1968, ApJ, 152, 971
  • Susa et al. (2015) Susa H., Doi K., Omukai K., 2015, ApJ, 801, 13
  • Tassis & Mouschovias (2007) Tassis K., Mouschovias T. C., 2007, ApJ, 660, 370
  • Tomida et al. (2013) Tomida K., Tomisaka K., Matsumoto T., Hori Y., Okuzumi S., Machida M. N., Saigo K., 2013, ApJ, 763, 6
  • Tomida et al. (2015) Tomida K., Okuzumi S., Machida M. N., 2015, ApJ, 801, 117
  • Tricco & Price (2012) Tricco T. S., Price D. J., 2012, Journal of Computational Physics, 231, 7214
  • Tricco et al. (2016) Tricco T. S., Price D. J., Bate M. R., 2016, Journal of Computational Physics, 322, 326
  • Troland & Crutcher (2008) Troland T. H., Crutcher R. M., 2008, ApJ, 680, 457
  • Tscharnuter (1987) Tscharnuter W. M., 1987, A&A, 188, 55
  • Tsukamoto et al. (2015a) Tsukamoto Y., Iwasaki K., Okuzumi S., Machida M. N., Inutsuka S., 2015a, MNRAS, 452, 278
  • Tsukamoto et al. (2015b) Tsukamoto Y., Iwasaki K., Okuzumi S., Machida M. N., Inutsuka S., 2015b, ApJ, 810, L26
  • Tsukamoto et al. (2017) Tsukamoto Y., Okuzumi S., Iwasaki K., Machida M. N., Inutsuka S.-i., 2017, PASJ, 69, 95
  • Turner & Sano (2008) Turner N. J., Sano T., 2008, ApJ, 679, L131
  • Umebayashi & Nakano (1981) Umebayashi T., Nakano T., 1981, PASJ, 33, 617
  • Umebayashi & Nakano (1990) Umebayashi T., Nakano T., 1990, MNRAS, 243, 103
  • Umebayashi & Nakano (2009) Umebayashi T., Nakano T., 2009, ApJ, 690, 69
  • Wardle (2007) Wardle M., 2007, Ap&SS, 311, 35
  • Wardle & Ng (1999) Wardle M., Ng C., 1999, MNRAS, 303, 239
  • Wurster (2016) Wurster J., 2016, PASA, 33, e041
  • Wurster et al. (2014) Wurster J., Price D., Ayliffe B., 2014, MNRAS, 444, 1104
  • Wurster et al. (2016) Wurster J., Price D. J., Bate M. R., 2016, MNRAS, 457, 1037
  • Wurster et al. (2017a) Wurster J., Bate M. R., Price D. J., Tricco T. S., 2017a, preprint, (arXiv:1706.07721)
  • Wurster et al. (2017b) Wurster J., Price D. J., Bate M. R., 2017b, MNRAS, 466, 1788
  • Wurster et al. (2018) Wurster J., Bate M. R., Price D. J., 2018, MNRAS, 475, 1859
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description