The density–magnetic field relation in the atomic ISM

Abstract

We present numerical experiments aimed to study the correlation between the magnetic field strength, B, and the density, n, in the cold atomic interstellar medium (CNM). We analyse 24 magnetohydrodynamic models with different initial magnetic field intensities (B₀ = 0.4, 2.1, 4.2, and 8.3 |$\mu$|G) and/or mean densities (2, 3, and 4 cm⁻³), in the presence of driven and decaying turbulence, with and without self-gravity, in a cubic computational domain with 100 pc by side. Our main findings are as follows: (i) For forced simulations that reproduce the main observed physical conditions of the CNM in the solar neighbourhood, a positive correlation between B and n develops for all the B₀ values. (ii) The density at which this correlation becomes significant (≲30 cm⁻³) depends on B₀ but is not sensitive to the presence of self-gravity. (iii) The effect of self-gravity, when noticeable, consists of producing a shallower correlation at high densities, suggesting that, in the studied regime, self-gravity induces motions along the field lines. (iv) Self-gravitating decaying models where the CNM is subsonic and sub-Alfvénic with β ≲ 1 develop a high-density positive correlation whose slopes are consistent with a constant β(n). (v) Decaying models where the low-density CNM is subsonic and sub-Alfvénic with β > 1 show a negative correlation at intermediate densities, followed by a high-density positive correlation.

MHD, ISM: magnetic fields, ISM: structure, local interstellar matter

1 INTRODUCTION

The relation between the hydrogen particle density and the magnetic field intensity in the interstellar medium has been studied during more than three decades through Zeeman effect measurements from different tracers corresponding to different density ranges. For the atomic interstellar gas, and more specifically for the cold neutral medium (CNM) traced through the 21 cm Hi line, a weak or null correlation has been reported by Troland & Heiles (1986). For molecular gas, in contrast, a positive correlation has been found from OH and CN data (Crutcher 1999). The emergence of a regime, where the magnetic field is correlated with density, has been traditionally interpreted as the result of self-gravity effects becoming important in generating gas motions perpendicular to the field lines that in turn produce a compression of the freezing-in magnetic field. However alternative explanations, based exclusively on the interactions between gas motions and magnetic field fluctuations, have also been proposed. On one hand, scaling relations between the density and the magnetic pressure produced by different magnetohydrodynamic (MHD) wave modes in the non-linear regime have been analytically found (McKee & Zweibel 1995; Passot & Vázquez-Semadeni 2003). On the other hand, numerical simulations of compressible magnetized turbulent isothermal gas have shown that a positive correlation between the density and the magnetic field intensity develops when supersonic motions are present, both in the super-Alfvénic as well as in the sub-Alfvénic case (Padoan & Nordlund 1999; Ostriker, Stone & Gammie 2001; Burkhart et al. 2009; Collins et al. 2011).

The most recent observational version of this relation has been given by Crutcher et al. (2010), who compiled data from the available surveys tracing gas with densities going from few tens to 10⁶ cm⁻³. They used these data to find the most probable maximum value for the total magnetic field intensity in the case where no correlation is assumed for densities below a transition density n_Ht, and a positive correlation, with an exponent α, is assumed for higher densities. In that work, Bayesian statistics is used to find n_Ht = 300 cm⁻³, α = 0.65, and the maximum magnetic field value for the no correlation regime, B_max = 10 |$\mu$|G. The value of n_Ht corresponds to the mean hydrogen particle density in giant molecular clouds, favouring the previously mentioned interpretation in terms of the effects of self-gravity. Nevertheless this value lies also in the density range corresponding to the CNM gas, where supersonic motions are thought to be common, cooling is non-negligible, and self-gravity is usually subdominant.

Additionally, recent observations have found a preference for the projected magnetic field on the plane of the sky to be parallel to matter structures in the CNM, traced by Hi (Clark, Peek & Putman 2014) as well as by dust emission (Planck Collaboration XXXII 2016). This alignment could be a consequence of the presence of super-Alfvénic turbulence, a regime in which gas compressions produced by supersonic motions are able to generate an enhancement of the magnetic field component perpendicular to motions (Planck Collaboration XXXII2016), suggesting the possible presence of a correlation between B and n_H in the CNM. The same trend in the alignment between the plane of the sky magnetic field and the dense structures has been recovered in numerical simulations of different types: colliding flows in a magnetized medium (Inoue & Inutsuka 2016) and driven turbulence models (Villagran & Gazol 2017). In the latter work we have found that the alignment does also hold for the 3D local magnetic field, reinforcing the possibility of magnetic field compression. This brings up the question of how the magnetic field intensity is related with density in a gas with thermal and dynamical properties similar to those reported for the CNM and where regions are predominantly magnetically subcritical (Heiles & Troland 2004).

