© ESO, 2011

1. Introduction

Modeling the gas chemistry is an important step towards correctly describing the growth of cosmological structures, the formation and evolution of galaxies, and star formation in general. For instance, the molecular hydrogen is one of the most efficient coolants, and its abundance eventually determines the total amount of stars in the Universe. Structure growth and galaxy formation and evolution are customarily investigated by means of large numerical simulations in which a wide set of chemical reactions taking place in the ISM should be considered to get and follow the key molecules (elemental species in general) eventually governing the efficiency of the star formation and gas cooling. However, we must face the growing standard complexity of a typical NB-TSPH model that includes particles of dark matter, particles of baryonic matter (this in the form of stars and gas, in turn divided into several thermal and chemical phases, such as (i) cold, warm and hot, (ii) atomic and molecular, (iii) neutral and ionized), sources of energy heating and cooling, energy feedback, and easily many other physical processes. For this reason, too, a detailed chemical description of the ISM would drastically reduce the computational performances of any numerical algorithm (code) that one may adopt to this purpose. This requires a strategy for optimizing the chemical accuracy of the ISM model and the computational speed.

In this paper we present a new model of the ISM and the associated code we have developed to explore the ISM properties over wide ranges of the physical parameters and, at the same time, to cope with the above difficulties. The model and companion code are named robo¹.

The model deals with an ideal ISM element of unit volume, containing gas and dust in arbitrary initial proportions, whose initial physical conditions are specified by a set of parameters, which is allowed to evolve for a given time interval. The history leading the element to that particular initial physical state is not of interest here. The ISM element is mechanically isolated from the host environment; i.e. it does not expand or contract under the action of large-scale forces. It can, however, be interested by the passage of shock waves caused by physical phenomena taking place elsewhere (e.g. supernova explosions). Furthermore, it neither acquires nor loses material, so the conservation of total mass applies, even if its chemical composition can change with time. It is immersed in a bath of UV radiation generated either by nearby or internal stellar sources and in a field of cosmic ray radiation. It can generate its own radiation field by internal processes and so it has its own temperature, density, and pressure, each related by an Equation of State (EoS). If observed from outside, it would radiate with a certain spectral energy distribution. For the aims of this study, we do not need to know the whole spectral energy distribution of the radiation field pervading the element, but only its UV component. Given these hypotheses and the initial conditions, the ISM element evolves toward another physical state under the action of the internal network of chemical reactions changing the relative abundances of the elemental species and molecules, the internal heating and cooling processes, the UV radiation field, the field of cosmic rays and the passage of shock waves. In view of the future applications of this model in dynamical simulations of galaxies, the integration time interval is chosen in such a way that (i) it is long enough to secure that the secular evolution of the gas properties is achieved, and (ii) it is short enough to secure that the physical properties of the ISM at each instant are nearly independent of the external variations in the large-scale properties of the system hosting the ISM element. In other words, robo associates a final state (another point in the same space) to any initial state (a point in the multidimensional space of the physical parameters) through a path. The model is like an operator determining the vector field of the local transformations of the ISM from an initial state to a final one in the space of the physical parameters. This is the greatest merit of this approach, which secures the wide applicability of the model.

The new ISM model stands on recipes of the internal physical processes falling in between those developed by Glover & Jappsen (2007) and Smith et al. (2008). The first one follows the thermal and chemical evolution of the low-metallicity gas in large numerical simulations. The chemical network includes a detailed treatment of H and He but neglects the molecules formed with heavy elements as the CO molecule. Dust is included to compute its contribution to the formation of molecular hydrogen, but its evolution is not calculated. The chemical code is an on-the-fly routine, running as part of a wider code following cosmological simulations of structure growth.

The model proposed by Smith et al. (2008) uses the non-equilibrium treatment for hydrogen-like species and the standard equilibrium approximation for all the remaining chemical species. It does not take any type of dust into account. To calculate the cooling rates, Smith et al. (2008) use CLOUDY (Ferland et al. 1998) and get the cooling rates as a function of temperature, density, and metallicity. By doing this, it is possible to include a large chemical network and a wide set of coolants, but the price to pay is that several oversimplifications of the problem are mandatory, e.g. the assumption of ionization equilibrium. A similar approach has been proposed by Hocuk & Spaans (2010) using the Meijerink & Spaans (2005) PDR model instead of CLOUDY.

Our ISM model and associated code robo can not only describe the gas evolution in great detail but also includes large chemical networks and the presence of various types of dust which follow the chemistry and the complex interplay between grain destruction and formation. Many gas and dust components are taken into account, among which we recall the molecular hydrogen and the metal coolants (Santoro & Shull 2006; Maio et al. 2007) or HD (McGreer & Bryan 2008). To track the evolution of these components, the ISM model and robo take various physical processes into account that may affect the behavior of the whole system. For instance, dust is very efficient in forming both H₂ and HD (Cazaux et al. 2008). Including dust formation and destruction is a formidable task. Among other processes, dust can be destroyed by shocks that deserve special care to be properly modeled. Here we have considered two different approaches. The first one makes use of a mean shock speed that is assumed to be the same for all the gas particles, neglecting the effects of the environment. In reality this is too crude a description. The second approach starts from the notions that the shocks develop when the motion of the gas particles becomes turbulent and that the distribution of turbulent velocities obeys the one predicted by the Kolmogorov law. This is significantly better than the previous case, so it is the approach we prefer.

Finally, we would like to include the gas model (and results of robo) in numerical simulations of galaxies in a simple way. We suppose that a numerical code calculates the formation and evolution of a model galaxy from an initial stage at high redshift using the standard NB-TSPH technique. At each time step, it requires an update of the chemical status of the gas particles. This is done for every gas particle of the simulation (typically from 10⁴ up to 10⁷ particles, depending on the simulation under consideration). We have two methods to our disposal. The first one is a real-time chemical updater. This approach must use simple physics and a powerful computer in order to save computing time. The second one is to use model grids, calculated in advance for a wide range of the input parameters, in such a way as to cover the plausible space of the initial conditions. Since increasing the parameters also increases the space dimensions of the grids, thus making data handling cumbersome, we make use of the ANNs technique to get rid of this difficulty. Once the ANNs are instructed to reproduce the (robo) results as a function of the parameters, they should replace robo in the NB-TSPH simulator of galaxies. In our case the NB-TSPH model is EvoL by Merlin & Chiosi (2006, 2007) and Merlin et al. (2010). The use of ANNs can greatly improve upon this point of difficulty. Here we briefly touch upon this problem and leave the detailed discussion of it to a forthcoming paper (Grassi et al. 2011).

The plan of the paper is as follows. In Sect. 2 we give a detailed description of the physics behind the ISM model and robo. This section is divided in three parts: the chemistry of the gas phase, the presence of dust grains of different types (including their formation by accretion and destruction), and finally the heating and cooling processes. In Sect. 3 we describe the general characteristics of the code. Section 4 is dedicated to describing the results of the calculations run to validate robo. Section 5 describes how to include the results of robo in the NB-TSPH simulations. Finally, some concluding remarks are presented in Sect. 6.

2. Physical model of the ISM

In this section we describe the physics of the ISM model that is used by robo. The problem and associated code are divided into three parts, mirrored by the structure of this section: gas chemistry, dust, and cooling. All these aspects are mutually coupled and overlapped.

2.1. The chemical network for the ISM

The kernel of the numerical code modeling the gas chemistry is the network of chemical reactions among different elemental species (atoms and molecules, neutral or ionized dust grains of different types) and free particles (such as the electrons), the network of photochemical reactions between the above elemental species, and the radiation field. In addition to this, we must include all the processes creating and destroying molecules and dust grains with particular attention to those that are most efficient for cooling the gas. For this reason our chemical network follows species like H₂, HD, and metals (C, O, Si, Fe, and their ions).

The number of reactions depends on how many species are tracked and how many details are included in their description. A number of codes study the behavior of the ISM (Cen 1992; Katz et al. 1996; Galli & Palla 1998; Anninos et al. 1997; Glover & Savin 2009), each of which has a different number of species to follow and a different degree of sophistication for the physical processes taken into consideration.

We keep track of the following 27 elemental species or molecules plus the free electrons: H, H⁺, H⁻, , , , , , , , , , , , , , , , , , , , , , , , , and .

These species are divided into four groups. The first one contains hydrogen-only based species. The second group is composed of deuterium-based species, in which HD plays the key role in the gas cooling. The third group lists helium and its ions. Carbon, oxygen, silicon, and iron with their ions and compounds form the fourth group. The free electrons link all the four groups together. The species from through are introduced to follow the formation and destruction of to be described below.

The reactions in which all the above species are involved are

• collisional ionization (A + e⁻ → A⁺ + 2e⁻),

• photo-recombination (A⁺ + e⁻ → A + γ),

• dissociative recombination (),

• charge transfer (A⁺ + B → A + B⁺),

• radiative attachment (A + e⁻ → A⁻ + γ),

• dissociative attachment (A + B⁻ → AB + e⁻),

• collisional detachment (A⁻ + e⁻ → A + 2e⁻),

• mutual neutralization (A⁺ + B⁻ → A + B),

• isotopic exchange ()

• dissociations by cosmic rays (AB + CR → A + B),

• neutral-neutral (AB + AB → A₂ + B₂),

• ion-neutral (),

• collider (AB + C → A + B + C),

• ionizations by field photons (A + γ → A⁺ + e⁻)

where A and B are two generic atoms. To these reactions we must add those regulating the abundance of to be described below.

The chemical network governing the ISM model is a classical system of differential equations in which each equation is a Cauchy problem of the form (1)where n_i(t) is the number density of the ith species with known initial value n_i(0). In Eq. (1) R_lm(T) is the rate of the reaction between the lth and the mth species expressed in units of cm³/s. Equation (1) is written as the sum of all the reactions forming the ith species (∑ _lmR_lm(T)n_l(t)n_m(t)) minus the sum of all the reactions destroying the ith species (∑ _jR_ij(T)n_i(t)n_j(t)). The indices i,j,l, and m run from 1 to 28 as in the list of elemental species and free electrons already mentioned. The one-to-one correspondence between the indices in the Eq. (1) and elemental species is given in Table 1.

