EDP Sciences
Free Access
Issue
A&A
Volume 600, April 2017
Article Number A114
Number of page(s) 12
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201629905
Published online 12 April 2017

© ESO, 2017

1. Introduction

It is widely accepted that molecular clouds are multiphase entities where the dense and cold gas, called the cold neutral medium (CNM), is mixed to a more diffuse and warm phase, known as the warm neutral medium (WNM). Interstellar turbulence not only leads to large density fluctuations (Falgarone et al. 1995), but it also mixes up these two phases. Using 1D simulations Lesaffre et al. (2007) have shown that the turbulent diffusion can transport molecules into warmer regions by spreading out the transition layers between the CNM and the WNM. More recently, in Valdivia et al. (2016, hereafter Paper I), we showed a similar effect in 3D magnetohydrodynamical (MHD) simulations of turbulent molecular clouds. In particular we showed that molecular hydrogen (H2) can be formed at intermediate densities in transient structures and subsequently carried toward warmer regions by the turbulent motions within the gas. This warm H2 is in good agreement with observations of excited populations of H2 in translucent molecular clouds (Rachford et al. 2002; Gry et al. 2002; Lacour et al. 2005), and it is reasonably successful at explaining these populations as collisionally excited states.

A possible related problem is the formation of the methylidine cation CH+ in the interstellar medium (ISM), which requires molecular hydrogen to be formed efficiently. CH+ was one of the first molecules observed in space (Douglas & Herzberg 1941), and it has been detected in a wide variety of environments, but its ubiquity and high abundance in the diffuse environments in the ISM has been a puzzling problem for more than 70 yr. In the Milky Way, 12CH+ has been detected in visible wavelengths in the local ISM (Crane et al. 1995; Gredel 1997; Weselak et al. 2008), while the isotopologue has been detected in infrared wavelengths deeper in the Galactic disc (Falgarone et al. 2005, 2010; Godard et al. 2012). This molecule has even been observed in the ISM of external galaxies (Rangwala et al. 2011; van der Werf et al. 2010; Ritchey et al. 2015; Spinoglio et al. 2012).

CH+ is easily destroyed by reactions with electrons and hydrogen atoms, but also by reactions with H2 molecules. Under these conditions, the only reaction efficient enough to counterbalance the fast destruction is (1)As this reaction is highly endothermic (Agúndez et al. 2010), it has been proposed that the observed abundances might be related to warm layers of gas resulting, for instance, from turbulent mixing or turbulent dissipation processes (Crane et al. 1995; Gredel 1997; Lambert & Danks 1986).

From the theoretical point of view, the attemps to explain the large amounts of CH+ in the diffuse ISM include shock heated gas (Draine & Katz 1986; Pineau des Forêts et al. 1986; Flower & Pineau des Forêts 1998; Lesaffre et al. 2013), turbulent dissipation in vortices (Joulain et al. 1998; Godard et al. 2009, 2014; Falgarone et al. 2010; Myers et al. 2015), and the turbulent mixing between the CNM and the WNM (Lesaffre et al. 2007). Several of these models have proven successful in reproducing not only the abundances of CH+, but also those of many other species (e.g. SH+, CO, HCO+) and their correlations observed in the diffuse medium under the constraint of the mechanical energy injected at large scale. However, they did so assuming idealised 1D steady-state structures and adopting either a single type of structure along the line of sight or very simplistic velocity, density, or magnetic field distributions with no link to the large scale dynamics of the gas.

In their recent work, Myers et al. (2015) were able to reproduce the observed CH+ abundances for a turbulent molecular cloud with typical physical conditions for the ISM. Using ideal MHD simulations, they showed that almost all of the CH+ molecules are produced in regions where the ion-neutral drift velocities are high (vd ≳ 3−4 km s-1). CH+ production is controlled by drift velocities higher than 3 km s-1 which, although rare, are present in a sufficient number to be statistically significant and to produce CH+ column densities comparable to observations. However, there are several caveats related to the simplifying assumptions that might be influencing their results. More specifically, the ionisation fraction is assumed to be constant and equal to the ionisation fraction of the dense phase, which might overestimate the drift velocities in the diffuse phase (where most of the CH+ is produced). Another incorrect physical assumption is that the H2 molecular fraction is assumed to be 10% throughout the entire domain. Taken all together, these assumptions likely lead to an artificially high CH+ abundance.

The goal of our paper is to study the formation of CH+ on a large scale in diffuse molecular clouds. We use a dynamically calculated abundance of molecular hydrogen (out of equilibrium H2) and include a more detailed description of the microphysical processes. We compare our data with available observations and with the previous theoretical work of Myers et al. (2015, hereafter MML15). First, we explore the role of the turbulent mixing, which transports molecular hydrogen to the warm gas, in a realistic 3D structure by assessing what role is played by the warm reactants. Second, we address the question of whether the approximated treatment of the ambipolar diffusion in the ideal MHD limit can explain the observed abundances of CH+ under realistic physical conditions.

The paper is structured as follows. In Sect. 2 we describe the numerical simulation and the chemical solver used to post-process it. We also describe how we estimate the ion-neutral slip velocity, and we analyse the validity of our approach. In Sect. 3 we present our results on the importance of warm reactants, and the importance of the ion-neutral drift, on the abundance of CH+. We conclude the paper in Sect. 4. We give a detailed description of our chemical solver and a detailed study of the chemical timescales for each individual species in the Appendix.

2. Methodology

2.1. MHD simulation

thumbnail Fig. 1

Slice showing the total number density and the H2 fraction. The arrows depict the velocity field projected onto the xy plane.

Open with DEXTER

The numerical simulation used in this work is an ideal1 MHD multiphase simulation of a realistic molecular cloud using the RAMSES code (Teyssier 2002); it is fully described in the appendix of Paper I.

