Modelling dust evolution in a multiphase, inhomogeneous ISM

Modelling dust evolution in galaxies with a multiphase, inhomogeneous ISM

Svitlana Zhukovska11affiliation: Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany , Clare Dobbs22affiliation: University of Exeter, Stocker Road, Exeter EX4 4QL, United Kingdom , Edward B. Jenkins33affiliation: Princeton University Observatory, Princeton, NJ 08544-1001, USA , Ralf Klessen44affiliation: Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany

We develop a model of dust evolution in a multiphase, inhomogeneous ISM including dust growth and destruction processes. The physical conditions for grain evolution are taken from hydrodynamical simulations of giant molecular clouds in a Milky Way-like spiral galaxy. We improve the treatment of dust growth by accretion in the ISM to investigate the role of the temperature-dependent sticking coefficient and ion-grain interactions. From detailed observational data on the gas-phase Si abundances measured in the local Galaxy, we derive a relation between the average and the local gas density which we use as a critical constraint for the models. This relation requires a sticking coefficient that decreases with the gas temperature. The relation predicted by the models reproduces the slope of of the observed relation in cold clouds. It is steeper than that for the warm medium and is explained by the dust growth. We find that growth occurs in the cold medium for all adopted values of the minimum grain size from 1 to 5 nm. For the classical cut-off of  nm, the Coulomb repulsion results in slower accretion and higher than the observed values. For  nm, the Coulomb interactions enhance the growth rate, steepen the slope of relation and provide a better match to observations. The rates of dust re-formation in the ISM by far exceed the rates of dust production by stellar sources. After the initial 140 Myr, the cycle of matter in and out of dust reaches a steady state, in which the dust growth balances the destruction on a similar timescale of 350 Myr.

1 Introduction

A fraction of metals in the interstellar medium (ISM) of galaxies in both the local and high-redshift Universe resides in tiny refractory solid particles or dust grains (Dorschner & Henning, 1995). Interstellar dust constitutes less than 1% of the ISM by mass, but it has manifold impact on the physics and chemistry of the ISM. Because dust locks some elements away from the gas phase, it reduces abundances of important gas coolants such as C and Si (Bekki, 2015b; McKinnon et al., 2016). One of the most important roles of interstellar grains is that they facilitate the formation of molecular hydrogen (H) on their surfaces (Hollenbach & Salpeter, 1971). The H molecule is the main component of molecular clouds, which are the cradle of star formation in most of the Universe (Klessen & Glover, 2016). Dust absorbs ultraviolet emission from young massive stars and re-emits it in the infrared, therefore the spectral energy distribution from dust is one of the primary indicators of star formation (Calzetti, 2013). Moreover, continuum emission from interstellar dust is a commonly used tracer of cold gas in galaxies (e.g., Santini et al., 2014).

Despite the utmost importance of interstellar dust for the ISM evolution, high-resolution numerical simulations of the ISM do not yet follow dust evolution. It is commonly assumed that dust abundance is scaled with metallicity and the scaling factor, grain size distribution and chemical composition have the average characteristics of dust in the local Galaxy (e.g. Walch et al., 2015, and references therein). Galactic-scale chemodynamical models consider the dependence of galactic evolution on dust by including H molecule formation on grain surfaces, but they follow the evolution of metallicity only and assume a fixed dust-to-metals ratio (e.g. Christensen et al., 2012; Forbes et al., 2016, and references therein). Hu et al. (2016) recently investigated how the choice of dust-to-metal ratio value affects the evolution of dwarf galaxies through photoelectric heating process, keeping other dust properties fixed.

There is a strong observational evidence that cosmic dust properties are not universal. Rémy-Ruyer et al. (2014) showed that there is large scatter in the observed dust-to-gas ratio vs. metallicity relation in local galaxies. Moreover, metal-poor dwarf galaxies tend to have lower dust-to-metal ratio compared to normal spiral galaxies. Recently, IR and far-IR surveys of dust emission opened a new venue to probe the spatial variations of dust properties on the small scale. Roman-Duval et al. (2014) analysed the IR emission maps of the Magellanic Clouds combined with various gas surveys and found large variations in the gas-to-dust ratio that are strongly correlated with the dust surface density distribution. Recently the far-IR emission survey with the Planck Satellite discovered similar variations of the gas-to-dust ratio both from cloud to cloud and within regions of individual clouds in the Milky Way (Reach et al., 2015). An estimate of the true values of the dust-to-gas ratio at high surface densities is a challenging task, since different effects may cause variations of the apparent dust-to-gas ratio measured from observations: possible underestimation of the molecular gas mass because of the presence of CO-dark molecular gas, variations of dust emissivity in far-IR due caused by grain coagulation in dense clouds, and actual increase of the dust-to-gas ratio due to gas-grain interactions (see Roman-Duval et al., 2014, for discussion). Variations of dust properties from one line of sight to another in the local galaxy are also supported by observations of extinction curves (Fitzpatrick & Massa, 2007), dust opacities (Roy et al., 2013), spectral characteristics (Dartois et al., 2004), and interstellar element gas-phase abundances (Savage & Sembach, 1996a; Jenkins, 2009). Thus, not only the dust-to-metal ratio, but also the composition and the size distribution of interstellar grains can differ from the standard. It is thus crucial to consider dust evolution in the context of modelling physics and chemistry of the ISM and using dust as a tracer to study the interstellar gas.

Evolution of interstellar dust is incorporated in one-zone models of the chemical evolution of galaxies, which consider dust properties averaged over the entire galaxy or over the vertical direction within concentric rings in spiral galaxies (Dwek & Scalo, 1980; Dwek, 1998; Hirashita, 1999; Zhukovska et al., 2008; Zhukovska & Gail, 2009; Zhukovska & Henning, 2013; Calura et al., 2008; Mattsson & Andersen, 2012; Asano et al., 2013b; de Bennassuti et al., 2014; Hirashita et al., 2015). These models demonstrate that the galactic dust content is mainly determined by the balance between dust destruction in interstellar shocks and dust production by evolved stars and growth by accretion of gas-phase species in the ISM. Given the relatively short timescales of dust destruction in the ISM of a few  yr (Draine & Salpeter, 1979; Seab & Shull, 1983; Jones et al., 1994, 1996; Serra Díaz-Cano & Jones, 2008; Bocchio et al., 2014), the dust growth in the ISM becomes the dominant dust source in galaxies on the timescale from a few  Myr to a few Gyr (Zhukovska, 2014). An accurate treatment of dust growth is therefore important for modelling of interstellar grains. Although simple dust evolution models can describe the average dust properties in these present-day galaxies, they do not consider the complex structure of the ISM and its impact on dust growth rate and are not able to describe the observed variations of dust properties across the ISM phases.

Recently, there have been a few attempts to take in to account the evolution of dust in numerical hydrodynamic galaxy simulations. Bekki (2013) included the time evolution of interstellar dust abundances in the chemodynamical model of disk galaxies and made a step forward by coupling it to the galactic evolution via H formation on grain surfaces. The hydrodynamics of interstellar gas is modelled by the Smoothed Particle Hydrodynamics (SPH) method. With assumption that dust is coupled to gas, the model considers variations of the dust-to-metal fraction in each gas particle. A disadvantage of the model of Bekki (2013) is the current implementation of dust processing in the ISM in which growth and destruction processes do not depend on the local conditions and occur on constant timescales. In following works, the dependence of the growth timescale on local conditions is included by scaling it with temperature and density of the particles (Yozin & Bekki, 2014; Bekki, 2015a). A new live dust particle model has been presented by Bekki (2015b), who decoupled gas and dust particles and implemented additional gas-grain interactions and radiation pressure on grains. Recently, dust evolution has been incorporated in the moving-mesh simulation code AREPO and used in zoom-in cosmological simulations of Milky Way-sized galaxies (McKinnon et al., 2016). Their assumptions for the dust model are similar to those outlined in Bekki (2015a), with the exception of the destruction timescale, which they related to the local SN rate. While their simulations agree with a number of observables, they over-predict the dust-to-metal ratio in the circumgalactic medium.

The hydrodynamic galactic simulations with evolving dust clearly demonstrate that grains influence evolution of galaxies, albeit not as strong as SN or active galactic nuclei feedback. However, the existing simulations should be improved in several ways. For example, they do not treat properly the dependence of the growth timecale on metallicity. The timescale of dust growth in the ISM is inversely proportional to the metallicity, resulting in the existence of the critical metallicity for the dust growth (e.g., Zhukovska, 2008; Zhukovska & Gail, 2009; Asano et al., 2013a). While some simulations include a dependence of the growth timescale on the local temperature and density (Bekki, 2015a, b; McKinnon et al., 2016), they assume a fixed sticking coefficient over the entire simulation volume. With this assumption, grains can grow by accretion even in the warm or hot gas, If they stay sufficiently long time there (1 Gyr). However, sticking of the impinging species to the grain surfaces in this case is unlikely because of their high thermal energies (D’Hendecourt et al., 1985). Moreover, the resolution of cosmological simulations () is not sufficient to investigate dense regions where dust growth is expected or to address the observed variations in dust abundances across the ISM phases in the local Galaxy.