To clarify the meaning of Eq. (1) we show the case of a two-reaction system (2)where A_i is the number density of the generic species in units of cm^-3. Looking at the density variation of the species A₀ over the time step dt the density n_A0 changes as (3)In this equation the positive term is due to the second reaction (between A₄ and A₅) that increases the total quantity of A₀ (and also A₃), while the negative term is due to the first reaction in which A₀ is destroyed because it reacts with A₁ to form A₂ and A₃.

We adopt here a minimal description containing only the 28 species and/or the molecular compounds listed in Table 1, the 64 reactions listed in Table 2 (where reactions among hydrogen, deuterium, and helium are considered), in Table 3 (where metals and cosmic rays are involved), in Table 4 (where the reactions of the CO mini-network are listed), and finally, the 12 photochemical processes listed in Table 5. In all the tables the gas temperature T is in Kelvin; the electron temperature T_e is in eV. In Table 3 the cosmic ray field is ζ_CR s^-1 and it is the rate of H₂ ionization by cosmic rays. For Table 4 we have T₃₀₀ = T/300. The expressions in Table 5 are ξ = E/E_th, , , x_f(a,b) = E/a − b, and x_Si = 1.672 × 10^-5.

The various types of reaction are described below in some detail. Tables 2–4 also sample the chemical reactions according to the same four groups in which we separated the elemental species depending on the type of reaction and the information at our disposal to derive the reactions rates. In doing this, the notation in use may appear rather complex. This is done on purpose, and the main motivation is that we want to keep the same formalism as adopted in the literature for easy comparison. First of all, there are 64 reactions among particles, whether atoms or molecules or free electrons. Each reaction is identified by a progressive number from 1 to 64 and so the associated reaction rate named k_n (with n from 1 to 64). Each reaction will be a term in the righthand of the 28 differential equations of system (1) and the associated reaction rate will coincide with one of the R_ij terms of Eqs. (1). To establish the exact correspondence is a matter of patient work, which is not the prime interest here. Most of these reaction rates are taken from Abel et al. (1997) and Glover & Savin (2009) to whom we refer for all the details. Some of these reactions are discussed in this section, and we pay attention to those rates that are vividly debated in literature.

Minimal reaction network for CO. Since this molecule plays a key role in determining the properties of the ISM, we have included a simplified network of reactions taken from Nelson & Langer (1997) to follow the creation/destruction of the CO in detail. These reactions are listed in Table 4 together with the corresponding rates of Woodall et al. (2007). Note that T₃₀₀ = T/300, with T in K.

Table 4 shows how CO is formed from C⁺ and O species, using CH, CH₂, and as intermediate products/reactants, and H₂ and e⁻ as main reactants.

Our model includes UV radiation, so we must consider the photodissociation of the molecules included in the network. The main effects of the radiation are on the carbon-based molecules and the corresponding photo-ionization rates are (4)and (5)where G₀ is the ratio of the adopted UV flux to that by Habing; i.e. I_Hab = 1.2 × 10^-4   erg/cm²/s/sr. The expression indicates that we are dealing with the CH, CH₂, , and molecules.

Basing on Eq. (5), we must consider the products of all the molecules present in reactions like CH_x^(+ ) + γ → products. However, at this stage we are not interested in the detail of these reactions; therefore, following Nelson & Langer (1997), we introduce the parameter β(G₀) controlling the efficiency of the reactions from 58 through 61 in Table 4. This parameter is defined as β(G₀) = exp [G₀·ln(ξ)]  with ξ = 5 × 10^-10 from Eq. (5). In the absence of the UV radiation (G₀ = 0), we get β = 1: this means that the existing photons cannot affect the formation of the various . Instead, when G₀ = 1, it follows β = ξ, which means that the formation reactions have reduced their efficiency. Finally, when G₀ → ∞, (unphysical) then β → 0, which means that the molecules are completely destroyed before they can interact. The coefficient of the ith reaction now is , with i ∈  [58,61]  and Eq. (5) is not included in the code.

Photochemistry. In addition to this, we consider the group of photochemical processes listed in Table 5 for which we provide the cross section of the reaction and the analytical expression to derive the reaction rate as a function of the existing radiation field.

In the present model of the ISM we do not include the chemical reactions with lithium and its compounds (e.g. LiH and LiH⁺) because according to Prieto et al. (2008) and Mizusawa et al. (2005), they are not important coolants. After we neglecting lithium and compounds, the number of species to follow, including metals and deuterium compounds, amounts to the 28 we have considered.

There is some interest in the CO molecule because it is an important coolant at very low temperatures. Among others, recipes for the formation of the interstellar CO are presented by Ruffle et al. (2002) and Glover et al. (2010). However, explicitly following the formation of the CO molecule would be too complicated for our aims. See below for the cooling rate by the CO molecule.

In our chemical network we also include the effect of the cosmic ray radiation field. For convenience, we list in Table 3 the ionization rates for cosmic rays given by Walmsley et al. (2004). Cosmic rays are important for the cooling because they can destroy molecules that contribute to it, thus influencing the temperature of the interstellar medium (see the reactions shown in Table 3). In particular, cosmic rays are able to destroy H₂ and HD and to ionize atoms. The role of the cosmic rays is still largely unknown because the strength of their radiation field in various regions and epochs of the Universe is not firmly determined. While their effect can be explored in the vicinity of a strong source, for a random sample of the Universe there are no exact measurements of the cosmic ray radiation field. Furthermore, most NB-TSPH simulators like EvoL do not contain the full description of cosmic rays, but simply consider their presence as a parameter. Therefore, even if these reactions have been considered and their rates presented in Table 3. they have not been included in practice in the system of Eqs. (1).

Photochemical reactions. For the photochemical reactions we adopt the model and rates of Glover & Jappsen (2007) and Verner & Ferland (1996). In general the reaction rate is given by (6)where J(E) = J(hν) is the energy flux per unit frequency and solid angle of the impinging radiation field, σ(E) is the cross-section, τ(E) the gas opacity at the energy E, and f(E) a numerical factor accounting for the effects of the secondary ionization, which are negligible if the UV radiation is not dominated by X-rays (Glover & Jappsen 2007). The rate is in .

The cross sections σ(E) are expressed in two different ways according to the reaction under consideration. The reactions and associated rates are listed in Table 5. For the reactions from 1 through 8 we note that (i) the quantity ; (ii) for H₂ and HD photodissociation, we use the model proposed by Glover & Jappsen (2007). We have (7)without considering self-shielding. For the HD we have (8)For the remaining reactions involving metals (i.e. from 9 through 12), the rates contain the fits given by Verner & Ferland (1996) using the expression (9)where . The meaning of the various symbols is given in Table 5. More details on the reactions and companion quantities and symbols are given by Glover & Jappsen (2007) and Verner & Ferland (1996). See also the references therein.

The UV radiation field is calculated as in Efstathiou (1992), Vedel et al. (1994), and Navarro & Steinmetz (1997). In particular we have (10)where z is the redshift, ν_H = E_H/h is the frequency corresponding to the hydrogen first-level energy threshold, α_UV = 1 and J₂₁(z) is (11)where in our case J₂₁ = 1. Finally, to derive the reaction rates, we integrate Eq. (6) from E = E_th to E = ∞ over the energy range of the radiation field photons.

The integrals must be calculated at each time step if the radiation field changes with time or once for all at the beginning of the simulation if the radiation field remains constant. The integrals are calculated using Romberg’s integration method.

Heating by photodissociation. Heating by photodissociation of molecular hydrogen and UV pumping, H and He photo-ionization, H₂ formation in the gas and dust phase, and finally ionization from cosmic rays are modeled as in Glover & Jappsen (2007). All these processes are listed in Table 6 where n_cr is the critical density and R_d is the photodissociation rate. To describe the heating by the photoelectric effect we adopt the model proposed by ? and discuss it separately in Sect. 2.4.1 below.

The H and He photo-ionizations increase the gas energy at the rate (12)where σ(E) is the reaction cross-section, J(E) is the photon flux, e^−τ(E) is the optical depth for a photon of energy E, η is the efficiency of the process (i.e. the fraction of energy converted to heat), and (E − E_th) is the difference between the energy of the photon and the energy of the atomic ground level. The rates Γ are given in Table 6.

2.2. Reaction by reaction: notes on a few cases

In this section we examine in some detail a few chemical reactions that are widely debated and are affected by large uncertainties. The reaction numbers are the same as those used in Tables 2, 3, and 5 for the sake of an easy identification.

(reaction 9). Each author has his own favored rate for this reaction so that the overall uncertainty is large. Our prescription is as follows. Up to 3 × 10⁴  K we adopt the data of Savin et al. (2004), based on the latest and most accurate quantum-mechanical calculations of the vibrationally resolved cross sections for the charge transfer at the center-of-mass collision energies from the threshold (~1.8 eV) up to 10 eV. For higher collision energies, up to approximately 10⁴ eV, we derive the cross section following the suggestions by Barnett et al. (1990) and Janev et al. (1987), which stand on the best known experimental (see for instance Gealy & van Zyl 1987) and theoretical data, as well as on the recent measurements by Kusakabe et al. (2003). These data were smoothly matched to those of Krstić (2002) at low energies, thus yielding the most updated cross section for the charge transfer from H to in the vibrationally ground state. This cross section is shown in the top panel of Fig. 1. Following Savin et al. (2004), from this cross section we calculate the rate coefficients for temperatures from the threshold up to 10⁸ K, as shown in the bottom panel of Fig. 1.