Zeeman observations from h i absorption are consistent with an uncorrelated magnetic field intensity with density, but they are not conclusive. Actually, the available data, analysed in Crutcher et al. (2010), make it hard to distinguish between the following three different cases: (i) all the observed clouds have the same total magnetic field, (ii) the total magnetic field intensity of each region is drawn from a uniform distribution, or (iii) the distribution of total magnetic field intensity leading to the observed line-of-sight field is the one obtained with the numerical model presented by Falceta–Gonçalves, Lazarian & Kowal (2008). This fact, along with the large uncertainties on the line-of-sight magnetic field measures and on the volumetric density determination – from the spin temperature of each component and a uniform pressure (3000 K cm⁻³) assumption –, and with the alignment results previously mentioned, constitute a strong motivation for a detailed investigation of the B–n_H relation in physical conditions similar to those of the CNM.

This relation has also been examined in the context of non-isothermal simulations. The simulations utilized to this end can be roughly divided in two groups; on one hand models considering a supernovae driven ISM (Balsara & Kim 2005; de Avillez & Breitschwerdt 2005; Hill et al. 2012; Hennebelle & Iffrig 2014; Pardi et al. 2017) and, on the other, models aimed to investigate molecular cloud formation (Hennebelle et al. 2008; Banerjee et al. 2009; Heitsch, Stone & Hartmann 2009; Inoue & Inutsuka 2012; Fogerty et al. 2016; Wu et al. 2017). In the first group, numerical requirements needed to include all the relevant physics and/or to take into account large scales (≥1 kpc) do not allow to use enough resolution to study the CNM gas. In the second group, the focus has been put on the study of the properties and the evolution of molecular gas that needs to resolve small scales in high-density regions. For this reason, these studies frequently use numerical schemes that give resolution priority to dense regions (i.e. smoothed particle hydrodynamics or adaptive mesh refinement) but treat more diffuse CNM-like gas with relatively low resolution. In the more recent works of this group, where the diffuse gas has enough resolution, no systematic attention has been paid to the B–n relation in this gas.

The aim of this work is thus to make a systematic numerical study of the relationship between the density and the magnetic field intensity in the CNM, by analysing models where the relevant physical ingredients of this kind of gas in the solar neighbourhood are taken into account, and where an adequate resolution is used.

The plan of the paper goes as follows. In the next section we briefly summarize some previous theoretical results concerning the relation between B and n, then, in Section 3 we describe the set-up of the models that we analyse in Section 4. These results are discussed in Section 5 and, finally, a conclusion is presented in Section 6.

2 THEORETICAL PREDICTIONS

In this section we briefly summarize some theoretical predictions for the B−n relation that are useful for the discussion and interpretation of our results.

The case of an isothermal spherical cloud dominated by gravity and contracting under its effects has been studied by Mestel (1965), who predicted B ∝ n^2/3. The same behaviour is expected whenever isotropic gravitational contraction is present in flux-freezing conditions, i.e. when the magnetic field is not strong enough to modify the geometry of the contracting region. When contraction happens under a constant ratio of thermal (P_th) to magnetic (P_B) pressure, β = P_th/P_B, and P_th ∝ n^γ, then B ∝ n^γ/2 (Mouschovias 1991). Thus, for the particular case of isothermal gas at constant β, B ∝ n^1/2. In this case the mass to magnetic flux ratio of the region increases as the contraction proceeds because this process is not isotropic. This behaviour is expected when the presence of the magnetic field has some influence on the contraction geometry.

When gravitational contraction is present, the magnetic field can also be amplified by the small-scale turbulent dynamo produced by the conversion of potential energy into turbulent kinetic energy (e.g. Federrath et al. 2011; Federrath 2016). As this process is highly dependent on the magnetic diffusivity, its contribution to numerical results is highly dependent on the resolution. Federrath et al. (2011) argued that a minimum number of 30 grid cells per Jeans length is required in order to observe some effects of this kind of dynamo.

The transition from an uncorrelated regime at low densities to a regime characterized by a correlation with a value of 2/3 for the exponent has also been interpreted as the result of diffusion caused by the turbulent magnetic reconnection (Lazarian, Esquivel & Crutcher 2012). In that work the authors present a simple reconnection diffusion model to show that this process can effectively remove magnetic field in molecular cloud envelopes and argue that, being more efficient at low densities, there is a scale range where reconnection diffusion can act on time-scales shorter than those of gravitational collapse.