The molecular cloud is formed by colliding streams of atomic gas and the general set-up is very similar to Valdivia & Hennebelle (2014) (see also Audit & Hennebelle 2005; Hennebelle et al. 2008). The simulation box is a cube of side L = 50 pc filled with atomic gas of density n = 1 cm-3 and temperature T = 8000 K, embedded in the interstellar radiation field (ISRF) of strength G0 = 1.7 in Habing units (Habing 1968), where two converging streams of WNM are injected from the x boundaries with a slightly turbulent stationary velocity field of module V0 = 15 km s-1. The corresponding injection of mechanical energy at the integral scale is Ėin ~ 3 × 1033 erg s-1, or equivalently , in agreement with the mean value of the kinetic energy transfer rate deduced from the observed velocity dispersions of molecular clouds seen in CO (Hennebelle & Falgarone 2012). The magnetic field is initially aligned with the x direction and it has an initial strength of 2.5 μG.

The simulation uses an adaptive mesh refinement (AMR) technique. The minimum and maximum resolution levels are min = 9 and max = 11, reaching an effective numerical resolution of 20483 or, equivalently, a spatial resolution of some 0.025 pc. The refinement criterion is density, with density thresholds at nthresh = 50 cm-3, and nthresh = 100 cm-3. The time step of the simulation is variable and smaller than ~ 180 yr, i.e. small enough to follow the propagation of dynamical perturbations in the CNM and WNM (~ 104 yr).

The simulation follows the formation and destruction of H2 and its thermal feedback. The effects of shielding from ultraviolet radiation by dust particles and H2 absorption (self-shielding) have both been included for the computation of the H2 photodissociation rate using our tree-based method, detailed in Valdivia & Hennebelle (2014) and in Paper I. The simulation produces turbulent and highly structured molecular clouds that exhibit a wide range of physical conditions (in density, temperature, and shielding parameters). Important fractions of out-of-equilibrium H2 are found in the diffuse and warm phase of the gas as a consequence of the turbulent mixing and the shielding against photodissociation provided by the multiphase structure. Molecules of H2 initially formed in transient dense regions are spread by turbulence into the interclump medium where they can survive thanks to the shielding provided by the cloud structure. The presence of H2 in a warm environment increases its excitation (Paper I) and leads to column densities of H2 rotational levels comparable to those observed with the Copernicus and FUSE Telescopes (e.g. Spitzer et al. 1974; Frisch 1980; Frisch & Jura 1980; Lambert & Danks 1986; Gry et al. 2002; Lacour et al. 2005).

The snapshot used in this work corresponds to an evolution time of 15 Myr in the simulation. Figure 1 shows the local number density in a cut through the middle plane, as well as a small region showing the local H2 fraction. This figure depicts the influence of dynamics on the transport of H2 molecules.

2.2. Method description

To understand the role of warm and out of equilibrium H2 and of the multiphase structure of molecular clouds on the chemistry of the ISM, we perform a post-treatment of the numerical simulation described above. The data are extracted with the python module PyMSES (Labadens et al. 2012). Once extracted, the chemical composition of every cell of the simulation is computed assuming chemical equilibrium for all species except for H and H2.

The chemical solver used for the post-treatment is taken from the Meudon photon dominated regions (PDR) code2 (e.g. Le Petit et al. 2006; Bron et al. 2014), stripped of surface reaction processes and detailed treatments of grain physics and radiative transfer, and modified to allow the abundance of H2 to be fixed beforehand. For each cell, the solver takes as input the local properties of the gas, i.e. the total hydrogen density nH, the kinetic temperature TK, the external UV radiation field χ, the visual extinction AV, the ion-neutral velocity drift vd, the abundance of H2n(H2), and the shielding of H2 from UV photons: (Paper I). The MHD simulation is ideal and hence the ion-neutral drift also has to be estimated by a post-processing treatment (see Sect. 2.4 below). In output, the solver returns the at-equilibrium abundances of the 147 variable species included in the chemical network. Detailed descriptions and tests of the solver and of the method used to fix the value of n(H2) are presented in Appendix A.

It is important to note that in this work we only consider the shielding of CO by dust, and neglect both the self-shielding and the shielding induced by H2 line absorption.

2.3. Timescales and transient chemistry

Computing the composition of the gas at chemical equilibrium with a fixed abundance of H2 is valid only if the timescales required to reach the equilibrium are smaller than that of H2, but also smaller than the typical timescales of variation of dynamical quantities (nH, v, TK). To check these assumptions, the chemical timescales have been estimated in Appendix A as functions of the physical properties of the gas, in particular the density and temperature. These timescales (see Fig. A.2) are found to vary between 102 and 106 yr depending on the species considered and the gas density. CH+, for instance, reaches its equilibrium abundance in ~ 2 × 106/nH yr, regardless of the gas temperature.

A comparison with the results obtained for H2 shows that molecular hydrogen has an evolution time longer than that of any other species over almost the whole range of physical conditions spanned in the simulation. This result not only proves the existence of an equilibrium solution when the abundance of H2 is fixed, it also indicates that this solution is a coherent description of the chemical composition of the gas over time yr for cm-3 and yr for cm-3.

The comparison with dynamical timescales is, unfortunately, less satisfactory. With the spatial resolution of the simulation, perturbations of physical conditions in the WNM and CNM propagates at sound speed over timescales of ~ 104 yr. Therefore, while computing the chemistry at equilibrium may be marginally valid for high-density gas ( cm-3), it is probably not appropriate for cells at lower density. In this case, an accurate treatment would require following the time-dependent evolution of the entire chemistry, which is still out of reach of any numerical simulation.

While the problem seems hopeless, we note that the chemical timescales presented in Appendix A.3.3 have been calculated assuming that only H2 is out of equilibrium in the molecular clouds. It is probable, however, that only a few other species evolve, like H2, over a long period, and need to be treated on the fly in the simulation. If so, computing the equilibrium abundances of all species except for a few critical ones would be a correct and very efficient method. Therefore, although our approach is not without fault, we argue that it gives a first estimation of the chemical state towards which the gas tends. It also provides a framework for future developments where other species could be considered out of equilibrium depending on the dynamical timescales, the properties of the gas, and the chemical history. All these developments are underway.

Another limitation of computing the chemistry at equilibrium is that it neglects any transient chemical event which can only be accounted for with time averaged abundances. However, we argue that taking into account such a transient chemical event within each fluid cell is not necessarily correct as long as the dynamics is not resolved over the same temporal scales.

2.4. Ion-neutral drift