We fit the reaction rate (in units of cm³   s^-1) with the analytical expression (13)where T is the temperature in Kelvin. Two intervals are considered for the temperature, i.e. T =  [10²,   10⁵]   K and T =  [10⁵,   10⁸]   K and two different fits are derived. The coefficients a_i for the two fits are given in Table 7. The rates obtained from this analytical fit are compared to those derived from the numerical data. No difference can be noticed as shown in the bottom panel of Fig. 1.

Associative detachments of H⁻ and D⁻ (reactions from 21 to 24). These processes include the following four reactions Following Glover et al. (2006) and using their notation, the reaction rate is expressed as k = const. × ξ where ξ is a parameter taking values in the interval  [0.65,5.0] . In our case we obtain k₂₁ = k₂₄ = 10^-9   ξ and k₂₂ = k₂₃ = 10^-9   ξ/2. In robo the parameter ξ can be varied in the above interval to investigate its effects on the overall results. In the present study we adopt ξ = 0.65.

H⁺ + H⁻ → H + H (reaction 29). For the mutual neutralization of H⁻ and H⁺ we adopt the cross section given by Croft et al. (1999) and Glover & Abel (2008). The mutual neutralization rate is related to the rate of the associative detachment H⁻ + H → H₂ + e⁻, in competition for the available H⁻ ions. Both rates are important for the formation of H₂ (see Glover & Abel 2008, for more details). It is worth recalling here that other estimates for this reaction rate have been reported by Glover & Abel (2008). Moseley et al. (1970) proposed the rate k = 5.7 × 10^-6T^-0.5 + 6.3 × 10^-8 − 9.2 × 10^-11T^0.5 + 4.4 × 10^-13T cm³   s^-1, whereas Dalgarno & Lepp (1987) gave k = 7.0 × 10^-7T^-0.5 cm³   s^-1.

He⁺ + e⁻ → He + γ(reaction 35). This process can occur either via direct radiative recombination or di-electronic recombination, followed by radiative relaxation. Therefore, the reaction rate is the sum of two terms (14)in units of cm³   s^-1. The first term is the di-electronic recombination rate, described by (15)where T_e is the temperature expressed in eV. The second term is the radiative recombination rate whose temperature dependence is (16)Both terms are taken from Abel et al. (1997).

Z⁺ + e⁻ → Z + γ(reactions 38, 40, 42, 44). An important process to include is the metal recombination. The metals considered by our ISM model and robo are C,Si,O, and Fe. The recombination rates of the metals are calculated according to the formalism proposed by Verner & Ferland (1996)(17)in units of cm³s^-1, where T is the temperature, and . In Table 3 we give the fit coefficients A, b, τ₀ and τ₁ for each metal.

Z + e⁻ → Z⁺ + 2   e⁻(reactions 39, 41, 43, 45). For the collisions between electrons and metals we use the fit proposed by Voronov (1997) expressed by (18)in units of cm³s^-1, where U = ΔE/T in which ΔE is the energy difference between the two atomic levels involved in the process, and T is the temperature. Both ΔE and T are in eV. The parameters A, P, ΔE, X, and K are given in Table 3.

2.3. The dust

Dust grains take part to the process of molecule formation, e.g. HD and H₂ form on the surface of dust grains; therefore, all the physical processes involving the dust must be described and treated with the highest accuracy. To understand how grains take part in the process of molecule formation, we need to know the mechanisms of grain formation and destruction, along with the distribution of the grain temperature and size.

Given an initial set of dust composition and abundances, our dust model follows the evolution of the dusty components during the history of the ISM due to the creation of new grains and the destruction of the existing ones. Creation of new dust grains is governed by several processes that have different efficiencies depending on the size of the dust particles. The same applies for the destruction that is mainly due to the thermal motion of the gas particles and shocks from supernovæ. Thermal destruction is quite easy to model, because the only parameter at work is the gas temperature. Shock disruption is more difficult to evaluate. The main uncertainty comes from discrete numerical hydrodynamical simulations only being able to follow shocks up to a given (often too coarse) resolution, yet insufficient for the microscopic description required here. To cope with this, we have followed a “statistical” approach. In brief, once identified the gas particles of the NB-TSPH simulations with their turbulent velocities (suggested by their velocity dispersion), we assume that all the shocks inside them follow the Kolmogorov law. See below for more details.

Furthermore, grain destruction may depend on their size. This is the case of the destruction by thermal motions, where the small grains are destroyed before the large ones. As a consequence of this, the size distribution function and the abundances of dust grains of different type are continuously changing with time. The ever-changing size distribution function plays an important role in the formation mechanism of HD and H₂, which are among the most efficient coolants. In brief, changing the distribution function of the dust grains means changing the quantity of key-role molecules produced by dust, as shown in Cazaux & Spaans (2009).

We analyze here the different aspects of the grain evolution and their role in the formation of coolant molecules. First, we focus on the distribution function of the dust grains, then we describe the formation of coolants on grains. Finally, we analyze the destruction and the formation of dust, as well as the effects of the grain temperature.

2.3.1. Size distribution function of dust grains

We adopt a simple power-law, MRN-like, grain size distribution function (Mathis et al. 1977; Draine & Lee 1984). This is given by dn(a)/da ∝ a^−λ, with a the grain dimension, n(a) the corresponding number density of grains with dimension a, and λ = 3.5; this distribution covers the range 5  Å < a < 2500  Å and is extended to the smallest grain dimensions in such a way as to include the typical polycyclic aromatic hydrocarbons (PAHs) region (Li & Draine 2001b). Even if more complicated distributions have been proposed (Weingartner & Draine 2001a), our simple choice is suitable as an initial condition for the purposes of this work.

Indeed, it is worth pointing out here that the initial size distribution function changes with time owing to dust destruction and formation, which in turn depend on the grain size. Consequently, the power law we started with does no longer apply. This will soon affect the formation of H₂ or HD. For this reason, the knowledge of the current size distribution function is paramount.

First of all, we split the interstellar dust in two main components: the carbon-based (thereinafter simply carbonaceous grains) and the silicon-based grains (thereinafter the silicates). These play different roles in the formation of molecular hydrogen, as well as in different formation and destruction rates. In principle, there should be an additional group to consider, namely the PAHs, because the three types of dust have different efficiencies in the H₂ and HD formation. However, there is the lucky circumstance that PAHs and graphite grains have similar efficiencies to those shown by Fig. 2 of Cazaux & Spaans (2009). Basing on this, we lump together graphite grains and PAHs and treat separately the silicates. Therefore, thereinafter we refer to the mixture graphite grains plus PAHs as “carbonaceous grains”.

2.3.2. Dust-driven H₂ and HD formation

The model we adopt to describe the so-called dust phase is the one proposed by Cazaux & Spaans (2009). In brief, the dust phase has its own network of reactions that establish the relative abundances of the various types of dust and govern the formation of H₂ and HD. We can neglect the presence of other molecules like CO. The newly formed H₂ and HD are fed to the gas phase and enter the system of equations governing the abundances of the ISM and vice versa. In reality, the two phases should be treated simultaneously, governed by a unique system of differential equations determining the number densities of the elemental species, dusty compounds, and free electrons at once. In practice this is hard to do because in general the time scales of the various processes creating and destroying the elemental species of the gas phase are much different from those governing the formation/destruction of dust grains. To cope with this difficulty, the Cazaux & Spaans (2009) model separates the two phases, provides the results for the dust phase, and more important provides a link between the two phases that allows determining the formation rate of H₂ and HD by dust as a function of the gas temperature T_g and dust temperature T_d. This link is named the “Sticking Function S(T_g,T_d)”. The sticking function somehow quantifies the probability that an hydrogen atom striking a grain sticks to the surface rather than simply bouncing off. The sticking function characterizes the probability for an atom to remain attached to a grain. Here we strictly follow this way of proceeding.

According to the Cazaux & Spaans (2009) model, the rate at which H₂ molecule forms over the grain surface is (19)where n(H) is the density of hydrogen in the gas phase, is the gas thermal speed (k is the Boltzmann constant and m_H is the hydrogen mass), n_d the dust number density, σ the grain cross section (i.e. πa²), ϵ_H2 the intrinsic efficiency of the process, and S(T_g,T_d) the sticking function. This latter is in turn a function of the dust and grain temperatures (20)where T_g and T_d are in Kelvin. With this function the probability for a gas molecule to remain stuck on a dust grain is higher in a cold gas than in a hot one. This probability is even higher for the cold grains. For HD the expressions are similar.

Equation (19) contains the intrinsic efficiency ϵ_H2 (or as ϵ_HD for HD), which is the probability for the process to occur. In our case it is the probability that atoms, which are stuck to a grain, react to form H₂ or HD.

For the carbonaceous grains, the efficiencies ϵ_H2 and ϵ_HD coincide and are equal to (21)The efficiency for the silicates ϵ_Si is (22)where β_d = 4 × 10⁹ for H₂ and β_d = 5.6 × 10⁹ for HD. The term ℱ is a function of the gas temperature and can be written as (23)where T_H is given by the expression (24)The various quantities appearing in the above relationships are listed in Table 8. For more details on these equations see Cazaux & Spaans (2009). Comparing Eqs. (21) and (22), we see that the efficiency is high when the dust temperature is low. For the silicates the efficiency window is shorter than for the carbonaceous grains. For the silicates the efficiency is high for temperatures up to 20  K and then falls by two orders of magnitude. The carbonaceous grains are efficient for temperatures up to 100  K, where the efficiency is still 0.1 (instead of 0.01 as for the silicates). Finally, it is worth noticing that the efficiency profile is smoother for the carbonaceous grains.

From all these considerations it follows that cold dust and warm carbon-dominated dust in a cold gaseous environment are the most efficient drivers for the formation of coolant molecules like H₂ and HD.

2.3.3. Grain formation

Here we briefly examine the formation of dust grains. We start with an initial number density with the size distribution given in Sect. 2.3.1 and with a given ratio between the silicates and the carbonaceous grains. The latter is a free parameter varying in the range  [0,1] , where zero stands for dust made of sole carbonaceous grains; one is for dust made of sole silicates.