In the current study, we describe a new model of dust evolution in the inhomogeneous ISM including dust destruction by SN shocks and dust growth in the ISM. We apply the model to study the evolution of the three-dimensional (3D) dust distribution in the local Galaxy using histories of physical conditions from hydrodynamic simulations of the lifecycle of giant molecular clouds (GMCs) (Dobbs & Pringle, 2013). These simulations have sufficiently high resolution compared to the previous works to investigate the changes experienced by grains as they cycle between molecular clouds and ambient ISM. We substantially improve the treatment of dust growth in the ISM in two ways: by including a temperature dependent sticking coefficient and Coulomb interactions in calculations of the growth timescale. An advantage of our post-processing approach is that we can run multiple models to investigate how different model assumptions affect the resulting distribution of dust. Our main goal is to use the analysis of the large amount of element depletion data measured in the local Galaxy (Jenkins, 2009) to constrain the uncertainties in microphysics of the growth process.

The paper is organised as follows. We describe the initial conditions and main assumptions of these simulations in Sect. 2. Histories of physical conditions of gas parcels from the hydrodynamic simulations provide input for the dust evolution model. Section 3 presents the formulation of the dust evolution model and our choice of model parameters. Section 4 introduces the data on element depletions in the local Milky Way providing observational constraints for the model. Results of the model calculations of dust evolution as traced by Si element depletions are presented in Sect. 5, where we discuss the timescales of dust destruction and formation, the rates of dust growth in the ISM as a function of ambient density and the distribution of the element depletion values. Sect. 5.4 compares the theoretical and observed trends of Si depletion with density. Our conclusions are presented in Sect. 6.

2 Simulations of lifecycle of giant molecular clouds

In this paper, we utilise a hydrodynamic simulation of the ISM in a spiral Milky Way-like galaxy described in Dobbs & Pringle (2013). The simulation is performed using the SPH code sphNG (Benz et al., 1990; Bate et al., 1995; Price & Monaghan, 2007). The gas in the simulation is subject to a galactic potential. We use a logarithmic galactic potential (Binney & Tremaine, 1987) which provides a flat rotation curve, and a two armed spiral potential from Cox & Gomez (2002). The spiral potential is fixed, the spiral arms rotating with a constant pattern speed of , and represents a perturbation of a few per cent compared to the logarithmic potential. GMCs form within the arms through a mixture of gravitational and thermal instabilities in the gas, and cloud-cloud collisions (Dobbs, 2008). There is a total of 8 million gas particles, and each one has a mass . The gas is situated within a radius of kpc, and the gas settles in approximate vertical equilibrium after Myr (see Dobbs et al. 2011). The average surface density of the gas is 8.

Cooling and heating of the ISM are incorporated according to Glover & Mac Low (2007), with temperatures representative of a multiphase medium, spanning a range K to slightly above K. The cooling processes include collisional cooling (mainly from C, O and Si, but H cooling is also included) and the main heating process is photoelectric heating. The cooling is metallicity dependent, with the simulation assuming solar metallicity. Although in reality the metallicity would vary with radius, solar metallicity is consistent with the metallicity used in the dust post-processing which reflects the properties of the ISM at the solar radius. Self gravity of the gas is also taken into account.

Stellar feedback is included in the form of supernovae feedback: once gas reaches a density of 500 and is gravitationally collapsing. Thermal and kinetic energy are added to particles within a radius of about 15 pc from the densest particle, which is based on the smoothing length at these densities. The energy is inserted according to a snowplough solution, for the exact equations see the appendix of Dobbs et al. (2011). The amount of energy added is calculated as


where is the star formation efficiency, chosen to be 0.05, is the mass of molecular hydrogen within , erg is the energy of 1 supernova. We adopt a Salpeter IMF such that 1 supernova occurs for every 160 of stars formed.

The hydrodynamic simulation then provides input for the dust evolution model such as the density, temperature and masses of the SPH particles at each stored snapshot with cadence of 1 Myr. The particles which have undergone a recent supernova feedback event can also be identified. In particular, for each supernova event, the dust evolution model uses the total mass of gas where feedback is injected,


where refers to the th feedback event, is the number of particles in the th feedback event (typically around 20). All other assumptions about how the dust evolves in the ISM, and how it is destroyed by feedback, are applied in a post process step according to the dust evolution model described in the next Section.

3 Dust evolution model

In evolved metal-rich systems such as the present-day Milky Way, the timescales of enrichment of the ISM with chemical elements from stars are significantly longer than the timescales of mass transfer between the ISM phases. This allows us to make a number of simplifying assumptions and focus on the dominant sinks and sources of dust. Firstly, we fix the total (dust+gas) element abundances, and, secondly, we consider the dust production by growth in the ISM and neglect dust input from stars (O’Donnell & Mathis, 1997). Additionally, we only include destruction of dust in the interstellar shocks and neglect destruction by star formation since its timescale of Gyr is much longer than the current estimates for the dust lifetimes against destruction in shocks. We check the validity of these two latter assumptions by comparing the corresponding dust formation rates in Sect. 5.

We model the evolution of dust grains in the ISM by solving numerical differential equations for the degree of condensation, i.e. the fraction of the key element in dust. The key element is usually a constituent of the considered dust species which determines the reaction rate for dust growth. It is usually the species with the lowest gas-phase abundance. The is related to the logarithmic depletion as The degree of condensation can be converted to the dust mass in a gas parcel as


where and are the key element abundance and its atomic weight, and is the mass fraction of the key element in the considered solid, is the mass of an SPH particle, and is the hydrogen mass fraction in the gas.

Dust evolution is governed by growth and destruction processes, so that the change of the degree of condensation is

Physical quantity Value
Key element Si
Solid material density () 3.13
Element abundance ()
Initial depletion
Threshold for destruction in diffuse gas () 1
Timescale for destruction in diffuse gas (Gyr) 0.1
Table 1: Basic data used in model calculations

3.1 Growth in the ISM

The mechanism of grain growth in the ISM is not fully understood. It has been suggested that accretion of gas-phase species on silicate and carbonaceous grains occurs selectively to keep them as two distinct populations (see discussion in Draine, 2009). We therefore model the evolution of silicate dust independently of carbon dust. Note that, in the present work, we investigate the accretion of refractory elements on the grains which proceeds prior to the accretion of complex ice mantles and coagulation in dense cores of molecular clouds (e.g. Joseph et al., 1986; Weingartner & Draine, 1999; Voshchinnikov & Henning, 2010).

One-zone dust evolution models usually assume that the dust growth occurs in molecular clouds based on their higher densities and consequently shorter collision timescales of gas-phase species with grain surfaces compared to the diffuse medium (Dwek, 1998; Zhukovska et al., 2008). We relax this assumption to include possible dust growth in the cold neutral medium (CNM) defined accordingly to the gas temperature range of 30 – 300 K (e.g. Mihalas & Binney, 1981).

The change of the degree of condensation owing to the dust growth in the ISM by collisions of gas phase species with the grain surfaces is (Zhukovska et al., 2008)


Here is the growth timescale of a given dust species


where is the sticking coefficient that expresses the probability of sticking of the growth species in collision with the surface, is the density of the solid. The thermal velocity and the number density of the gas are determined by the local conditions recorded in the particle histories. The average grain radius is




is the th moment of the grain size distribution defined so that is the number of grains with radii from to .

In the diffuse ISM, most of the key species are singly ionised, therefore electrostatic interactions can enhance or reduce accretion rates, depending on the grain charge. To include this effect of Coulomb interactions in our model, we adopt the enhancement factors from work of Weingartner & Draine (1999) calculated for the standard interstellar UV radiation field, collisions with singly ionsed gas species and the grain charge distribution described in Weingartner & Draine (2001). We include enhancement effects in the CNM using a modified average grain radius in Eq. (6)


where the electrostatic factor accounts for the change in the cross section of an interaction between ion and grain. For neutral particles in MCs, .

The timescale of dust re-formation in the ISM is defined as


which for the simulations with equal mass particles can be written using summations over all particles


We assume that dust mass grows in the quiescent ISM and exclude the gas particles which undergo SN feedback event described in the next section.

3.2 Dust destruction in SN shocks

Sputtering of grains in high velocity shocks () is the main mechanism of dust destruction in the ISM. The large scale simulations used in this work do not resolve the scale on which gas is shocked to high velocities, therefore our approach to calculations of dust destruction utilises the existing SN energy feedback scheme described in Sect. 2.

