Simulating the atomic and molecular content of molecular clouds using probability distributions of physical parameters

TL;DR

This paper introduces a fast, scalable method for modeling the atomic and molecular composition of interstellar clouds using pre-calculated PDR calculations and probability distributions of physical parameters, bridging the gap between detailed simulations and large-scale observations.

Contribution

It presents a novel approach that models the ISM at large scales with an ensemble of 1D PDR calculations linked to hydrodynamical simulation data, enabling efficient analysis of diverse environments.

Findings

Low density medium remains HI- and CII-dominated under all conditions.

Denser clouds are mostly molecular, but carbon phase varies with environmental parameters.

The framework is fast, flexible, and suitable for large-scale ISM studies.

Abstract

Modern observations of the interstellar medium (ISM) in galaxies detect a variety of atomic and molecular species. The goal is to connect these observations to the astrochemical properties of the ISM. 3D hydro-chemical simulations attempt this but due to extreme computational cost, they have to rely on simplified chemical networks and are bound to individual case studies. We present an alternative approach which models the ISM at larger scales by an ensemble of pre-calculated 1D thermo-chemical photodissociation region (PDR) calculations that determine the abundance and excitation of atomic and molecular species. We adopt lognormal distributions of column density (\avpdf s) for which each column density is linked to a volume density as derived by hydrodynamical simulations. We consider two lognormal {\avpdf}s: a diffuse, low density medium with average visual extinction of $\overline{{\rm A}_V}=0.75\,{\rm mag}$ and dispersion of $\sigma=0.5$ and a denser giant molecular cloud with $\overline{{\rm A}_V}=4\,{\rm mag}$ and $\sigma=0.8$. We treat the UV radiation field, cosmic-ray ionization rate and metallicity as free parameters. We find that the low density medium remains fully HI- and CII-dominated under all explored conditions. The denser cloud remains almost always molecular (i.e. H$_2$-dominated) while its carbon phase (CO, CI and CII) is sensitive to the above free parameters, implying that existing methods of tracing H$_2$-rich gas may require adjustments depending on environment. Our numerical framework can be used to estimate the PDR properties of large ISM regions and quantify trends with different environmental parameters as it is fast, covers wide parameter space, and is flexible for extensions.

Tables

Table 1. Table 1: Initial gas-phase abundance of species used for the grid of PDR simulations. The gas-to-dust ratio is assumed to be 100.

H	$4.00\times {10}^{-1}$	C+	$1.00\times {10}^{-4}$
H2	$3.00\times {10}^{-1}$	O	$3.00\times {10}^{-4}$
He	$1.00\times {10}^{-1}$

Table 2. Table 2: Fractional abundances of species as calculated for the ‘Reference ISM’ defined in § 3.3 adopting A V subscript 𝐴 𝑉 A_{V} -PDF and A V − n H subscript 𝐴 𝑉 subscript 𝑛 H A_{V}{-}n_{\rm H}  relations described in § 3.2 . Columns 1 − 6 1 6 1{-}6 refer to the shapes of A V subscript 𝐴 𝑉 A_{V} -PDF, mean ( M ≡ A V ¯ 𝑀 ¯ subscript 𝐴 𝑉 M\equiv\overline{A_{V}} ), width ( σ 𝜎 \sigma ), median ( m 𝑚 m ), mode ( ℳ ℳ \cal M ) and A V − n H subscript 𝐴 𝑉 subscript 𝑛 H A_{V}{-}n_{\rm H} , respectively. Columns 7 − 11 7 11 7{-}11 correspond to the average fractional abundances of species of CO, C i , C ii , H i and H 2 .

$A_V$-PDF	$M$	$\sigma$	$m$	$ℳ$	$A_V-n_H$	$f_{\mathrm{CO}}$	$f_{\mathrm{CI}}$	$f_{\mathrm{CII}}$	$f_{\mathrm{HI}}$	$f_{H_2}$
						$\times {10}^{-5}$	$\times {10}^{-5}$	$\times {10}^{-5}$
Case-1	0.75	0.5	0.66	0.52	Upper bound	0.0336	0.898	9.094	0.694	0.153
					Average	0.0900	1.855	8.051	0.548	0.226
					Lower bound	0.1569	2.616	7.210	0.461	0.269
Case-2	4.0	0.8	2.90	1.53	Upper bound	4.471	2.981	2.536	0.114	0.442
					Average	5.317	2.842	1.827	0.080	0.460
					Lower bound	5.800	2.714	1.470	0.065	0.467

Equations

$${\mathrm{C}}^{+}+{\mathrm{H}}_2⟶{\mathrm{CH}}^{+}+\mathrm{H}.$$

$$\mathrm{PDF}(A_V;\mu ,\sigma )=\frac{1}{A_V\sigma \sqrt{2\pi }}\exp \left[-\frac{(\ln A_V-\mu {)}^2}{2{\sigma }^2}\right],$$

$$M\equiv \exp \left(\mu +\frac{{\sigma }^2}{2}\right)=\overline{A_V},$$

$$m\equiv \exp (\mu )=\overline{A_V}\exp \left(-\frac{{\sigma }^2}{2}\right),$$

$$\mathcal{M}\equiv \exp (\mu -{\sigma }^2)=\overline{A_V}\exp \left(-\frac{3{\sigma }^2}{2}\right).$$

$$f_{\mathrm{sp}}=\frac{{\sum }_{i=1}^q{N}_i(\mathrm{sp})\cdot {\mathrm{PDF}}_i}{{\sum }_{i=1}^q{N}_i(\mathrm{H},\mathrm{tot})\cdot {\mathrm{PDF}}_i},$$

$$⟨f_{\mathrm{sp}}⟩=\frac{{\sum }_j^{\ell }{f}_{\mathrm{sp},j}\cdot {\mathrm{PDF}}_j}{{\sum }_j^{\ell }{\mathrm{PDF}}_j},$$

$$⟨f_{\mathrm{sp}}⟩$$

$${\left[f_{\mathrm{sp}}\right]}_{n_{\mathrm{H}}}$$

$$A_V(n_{\mathrm{H}})=g(A_V,n_{\mathrm{H}})\left\{n_{\mathrm{H}}+\lambda n_{\mathrm{H}}\left(1-\frac{\log (n_{\mathrm{H}}/n_{\min })}{\log (n_{\max }/n_{\min })}\right)\right\}\left[\frac{\mathrm{mag}}{{\mathrm{cm}}^{-3}}\right]$$

Full text

††pubyear: 2018††pagerange: Simulating the atomic and molecular content of molecular clouds using probability distributions of physical parameters–B

1 Introduction

Astrophysical and cosmological calculations are nowadays becoming increasingly advanced, with the general trend to couple the (magneto-) hydrodynamical evolution of the interstellar medium (ISM) with detailed chemistry in all gas phases. The goal of these new codes is to understand the evolution of molecular clouds and how this links to galaxy evolution (e.g. Inoue & Inutsuka, 2012; Walch et al., 2015; Kim & Ostriker, 2017), to determine how metallicity, ultraviolet radiation and cosmic-rays affect the properties of the ISM (e.g. Girichidis et al., 2016; Richings & Schaye, 2016a, b), to identify H2-rich but CO-poor gas in galaxies (e.g. Clark et al., 2012; Smith et al., 2014; Bournaud et al., 2015), and to study the star formation process (e.g. Hu et al., 2017). Although such codes can simulate the ISM to the point of a direct comparison with observations, they are computationally very demanding and expensive, requiring often an interdisciplinary synergy of groups.

When it comes to chemistry, a key challenge of such hydro-chemical codes is to determine as detailed as possible the atomic-to-molecular transition and the transition of the different carbon phases (e.g. Glover et al., 2010; Offner et al., 2013). This problem has been the main focus of photodissociation region (PDR) studies over the past 30 years or so (Hollenbach & Tielens, 1999; Röllig et al., 2007). In general, PDR codes do not attempt to co-evolve hydrodynamically the ISM. Instead, they apply detailed micro-physics that are sensitive to many user-defined parameters and that together control the overall chemistry, excitation and thermal balance (e.g. van Dishoeck & Black, 1988; Ferland et al., 1998; Bell et al., 2006; Le Petit et al., 2006; Bisbas et al., 2012). These codes are able to estimate the abundances of hundreds of species and the line emission of various coolants. This information can be used to create synthetic images to be contrasted against real observations (see Haworth et al., 2018, for a review). However, a fundamental limitation of these codes is that they have to adopt many ad hoc user-defined parameters describing the ISM, such as the density distribution, gas velocity, the structure of the ultraviolet radiation field and others. These ‘free’ parameters play a crucial role when determining the abundance and emissivities of species and thus affect the physical interpretation of observations.