According to Dwek (1998), the temporal variation in the size distribution function of grains caused by accretion is given by (25)where all the quantities have the same meaning as in Eq. (19), but for α(T_g,T_d) and c_d.

The quantity α(T_g,T_d) is a sort of sticking coefficient depending on the gas and dust temperature and the type of dust. This coefficient is therefore forced to change in the course of the evolution. Furthermore, n_d (and consequentially n) are functions of a. Equation (25) is similar to Eq. (19) because the process is similar, except that now this mathematical description is applied to the carbon atoms that remain stuck to the carbon lattice of the grain. Our expression differs slightly from the original one of Dwek (1998) because the parameter c_d is introduced in Eq. (25) to take some considerations made by Dwek (1998) himself into account. Inserting c_d = 1 we can obtain the time scale of the process. In a standard cold gas this time scale is τ ≈ 2 × 10⁴  yrs, which is significantly smaller than the observational estimates. Normal evaporation, caused by cosmic rays or UV heating and grain-grain collisions can halt the growing of the dust grains. The factor c_d somehow takes this phenomena into account. It is estimated to be on the order of c_d = 10^-3. We name c_d the delay factor.

For the sticking coefficient α(T_g,T_d), we make use of the data by Leitch-Devlin & Williams (1985) and consider a carbon atom of the gas phase as impinging on a carbon lattice. The equation fitting the data has the form (26)where T_g and T_d are the gas and dust temperatures, respectively (see Fig. 2). The sticking coefficient α(T_g,T_d) has no dimensions. The data of Leitch-Devlin & Williams (1985) provide a good dependence on the gas temperature T_g, but a poor one on the dust temperature T_d. Therefore an accurate fit is not possible. The above relationship can be used in the temperature intervals 10  K ≤ T_g ≤ 1000  K and 3  K ≤ T_d ≤ 300  K. At present we also use the same model for silicate grains even though this approximation might not be accurate. We plan to improve upon this point in the future.

With this model the formation efficiency is higher when the interstellar medium has a temperature T ≈ 100  K and when the dust grains have a temperature of about 300  K. For dust grains with temperature higher than 300  K or lower than 3 K, the formation efficiency is unknown. Though the fits from Leitch-Devlin & Williams (1985) are not very accurate, they are accurate enough for our purposes because we are only interested in the shape and maxima of the curves in Fig. 2.

Finally, we note that we have already described the formation of dust by means of a general process of accretion in the ISM, leaving aside other sites of dust formation like the envelopes of obscured AGB stars, Wolf-Rayet stars, and remnants of supernovæ (Dwek 1998). The reason behind this is that all of these are external sources of dust that eventually enrich the ISM in dust content and therefore determine a different initial dust content as input for ROBO. Indeed, for our purposes, only the internal sources of dust in the unit volume are important. The information about the initial conditions of the dust mixture should be provided to ROBO by the companion NB-TSPH code EvoL, of properly taking the dusty yields coming from the stars into account.

2.3.4. Grain destruction

In our model we assume that there are three processes destroying the grains of the ISM, i.e. shocks, in particular supernovae induced shocks, vaporization by high-velocity shocks, and thermal sputtering, in which dust is destroyed by the thermal motions of the gas. We need to describe the three processes in a way suited to the NB-TSPH formalism.

Destruction by shocks. This phenomenon is difficult to model. First of all shocks may easily induce turbulence in the ISM and to be suitably described one needs some assumptions about the velocity fields of the fluid elements. The main difficulty arises from the fact that within an SPH gas particle, because of the typical mass resolution, unresolved shocks may take place at different velocities. To avoid this difficulty, we treat the gas as a turbulent fluid, assuming that even inside a gas particle the turbulent nature of the fluid is preserved, and use the Gaussian velocity distribution φ(v) derived from the Kolmogorov law for the power spectrum E(k) ∝ k^−αdk with α = 5/3. Second, grains of different masses and sizes move at different velocities: the large, more massive grains moving slower than the smaller less massive ones. This suggests adopting the point of view in which the small fast grains are considered as projectiles impinging on the massive slow grains considered as targets. To summarize, the ISM we are dealing with is turbulent because of the shocks crossing it and rich in dust grains of all possible dimensions and masses interacting with each other and with the shock fronts.

Given these premises, we model the destruction of dust grains by shocks according to the picture by Hirashita & Yan (2009). The dust grains are grouped in N = 100 bins of number density according to their mass. The mass of a grain is in turn expressed as the product of a reference density, here assumed to be the density of the graphite grains ρ_gr = 2.3   g/cm³ times the volume of the grain V_gr = 4/3πa³ where a is the radius of the grain (assumed to be spherical for simplicity). Therefore each density bin contains grains whose mass goes from m_low to m_up, where both mass limits vary with the bin. The corresponding radii are given by and . The grain radii follow the distribution law given by dn(a) = C_nor × a^-3.5da, where C_nor is a suitable normalization factor to be determined. Therefore, the number density of the grains in each mass interval is (27)The total number density of dust grains of any size and hence mass is (28)When a projectile of mass m_j hits a target of mass m_i, the target loses a fraction of mass f_shm_i if m_j < m_i/2, or it is totally destroyed if m_j ≥ m_i/2 (i.e. f_sh = 1). In both cases the projectile is destroyed. The lost mass, i.e. a fragment of lower mass, so the dimension is assigned to the mass (size) bin according to the adopted power law. The remaining part of the target (remnant) of mass (1 − f_sh)m_i is added to the appropriate mass, hence size bin, so the number density variation of the ith bin is (29)where the subscript L indicates the term representing the mass lost in fragments, F the mass gained from the fragmentation, and finally R the mass of the remnants moving from other bins to the ith one. The first term must be negative.

To estimate the various contributions to the total number density variation, we must consider the probability of an impact. This will be proportional to the number density of the targets n_i, the density of the projectiles n_j, their relative speed v_ij, and their sizes (a_i and a_j). The assumptions on the speed will be discussed later in this section. For the targets belonging to the ith bin, we find that the mass lost in fragments corresponds to the density change of (30)where (31)and the impact coefficient α_ij is (32)Consequently, the total density variation of newly created fragments is (33)They will be distributed among the bins with a power law to obtain the (dn_i/dt)_F of each bin.

The variation in the remnant number density is given by (34)where k is the index of the bin receiving the remnant given by (35)with i the index of the initial grain that suffered fragmentation and Δm the mass range of each bin, (m_up − m_low)/N.

The shock velocities obey a Gaussian distribution normalized to (36)where v_low and v_up are the limits of the shock velocities over which the normalization of φ(v) applies, thus fixing the normalization constant. The relative speed v_ij of the grains is obtained from the velocity distribution of Eq. (36) with v_low = 1 km s^-1 and v_up = 200 km s^-1 and a total of 30 velocity bins. The velocity φ(v) is a Gaussian centered at v = 100 km s^-1 with a dispersion of 30 km s^-1. The integral (37)yields the relative weight of each velocity bin. The relative velocities of the projectiles are distributed according to these weights and the total density of targets in the ith bin changes as described by Eqs. (30) and (32).

Vaporization. Since we also deal with high-velocity shocks, vaporization of dust grains into the gas phase may become important. We use the formalism of Tielens et al. (1994) for the impact of carbon grains. The fraction of vaporized material is (38)with f_vap given by (39)where v_t = 23 km s^-1 is the shock threshold velocity, m_i and m_j are the masses of target and projectile grains, respectively. Equation (39) must satisfy the condition v ≥ v_t. The quantities f_v1 and f_v2 are taken from Tielens et al. (1994) and are given by the relations (40)and (41)where (42)with ℳ the Mach’s number and s = 1.9.

Destruction by thermal sputtering. To evaluate the fraction of grains destroyed by thermal sputtering, we adopt the approximation by Draine & Salpeter (1979). A grain of dust in a medium with temperature T ≤ 10⁶  K has a destruction time (i.e. its lifetime) with a in nm. From this we can obtain the destruction rate per second. The above relation is only valid for high temperatures.

To include the temperature dependence of the destruction rate, we refer to Tielens et al. (1994), with the aid of which we model the dependence on the gas temperature of the lifetime of the dust. As expected, dust grains with lower temperature have a longer lifetime. Finally, according to Draine & Salpeter (1979), the small grains have a shorter lifetime than the large ones in the same environment.

The grain temperature depends on the radiation field generated by all the stellar sources. However, for the sake of simplicity, in this study we used a fixed value for the grain temperature. Therefore, a big improvement would be given by including a description of the photon diffusion (e.g. Mathis et al. 1983). Our choice is partially justified by the fact that the companion NB-TSPH code EvoL does not yet include photon diffusion. This leads us to postpone the implementation of different grain temperatures to when photon diffusion is included in EvoL.

2.4. Heating and cooling

At this stage, it is worth discussing in some detail the role played by the heating and cooling during the history of star formation in a galaxy or a cosmological simulation. We have already emphasized that coolants are the key elements for the formation scenario. As the gas cools down, more and more spatial structures and stars (with a certain mass function, hence mass-luminosity law) are born. Without coolants neither structures nor stars would form. In this context, the grains and their temperature in turn play the key role in the formation of efficient coolant molecules like H₂ and HD. Cold grains form more molecules than the warm ones. Since the dust temperature depends on the surrounding photon flux, it means that stars in the neighborhood are crucial for heating the dust particles. All this forms a closed loop of interwoven physical processes: dust forms coolants – coolants induce star formation – stars heat the dust. Understanding the details of this mutual interaction will allow us to get clues on the star formation history in general.

In this scene a starring actor is the gas cooling. Chemical reactions are sensitive to the gas temperature, hence to the gas cooling; indeed, the formation of coolants depends on the gas temperature, which in turn depends on the cooling process. We split the cooling of the gas in several sources characterized by the physics of the dominant process. Above 10⁴   K, there are two very popular descriptions of the cooling, namely Cen (1992) and Sutherland & Dopita (1993). At lower temperatures we have the metal cooling (Maio et al. 2007), the molecular hydrogen cooling (Glover & Abel 2008; Galli & Palla 1998), and finally, the HD cooling (Lipovka et al. 2005).

