The role of molecular gas formation on chemical evolution

# Galaxy chemical evolution models: The role of molecular gas formation

Mercedes  Mollá, Ángeles I. Díaz, Yago Ascasibar, and Brad K. Gibson
Departamento de Investigación Básica, CIEMAT, 28040, Madrid. Spain
E.A. Milne Centre for Astrophysics, University of Hull, Hull, HU6 7RX, United Kingdom
E-mail:mercedes.molla@ciemat.es
Accepted Received ; in original form
###### Abstract

In our grid of multiphase chemical evolution models (Mollá & Díaz, 2005), star formation in the disk occurs in two steps: first, molecular gas forms, and then stars are created by cloud-cloud collisions or interactions of massive stars with the surrounding molecular clouds. The formation of both molecular clouds and stars are treated through the use of free parameters we refer to as efficiencies. In this work we modify the formation of molecular clouds based on several new prescriptions existing in the literature, and we compare the results obtained for a chemical evolution model of the Milky Way Galaxy regarding the evolution of the Solar region, the radial structure of the Galactic disk, and the ratio between the diffuse and molecular components, HI/H. Our results show that the six prescriptions we have tested reproduce fairly consistent most of the observed trends, differing mostly in their predictions for the (poorly-constrained) outskirts of the Milky Way and the evolution in time of its radial structure. Among them, the model proposed by Ascasibar et al. (2017), where the conversion of diffuse gas into molecular clouds depends on the local stellar and gas densities as well as on the gas metallicity, seems to provide the best overall match to the observed data.

###### keywords:
Galaxy: abundances–Galaxy: molecular gas–Galaxy: star formation
pagerange: Galaxy chemical evolution models: The role of molecular gas formationLABEL:lastpagepubyear: 2011

## 1 Introduction

Chemical elements appear in the Universe as a consequence of three main processes of production: the Big-Bang nucleosynthesis, the spallation process and the stellar nucleosynthesis. The first one only explains H and He and some traces of other light elements, while the second one transforms, by the collision of cosmic rays, some chemical elements in other lighter by the fragmentation of the weightier ones. The third one is actually the most important in forming the majority of elements of the Universe.

Chemical evolution models are the classical tool to interpret observed elemental abundances, and associated quantities such as gas and stellar surface densities, star formation histories, and the distribution of stellar ages. Elemental patterns carry the fingerprint of star formation timescales from their birth location, regardless of a star’s present-day position. Chemical evolution codes solve a system of first order integro-differential equations, assuming an analytical star formation (SF) law, an initial mass function (IMF), stellar lifetimes, and nucleosynthetic yields.

In Mollá & Díaz (2005, hereafter MD05) we calculated a grid of 440 theoretical galaxy models, applied to 44 different radial distributions of galaxy mass and with 10 different values between 0 and 1 for the molecular gas and star formation efficiencies, which implies 10 different evolutionary histories for each one of the different sizes/masses of galaxies. This grid was calibrated on the Milky Way Galaxy (MWG), by using the data of surface densities of HI, H, stars and SF rate (SFR), and also of C, N, and O abundances, and with time evolution data for the solar region such as the metal-enrichment or the SFR histories. Our model has the advantage over similar ones in the literature of assuming that stars form in two steps: first the molecular clouds form from the diffuse gas following a classical Schmidt law, then cloud-cloud collisions lead to the formation of stars. Therefore, we derive radial distributions for both gas phases in a galaxy. This allows a certain delay in the star formation, doing smoother star formation histories. Results from that work show that, as expected, the SFR radial distributions follow those of the molecular gas; however both showing a drop in the inner regions of the disks, at variance with observations.

In order to explore those differences further, we have started the computation of a new grid of models including updates to all the relevant ingredients as the stellar yields in Mollá et al. (2015), or the infall rate law in Mollá et al. (2016). In the present work, our aim is to improve the predicted and SFR radial profiles, while maintaining abundance radial gradients in agreement with those observed. For that we will study the changes in the molecular gas radial distributions when using different prescriptions existing in the literature about the formation of molecular clouds from the diffuse gas and select which of them is the best one to reproduce the existing observations on the MWG taken as the calibration galaxy.

Section 2 of this work presents some of these new prescriptions taken from Blitz & Rosolowsky (2006); Krumholz et al. (2008, 2009); Gnedin & Kravtsov (2011); Ascasibar et al. (2017), describing the different models we compute. In Section 3 we show our model results for the MWG, comparing them with recent observational data as given in Mollá et al. (2015, Appendix). In Section 4 we check if our results reproduce the universal radial profile H/HI given by Bigiel et al. (2008), analyzing the dependencies of this ratio on the galaxy morphological type and on the dynamical mass. Finally, in Section 5 we draw our Conclusions.

## 2 Chemical evolution models: HI to H2 conversion prescriptions

In this section we are going to summarize our basic multiphase chemical evolution model (hereinafter MULCHEM) and describe the different prescriptions used in this work, to be included in it, to form molecular clouds, H, from the diffuse gas, HI.

### 2.1 Description of the basic model