Coupling hydrodynamics with detailed chemistry has the potential to simulate the ISM quite realistically. Indeed, first promising efforts have been made towards this direction (Bisbas et al., 2015b; Motoyama et al., 2015; Haworth et al., 2015), although the computational expense is substantial. To reduce the computational cost of hydro-chemical simulations, simpler (and therefore faster) chemical networks have been implemented (Nelson & Langer, 1997; Glover et al., 2010; Grassi et al., 2014). These networks provide an estimate of the local gas temperature, leading to a more accurate estimation of the thermal pressure and therefore the overall ISM hydrodynamics. They are nowadays frequently used to study isolated turbulent molecular clouds in different conditions (from the solar neighbourhood to the Galactic Centre and the distant Universe) and star formation therein when only the most important chemical reactions are accounted for. However, for now extensive parameter studies or detailed simulations of entire galaxies remain prohibitive due to their computational cost.

One important task of chemical calculations is to determine how commonly used methods for tracing H2-rich gas, e.g. via CO, depend on various factors. The most important of them are (i) the average ISM metallicity (Pak et al., 1998), (ii) the ambient FUV radiation field (Pak et al., 1998; Wolfire et al., 2003), (iii) the average temperature, density and dynamical state of the gas (Dickman et al., 1986; Young & Scoville, 1991; Bryant & Scoville, 1996; Papadopoulos et al., 2012) and (iv) the CO destruction due to cosmic rays (Bisbas et al., 2015a, 2017a). These factors are even more important in studies of low-metallicity or high-$z$ galaxies due to their extreme ISM conditions. We thus require predictions of how the partitioning between the carbon phase scales with critical ‘environmental’ parameters.

To connect hydrodynamical simulations, chemical modeling and observations, the visual extinction, $A_{V}$, probability density function ($A_{V}$-PDF) is considered. At solar metallicity, $A_{V}$ is related to the total H-nucleus column density, $N_{\rm H}$, through $A_{V}=N_{\rm H}\cdot 6.3\times 10^{-22}$ (Weingartner & Draine 2001; Röllig et al. 2007; see also Bohlin et al. 1978; Rachford et al. 2009). An $A_{V}$-PDF is defined as the probability of finding ISM gas within a visual extinction in the range of [$A_{V},\,A_{V}+{\rm d}A_{V}$]. Various observational studies have determined $A_{V}$-PDFs of Galactic molecular clouds (e.g. Goodman et al., 2009; Kainulainen et al., 2009; Froebrich & Rowles, 2010; Abreu-Vicente et al., 2015; Schneider et al., 2015a, b, 2016), the main characteristic of them being that they consist of a log-normal component for $A_{V}/\overline{A_{V}}\lesssim 1$ (where $\overline{A_{V}}$ is the mean observed visual extinction) and a power-law component for $A_{V}/\overline{A_{V}}\gtrsim 1$. Several theoretical studies explored the nature of these two components (Tassis et al., 2010; Kritsuk et al., 2011; Burkhart et al., 2013; Girichidis et al., 2014; Burkhart et al., 2015) and concluded that the power-law tail is a result of self-gravity possibly leading to star formation. Other groups have suggested that the power-law tail may be a result of the external pressure to gravitationally unbound entities (Kainulainen et al., 2011), that $A_{V}$-PDFs of collapsing protoclusters may not exhibit a power-law tail at all (Butler et al., 2014) or even that an $A_{V}$-PDF may be a superposition of various log-normal components (Brunt, 2015). On the other hand, recent studies propose that the distinction between log-normal and power-law components may naturally result from statistical errors connected with the consideration of the last (outer) closed contour (Lombardi et al., 2015; Ossenkopf-Okada et al., 2016; Alves et al., 2017; Körtgen et al., 2018) and that all PDFs are in fact power-laws.

Whatever their precise form, it is these $A_{V}$-PDFs that can be used as inputs to estimate the astrochemical properties of the ISM with much simpler models. Indeed, such an approach has been used to connect ensembles of column densities with spectral energy distributions of different cooling lines (e.g. Krumholz & Thompson, 2007; Narayanan et al., 2008). One such example is the despotic code (Krumholz, 2014) which uses one-zone models to represent optically thick molecular clouds and calculates line luminosities and cooling rates based on an escape probability formalism. Complementary, Leroy et al. (2017) use radex (van der Tak et al., 2007) to model an ensemble of (one-zone) columns to analyse line emissions in extragalactic observations. These studies move in the direction of constructing a general framework to estimate average ISM properties accounting for complex (unresolved) substructure which is particularly useful when dealing with observations of the turbulent ISM over large spatial scales. However, the previous methods lack of detailed PDR chemistry which could potentially provide better insights into the physical state of the observed ISM.

The motivation of the present work is to further improve this framework to estimate the atomic and molecular content of ISM regions from their substructure by using more detailed and self-consistent thermo-chemical calculations. The proposed method considers user-specified $A_{V}$-PDFs as input to parametrise the ISM at large scales and links them with a grid of pre-run 1D thermo-chemical calculations of uniform-density providing average PDR properties, such as abundances of species and line emissivities. To do this, each $A_{V}$ value is connected to a most probable local H-nucleus number density, $n_{\rm H}$ (and therefore a most probable size of a uniform-density sphere) as determined by hydrodynamical simulations. Each individual PDR is integrated over depth up to this $A_{V}$ value (also sometimes called $A_{V}^{\rm tot}$, see van Dishoeck & Black 1988). The $A_{V}$-PDF is the key element to perform chemical modeling of observed ISM regions with complex substructure, without the necessity of performing time-consuming 3D hydro-chemical calculations, and to be able to assess quickly how the ISM characteristics change with critical parameters such as the radiation field ($\chi$), cosmic ray ionization rate ($\zeta_{\rm CR}$) and the metallicity ($Z$). The above is of particular interest in studies of unresolved systems such as extragalactic objects.

In this paper we focus on how the abundances of hydrogen (H2, Hi) and carbon phase (CO, Ci, Cii) depend on $\chi$, $\zeta_{\rm CR}$ and $Z$ in ISM distributions corresponding to ambient low density atomic gas in galaxies and denser clouds that make up Giant Molecular Clouds. These media are parametrised by lognormal $A_{V}$-PDFs (see §3.2). The low density medium has mean extinction along the line-of-sight of the observer of $\overline{A_{V}}=0.75$ mag, a dispersion of $\sigma=0.5$ in logarithmic intervals of $\overline{A_{V}}$ and mean hydrogen density of $n_{\rm H}\sim 80$ cm*-3*. We will find that such a medium consists predominantly of Hi and Cii (Figure 1). The denser molecular cloud has $\overline{A_{V}}=4$ mag, $\sigma=0.8$, $n_{\rm H}\sim 600$ cm*-3* and consists of H2 and a mixture of CO, Ci, Cii.

The paper is organized as follows. Section 2 gives a brief overview on the carbon cycle and the atomic-to-molecular transition. Section 3 describes the numerical methodology including the PDR calculations. Section 4 presents the results of our approach. Section 5 examines how various other approaches compare to our presented one. Section 6 discusses the limitations of this methodology and provides directions for further improvements. We conclude in Section 7. Appendix A presents how the atomic-to-molecular transition is affected due to the consideration of the suprathermal formation of CH+, which in turn leads to enhanced CO at low extinctions. Appendix B presents the here adopted relation between the local volume density and the effective (local) visual extinction.

In a follow-up paper we will consider the radiative transfer to translate the here derived abundances and excitations to observable emission line fluxes. However, as abundances and excitations are the key elements of modeling the ISM, we will concentrate in this first paper on a detailed investigation of their properties as function of environmental conditions.

2 The carbon cycle and atomic-to-molecular transition