2.4.1. Heating by photoelectric ejection of electrons from dust grains

The photoelectric ejection of electrons by dust grains is an important source of heating. The model proposed here is based on Bakes & Tielens (1994) and Weingartner & Draine (2001b). The photoelectric heating is given by (43)where W is the FUV dilution factor, σ_abs the photon absorption cross section, Y_ion the photoelectric ionization yield, F(ν) the UV radiation flux, and g(N_C,IP_Z) the kinetic energy partition function. Here, N_C is the number of carbon atoms that form the dust (N_C ~ 0.5   a³ with a the radius of the grain in Å, according to Li & Draine 2001b) and IP_Z the ionization potential of a grain of charge Z, and ν_H is the Lyman frequency and ν_Z the frequency corresponding to IP_Z. In our model the range of grain sizes goes from 5 Å  to 100 Å, since the total heating is mostly due to small grains (see Bakes & Tielens 1994,   for details).

The total photoelectric heating is (44)where f(N_C,Z) is the probability to find a grain composed of N_C carbon atoms at a certain charge Z. This can be computed by considering the collisions with electrons and ions. The detailed balance yields (45)with J_pe the rate of photoelectron emission, J_ion, and J_e the accretion rates of ions and electrons. For the detailed calculation of f(Z) see Bakes & Tielens (1994), while for all the details about J_pe, J_ion and J_e see Weingartner & Draine (2001b).

As the FUV absorption cross section σ_abs depends on the absorption efficiency of FUV photons by the grains, we use the data of Draine & Lee (1984) and Laor & Draine (1993) for the graphite grains and Li & Draine (2001a,b,c) for the PAHs. We also use the model proposed by Li & Draine (2001a,b,c) to create a mixed population of PAHs and graphite grains. Following their paper the total absorption efficiency is (46)defining (47)by q = 0.01 and a₀ = 50 Å (see Li & Draine 2001a,b,c, for more details).

In the same way we can also obtain the cooling associated with the electron recombination. This is given by (48)where n_i is the density of the charged particle with sticking coefficient s_i = 1 and mass m_i and a is the grain size. The term is described in detail in Bakes & Tielens (1994). The total cooling is obtained as in Eq. (44). In Fig. 3 we show the results from a simple model with a grain distribution with n_dust = 10^-5  cm^-3 and a gas with a temperature of 10² K.

2.4.2. Which are the sources of cooling to consider?

Cooling at high temperatures (T ≥ 10⁴). Two main sources of cooling can be found in the literature: (i) the formulation proposed by Cen (1992) which includes the following mechanisms:

• collisional ionization – H, , , ,

• recombination – H⁺, , ,

• dielectronic recombination – ,

• collisional excitation – H(all n), (n = 2), (n = 2,3,4),

• bremsstrahlung – all ions,

which is particularly suited in presence of non equilibrium reactions for the H and He chemistry; (ii) the description of Sutherland & Dopita (1993, thereinafter SD93), which is particularly suited to treat cooling in presence of metals and to describe it as function of the gas metallicity. Their case for Z = 0 can be left aside. Both cooling models are used for temperature higher than 10⁴   K. We refer to them by the acronyms C92 or SD93. In our models of the ISM we adopt both sources as appropriated to the current physical conditions.

A difference between C92 and SD93 is that the cooling rate is described in the former by numerical fits (no further interpolations are needed), whereas in the latter one has to interpolate huge tabulations of data in temperature and metallicity. To this aim we used a surface fit (routine SFIT in IDL) where the plane surface passes through four points that correspond to two discrete values in temperature and two discrete values in metallicity. The analytical expression of this plane yields the cooling rate. Given the two pairs of interpolating points, (T₀,Z₀), (T₀,Z₁), (T₁,Z₀), and (T₁,Z₁), the analytical function that describes the surface can be written in the form (49)with a_i the fit coefficients. So, for a generic point (T,Z), for which T₀ ≤ T ≤ T₁ and Z₀ ≤ Z ≤ Z₁, the rate cooling is given by Eq. (49).

H₂cooling. The H₂ cooling requires a different description. Hollenbach & McKee (1979) first evaluated the molecular hydrogen cooling. A modern, widely used source of H₂ cooling is by Galli & Palla (1998). Although this cooling function is quite accurate, we prefer to add some more details by including the function found by Glover & Abel (2008), which takes not only the H₂ − H interaction into account, but also the collisions with He, H⁺, H₂, and free electrons. It is described by (50)where T₃ = T/10³  K with T the gas temperature, a_ik is the ith fit coefficient of the kth species (k =  { H,H⁺,H₂,e⁻,He } ), and n are the number densities. The orto-para ratio is assumed to be 3:1. Outside the temperature range of the fits we use the molecular hydrogen cooling by Galli & Palla (1998).

HD cooling. To describe the cooling by the deuterated molecular hydrogen we use the model proposed by Lipovka et al. (2005). It includes HD roto-vibrational structures, radiative and collisional transitions for J ≤ 8 rotational levels, and the vibrational levels v = 0,1,2,3. It has been found that including the roto-vibrational transitions increases the cooling efficiency of the HD. The fit provided by the authors depends on the gas density and temperature. It can be parameterized as (51)where c_ij is the matrix whose elements are given in Table 9. Using Eq. (51) is particularly convenient from a numerical point of view as it provides fast evaluations of the cooling by the HD molecule.

Cooling by the CO molecule. Owing to founding hypotheses of our model of the ISM presented in the introduction, the ISM is optically thin so that it is worth considering the straight cooling of the CO molecules by radiative transitions among rotational and vibrational energy levels.

To this aim, we may adopt the steady-state equations of McKee et al. (1982) and the rates calculated by Schinke et al. (1985) by applying the L → 0 coefficients to the Goldflam et al. (1977) method and including the Depristo et al. (1979) correction (see these studies for more details).

To calculate the total rotational cooling, we use the equation of statistical equilibrium from McKee et al. (1982): (52)where n_j is the number density of the CO molecules in the state j = J, n_H2 the number density of the molecular hydrogen, γ_ij the rate coefficient of the transition, and A_j the A-coefficient from the state j to j − 1. The system of equations is completed by the condition  ∑ _jn_j = n_CO. The solution of this linear system is straightforward. The emission due to the transition from J to J − 1 is (53)from which we get the total cooling rate Λ_rot = n_CO ∑ _JI_J.

The vibrational component is taken from Hollenbach & McKee (1979), considering collisions with H₂ whose rates, are (54)where E₁₀/k = 0.5E₂₀/k = 3080   K, T₃ = T^1/3, and k is the Boltzmann’s constant.

Cooling by the metals. The cooling process by the metals is included using the list of the metals adopted by Maio et al. (2007), namely CII, SiII, OI and FeII (see Appendix B in Maio et al. 2007). This source of cooling is particularly important for temperatures lower than 10²   K. The cooling is due to the fine structure level transitions of the ionized carbon 2P(3/2) → 2P(1/2), and similarly for the ionized silicates. There are two other species that take part in this process: neutral oxygen (nine transitions between levels , , , , and ) and ionized iron (five transitions between levels with i odd in the range  [1,9]  ∈ N). The cooling from each transition is additive and is given in erg/cm³/s by (55)where is the reaction rate for the hydrogen, and is the same but for the electrons, A_ij is the Einstein’s coefficient between the ith and jth levels, ΔE_ij is the energy difference between the levels, n_tot is the total number density of the species considered, and finally, n_H and n_e are the hydrogen and electron number density, respectively. To complete the equation we need to know (56)both for hydrogen and electrons; g_i and g_j are the level statistical weights, k_B is the Boltzmann’s constant, and T_gas is the gas temperature.

If the model is calculated at high redshift, the CMB pumping must be considered. This process is described by including the stimulated emission coefficient in the treatment of the cooling. The black-body radiation field of the CMB has a temperature that depends on the redshift as T_CMB(z) = (1 + z)T_CMB(0) K. However, as the models of the ISM refer to the current age, the temperature is kept constant and equal to the present day value (see below). In any case, owing to the low value of the present-day CMB temperature, this effect is neglected. To obtain the total cooling by the metals we use Λ_metals =  ∑ _iΛ_i where i indicates a generic level transition for any of the four metals we have included.

The cooling proposed by Maio et al. (2007) is improved implementing the cooling rates by Glover & Jappsen (2007); in particular, we use the de-excitation rates calculated for the collisions with ionized and molecular hydrogen when available. We also include the cooling from neutral C, Si (as in Glover & Jappsen 2007), Fe, and ionized O (as in Hollenbach & McKee 1989).

Total cooling. The cooling rate by all the above processes is additive and can be described by (57)where Λ_X is either Λ_C92 or Λ_SD93 as appropriate, and all the Λs are functions of the temperature. Both Λ_H2 and Λ_HD are the cooling functions of the two molecular species, and the terms in brackets are the cooling of the metals.

The overall rate of temperature change due to the heating and cooling is given by the following equation (58)where Γ is the heating source (if present), Λ the cooling source, n_i the number density (the sum is over all the elements), k_B the Boltzmann constant, and γ the adiabatic index defined in Merlin et al. (in prep.) as (59)where x is the fraction of the element indicated by the subscript. For an ideal gas of pure hydrogen this value is 5/3. If the gas is made only of molecular hydrogen we have γ = 7/5. In Eq. (59) we use a linear fit of the data proposed by Boley et al. (2007) for the adiabatic index of H₂ under log (T) ≤ 2.6 (see also erratum), considering a 3:1 ortho/para ratio mix.

The generic heating term Γ can be used to introduce heating phenomena like SNæ feedback, cosmic rays, and others. In the current version of robo, Γ is used as a free parameter that is kept constant during a single run. In any case, the Γ term can be specified by the user according to his scientific aims.