The ion-neutral drift, or ambipolar diffusion, can help to increase the rates of highly endothermic reactions which scale as exp(−ΔE/kT), where ΔE is the endothermicity. The relative motion between neutrals and ions increases the effective velocity dispersion, which results in the increase of the effective temperature at which the reaction occurs (Draine 1980; Draine et al. 1983; Flower et al. 1985): (2)where mi (mn) and Ti (Tn) are the mass and temperature of ions (neutrals), and ΔT is (3)In this expression vd is the ion-neutral drift velocity, μ is the reduced mass of the reaction, and k is the Boltzmann constant. When ions and neutrals have the same temperature (equal to the gas temperature, T), the effective temperature can be approximated as TeffT + ΔT. But when T< ΔTβ, where β = ΔE/k is the endothermicity expressed in kelvins, this approximation overestimates the reaction rate (Pineau des Forêts et al. 1986), and the approximation is no longer valid. A more accurate result is obtained by taking the ion-neutral reaction rate to be proportional to exp(−max{β/Teff,(β−3ΔT) /T}).

Although our simulation is in the ideal MHD approximation, we estimate the ion-neutral drift by neglecting the ion inertia in the ion equation of motion. The resulting balance between Lorentz forces and the ion-neutral drag force yields (Shang et al. 2002; Glassgold et al. 2005) (4)where B is the magnetic field, ρn and ρi are respectively the densities of neutrals and ions, γ is the coefficient for ambipolar diffusion, nj and nk are the number densities of the ionic and neutral species, vd,jk = | vjvk | is the specific drift velocity for species j and k, μjk is the reduced mass, and Kjk is the momentum transfer rate coefficient. Assuming the same drift velocity for all the species, vd can be approximated as follows: (5)

thumbnail Fig. 2

Volume-weighted (left panel), m(H2)-weighted (middle panel), and m(CH+)-weighted (right panel) 2D probability density functions (PDF) of the gas temperature (Tgas) versus the H2 fraction (f(H2)).

Open with DEXTER

Table 1

Momentum transfer coefficients.

It is important to note that Kjk, and thus the drift velocity vd, depends on the rms velocity vrms (in km s-1). Assuming that all ions and neutrals have the same temperature, this velocity can be written as (6)which is dependent on the ion-neutral drift velocity itself, and on the thermal velocity (Pinto & Galli 2008).

Since the effective temperature varies as the squared value of the drift velocity (see Eq. (3)), errors can be greatly amplified. For this reason we use an iterative method to give an accurate value of vd.

For the momentum transfer rate coefficients we use the accurate expressions given by Pinto & Galli (2008), and the Langevin rates when a more accurate rate was not available. For the particular case of collisions between C+ and H2, we assumed a fraction of the value used for the collisions between C+ and H, which corresponds to the ratio between the corresponding Langevin rates. In Table 1 we summarise the interactions included in our calculations3.

3. Results

In order to compare our data with available observations, we selected a uniform grid of 1024 lines of sight for which we solved for the chemistry and integrated for column densities. Each line of sight has been analysed using a constant numerical resolution of 1024 cells in the line of sight, or equivalently ~ 0.05 pc, which corresponds to an intermediate resolution level of = 10 for the native AMR resolution of the simulation.

3.1. Role of warm H2

thumbnail Fig. 3

Line of sight showing the influence of warm H2. Each panel shows the total number density and the gas temperature (top), the number density of H2 calculated at equilibrium and dynamically (centre), and the number density of CH+ (bottom).

Open with DEXTER

In Paper I we showed the presence of H2, which is not at chemical equilibrium, in warm regions. As H2 is a prerequisite to the formation of other molecules, and as warm reactants can boost the formation of molecules with high reaction barriers, we analyse the role of warm H2 on the formation of CH+. For this purpose we calculate the formation of CH+ in two different ways. In the first case we calculate the chemical abundances of all the species in our chemical network with the H2 abundances obtained in our simulation (out of equilibrium H2), while in the second case we calculate the full chemical equilibrium for all the species including H2.

Figure 2 shows the 2D probability density functions (PDF) in the form of a 2D histogram of the gas volume, the H2 mass (m(H2)), and the CH+ mass (m(CH+)) distributions in the simulation box. The left and centre panels show that a large number of cells are characterised by low densities of H2, while a reduced number of cells concentrate most of the mass in H2 form. The right panel shows that in spite of this distribution, most of the CH+ is produced in regions with intermediate H2 fractions (f(H2) ~ 0.3−30%), which do not correspond to the regions that dominate the volume or the mass of H2, and with gas temperatures as high as several 102−103 K.

thumbnail Fig. 4

Comparison of column densities of CH+ as a function of the total column density for the case where H2 is fixed from the simulation (solid circles), and for the case where H2 is calculated at equilibrium (open circles). The crosses are the observational data from Crane et al. (1995), Gredel (1997), and Weselak et al. (2008).

Open with DEXTER

We present a single line of sight in Fig. 3 to shed some light on the physical conditions that give rise to an enhancement on the CH+ abundance. This figure shows the local physical condition of the gas, as well as a comparison of the H2 density calculated dynamically in our simulation and that expected at equilibrium for the same physical conditions (total density, temperature, dust shielding, and H2 self-shielding). The bottom panel of this figure shows that most of the CH+ is produced in regions that present specific characteristics favourable to CH+ formation: the fraction of H2 is higher than is predicted at equilibrium, and temperatures are of the order of several 100 K. In both cases most of the CH+ is produced near the edge of the clumps, but the abundance of CH+ is up to three orders of magnitude larger when the warm H2 is used.

The resulting column densities for the grid of 1024 lines of sight are shown in Fig. 4. The column densities of CH+ are 3−10 times higher than those obtained with H2 at equilibrium and closer to the observed ones. However, the abundances are still underpredicted by a factor of ~ 6. This means that the warm H2 plays an important role in the production of CH+, although it is not enough to explain the observed abundances.

3.2. Role of the ion-neutral drift

thumbnail Fig. 5

Normalised volume-weighted drift velocity probability distribution (P), calculated as described in Sect. 2.4 (solid blue line) compared to what would be obtained using the γAD and x(e) prescriptions from MML15 (dashed red line). The analytic fit of MML15 is given for comparison (solid light blue line).