As described in Mollá et al. (2016), we start computing the radial mass distributions for 16 theoretical galaxies following Salucci et al. (2007) equations, defined in terms of the virial mass of the dark matter haloes, , and their associated rotation curves for the galaxy disk and halo. The virial masses are in the range  M, with associated barionic disk masses in the range  M. The initial gas in each model collapses onto the disk on timescales based in the Shankar et al. (2006) relationship, which gives the ratio between the galaxy disk and the dark matter halo masses, . From the rotation curves, we calculate the radial distribution of the galaxy dynamical mass, as well as the predicted disk mass, , at the present time. Thus, we obtain the infall rate necessary to produce, at the end of a given model evolution, the appropriate disk for each galaxy dynamical mass. These infall rates are defined by a collapse time scale given by these two accumulated masses distributions, an through the corresponding mass in each radial region, which we define as 1 kpc wide:

 ΔM\scriptsize tot(R)=M\scriptsize tot(
 ΔM\scriptsize D(R)=M\scriptsize disk(
 τ(R)=−13.2ln(1−ΔM\scriptsize D(R)ΔM\scriptsize tot(R))[Gyr] (3)

The infall rates inferred from expression (3) vary among different galaxies and among the different radial regions within them, and they imply a modest evolution with time of the infall rate for the disks but a strong evolution for bulges. The radial regions in a disk for a galaxy with  M(i.e., a MWG-like analog) show present day infall rates of  M yr for galactocentric radii  kpc, in agreement with Sancisi et al. (2008)’s data, being much lower for the outer regions,  kpc, (see Mollá et al., 2016, for a detailed discussion).

Once defined our initial scenario with these radial distributions of masses, it is necessary to solve the equation system which solves the evolution along the time of each radial region located at a galactocentric distance R, in the halo and in the disk. In MULCHEM we solve the following system:

 dgH(R)dt = −(κh,1(R)+κh,2(R))gnH(R)−f(R)gH(R)+WH(R) (4) ds1,H(R)dt = κh,1(R)gnH(R)−D1,H(R) (5) ds2,H(R)dt = κh,2gnH(R)−D2,H(R) (6) dgD(R)dt = −κc(R)gnD(R)+κ′a(R)c(R)s2,D(R)+κ′s(R)c2(R) (7) +f(R)gH(R)+WD(R) dc(R)dt = κc(R)gnD(R)−(κa,1(R)+κa,2(R)+κ′a(R))c(R)s2,D(R) (9) −(κs,1(R)+κs,2(R)+κ′s(R))c2(R) ds1,D(R)dt = κs,1(R)c2(R)+κa,1(R)c(R)s2,D(R)−D1,D(R) (10) ds2,D(R)dt = κs,2(R)c2(R)+κa,2(R)c(R)s2,D(R)−D2,D(R) (11) drH(R)dt = D1,H(R)+D2,H(R)−WH(R) (12) drD(R)dt = D1,D(R)+D2,D(R)−WD(R) (13)

These equations predict the time evolution of the different phases of the model. More importantly, we consider in our model MULCHEM (Ferrini et al., 1992, 1994; Mollá & Díaz, 2005) two phases of gas in the disk: diffuse gas, , and molecular gas, , which are well known to be essential ingredients in the process of star formation, but have only recently been explicitly included in other chemical evolution models (e.g. Kubryk, Prantzos, & Athanassoula, 2015b, a). MULCHEM divides stars in two ranges, , and , denoting low-mass and intermediate, and massive stars, respectively, and stellar remnants, . In all cases subscripts and correspond to disk and halo, respectively. The infall rate, , is the inverse of the collapse time , given by the expression 3, and the death rates are computed with:

 D1,H,D(R,t) = ∫m∗mminΨH,D(R,t−τm)mϕ(m)dm (14) D2,H,D(R,t) = ∫mmaxm∗ΨH,D(R,t−τm)mϕ(m)dm, (15)

being the mass limiting the low mass stars from intermediate mass and massive stars, and the lower and upper mass limits of the initial mass function , and the main sequence lifetime of a star of mass .

As explained in MD05, in our model it is assumed that the star formation takes place following a Schmidt law in the halo regions; in the disk, however, it occurs in two steps: first molecular clouds, , form from diffuse gas, , then stars form through cloud-cloud collisions. Besides, a second star formation process appears resulting from the interaction of massive stars, with the molecular clouds surrounding them. Therefore, we have different processes defined in the galaxy:

1. Star formation by spontaneous fragmentation of gas in the halo: , where we use

2. Disk formation by gas accretion from the halo or protogalaxy:

3. Clouds formation by diffuse gas. In our standard models this process is . (with ). This cloud formation law will be modified in the present models as we will describe in next subsections

4. Star formation due to cloud collision:

5. Diffuse gas restitution due to cloud-cloud collision:

6. Induced star formation due to the interaction between clouds and massive stars:

7. Diffuse gas restitution due to the induced star formation 111 Massive stars induce SF in the surrounding molecular clouds but their radiation also destroys a proportion of them.: ,

where , , and are the proportionality factors of the star formation in the halo, the cloud formation, the cloud-cloud collision and the cloud-massive stars interactions (last two create stars from molecular clouds)222For the sake of avoiding misinterpretations with other quantities, we have changed the letters , , and , used in our previous works to denote the parameters defining the star formation in the halo, the molecular clouds formation in the disk, the cloud-cloud collisions and the cloud-massive stars interactions, by , , , and , respectively. Thus, the old efficiencies , , and are called now , , , and .. Since stars are divided in two groups, , and , the parameters involving star formation are divided in the two groups too, thus: , , and , where terms and refer to the restitution of diffuse gas due to the cloud-cloud collisions and massive stars-cloud interaction processes.

Thus, the star formation law in halo and disk is:

 ΨH(R,t) = (κh,1(R)+κh,2(R))gH(R)n (16) ΨD(R,t) = (κs,1(R)+κs,2(R))c(R)2+(κa,1+κa,2)c(R)s2,D(R) (17)

The factors , , and have a radial dependence, as discussed in previous studies (Ferrini et al., 1994; Mollá, 2014):

 κh(R) = ϵh(G/VH(R))1/2 (18) κs(R) = ϵs(3./VD(R)) (19) κa(R) = ϵa(Gρc)1/2/, (20)

where is the universal gravitational constant, and are the halo and the disc volumes of each radial region, is the average cloud density and is the average mass of massive stars. Thus, we converted these parameters in relationships with the volume and with other quantities that we called efficiencies, , and , which represent probabilities (in the range ) associated to the processes of conversion among the different phases and that we assume constant for all radial regions within a galaxy. The efficiency to form stars in the halo, , is obtained through the selection of the best value , able to reproduce the SFR and abundances of the Galactic halo (see Ferrini et al., 1994, for details). We assumed that it is approximately constant for all halos. The last efficiency, , was also obtained from the best value for the MWG, and it is also assumed constant for all galaxies, since these interactions between clouds and massive stars are local processes. Thus, we have only one free parameter: . In the present work we are going to modify the formation of H, keeping this efficiency as a free parameter with values in the range [0–1] as in MD05 and Mollá (2014). We compute 10 possible values using:

 ϵs=exp−NT2/8, (21)

where NT is a free parameter (although related to the morphological type, see details in Sections 2.2 and 2.3), with 10 assigned values between 1 and 10, as given in column 2 of Table 1.

The equations of the chemical abundances are:

 Xi,H(R)dt=Wi,H(R)−Xi,H(R)WH(R)gH(R) (22)
 Xi,D(R)dt=Wi,D(R)−Xi,D(R)WD(R)+f(R)gH(R)[Xi,H(R)−Xi,D(R)]gD(R)+c(R) (23)

where are the mass fractions of the 15 elements considered by the model: H, D, He, He, C, O, N, C, Ne, Mg, Si, S, Ca, Fe, and the rich neutron isotopes created from C, O, N and from C; and the restitution rates are:

 Wi,H,D(R,t)=∫mmaxmmin[∑j~Qij(m)Xj(t−τm)]ΨH,D(t−τm)dm (24)

To compute the elemental abundances we use the technique based in the matrices formalism (Talbot & Arnett, 1973; Ferrini et al., 1992; Portinari, Chiosi, & Bressan, 1998). Following this formalism, each element of a matrix, gives the proportion of a star which was element and is ejected as when the star dies. Thus,

 Qi,j(m)=mi,j,expmj, (25)

and

 Qi,j(m)Xj=mi,j,expm. (26)

If we take into account the number of stars of each mass , given by the initial mass function, that eject this mass , we have:

 ~Qi,j(m)=Qi,jmϕ(m), (27)

and, therefore, the term in Eq. (22) represents the mass of an element ejected by all stars of mass . This method allows to relax the hypothesis of solar proportions in the ejection, since each element relates with its own sources. It was originally introduced to compensate for the lack of stellar models of different metallicities, when the dependence on Z was not included in the stellar yields calculations Now that stellar yields for different metallicities are available, the use of matrices allows to us to take into account possible differences of chemical composition within a given Z. Stellar yields have usually been computed assuming only solar relative abundances among the different elements at a given Z. However, the relative abundances of elements are not always solar nor constant along the evolutionary time. Therefore, as Portinari, Chiosi, & Bressan (1998) explained, "The matrix links any ejected species to all its different nucleosynthetic sources, allowing the model to scale the ejecta with respect to the detailed initial composition of the star through the Xj’s".

Following Mollá et al. (2015), we use the stellar yield sets from Limongi & Chieffi (2003); Chieffi & Limongi (2004) for massive stars, together with yields from Gavilán et al. (2005, 2006) for low and intermediate mass stars, combined with the IMF from Kroupa (2001), which is one of the best combinations able to reproduce the MWG data. For supernovae type Ia (SNe-Ia) we use the rates given by Ruiz-Lapuente et al. (2000), who provides a numerical table (private communication) with the time evolution of the supernova rates for a single stellar population, computed under updated assumptions about different scenarios and probabilities of occurrence. The stellar yields for SNe-Ia are those given by Iwamoto et al. (1999). As shown in Mollá et al. (2015), the observational data for the MWG are well reproduced with our MULCHEM model. The set of data used for comparison with the results of our models were included as Appendix A in Mollá et al. (2015) and refers to a) the time evolution of the Solar Region –located at a galactocentric distance of 8 kpc–, (star formation and enrichment histories, as and , alpha-element over-abundances , as a function of the metallicity , and the metallicity distribution function, MDF, or proportion of stars in each [Fe/H] bin); b) the radial distributions within the Galactic disk for the present time ( and surface densities, and , the stellar profile , and the radial distribution of the star formation rate surface density, ); and c) the radial gradients of C, N and O elemental abundances.