3. Code description

The code has been developed with IDL². Its user friendly features help the development of applications that are otherwise difficult to build in FORTRAN. The code is divided in self-explanatory procedures (routines) that are grouped in four classes (gas chemistry, dust, cooling, and general code behavior). The main routine calls all other subroutines that are needed to calculate the gas evolution. The first group of routines calculates the reaction rates and updates the density of each species.

Mass conservation. The total mass per unit volume of the ISM at time t is M(t) =  ∑ _in_i(t)   m_i, where n_i(t) is the current number density of the ith species, m_i its molecular or atomic mass, and the sum is over all the species. While the number densities can vary with time, the total mass must be conserved. In the numerical integration of the system of differential Eqs. (1), the conservation of the total mass is not always guaranteed, because it depends on the physical processes, the choice of initial parameters and the time step. Particular care is paid in our test computations to securing and checking run-by-run the conservation of the total mass, i.e. where M(t₀) is the initial total mass of the system per unit volume.

Time steps. Even with short time steps, it may happen that some species reach negative values in a single time step. This is because the differential of a species could be negative, and its absolute value higher than the relevant number density. This problem could be solved by forcing the species to be positive. This would cause difficulties with the mass conservation, and then in the subsequent time step the solver would again change the sign of the species abundance, and finally, the cooling would produce some nonphysical effects. In brief, a negative value of the abundance of the species would artificially turn the cooling into heating, which is clearly impossible. The obvious way out is to suitably choose the time step.

Numerical Solvers. In our system of differential Eqs. (1) there are 76 reactions to deal with (64 reactions plus 12 photochemical processes), so a fast and accurate solver is needed. We use the routine LSODE of IDL. LSODE uses adaptive numerical methods to advance the solution of a system of ordinary differential equations by one time step (Hindmarsh 1983; Petzold 1983), optimizing the process for user-defined time steps. In our case LSODE absolute tolerance is 10^-40, and the relative one is 10^-12.

Integrators. The number densities of all the particles (atomic species, molecules, dusty components, and free electrons) are governed by the balance between creation and destruction processes and, even more important, span very wide ranges of values.

Ranges of applicability. robo can be safely used in the following intervals for temperature, density, and metallicity: 10 ≤ T ≤ 10⁷ K, 10^-12 ≤ n ≤ 10³, and 10^-12 ≤ Z ≤ 10. The density range is somewhat limited towards the high value end: observational estimates of the density in molecular clouds can indeed reach higher values. Work is in progress to extend the density range. Furthermore, let us remind the reader that Z = 1 means  [Fe/H]  = 0, the solar case, and Z = 10 means  [Fe/H]  = 1; i.e., the ratio of the number densities (Fe/H) is equal to ten times solar. Finally, as not all the reaction rates cover the whole range of values, some extrapolations are required.

The aforementioned ranges of applicability are chosen in such a way as to guarantee the numerical stability of the system. Indeed, the core of the model is the system of differential equations that describe the chemical evolution of the gas according to the reactions listed in Tables 2 through 5. The stability of the system is measured by the conservation of the mass of the ISM elements. Therefore, we carefully checked whether robo satisfied this conditions for values of the parameters spanning over large volumes of the parameter space. The results show that this is always the case.

3.1. Free parameters

robo contains 47 parameters governing the physics, the mathematics, and the numerical procedure. Here we briefly comment on the most important ones.

• Ionized fraction of metals: it fixes the fraction of C, O, Si, and Fe at the firstionization level.

• Metallicity: the metallicity of the gas is defined as Z = 10 ^[Fe/H] . It is worth recalling here that Z = 1 is the solar case corresponding to  [Fe/H]  = 0 (see below for more details).

• Number densities: the number densities of the 28 elemental species and/or molecules, free electrons, and dusty components are all in cm^-3.

• Fraction of carbon dust: the percentage of carbonaceous grains. If it is equal to 1, all the dust is made by graphite grains and PAHs; if this parameter is equal to 0, all the dust is made by silicates. Intermediate values are also possible.

• Gas temperature: the temperature of the gas depends on the cooling and heating processes and changes during the gas evolution.

• Dust temperature: the temperature of the dust is kept constant during each run. It controls the formation of the molecules on the surface of the grains and also the accretion of the dust grains themselves.

• CMB temperature: the gas temperature cannot go below the temperature of the cosmic microwave background (CMB). At present the CMB temperature is kept constant during each run.

• Cosmic ray field: the field formed by the cosmic rays that populate the gas. It destroys or ionizes molecules like H₂. This field is expressed in s^-1, thus corresponding to the rate of ionization of the H₂ by the cosmic rays.

• Total integration time of a model: the total time of each run in s.

• Time step: the time step used in the models. A typical value is 10³ years. However, it can be changed by the routine LSODE.

• Finally, there is a number of flags activating or switching off the different sources of heating and cooling and the mechanisms of dust formation and destruction we have described in the previous sections.

4. Models and discussion

To validate robo we calculated two groups of models of gas evolution: the first one consists of 600 dust-free cases. Each case corresponds to a different combination of the parameters. The second group consists of 32 models of gas evolution, but the presence of dust grains is taken into account in them.

The aim here is to investigate the evolution of the same elemental species in the presence or absence of the dust. In addition to this, the large number of cases per group allow us to investigate the model response at varying the key parameters. Indeed, owing to the large number of parameters at our disposal, it is not possible to explore the whole parameter space.

Finally, it is worth noticing that the choice of the parameters for the two groups of models is somehow guided by robo mainly being designed to become a sort of ancillary tool for NB-TSPH codes. Therefore, in view of this, the two groups of models use the same parameter space adopted by EvoL.

4.1. Dust-free models

Each model is a simulation of a unitary volume of gas in the absence of dust. The results consist of 600 (see below) evolutionary models of the gas temperature and number densities of the 28 species under consideration plus the free electrons. As already emphasized, in each model special care is paid to secure the mass conservation and to avoid unphysical negative number densities. The models of this set are calculated neglecting the presence of any type of dust, and they are meant to check the internal consistency of robo.

4.1.1. Set up of the parameters

Each model gives the thermal and density history of the gas for an assigned set of parameters. The number of parameters and the values assigned to each of them determine the total number of models to calculate. This is given by , where N_m is the total number of models, p_i the number of different values for the ith parameter, and N the total number of the parameters, which in our case amounts to 47 (see Sect. 3.1).

For all the models, we must specify the initial value of the metallicity, the number densities n_H2, n_H⁺, n_e⁻, and the gas temperature. We adopt four values for the metallicity Z = 10[Fe/H]:  { 0,10^-6,10^-3,1 } ; five values for the H₂ number density:  { 10^-10,10^-6,10^-2,10^-1,1 }  cm^-3; three values for the H⁺ number density:  { 10^-10,10^-1,1 }  cm^-3; two values for the electron number density:  { 10^-10,10^-1 }  cm^-3; five values for the gas temperature:  { 10,10³,10⁴,10⁶,10⁸ }  K. More details on the metallicity are given below when we discuss the number densities of the metals. All the other parameters of the models have a constant value.

In all the models, the following quantities have a fixed initial value. They are n_H = 1.0 for the neutral hydrogen, n_H⁻ = 10^-9 for the hydride, for the molecular hydrogen and n_He = 0.08 for the helium. No helium ions are present at the beginning, so we have n_He⁺ = n_He ^+ +  = 0. The initial values of all the deuteroids are set to 10^-25: these are the deuterium atom, its ions (D⁺ and D⁻), the molecules HD and D₂, and the ion HD⁺.

Metals (i.e. C, O, Si, Fe and their ions C⁺, O⁺, Si⁺, and Fe⁺) are computed from the metallicity as follows³. We start from the general (60)where n_Fe and n_H are the number densities of iron and hydrogen, and the definition of metallicity we have adopted, namely Z = 10[Fe/H]. It must be emphasized that this notation for Z is different from the commonly used definition of metallicity, that is, Z = 1 − X − Y with X and Y indicating the abundances by mass of hydrogen and helium. In the usual meaning, Z is therefore the mass fraction of all the species heavier than helium. In our definition Z is simply related to the iron content [Fe/H] by the expression Z = 10 ^[Fe/H] . With this notation, Z can be higher than one: Z = 1 corresponds to the solar iron abundance. Considering now a generic metal indicated by X we can write (61)and (62)Therefore if we know the ratio , the total metallicity Z, and the number density n_H of hydrogen in the model, we can derive the number density of the generic metal X. By doing this we assume that for any metallicity the relative abundances of the elements follow the solar partition. In other words we do not consider here the possibility that the gas may have a different distribution of the heavy elements from that in the Sun with respect to the hydrogen (α-enhancement problem).

We use the same method to derive the number density of metal ions given by (63)where κ represents the ionized fraction n_X⁺/n_X of the metal X. In the present models we set κ = 0, so that all the metals start as neutral species. If the initial temperature is high enough, the hot gas can ionize the metals very early on in the course of the ISM evolution.

The temperature of the gas is a free parameter, but cannot be lower than the temperature of the CMB. For the present models we set T_CMB = 2.73  K, the present day measured value (Boggess et al. 1992; Fixsen 2009).

Each model of the ISM is followed during a total time of 3.15 × 10¹⁴  s, which approximately corresponds to 10⁷  yrs. This choice stems from the following considerations. Each model of the ISM is meant to represent the thermal-chemical history of a unit volume of ISM whose initial conditions have been established at a given arbitrary time and whose thermal-chemical evolution is followed over a time scale long enough for secular effects to develop but short enough to closely correspond to a sort of instantaneous picture of large-scale evolution of the whole system hosting the ISM unit volume. The initial physical conditions are fixed by a given set of parameters each of which can vary over wide ranges. One has to solve the network of equations for a time scale that is long enough to reveal the variations due to important phenomena such as star formation, cooling, and heating, but not too long to let the system depart from the instantaneous situation one is looking at. The value of 10⁷ yr resulted in a good compromise.