Open with DEXTER

The relative velocity between ions and neutrals adds a non-thermal component to the gas temperature that increases the effective temperature at which ion-neutral reactions occur. To assess the role of the ambipolar diffusion on the production of CH+, we first analyse the distribution of ion-neutral drift velocities and then its impact on the effective temperature distributions.

3.2.1. Drift velocity distributions

Figure 5 shows the distribution obtained by MML15 represented by their analytic fit (light blue line); the distribution of the drift velocity that we obtain using their prescription, which assumes that the dominant ion is C+, with a constant abundance of x(C+) = 1.6 × 10-4 and an ambipolar diffusion coefficient given by γ = 8.47 × 1013 cm3 s-1 g-1 (Draine 1980), applied to our data (dashed red line); and the distribution obtained for our data using our prescription as described in Sect. 2.4 (solid blue line).

When we apply the prescription of MML15 to our data (dashed red line), the distribution we obtain is different to their distribution. We obtain an asymmetric distribution, narrower around the maximum, with a flatter region at low drift velocities and displaying a sharply decaying tail at higher drift velocities. The maximum, when using our simulation data, is located at ~ 105 cm s-1, which is one order of magnitude higher than the MML15 distribution. Nevertheless, our distribution predicts fewer events (roughly one order of magnitude) of velocities vd ≳ 10 km s-1. The difference in shape might be a consequence of the biphasic nature of our simulation, while the higher peak can be easily explained by the different mean densities of the simulated clouds. As in a first approximation the drift velocity is inversely proportional to the squared density (see Eq. (4)); a mean density of in the case of MML15 and in our case leads to a difference of roughly a factor of 10.

When we use our prescription (using the iterative method and the ion densities calculated with our chemical solver) we obtain a distribution shifted towards lower drift velocities. The maximum of our distribution is located around ~ 4 × 103 cm s-1, and most of the points in the simulation show very weak drift velocities. Points with drift velocities higher than 1 km s-1 are extremely rare. These differences can be explained by the different assumptions we used. In the prescription of MML15 the electronic fraction is assumed to be constant and dominated by the C+ abundance, which is a valid assumption only in the dense and cold gas. This approximation underestimates the electronic fraction by one to two orders of magnitude in the more diffuse and warm gas where the electronic fraction is indeed dominated by the ionised hydrogen H+ (xe ~ 10-2−10-3, as shown in Fig. 8). It is important to remember that vd ∝ 1 /ni, and a difference of two orders of magnitude is amplified to four when calculating the effective temperature. In addition, the ambipolar coefficient γ is assumed to be determined solely by interations between C+ and H2.

thumbnail Fig. 6

Normalised mass-weighted probability distribution (P) of the gas temperature, log Tgas (dashed light blue line), and the effective temperature, log Teff (solid magenta line), showing the two phases of the ISM.

Open with DEXTER

In summary, this figure shows that the simplified prescription of MML15 applied to our data greatly overestimates the ion-neutral drift velocity distribution.

Figure 6 shows the distribution of the gas temperature and the distribution of the effective temperature. It reveals the presence of the two distinct phases of the ISM, the WNM and the CNM (Field et al. 1969), but also a transition phase that correspond to gas transiting through the thermally unstable phase between the two main phases. As the effect of the ambipolar diffusion is local and intermittent, it affects only a small fraction of the gas, and therefore the two phases are not expected to differ significantly in the effective temperature distribution. To noticeably change the relative distribution of the CNM and WNM, a mean drift velocity of at least 1 km s-1 over the bulk of the CNM is necessary. Differences should only appear in the intermediate temperature domain where the effective temperature is expected to vary from the gas temperature for sufficiently high drift velocities (vd ≳ 1 km s-1). However, Fig. 6 shows that the two distributions are almost indistinguishable, revealing that the effect of the ambipolar diffusion is negligible, at least under the conditions analysed in this work.

3.2.2. Role of ambipolar diffusion on the production of CH+

thumbnail Fig. 7

Distribution of the abundance of CH+ as a function of density and gas temperature (panel a)), and as a function of density and drift velocity (panel c)). Panels b) and d): contribution from the ambipolar diffusion.

Open with DEXTER

As we recalled in the previous section, both the gas temperature and the ion-neutral drift can increase the reaction rates. To disentangle the role of the ion-neutral drift on the production of the CH+ molecule, we calculated the equilibrium abundances of all the species (except for H and H2) in two different ways. For the first case we calculated the chemical abundances using only the gas temperature, while in the second case we computed it taking into account the fact that the reaction rates depend on the ion-neutral drift velocity vd, as defined in Sect. 2.4. With these two distributions we compute the excess produced by the ion-neutral drift as Δn(CH+) = n(CH+)(Teff)−n(CH+)(Tgas).

thumbnail Fig. 8

Example of line of sight, showing a case of strong ion-neutral drift.

Open with DEXTER

thumbnail Fig. 9

Column densities of CH+ as a function of the total column density. Including the ion-neutral drift (solid circles), and without including this effect (open circles). The crosses are the observational data as in Fig. 4.

Open with DEXTER

Figure 7 shows the distribution of CH+ and the contribution from the ambipolar diffusion Δn(CH+) as 2D histograms. The two leftmost panels show the distributions as a function of the gas number density and gas temperature, while the two rightmost panels show the same distributions as a function of the gas number density and the ion-neutral drift velocity. This figure shows that the ion-neutral drift produces a small but not always negligible contribution to the total production of CH+. A close comparison between the second and fourth panels shows that the production of CH+ increases at temperatures ranging from several hundreds to thousands K and at low densities up to some ~ 10 cm-3 due to the ion-neutral drift, consistent with high values of the drift velocity at the same density regime.