### 2.2 Our standard model

The molecular cloud phase in our standard model follows the equation:

 dc(R)dt=κc(R)gD(R)n−κs(R)c2(R)−κac(R)s2,D(R), (28)

where the first term refers to the formation of molecular cloud from diffuse gas, while the next two others describes how these clouds disappear due to the cloud-cloud collisions (which create stars and also restitute diffuse gas to the ISM) and as the consequence of the interaction of massive stars with the cloud surrounding them (this process also create stars and new diffuse gas which return to the ISM).

The proportionality factor depend on the geometry of the regions, in a similar way that the explained for the three others in Eq. 18, 20, and 21:

 κc(R)=ϵc(G/VD(R))1/2 (29)

In MD05 we used the data from Young et al. (1996) corresponding to and – as calculated from the masses of HI and H–, and the SFR, , to obtain a relation between both efficiencies, and , as a function of the galaxy morphological type (see Fig. 4 from MD05), finding that . Using Eq. 21 for , it implies:

 ϵc=exp−NT2/20 (30)

We computed this way both efficiencies simultaneously as a function of the free parameter , that may be associated to the morphological type index (see MD05 for details about this relationship). The resulting efficiencies are given in column 3 of Table 1. The model using these prescriptions is called STD.

### 2.3 Our modified model