In each feedback event, the energy from SNe is equally distributed between the gas particles with the total mass given by Eq. (1). In the same feedback event, dust is completely destroyed in a mass of gas which for the homogeneous ambient medium is defined as (McKee, 1989)


where and are the initial and final velocities of the SNR expanding into an ambient medium of density , respectively, is the mass of gas swept up by a shock with velocity in the range of , is the degree of dust destruction in a gas parcel shocked to the velocity . To calculate for each SN remnant, we use Eq. (12) together with an analytical expression for given by Dwek et al. (2007) which describes the adiabatic expansion and pressure-driven snow-plough stages of the SNR evolution.

With the resolution of our hydrodynamic simulations, is lower than . We therefore make two assumptions which allow us to simplify the treatment of dust destruction. We approximate the total dust mass destroyed in the th event by and assume the dust content is reduced equally in all affected gas particles by a fraction


The change in the condensation fraction in each gas particle in the th feedback event is . The change of condensation degree in a gas parcel which experienced the th feedback event is


where is the time step following the SN explosion. By using Eq. (3) and summing over all particles affected by the total feedback events in GMCs, we derive the total dust mass destruction rate


where is the total number of feedback events that occurred between two subsequent snapshots.

The current implementation of stellar feedback in our simulations neglects such forms of stellar feedback as stellar winds and radiation, prostellar jets and SNe exploding in low-density regions. The impact of different mechanisms of stellar feedback on the dynamics of the ISM is discussed by Walch et al. (2015). Among different processes, blast waves in the low-density gas are most important for dust destruction (McKee, 1989; Draine, 1990). Observationally it has been determined that 20–25% of OB stars belong to the field stars and will likely explode in the diffuse ISM (Oey & Lamb, 2012). Without account of these SNe we may underestimate the overall dust destruction in the ISM. To investigate their effect on the outcome of simulations, we also run models with an additional term in the total destruction rate for all particles with density below adopted as an upper limit for the diffuse ISM


where is the timescale to destroy all dust in the diffuse medium. This destruction rate is added on each time step only to those particles which are not affected by current stellar feedback events in GMCs. Note that in order to derive the timescale to destroy all grains in the ISM by this process has to be divided by the mass fraction of diffuse gas.

Figure 1: Sticking coefficient for H chemisorption on the silicate surface derived by Chaabouni et al. (2012) (solid line) and that for physisorption derived by Leitch-Devlin & Williams (1985) for dust temperatures 10, 20 and 50 K.
(nm) (nm)
normalaaThe average grain radius derived from Eq. (7). enhancedbb derived from Eq. (9) for the enhancement factor for silicate grains in the CNM.
1 15.8 0.8
3 27.4 7.7
5 35.4 48.2
Table 2: Average grain radius adopted in model calculations for the selected minimum sizes of the MRN grain size distribution

3.3 Model parameters

To solve equations for the condensation degree numerically, we need to make assumptions about the main properties of dust. For the total (gas+dust) chemical composition, we adopt solar metallicity and element abundance ratios from Lodders et al. (2009). In this work, we focus on growth of silicate dust, assuming that Si is the key species in silicate growth. For the adopted olivine-type composition (MgFeSiO) of silicates, the mass fraction of Si in dust is . All physical quantities used in the calculations are summarised in Table 1.

For the initial dust abundance, we take a constant value for all gas particles (). It is on the lower limit of the general depletion level in the WNM derived by Savage & Sembach (1996b) and somewhat lower than the value suggested by the depletion data from work by Jenkins (2009) (see Sect. 4). This value is consistent with the Si abundance required by models for extinction, emission, and polarisation of light by dust in the diffuse ISM in the solar neighbourhood (Siebenmorgen et al., 2014). We intentionally take low initial dust abundances for all particles with higher density, to explore how the system reaches the balance between destruction and production processes.

Efficiencies of grain destruction for different dust species are obtained from extensive numerical calculations of grain evolution in the shocked gas. For , we adopt the values from works by Jones et al. (1994) and Jones et al. (1996) derived for steady shock models.

The dust growth rate depends on the choice of the grain size distribution through Eq. (6). For simplicity, we apply a commonly used power law distribution (Mathis et al., 1977, , hereafter called MRN), which runs from  nm to  nm. The value  nm is used in our reference models. To investigate how the minimum grain size affects our results, we perform additional model calculations for  nm and 5 nm and discuss them in Sect. 5.6. The corresponding values of the average grain size are listed in Table 2.

3.3.1 Description of the models

The basic name of our models is ”MRNnm”, where is the minimum grain size in nm and MRN indicates that all models assume MRN grain size distribution. We add ”E” at the beginning of the name of models in which the average timescale of accretion is calculated with account of electrostatic interactions due to ion-grain collisions in the CNM, otherwise a fixed value of is used in all phases (Table 2). The models with additional destruction in the diffuse medium as described in Sect. 3.2 are denoted with an additional prefix ”C”. The model ECMRN3nm is our best fit model for which most of the analysis is given below.

3.3.2 Sticking coefficient

The sticking coefficient is the major source of uncertainty of our model. One-zone dust evolution models usually assume that growth occurs in molecular clouds and key species colliding with the grain surface stick perfectly, so that (O’Donnell & Mathis 1997; Zhukovska et al. 2008; Asano et al. 2013b; see also Weingartner & Draine 1999). Hirashita & Kuo (2011) adopt a more conservative value of . A constant value of has been also adopted in recent numerical simulations of galaxies with dust evolution (e.g., Bekki, 2015b; McKinnon et al., 2016). This assumption is justified at low gas temperature  K and dust temperatures  K because the kinetic energy of incident Si atoms is significantly lower than their binding energy on the surface (e.g., D’Hendecourt et al., 1985). It may however overestimate the dust growth rate in the warm medium in hydrodynamical simulations.

In this work, we include gas-grain interactions in the diffuse ISM where most of the silicon is singly ionised. The sticking coefficient for Si on silicate surface is not known. Watson & Salpeter (1972) discuss possible outcomes of the interaction between a positive ion and a negatively charged grain. They point that when an ion approaches the grain surface, one of the electrons in the grain may tunnel through the work function and neutralise the ion before it reaches the surface. The probability that the atom will remain on the surface upon collision increases dramatically if it can be chemisorbed (Watson & Salpeter, 1972). This possibility is supported by recent experiments in Jena which demonstrate that adsorbed species can form chemical bonds characteristic for silica and silicates at substrate temperatures of 10 K without an activation barrier (Krasnokutski et al., 2014; Rouillé et al., 2014).

The sticking coefficient depends on many parameters including poorly understood surface properties of interstellar grains. With these uncertainties, the observed relation between interstellar depletions and gas density may provide an observational constraint for the choice of . To this end, we investigate various choices of the sticking coefficient: a) a fixed value for  K and for  K; b) a fixed value for all temperatures; c) experimentally measured for chemisorption of H molecules on silicate surface by Chaabouni et al. (2012); and d) from theoretical calculations of physisorption performed by Leitch-Devlin & Williams (1985). In absence of estimates of for physisorption of on silicate grain surfaces, we provisionally adopt the data for carbon atoms arriving on graphite lattice from work by Leitch-Devlin & Williams (1985), in the functional form derived by Grassi et al. (2011). The dust temperature is fixed to 20 K.

We deem the temperature dependence of the sticking coefficient measured by Chaabouni et al. (2012) more plausible than that of the derived by Leitch-Devlin & Williams (1985) (Fig. 1). The former increases at low gas temperatures and approaches 1 at  K, while the latter peaks at and decreases at lower gas temperatures. However, after initial testing and comparison of model predictions with the observational constraints discussed in Sect. 4, we assume the case a) with a simple step function dependence of on for the reference model. It restricts the growth with temperature and does not depend on the exact shape of that has yet to be determined by future experiments.

As will be shown later in the manuscript the assumption of for all temperatures results in too high average depletions compared to the range of values derived from depletion studies and should not be used.

Figure 2: Multiple points with error bars: Gas-phase abundances (Jenkins, 2009) versus the corresponding logarithms of the local gas densities , as inferred from absorption features arising from the C i fine-structure levels in various directions (Jenkins & Tripp, 2011). These measurements are indicated with colors that correspond to estimates for the fractions of the total material that were sampled by the fine-structure population ratios. A solid red line indicates the least-squares linear fit to the data. Single black pentagon: An estimate for the representative abundance for the warm neutral medium (WNM) for a value of equal to 0.12 that corresponds to depletions measured by Savage & Sembach (1996a) for this phase.

4 Depletion of elements in the ISM

In this Section, we discuss observational data used to constrain our model. Interstellar dust abundances can be studied by analysing the gas-phase abundances of refractory elements, assuming that the elements missing from the gas are locked in dust. The abundances of free atoms and ions in space can be determined by analysing absorption features appearing in the ultraviolet spectra of background stars (Spitzer & Jenkins, 1975). The logarithmic depletion of an element in the ISM is the gas-phase element abundance relative to a reference abundance, for which we adopt the element abundance in proto-Sun reflected by the present-day photospheric abundance with an increase of 0.07 dex to account for gravitational settling (Lodders, 2003)