As shown in Figs. 5 and in 6 the effect of the ion-neutral drift is in general minor, but some regions of high drift velocities can arise. This is the case in Fig. 8, where we selected a line of sight that displays a case of strong ion-neutral drift. This line of sight shows that high drift velocities arise as narrow features in regions with low density, where the electron abundance is in general higher than 10-3, and dominated by H+ rather than by C+. The production of CH+ is locally strongly boosted (by ~ 1 order of magnitude), and although punctual and rare, such local events might have a great influence on the integrated column density. To check whether the ambipolar diffusion can have a significant impact on the global picture or if the effect is statistically unimportant, we integrate total column densities and CH+ column densities along the selected lines of sight for the two cases analysed in this section. The global effect of the drift velocity is shown in Fig. 9. This figure shows that, statistically, the effect of the ambipolar diffusion is almost imperceptible. Nevertheless, this effect is resolution-dependent, so probably the relevant scales for the ambipolar diffusion are underdescribed. This effect is explored in the following section.

3.3. Effect of the resolution on the vd distribution

Our estimate of vd from the ideal MHD simulations depends on, among other variables, the gradient of the magnetic field and the ion-neutral momentum transfer rate (sensitive to both the ionisation degree and the temperature). Both of these quantities are likely to be sensitive to the resolution, so it is important to analyse the role of the numerical resolution on the distribution of drift velocities.

To shed some light on the role of resolution, we present in Fig. 10 the distribution of drift velocities for three identical simulations, at 15 Myr, performed using three different numerical resolutions. Each simulation is characterised by the minimum (min) and maximum (max) resolution levels allowed in the AMR4, and by the constant intermediate resolution (dxa) used for the analysis. The low-resolution simulation has min = max = 7, and dxa ≃ 0.4 pc; the intermediate-resolution simulation has min = 8, max = 10, and dxa ≃ 0.1 pc; and the high-resolution simulation has min = 9, max = 11, and dxa ≃ 0.05 pc.

Figure 10 shows the distribution of the drift velocity obtained for each of the simulations. Since in our MHD simulations the only dissipative process is numerical truncation, the magnetic current (∇ × B) power will pile up at smaller scales for higher resolution, which can increase our estimate of vd. On the other hand, the numerical dissipation heating will also occur at smaller scales as the resolution increases, which will impact the temperature and the ionisation degree and may alter the momentum transfer rates in a way which is not easy to predict. The current resolution study indicates that the net effect is to shift the drift velocity distributions towards larger values for higher resolution simulations.

If the ambipolar diffusion is included in the MHD simulations, this implies accurately resolving the local ambipolar diffusion scale λd = πvA/ (γρi), where vA is the Alfvén speed (Ntormousi et al. 2016; Hennebelle & André 2013). This can be very stringent, although Momferratos et al. (2014) has shown that the largest drift velocities can be obtained for scales much larger than λd. In addition, observational constraints, based on kinematical observations of ionic and neutral species, set the drift scales at mpc-scales (Hezareh et al. 2014; Li & Houde 2008). Future studies should therefore strive to treat the dissipative processes more accurately, but it is not clear how this will affect our results.

thumbnail Fig. 10

Normalised probability distribution (P) of the ion-neutral drift velocity estimated from ideal MHD simulations at different numerical resolutions. Numbers in the legend indicate the lowest and highest resolution level () used for the AMR.

Open with DEXTER

4. Conclusions

We assessed the importance of two key factors that contribute to producing the methylidine cation CH+ in diffuse molecular clouds, namely the role of warm reactants and the role of ambipolar diffusion. To this end, we post-processed an ideal MHD simulation of a realistic turbulent molecular cloud, which includes the evolution of the thermal state of the gas and the formation of molecular hydrogen. For the post-processing we developed a chemical solver able to calculate the chemical equilibrium for a given set of physical conditions and able to fix the abundance of H2 as an input parameter if required. The chemical solver provides the abundances of all the species in the chemical network (149 species in our case) very efficiently.

We compared the expected abundances of CH+ obtained using a full equilibrium with those obtained when fixing the abundances of H2 to the value obtained dynamically in the simulation. We show that molecular hydrogen plays a fundamental role on the chemistry of the ISM. Since the evolution time of molecular hydrogen is long, its formation and evolution must be followed carefully in numerical simulations. In particular, we show that the excess of H2 found in warm gas as a consequence of the multiphase structure and the turbulent mixing within interstellar gas is a very important ingredient in the understanding of the warm chemistry.

More specifically, the formation of CH+ seems to be more efficient in regions where H2 is not expected at equilibrium. CH+ is formed mainly in regions characterised by H2 fractions close to ~ 1−10% and by temperatures higher than 300 K, most likely of the order of 103 K. These specific regions, such as the external layers of clumps, meet the necessary conditions to allow endothermic reactions to produce CH+ efficiently. Nevertheless, the abundances of CH+ are still underpredicted compared to observations (Crane et al. 1995; Gredel 1997; Weselak et al. 2008) by a factor of the order of 6, suggesting that missing elements might contribute in the same measure. A possible clue to solving this puzzle, other than the dissipation of turbulence, is the excitation of rotational and vibrational levels of H2 pumped by the surrounding UV field. Molecules of H2 transported towards regions where the shielding is mild are usually found in excited states, and the internal energy available in the molecule can be used to overcome the reaction barrier (Zanchet et al. 2013; Herráez-Aguilar et al. 2014). It is also worth noting that given the timescales involved, the excited states are likely to be out of equilibrium and therefore need to be followed in real time in the simulation. A different but related issue would be the initial and boundary conditions used in our simulation. Since real molecular clouds are likely formed by multiphase flows made up of a mixture of WNM and CNM gas already containing molecular gas, the abundance of H2 – mainly in low-density regions – might be underestimated. A larger fraction of H2 in such regions might lead to an increase in the CH+ abundance. These effects are not treated in this paper, but we plan to address them in future works.

We calculated the abundances of the dominant ionic species using our full chemical network. This, along with accurate momentum transfer rates, leads to a distribution of drift velocities dominated by low values (peak around ~ 4 × 103 cm s-1), where drift velocities higher than 1 km s-1 are extremely rare events. Our calculation of ion-neutral drift velocities does not produce a high-velocity tail, as has been shown in previous works (Myers et al. 2015). This highlights the importance of a good description of small-scale physical processes, such as the calculation of the electron density, as well as the description of the main ion-neutral interactions in order to avoid unrealistic estimates of the ion-neutral drift velocities. It is worth noting that since the corrective term ΔT depends quadratically on the drift velocity, and at the same time the formation rates are proportional to exponentials involving ΔT, overestimated drift velocities can lead to an artificial overabundance of CH+.