Looking at Fig. 4 from MD05, it can be actually seen that, although a constant value for the ratio is statistically significant, a straight line fit is also reasonable:

 lnϵclnϵs=0.12+0.07NT (31)

As for , the efficiency to form stars from molecular clouds, we have kept their values as in MD05, we now modify accordingly:

 ϵc=exp−NT28(0.12+0.07NT) (32)

This changes slightly the values for both efficiencies as shown in Table 1 where, for each value of (column 1), we give the efficiency to form stars, (column 2), and the old and new values for the molecular gas formation efficiency, (columns 3 and 4, respectively).

In Fig. 1 we plot as a function of both for the STD model (Eq.  30) and our modified model (Eq.  32) that we call MOD. It is evident that according to the newly adopted function, late type galaxies will form a smaller amount of molecular gas, while this phase will be increased for the earlier types, as compared to the STD model. The equation defining the molecular cloud mass, as the corresponding parameter are the same for both models, STD and MOD. Only the efficiencies change among them.

### 2.4 Ascasibar et al. (2017) prescriptions

In this case, we follow a prescription based on Ascasibar et al. (2017) to calculate the value of the parameter . The most relevant feature of this model is that the time scale for the formation of molecular clouds:

 κc,ASC(R)=2ndust(R)⟨σv⟩dust, (33)

which depends on the local number density of dust grains and the thermally-averaged cross-section for the condensation of hydrogen molecules onto their surface. We assume that the number of dust particles is proportional to the metal content of the gas, i.e. , and we set the normalization  cm s to the reaction rate estimated by Draine & Bertoldi (1996) for a gas of solar composition (, Asplund et al., 2009) at a temperature  K. This way:

 κc,ASC(R)∼2nH(R)⟨σv⟩⊙Z(R)+ZeffZ⊙ (34)

Other physical processes unrelated to the dust grains (e.g. the channel) also contribute to the formation of hydrogen molecules. Although such alternative mechanisms are only relevant for the lowest metallicities, they are extremely important in the early universe, and we represent their combined action by an effective term that becomes dominant when (see e.g. Glover & Jappsen, 2007).

The gas density includes the contribution of both diffuse gas and molecular clouds, and, using the classical gas equation, we set each one as:

 n{g,c}(R)∼P(R)kB(T{g,c}+Teff) (35)

where denotes the Boltzmann’s constant, and we adopt  K. In addition to thermal pressure, non-thermal pressure support from turbulent motions, cosmic rays, and magnetic fields is accounted for by an extra term with  K (implying, under the assumption of equipartition, a velocity dispersion of the order of  km/s for the diffuse gas and  km/s for the molecular clouds). is the mid-plane pressure which depends, in turn, on the total and gas surface densities (see e.g. Elmegreen, 1989; Leroy et al., 2008):

 P(R)≈π2GΣgas(R)[Σgas(R)+Σ⋆(R)] (36)

Substituting in equation (34), one finally arrives to:

 κc,ASC(R)=2.67Σgas(R)[Σgas(R)+Σ⋆(R)][Z(R)+1.4×10−5] (37)

with surface densities in units. The evolution of the molecular mass is given by:

 dc(R)dt=κc,ASC(R)gD(R)−κs(R)c2(r)−κac(R)s2,D(R), (38)

in a similar way to models STD and MOD with , (of course, it is also necessary to modify this index in the term associated to the formation of molecular clouds in the equation for the diffuse gas). We will refer to this prescription as ASC.

The parameter is represented in Fig. 2 as a function of the total gas (atomic molecular) surface density for different stellar densities, as labelled. As this ratio also depends on the metallicity, we have used the local relation between metallicity and stellar-to-gas fraction proposed by Ascasibar et al. (2015), to compute the parameter drawn in that Fig.2:

 Z(S)Zmax=S1+S/Scrit (39)

where denotes the surface density ratio between stars and total gas in each radial region, is the asymptotic metallicity attained in the chemically-evolved limit , and marks the transition between the chemically-young regime, where , and the flat asymptotic behaviour (see Ascasibar et al., 2015, for more details). This equation (39) is, however, used here only for purposes of doing the figure. In the ASC model the metalllicity is obtained as a result in each radial region and time step. In practice, is roughly proportional to when and (i.e. in the very first stages of galactic evolution) and to the product when (i.e. over most of the life of the galaxy).

### 2.5 Gnedin and Kravtsov (2011) prescriptions

In this case, as in the next two cases, we are going to use the prescriptions given by some authors which calculate the relation between the molecular gas and the total gas, defined through the ratio:

 fH2=ΣH2Σgas, (40)

which is equivalent to in our notation (since the area is the same for both surface densities) . This implies that:

 c(R)=fH21−fH2gD(R), (41)

and consequently,

 dc(R)dt=fH21−fH2dgD(R)dt+gD(R)ddt(fH21−fH2) (42)