There are significant variations of interstellar element depletions between different lines of sight, which are probably related to dust evolution driven by local conditions and the recent history of the gas. The observed determinations of element depletions can therefore be used to validate and provide useful constraints for our model.

One simple quantity commonly used in the interpretation of interstellar element depletions is the average gas density along the sight line, defined as the ratio of the column density of hydrogen (both atomic and molecular) to the distance of the background star


This is a volume-weighted determination of local densities. It has been well established that strongly correlates with element depletions (e.g., Joseph et al., 1986; Joseph, 1988; Jenkins, 1987; Savage & Sembach, 1996b; Jenkins, 2009; Haris et al., 2016). However, this correlation may be governed more by the relative proportions of some representative high density material compared to very low density gas (Spitzer, 1985; Jenkins et al., 1986), rather than an average for the local densities over regions containing mostly neutral hydrogen along the sight line. In reality, the value of provides a lower limit for the gas densities associated with measured element depletions, with an understanding that the true average hydrogen density for the measurements may be significantly higher if much of the sight line is filled with fully ionised (hence invisible) gas.

The outputs from our dust formation models show details on the distributions of depletions as a function of local density. The lack of an exact correspondence between and representative local densities presents a challenge in our making meaningful comparisons of the observed depletions to the dust models.

To derive the relation between and the mean sight-line density , we start with the data compiled by Jenkins (2009) which included over 243 lines of sight probing a large range of physical conditions. The sight lines with extraordinary high velocities typical for shocked gas were excluded from the analysis. Only stars situated within about 1.5 kpc from the solar Galactocentric radius were used for quantitative analysis, which avoids the possible consequences of the Galactic radial abundance gradient. In the work of Jenkins (2009), the depletions were expressed in terms of a factor designated as , a parameter that characterised the relative strengths of some available element depletions for a given line of sight. Owing to the fact that all of the different element depletions behave in well defined manners with respect to , this parameter can be used as a proxy to determine the depletion of any particular element (in our case ), even if that element was not observed for the sight line under consideration. We may express how relates to using a linear fit to the data plotted Fig. 16 in Jenkins (2009) and the abundance as a function of (by the use of Eq. (10) in Jenkins (2009)), which allows us to derive the following relation between the gas-phase abundance and the average density


There is a more direct way to probe the local gas density in diffuse cold clouds based on an analysis of the collisional excitation of the fine-structure levels of interstellar neutral carbon. Jenkins & Tripp (2011) investigated the C i absorption lines seen in the spectra of stars recorded in the highest resolution echelle modes of the Space Telescope Imaging Spectrograph on board the Hubble Space Telescope. This investigation included 89 out of the 243 lines of sight from the data sample from Jenkins (2009). Using the estimates of thermal gas pressure and temperatures evaluated by Jenkins & Tripp (2011), we calculate the gas densities in the interstellar cold clouds and paired them with the gas-phase abundances along the same sight lines. The outcomes are shown in Fig. 2. The best fit to these results is represented by the equation


and it is shown as a solid red line in the figure. This outcome is based on minimizing the for both the errors in and using the routine fitexy (Press et al., 2007).

As a result of the fact that the ionization equilibrium between C i and C ii is governed by the local electron density, the measurements of the fine-structure level populations for C i are biased in favor of denser gas, which means that the derived hydrogen densities are restricted to the range from about to . It follows that the WNM is not sampled by the C i because the characteristic densities are much lower. For this reason, we must try to understand empirically how influential this WNM contamination is for our determinations of representative local densities of neutral hydrogen that appear in Fig. 2.

For every sight line, Jenkins & Tripp (2011) estimated the fraction of gas that was sampled by the C i fine-structure excitations by examining the absorption profiles of other species in their dominant ionization stages, which could serve as proxies for C II. Based on this information, they stated this sampling fraction in terms of a quantity that they called ”fraction of C II observed”(see their Table 3). If this percentage is low, one must conclude that the WNM dominates over the CNM, and this effect could lead to an outcome for the silicon depletion that is less severe than that for the dense region (Savage & Sembach, 1996a).

To evaluate the importance of the WNM contributions, we re-evaluated the best-fit relationship only for cases where the fraction of C ii in the measurable C i velocity range exceeded 50%. The change of the coefficients in Eq. (20) was inconsequential.

In their review of interstellar abundances, Savage & Sembach (1996a) stated general values for depletions in the WNM in the Galactic disk, which Jenkins converted to a value for the parameter equal to 0.12, which in turn yields . We include in Fig. 2 a point that shows this level of Si depletion at a location for a typical density of the WNM, . The value for the representative density is justified by Jenkins (2013) – see Section 7.4 of that paper. Gas in the Galactic halo probably represents material that has been shocked and accelerated, and here the abundances may be even higher because the grains have been eroded even further.

Figure 3: Evolution of rates of dust destruction in the interstellar shocks and production by growth in the ISM (red thin and blue thick lines, respectively). The rates of dust production in the ISM, destruction by SN shocks and stardust injection rate from 1D dust evolution model (Zhukovska et al., 2008) are shown for comparison (lines with filled triangles, squares and circles, respectively). Top and bottom panels show models with and without additional destruction in the diffuse ISM (ECMRN3nm and EMRN3nm, respectively).
Figure 4: Evolution of the total supernova rate in the simulation volume.

5 Results

In the following, we present the results of the calculations of dust evolution as traced by Si abundance in the Milky-Way-like disk galaxy simulated for 270 Myr. As our goal is to compare the model predictions with the observed depletion data limited to a few kpc around the Sun (Jenkins, 2009), we select the particles in a ring with galactic radii from 6 to 9 kpc representing the conditions similar to those in the solar neighbourhood.

5.1 Dust production rates

Our model predicts the rates of destruction and production of silicate dust with a fixed olivine-type composition. If we assume that the carbon-to-silicate dust mass ratio does not vary significantly across the ISM phases, we can estimate the total rates of dust production and destruction in the ISM. We adopt the carbon-to-silicate ratio of 0.5 inferred for interstellar dust in the local diffuse ISM (Dwek, 2005).

Figure 3 shows the variations of the interstellar dust production and destruction rates per unit of the surface area as a function of time after the start of our calculations. On two panels, we show the results for models EMRN3nm and ECMRN3nm corresponding to the case of dust destruction only by feedback from massive stars in GMCs and the case with the additional destruction in the diffuse medium.

The total dust growth rate is higher at the beginning because of the low values of the initial element depletions in the simulation volume. The initial peak is followed by the dip in the production rate at 90 Myr resulting from the enhanced SN rate (Fig. 4) caused by the initial conditions of the hydrodynamic simulations.

Figure 3 illustrates that, after the initial 140 Myr, the cycle of matter in and out of dust reaches a steady state, in which interstellar dust is distributed over the ISM phases in such a way that the dust production rate balances the destruction rate. The destruction and production rates converge to a similar value in the steady state controlled by the destruction timescale. It is about for models EMRN3nm, which determines a lower limit for the dust destruction rate and corresponds to destruction only by feedback from massive stars in GMCs. For ECMRN3nm model, the value of destruction/production rates in the steady state is by a factor of 2.5 higher and attains the value of .

The resulting total rate of dust destruction/ISM growth for ECMRN3nm model matches remarkably well the value from a simple one-zone model of chemical evolution of the solar neighbourhood with dust derived by Zhukovska et al. (2008). The one-zone model considers the evolution of dust surface densities averaged over the vertical direction as well as in a 1 kpc wide ring with the solar galactocentric radius. It follows the evolution of dust and gas chemical abundances for 13 Gyr, starting from the primordial chemical abundances to the present day values. Time variations of the destruction and ISM dust growth rates predicted by the one-zone model are also displayed in Fig. 3. For comparison with our results, we take the rates from the one-zone model for the last 270 Myr before the present time. These rates from the one-zone model do not noticeably vary over the considered time span because it is much shorter than the present-day timescale of chemical enrichment.

1D dust evolution model also predicts the dust injection rate by AGB stars and SNe based on the star formation history of the the solar neighbourhood and stellar dust yields described in Zhukovska et al. (2008) and Ferrarotti & Gail (2006). The current rate of dust input from stars is (Fig. 3). This is 33 times lower than the dust production rate by ISM growth in our reference model ECMRN3nm. This corroborates our model assumption that growth in the ISM is the dominant dust source in the present-day Milky Way and dust production in stellar winds can be neglected. Note that this assumption does not hold during the early evolution of the Milky Way or in young metal-poor dwarf galaxies (Zhukovska et al., 2008; Zhukovska, 2014).