On the other hand, the relationship between the magnetic pressure produced by undamped Alfvén waves and the density has been studied by McKee & Zweibel (1995). They found that, in hydrostatic equilibrium, this pressure, P_A, is proportional to n^1/2; while in a slow-varying gas (moderate Alfvénic Mach number, M_A, and long wavelength), P_A ∝ n^3/2 for three cases: (i) isotropic Alfvén waves, (ii) spatially uniform density, and (iii) self-similar contraction. Finally they argue that in the presence of a shock, the exponent relating P_A with n could be as large as two , depending on the shock direction. These results have been recovered by Passot & Vázquez–Semadeni (2003), who performed a perturbation analysis on top of circularly polarized Alfvén waves. They also investigated the scaling of the magnetic pressure produced by magneto-sonic modes with density. To this end they used simple waves (Boillat 1970; Webb et al. 1996) in a slab geometry, and found that the magnetic pressure associated with the fast mode is P_B ∝ n², while for the slow mode P_B ∝ c₁ − c₂n (with c₁ and c₂ being constants). The latter expression implies that when a slow mode is dominant, P_B is independent of density at low enough densities and it decreases n for higher values. Comparing the different scalings resulting from their analysis Passot & Vázquez–Semadeni (2003) concluded that at low densities the B–n relation is dominated by the slow mode and no correlation is expected, while at high enough densities, when the contribution of the slow mode disappears, a positive correlation develops due to the domination of the Alfvén and the fast modes. Note however, that the application of this analysis to a gas with physical properties similar to those of the CNM remains unclear. Mainly because magneto-sonic scalings are derived in a limit where important density fluctuations with weak magnetic field perturbation, a situation expected in the presence of supersonic motions and a strong magnetic field, cannot be considered. Note that a negative correlation has been numerically observed for 3D isothermal subsonic and sub-Alfvénic turbulence (Burkhart et al. 2009), although its presence seems to depend on the driving scheme, and more specifically on the correlation time associated with successive forcing events, whose effects become relevant at transonic velocities (Yoon, Cho & Kim 2016).

3 THE SIMULATIONS

In this paper, we analyse the results of 24 simulations that correspond to six groups with four different configurations: with/without self-gravity and forced/decaying turbulence. Each group has different initial magnetic field intensity, B₀, and/or different initial density, n₀ (see Table 1). Forced simulations are part of models analysed in Villagran & Gazol (2017, we keep the same labels for these models), which have been restarted with turbulent forcing turned off in order to produce the decaying turbulence models that we present in this work (labelled with an additional D). The main physical properties of the CNM in the solar neighbourhood are reproduced by forced models, particularly by those with n₀ = 2 cm⁻³, B₀= 4.2, and 8.3 |$\mu$|G. In fact, as discussed in detail in Villagran & Gazol (2017), the Alfvénic Mach number histogram for the CNM peaks at 1.7 and 1, for 4.2 and 8.3 |$\mu$|G, respectively. These values encompass the observationally determined estimation of 1.4 for a mean magnetic field intensity of 6 |$\mu$|G (Heiles & Troland 2003). The corresponding sonic Mach number histograms for the CNM resulting from all the forced models have peaks at moderately supersonic values (≈3.6) thatare also in agreement with the observational estimation by Heiles & Troland (2003).

Table 1.

Run parameters.

Model	n0	T0	B0	Forcing	Self-gravity
	cm−3	K	|$\mu$|G
B01S	2	1500	0.4	On	off
B01G	2	1500	0.4	On	on
B01SD	2	1500	0.4	Off	off
B01GD	2	1500	0.4	Off	on
B05S	2	1500	2.1	On	off
B05G	2	1500	2.1	On	on
B05SD	2	1500	2.1	Off	Off
B05GD	2	1500	2.1	Off	On
B10S	2	1500	4.2	On	Off
B10G	2	1500	4.2	On	On
B10SD	2	1500	4.2	Off	Off
B10GD	2	1500	4.2	Off	On
B20S	2	1500	8.3	On	Off
B20G	2	1500	8.3	On	On
B20SD	2	1500	8.3	Off	Off
B20GD	2	1500	8.3	Off	On
B10Sn3	3	733	4.2	On	Off
B10Gn3	3	733	4.2	On	On
B10SDn3	3	733	4.2	Off	Off
B10GDn3	3	733	4.2	Off	On
B10Sn4	4	518	4.2	On	Off
B10Gn4	4	518	4.2	On	On
B10SDn4	4	518	4.2	Off	Off
B10GDn4	4	518	4.2	Off	On