Figure 1 illustrates a PDR structure illuminated from one side for conditions typically found in the Milky Way. In order to understand the change in thermal and chemical structure of the ISM, we need to understand when each chemical element transitions from one phase to another. The atomic-to-molecular (Hi-to-H2) transition is vital for understanding the evolution of the ISM and its chemistry (Jura, 1974; Black & Dalgarno, 1977; van Dishoeck & Black, 1986; Glover et al., 2010; Offner et al., 2013; Sternberg et al., 2014; Bialy et al., 2015b). This transition is determined and controlled by many different parameters (e.g. Bialy et al., 2017), the most important of which are the far-ultraviolet radiation (FUV; Jura, 1974; Hollenbach & Tielens, 1999), the cosmic-ray ionization rate (e.g. Meijerink et al., 2011), and the gas phase metallicity (e.g. Schruba et al., 2018). Other parameters such as shock (mechanical) heating (Meijerink et al., 2011) and X-ray heating (Maloney et al., 1996; Meijerink et al., 2006) may play significant role in special environments, however these are not examined in this work.

Carbon is the most important species for observing the properties of the ISM. It can be found in three major phases: ionized (Cii), atomic (Ci) and molecular in the form of carbon monoxide (CO). All three species are major coolants in the ISM through their spectral line emission (Hollenbach & Tielens, 1999). CO is the most abundant molecule after H2 and is frequently used to trace cold H2 gas (Dickman et al., 1986; Solomon et al., 1987; Bolatto et al., 2013), since H2 is a very poor emitter in most ISM conditions. The carbon phase pathway (Cii/Ci/CO) is connected with different ISM evolutionary stages (e.g. van Dishoeck & Black, 1988; Sternberg & Dalgarno, 1995; Beuther et al., 2014; Walch et al., 2015), and studying it is thus of significant importance.

In the initial stages of the formation of clouds, the ISM is found in a diffuse and ionized state, thus rich in Cii and bright in the [Cii] $158\,\mu$m fine-structure line (Röllig et al., 2006; Nordon & Sternberg, 2016). Here, hydrogen is mainly found in atomic form (Hi), unless a source of Lyman-$\alpha$ photons is present ($h\nu>13.6\,{\rm eV}$, i.e. an OB star or cluster of massive stars) which ionizes the ISM and creates Hii regions. Turbulence and self-gravity can increase locally the density of the ISM leading to larger total column densities (Elmegreen & Scalo, 2004; Hennebelle & Falgarone, 2012). This creates conditions in which Cii recombines forming Ci and subsequently CO. Due to its higher abundance and self-shielding, hydrogen is transformed earlier than carbon (i.e. at lower visual extinctions) and forms H2 molecules. CO can also self-shield, but to a lesser degree than H2, so dust absorption is also important. Thus, if the ISM density increases even more, to conditions typically found in giant molecular clouds, the dust shields the propagation of the FUV ($6<h\nu<13.6\,{\rm eV}$) radiation and this creates vast regions of molecular gas rich in CO and other molecules (e.g. van Dishoeck & Black, 1988; Bergin et al., 2004; Langer et al., 2010; Glover & Mac Low, 2011).

Cosmic-rays penetrate the ISM even at high column densities and change their chemistry in places where (external) FUV radiation cannot reach (see Strong et al., 2007; Grenier et al., 2015, for a review). The effect of elevated cosmic-ray energy densities on the chemistry related to carbon phases and the atomic-to-molecular transition was studied numerically by Meijerink et al. (2011) and analytically by Bialy & Sternberg (2015). Bisbas et al. (2015a, 2017a) investigated in detail the consequences on the traceability of molecular gas in star-forming galaxies at low and high redshift, since higher $\zeta_{\rm CR}$ values scale proportionally to the star formation rate (Papadopoulos, 2010). As the $\zeta_{\rm CR}$ ionization rate increases, cosmic-rays produce a high amount of He+ which then reacts with CO creating Cii. The latter then quickly recombines with free electrons creating Ci. At the same time, the H2 molecule remains remarkably unaffected for densities expected in diffuse molecular clouds. This effect of cosmic-ray induced CO destruction has the potential to create vast amounts of CO-poor H2 gas (Bisbas et al., 2015a), much higher than those created due to the presence of FUV radiation. The term ‘CO-poor’ is used to describe the molecular gas state in which the CO abundance with respect to H2 is significantly lower than normally expected, whether caused by enhanced cosmic rays or any other parameter (van Dishoeck, 1992).

Low metallicity gas as typically found in dwarf galaxies (e.g. Stark et al., 1997; Tafelmeyer et al., 2010; Requena-Torres et al., 2016; Pineda et al., 2017), the outer parts of large galaxies ($R\gtrsim 15\,{\rm kpc}$, Balser et al., 2011; Hayden et al., 2014) and galaxies at high redshift (e.g. Carniani et al., 2018) have low dust-to-gas ratios which reduce the H2 formation rate (which happens on the surfaces of dust grains) as well as the shielding of CO photodissociation by dust. In addition, the reduced C and O abundances imply that the CO column builds up more slowly so that its self-shielding sets in later. This has profound effects on all chemical processes of these systems (Maloney & Black, 1988; Wolfire et al., 1995; Pak et al., 1998; Madden et al., 2006; Bialy & Sternberg, 2015) and modifies the stratification of PDR species (Fig. 1). Observations of low metallicity clouds show that CO is restricted to small, well shielded dense gas peaks surrounded by Hi- and H2-rich gas (e.g. Schruba et al., 2017), and that CO directly scales with the line-of-sight extinction $A_{V}$ (Lee et al., 2015, 2018).

3 Numerical method

3.1 PDR astrochemical calculations

For the purposes of this work, the publicly available code 3d-pdr111https://uclchem.github.io/3dpdr.html is used, which treats the astrochemistry of photodissociation regions (Bisbas et al., 2012). Although the code is able to perform PDR calculations for 3D density distributions, here we apply it to ensembles of 1D columns. The gas column is illuminated from one side by a plane-parallel ultraviolet (UV) radiation field $\chi$ (normalized to the spectral shape of Draine, 1978) with UV radiation decreasing with depth into the cloud. We consider UV photons belonging to the Lyman-Werner band ($11.2-13.6\,{\rm eV}$) which can dissociate molecules but do not ionize hydrogen. The attenuated UV radiation is a function of the optical depth $\tau$ and is given by $\chi=\chi_{0}\exp(-\omega_{\lambda}\tau)$ with $\omega_{\lambda}=3.02$ at $\lambda=1000\,\textup{\AA}$ (see Röllig et al., 2007). The code performs full thermal calculations by balancing various heating and cooling processes as function of depth along the column, and adopting an escape probability approach determines the local gas temperature as a self-consistent solution (see Bisbas et al., 2012, 2017a, for full details). The abundances of species and their line emissivities are then calculated for the local (attenuated) UV field and gas density. Here, a subset of the UMIST2012 network is used (McElroy et al., 2013) consisting of 33 species and 330 reactions.

In this work, CO freeze-out on dust grains is not included. However, the route of suprathermal formation of CO via CH+ is included, following the methodology described by Visser et al. (2009). This route is found to give a better agreement between the column densities of H2 and CO at low visual extinctions. Federman et al. (1996) argued that the observed enhancement of CO column densities at low visual extinctions ($A_{V}<1\,{\rm mag}$) is possibly due to non-thermal motions between ions and neutrals as a result of Alfvén waves interacting with the cloud in its outermost parts, thus helping to overcome the energy barrier of the reaction

[TABLE]

This interaction increases the gas (kinetic) temperature by the amount of $\mu v_{A}^{2}/3k_{\rm B}$, where $\mu$ is the reduced mass of reactants from all ion-neutral reactions, $v_{A}$ is the Alfvén speed and $k_{\rm B}$ is the Boltzmann constant. This increment in gas temperature is able to change the pathway normally followed to form CO, the latter being formed now due to reactions between CH+ and O in addition to C+ and OH (Federman et al., 1996). Interestingly, such a turnover in the formation of CO via the OH channel has also been found to result from elevated cosmic-ray energy densities (Bisbas et al., 2017a; Bialy & Sternberg, 2015). This means that cosmic-rays can initiate an alternative chemical pathway of CO formation. In addition, the above reaction shifts the Hi-to-H2 transition at higher $A_{V}$, since it dissociates H2 in the outermost parts of the cloud where it could otherwise remain molecular, i.e. if the energy barrier cannot be overcome (see Appendix A).