Figure 5: Evolution of dust production rates in the ISM in model ECMRN3nm separated into gas density bins of width of 1 dex from 0.05 to 500.
Figure 6: Top panel. Evolution of timescales of silicate dust destruction and production by growth in the ISM for model ECMRN3nm (red thin and blue thick lines, respectively). Production rate for model CMRN3nm is shown with dashed line. Bottom panel. The same for the EMRN3nm and MRN3nm models without additional destruction in the diffuse ISM.

5.1.1 Where do interstellar grains grow?

To investigate the conditions that favour dust growth in our hydrodynamical simulations, we compare the dust production rates within the logarithmic gas density intervals from 0.05 to 500 with a stepsize of 1 dex. Time evolution of the rates in each density bin for ECMRN3nm model is shown in Fig. 5. It demonstrates that most of growth occurs at densities , which is 49% of the total ISM growth rate. Gas with densities contributes 30% and more diffuse clouds with 17% to the total growth rate, correspondingly. The contribution from the diffuse medium with constitutes only 3% of the total rate. This is not surprising, since the growth is limited to  K in the reference model.

Our results do not drastically change for models with different temperature-dependent sticking coefficient tested in this work (Chaabouni et al., 2012; Leitch-Devlin & Williams, 1985). Compared to a simple step function dependence on temperature in the reference models, the sticking coefficient usually decreases with gas temperatures from the maximum value close to 1 similar to the behavior of the data from Chaabouni et al. (2012) shown in Fig. 1. Such dependence somewhat reduces the growth rate in the density bin and increases it in the next density bin () resulting in very similar values of growth rates in these two bins.

Our simulations do not resolve the dense gas well, therefore its actual contribution can be higher than the present value of a few per cent.

Figure 7: Probability density functions of the depletions per unit log depletion for the final simulation snapshot split in different gas density intervals spanning from 0.05 to 500(solid lines with various symbols, as indicated). Black line shows the PDF for all densities. All PDFs are normalised to the total gas mass in the considered volume. Top and bottom panels display models CMRN3nm and ECMRN3nm, respectively.

5.2 Timescales of dust destruction and re-formation in the ISM

Figure 6 shows the evolution of the average timescales of silicate dust destruction and re-formation in the ISM derived using the dust production and destruction rates for models ECMRN3nm and EMRN3nm (Fig. 3). The early evolution of the dust formation timescale is strongly affected by the initial conditions for both the hydrodynamical simulations and dust abundances. The adopted constant depletion at the onset of the simulations causes a rapid growth phase with the formation timescale of a few hundred Myr. Later on, the growth is slowed down as a consequence of a temporal reduction of the GMC mass fraction by feedback from the initial burst of star formation. This reduction in the star formation activity manifests in Fig. 4 as a dip in the SN rate following the peak corresponding to the initial starburst. Strong feedback from the starburst suppresses the growth in clouds and increases to over 1 Gyr by a time of 90 Myr in both models. However, the steady state values of the destruction and formation timescales, reached after the balance between disruption and formation of GMCs is established, differ significantly for the two models.

The average timescale of dust destruction calculated for model ECMRN3nm is 3 times shorter than for model EMRN3nm and reaches down to 350 Myr. For the adopted low value of the timescale of the destruction rate in the diffuse medium  Myr, is determined by the and the mass fraction of the diffuse ISM in our simulations which is about 20%. Efficient destruction in the diffuse phase with a short , resulting in the short lifetime of grains in the ISM, are required for our model to reproduce the mean depletion in the diffuse medium. Our results thus re-enforce the assumption proposed in earlier studies that sputtering of grains in the diffuse ISM by single SNe is the dominant mechanism of dust destruction (McKee, 1989).

When only destruction by SN feedback in GMCs is included (model EMRN3nm), the lifetime of grains in the ISM at the end of simulations reaches 1 Gyr. Model EMRN3nm overpredicts the high depletions of Si in the WNM compared to the observed value, and , respectively. The value of 350 Myr predicted by model ECMRN3nm at the end of simulations is 50 Myr lower that the lifetime of silicate grains in the Milky Way estimated by Jones et al. (1996). Following McKee (1989), they combined the observed mass of the ISM and galactic SN rate with their estimate of the given by Eq. (12) resulting from extensive numerical calculations of the dust destruction in the blast wave.

We can also compare the timescale from our simulations with other studies which improved various aspects of the physics of grain destruction in a shock, but kept the same simple approach to the grain cycle in the homogeneous ISM. Bocchio et al. (2014), for example, re-evaluated the dust lifetimes against destruction using models with a better treatment of the dust dynamics in the shock. With consideration of uncertainties in the observed mass of the ISM and SN rate, they derived the lifetime of  Myr for silicate grains, which is very close to the value predicted by our best fit model. Slavin et al. (2015) proposed much longer lifetimes of silicate grains of 2–3 Gyr. Although they included the hydrodynamical shock evolution in calculations of the destruction efficiency , the new value of from their work is not very different from the previous works. The main reason for their longer is the new estimates of the total gas mass and SN rate from recent observations.

The average dust formation timescale in the ISM for models with enhanced collision rates of cations with grains in the CNM is not very different from those without it (Fig 6). The differences are larger for models with additional destruction in the diffuse medium. They can be explained by differences in and element depletion levels in these models. The growth timescale in the CNM is 4 times shorter for model ECMRN3nm which results in rapid depletion of from the gas phase to higher levels than in model CMRN3nm. Since increases with in individual gas particles and decreases with their element condensation fractions , see Equations (5) and (6), the interplay between these two factors determines the ratio of predicted by the models with and without enhanced collision rates.

5.3 Distribution of Si depletion with density

We analyse the complex distribution of Si depletion values in the simulation volume by means of mass-weighted probability density functions (PDFs) calculated in logarithmically equal gas density intervals. Figure 7 displays the PDFs of the final distribution of Si depletions for the logarithmic gas density intervals from 0.05 to 500 with steps of 1 dex. All PDFs are normalised to the total mass of gas. The main feature of the derived distributions is broad asymmetric profiles which become shallower in higher gas density bins as a consequence of the dust growth. Increase is commensurate with the observed depletion trends discussed in detail in Sect. 5.4.

The broadest depletion distribution derived for the densities reflects rich dynamic histories of particles at these densities. Particles with the highest depletions belong to clouds which were formed or evolved from the disruption of larger clouds. They have spent longer time at larger densities. On the other hand, the low density tail of the PDF corresponds to gas which has recently undergone compression from the diffuse gas so that the depletion values have not yet substantially increased. The highly asymmetric profiles of PDFs shown in Fig. 7 with a steady increase towards the less severe depletions are a direct consequence of the additional dust destruction in the diffuse phase. The mean depletions for these PDFs nevertheless provide the best match to the observed depletion values.

For the selected minimum grain size of 3 nm, the enhanced collision rates due to electrostatic focusing modestly affects the PDFs of Si depletion (Fig. 7, bottom panel). This process is responsible for the shallower slopes of PDFs for the density range owing to the enhanced growth under these conditions. The PDF for looks lower for the ECMRN3nm model, because a larger fraction of gas at these densities has depletion values below .

Figure 8: Relation between mean Si depletion and gas density derived using the final PDFs of gas-phase abundances for models CMRN3nm and ECMRN3nm (solid and dashed black lines, respectively). The filled areas around the synthetic relations indicate the corresponding standard deviations. The pentagon symbol indicates a generalised depletion level at a location for a typical density of the WNM from observations. The thick red line shows the linear fit to the observed data given by Eq. (20) and the thin red line shows the lower limit for this relation derived for the mean gas density on the line of sight given by Eq. (19) (see explanation in Sect. 4). The yellow area shows the range of values between two relations. Typical values for the warm and cool galactic disk (open and filled diamond, respectively) derived by Savage & Sembach (1996a) are shown for comparison.
Figure 9: The same as in Fig. 8 compared to the relations between mean Si depletion and gas density derived for models MRN3nm and EMRN3nm (blue solid and dashed lines with circles, respectively) without additional dust destruction in the diffuse ISM. We do not show the standard deviations for the synthetic relations for clarity of the figure.

5.4 Si depletion – gas density relation

In order to compare our model predictions with the observed relation between Si depletion and gas density (Sect. 4), we compute the mean values and standard deviations of the interstellar Si depletion using the PDFs for the final simulation snapshot. For this analysis, the PDFs are computed in the same way as described in the previous section, but within much finer logarithmic gas density intervals. The resulting relations between the gas phase Si abundance and gas density are shown in Fig. 8 for reference models ECMRN3nm and CMRN3nm, in which a simple step function is adopted for the sticking coefficient ( for  K and for higher ).

Figure 8 also shows the relations between the mean values and the gas density derived from observations. The relation based on the gas densities probed by C i fine structure lines (Eq. (19)) is limited to a narrow density range of . A generalised Si depletion level in the WNM inferred from observations provides a constraint for dust evolution models in the WNM density regime. Additionally, Figure 8 shows the relation derived for the mean gas density on the lines of sight given by Eq (19), which we consider as a lower limit to the observed relation.