Gnedin & Kravtsov (2011, hereinafter GNE) presented a detailed description of a phenomenological H formation model. Their results give the atomic-to-molecular transition by means of the ratio between the molecular and the total hydrogen, which they give as a function of gas density or column density with dependencies on the dust-to-gas ratio, and on the radiation flux at 1000 Å, . They parametrize these dependencies in convenient fitting formula useful for their inclusion in semi-analytic models and cosmological simulations that do not include radiative transfer and H formation.

Values for range from 10 to 3 relative to the MWG value, and scaling to metallicity, while those for , are between 0.1 and 100 times the MWG value. The ratio is approximated by the equation:

 fH2,GNE=11+exp(−4xGNE−3x3GNE), (43)

where

 xGNE = Λ3/7ln(DMWGnHΛn∗), (44) n⋆ = 25cm−3, (45) Λ = ln(1+gGNED3/7MWG(UMWG/15)4/7) (46)

The function is a function giving the transition between self-shielding and dust shielding, defined as:

 gGNE = 1+αGNEsGNE+s2GNE1+sGNE, (47) with: (48) sGNE = 0.04D∗+DMWG, (49) αGNE = 5UMWG/21+(UMWG/2)2 (50) D∗ = 1.5×10−3×ln(1+(3UMWG)1.7) (51)

In these equations refers to the hydrogen particle total density in their different phases and denotes a fiducial value included only for normalization purposes in (32).

We have computed these functions for 4 values of ,1,10 and 100 and 8 values of , which is given by its relation with the metallicity: , from to . The results are shown in Fig. 3 for a value of in the top panel, only modifying , and for a value of for variable values of in the bottom one. With these equations, the ratio depends mainly on the gas surface density with the two other parameters modifying the final results: changing slightly the slope, and shifting the curve to the right or the left. In fact, in STD and MOD models this dependence on the gas surface density also appears, through the disk volume in Eq. 29, which finally translates into a radial dependence, i.e. the diffuse to molecular cloud conversion depends on the gas density in each radial region. In the Gnedin & Kravtsov (2011) case, it also depends on the metallicity (by means of ) and the parameter .

Such as we may see in that Fig. 3, the function may be approximated as a function with the adequate transformation of the x-axis, using (where is the ). For instance, for and , we find:

 fH2=12tanh[3×(x−1.7)]+1 (52)

We plot this fitted function in Fig. 3b) as a red short-long-dashed line.

By assuming this dependence we may calculate:

 dfH2dt = dtanh(u)dududt (53) = 12sech2(u)A(dlogΣgD+cdt) (54) = A2sech2u1gD+c(dgDdt+dcdt) (55)

Using this last expression and Eq.52 within Eq.42, and with the adequate mathematics, we obtain this equation for the formation of molecular clouds from diffuse gas:

 dcfordt=κcdgDdt, (56)

where:

 κc = fH2+(A/2)sech2(u)1−fH2−(A/2)sech2(u) (57) = 1+tanhu+Asech2(u)1−tanhu−Asech2(u) (58)

The function tends to zero for lower and higher gas densities (see Fig.3), that is, the derivative of the function when this is almost flat, is practically zero, what it occurs at the beginning and at the end of each curve. In this case . In the middle of the function, when this increases rapidly from 0 to 1, the function has a significant value which may reach a value of 1. We have estimated the value of for the case and for the region defined by , or . It results that the expression began to be negative when or , that is, . We could interpret that the molecular clouds start to be destroyed when . Since we are interested only in the formation processes, we prefer to assume that the formation of molecular clouds is proportional to , and computing this destruction with the usual equations used in our code. This way, our final equation for the evolution of both phases of gas in the disk will be given by:

 dc(R)dt = κc(R)dgD(R)dt−κs(R)c2(r)−κac(R)s2,D(R) (59) dgD(R)dt = κ′s(R)c2(r)+κ′ac(R)s2,D(R)+f(R)gH(R)+WD1+κc(R) (60)

with

 κc,GNE(R)=fH21−fH2, (61)

calculated with the expressions given by Eq. 43 and the other terms defined as in the STD and MOD models.

In order to include the above equations into our chemical evolution code, we use, as said, the metallicity to calculate the parameter at each time, and the total gas density to compute the hydrogen density . To compute the we need to use the density . We have assumed that this density may be calculated with the surface density of the total gas dividing by the width of the disk and doing the adequate change of units from to , we obtain that . Regarding the UV flux, we have assumed for all our simulations since we are modeling MWG. Moreover, it is difficult to evaluate its value in our chemical evolution models prior to the computation of the photometric evolution of each theoretical galaxy (only masses in the different phases are calculated in this step). In panel a) of Fig.4 we represent the as a function of the surface gas density for different values of the metallicity which, in our code, as also does the gas density, varies continuously when the star formation takes place. We label models following the prescriptions described in this section, GNE.

### 2.6 Blitz & Rosolowsky (2006) prescriptions

Blitz & Rosolowsky (2006) proposed that the ratio of atomic to molecular gas depends basically on the total pressure at each radius by the equation:

 Rmol(R)=MH2(R)MHI(R)=[P(R)P0]αP, (62)

being a parameter.

Using the expression:

 P(R)=π2GΣgas(R)[Σgas(R)+fσΣ⋆(R)], (63)