To demonstrate how the suprathermal formation of CO via CH+ operates, we have performed a full thermochemical PDR calculation of an one-dimensional uniform density distribution with $n_{\rm H}=100\,{\rm cm}^{-3}$ interacting with a plane-parallel radiation field with strength $\chi/\chi_{0}=1$. Here, the cosmic-ray ionization rate is taken to be $\zeta_{\rm CR}=1.3\times 10^{-17}\,{\rm s}^{-1}$, the Alfvén speed is taken to be $v_{A}=3.3\,{\rm km}\,{\rm s}^{-1}$ and the calculation is terminated once steady-state thermal balance has been reached. The initial gas phase elemental abundances considered (normalized to hydrogen), where He $=0.1$, C${}^{+}=1.4\times 10^{-4}$ and O $=3.1\times 10^{-4}$. All aforementioned values are chosen to allow for a direct comparison with the results of Visser et al. (2009).

Figure 2 shows the resultant correlation of $N$(H2) and $N$(CO) when the suprathermal formation of CO via CH+ has been considered (green solid line) and when it has not (dashed line). As can be seen, $N$(CO) is approximately one order of magnitude higher in the first case when $N$(H2)$\lesssim 10^{21}\,{\rm cm}^{-2}$, the latter of which corresponds to a visual extinction of $A_{V}\sim 1.5\,{\rm mag}$. The red squares are observations taken by Sheffer et al. (2008) in Galactic diffuse molecular clouds. The blue solid line is the correlation estimated by Visser et al. (2009) obtained by averaging 100 isothermal PDR calculations. Since a significant fraction of the ISM can be in this low column density phase, it is important to treat this region accurately.

3.2 Abundances of species for given $A_{V}$-PDFs

The aim of this work is to construct a method that determines the average chemical properties of the ISM, i.e. the abundances and excitation of species, by using as input probability density functions (PDFs) of the main physical parameters that describe interstellar clouds. The most fundamental parameter for this work is the total hydrogen column density $N_{\rm H}$-PDF. Hereafter, this parameter will be referred to with the more commonly used relation of $A_{V}$-PDF, where $A_{V}$ is the extinction related to $N_{\rm H}$ through $A_{V}=N_{\rm H}\cdot 6.3\times 10^{-22}$ for solar metallicity galactic clouds (Weingartner & Draine, 2001; Röllig et al., 2007). When modelling lower metallicity regimes, the same $N_{\rm H}$-PDF is adopted but the associated $A_{V}$ varies proportional to the metallicity, $Z$. In our 1D-PDR models, the visual extinction in the $A_{V}$-PDF function is taken to be the extinction integrated along the line-of-sight of the observer. The latter is also denoted as $A_{V}^{\rm tot}$.

As a first approach towards this aim, simple lognormal $A_{V}$-PDF functions of the following form are considered:

[TABLE]

where $\mu\equiv\overline{\ln A_{V}}$ and $\sigma$ are the mean and standard deviation of the variable’s natural logarithm, respectively. Three different statistical characteristic values can be then defined. The mean value, $M$, which corresponds to the average value of all $A_{V}$ (i.e. $\overline{A_{V}}$), defined as

[TABLE]

the median value, $m$, which corresponds to the value of $A_{V}$ that divides the PDF area in two equal parts, defined as

[TABLE]

and the mode, ${\cal M}$, which corresponds to the peak value of the PDF (where $dp/dA_{V}=0$), defined as

[TABLE]

By adopting appropriate values of $\overline{A_{V}}$ and $\sigma$, two different $A_{V}$-PDFs are created corresponding to a hypothetical low density, atomic-dominated medium (with $\overline{A_{V}}=M=0.75\>{\rm mag}$ and $\sigma=0.5$ giving $m=0.66$ and ${\cal M}=0.52$; hereafter ‘Case-1’) and a denser giant molecular cloud (with $\overline{A_{V}}=M=4\>{\rm mag}$ and $\sigma=0.8$ giving $m=2.90$ and ${\cal M}=1.53$; hereafter ‘Case-2’). We note that we do not attempt to model any specific region reported in observations and that we ignore any power-law tail that may appear in a region similar to the Case-2 (although its contribution can be considered by the presented method, as discussed in §4.5). Instead, the mean and the width are chosen to roughly mimic typical structures of the ISM to test the proposed method and determine trends with physical parameters. Such PDFs have been found in the local ISM (e.g. Kainulainen et al., 2009). The shapes of these PDFs are illustrated in Fig. 3.

To estimate the average column-integrated fractional abundances, $f_{\rm sp}$, of each different species (sp), the following formula is applied

[TABLE]

where $N_{i}$ is the column density multiplied by the frequency ${\rm PDF}_{i}$, which is determined by the $A_{V}$-PDF function given by Eqn. 2 and $N({\rm H,tot})=N({\rm H\textsc{i}})+2N({\rm H_{2}})$ is the total hydrogen column density. The summation over the index $i$ considers the PDR calculations at a linear grid of $A_{V}$ values.

To estimate the column density, $N$, a relation connecting the H-nucleus number density, $n_{\rm H}$, with the most probable value of visual extinction, $A_{V}$, is required. It is reasonable to expect that low $n_{\rm H}$ are most likely found at low $A_{V}$, therefore being more strongly affected by UV radiation. On the contrary, high $n_{\rm H}$ are expected at high $A_{V}$. Three-dimensional hydrodynamical simulations from parsec (Glover et al., 2010; Offner et al., 2013) to kiloparsec scales (Van Loo et al., 2013; Seifried et al., 2017; Safranek-Shrader et al., 2017; Gong et al., 2018) remarkably demonstrate this correlation. In this work, the ‘default’ $A_{V}{-}n_{\rm H}$ relation is an average over four different such relations found in the literature (see Appendix B for details). In the upper panel of Fig. 4, the default $A_{V}{-}n_{\rm H}$ relation is plotted as solid line along with the adopted upper (dot-dashed line) and lower bounds (dashed line) to mimic the relations derived from hydrodynamical simulations. In the following, each upper and lower bound is treated independently as an additional $A_{V}{-}n_{\rm H}$ relation. This is to demonstrate how the estimated abundances of species depend on the choice of such relation.

By relating the line of sight integrated visual extinction, $A_{V}$, with the number density $n_{\rm H}$, we are in fact connecting each cloud with a specific uniform-density PDR calculation. Strictly speaking, this is an approximation since $A_{V}$ is a local value in the simulations used to derive this relation. Our constant density clouds are used as a simplified proof of concept for the statistics presented here. A more realistic assumption would be to use a given $A_{V}{-}n_{\rm H}$ relation as an input directly to the PDR code, so that the edge of the cloud has a lower density than the center, and calculate the thermal balance and abundances of species accordingly. The latter method is to be developed in a subsequent paper.

The $A_{V}{-}n_{\rm H}$ relations described above imply an average size, $L$, corresponding to the ‘depth’ or the ‘line-of-sight extent’ of a hypothetical uniform density sphere of $L=A_{V}/6.3\times 10^{-22}n_{\rm H}$. The bottom panel of Fig. 4 shows the scaling between $L$ and $n_{\rm H}$. The mode of the ‘Case-1’ $A_{V}$-PDF corresponds to a medium with $n_{\rm H}\sim 20\,{\rm cm}^{-3}$, $L\sim 13\,{\rm pc}$ and $M\sim 4.5\times 10^{3}\,{\rm M}_{\odot}$ and the mode of the ‘Case-2’ corresponds to $n_{\rm H}\sim 600\,{\rm cm}^{-3}$, $L\sim 1.3\,{\rm pc}$ and $M\sim 136\,{\rm M}_{\odot}$.

In our analysis we consider two density limits: a minimum of $n_{\rm H}=1\,{\rm cm}^{-3}$ and a maximum of $n_{\rm H}=10^{6}\,{\rm cm}^{-3}$, below and above which no PDR modelling is performed. Regions corresponding to densities (and thus $A_{V}$ values) outside these bounds are shown as grey shaded regions in Figs. 3 and Fig. 4. For these regions we adopt the 1D-PDR calculation performed at the respective density limit and integrate the PDR model to the total $A_{V}$ of that cloud. For $A_{V}>30\,{\rm mag}$ (and consequently for the models with $n_{\rm H}\gtrsim 10^{5}\,{\rm cm}^{-3}$), we adopt the chemical composition calculated at the maximum extinction of $A_{V}=30\,{\rm mag}$ but integrate the PDR models to the $A_{V}$ as dictated by the $A_{V}{-}n_{\rm H}$ relation. This inserts negligible errors both due to the fact that the UV radiation has been severely attenuated and thus chemistry remains unchanged, and also because the probability of such a high $A_{V}$ is typically negligible ($\lesssim 10^{-4}$ for the Case-2 $A_{V}$-PDF).