The initial value of the time step is 3.15 × 10¹⁰  s ≈ 10³  yrs. This time step determines the minimum number of steps required to cover the time spanned by a model. It means that each simulation needs at least 10⁴ iterations to be completed. The LSODE integrator may introduce shorter time steps depending on the complexity of the problem, so 10⁴ is the minimum number of required steps. This value for the time step seems to keep the system stable during the numerical integration.

While integrating, all the sources of cooling are kept active: metal cooling, H₂ cooling according to the prescriptions by Galli & Palla (1998) and Glover & Abel (2008), and finally, the cooling from deuterated hydrogen.

In general, the initial values of the number densities fall into three groups: (i) the elemental species with constant initial values, the same for all the models (namely H⁻, , He⁺ and He ^+ + ); (ii) the elemental species whose initial values are derived from other parameters (namely He, all the deuteroids and the metals, such as C and O, which depend on the choice the total metallicity Z), and finally, (iii) the elemental species with free initial conditions (namely H, H₂, H⁺, and e⁻).

Hydrogen group: H, H⁺, H⁻, H₂, and . The initial values of the species H⁻ and are n_H⁻ = 10^-9 cm^-3,  cm^-3 according to Prieto et al. (2008), while the three other hydrogen species have free initial values.

Deuterium group: D, D⁺, D⁻, D₂, HD, and HD⁺. The number densities of the deuteroids are calculated from their hydrogenoid counterparts. For the single atom species we have n_D = f_D   n_H, n_D⁺ = f_D   n_H⁺, and n_D⁻ = f_D   n_H⁻, where f_D = n_D/n_H. For the molecules, we can consider the ratio f_D as the probability of finding an atom of deuterium in a population of hydrogen-deuterium atoms. This assumption allows us to calculate the HD, D₂, and HD⁺ number densities as a joint probability. For HD and HD⁺ we have n_HD = f_D   n_H2 and , but for D₂ is as the probability of finding two deuterium atoms is . This is valid as long as f_D ≪ 1.

Helium group: He, He⁺, He ^+ + . The ratio n_He/n_H is set to 0.08, thus allowing the initial value of n_He to vary according to the initial value for n_H. The initial number densities of the species He⁺, He ^+ +  are both set equal to zero in all the models.

Metals group: C, C⁺, O, O⁺, Si, Si⁺, and Fe, Fe⁺. The Fe number density of the ISM is (64)where is the iron-hydrogen ratio for the Sun. To retrieve the number density of a given metal X we use n_X = n_Fe·f_X, where f_X is the metal-iron number density ratio in the Sun.

The list of the species whose initial number densities are kept constant in all the models of the ISM is given in Table 10. The values listed here are either fixed to a constant value or based upon the number density of one of the free hydrogen species n_H, n_H2 and n_H⁺ via the f_D factor. Since H⁻ and are constant, then also D⁻ and HD⁺ are fixed. Values are indicated as a(b) = a × 10^b.

The chemical composition of the ISM is typically primordial with a mass abundances of hydrogen X = 0.76, and helium Y = 0.24 and all metals Z ≃ 0. The helium-to-hydrogen number density ratio corresponding to this primordial by mass abundances is n_He/n_H ≃ 0.08. The adopted cosmological ratio for the deuterium is f_D = n_D/n_H ≃ 10^-5. With the aid of these numbers and the above prescriptions we get the number density ratios listed in Table 10 and the initial values of the number density in turn.

4.1.2. Results for dust-free models

In this section we discuss the results we have obtained for the dust-free models. Together with this we also consider a side group of models calculated with some specific assumptions (e.g. without metal cooling and others to be explained in the course of the presentation) to underline some effects and to better understand the physical processes taking place in the ISM.

Special attention is paid to check whether there are some unexpected effects or drawbacks in the models and in robo. We also describe how different physical quantities affect the overall behavior of the gas. This is achieved by varying a single parameter at a time and keeping constant all the others.

Temperature and H₂ density. In Fig. 4 we compare the temperature evolution for different values of the initial molecular hydrogen densities in two cases with very different metallicity, i.e. Z = 10^-10 and Z = 1. The important role played by the metal cooling is soon visible comparing the two panels. In the bottom panel (high-metallicity Z = 1), cooling is so strong that all the models simultaneously reach in a short time interval nearly the same final temperature of about 10  K. In the top panel (low-metallicity case) the cooling time scale is much longer than in the bottom panel because the contribution to cooling by the metals is negligible. Comparing each curve with the others in both panels, the effect of the molecular hydrogen is evident: higher H₂ densities imply steeper cooling curves and a shift toward lower temperatures. In the metal-free models, the temperature systematically shifts to lower values at increasing the density of the molecular hydrogen (top panel of Fig. 4). The effect is similar but enhanced in the case of metal-rich models (bottom panel of Fig. 4).

Metallicity. Figure 5 shows the evolution of the temperature (top panel) and H₂ number density (bottom panel) of models with different metallicity Z. Each curve is characterized by a different initial metallicity. The sources of UV radiation are at work in these models. As in the previous figure, the cooling strongly depends on the metallicity. The marked knee in the temperature-age relationship is caused by the cooling via the metallicity. Looking at the H₂-age relationship, we note that the formation of H₂ is favored when more free electrons are present (an effect that becomes clear when considering Figs. 5–7).

UV radiation. The UV radiation plays a key role in the H₂ density evolution. To cast light on the issue, we calculated models with no UV sources. They are displayed in Fig. 6. The top panel shows the temperature - age relationship, whereas the bottom panel shows the H₂ number density - age relationship. Comparing the bottom panel of Fig. 5 to that of Fig. 6, it is clear the effect of the UV radiation. When the UV radiation is present (Fig. 5), H₂ is destroyed very early on, whereas when the UV radiation is absent (Fig. 6) the number density of H₂ remains nearly constant during the whole evolution. We also note that the effect of the cooling is stronger in the first case. This happens because the UV flux increases the density of free electrons, and since they act as colliders, the eventual higher number of collisions determines a marked improvement in the cooling process. This effect of the UV flux on the number density of free electrons is displayed in Fig. 7, where we show the time evolution of this quantity at varying the metallicity and switching on (top panel) and off (bottom panel) the UV radiation.

Ionized versus neutral metals. The difference in the cooling between Fig. 5 (top panel) and Fig. 6 (top panel) could be attributed to the different efficiencies of the cooling by ionized and neutral metals. Indeed, it is worth noticing that, in order to make the effects of the metals evident, the initial number density of H₂ in Figs. 5 and 6 is set to 10^-6. The UV radiation affects the ionization state of atoms and molecules, thus varying the partition between ionized and neutral metals. To clarify this issue, we first compare models in the presence of UV flux in Fig. 8, but the metal cooling by ionized species (top panel) or neutral species (bottom panel) are alternatively switched off. The initial number density of H₂ is always set to 10^-6. For each metallicity, the values shown in Fig. 5 (top panel), where both neutral and ionized metals are included, are roughly the sum of the values presented in Fig. 8 (both panels). In the top panel of Fig. 8, the cooling by ionized metals is switched off so that we would expect a behavior similar to that of Fig. 6 (top panel), where the UV radiation is switched off and the neutral metals are favored. It must be pointed out, however, that when the UV is switched off (Fig. 6) most metals are neutral. Even if the ionized metals are switched off (Fig. 8), the fraction of neutral metals will be much lower. Indeed, the shape of the curves displayed in Fig. 8 (top panel) remains similar to those shown in Fig. 5 (top panel), and the long cooling time scales of Fig. 6 cannot be reproduced.

For the same effect of the UV flux on metals, in the top panel of Fig. 8 (where the UV flux is on and the ionized metals are off) cooling does not depend on the metallicity. The UV radiation easily ionizes metals like silicon, carbon, and iron: UV radiation below 13.6 eV from various astrophysical sources generates a UV background that ionizes atoms with a first ionization potential lower than 13.6 eV (Maio et al. 2007). Consequently, with the UV flux switched on, the total amount of neutral metals becomes negligible. If the cooling by ionized metals is switched off and the contribution to the total cooling by neutral metals is negligible, it follows that the metals almost do not contribute to cooling. When the ionized metals are put back and the negligible amount of neutral metals is dropped (Fig. 8 – bottom panel), we recover the situation of Fig. 5 as expected.

Now the UV flux is switched off as in Fig. 6 and the metal cooling by ionized species (top panel) or neutral species (bottom panel) is alternatively removed. As we can see in Fig. 9, with no cooling by ionized metals and only with the contribution by neutral metals, we are able to reproduce the behavior observed in Fig. 6 (top panel). The case without neutral metals and only with cooling by ionized metals is not shown because it would lead to an almost constant temperature in all the cases, consistent with the picture we have just described. In brief, in this case, both H₂ (the input number density of H₂ is set to a lower value, thus contributing less to the cooling) and the ionized metals (without UV flux, as expected, most of the metals are neutral) little contribute to the cooling process. We can finally conclude that the partition between ionized and neutral metals fully explains the difference between Figs. 5 and 6.

Metallicity and H₂ vs. Cooling. In Fig. 10 are plotted four different models with different combinations of the metallicity (very high or very low) and cooling by metals (included or switched off)⁴. Furthermore, to examine the effect of the H₂ cooling, we set the input number density of H₂ equal to 10^-1 cm^-3 (see Sect. 4.1.1 for more details on the selected input values for H₂). This plot highlights how the cooling by metals and H₂ affects the temperature and thus the physical status in the gas. The steep decrease of the temperature from the beginning of the simulations and common to all four models is due to the strong H₂ cooling. As expected, the model with high-metallicity and cooling by metals (top panel of Fig. 10, dashed line) is the one experiencing the strongest cooling. The two models with very low-metallicity have a similar behavior independent of the presence or absence of the cooling by metals. The last case, with high-metallicity and no cooling by metals, is the one with the lowest temperature decrease. The reason for it is explained in Fig. 11, where the evolution of H₂ for the four cases under consideration is shown.