All the simulations result from the integration of the ideal MHD equations in the presence of cooling, using 512³ grid cells to represent a periodic cubic box of 100 pc by side, filled with gas at mean density and mean pressure lying in the thermally unstable regime (Field 1965) of the atomic interstellar gas. For this, we use a Total Variation Dismishing (TVD) code based on the one presented by Kim et al. (1999). In the forced models, the initial uniform density, n₀, and temperature, T₀, have thermal equilibrium values according to our cooling function. The cooling function is a fit, based on Wolfire et al. (2003) results for the solarG alactocentric radius (Gazol & Villagran 2016), that assumes a constant heating rate Γ = 22.4 × 10⁻²⁷ erg cm³s⁻¹. Driven models start at rest with a uniform magnetic field, |$\boldsymbol {B_0}=B_0\boldsymbol {x}$|⁠, and are forced in the Fourier space with solenoidal forcing applied at a fixed wavenumber corresponding to a physical scale of 50 pc. The amplitude of velocity perturbations is fixed by a constant injection rate of kinetic energy that we chose so as to approximately get the same rms velocity for all the driven models. The self-gravitating models are very gravitationally stable, having Jeans lengths of 552, 382, and 278 pc, for n₀ = 2, 3, and 4 cm⁻³, respectively. For the resolution we use, the Truelove criterion (Truelove et al. 1997) to avoid artificial fragmentation is satisfied for densities below ≈6350 cm⁻³, while the maximum density to resolve small-scale dynamo produced by gravitational-induced motions (Federrath et al. 2011) is ≈5000 cm⁻³. For further details concerning the model and the code see Villagran & Gazol (2017), Gazol & Kim (2010), and Kim et al. (1999).

4 RESULTS

In this section, we analyse the relation between the magnetic field intensity and the density for data resulting from the models described in the previous section. As a first approach, we take into account all the gas in the box and we analyse the B–n relation at densities comparable with those of the CNM, making no distinction between individual dense structures. In order to characterize the behaviour of B, we use the median magnetic field at each density value, B_m, computed from the 2D histogram. Note that we have tested the effects of using the mean value of B instead of its median value and we did not find important differences for the results presented below. For the sonic (M_s) and the Alfvénic Mach numbers, we also utilize their median value at each density resulting from the corresponding 2D histogram.

4.1 The effect of varying B₀

4.1.1 Forced turbulence

We first present results from forced models with n₀ = 2 cm⁻³ and different initial magnetic field intensities with and without self-gravity.

From Figs 1 and 2, where the median value of M_s and M_A as a function of density is displayed for non-self-gravitating and self-gravitating models, respectively, it can be seen that for all the models in this set, the gas with densities similar to those of the CNM (n ≳ 10 cm⁻³) is moderately supersonic and super-Alfvénic, in agreement with observations (Heiles & Troland 2003, 2004). As the β plasma parameter is related with M_s and M_A by |$\beta =2M_\mathrm{ A}^2/M_\mathrm{ s}^2$|⁠, the same figures imply that in our models, the dense gas has β ≲ 1 that is also consistent with observations (Heiles & Crutcher 2005).

Figure 1.

Sonic (triangles) and Alfvénic (squares) Mach numbers for forced models without self-gravity. The symbols represent the median value of M_s and M_A at each density. Black straight lines at log n = 1 and log M = 0 are included as references.

Figure 2.

Sonic (triangles) and Alfvénic (squares) Mach numbers for forced models with self-gravity. The symbols represent the median value of M_s and M_A at each density. Black straight lines at log n = 1 and log M = 0 are included as references.

In Fig. 3, we display B_m as a function of density for models with (empty circles) and without (filled circles) self-gravity for the four B₀ values that we consider. These plots have been obtained from averaging over the last 2.9 turbulent crossing times. First of all, it can be noticed that this set of models covers a region in the n−B plane in which a large number of the available Hi Zeeman data (black crosses; Heiles & Troland 2004) are included. The second noticeable feature is that the magnetic field intensity grows with density for logn ≥ 1.5 for all the B₀ values we consider. For low B₀ models, the rise in B_m starts at smaller densities, where the gas motions are already pre-dominantly super-Alfvénic (see Figs 1 and 2). As B₀ increases, the density at which B_m starts to grow with density, as well as the density at which gas motions become mainly super-Alfvénic, shifts to larger values. Finally, Fig. 3 shows that the effect of self-gravity is not appreciablefor high B₀ models but, as expected, it becomes stronger as B₀ decreases. However, in this case the increment of B_m with n is shallower when self-gravity is included, which is an unexpected behaviour, suggesting the presence of density enhancements without the corresponding increase of B. This could be a consequence of gravitational force producing motions parallel to the field lines, i.e. motions that do not compress the magnetic field lines. Unfortunately, our simulations do not have enough high-density points in order to quantify this effect. Although the growth of logB_m with log n does not seem to follow a linear relation, for reference, we computed the slope of a linear least-squares fit for 2 ≤ logn ≤ 2.6. For this density range, the mean polytropic index according to our thermal equilibrium curve is 0.7. For the non-self-gravitating runs we find a slope α of 0.43, 0.38, 0.28, and 0.27 for B₀ = 0.4, 2.1, 4.2, and 8.3 |$\mu$|G, respectively; while for the corresponding self-gravitating models α = 0.30, 0.28, 0.32, and 0.32. These values indicate that, for the low B₀ cases, the presence of self-gravity changes values of α from >γ/2, for which the magnetic pressure increases faster with density than thermal pressure (see Section 2), to values <γ/2, where the opposite is true, i.e. the presence of self-gravity causes a qualitative change in the B−n relation for low B₀ values.