A conspicuous feature of the synthetic mean relation for all models is the large standard deviation of about 0.5 dex at the lowest gas densities up to one dex at higher densities. These large variations predicted by the dust evolution models agree well with the dispersion in the real values from observations (Fig. 2). The dispersion arises from different dynamical histories of gas parcels with the same density that undergo various stages of the lifecycle of GMCs. Gas that previously resided in GMCs on its way to the WNM has higher Si depletions, while matter from the WNM has lower element depletions because of dust destruction process. This effect can be also seen in the widths of PDFs calculated within the larger density bins discussed above.

Another characteristic feature of the synthetic mean relation is its different slopes at the low and high gas densities seen as the double power law in Fig. 8. The slope in the dense gas regime is steeper and its value of 0.5 is in excellent agreement with the slope of the observed relation. This agreement provides the evidence of dust growth by accretion in the ISM. With enhanced collision rates due to ion-grain interactions (model ECMRN3nm), the accretion timescale is shorter and the slope in the diffuse gas is therefore steeper compared to the CMRN3nm model.

Comparison of the slopes at low densities for models ECMRn3nm and CMRN3nm reveals that this slope also depends on the accretion timescale, although the growth is limited to temperatures below 300 K in these models. This is also supported by the fact that dust destruction in the diffuse phase does not depend on the local density in our model. The slope at low density regime is determined by the dust depletion at high densities, which is controlled by the growth timescale. Additional important factor controlling the distribution of dust abundances at low densities is the mass circulation between the GMCs and diffuse medium. In this work, it is set by numerical simulations of GMC evolution and is the same in both dust models, hence the differences in the slope at low densities between models CMRN3nm and EMRN3nm as well as in other models discussed in the remainder of this paper are the ”memory” of dust growth by accretion at higher densities.

In the reference models described above we include dust destruction in the diffuse ISM in addition to destruction by type II SNe in GMCs. We compare the synthetic relations between and obtained from the reference models with those from models MRN3nm and EMRN3nm in Fig. 9. All parameters are the same with the exception of dust destruction, which in models MRN3nm and EMRN3nm occurs only through feedback from massive stars in GMCs. Models without additional destruction in the diffuse ISM overpredict Si depletion by 0.5 dex in the case without ion-grain interactions included and by 1 dex in the model EMRN3nm with enhanced collision rates due to ion-grain interactions.

The observed distribution of gas phase element abundances thus implies that some fraction of stellar feedback energy is injected in the diffuse medium. SN energy input only in the dense phase results in too low values compared to the observations. The importance of stellar feedback energy injection in the diffuse ISM is also stressed in work by Walch et al. (2015), who analyse how location of SN explosions affects the structure of the ISM by means of high resolution hydrodynamical simulations of the ISM.

The present dust evolution model does not include dust and gas input from stars. Since the volume filling factor of the diffuse medium is much larger than that of the molecular clouds, evolved stars inject matter mainly in the diffuse phase. For stardust evolution model of Zhukovska et al. (2008), the dust input rate from stars is significantly lower than the rate of dust mass growth in the ISM (see Sect. 5.1) and in stellar ejecta is higher than the observed value in the warm phase. Therefore, adding stellar dust sources self-consistently would result in slightly higher values of in the diffuse gas, but would not change the slope of the relation.

In this work, we focus on the matter cycle at the solar galactocentric radius. Dobbs & Pringle (2013) demonstrate that the dominant mechanisms of disruption and formation of the GMCs are different in the outer and inner disk. Clouds in the inner disk are more likely to get destroyed by stellar feedback and sheer, while large scale gas dynamics is more important on the larger radii. The radial distribution of dust and, correspondingly, element depletions in the gas phase are affected by the changes in dynamical evolution of GMCs with distance from the Galactic centre, in addition to the metallicity radial gradients. Our preliminary study indicates that shorter lifetimes of the clouds in the outer disk result in lower element depletions. We investigate how the radial variations in the dynamical evolution of clouds impact the dust abundances and the relation in a forthcoming paper.

Figure 10: Relation between Si depletion and gas density derived using the PDFs of gas-phase abundances for the final snapshot in the simulations for different sticking coefficients: from work of Leitch-Devlin & Williams (1985), Chaabouni et al. (2012) and (lines with triangles, squares and crosses, respectively). Solid and dashed lines show models CMRN3nm and ECMRN3nm, respectively. The observational data are displayed in the same way as in Fig. 8.

5.5 Observational constraints for the sticking coefficient

The sticking coefficient affects the final distribution of the interstellar depletions by limiting the densities of the ISM at which dust mass can grow by accretion of gas-phase species. In order to analyse how the modelled depletions depend on the choice of the sticking coefficient, we perform the simulations for models with the same parameters as ECMRN3nm and CMRN3nm, but assuming different as discussed in Sect. 3.3.2. The resulting relations for the mean depletion as a function of gas density are shown in Fig. 10.

The reference model includes the dust growth only for  K with the maximum value of . If the temperature restriction is lifted, the theoretical relation shifts by at low densities and by at high densities. To balance dust growth at high temperatures and to lower depletions in the diffuse ISM in models MRN3nm and EMRN3nm and to match the observational data would require an unreasonably short timescale of destruction in the diffuse medium in Eq (16). We therefore conclude that the assumption of the maximum coefficient for all temperatures, commonly used in recent hydrodynamical simulation with dust (e.g., McKinnon et al., 2016; Bekki, 2015a, b), tends to overestimate the dust production rates in hydrodynamic numerical simulations of galactic evolution.

The sticking coefficients derived by Chaabouni et al. (2012) and by Leitch-Devlin & Williams (1985) have different dependences on the gas temperature. The latter has the maximum at about 200 K, while the former peaks at  K. The fact that both prescriptions result in the similar relations reveals that the mean depletions are not very sensitive to the choice of , as long as it decreases with gas temperature. CMRN3nm models with from Chaabouni et al. (2012) and Leitch-Devlin & Williams (1985) yield too high values and a shallower slope than the values from observational data. Calculations for EMRN3nm models provide a better fit to the observations than the CMRN3nm model.

Figure 11: Relation between Si depletion and gas density derived from the final distribution of Si gas-phase abundances in the simulations for minimum grain size , 3, and 5 nm (lines with squares, circles, and triangles, respectively) for CMRN and ECMRN models (solid and dashed lines, respectively). The observed relation is shown in the same way as in Fig. 8.

5.6 Dependence on the grain size distribution

In the previous sections we discussed various results of dust modelling for a fixed grain size distribution. The size distribution for silicate grains, in particular its lower limit, is highly uncertain parameter and varies among different models. It is determined by fitting the observed spectral energy distribution of interstellar dust in the diffuse medium (Zubko et al., 2004; Siebenmorgen et al., 2014). The grain size distribution enters the rate of the dust growth in the ISM via Eqs. (59). Variations in the grain sizes affect the final element depletions through the change in the total surface area, and additionally, through enhanced collision rates with smaller grains (Weingartner & Draine, 1999).

In the following, we inspect how the variations in the lower limit of the grain size distribution affects the relation between the average depletion and gas density predicted by our models. We use the same value m for the upper limit for the grain sizes, as it has has little effect on the total surface area. All other parameters are taken the same as for the reference model. The results for models CMRN and ECMRN with , 3 and 5 nm for the final simulation snapshot are shown in Fig. 11.

The synthetic average vs. relations show double power law behaviour noted above for the reference models. With the decrease of the minimum grain size, the total surface area of grains increases and the accretion in the ISM occurs faster resulting in steeper slopes of the relation and transition to the efficient growth phase at lower densities. This effect is additionally enhanced by the presence of negatively charged small grains if the collision rates of grains with ions are included (see Fig. 1 in Weingartner & Draine, 1999). The differences in depletion values between ECMRN and CMRN models may reach 1 dex at high densities for  nm and  nm (Fig. 11). The ratio of the gas densities at which the transition to effective growth takes place roughly corresponds to the ratio of the growth timescales that for models shown in Fig. 11 is determined by the average grain radius (Table 2).

We do not confirm existence of the so-called Coulomb barrier for dust growth suggested recently by Ferrara et al. (2016) as a process to prevent dust growth. In contrary, we find that growth by accretion in the CNM occurs even for assumption of the “classical” cut-off of the MRN size distribution  nm due to the fact that only silicate grains larger than 20 nm are positively charged. In this case, however, repulsion between ions and positively charged grains results in a somewhat longer accretion timescale and shallower slopes of the depletion–gas density relation than in model CMRN5nm. The mean values in the dense gas for models CMRN5nm and ECMRN5nm are higher than the observational data (Fig. 11), we therefore favour the models with  nm.