Our distribution of vd does not impact the effective temperature distribution, and thus is not able to explain the observed abundances of CH+. Even though it can increase locally the abundance of CH+, the role of the ambipolar diffusion is very limited because high drift velocities are rare events. However, we show that resolution effects could have a strong influence on the high-velocity tail. Therefore, the ambipolar diffusion physics and its scales need to be taken into account in future works in order to describe this process and its impact on the production of interstellar CH+.


1

It should be remembered that non-ideal MHD simulations are out of reach of current computational capabilities, requiring complicated implicit or semi-implicit schemes, or extremely small time steps.

2

Version 1.5.2 available at http://ism.obspm.fr

3

We do not include collisions with He, which further reduce the drift velocity distribution only slightly.

4

The associated resolution for a level is dx = L/ 2.

6

Fractional abundances of PAHs and grains are computed assuming a dust-to-gas mass ratio of 0.01, a PAH mass fraction of 4.6% (Draine & Li 2007), an MRN size distribution of grains ranging from 0.005 to 0.3 μm, and a log-normal distribution of PAHs centred on 6.4 Å.

Acknowledgments

We are grateful to Alexandre Faure, who kindly provided us with the Langevin rates, and to Edith Falgarone and Eva Ntormousi for valuable discussions. We acknowledge the financial support of the Agence Nationale pour la Recherche through the COSMIS project. This research has received funding from the European Research Council under the European Community Seventh Framework Programme (FP7/2007–2013 Grant Agreement No. 306483). We thank the French Programme Physique Chimie du Milieu Interstellaire (PCMI). This work was granted access to HPC resources of CINES under the allocation x2014047023 made by GENCI (Grand Equipement National de Calcul Intensif). This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

References

Appendix A: The chemical solver

The chemical solver built in this work is an adaptation of the solver used in the Meudon PDR (photodissociation regions) code (Le Petit et al. 2006), modified to simplify the treatment of a few chemical processes (e.g. surface reactions, H2 self shielding), and optimised to improve the computational speed. The solver, which we describe in greater detail below, is available on the ISM numerical platform of the Paris Observatory5.

Appendix A.1: Description

The solver considers a static fluid cell of density nH and temperature TK, irradiated by an external radiation field χ (expressed in Mathis’ unit, Mathis et al. 1983) attenuated by a visible extinction AV, and pervaded by cosmic ray particles. The cell contains NX species (resulting from the combinations of NA different atoms) that interact with each other through a given network of chemical reactions. In this configuration the code computes the chemical state of the fluid cell at equilibrium. Assuming that chemical reactions obey the law of mass action kinetics, this equilibrium is defined by the following set of algebraic equations (A.1)where n(X) and nA are the density of species X and the number of atoms A (in cm-3), N the number of reactions, Rj,k and κj the reactants and reaction rate of reaction j, sj(X) the stoichiometric coefficient of X in reaction j, m(A,X) the multiplicity of atom A in species X, and c(X) the charge of X. The chemical networks used in the solver are conservative; therefore, the first line above provides a system of NXNA−1 independent equations which is completed by conservation equations for each atom and for the electrons, chosen to replace the evolution equation of the most abundant atom carrier and charge carrier.

As described in Le Petit et al. (2006), the system of equations is solved with a Newton-Raphson scheme modified to prevent exploration of solutions with negative abundances. Iterations stop when the relative variation of abundances falls below a given threshold and if this solution corresponds to an equilibrium for the chemistry, i.e. when (A.2)where F(X) and D(X) are the total formation and destruction rates of X and ε1 and ε2 are two parameters set to 10-6 and 10-3, respectively. When called for the first time, the solver starts with initial conditions representative of the high-ionisation phase (Le Bourlot et al. 1995): all carbon and sulfur in C+ and S+, and all oxygen and nitrogen in O and N. When called in sequence, the initial conditions are the set of abundances corresponding to the last solution found by the solver. This method is done to favour continuity of chemical states computed on adjacent fluid cells. However, it does not preclude the existence of one or several other chemically stable solutions for a given set of physical conditions (e.g. low-ionisation phase, Le Bourlot et al. 1995).

Table A.1

Parameters of the chemical solver.

The main input parameters of the solver and the ranges of values they span in the simulation are given in Table A.1. In addition, the code accepts four optional input quantities: the ion-neutral velocity drift (see Sect. 2.4), the self-shielding factors for the photodissociation of H2 and CO (see Eq. (11) of Paper I), and the fractional abundance of H2 (see Sect. A.3).

Appendix A.2: Chemical network