Figure 3.

B–n relation for forced models with n₀ = 2 cm⁻³ as traced by the median of the 2Dhistogram. Filled(empty) circles are for models without(with) self-gravity. For each initial magnetic field, continuous(dashed) lines are for B values containing 25|$\hbox{ per cent}$| of points below and above the median for models without(with) self-gravity. Black lines are the most probable maximum values found by Crutcher et al. (2010), while black crosses are the H i data from Heiles & Troland (2004).

4.1.2 Decaying turbulence

Although the models discussed in the previous section satisfactorily reproduce a number of observed CNM properties (see also Villagran & Gazol 2017), it seems natural to raise the questions of how B is related to n in the case of decaying turbulence, and which are the effects of self-gravity in this case. The B–n relation, as measured again by the median value of the 2D histogram, for unforced restarts of models presented in Section 4.1.1 are shown in Figs 4 and 5 for the non-self-gravitating and the self-gravitating group, respectively; while the corresponding plots for M_s and M_A as a function of n are displayed in Figs 6 and 7. These figures have been done for times at which more than 85 per cent of the initial kinetic energy E₀ has been dissipated. In order to improve the statistics we average over 10 snapshots during which the fraction of dissipated energy varies among runs between 0.2 and 2 per cent of E₀. It can be seen that the models with a relatively high magnetic field, for which the CNM remains supersonic but becomes sub-Alfvénic for the low-density zone, show an almost flat B_m with small (≲100 cm⁻³) maximum densities, and their behaviour does not seem to be affected by the presence of self-gravity. In contrast, in simulations with low B₀, for which the CNM is predominantly subsonic and sub-Alfvénic at low densities, B_m varies with n in a more complex way that is affected by self-gravity. In fact, at densities below the CNM minimum density these two models show a slight (B₀ = 2.1 |$\mu$|G) or moderate (B₀ = 0.4 |$\mu$|G) increase of B_m with n, but in a density range starting at n ≈ 10 cm⁻³B_m clearly decreases with n. The width of this range is smaller for the lowest B₀ value and seems to be independent of self-gravity. In both cases, with and without self-gravity, B_m starts to grow with n again at large enough densities. Non-self-gravitating low B₀ models reach maximum density values even lower than those attained by models with a larger initial magnetic field. In contrast, self-gravitating low B₀ models show a sustained growth of the magnetic field intensity with density over a density range that, for B₀ = 0.4 |$\mu$|G (green symbols in Fig. 5) goes up to n ≈ 10³ cm⁻³. From Fig. 7 it can be seen that in this case, for densities at which B_m increases with n, the median value of the sonic Mach number slowly increases with n while M_A remains relatively constant at trans-Alfvénic values, implying a decline in the β value. For this model, the slope of a linear fit for the B–n relation at 10² ≤ n ≤ 10³ cm⁻³ is ≈0.41, which is larger than γ/2, consistent with a magnetic pressure growing faster with n than thermal pressure.

Figure 4.

B–n_H relation for decaying models with n₀ = 2 cm⁻³ without self-gravity. For each initial magnetic field continuous (dashed) lines are for B values containing 25|$\hbox{ per cent}$| of points below and above the median for models without(with) self-gravity. Black lines are the most probable maximum values found by Crutcher et al. (2010), while black crosses are the Hi data from Heiles & Troland (2004).

Figure 5.

Same as Fig. 4 for decaying self-gravitating models with n₀ = 2 cm⁻³.

Figure 6.

Sonic (triangles) and Alfvénic (squares) Mach numbers for decaying models with n₀ = 2 cm⁻³ and without self-gravity.

Figure 7.

Sonic (triangles) and Alfvénic (squares) Mach numbers for self-gravitating decaying models with n₀ = 2 cm⁻³.

4.2 The effect of varying the initial density

In this section, we present results from the models with higher mean density and B₀ = 4.2 |$\mu$|G both in the forced and in the decaying case. These models are useful to further study the effects of self-gravity.

4.2.1 Forced turbulence

Fig. 8 shows the B–n relation for forced high-density models compared to models with n₀ = 2 cm⁻³. For these figures we used the same time interval as for the figures presented in Section 4.1.1. As expected, the effects of self-gravity are noticeable at high densities (n ≳ 100 cm⁻³) where, unexpectedly, for the more massive model the magnetic field intensity grows slower with density, compared with the non-self-gravitating case. In fact the slopes of a linear least-squares fit for the same density range as the one used in Section 4.1.1 are α = 0.28, 0.26, and 0.30 for n₀ = 2, 3, and 4 cm⁻³, respectively, for the non-self-gravitating cases; while the corresponding values for the self-gravitating models are α = 0.32, 0.28, and 0.24. Thus, as in Section 4.1.1, when the relative relevance of self-gravity increases, B_m presents a slower growth with n. Note that for the models with increased mass, the behaviour of M_s and M_A does not considerably differ from the one we have shown in Figs 1 and 2 for B₀ = 4.2 |$\mu$|G.