We conclude that silicate grains with the sizes that are smaller than the cut-off of the MRN size distribution of  nm are necessary to reproduce the observed slope of the relation in our model. Underlying this result is the assumption that some chemical selection process prevents sticking of on very small carbonaceous grains such as desorption upon photoexcitation by UV irradiation (Tielens, 1998; Draine, 2009).

We do not take into account that the grain size distribution evolves in dense cores of molecular clouds owing to coagulation process. Hirashita (2012) showed that because element depletion by accretion occurs at lower gas densities than grain coagulation, the reduction of total grain surface area by coagulation only slightly affects the dust mass growth in the ISM. Indeed, we find that efficient dust growth occurs already in the CNM and consumes almost all gas-phase at much higher densities than the coagulation regime ().

6 Conclusions

We develop a three-dimensional model for interstellar dust evolution as traced by abundance based on numerical hydrodynamical simulations of giant molecular clouds in a Milky-Way-like galaxy. This work focuses on the model predictions in a ring centred on the solar galactocentric radius so that they can be tested against the gas-phase abundances measured along numerous sight lines and probing various interstellar conditions in the local Milky Way. Using the depletion data and measurements of the local gas density probed by fine-structure lines of neutral carbon, we derive a relation between the average Si abundance in the gas and the local gas density which we use as a critical constraint for the models.

We demonstrate that the present rate of dust re-formation in the ISM of in the solar neighbourhood is substantially larger than the dust injection rate by stars (). Most of dust growth occurs in the density range accounting for a half of the total growth rate. We find that the dust mass growth takes place in the cold neutral medium for all considered values of the minimum grain size from 1 to 5 nm, despite the Coulomb repulsion between positive impinging ions and large positively charged grains. Because of negatively charged small grains, the accretion timescale for  nm is shorter with account of Coulomb interactions.

The model includes destruction of grains by energetic feedback from SNe that occur mostly in the GMCs in the hydrodynamical simulations used for the dust modelling. This destruction process has a relatively long timescale of 1 Gyr and tends to over-predict the dust abundances in the WNM. The latter constraint can be only reproduced with the implementation of an additional destruction process in the diffuse ISM representing non-correlated SNe from run-away massive stars and type Ia SNe. We derive an average lifetime of silicate grains against destruction of 350 Myr for our best-fit model. Our findings re-enforce the assumption that destruction of grains in the diffuse ISM by single SNe is a dominant pathway of dust destruction.

We calculate the probability density distributions of the and their corresponding mean values and standard deviations for different gas densities. The resulting relation between the average and the gas density derived for dust distribution at the end of the simulations is used to constrain the dust growth model by comparison with the observed relation. We find that the assumption of a constant sticking coefficient adopted in recent hydrodynamics models with dust significantly overestimates the dust growth rates. The observed relation requires a sticking coefficient that decreases with temperature. A simple assumption that the growth is limited to the CNM and clouds with  K provides a slightly better fit in the CNM density range than the provisional models with temperature-dependent sticking coefficients from Chaabouni et al. (2012) and Leitch-Devlin & Williams (1985). Including the enhanced collision rates of cations with grains in the CNM and adopting the lower limit for the grain size of  nm in the dust growth model are crucial to reproduce the observed slope of the relation.

The dust evolution model proposed in this work does not require the presence of a population of very small silicate grains like the one known to exist for carbonaceous grains (Draine & Anderson, 1985; Draine & Li, 2001). Including the 1 nm grains in the model leads to overly high depletions in the entire simulation volume. If such a population exists, then an additional mechanism, other than sputtering in SN shocks, has to be included in the model to remove accreted from the grain surfaces and provide the observed levels of depletion in the ISM.