The elemental abundances adopted in this work are set to the values observed in the solar neighbourhood and compiled by Flower & Pineau des Forêts (2003), assuming no mantles on grain surfaces. The chemical network is the most recent version of the network used in the Meudon PDR code, slightly adapted to simplify the treatment of a few chemical processes. This network is available online (http://ism.obspm.fr, network_valdivia_2016.dat) and contains 149 species interacting with each other through 2692 reactions.

Formation of H2 on grain surfaces and destruction by photodissociation are computed according to Paper I. Contrary to the PDR model, which performs a coherent calculation of grain surface chemistry including both the Eley-Rideal and the Langmuir-Hinshelwood mechanisms (Le Bourlot et al. 2012; Bron et al. 2014), the formation rate of H2 is computed via a simple function, scaled to the mean value of the formation rate observed in the ISM (Gry et al. 2002), and tuned to depict the dependence of the sticking coefficient of H on grains on the kinetic temperature. The photodissociation rate of H2 is set to 3.3 × 10-11χfsh, H2 s-1, where fsh, H2 is a parameter (see Table A.1) used to include the effects of shielding of the radiation field by H2 line and dust continuous absorptions (e.g. Draine & Bertoldi 1996). All other gas phase photoreaction rates are computed using the exponential fits of van Dishoeck (1988) and van Dishoeck et al. (2006) at a given value of the visible extinction.

The treatment of the charges of polycyclic aromatic hydrocarbons (PAHs) and very small grains and of the processes of electron transfer between those particles and the ions is similar to the treatment performed in the Turbulent Dissipation Region (TDR) model (Godard et al. 2014). PAHs and grains are described as spherical particles with radius of 6.4 Å and 0.01 μm, respectively, and fractional abundances6 of 4.2 × 10-7 and 4.8×10-10. Each particle is supposed to exist in charge states ranging from − 1 to 5. Photoionisations and photodetachments of all the corresponding species are modelled as photoreactions with rates deduced from detailed computations of the photoelectric effect with the Meudon PDR code. Electronic attachments and ion recombinations on dust (in their various charged states) are treated using the prescription of Draine & Sutin (1987) and Weingartner & Draine (2001a,b).

Appendix A.3: Fixed H2 abundance

The 3D MHD simulation of the diffuse medium used in this work was built to compute the time dependent evolution of the abundance of molecular hydrogen. In order to estimate the impact of out-of-equilibrium H2 on the chemistry, the solver offers the possibility of fixing the abundance of H2 in the cell. In this case, the code searches for a chemical state of the gas where almost all the species are at equilibrium.

This option is neither straightforward nor inconsequential. Fixing the abundance of H2 is equivalent to removing an equation from the first line of system A.1. System A.1 therefore becomes an overdetermined problem composed of NX−1 variables linked by NX independent equations. Consequently, if the abundance of H2 is fixed and out of equilibrium, either the system is not conservative or there is at least one other species out of equilibrium.

To simplify the problem, the solver works on the assumption that only one species, SH, is out of equilibrium in addition to H2. While this hypothesis is a good approximation in some cases, it raises several issues: What criteria should be used to select SH? Does this choice guarantee the existence of a solution over a wide range of physical conditions (temperature, density, and molecular fraction)? If a solution exists, is it a correct representation of the ISM chemistry induced by out of equilibrium H2?

Appendix A.3.1: Selection criterion

There are several arguments that constrain the selection of the out of equilibrium species SH. First, a chemical network used with a constant abundance of molecular hydrogen is conservative for all elements except H. A necessary condition to invert system A.1 is therefore to remove an evolution equation of an H-bearing species. Secondly, the chemical equilibrium state of the ISM is driven by its most abundant constituents. In order to find a stable solution, SH should thus be abundant enough so that its density can be slightly adjusted to conserve the number of H atoms without modifying the at-equilibrium abundances of all the other species. Finally, since the timescale allowed to reach an equilibrium is limited by the resolution of the simulation, SH should be among the species which react the most quickly to any variation in the abundance of H2, i.e. those most likely to be out of equilibrium if H2 is.

We therefore define SH as the most abundant H-bearing species besides H2. With the parameters explored in the Appendix and the physical conditions spanned in the simulation (see Table A.1), we find that such a criterion always leads to atomic hydrogen (H), both in the most diffuse and ionised phases of the ISM (nH ~ 0.1 cm-3, AV< 10-6) and in the fully molecular gas. For intermediate or large molecular fractions (f(H2) ≳ 10-2), the selection method fulfils all the constraints. In particular, H is not only out of equilibrium, but 1D time-dependent chemical models also show that its abundance varies oppositely to that of H2: dn(H)/dt = −dn(H2)/dt. In contrast, the method is found to be frailer at low molecular fraction (f(H2) ≲ 10-2). Indeed, in this case, the abundance of H is less dependent on that of H2, hence several species are more likely than H to be out of equilibrium. Despite this limitation, we find that selecting H is paramount even in these conditions in order to ensure the numerical stability of the solver.

It is worth noting that H is the main species out of equilibrium in a large majority of cases, which strongly supports the method used in Paper I: reducing the chemical network to H and H2 to compute the abundance of molecular hydrogen.

Appendix A.3.2: Existence of a solution

thumbnail Fig. A.1

Percentage of convergence failures of the chemical solver as functions of the molecular fraction and the density (left), the temperature (middle), or the visible extinction (right).

Open with DEXTER

Dynamical systems such as chemical networks under mass action kinetics are known to display a variety of asymptotic behaviours depending on their control parameters (e.g. Vidal & Lemarchand 1988; Feinberg 1987). In particular, fixing the abundance of one species, which is equivalent to adding a reservoir and a draining route to the network, can prevent the existence of a stable equilibrium state. For instance, if the abundance of H2 is too high, the steady state abundances of all other H-bearing species produced directly or indirectly via H2 cannot be small enough to guarantee the conservation of H atoms. In some cases, fixing the abundance of H2 can thus maintain the entire system out of equilibrium.

To identify such cases, we explored the convergence of the solver over a grid of 2.56 million models. The results displayed in Fig. A.1 reveal two regions in the 4D parameter space.

The most pre-eminent region corresponds to high-density (nH ≳ 103 cm-3), shielded (AV ≳ 1) environments with moderate molecular fraction (10-1f(H2) ≲ 1) and kinetic temperature (300 K ≲ TK ≲ 103 K). Testing those models without fixing the abundance of H2 shows similar results, which indicates that convergence failures are due to the Newton-Raphson algorithm and the topology of the Jacobian of system A.1 rather than the absence of fixed point for the chemistry. When such cases occur, we therefore choose to guide the algorithm by computing the solution at high temperature and progressively decrease the temperature towards the target value.

The least pre-eminent region in Fig. A.1 corresponds to low-density (nH ≲ 1 cm-3), illuminated (AV ≲ 0.5) environments with a high molecular fraction (f(H2) ≳ 0.95) and kinetic temperature (TK ≳ 103 K). As predicted above, convergence failures happen here because the abundance of H2 is fixed and higher than a critical value. However, since such physical conditions are never reached in the simulation, we find that fixing the abundance of H2 and SH never precludes the existence of an equilibrium state for the rest of the chemistry.

Appendix A.3.3: Chemical timescales

thumbnail Fig. A.2

Timescales required by the species to reach their equilibrium abundances as a function of the density of the gas for a temperature fixed at 50 K (left panel) and as a function of the temperature of the gas for a density of 10 cm-3 (right). The H2 timescales (coloured curves) are computed in configuration (a) (see main text) for fsh,H2 = 104 (solid), 102 (dotted), and 1 (dot-dashed). Each curve results from an average of timescales computed for ten random initial conditions and two values of AV (0.1 and 1). The timescales of single molecules (coloured points) and the range of timescales of all other species (black segments) are computed in configuration (b) (see main text). These timescales are averaged over ten random initial conditions, two values of AV (0.1 and 1), and four values of f(H2) (10-3, 10-2, 10-1, and 1).

Open with DEXTER

Finally, fixing H2 presumes that the timescales required for each species to reach their equilibrium state in a medium at constant molecular fraction is smaller than that required for H2. To test this hypothesis we have computed the equilibrium timescales τ(X) of every species X using the time-dependent integrator of the TDR code (Godard et al. 2014) set up to follow the chemical evolution of a fluid cell with constant density and temperature. We define the timescale τ(X) as the longest time for which the abundance of species X differs from its asymptotic value by more than 10%. With this definition, the integrator was run over grids of models of different density, temperature, and extinction in two different configurations: (a) considering all species as variables and adopting the H2 self-shielding factor fsh, H2 as a parameter and (b) assuming a constant molecular fraction f(H2). All models were finally run with ten different initial conditions randomly drawn on a logarithmic scale and rescaled to match the input elemental abundances (Glover et al. 2010 Appendix A).

The chemical timescales averaged over initial conditions, visual extinction, and H2 molecular fraction (for configuration b only) are displayed in Fig. A.2 as functions of the density and the temperature of the gas. Dispersion analysis shows that these equilibrium timescales are computed with an uncertainty of about a factor of 3. As expected, the equilibrium timescale of H2 derived in configuration (a) depends strongly on the self-shielding parameter, but only weakly on the density and the temperature of the gas. In contrast, the chemical timescales derived in configuration (b) are found to span a wide range of values (three orders of magnitude), to strongly depend on the gas density, and weakly depend on the gas temperature.

Figure A.2 shows that H2 is among the species with the longest evolution times. The equilibrium timescale of H2 is found to be larger than that of any other species over a wide range of physical parameters, as long as cm-3 or . These wide conditions are fulfilled in more than 90% of the mass of the gas (Paper I, Fig. 7), but also in more than 92% of the environments where CH+ is produced (see Fig. 7). It therefore not only validates the existence of a chemical equilibrium in a medium with a fixed molecular fraction, but also justifies our approach to single out H2 in the chemistry.

It is, finally, important to note that these considerations justify treating the chemistry at equilibrium only with respect to H2, but not necessarily with respect to the gas dynamics. A more detailed comparison between the chemical timescales and the dynamical resolution of the simulation and the importance of transient events are discussed in Sect. 2.3.

All Tables

Table 1

Momentum transfer coefficients.

Table A.1

Parameters of the chemical solver.

All Figures

thumbnail Fig. 1

Slice showing the total number density and the H2 fraction. The arrows depict the velocity field projected onto the xy plane.

Open with DEXTER
In the text
thumbnail Fig. 2

Volume-weighted (left panel), m(H2)-weighted (middle panel), and m(CH+)-weighted (right panel) 2D probability density functions (PDF) of the gas temperature (Tgas) versus the H2 fraction (f(H2)).

Open with DEXTER
In the text
thumbnail Fig. 3

Line of sight showing the influence of warm H2. Each panel shows the total number density and the gas temperature (top), the number density of H2 calculated at equilibrium and dynamically (centre), and the number density of CH+ (bottom).

Open with DEXTER
In the text
thumbnail Fig. 4

Comparison of column densities of CH+ as a function of the total column density for the case where H2 is fixed from the simulation (solid circles), and for the case where H2 is calculated at equilibrium (open circles). The crosses are the observational data from Crane et al. (1995), Gredel (1997), and Weselak et al. (2008).

Open with DEXTER
In the text
thumbnail Fig. 5

Normalised volume-weighted drift velocity probability distribution (P), calculated as described in Sect. 2.4 (solid blue line) compared to what would be obtained using the γAD and x(e) prescriptions from MML15 (dashed red line). The analytic fit of MML15 is given for comparison (solid light blue line).

Open with DEXTER
In the text
thumbnail Fig. 6

Normalised mass-weighted probability distribution (P) of the gas temperature, log Tgas (dashed light blue line), and the effective temperature, log Teff (solid magenta line), showing the two phases of the ISM.

Open with DEXTER
In the text
thumbnail Fig. 7

Distribution of the abundance of CH+ as a function of density and gas temperature (panel a)), and as a function of density and drift velocity (panel c)). Panels b) and d): contribution from the ambipolar diffusion.