where is the gravitational constant, is the galactocentric radius and , is the ratio of the vertical velocity dispersions of gas and stars. Assuming an exponential function for the disk density, and a mean value for along the disk, they arrived to the final equation:

 Rmol(R)=1.3810−3[Σgas(R)+0.1Σgas(R)√(Σ⋆(R)Σ⋆,0]0.92 (64)

The equations giving the time evolution of the diffuse and molecular phases of gas are in this case the same as used in the GNE model. Since in our notation, we have a similar , calculated with Eq.64, while .

In this case the molecular ratio depends primarily on the gas surface density but also on the stellar surface density, which may be important when the disk galaxies are gas-poor. At any rate this equation refers to the total pressure averaged over a particular galactocentric distance, which may present problems for its application to cosmological simulations since no correction for clumpliness or the presence of a warm phase are taken into account, but it is totally valid for our purposes.

In panel b) of Fig.4 we represent for different values of the stellar surface density. Obviously in our code this density, as well as the gas surface density, varies continuously when the star formation takes place. Models calculated following this prescription are labelled BLI.

### 2.7 Krumhold et al. (2008, 2009) prescriptions

Krumholz et al. (2008, 2009) studied the transition of the diffuse into molecular gas finding that the conversion depends on the interstellar radiation field, through a dimensionless strength , here renamed , which includes the dependence on the dust properties (basically through the metallicity ), the radiation field and the atomic gas density as follows:

 fH2,KRU=1−34(sKRU1+0.25sKRU), (65)

where

 sKRU=ln(1+0.6χ+0.01χ2)0.04(Z/Z⊙)(ΣHIM⊙pc−2), (66)

where the function , in turn, is basically dependent on the metallicity:

 χKRU=3.11+3.1Z/Z0.365⊙4.1 (67)

These expressions are valid for and actually represent functions depending on metallicity and gas density. Models using these equations are called KRU.

It should be mentioned that in Krumholz (2013) the author performs a careful analysis of the formation of stars (and consequently of the molecular cloud formation) in the low molecular gas regime. According to his work, it is necessary to compute a minimum gas depletion time (in Gyr), which the author obtained from his Eqs. (21), (22) and (27), and then use this value to calculate . This technique produces different curves for very low , as appropriate for the early evolution of galaxies, and for the normal regime, when this parameter reaches values higher than 0.2–0.3 (see his Fig.1 where this behavior is shown). We have also included this method in our code to choose a different value of in the low molecular gas regime. As seen in our Fig. 4, it also makes the agreement with the GNE prescriptions improve. As in BLI model, the equations of the evolution for both phases of gas are similar to the ones given for GNE, with a similar value of in which is calculated as explained.

Panel c) of Fig.4 shows the evolution of with gas surface density for different values of the metallicity. This metallicity, as in the case of the GNE model, varies continuously in our chemical evolution code. However, the prescriptions of GNE and KRU result in rather different functions, as shown.

## 3 Results

### 3.1 Comparison with the Milky Way Galaxy

In this subsection we compare the results of MULCHEM applied to the MWG using the different prescriptions to form molecular clouds. The basic model is the same as described in Mollá et al. (2015) and in the above section 2. The scenario starts with a protogalaxy with a virial mass of M assumed spherically distributed following a Burkert (1995) profile. This mass, initially in gas phase, falls over the equatorial plane forming out the disc. The infall rate is assumed as the necessary to create, along a Hubble time, a disk with a radial profile as observed. Details about the evolution of this infall rate and comparison with existing data and with cosmological simulations results are given Mollá et al. (2016). Once this gas in the disk, molecular clouds are created from the diffuse gas, and then, by cloud-cloud collisions and by the interaction of massive stars with these clouds, stars form.

We run 6 models, summarized in Table 2. The only differences among the six MWG models are the different methods to convert HI in H and in the values of , which are computed from the corresponding efficiencies and volumes in each radial region for STD and MOD models, while are obtained as explained in the above section for the other four ones. In Table 2 we give for each model the main characteristics on which the conversion of HI to H depends and also the color of the line used in our figures for the MWG models and the references in which we base our work. In the first two models, the conversion of diffuse gas, , into molecular gas , occurs by assuming a Schmidt law with a power slope . In the last 4 models, however, the transition from diffuse to molecular phase is assumed as a linear process defined by . The equations given this ratio as a function of the total density of gas or as a function of other quantities are given in the previous section for the different authors.

In Fig. 5 we represent the time evolution of the Solar Region showing the SFR and the enrichment histories as a function of time (left panels) and the alpha-element abundance, [/Fe], and the MDF as a function of the Iron abundance, , taken as a proxy for time (right panels). In panels a), b) and c), it is rather evident that results for the final times, including the present one, are very similar for the six models. However, some differences appear at early times. Thus, all models show a maximum in the SFR, but for models GNE, BLI, and KRU this maximum is reached before the other three. In particular, ASC model has this maximum shifted towards slightly later times than the rest of models ( Gyr instead  Gyr) and this maximum is higher.

These differences translate for the enrichment history during the first 3  Gyr shown in panel a) into a slower increase in ASC, as expected from a later SFR, while in BLI , the Fe abundance increases very rapidly, as corresponds to the earliest SFR. On the other hand, all the other models, GNE, KRU, STD and MOD behave similarly. After 4 Gyr all models show practically the same results, with ASC reaching the highest Fe abundance despite starting to increase later than them all. Differences among models are, however, smaller than errors. In panels c) and d) differences appear as a consequence of these distinct early SFR, with BLI showing slightly higher [/Fe] at the poor-metal end of the correlation, while ASC showing the lowest. Above  dex, lines cross, ASC being the closest to data at face value. The corresponding MDFs are not equal, with ASC showing a larger proportion of metal-rich stars and all the other models showing a maximum higher than observed around the solar abundance. ASC is the only one that does not show a high peak at solar as the other models do. However, some of them, as BLI, reproduce better the metallicity distribution in the range .