S.Z. acknowledges support by the Forschungsgemeinschaft through SPP 1573: “Physics of the Interstellar Medium”. We are grateful to Nikolai Voshchinnikov and Simon Glover for reading the manuscript. We thank the referee for useful comments that helped to improve the manuscript. We gratefully acknowledge the Max Planck Computing and Data Facility for providing their user support and computing time on the Odin and Hydra clusters. C.D. acknowledges support by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2011-2016 grant agreement no. 280104, LOCALSTAR). The findings reported in Section 4 of this paper arose from conclusions supported by HST archival program numbers 09534 and 10279, which were provided by NASA through grants from the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.


  • Asano et al. (2013a) Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013a, Earth, Planets and Space, 65, 213
  • Asano et al. (2013b) Asano, R. S., Takeuchi, T. T., Hirashita, H., & Nozawa, T. 2013b, Monthly Notices of the Royal Astronomical Society, 432, 637
  • Bate et al. (1995) Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, Monthly Notices of the Royal Astronomical Society, 277, 362
  • Bekki (2013) Bekki, K. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2254
  • Bekki (2015a) —. 2015a, The Astrophysical Journal, 799, 166
  • Bekki (2015b) —. 2015b, Monthly Notices of the Royal Astronomical Society, 449, 1625
  • Benz et al. (1990) Benz, W., Cameron, A. G. W., Press, W. H., & Bowers, R. L. 1990, The Astrophysical Journal, 348, 647
  • Binney & Tremaine (1987) Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton: Princeton University Press)
  • Bocchio et al. (2014) Bocchio, M., Jones, A. P., & Slavin, J. D. 2014, Astronomy and Astrophysics, 570, A32
  • Calura et al. (2008) Calura, F., Pipino, A., & Matteucci, F. 2008, Astronomy and Astrophysics, 479, 669
  • Calzetti (2013) Calzetti, D. 2013, in Secular Evolution of Galaxies, ed. J. Falcón-Barroso & J. H. Knapen (Secular Evolution of Galaxies), 419
  • Chaabouni et al. (2012) Chaabouni, H., Bergeron, H., Baouche, S., et al. 2012, Astronomy and Astrophysics, 538, A128
  • Christensen et al. (2012) Christensen, C., Quinn, T., Governato, F., et al. 2012, Monthly Notices of the Royal Astronomical Society, 425, 3058
  • Cox & Gomez (2002) Cox, D. P., & Gomez, G. C. 2002, THE ASTROPHYSICAL JOURNAL SUPPLEMENT SERIES, 142, 261
  • Dartois et al. (2004) Dartois, E., Mu Oz Caro, G. M., Deboffle, D., & D’hendecourt, L. 2004, Astronomy and Astrophysics, 423, L33
  • de Bennassuti et al. (2014) de Bennassuti, M., Schneider, R., Valiante, R., & Salvadori, S. 2014, Monthly Notices of the Royal Astronomical Society, 445, 3039
  • D’Hendecourt et al. (1985) D’Hendecourt, L. B., Allamandola, L. J., & Greenberg, J. M. 1985, Astronomy and Astrophysics, 152, 130
  • Dobbs (2008) Dobbs, C. L. 2008, Monthly Notices of the Royal Astronomical Society, 391, 844
  • Dobbs et al. (2011) Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, Monthly Notices of the Royal Astronomical Society, 417, 1318
  • Dobbs & Pringle (2013) Dobbs, C. L., & Pringle, J. E. 2013, Monthly Notices of the Royal Astronomical Society, 432, 653
  • Dorschner & Henning (1995) Dorschner, J., & Henning, T. 1995, The Astronomy and Astrophysics Review, 6, 271
  • Draine (1990) Draine, B. T. 1990, in Astronomical Society of the Pacific, The evolution of the interstelllar medium, ed. L. Blitz, 193–205
  • Draine (2009) Draine, B. T. 2009, in ASP Conf. Ser. 414, Cosmic Dust - Near And Far, ed. T. Henning, E. Grün, & A. Steinacker, 453
  • Draine & Anderson (1985) Draine, B. T., & Anderson, N. 1985, The Astrophysical Journal, 292, 494
  • Draine & Li (2001) Draine, B. T., & Li, A. 2001, The Astrophysical Journal, 551, 807
  • Draine & Salpeter (1979) Draine, B. T., & Salpeter, E. E. 1979, Astrophysical Journal, 231, 438
  • Dwek (1998) Dwek, E. 1998, Astrophysical Journal, 501, 643
  • Dwek (2005) Dwek, E. 2005, in AIP Conference Proceedings, ed. C. Popescu & R. Tuffs (AIP), 103–122
  • Dwek et al. (2007) Dwek, E., Galliano, F., & Jones, A. P. 2007, The Astrophysical Journal, 662, 927
  • Dwek & Scalo (1980) Dwek, E., & Scalo, J. M. 1980, Astrophysical Journal, 239, 193
  • Ferrara et al. (2016) Ferrara, A., Viti, S., & Ceccarelli, C. 2016, 1606.07214
  • Ferrarotti & Gail (2006) Ferrarotti, A. S., & Gail, H.-P. 2006, Astronomy and Astrophysics, 447, 553
  • Fitzpatrick & Massa (2007) Fitzpatrick, E. L., & Massa, D. 2007, The Astrophysical Journal, 663, 320
  • Forbes et al. (2016) Forbes, J. C., Krumholz, M. R., Goldbaum, N. J., & Dekel, A. 2016, eprint arXiv:1605.00650, 1605.00650
  • Glover & Mac Low (2007) Glover, S. C. O., & Mac Low, M.-M. 2007, THE ASTROPHYSICAL JOURNAL SUPPLEMENT SERIES, 169, 239
  • Grassi et al. (2011) Grassi, T., Krstic, P., Merlin, E., et al. 2011, 533, A123
  • Haris et al. (2016) Haris, U., Parvathi, V. S., Gudennavar, S. B., et al. 2016, The Astronomical Journal, 151, 143
  • Hirashita (1999) Hirashita, H. 1999, The Astrophysical Journal, 522, 220
  • Hirashita (2012) —. 2012, Monthly Notice of the Royal Astronomical Society, 422, 1263
  • Hirashita & Kuo (2011) Hirashita, H., & Kuo, T.-M. 2011, Monthly Notice of the Royal Astronomical Society, 416, 1340
  • Hirashita et al. (2015) Hirashita, H., Nozawa, T., Villaume, A., & Srinivasan, S. 2015, Monthly Notices of the Royal Astronomical Society, 454, 1620
  • Hollenbach & Salpeter (1971) Hollenbach, D., & Salpeter, E. E. 1971, The Astrophysical Journal, 163, 155
  • Hu et al. (2016) Hu, C.-Y., Naab, T., Walch, S., Glover, S. C. O., & Clark, P. C. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3528
  • Jenkins (1987) Jenkins, E. B. 1987, Interstellar Processes, 134, 533
  • Jenkins (2009) —. 2009, The Astrophysical Journal, 700, 1299
  • Jenkins (2013) —. 2013, The Astrophysical Journal, 764, 25
  • Jenkins et al. (1986) Jenkins, E. B., Savage, B. D., & Spitzer, L. J. 1986, The Astrophysical Journal, 301, 355
  • Jenkins & Tripp (2011) Jenkins, E. B., & Tripp, T. M. 2011, The Astrophysical Journal, 734, 65
  • Jones et al. (1996) Jones, A., Tielens, A., & Hollenbach, D. 1996, The Astrophysical Journal, 469, 740
  • Jones et al. (1994) Jones, A. P., Tielens, A. G. G. M., Hollenbach, D. J., & McKee, C. F. 1994, Astrophysical Journal, 433, 797
  • Joseph (1988) Joseph, C. 1988, The Astrophysical Journal, 335, 157
  • Joseph et al. (1986) Joseph, C. L., Snow, T. P., Seab, C. G., & Crutcher, R. M. 1986, Astrophysical Journal, 309, 771
  • Klessen & Glover (2016) Klessen, R. S., & Glover, S. C. O. 2016, Star Formation in Galaxy Evolution: Connecting Numerical Models to Reality, Saas-Fee Advanced Course, Volume 43. ISBN 978-3-662-47889-9. Springer-Verlag Berlin Heidelberg, 2016, p. 85, 43, 85
  • Krasnokutski et al. (2014) Krasnokutski, S. A., Rouillé, G., Jäger, C., et al. 2014, The Astrophysical Journal, 782, 15
  • Leitch-Devlin & Williams (1985) Leitch-Devlin, M. A., & Williams, D. A. 1985, Monthly Notices of the Royal Astronomical Society, 213, 295
  • Lodders (2003) Lodders, K. 2003, The Astrophysical Journal, 591, 1220
  • Lodders et al. (2009) Lodders, K., Palme, H., & Gail, H.-P. 2009, in Solar System, ed. J. E. Trümper (Berlin, Heidelberg: Springer Berlin Heidelberg), 712–770
  • Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, Astrophysical Journal, 217, 425
  • Mattsson & Andersen (2012) Mattsson, L., & Andersen, A. C. 2012, Monthly Notices of the Royal Astronomical Society, 423, 38
  • McKee (1989) McKee, C. 1989, Interstellar Dust: Proceedings of the 135th Symposium of the International Astronomical Union, 135, 431
  • McKinnon et al. (2016) McKinnon, R., Torrey, P., & Vogelsberger, M. 2016, Monthly Notices of the Royal Astronomical Society, 457, 3775
  • Mihalas & Binney (1981) Mihalas, D., & Binney, J. 1981, Galactic astronomy: Structure and kinematics, 2nd edn., Vol. 1 (San Francisco, CA, W. H. Freeman and Co., 1981. 608 p.)
  • O’Donnell & Mathis (1997) O’Donnell, J. E., & Mathis, J. S. 1997, Astrophysical Journal, 479, 806
  • Oey & Lamb (2012) Oey, M. S., & Lamb, J. B. 2012, Four Decades of Research on Massive Stars ASP Conference Series, 465, 431
  • Press et al. (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (Oracle Corporation: Cambridge University Press)
  • Price & Monaghan (2007) Price, D. J., & Monaghan, J. J. 2007, Monthly Notices of the Royal Astronomical Society, 374, 1347
  • Reach et al. (2015) Reach, W. T., Heiles, C., & Bernard, J.-P. 2015, The Astrophysical Journal, 811, 118
  • Rémy-Ruyer et al. (2014) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, Astronomy and Astrophysics, 563, A31
  • Roman-Duval et al. (2014) Roman-Duval, J., Gordon, K. D., Meixner, M., et al. 2014, The Astrophysical Journal, 797, 86
  • Rouillé et al. (2014) Rouillé, G., Jäger, C., Krasnokutski, S. A., Krebsz, M., & Henning, T. 2014, Faraday Discussions, 168, 449
  • Roy et al. (2013) Roy, A., Martin, P. G., Polychroni, D., et al. 2013, The Astrophysical Journal, 763, 55
  • Santini et al. (2014) Santini, P., Maiolino, R., Magnelli, B., et al. 2014, Astronomy and Astrophysics, 562, A30
  • Savage & Sembach (1996a) Savage, B. D., & Sembach, K. R. 1996a, Annual Review of Astronomy and Astrophysics, 34, 279
  • Savage & Sembach (1996b) —. 1996b, Astrophysical Journal, 470, 893
  • Seab & Shull (1983) Seab, C. G., & Shull, J. M. 1983, The Astrophysical Journal, 275, 652
  • Serra Díaz-Cano & Jones (2008) Serra Díaz-Cano, L., & Jones, A. P. 2008, Astronomy and Astrophysics, 492, 127
  • Siebenmorgen et al. (2014) Siebenmorgen, R., Voshchinnikov, N. V., & Bagnulo, S. 2014, Astronomy and Astrophysics, 561, A82
  • Slavin et al. (2015) Slavin, J. D., Dwek, E., & Jones, A. P. 2015, The Astrophysical Journal, 803, 7
  • Spitzer (1985) Spitzer, L. J. 1985, The Astrophysical Journal, 290, L21
  • Spitzer & Jenkins (1975) Spitzer, L. J., & Jenkins, E. B. 1975, Annual Review of Astronomy and Astrophysics, 13, 133
  • Tielens (1998) Tielens, A. G. G. M. 1998, The Astrophysical Journal, 499, 267
  • Voshchinnikov & Henning (2010) Voshchinnikov, N. V., & Henning, T. 2010, Astronomy and Astrophysics, 517, A45
  • Walch et al. (2015) Walch, S., Girichidis, P., Naab, T., et al. 2015, Monthly Notices of the Royal Astronomical Society, 454, 238
  • Watson & Salpeter (1972) Watson, W. D., & Salpeter, E. E. 1972, The Astrophysical Journal, 174, 321
  • Weingartner & Draine (1999) Weingartner, J., & Draine, B. 1999, The Astrophysical Journal, 517, 292
  • Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, THE ASTROPHYSICAL JOURNAL SUPPLEMENT SERIES, 134, 263
  • Yozin & Bekki (2014) Yozin, C., & Bekki, K. 2014, Monthly Notices of the Royal Astronomical Society, 439, 1948
  • Zhukovska (2008) Zhukovska, S. 2008, PhD thesis, Ruperto-Carola-Univesity, Heidelberg
  • Zhukovska (2014) —. 2014, Astronomy and Astrophysics, 562, A76
  • Zhukovska & Gail (2009) Zhukovska, S., & Gail, H.-P. 2009, Astronomical Society of the Pacific Conference Series, 414, 199
  • Zhukovska et al. (2008) Zhukovska, S., Gail, H.-P., & Trieloff, M. 2008, Astronomy and Astrophysics, 479, 453
  • Zhukovska & Henning (2013) Zhukovska, S., & Henning, T. 2013, Astronomy and Astrophysics, 555, A99
  • Zubko et al. (2004) Zubko, V., Dwek, E., & Arendt, R. G. 2004, THE ASTROPHYSICAL JOURNAL SUPPLEMENT SERIES, 152, 211
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description