Figure 8.

B–n_H relation for forced models with n₀ = 2 cm⁻³ as traced by the median of the 2Dhistogram. Filled(empty) circles are for models without(with) self-gravity. For each initial magnetic field continuous (dashed) lines are for B values containing 25|$\hbox{ per cent}$| of points below and above the median for models without (with) self-gravity. Black lines are the most probable maximum values found by Crutcher et al. (2010), while black crosses are the Hi data from Heiles & Troland (2004).

4.2.2 Decaying turbulence

In the decaying case, self-gravitating massive models develop a B–n relation qualitatively different from the one resulting from the n₀ = 2 cm⁻³ simulation. This can be observed in Fig. 9, where a very slight decay in B_m can be seen for all the models at densities below ≈30 cm⁻³. For non-self-gravitating runs and for the self-gravitating standard density model, this decrement continues at higher densities (note that this decay is similar to the one in Figs 4 and 5 for the corresponding B₀, but it looks stronger due to the different scales we use), but for massive self-gravitating models B_m starts to grow with density at ≈30 cm⁻³. This increase remains for more than an order of magnitude in density.

Figure 9.

B–n_H relation for decaying models with n₀ = 2 cm⁻³, n₀ = 3 cm⁻³, and n₀ = 4 cm⁻³ with (empty circles) and without (filled circles) self-gravity.

From Figs 10 and 11 where the median values of the sonic and the Alfvénic Mach numbers for each density bin resulting from non self-gravitating and self-gravitating models are plotted, it can be seen that in all the models in this set M_A ≲ M_s, implying that the magnetic pressure is dominant over or balances the thermal pressure. Also, as in the models discussed in Section 4.1.2, the simulations for which self-gravity leads to a regime where B_m increases with n are those with predominately subsonic motions at CNM densities. For the most massive case, it also can be seen that (i) at n ≳ 100 cm⁻³M_s ≈ M_A, i.e. magnetic and thermal pressures are approximately equal and (ii) at n ≥ 200 cm⁻³ the gas becomes supersonic and super-Alfvénic. The slope of the B_m−n relation in this case (cyan points in Fig. 9) is ≈0.37 for 10² ≤ n ≤ 10³ cm⁻³ that is consistent with α ≈ γ/2. For the n₀ = 3 cm⁻³ model, the high-density gas is near the constant β regime, but does not reach it and remains subsonic and sub-Alfvénic, suggesting that the effects of self-gravity are not strong enough to produce an important acceleration of the gas.

Figure 10.

Sonic (triangles) and Alfvénic (squares) Mach numbers for decaying models with different n₀ without self-gravity.

Figure 11.

Sonic (triangles) and Alfvénic (squares) Mach numbers for decaying self-gravitating models with different n₀.

5 DISCUSSION

The results described in the previous section show that the magnetic field intensity could have a variety of scalings with the density in a gas with thermal properties similar to those of the local CNM, and formed by the development of thermal instability. Despite this diversity in the B–n relation, all the models that we presented have B and n values consistent with at least part of the observed Hi data reported by Heiles & Troland (2004) and used by Crutcher et al. (2010). The low-density region of the most probable maximum total B value found by the latter authors appears in fact as an upper limit to our models. However, in order to further compare them with observations, a detailed analysis of the physical properties in the resulting individual dense structures is needed.

5.1 Forced turbulence

When the gas is forced to velocities leading to sonic Mach numbers consistent with those derived from observations, magnetic field compression produces median M_A values for the dense gas in the trans-Alfvénic or in the slightly super-Alfvénic regime (M_A < 5), which are also consistent with observationally obtained values. At the minimum CNM density, we get M_A ranging between ≈0.8 and ≈3 (see Figs 1 and 2). Thus, at densities corresponding to the CNM all our forced models are moderately supersonic and slightly super-Alfvénic. Even for this relatively narrow range in the dense gas M_A, we obtain different behaviours for the B–n relation in gas with densities similar to those of the CNM.