The model with high-metallicity and no cooling by metals has the lowest H₂ density compared to the other three and consequently has low cooling. From Fig. 11 we also see that the model with high-metallicity and cooling by metals (dashed line) has the highest H₂ density. This can explain the results in the top panel of Fig. 10, in the sense that the cooling process here could be mainly due to the contribution coming from the molecular hydrogen. The bottom panel of Fig. 10 shows the same models with no cooling by H₂. In this case only the model with cooling by metals included and high-metallicity undergo a significant cooling. Indeed, if the effect of H₂ is neglected, we need metals in significant amounts to have strong cooling.

4.2. Models with dust: parameter set up

In this section we consider models in which the effect of the dust grains is taken into account. We adopt here the same values for the parameters as in the dust-free models.

The minimum initial metallicity is Z = 10^-6, and the initial number density of the molecular hydrogen is n_H2 = 10^-6  cm^-3. The number densities of electrons and H⁺ are free parameters; they are n_H⁺ =  { 10^-10,10^-1 }  and n_e⁻ =  { 10^-10,10^-1 } , both in units of cm^-3 as usual. The initial temperature is set to T = 10⁴  K.

The only difference with respect to the previous models is the presence of dust. We adopt four values for the number density of dust grains, namely n_dust =  { 0,10^-3,10^-2,10^-1 }  cm^-3. In these models the composition of the dust mixture is 50% carbonaceous grains and 50% silicates.

We also calculate models with or without dust sputtering by shocks, in order to describe the behavior of a turbulent gas particle with and without the grain depletion due to the shocks.

All the other parameters remain the same as in the dust-free models, such as the number densities of different elements and the cooling processes. In all these models thermal sputtering and dust formation are active, so the dust properties are let change during the evolution of the interstellar medium.

4.2.1. Results for models with dust

The series of plots going from Fig. 12 through Fig. 14 show the temporal evolution of three important quantities of the models, namely the gas temperature, the number densities of dust, and molecular hydrogen for different amounts of initial dust and different initial metallicities. Clearly there is a tight relationship between the initial amount of dust, the temporal behavior of temperature, and the number density of H₂. First, at a given low-metallicity and by increasing the dust fraction, the temperature decreases earlier and faster (top panel of Fig. 12), whereas if the metallicity is high there is no remarkable effect of the increased dust content (bottom panel of Fig. 12). Also, for the low-metallicity, the trend is anticipated as the dust content increases.

Looking at the temporal evolution of the number density of H₂ shown in the panels of Fig. 13 we note that, in coincidence with the temperature fall off and subsequent gentle decrease, the H₂ number density first decreases and then increases, forming a local minimum. This happens because the dust drives the formation of H₂: models with more dust suffer more cooling and hence form more H₂. Only if the metallicity is very high, does the cooling effect of dust grains lose importance as shown by the bottom panel of Fig. 12. To conclude, the curves displayed in the top and bottom panels of Fig. 13 nearly have the same shape, but a different vertical offset. The position and depth of the minima depends on the temperature variation, and the vertical offset clearly depends on the amount of dust present in the gas. This means that amount of dust and H₂ are closely related.

Finally, we calculate and present six models (three for each initial amount of the dust density) that include now the dust destruction by shock sputtering and vaporization. Their temporal evolution is limited to the first 10⁵ years only. The results are shown in Fig. 15 (the gas temperature) and Fig. 16 (the dust grains number density as a function of the size of dust grains). Looking at Fig. 15, the temporal evolution of the gas temperature is the same as in the previous case with the sole thermal sputtering at work. In contrast, the dust number density undergoes big changes as shown in Fig. 16, because the shock sputtering is very efficient. Figure 16 shows how the dust distribution changes with the time for the n_dust = 10^-3 cm^-3 case. After plotting the data of Fig. 16, the dust grains have been grouped in bins according to their size⁵. The lines represent different distributions at different ages, namely 10, 10², 10³, and 10⁴ years from top to bottom. As expected, the shock first destroys the large size grains and then the small ones.

5. Including ROBO results in an NB-TSPH code

In this section we briefly describe how robo is planned to be used as an ancillary tool for the NB-TSPH code EvoL to calculate the thermal and chemical properties of the gas particles. In the following we present some preliminary results. The method and its results will be widely discussed in a forthcoming paper (Grassi et al. 2011).

We already mentioned that two different techniques can be used. The first one is a real-time method, in which the elemental abundances are calculated solving the Cauchy problem of the chemical network. It means that for every time step and for every particle, EvoL needs to perform this computation using robo as a routine dedicated to this purpose. The advantage is that elemental abundances, thermal properties of the gas particles, and cooling can be calculated with high precision. However, this would require large computational resources, and unfortunately, the more elements it tracks, the more CPU time is needed. Furthermore, the network complexity and additional physical processes (e.g. detailed cooling) would increase the computational cost of it.

To avoid these problems an NB-TSPH code can use the s- called “grid” method. It consists of calculating a large number of models exploring different scenarios produced with different sets of the input parameters. The results are then tabulated. When a particle in the cosmological or galaxy simulation needs to know its evolution after a time step, the particle looks for the closest parameters set in the grid, and interpolate the evolved parameters. This a lightweight method, but the main drawback is that the number of free parameters cannot be too large, mainly because interpolating multidimensional grid can be difficult and inaccurate (depending on the grid coarseness).

Our method is similar to the latter one. First of all, we select a set of parameters. Second, we run robo for each parameter combination over the total integration time we have chosen (3.15 × 10¹⁴ s). The total number of simulations depends on the number of parameters and on the multiplicity of each parameter (see Sect. 4.1.1). This part can easily be achieved using for example a cluster of computers.

Now we have a huge number of models characterized by a large set of parameters. With the standard grid method, interpolating the requested data would be complicated. To solve this problem we use a powerful tool: the ANNs. There are different kinds of ANN architectures depending on the task it needs to perform. We choose the so called back-propagation algorithm (Rumelhart et al. 1986). It is one of the most used and versatile methods for solving a wide range of problems. This mathematical method consists of an algorithm that behaves like the gas model itself. Our ANN learns how the gas evolves using the models produced with robo. We first perform a training stage, in which a subset of randomly chosen models is shown to the ANN. After some iterations the ANN converges to a stable state. At this point we look at the complementary of the subset of models: its aim is to calculate the ANN error on unknown data. If the error is small enough, we can assert that the ANN behaves like the gas model. When this happens, in EvoL robo can be safely replaced by the ANN.

The great advantage is that the ANN can give the requested solution with ease. It is very light, since it consists of three small matrices of so-called synaptic weights. Each particle needs only to perform a matrix product, which is (in this scenario) a fast algebraic operation. In this way retrieving chemical evolution for each particle is fast. To illustrate the ability of the ANN to provide good representations of the physical state of the ISM in Fig. 17, we show the comparison between the theoretical data calculated by robo for the free electron number density and fed to the ANN (squares) and the same data retrieved by the ANN (crosses). The agreement is remarkably good. The same technique can be applied to any other variable of interest here.

6. Conclusions

We have presented a model of the ISM that provides a detailed description of the gas chemistry and evolution, the formation and destruction of dust grains of different types, and finally, a thorough description of the cooling process over a wide range of physical parameters and initial conditions. The way the model is conceived corresponds to an instantaneous picture of the physical state of an elementary volume of the ISM characterized by a set of physical parameters assumed here as the initial conditions of a given volume element. Under the action of the ISM models, the initial physical state evolves on a secular time scale. This provides us a sort of vector field telling how a given physical state will evolve (how much and in which direction in the multidimensional space of the physical conditions). The integral of the elementary volumes over the underlying evolutionary path of the grand physical quantities like density and temperature (all of these functions of space and time) of the host system (a galaxy or a cosmological simulation) will give us the detailed evolution of the ISM. This is the big advantage offered by the model, securing it a wide range of applicability.

We have presented here the results for dust-free and dust-rich ISM at varying the key parameters. The first group of models for a dust-free medium is meant to understand how the ISM behaves in the absence of dust grains. These models highlight the importance of the different kinds of cooling that are dominant in different kinds of environments. In particular, we call attention to the role of metals and free electrons in driving the physical behavior of the ISM via their effect on the gas cooling during its evolution. The dust-rich ISM allows us to understand how the ISM responds to the presence of the dust. In particular, we analyzed the temperature variations caused by the presence of dust in different amounts. We have also explored how the creation and destruction of the dust grains (the latter induced by shock and thermal sputtering) affects the evolution of the ISM.

The ISM model and companion code were created as auxiliary tools for NB-TSPH simulations in the context of galaxy cosmological simulations of the Universe and models of galaxy formation, structure, and evolution. Our specific aim is to give a more accurate description of the gas component in EvoL or in similar codes in literature. Finally, robo is also designed to run in small and middle-size computers. Thanks to robo, detailed gas physics can be inserted in NB-TSPH simulations at low computational costs.

To include the results of robo in our NB-TSPH code, we plan to use the ANNs that are more accurate, faster, and easier to implement than the standard fits on multidimensional grids. A complete account of this will be made public soon (Grassi et al. 2011).

Future implementations of robo are planned, among which we mention the inclusion of the photo-ionization by single stellar populations of different age and chemical composition in the chemical network and a better determination of the grains temperature that is tightly related to the local stellar radiation field.

---

Acknowledgments

The authors would like to thank the referee, Dr. Simon Glover, for useful and constructive remarks that greatly helped improve the first version of the manuscript. T. Grassi is grateful to Dr. F. Combes for the kind hospitality at the Observatoire de Paris – LERMA under EARA grants, where part of the work was developed and for the many stimulating discussions. P. Krstic acknowledges support from the US DOE Office of Fusion Sciences through ORNL, under contract No. DE-AC05-00OR22725 with UT-Battelle, LLC.

References

All Tables

All Figures