In Fig. 6 we show the corresponding radial distributions of the surface density of both phases of gas, stars and SFR. These radial distributions result very similar for the six models within a radius  kpc (just where the optical radius is defined) with all of them showing a good fit to the radial profiles and it is not possible to distinguish the different lines in panel b), c) and d). The reason for this is that the same total mass is chosen for the six models and the infall rate is also the same for all of them; simultaneously, the SFR is quite similar too. There are, however, important differences out of this galactocentric distance. This occurs because, at variance with our classical models STD and MOD in which the SFR begins at the initial time , the other four models require a threshold value for the total or stellar density or the metallicity to initiate the formation of molecular clouds and the subsequent SFR. In fact, we have added a threshold metallicity to the Z term in GNE, KRU and ASC models, to avoid the parameter to be always zero (and the SFR too) thus preventing any evolution. This threshold translates into an abrupt increase of in the last 4 models when the required conditions are reached and, simultaneously, a strong decrease of the star and gas surface densities in the outer regions of disks in different degrees, due to low densities or metallicities. Panel a) showing the diffuse gas surface density is where the largest differences arise among the different models. This is expected since they differ in the prescriptions to convert this diffuse gas into molecular gas. At any rate, all of them show a flatter radial distribution than the other three panels, with ASC showing the best agreement with data up to a galactocentric distance of 14-15 kpc, while the others look slightly high, in the limit of the error bars. The radial distributions resulting from models KRU, GNE and ASC, on the other hand, show a maximum just at  kpc, that we interpret as due to a strong decrease in the parameter and that can also be seen in the decreasing of the distributions in panels b), c) and d). This way, outer than 15-16 kpc, BLI shows a better agreement with data than ASC and also better than KRU, and GNE. Nevertheless, we should bear in mind that the large errors in the data corresponding to the outer regions do not allow to discriminate among the different models. In panel b), the comparison is clearer: all models behave basically in the same way within the optical disc. At face value, BLI and ASC fit better the observations in the outer disk, while STD and MOD show values higher than expected, and KRU and GNE show a very steep decrease. In any case, it is again necessary to take into account the large error bars of these estimates. In c) all models are plausible up to a galactocentric distance of 14-15 kpc and differences arise out of this radius. As there are no data out of 16 kpc any model seems valid. In panel d) the data points lie between BLI and ASC models and STD and MOD. Summarizing, GNE and KRU yield the worst agreement with data, while BLI and ASC seem to provide the best one.

It is necessary to take into account, when analyzing these results, that all the new prescriptions, that is, GNE, BLI, KRU and even ASC models, have some parameters taken as constant which actually might be changed, thus modifying the results. Thus, GNE uses a parameter which we take as 1, just because we are calculating a model for MWG. This parameter, however, could be slightly different along the galactocentric radius, since it depends basically on the ultraviolet stellar radiation at 1000 Å  which, in turn, would change with the number of massive stars born in each time and region. With a lower value, , (see Fig. 3), a high value of may be reached at low densities of gas. Therefore, the SFR will be higher, since it depends on the molecular gas, and also the stars; simultaneously, a strong decrease of the diffuse gas will be produced. On the contrary, with a higher value, as , the molecular clouds and stars formation will decrease and, therefore, HI density will maintain high. Using a parameter variable with radius would do possible to fit the outer regions of the Galactic disk better than with a constant value. We have checked this idea with in  kpc and a straight line from kpc and . We have seen that effectively the fit improves in the outer disk beyond 14 kpc. It will be necessary to analyze carefully this subject and its consequences in a future work, after calculating the spectral energy distribution of each time and region, and the corresponding variation of the intensity of 1000 Å  in order to include a more precise dependence on this parameter. Similarly it occurs with other quantities included in the other prescriptions, as the local number density of dust grains in ASC, calculated with a fiducial value while other authors use the half of this value (Wolfire et al., 2008), or the functions or included in KRU, which are obtained through approximations and simplifications by using typical values of the number density, dust cross section, and H formation rate coefficient, to compute a value , which may also be different in other times or other radial regions of the MWG. These considerations do that the proportionality factors included in Eq.(37), (43), (63) and (64), might vary a factor between 2 and 5. This implies that some uncertainties are associated to these model results, and therefore to select the best one is not a straightforward task.

### 3.2 The evolution of the radial distributions

We plot in Fig. 8, Fig. 9 and Fig. 10 the radial distributions of star formation rate surface density, stellar profile surface density and the oxygen abundance , obtained for the six models in different redshifts/evolutionary times. To transform our evolutionary times, in redshifts, , we have used the relationship between and as given by MacDonald (2006), as explained in Mollá (2014), using a flat cosmological model with and (Planck Collaboration et al., 2014). Moreover, we calculate the evolution for 13.2 Gyr, by assuming an Universe age of 13.8 Gyr, that is, there is a time shift to start the galaxy formation. We see in Fig. 8 to Fig. 10 that the models differ mainly at the outer disk, , for . However, when we look at the higher redshifts distributions, differences are larger at shorter radius. In particular, at each radial distribution seem different of the others. This is particularly clear for the O/H abundance radial distributions.

However, if we represent as a function of the normalized radius , Fig 11, where the effective radius, , is defined as the radius which encloses the half stellar mass of the galaxy, we see that all models coincide until , where, by definition, . This occurs for all redshifts until . For redshifts higher that this value it is difficult to estimate the effective radius, which is smaller than 1 kpc in most of cases. Since we compute our models with a radius step or even 2  kpc, our results are not longer accurate enough.