The combination of a flat low-density regime with an approximate power law at higher densities, hypothesized by Crutcher et al. (2010) as the most probable maximum B, is recovered for the higher B₀ models (whose Alfvénic Mach numbers are closer to the observed values), while a continuous increase of the median B with density is observed for the lower B₀ models. Thus, at large enough densities, an increase of B with n is observed for all the initial magnetic field values that we considered. This trend with B₀ has also been obtained by Mocz et al. (2017), who followed the collapse of gas with physical properties similar to those of pre-stellar cores, (i.e. self-gravitating, super-critical, isothermal, and highly supersonic gas) with initial magnetic fields corresponding to Alfvénic Mach numbers of 35, 3.5, 1.2, and 0.35. These values are in fact comparable to those obtained for our forced models using the values of the initial Alfvén speed and the rms velocity (see table 2 in Villagran & Gazol 2017) that are 11.25, 1.12, 1.05, and 0.50, respectively. We observe the same trend with and without self-gravity. Then, the emergence of a combined regime and particularly the appearance of a positive correlation, seems to be independent of the presence of self-gravity. In fact, in our models the density from which B starts to show an important rise is independent of the presence of self-gravity and does not seem to be modified when more mass is available, at least for B₀ = 4.2 |$\mu$|G. Even for models with B₀ = 8.3 |$\mu$|G, in which turbulent compression drives the median B to values close to the most probable maximum found by Crutcher et al. (2010) of 10 |$\mu$|G, the rise of B with density starts at n ≈ 30 cm⁻³ that is approximately one order of magnitude below the value found by these authors for the switch between the constant B regime and the power-law regime. This is consistent with the B–n relation presented in previous numerical works studying molecular cloud formation, where the diffuse gas is sufficiently resolved. Inoue & Inutsuka (2016) consider the formation of molecular clouds due to converging magnetized Hi flows in the absence of self-gravity, finding that the magnetic field intensity increases with n from relatively low-density values. Also Fogerty et al. (2016), who investigate the role of shear and magnetic fields in molecular cloud formation, report B–n distributions consistent with a power-law growth starting ≈10 cm⁻³. In Wu et al. (2017), where the collision of two GMC is addressed, a positive correlation can be observed at densities of about 10 cm⁻³. Note however that Heitsch et al. (2009), who study the effects of magnetic field on the flow-driven formation of molecular clouds in non-self-gravitating models, find a very weak correlation between B and n for different magnetic field strengths and orientations. The discrepancy with our results could be due to the low |$v$|_rms(M_s) they get for the cold gas, which are ≈0.4–0.15 times the ones we obtain in the forced models, that are consistent with the CNM observations.

Without self-gravity, the slope of the logarithmic B–n relation for 100 ≲ n ≲ 300 cm⁻³ increases with M_A. This is consistent with high M_A gas being more prone to magnetic field compressions and with numerical results for a turbulent isothermal gas (Padoan & Nordlund 1999; Ostriker et al. 2001; Burkhart et al. 2009; Collins et al. 2011). In the presence of self-gravity, however, the slope is less affected by B₀ differences, with low B₀ models showing a shallower growth of B with n. In fact, for the high M_A models, the presence of self-gravity switches α from values >γ/2 to values <γ/2, that implies that the gas switches from a regime where the magnetic pressure grows with density faster than the thermal pressure to a regime with the opposite behaviour. In the presence of self-gravity, the slope is around 0.3, for the four forced models with n₀ = 2 cm⁻³. A shallower power law in the presence of gravity is also observed for the model with the highest mean density and B₀ = 4.2 |$\mu$|G but in this case, where for the non-self-gravitating model α < γ/2, the change does not cause a qualitative switch. It is worth noting that the models for which self-gravity produces a shallower B–n relation are those for which Villagran & Gazol (2017) find the highest level of alignment between |$\mathbf {B}$| and ∇n. The origin of the reduction in the correlation could be the possible presence of density enhancements not associated with dynamic compressions, presumably associated with motions parallel to the local magnetic field, which we expect to be more relevant when magnetic fields are relatively less important with respect to gravity, i.e. for small B₀ or large n₀.

Note that the opposite effect of self-gravity, i.e. a steepening of the B–n relation, has been reported by Collins et al. (2011), who found that when no self-gravity is included, the magnetic field scales with density with an exponent 0.4, while after the gravity is turned on, the value becomes 0.48. However, our physical setup (see Section 3) is different from theirs. In fact, they follow the evolution of an isothermal highly supercritical, highly supersonic (M_s ≈ 10), and super-Alfvénic region because the scope of their work was to study the properties of a turbulent collapsing region. Their setup is similar to the one recently used by Mocz et al. (2017), who find a decrease of the scaling index for the strongest initial magnetic field (120 |$\mu$|G), for which α switches from ∼2/3 to ∼1/2.

5.2 Decaying turbulence