A large grid of constant density simulations is then constructed (see §4) and a look-up table is used to recall a particular PDR simulation and to calculate the column densities of CO, Ci, Cii, Hi and H2 (the species examined in this paper), as well as the total H-nucleus column density, up to that value of visual extinction. Equation 6 is then applied for both $A_{V}$-PDF cases. Chemical iterations are converged once equilibrium is reached. 3d-pdr performs these calculations for a chemical time of $8\,{\rm Myr}$, although equilibrium is typically reached after the first couple Myrs.

3.3 Abundances of individual uniform-density PDRs

As reference for our chemical modeling, we consider the following set of one-dimensional, uniform density PDR calculations. We consider an H-nucleus number density that spans 6 orders of magnitude222Densities higher than $10^{6}\,{\rm cm}^{-3}$ correspond to very small probabilities in both PDF cases explored and, as we consider only relatively low critical density molecules, are not considered further. ($n_{\rm H}=10^{0}-10^{6}\,{\rm cm}^{-3}$) at a resolution of $0.01$ per logarithmic dex (therefore $q=600$) and assume a maximum visual extinction of $A_{V}^{\rm max}=30\,{\rm mag}$ at all times. The other ‘environmental’ parameters are kept fixed at a plane-parallel FUV radiation field with strength $\chi/\chi_{0}=1$ (Draine, 1978), a cosmic-ray ionization rate per H atom of $\zeta_{\rm CR}=10^{-16}\,{\rm s}^{-1}$ representing the average value observed in the Milky Way (Dalgarno, 2006; Indriolo et al., 2015; Neufeld & Wolfire, 2017), and adopting solar metallicity. The column density of each species is derived by integrating the abundance at each depth into the cloud up to the $A_{V}^{\rm tot}$ as sampled by the $A_{V}$-PDF of either Case-1 or Case-2. The initial gas-phase chemical abundances are listed in Table 1. These abundances take into account that some fraction of carbon and oxygen is locked up in solids.

Figure 5 shows the column densities of Hi and H2 (upper panel) and carbon phase (lower panel) for individual uniform density, one-dimensional PDRs as a function of $A_{V}^{\rm tot}$ when adopting the ‘average’ $A_{V}{-}n_{\rm H}$ relation (solid lines) or the upper or lower bounds (dot-dashed or dashed lines). For these PDRs, the atomic-to-molecular transition occurs at $A_{V}\sim 0.9\,{\rm mag}$ (we note that we define this transition by $N_{\textsc{Hi}}=N_{\rm H_{2}}$) which corresponds to a column with density $n_{\rm H}\sim 130\,{\rm cm}^{-3}$ and length $L\sim 3.2\,{\rm pc}$. The PDRs are Cii-dominated for $A_{V}^{\rm tot}\lesssim 1.5\,{\rm mag}$ (with $n_{\rm H}\sim 600\,{\rm cm}^{-3}$ and $L\sim 1.3\,{\rm pc}$), Ci-dominated for $1.5\lesssim A_{V}^{\rm tot}\lesssim 3.2\,{\rm mag}$ (with $600\lesssim n_{\rm H}\lesssim 3700\,{\rm cm}^{-3}$ and $0.4\lesssim L\lesssim 1.3\,{\rm pc}$) and CO-dominated for $A_{V}^{\rm tot}\gtrsim 3.2\,{\rm mag}$ ($n_{\rm H}\gtrsim 3700\,{\rm cm}^{-3}$ and $L\lesssim 0.4\,{\rm pc}$).

3.4 Convergence test

The accuracy of the method is examined by performing a convergence test for the $q$ number of total PDR simulations needed to sample a given $A_{V}$-PDF. First, we lower this number by one half ($q=300$), implying a 0.02 dex logarithmic resolution in the $n_{\rm H}$ space. Second, we double the initial $q$ number ($q=1200$), corresponding to 0.005 dex resolution, at which point additional PDR simulations are run.

In all cases, the changes in abundance stated in Table 2 (see §4) between $q=300$ and $q=1200$ are always $\lesssim 2\%$, and between $q=600$ and $q=1200$ are $\lesssim 1\%$, indicating convergence. It is interesting to note that these changes are related with the different carbon phases, while variations in the Hi and H2 abundances for different $q$ are negligible.

4 Results

4.1 Abundances for reference environmental parameters

Table 2 lists the derived fractional abundances of species, $f_{\rm sp}$, for the two different $A_{V}$-PDFs considered here for the constant ‘environmental’ parameters defined in §3.3. The Case-1 $A_{V}$-PDF describes a diffuse, low density medium which is found to be rich in Hi and Cii. When varying the $A_{V}{-}n_{\rm H}$ relation within its bounds, the abundance of the more rare species CO, Ci, and H2 are affected at a factor of $\lesssim 2{-}3$ level. The Case-2 $A_{V}$-PDF describes a denser molecular cloud and is found to be rich in H2 and contains a mix of Cii, Ci and CO. We note that in all cases $f_{\textrm{CO}}+f_{\textsc{Ci}}+f_{\textsc{Cii}}\approx 10^{-4}$, which is the carbon abundance adopted in the PDR simulations (Table 1), and $f_{\textsc{Hi}}+2f_{\rm H_{2}}\approx 1$. Within the boundary of the $A_{V}{-}n_{\rm H}$ relation, the abundances of the rare species Hi and Cii change at a factor $\lesssim 1.5$ level.

With this ‘reference’ as a starting point, three suites of PDR calculations are performed in which the incident FUV radiation field $\chi$, the cosmic ray ionization rate $\zeta_{\rm CR}$, and the metallicity $Z$ are varied as free parameters. In particular, we consider four different intensities of the incident FUV radiation field ($\chi/\chi_{0}=10^{0}-10^{3}$), four different cosmic-ray ionization rates ($\zeta_{\rm CR}=10^{-17}-10^{-14}\,{\rm s}^{-1}$), and four different metallicities ($Z/Z_{\odot}=0.1,\,0.2,\,0.5,\,1$). Figure 6 shows the results for the Case-1 $A_{V}$-PDF and Figure 7 those for Case-2.

4.2 Varying the FUV radiation field intensity

In the first suite of calculations, the external FUV radiation field $\chi$ is varied as free parameter while $\zeta_{\rm CR}=10^{-16}\,{\rm s}^{-1}$ and $Z=1\,Z_{\odot}$ are kept constant. The left columns of Figs. 6 and 7 show the average abundance of Hi, H2, Cii, Ci and CO as a function of $\chi/\chi_{0}$ for the Case-1 and Case-2 PDF models, respectively. Since the Case-1 PDF corresponds to a diffuse, low density medium, it is expected to be predominantly in atomic form with high abundances of Cii. Indeed, as can be seen in Fig. 6, for $\chi/\chi_{0}=1$, $f_{\textsc{Hi}}$ is already $\gtrsim 0.5$ and $f_{\textsc{Cii}}\gtrsim 8\times 10^{-5}$ (which is $\sim 80\%$ of carbon) for all $A_{V}{-}n_{\rm H}$ relations. The abundances of H2, Ci and CO differ by a factor between two to four as a result of the adopted upper and lower $A_{V}{-}n_{\rm H}$ relations for the low $A_{V}$ regime. However, the abundances of these species are already much smaller when compared to the more dominant Hi and Cii, respectively, and such a difference can therefore be considered negligible. H2, Ci and CO decrease monotonically with increasing FUV radiation field intensity. For the case of $\chi/\chi_{0}=10^{3}$, this low density medium is Hi- and Cii-dominated by $\gtrsim 99\%$.

The much denser molecular cloud described by the Case-2 PDF (left column of Fig. 7) is dominated by H2 for FUV radiation fields up to $\chi/\chi_{0}\approx 300$, above which it becomes Hi-dominated. For $\chi/\chi_{0}\lesssim 7$ the cloud is dominated by CO, above which is dominated by Cii. It is interesting to note that the cloud never becomes Ci-dominated as the FUV field increases. Instead, it is predicted to be dominated by warm gas rich in H2 and Cii. The largest fraction of Ci is $\sim 2.8\times 10^{-5}$ and for $\chi/\chi_{0}=1$. Due to the overall high density, the shaded region due to the decreasing $A_{V}{-}n_{\rm H}$ bounds is thinner than in Case-1 described above.