The fact that the radial gradient be the same when it is measured as a function of a normalized radius related with the stellar population existing in each time is reasonable, since stars are who produce the elements and eject them to the interstellar medium. Therefore, a similar profile for the region where the stars have been created in each time is expected for all models. Differences must be appears when the star formation is not longer the same: in the outskirts of galaxies when the effect of a threshold, which do not appear equally in all models, is evident.

The effective radius, as shown in Fig. 12, is different for each model, showing the fact that the SFR history is not the same for every one, as due to the distinct time scale to create the molecular clouds obtained from the different prescriptions described in Section 2. In this figure, it is evident that each model creates the stellar disk to a different rate, the ASC model being the one that does it more later, compared with the others, and being also the one more in agreement with the observational data from Trujillo et al. (2007); Buitrago et al. (2008), shown as dots and squares. These authors obtain their data from CANDELS for galaxies with stellar masses  M. Although they divide their sample in two bins with disks and spheroids, respectively, and we use the ones for the disk galaxies, the effective radius is computed including the bulges, and, moreover, are biased towards massive early galaxies, and maybe for this reason their averaged effective radii are shorter than those we find for a MWG-like model (whose stellar mass is within their data range). The use of the new prescriptions in models GNE, BLI, KRU and ASC produces disks more slowly than the use of a free parameter, as in models STD and MOD in which the formation of the disk occurs very quickly (); in GNE and KRU models this happens at and in ASC and BLI at . On the other hand, GNE and KRU cross at and GNE also crosses with two other: ASC and BLI, at . GNE reaches the smallest at . BLI and ASC show a very similar behaviour in this figure.

### 3.3 The relation H2/HI as function of the radius

In this section we are going to compare our models with the correlation found by Bigiel et al. (2008); Bigiel & Blitz (2012). In the first work, a correlation between both components H and HI as a function of the normalized radius is found for seven galaxies. Following their findings:

 logΣH2ΣHI=1.0−0.977×R/Reff, (68)

And a similar relationship was also found by Leroy et al. (2008), for which

 logΣH2ΣHI=1.025−0.789×R/Reff, (69)

In Fig. 13, panel a), we show this radial distribution of the ratio , using the normalized radius , for the six models compared with the above expressions. In panel b) we show a similar figure where this ratio is represented as a function of the stellar surface density, compared with the correlation also found by Leroy et al. (2008): . In both case we have include the observational data corresponding to MWG as full black dots.

Clearly our old prescriptions as included in STD and MOD models are not as good as the new ones, since they produce very flat radial distributions, even increasing at the outer regions, in disagreement with the MWG data at face value (although still within errors, which are quite large in these outer regions of the disc). They also differ in the observational trends shown in Bigiel et al. (2008); Leroy et al. (2008). Models GNE,BLI, and KRU have a better behavior, but actually only ASC shows a shape closer to the correlations found by these authors. The same can be seen in the figures at the right panel. At variance with models STD and MOD using the old efficiency to form molecular clouds, models computed with the new prescriptions show an increase in the ratio with the stellar surface density. However, the slope of these relationships are smaller than the observed one, ASC being the closest model to that line, although it does not present a straight line in the logarithmic scale but a curve.

## 4 Conclusions

1. We have computed six different Galactic chemical evolution models in which different prescriptions to form molecular gas from the diffuse gas are used.

2. We have checked that the Solar Region and the Galaxy disk data (referring to the radial distributions of diffuse and molecular gas, stellar profile and star formation rate surface densities) are reproduced as shown in Fig. 5and 6. However, it is evident that the ASC prescriptions give a model in a better agreement with the data.

3. For what refers to the C, N and O abundances (Fig. 7), all models are able to fit the observations until the optical radius  kpc. After this radius, in the outer regions of disk, differences among models arise clearly, with STD and MOD showing more extended disks, while GNE shows the smallest ones. ASC, as BLI, shows an intermediate behavior. Given the data for these regions, scarce for C, and reaching far for O, it is quite difficult to determine which of these models is the best one.

4. We present the evolution along redshift of these quantities obtaining the same conclusion about the size of the disk: STD and MOD show more extended disks for all redshifts, while ASC and BLI have the smallest for . However, for GNE has the smallest disk. T his implies a different grow of the stellar disk in each model, being the ASC model the one which starts to increase later, although it arrives to a similar to the other models (except GNE) at (Fig. 12).

5. The Oxygen abundances show a similar slope in their disk distributions when only radial regions within are considered. It is necessary to take into account that the disk grows at a different rate in each model. In consequence the radial gradient of oxygen abundance measured within the optical disk is the same for all models and all redshifts, independently of the differences of grow existing among models (Fig. 11. This point will need more detailed analysis in a future work.

6. The ratio is better reproduced now with the new prescriptions to form molecular clouds from diffuse gas than before (Fig. 13). This occurs in the radial distribution of this ratio as a function of the normalized radius , as in the relationship of this quantity with the stellar surface density.

7. This implies that the ASC method to compute the conversion between gas phases, including a dependence on the stellar and gas densities and on the global metallicity, seems to produce good results when compared to observational constraints. This model, however, produces a worst fit to the SFR at early times. Furthermore, the N and O radial gradients seem to be better reproduced by the STD and MOD models at the outer disk. However, if we want to use a more realistic prescription to form molecular clouds than a free parameter, ASC may be considered the best choice to be used in modern numerical chemical evolution models as well as in cosmological simulations.

## 5 Acknowledgments

We acknowledge the anonymous referee for very helpful comments which have improved this manuscript.This work has been supported by DGICYT grant AYA2013-47742-C4-4-P and AYA2013-47742-C4-3-P and