In the decaying case, the CNM becomes sub-Alfvénic at the CNM low-density limit for all our models, and we get three different behaviours. (i) No significant correlation between B_m an n, besides a very slight decrease of B_m, and low (<100 cm⁻³) maximal densities. This occurs for two types of models. On one side, the high B₀ and n₀ = 2 cm⁻³ models both, with and without self-gravity, where the CNM is supersonic by the time we use for the analysis and for all the density values. On the other side, models with n₀ = 3 and 4 cm⁻³ without self-gravity, where the dense gas is predominantly subsonic at the CNM lower transition density but becomes supersonic at larger but still relatively low densities. For the two types of models with this behaviour the CNM becomes super-Alfvénic at large enough densities, but M_A remains smaller than or comparable with M_s, indicating that magnetic pressure dominates over or balances thermal pressure. (ii) A very slight decrease of B_m followed by a positive correlation is observed for high-density models including self-gravity, for which the CNM remains mainly subsonic and sub-Alfvénic (with M_A ≲ M_s) for n ≲ 200 cm⁻³, and a trend towards constant β as n increases is observed. (iii) A clear negative correlation for the low-density CNM followed by a positive correlation is observed for low B₀ models. By the time we analyse these models, the CNM has become subsonic at the low CNM density limit but, in contrast with the models with the two previous behaviours, the CNM has β ≥ 1. Interestingly, the maximum density with a negative correlation seems to depend more on the B₀ value than on the presence of self-gravity, suggesting an MHD origin for the switch on the B–n relation. In these models, the presence of self-gravity produces a high-density B growth that is better defined and expands to larger densities.

The first behaviour is in contrast with the results reported by previous works for isothermal gas, where supersonic motions do always lead to a positive correlation between B and n. However, this difference is somehow expected, because in a cooling gas M_s> 1 results from the combined effects of motions, which can potentially compress the magnetic field lines, and the lowering temperature, that is not accompanied by a contraction so that it does not affect the magnetic field intensity.

The second behaviour does also differ from what has been found in isothermal gases but in this case the difference can be attributed to the presence of self-gravity in a region with a sufficiently large amount of mass. It is interesting, however that the positive correlation does also develop from n ≈ 30 cm⁻³.

The interpretation of the B_m decay observed in (ii) and (iii) is not straightforward. On one hand, accordingly with Passot & Vázquez–Semadeni (2003), this negative correlation suggests the domination of the slow magneto-sonic mode, and the transition to a positive correlation could be due to the transition towards a regime dominated by the fast or the Alfvén mode. However, in order to confirm this origin a detailed analysis of the relative contribution of different MHD modes should be performed. On the other hand, the decay of B_m could also be due to fast turbulent reconnection (Lazarian & Vishniac 1999; Lazarian et al. 2012). In this case, in order to verify this hypothesis, a detailed study of the parameters characterizing the turbulence that is present in the models where we observe the decay (i.e. injection scale and spectrum scaling) isneeded.

Note that the fact that the second and the third behaviour appearing at M_s< 1 makes them less plausible to occur in the real CNM, at least in the solar neighbourhood. Whether or not they can develop in cold gas with subsonic velocity dispersion depends also on the particular thermal and magnetic properties of that gas. For the third behaviour, there is the additional requirement of β > 1 that is also inconsistent with observations characterizing the local CNM as a β < 1 gas. It is interesting, however, that the negative correlation that we report appears in an unforced case, where the effects of the driving mechanism (Yoon et al. 2016) do not affect.

6 CONCLUSIONS

From the results presented in this work we can conclude the following:

• At realistic rms velocities and magnetic field strengths, a cooling gas with thermal pressure and temperature similar to those of the local CNM develops a positive correlation between the magnetic field intensity and the density starting at n ≲ 30 cm⁻³. The emergence of this positive correlation does not depend on the presence of self-gravity. The presence of this positive correlation in a gas with supersonic and super-Alfvénic motions is consistent with previous results from isothermal simulations.

• In this velocity regime, a B–n relation constituted by the combination of a flat B at low densities followed by a power law at high densities develops for a strong enough initial magnetic field intensity. This is consistent with the assumption made by Crutcher et al. (2010) for the most probable maximum magnetic field intensity according to observational data.

• The slope of this positive correlation depends on the relative relevance of self-gravity, but the density at which it emerges seems to be independent of it. For our – highly gravitationally stable – models, the effect of self-gravity is to produce a shallower B–n relation in models with either, a relatively weak initial magnetic field or a relatively high mean density. This effect suggests that, for this parameter choice, self-gravity induces motions parallel to the field lines.

• For decaying models with the CNM reaching subsonic and sub-Alfvénic regimes with β ≲ 1, a positive correlation develops only in the presence of self-gravity. The slope of this correlation at large densities is consistent with a constant β(n) regime.

• For decaying models where the low-density CNM is subsonic and sub-Alfvénic with β > 1, we observe a narrow density range with a negative correlation followed by a positive one, presumably both originated by an MHD process.

ACKNOWLEDGEMENTS

We acknowledge the support of UNAM-DGAPA-PAPIIT (Universidad Nacional Autóonoma de Méexico-Direccióon General de Asuntos del Personal Acadéemico-Programa de Apoyo a Proyectos de Investigacióon e Inovacióon Tecnolóogica) through the grant no.: IN110316. This work has made extensive use of the NASA’s Astrophysics Data System Abstract Service.