For the Case-1 PDF, the average abundances of H2 and of the carbon phases change considerably with varying the radiation field. This is because the low density medium is translucent at low values of $A_{\rm V}$, and it implies drastic changes in the abundance ratios of CO/H2, Ci/H2 and Cii/H2. For the Case-2 PDF, none of the abundances experience such large changes and the abundance of H2 scales closely with CO and Ci (but not Cii).

4.3 Varying the cosmic-ray ionization rate

As a next application, we study the change in abundance of key species as a function of $\zeta_{\rm CR}$. The middle columns of Figs. 6 and 7 show the results of these calculations. For the Case-1 PDF, the abundance of Cii remains very high regardless of $\zeta_{\rm CR}$. Likewise, Hi increases and this makes the low density medium to be $\gtrsim 70\%$ in the atomic phase, particularly for $\zeta_{\rm CR}\gtrsim 10^{-15}\,{\rm s}^{-1}$. For $\zeta_{\rm CR}\approx 10^{-15}\,{\rm s}^{-1}$, CO has a local maximum (although its abundance is overall very low), which is most likely attributed to the formation of CO via the OH channel as described in Bisbas et al. (2017a). It is interesting to note that Ci remains nearly constant as $\zeta_{\rm CR}$ increases and that Ci is always one to two orders of magnitude more abundant than CO.

In contrast, the denser molecular cloud described by the Case-2 PDF remains predominantly molecular for all $\zeta_{\rm CR}$. The carbon phase sensitively depends on the cosmic-ray ionization rate. Above $\zeta_{\rm CR}\gtrsim 3\times 10^{-16}\,{\rm s}^{-1}$, CO is effectively destroyed and carbon is either in Ci or Cii. For $3\times 10^{-16}<\zeta_{\rm CR}<5\times 10^{-15}\,{\rm s}^{-1}$ the carbon is in Ci, while for higher $\zeta_{\rm CR}$ it is in Cii. As already noted by Bisbas et al. (2015a), for moderately enhanced cosmic-ray ionization rates, Ci becomes a better tracer of H2-rich gas than CO. In this regime, Cii is also very abundant, however it may well originate from other heating sources such as FUV radiation (see §4.2), shocks or X-rays (Mackey et al., 2018) and, therefore, is not always connected to the presence of H2.

4.4 Varying the metallicity

As a third test, we investigate how the abundances of species change with lower metallicity ($Z$), while keeping FUV intensity and the cosmic-ray ionization rate fixed. In this work, we take the depths of clouds at low metallicity to be the same as those for solar metallicity for a given density.

The right columns of Figs. 6 and 7 show the corresponding calculations. The low density medium represented by the Case-1 PDF is always found to be atomic dominated and as $Z$ decreases the abundance of Hi increases. Almost all carbon is in form of Cii. As metallicity decreases, the abundances of all carbon phases also decrease in parallel and no transition from one phase to another is observed. At all metallicities, the fractional abundance of CO relative to hydrogen is very small ($\lesssim 5\times 10^{-7}$).

The denser molecular cloud (Case-2) remains fully molecular for all metallicities. At solar metallicity, the carbon is two to three times more in CO, than in Ci or Cii. For $Z<0.5\,Z_{\odot}$ it is dominated by Ci. For $Z\sim 0.1\,Z_{\odot}$, Ci has two times the abundance of Cii and six times the abundance of CO. This is because H2 self-shields from FUV radiation, whereas CO is shielded by dust, which is reduced in low metallicity environments. In addition, the column for CO self-shielding builds up deeper into cloud due to the overall lower abundances of C and O. Therefore, CO photodissociates creating a surplus of Ci and Cii (Maloney & Black, 1988; Pak et al., 1998; Bolatto et al., 1999).

4.5 Considering an $A_{V}$-PDF similar to Taurus

Lastly, we consider an $A_{V}$-PDF with $\overline{A_{\rm V}}=1.8\,{\rm mag}$ and width $\sigma=0.49$ which parametrises the column density distribution of the nearby Taurus star-forming cloud as observed by Kainulainen et al. (2009) at $\sim 0.1\,{\rm pc}$ resolution. For our ‘average’ $A_{V}{-}n_{\rm H}$ relation, this corresponds to a density of $n_{\rm H}\sim 300\,{\rm cm}^{-3}$ and a size of $L\sim 2\,{\rm pc}$. Furthermore, we consider the impact of the power-law tail that has been observed for this region (Kainulainen et al., 2009) on the overall chemical composition. The power-law tail is shown in the left panel of Fig. 8 (dashed line) and has a slope of about $-2.3$ (in the units of $dp/dA_{V}$) at $A_{V}>3\,{\rm mag}$. While we do not attempt to examine in detail the physical conditions prevailing in the ISM of Taurus, it is interesting to explore how the average abundances for an observed $A_{V}$-PDF change as a function of the free parameters considered in this work.

The middle and right columns of Fig. 8 show how the abundances scale with varying radiation field or cosmic ray ionization rates while keeping the metallicity fixed at solar. Solid lines along with their associated shadowed regions correspond to the log-normal component only without the contribution of the power-law tail. Dashed lines show results of the entire $A_{V}$-PDF when the power-law tail has been considered (for which the ‘average’ $A_{V}{-}n_{\rm H}$ relation has been taken into account only, to avoid confusion). From the $A_{V}$-PDF it is seen that this cloud remains dominated by H2 for moderate-to-high values of cosmic-ray ionization rates ($\zeta_{\rm CR}\lesssim 3\times 10^{-15}\,{\rm s}^{-1}$), but is CO-poor particularly for $\zeta_{\rm CR}>10^{-16}\,{\rm s}^{-1}$ (Fig. 8 middle column).

Cosmic-rays attenuate as a function of column density which has consequences in the CO/H2 abundance ratio. For example, Padovani et al. (2018) find that the cosmic-ray ionization rate for $A_{V}>4\,{\rm mag}$ may be $\zeta_{\rm CR}\lesssim 10^{-17}\,{\rm s}^{-1}$. Thus, the densest part of the Taurus star-forming cloud may have lower $\zeta_{\rm CR}$ values than those examined here. This in turn implies higher CO abundances than those predicted in Fig. 8. It is interesting to note that the molecular gas is equally rich in Ci and Cii until it transits to the atomic phase, where it becomes Cii-rich. When the power-law component is taken into consideration, it is found that the abundance of H2 increases by a small amount which is reflected by the small decrease in Hi abundance. The most striking feature, however, is the sudden increase in CO abundance for all $\zeta_{\rm CR}$ and $\chi/\chi_{0}$ values explored. In particular, for $\zeta_{\rm CR}\lesssim 3\times 10^{-17}\,{\rm s}^{-1}$ it is found that the molecular gas is CO-rich contrary to the findings of the log-normal component-only.

The cloud remains also H2-rich for FUV intensities up to $\chi/\chi_{0}\sim 10-40$, as shown in the right column of Fig. 8. For values of $\chi/\chi_{0}\sim 1$, Ci and Cii (and CO when the power-law tail is considered) have comparable abundances but for higher intensities it becomes (and stays) very Cii-rich. There is thus an extended range of $\zeta_{\rm CR}$ ($\sim 10^{-17}-10^{-15}\,{\rm s}^{-1}$) and $\chi/\chi_{0}$ ($\sim 1-10$) in which the abundance ratio of Cii/H2 is high. When compared to the Case-2 $A_{V}$-PDF (which has $\sim 2.2$ times higher $\overline{A_{V}}$ and $\sim 1.6$ times larger $\sigma$) Cii-dominated H2-rich gas is observed for a more extended range of $\chi/\chi_{0}$ ($\sim 10-300$), although the molecular gas is either CO-, Ci or Cii-dominated depending on the value of $\zeta_{\rm CR}$.

5 Difference to uniform density calculations