Open with DEXTER
In the text
thumbnail Fig. 8

Example of line of sight, showing a case of strong ion-neutral drift.

Open with DEXTER
In the text
thumbnail Fig. 9

Column densities of CH+ as a function of the total column density. Including the ion-neutral drift (solid circles), and without including this effect (open circles). The crosses are the observational data as in Fig. 4.

Open with DEXTER
In the text
thumbnail Fig. 10

Normalised probability distribution (P) of the ion-neutral drift velocity estimated from ideal MHD simulations at different numerical resolutions. Numbers in the legend indicate the lowest and highest resolution level () used for the AMR.

Open with DEXTER
In the text
thumbnail Fig. A.1

Percentage of convergence failures of the chemical solver as functions of the molecular fraction and the density (left), the temperature (middle), or the visible extinction (right).

Open with DEXTER
In the text
thumbnail Fig. A.2

Timescales required by the species to reach their equilibrium abundances as a function of the density of the gas for a temperature fixed at 50 K (left panel) and as a function of the temperature of the gas for a density of 10 cm-3 (right). The H2 timescales (coloured curves) are computed in configuration (a) (see main text) for fsh,H2 = 104 (solid), 102 (dotted), and 1 (dot-dashed). Each curve results from an average of timescales computed for ten random initial conditions and two values of AV (0.1 and 1). The timescales of single molecules (coloured points) and the range of timescales of all other species (black segments) are computed in configuration (b) (see main text). These timescales are averaged over ten random initial conditions, two values of AV (0.1 and 1), and four values of f(H2) (10-3, 10-2, 10-1, and 1).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.