In this section, we investigate whether the abundances of lognormal $A_{V}$-PDFs can be adequately estimated by single uniform density PDR calculations as compared to the methodology presented in §3.2 of models that take into account the full distribution of extinction columns and associated volume densities. In particular, we consider lognormal $A_{V}$-PDFs and determine the chemical abundances at the mean, median and mode extinction columns (which are defined by Eqns. (3)–(5)). We keep the median constant ($m=3$ mag) and consider six different widths ($\sigma=0.01,0.1,0.2,0.5,1.0,1.2$), which also implies variations in the mean ($\overline{A_{V}}$). We thus construct six different PDF functions, namely $p_{1}-p_{6}$ corresponding to each different width, respectively. In the special case where $\sigma$ tends to zero, the PDF tends to a $\delta$-function and thus the aforementioned quantities become identical333In symmetrical distributions, such as a Gaussian one, the mean, median and mode are always equivalent..

The top left panel of Fig. 9 shows the shapes of the considered $A_{V}$-PDFs. For each $\sigma$, the abundances of Hi, H2, Cii, Ci and CO are calculated using uniform density PDRs corresponding to the mean, median and mode, as well as using the §3.2 method (referred to as ‘PDF’ in Fig. 9). All these results are shown in the remaining panels of Fig. 9. For the distributions with $m=3$ mag, values of $\sigma\lesssim 0.3$ show no appreciable difference in the abundances of species as calculated using a particular method.

As $\sigma$ increases, the difference between the derived abundances based on the above methods becomes significant. It is interesting to note that the abundances calculated with uniform density PDRs based on the values of the mean and the median always overestimate the abundances of H2 and CO, whereas they always underestimate the abundances of atomic species (Hi, Ci) and the abundance of Cii. Although the mode and ‘PDF’ approaches follow the same trend for Hi, H2 and Cii with increasing $\sigma$, there is a discrepancy when calculating the abundances of Ci and CO. This is because the column density of this species peaks only for a particular range of $A_{V}$ (see §2) and it is thus more strongly dependent on how this range of $A_{V}$ compares to the $A_{V}$ corresponding to the aforementioned approaches. Overall, the mode $A_{V}$ provides the best approach for single PDR calculations, but significant differences with the full PDF recommending using the here presented method that takes full account of the $A_{V}$-PDF.

6 Discussion

The method presented is quick and robust. Since it considers an $A_{V}$-PDF distribution as input, it is able to give a more realistic estimation of the average abundances of species compared to uniform (1D) PDR calculations with densities representing the average number density of a particular object. Thus, it can be used to provide fast and reasonable estimations of ISM structures spanning from a few parsec to kiloparsec scales. However, there are several assumptions and limitations made at this stage, which can be improved in future works.

In the most basic framework adopted here, it has been assumed that the chemical evolution time, $t_{\rm chem}$, is significantly smaller than the hydrodynamical evolution time, $t_{\rm hydro}$. The condition of $t_{\rm chem}\ll t_{\rm hydro}$ corresponds to quiescent (non-star-forming) molecular clouds the ISM of which is considered to be in statistical equilibrium. However, it is known that turbulence, which is always present in real clouds, can affect the chemical evolution of the ISM by cycling material from low to high density gas and back, and thereby altering the abundances of species. The H2 formation time is $t_{\rm H_{2}}\!\sim\!10^{9}\,{\rm yr}/(n_{\rm H}/{\rm cm}^{-3})$ (Hollenbach & McKee, 1979) which for a gas cloud of density $n_{\rm H}\sim 100\,{\rm cm}^{-3}$ gives $t_{\rm H_{2}}\sim 10^{7}\,{\rm yr}$. This time is needed for the atomic-to-molecular transition to reach equilibrium (Hollenbach & Tielens, 1999). On the other hand, for quiescent molecular clouds, the turbulent diffusion time (over which turbulence has significantly decayed) is a few $10^{6}\,{\rm yr}$ (Xie et al., 1995) and thus comparable to the H2 formation time for this specific density. Thus in a turbulent medium, equilibrium may be reached at different times, as discussed below.

Modern 3D numerical simulations of turbulent clouds that include time dependent chemistry find that turbulence can speed up the H2 formation time at a respective $n_{\rm H}$ when compared to a quiescent case (Glover & Mac Low, 2007a, b; Glover et al., 2010). Recent simulations of Milky Way-type molecular clouds by Seifried et al. (2017) indicate that such clouds may contain molecular hydrogen in regions with number densities of $n_{\rm H}\sim 30\,{\rm cm}^{-3}$ as a result of turbulent mixing, despite the fact that the local H2 formation time there is $t_{\rm H_{2}}\sim 10^{8}\,{\rm yr}$. This may have secondary effects in the connection of $A_{V}$-PDFs with the average abundances of species as presented here, since our PDR calculations do not account for such turbulent mixing and are thus not suitable to study turbulent media. For example, Bialy et al. (2017) find that once turbulence becomes supersonic, the development of strong density fluctuations distort the atomic-to-molecular transition. They further find that the mean and the median of the Hi column densities are affected by supersonic turbulence, but since turbulence decays more quickly as density increases, any changes to the $A_{V}{-}n_{\rm H}$ relation should predominantly affect low values of $A_{\rm V}$. Therefore, modifications of our methodology may be required only at the low $A_{V}$ end to account for turbulent mixing in regions of supersonic turbulence.

In all applications discussed in §4, the input parameters ($\chi$, $\zeta_{\rm CR}$, $Z$) are assumed to be uniform. This condition may hold only for small ISM patches and individual molecular clouds, unless there are strong, embedded UV sources. Despite this caveat, uniform boundary conditions are frequently adopted in 3D simulations of quiescent or star-forming clouds (e.g. Glover et al., 2010; Clark et al., 2013; Wu et al., 2017; Bisbas et al., 2017a, b, 2018). Likewise, the cosmic-ray ionization rate is considered constant with no attenuation as a function of $A_{V}$ (e.g. Clark et al., 2013; Bisbas et al., 2017a), and metallicity is assumed constant (e.g. Glover & Clark, 2012; Papadopoulos et al., 2018). However, when studying the ISM at larger (kpc) scale, the above assumptions may not hold as all three aforementioned input parameters may vary spatially. For example, Wolfire et al. (2003) provide an analytical model of how the FUV radiation field, cosmic ray ionization rate and metallicity vary with galactic radius in the Milky Way. Their model shows that all these quantities vary significantly (e.g. up to 20 times) in different places of the Galaxy and therefore should be accounted for when studying large ISM patches. Recent numerical simulations by Smith et al. (2014), Girichidis et al. (2016) and Li et al. (2018) also highlight the complexity of the ISM in this regard. Thus it can be expected that the ‘environmental’ parameters considered in our framework, need not always be assumed as constant values, but rather vary according to some distribution as well, which should be combined with the respective $A_{V}$-PDF. In this situation, Eqn. (6) should be modified as follows:

[TABLE]

where PDFj is the probability distribution function of the ‘environmental’ parameters ($\chi$, $\zeta_{\rm CR}$, $Z$) and $\ell$ is the number of different PDR suites, each suite consisting of $q$ simulations for a particular input parameter.

An important advancement in the present method will be the inclusion of radiative transfer to calculate the line emission of key tracers such as [CII], [CI], different CO transitions and others that are captured by observational studies. Several algorithms are now moving into this direction, as it provides the most direct way of comparing observational data with numeric simulations. Such tools include radex (van der Tak et al., 2007), lime (Brinch & Hogerheijde, 2010), radmc-3d (Dullemond et al., 2012) and more recently PyRaTe (Tritsis et al., 2018). However, calculation of abundances using PDR chemistry resulting in more accurate level populations and therefore optical depths, is not explicitly done in the aforementioned codes. A future aim of this work is to couple the presented method with a radiative transfer scheme. This will allow the determination of average line intensities for complex ISM distributions and predict how fundamental quantities—such as CO or Ci conversion factors—change as function of various ‘environmental’ parameters (e.g. Papadopoulos et al., 2004; Bolatto et al., 2013). For example, the middle panels of Fig. (7) show that as $\zeta_{\rm CR}$ increases, the carbon phase is more sensitive to the increment of $\zeta_{\rm CR}$ than the atomic-to-molecular transition is. As already highlighted by Bisbas et al. (2015a, 2017a), this is a strong indicator that the common technique of using CO as a tracer of H2 might be biased. In particular in high-redshift galaxies, that appear to have elevated $\zeta_{\rm CR}$ values, additional and more detailed examination is needed to establish the corresponding values of CO- and Ci-to-H2 conversion factors. At the same time, observations of low metallicity environments, where CO emission is traditionally found to be very weak (e.g. Leroy et al., 2005; Schruba et al., 2012, 2017), would greatly benefit from detailed knowledge of how (non-) detected line intensities relate to ISM properties.

7 Conclusions

We present a new numerical algorithm that uses $A_{V}$-PDF and $A_{V}{-}n_{\rm H}$ relations as inputs to estimate the average abundances of species with look-up tables from a suite of pre-calculated 1D thermo-chemical photodissociation region (PDR) calculations. This approach is much faster to full hydrochemical calculations where the computational cost is very high. It is thus a quick and alternative tool to estimate the average conditions of large ISM scales and assess trends with parameters, suitable for extragalactic studies.

Two hypothetical, lognormal $A_{V}$-PDF distributions are examined: a first one corresponding to a low density, atomic-dominated medium (with $\overline{A_{V}}=0.75$ mag and width $\sigma=0.5$) and a second one corresponding to a dense molecular cloud (with $\overline{A_{V}}=4$ mag and $\sigma=0.8$). For each $A_{V}$-PDF, the impact of different FUV radiation fields, cosmic-ray ionization rates and metallicities on the abundance of species was investigated. A third distribution (with $\overline{A_{V}}=1.8$ mag and $\sigma=0.49$) was considered describing the observed $A_{V}$-PDF of the Taurus molecular cloud.

In the case of the low density medium, the gas remains fully Hi- and Cii-dominated. However, in the case of the denser molecular cloud, we find that, although the cloud remains entirely molecular, its carbon phase can either be CO, Ci or even Cii-dominated, depending on the radiation field, cosmic ray rates and metallicity. For the particular case of the Taurus-like cloud and when considering the log-normal component only, we find that its molecular phase is equally Ci and Cii for a wide range of cosmic-ray ionization rates at a range of low UV intensities. It is Cii-rich for UV intensities up to a few for low cosmic-ray ionization rates, after which the medium becomes atomic. When considering the power-law tail in addition to the log-normal component, we find a small increase in H2 but a more striking increase in the CO abundances. The here presented algorithm is flexible and computationally efficient, and thus can be a powerful tool to determine the chemical composition of large ISM regions without the necessity to perform detailed and expensive astrochemical calculations.

We note that the algorithm may need to be modified to model kiloparsec-scale ISM regions. One reason for this is that supersonic turbulence can dynamically mix the gas so that a non-equilibrium approach is required. In addition, the input parameters of the external FUV radiation field, cosmic-ray ionization rate and metallicity are considered to be uniform for now. While this assumption seems reasonable for small ISM patches, it becomes challenged when considering the ISM at large scales. In principle, large scale variations can be accounted for by assuming distribution functions of the ISM properties and external parameters, and applying our algorithm ensembles of smaller ISM patches. Finally, coupling the present algorithm with radiative transfer will provide a powerful tool for studying the average intensities of emission lines of large ISM regions with complex substructure. This may be of particular importance for studies of high-redshift galaxies that often remain marginally resolved by observations and one-zone PDR calculations may fail to deliver reliably characterizations of their physical properties.

Acknowledgements

The authors thank the referee for the useful comments which improved the clarity of this work. This work is supported by a Royal Netherlands Academy of Arts and Sciences (KNAW) professor prize, and by the Netherlands Research School for Astronomy (NOVA). TGB acknowledges funding by the German Science Foundation (DFG) via the Collaborative Research Center SFB 956 “Conditions and impact of star formation”. The authors thank Alyssa Goodman, Desika Narayanan, Blakesley Burkhart, Shmuel Bialy, Amiel Sternberg, Daniel Seifried, Sebastian Haid, Annika Franeck, Stefanie Walch, Nicola Schneider, Volker Ossenkopf and Markus Röllig for insightful discussions. This research has made use of NASA’s Astrophysics Data System.

Appendix A The Hi-to-H2 transition

The network of chemical reactions following the inclusion of the suprathermal formation of CH+, described in 3.1, destroys the H2 molecule at low column densities, thereby shifting the atomic-to-molecular transition to higher $A_{V}$. To investigate this further, an isothermal PDR simulation of a uniform density cloud is performed ($T_{\rm gas}=50\,{\rm K}$, $n_{\rm H}=10^{3}\,{\rm cm}^{-3}$, $\chi/\chi_{0}=1$, $\zeta_{\rm CR}=10^{-16}\,{\rm s}^{-1}$) where the ‘suprathermal’ routine is switched ON and OFF.

Figure 10 shows the results of the fractional abundances (top) and the corresponding destruction rates (bottom) of the H2 molecule as a function of $A_{V}$. Once the suprathermal routine is switched ON, H2 is destroyed by Reaction 1 (see discussion in §3.1) which dominates over the H2 photodissociation reaction, therefore delaying the formation of H2-rich gas.

Appendix B The $A_{V}{-}n_{\rm H}$ relation

A collection of four different $A_{V}{-}n_{\rm H}$ relations has been considered in this work. In particular, the relations by Glover et al. (2010); Van Loo et al. (2013); Safranek-Shrader et al. (2017) and Seifried et al. (2017) have been taken into account. Glover et al. (2010) (blue solid line) performed magnetohydrodynamical (MHD) models of a periodic box of side length $L=20\,{\rm pc}$ filled with a $\langle n_{\rm H}\rangle=300\,{\rm cm}^{-3}$ gas and used a six-orthogonal ray approximation to estimate the effective visual extinction at all computational elements in the box. A six-orthogonal ray approximation was also used by Van Loo et al. (2013) (yellow solid line) to obtain the $A_{V}{-}n_{\rm H}$ relation over a galactic disk with $20\,{\rm kpc}$ in diameter and a spatial resolution of $\sim 8\,{\rm pc}$ or less (down to $\sim 0.5\,{\rm pc}$ over a smaller, $1\,{\rm kpc}$, region to study the GMC properties). Safranek-Shrader et al. (2017) (green solid line) performed MHD simulations of a $512\times 512\times 1024\,{\rm pc}^{3}$ box with a uniform resolution of $2\,{\rm pc}$ and they estimated the effective $A_{V}$ by considering the fraction between the unattenuated radiation field and the attenuated one due to dust. Finally, SILCC-Zoom simulations by Seifried et al. (2017) (red solid line) used a HEALPix based (Górski et al., 2005) 48-ray approximation to calculate the effective $A_{V}$ of an $88\times 78\times 71\,{\rm pc}^{3}$ region at a $0.12\,{\rm pc}$ resolution. As can be seen in Fig. 11, all aforementioned works find a very similar relation for $A_{V}{-}n_{\rm H}$ for ISM regions exceeding a few tens of parsecs in size up to a whole galaxy.

The characteristic of all aforementioned $A_{V}{-}n_{\rm H}$ relations is that they represent the most probable $A_{V}$ for a given $n_{\rm H}$ and thus they should not be considered as an one-to-one relation. As can be seen from detailed hydrodynamical models, for each $n_{\rm H}$ an $A_{V}$-PDF is associated with and needs be taken into account for a more accurate estimation of the average ISM property of that particular number density. In the case that this $A_{V}$-PDF at constant $n_{\rm H}$ is known (denoted as $\left[~{}\right]_{n_{\rm H}}$), the following equation (see also Eqn. 7) can be applied

[TABLE]

where $Q$ is the resolution (discretization) of the $\left[A_{V}-{\rm PDF}\right]_{n_{\rm H}}$. In this work, the above set of equations are not applied. Instead, a bounded average $A_{V}{-}n_{\rm H}$ is considered as explained below.

The solid black line represents the $A_{V}{-}n_{\rm H}$ relation that is assumed as the default in the present work. This is the average function when all aforementioned relations are taken into account. Since the spread in $A_{V}$ for constant $n_{\rm H}$ is not known, a lower and an upper value of $A_{V}{-}n_{\rm H}$ is assumed as derived using the equation

[TABLE]

where $g(A_{V},n{\rm H})$ is the ‘average’ $A_{V}{-}n_{\rm H}$ function (solid black line of Fig. 11), $\lambda=1$ is a fudge factor defining the upper bound (dot-dashed line), $\lambda=-1/2$ the lower bound (dashed line) and $n_{\rm min}$, $n_{\rm max}$ are the minimum and maximum values respectively of the H-nucleus number density used in $g(A_{V},n{\rm H})$. These bounds are decreasing with increasing $n_{\rm H}$ to approximate the trend of $A_{V}{-}n_{\rm H}$ spread in accordance to the above hydrodynamical models.