Shiva: the dust destruction model

TL;DR

Shiva is a numerical tool that simulates dust destruction and evolution in various astrophysical environments, accounting for physical processes like sputtering, shattering, and photo-processing, and models changes in grain properties over time.

Contribution

The paper introduces Shiva, a novel computational tool that models dust grain evolution and destruction, including chemical and physical transformations, in different interstellar media.

Findings

Shiva can simulate dust evolution in diverse environments.

It predicts dust grain lifetimes under various conditions.

The tool models spectral changes due to dust processing.

Abstract

We present a numerical tool Shiva designed to simulate the dust destruction in warm neutral, warm ionized, and hot ionized media under the influence of photo-processing, sputtering, and shattering. The tool is designed primarily to study the evolution of hydrogenated amorphous carbons (HACs), but options to simulate polycyclic aromatic hydrocarbons (PAHs), silicate and graphite grains are also implemented. HAC grain photo-processing includes both dehydrogenation and carbon atom loss. Dehydrogenation leads to material transformation from aliphatic to aromatic structure. Simultaneously, some other physical properties (band gap energy, optical properties, etc.) of the material change as well. The Shiva tool allows calculating the time-dependent evolution of the dust size distribution depending on hydrogen, helium, and carbon number densities and ionization state, gas temperature, radiation flux, relative gas-dust and grain-grain velocities. For HAC grains the evolution of band gap energy distribution is also computed. We describe a dust evolution model, on which the tool relies, and present evolutionary time-scales for dust grains of different sizes depending on external conditions. This allows a user to estimate quickly a lifetime of a specific dust grain under relevant conditions. As an example of the tool usage, we demonstrate how grain properties and corresponding infrared spectra evolve in photo-dissociation regions, H II regions, and supernova remnant shocks.

Equations

$$\begin{array}{rl}\frac{d{N}_{ij}}{dt}& =\underset{\text{aromatisation by photons}}{\underbrace{A_{ij+1}^{(1)}{N}_{ij+1}-A_{ij}^{(1)}{N}_{ij}}}+\underset{\text{photo-destruction}}{\underbrace{B_{i+1j}^{(1)}{N}_{i+1j}-B_{ij}^{(1)}{N}_{ij}}}\\ & +\underset{\text{aromatisation by ions in inelastic collisions}}{\underbrace{A_{ij+1}^{(2)}{N}_{ij+1}-A_{ij}^{(2)}{N}_{ij}}}+\underset{\text{sputtering due to inelastic collisions}}{\underbrace{B_{i+1j}^{(2)}{N}_{i+1j}-B_{ij}^{(2)}{N}_{ij}}}\\ & +\underset{\text{aromatisation by ions in elastic collisions}}{\underbrace{A_{ij+1}^{(3)}{N}_{ij+1}-A_{ij}^{(3)}{N}_{ij}}}+\underset{\text{sputtering due to elastic collisions}}{\underbrace{B_{i+1j}^{(3)}{N}_{i+1j}-B_{ij}^{(3)}{N}_{ij}}}\\ & +\underset{\text{aromatisation of big grains due to sputtering}}{\underbrace{A_{ij+1}^{(4)}{N}_{ij+1}-A_{ij}^{(4)}{N}_{ij}}}+\underset{\text{sputtering of big grains}}{\underbrace{B_{i+1j}^{(4)}{N}_{i+1j}-B_{ij}^{(4)}{N}_{ij}}}\\ & +\underset{\text{shattering}}{\underbrace{\frac{1}{2}\sum _{i_1=1}^{N_{\mathrm{m}}}\sum _{j_1=1}^{N_{\mathrm{eg}}}\sum _{i_2=1}^{N_{\mathrm{m}}}\sum _{j_2=1}^{N_{\mathrm{eg}}}{C}_{i_1{j}_1{i}_2{j}_2}^{ij}{L}_{i_1{i}_2}{N}_{i_1{j}_1}{N}_{i_2{j}_2}-N_{ij}\sum _{i_1=1}^{N_{\mathrm{m}}}\sum _{j_1=1}^{N_{\mathrm{eg}}}{L}_{i_{1}i}{N}_{i_1{j}_1}}}\end{array}$$

$$L_{i_1{i}_2}=\pi (a_{i_1}+a_{i_2}{)}^2{v}_{i_1{i}_2}^{\mathrm{col}}{F}_{\mathrm{C}}$$

$$F_{\mathrm{C}}=1-\frac{2Z_{i_1}{Z}_{i_2}{e}^2}{m_{i_1{i}_2}(a_{i_1}+a_{i_2})(v_{i_1{i}_2}^{\mathrm{col}}{)}^2},$$

$$A_{ij}^{(1,2,3,4)}=\frac{{\epsilon }_{ij}^{(1,2,3,4)}}{E_{\mathrm{gb}}^{j+1}-E_{\mathrm{gb}}^j}.$$

$$B_{ij}^{(1,2,3,4)}=\frac{{\mu }_{ij}^{(1,2,3,4)}}{m_{\mathrm{b}}^{i+1}-m_{\mathrm{b}}^i}.$$

$$Y_{\mathrm{diss}}^{\mathrm{CC}}(a_i,Z_i)=\underset{0}{\overset{\infty }{\int }}{Y}_{\mathrm{diss}}^{\mathrm{CC},\mathrm{sp}}(a_i,E,Z_i)p(E)dE.$$

$${\mu }_{ij}^{(1)}=\frac{{\mu }_{\mathrm{d}}}{N_{\mathrm{A}}}{Y}_{\mathrm{diss}}^{\mathrm{CC}}(a_i,Z_i)\underset{0}{\overset{\infty }{\int }}{Q}_{\mathrm{abs}}(a_i,E,Z_i)a_i^2\frac{F(E)}{E}dE,$$

$$v_{i_1{i}_2}^{\mathrm{col}}=({\sigma }_{\mathrm{d},i_1}^2+{\sigma }_{\mathrm{d},i_2}^2{)}^{1/2},$$

$$\begin{array}{rl}{\epsilon }_{ij}^{(2)}& =\frac{\delta }{N_{\mathrm{atoms}}}[R_{\mathrm{s}}^{\mathrm{H}}(\mathrm{H})+R_{\mathrm{s}}^{\mathrm{He}}(\mathrm{H})+R_{\mathrm{s}}^{\mathrm{C}}(\mathrm{H})\\ & +R_{\mathrm{s},\mathrm{T}}^{\mathrm{H}}(\mathrm{H})+R_{\mathrm{s},\mathrm{T}}^{\mathrm{He}}(\mathrm{H})+R_{\mathrm{s},\mathrm{T}}^{\mathrm{C}}(\mathrm{H})+R_{\mathrm{s},\mathrm{T}}^{\mathrm{elec}}(\mathrm{H})].\end{array}$$

$$\begin{array}{rl}{\mu }_{ij}^{(2)}& =\frac{{\mu }_{\mathrm{d}}}{N_{\mathrm{A}}}[R_{\mathrm{s}}^{\mathrm{H}}(\mathrm{C})+R_{\mathrm{s}}^{\mathrm{He}}(\mathrm{C})+R_{\mathrm{s}}^{\mathrm{C}}(\mathrm{C})\\ & +R_{\mathrm{s},\mathrm{T}}^{\mathrm{H}}(\mathrm{C})+R_{\mathrm{s},\mathrm{T}}^{\mathrm{He}}(\mathrm{C})+R_{\mathrm{s},\mathrm{T}}^{\mathrm{C}}(\mathrm{C})+R_{\mathrm{s},\mathrm{T}}^{\mathrm{elec}}(\mathrm{C})].\end{array}$$

$$\begin{array}{rl}{\epsilon }_{ij}^{(3)}& =\frac{\delta }{N_{\mathrm{atoms}}}[R_{\mathrm{n}}^{\mathrm{H}}(\mathrm{H})+R_{\mathrm{n}}^{\mathrm{He}}(\mathrm{H})+R_{\mathrm{n}}^{\mathrm{C}}(\mathrm{H})\\ & +R_{\mathrm{n},\mathrm{T}}^{\mathrm{H}}(\mathrm{H})+R_{\mathrm{n},\mathrm{T}}^{\mathrm{He}}(\mathrm{H})+R_{\mathrm{n},\mathrm{T}}^{\mathrm{C}}(\mathrm{H})].\end{array}$$

$$\begin{array}{rl}{\mu }_{ij}^{(3)}& =\frac{{\mu }_{\mathrm{d}}}{N_{\mathrm{A}}}[R_{\mathrm{n}}^{\mathrm{H}}(\mathrm{C})+R_{\mathrm{n}}^{\mathrm{He}}(\mathrm{C})+R_{\mathrm{n}}^{\mathrm{C}}(\mathrm{C})\\ & +R_{\mathrm{n},\mathrm{T}}^{\mathrm{H}}(\mathrm{C})+R_{\mathrm{n},\mathrm{T}}^{\mathrm{He}}(\mathrm{C})+R_{\mathrm{n},\mathrm{T}}^{\mathrm{C}}(\mathrm{C})].\end{array}$$

$$\begin{array}{rl}{\epsilon }_{ij}^{(4)}& =\frac{\delta }{N_{\mathrm{atoms}}}[R_{\mathrm{n}}^{\mathrm{H}}(\mathrm{H})+R_{\mathrm{n}}^{\mathrm{He}}(\mathrm{H})+R_{\mathrm{n}}^{\mathrm{C}}(\mathrm{H})\\ & +R_{\mathrm{n},\mathrm{T}}^{\mathrm{H}}(\mathrm{H})+R_{\mathrm{n},\mathrm{T}}^{\mathrm{He}}(\mathrm{H})+R_{\mathrm{n},\mathrm{T}}^{\mathrm{C}}(\mathrm{H})].\end{array}$$

$$\begin{array}{rl}{\mu }_{ij}^{(4)}& =\frac{{\mu }_{\mathrm{d}}}{N_{\mathrm{A}}}\left[R_{\mathrm{big}}^{\mathrm{H}}(\mathrm{C})+R_{\mathrm{big}}^{\mathrm{He}}(\mathrm{C})+R_{\mathrm{big},\mathrm{T}}^{\mathrm{H}}(\mathrm{C})+R_{\mathrm{big},\mathrm{T}}^{\mathrm{He}}(\mathrm{C})\right].\end{array}$$

$${\epsilon }_{ij}^{(1)}=\delta Y_{\mathrm{diss}}^{\mathrm{CH}}{\sigma }_{\mathrm{loss}}^{\mathrm{CH}}\underset{10\mathrm{eV}}{\overset{13.6\mathrm{eV}}{\int }}{Q}_{\mathrm{abs}}(a_i,E_{\mathrm{g}}^j,Z_i,E)\frac{F(E)}{E}dE,$$

$$Y_{\mathrm{diss}}^{\mathrm{CC},\mathrm{sp}}(a_i,E,Z_i)=\frac{2k_{{\mathrm{C}}_2{\mathrm{H}}_2}[1-Y_{\mathrm{ion}}(a_i,E,Z_i)]}{k_{\mathrm{diss}}(a_i,E)+k_{\mathrm{IR}}(a_i,E,Z_i)},$$

Full text

††pubyear: 2018††pagerange: Shiva: the dust destruction model–B

1 Introduction

A large number of infrared (IR) observations available to the date indicate that dust evolves in the interstellar medium (ISM). This evolution is related to both growth and destruction. Many studies have already been performed to show how these processes manifest themselves in various astrophysical environments. Dust grains are mostly destroyed in extreme conditions with enhanced ultraviolet (UV) radiation field (RF), shock waves, hot gas, or turbulent motions. At these conditions grains can be eroded by high energy photons (Jochims et al., 1994; Allain et al., 1996a, b), via sputtering by fast ions and electrons (Draine & Salpeter, 1979; Tielens et al., 1994), or by shattering in grain-grain collisions (Jones et al., 1996; Hirashita & Yan, 2009). Rates of these processes depend on the dust material. While it is generally assumed that interstellar dust grains are most likely silicate and carbonaceous (Mathis et al., 1977; Desert et al., 1990; Li & Draine, 2001; Zubko et al., 2004; Jones et al., 2013), the specific nature of carbonaceous material in cosmic dust is still debated. In some models, which successfully explain observed dust extinction/IR emission properties, it is assumed that big carbonaceous grains are graphite-like and small carbonaceous grains are PAH-like (Draine & Li, 2007; Siebenmorgen et al., 2014). In other models, it is assumed that carbon dust grains are hydrogenated amorphous carbons (Compiègne et al., 2011; Jones et al., 2013), with PAHs being their special case.

While non-evolving dust models are still successfully used to interpret available observations, various destructive processes for specific interstellar conditions and dust materials are also considered on a case-by-case basis. Given an ever increasing volume of observational data, it would be convenient to have a tool for a quick theoretical estimation of evolutionary time-scales and an analysis of a varying dust size distribution. With this in mind, we present a tool Shiva that combines primary destructive processes into a single model and allows calculating the dust evolution under various conditions. The Shiva tool can be used as a stand-alone feature or be included as a module into a more detailed hydrodynamical model, like the MARION model of an expanding Hii region presented in the works of Akimkin et al. (2015, 2017). Shiva shares many elements with the recently published THEMIS model (Jones et al., 2017), however, in addition to processes presented in THEMIS (photo-aromatisation, sputtering), the Shiva model also accounts for photo-destruction via carbon atom loss, aromatisation by ions and electrons, and shattering processes. Our numerical tool is publicly available111http://www.inasan.rssi.ru/~khramtsova/SHIVA.html.

The paper is organised as follows. In Sect. 2 we describe the main elements of the dust evolution model and present equations being solved. In Sect. 3 we demonstrate some results obtained with the Shiva tool. Finally, conclusions are given in Sect. 4. In Appendix A and B we describe rate coefficients for photo-destruction and sputtering processes in more detail.

2 The model

A model of cosmic dust evolution implemented in Shiva includes photo-processing, i.e. hydrogen atom loss (aromatisation) and carbon atom loss, sputtering by ions and electrons, and shattering in grain-grain collisions. The model can be applied in a wide range of astrophysical objects except for those, where accretion/coagulation processes play an important role (i.e. molecular clouds, protoplanetary discs). The main concepts of the model were described in the works of Murga et al. (2016a, hereafter Paper I) and Murga et al. (2016b). In this work we present an advanced version of the model which now includes aromatisation by ions and treatment of grain charge, employs a dedicated method for solution of stiff kinetic equations, and is publicly available. Below we recapitulate the features of the model, paying more attention to its modifications. Processes of sputtering and shattering are treated analogously to Paper I, so details can be found there. For reference, expressions for sputtering rates are given in Appendix A.

The key element of the model is consideration of aromatisation process of carbonaceous grains, that is, restructuring of dust grains from the aliphatic state to the aromatic state as a result of dehydrogenation. Alteration of the material structure leads to a change of its band gap energy, $E_{\rm g}$, which decreases as a grain becomes more aromatic. Other considered processes change grain mass. Thus, each grain can be described by two parameters, its mass $m$ (or an equivalent radius $a$, which we calculate assuming a compact grain with solid density of 1.5 g cm*-3* as appropriate for HAC material) and the band gap energy $E_{\rm g}$. These parameters evolve under the influence of external factors. We divide the mass interval [$m_{\rm{min}}$, $m_{\rm max}$] into $N_{\rm m}$ bins with $N_{\rm m}+1$ borders, $m_{\rm b}^{i}$. The band gap interval [$E_{\rm g}^{\rm min}$, $E_{\rm g}^{\rm max}$] is divided into $N_{\rm eg}$ bins with $N_{\rm eg}+1$ borders, $E_{\rm gb}^{j}$. Thus, the total number of bins is $N_{\rm m}\times N_{\rm eg}$. Each bin can be characterised by the grain number density $N_{ij}$, cm*-3*, which is the number of dust grains within the bin, so that the total dust number density is $\sum_{i=1}^{N_{\rm m}}\sum_{j=1}^{N_{\rm eg}}N_{ij}$. We assume that all dust grains are initially in the same aromatisation state. If the initial mass distribution $dn/dm$ is known, $N_{ij}$ can be expressed as $dn/dm(m_{\rm b}^{i+1}-m_{\rm b}^{i})$. The evolution of $N_{ij}$ can be written in a form of the discrete Smoluchowski equation (Smoluchowski, 1916) that is often used for the description of coagulation process (Dullemond & Dominik, 2005; Akimkin, 2015):

{strip}

[TABLE]

Here $A_{ij+1}^{(1,2,3,4)}$ are rate coefficients that describe a hydrogen atom loss and a corresponding increase in aromatisation in the direction from a $(j+1)$th bin to a neighbouring $j$th bin by $E_{\rm g}$ (or from $j$th to $(j-1)$th bin in case of $A_{ij}^{(1,2,3,4)}$). We assume that only a loss of carbon atoms leads to the change in grain mass, so that the mass bin of a grain is not changed due to aromatisation process. $B_{i+1j}^{(1,2,3,4)}$ are rate coefficients that describe the carbon atom loss by photo-destruction and sputtering and a corresponding change in grain mass from an $(i+1)$th bin to a neighbouring $i$th bin (or from $i$th to $(i-1)$th bin in case of $B_{ij}^{(1,2,3,4)}$). We assume that the loss of carbon atoms does not change the grain aromatisation state. $C_{i_{1}j_{1}i_{2}j_{2}}$ is a ratio between the mass transferred to $i$th bin due to shattering and the average mass of the bin (see other details in Paper I), and $L_{{i_{1}}{i_{2}}}$ is the shattering kernel:

[TABLE]

with $v^{\rm col}_{{i_{1}}{i_{2}}}$ being the velocity of a grain-grain collision and $F_{\rm C}$ being the Coulomb factor, which is estimated as

[TABLE]

where $Z_{i_{1}}$ and $Z_{i_{2}}$ are charge numbers of grains, $e$ is electron charge, $m_{i_{1}i_{2}}$ is a reduced mass of colliding grains. The value of $F_{\rm C}$ is set to zero if the above expression gives negative values.

The rate coefficients of aromatisation are calculated as

[TABLE]

Boundary values of $A^{(1,2,3,4)}_{ij}$ for $j=1$ and $j=N_{\rm eg}+1$ are set to zero.

The rate coefficients for carbon atom loss are calculated analogously:

[TABLE]

Boundary values of $B^{(1,2,3,4)}_{ij}$ for $i=N_{\rm m}$ are set to zero, but we allow grains to leave the bin with $i=1$. We define $\varepsilon^{(1,2,3,4)}_{ij}$ and $\mu^{(1,2,3,4)}_{ij}$ in the following sections.

The system of ordinary differential equations (ODE) (1) is stiff and should be solved using an appropriate method. We chose the solver CVODE222https://computation.llnl.gov/projects/sundials/cvode from the SUNDIALS software package (Hindmarsh et al., 2005) that allows solving stiff systems of ODE.

The Shiva code is implemented as a C*++* computer code. Detailed instructions for users can be found in the web link specified above.

2.1 Photo-processing

We consider photo-destruction processes that include the loss of both carbon and hydrogen atoms. In general, the absorption of a UV photon by a grain can lead to several outcomes: ionization (charging), IR photon emission, dissociation of C–H or C–C bonds. These processes occur with probabilities $Y_{\rm ion}$, $Y_{\rm IR}$, $Y_{\rm diss}^{\rm CH}$, $Y_{\rm diss}^{\rm CC}$, respectively, which sum up to unity. If we assume that photon energy is spent to dissociation, then a C–H bond is more likely to be destroyed as it is less strong than a C–C bond. Thus, a dust grain first loses its hydrogen atoms and becomes partially or fully dehydrogenated, so that $Y_{\rm diss}^{\rm CC}=0$ at this stage. Simultaneously the structure of a grain material changes from aliphatic to aromatic333This applies only to HAC material. as is described in detail in the works of Jones (2012a, b, c). In these papers it is mentioned that the maximum band gap energy value is 2.67 eV, and the minimum value depends on a grain radius, but we set this value equal to 0.1 eV for all dust grains as it was adopted for aromatic material in the work of Jones et al. (2013). The rate of change of the band gap energy $\varepsilon^{(1)}_{ij}$ due to aromatisation by photons can be found in Appendix B.

When a grain is dehydrogenated down to $E_{\rm g}=0.1$ eV (or is a PAH particle), another destruction scheme with carbon atoms loss starts to operate, and at this stage we assume $Y_{\rm diss}^{\rm CH}=0$. This process is relevant only for small particles, so we limit its action to grains with the number of carbon atoms less than 1000. The loss of carbon atoms requires more energy and apparently is only possible near a strong UV radiation source (Jochims et al., 1994; Allain et al., 1996a; Zhen et al., 2015). In the model we assume that subsequent dissociation proceeds analogously to PAHs as described in the work of Allain et al. (1996a)444It should be mentioned that small HAC particles described by Jones et al. (2013) evolve with time, and are very likely aromatic in the ISM. However, even small aromatic grains are not identical to PAHs, as they are volumetric in a general case, though in some specific cases they can be PAHs (Jones, 1990). But we assume that dissociation properties (like the binding energy) for small grains are always the same as for PAHs.. Expression for dissociation probability of dehydrogenated grains $Y_{\rm diss}^{\rm CC,sp}$ can be found in Appendix B.

In the model we consider the multi-photon heating mechanism described in the work of Guhathakurta & Draine (1989). It implies that a dust grain can accumulate energies of several absorbed photons before this energy is spent on the grain IR emission. The internal energy of a grain can grow several times higher in this case than in a single-photon approach, and, thus, dissociation occurs more effectively. We calculate the energy history for a grain with a radius $a$ immersed in the radiation field with the flux $F(E)$ and then derive the probability distribution of its internal energy $p(E)$ using the approach described in the work of Pavlyuchenkov et al. (2012). Importance of the multi-photon approach for small grain destruction has been noted in the work of Visser et al. (2007). Thus, the dissociation probability $Y_{\rm diss}^{\rm CC,sp}$ should be integrated over all possible internal energies:

[TABLE]

The rate of a grain mass decrease as a result of photo-destruction is expressed as

[TABLE]

where $Q_{\rm abs}$ is a grain absorption efficiency coefficient, $F(E)$ is the incident radiation flux, erg cm*-2* s*-1* eV*-1*, $E$ is photon energy, $\mu_{\rm d}$, g mol*-1*, is molar weight of the grain material, $\mu_{\rm d}=12-11X_{\rm H}$, $X_{\rm H}$ is hydrogen atom fraction, and $N_{\rm A}$ is Avogadro constant. The ratio of $\mu_{\rm d}$ and $N_{\rm A}$ is used to convert the rate from the number of atoms per second to gram per second. In the current study we assume that the radiation field is represented by the scaled interstellar radiation field (ISRF), taken from Mathis et al. (1983), which we refer to as MMP83. The scale factor is $U$, so that the local radiation flux is $U\times{\rm MMP83}$. Other spectra $F(E)$ can be used as well. The adopted optical properties (in particular, $Q_{\rm abs}$) are discussed in subsection 2.2.

Photo-destruction of carbonaceous grains is in fact a blurred question. Partly this is related to the lack of knowledge of which types of grains are encountered in space. Since the unidentified infrared (UIR) bands had been discovered (e.g. Gillett et al., 1973), further investigations have constrained possible assumptions, but the final clarification is yet to come. While it is more or less reliably established that small carbonaceous grains should contain at least some aromatic rings to explain the UIR bands (Allamandola et al., 1985), their actual structure, size, chemical formulae are debated. Several candidates that are considered to date are PAHs (Leger & Puget, 1984; Allamandola et al., 1985), HAC grains (Duley, 1985; Jones et al., 2013), mixed aromatic-aliphatic organic nanoparticles (Kwok & Zhang, 2011). Even if we adopt the most studied candidates, PAHs, as typical representatives of the smallest grains, their exact shape is also far from being fully understood. While some studies indicate that PAHs should be compact, pericondensed and consist of no less than 30–40 atoms (Allamandola et al., 1989; Jochims et al., 1994; Langhoff, 1996), other works offer different conclusions (Ekern et al., 1998), stressing that the photo-destruction rate strongly depends on many parameters and not only on grain size and shape. It has been shown that biphelynene (C12H8) is a stable molecule in spite of consisting of only two aromatic rings. Coronene (C24H12) loses hydrogen atoms under the influence of UV radiation, but catacondensed chrysene or tetracene (both C18H12) lose one or two hydrogen atoms and acetylenic group and then become photostable. So, apparently, there is no universal approach to the calculation of photo-destruction rates and products for any PAHs, and it is necessary to consider each molecule individually.

The situation is not much clearer with the photo-destruction of HAC grains. Many experiments prove that their dehydrogenation occurs under UV irradiation (Mennella et al., 2001). Rates and products of HAC photo-destruction were studied in recent laboratory experiments (Alata et al., 2014; Alata et al., 2015), but perhaps these rates cannot be applied directly to the ISM. Interstellar HAC grains are probably more aromatic, at least on the surface, and consequently are more stable. So in the absence of a general theory for the photo-destruction of HAC grains, at the first destruction step we apply experimental data, if the grains are predominantly aliphatic, and the RRK (Rice-Ramsperger-Kassel) theory, adopted to PAHs by Allain et al. (1996a), if the grains are aromatic.

2.2 Grain charge

Supplementing our previous work, we include a charge of a grain $Z$ in the calculations of the destruction rates. It was shown that charged and neutral small grains have different absorption properties (Bohren & Hunt, 1977; Kocifaj & Klačka, 2012; Kocifaj et al., 2012). Moreover, the fate of an absorbed photon (ionization, IR photon emission, dissociation of C–H and C–C bonds) also depends on a grain charge. In this paper we study how the charge influences branching over these possible channels. We consider only variations in $Q_{\rm abs}$ in the mid-IR range and not in the UV range, so the main factor is a grain cooling rate, not a grain heating rate.

The charge state of dust grains may affect the efficiency of the processes considered in this work. Quite naturally, even large charge numbers are not important for shattering as at high collision velocities needed for the fragmentation the Coulomb factor is close to unity. Sputtering of big grains is also insensitive to their charge. However, charged and neutral small grains do respond differently to photo-destruction and sputtering. First, small ionized grains have different optical properties that affect their IR-emission rate. Under non-equilibrium conditions, when single-photon heating is important, IR opacities define the rate, with which the grain sheds the extra energy, and destruction probability. Second, it is more difficult to ionize grains that are already positively charged. This is why the dissociation yields are different for ionized and neutral grains. Here we consider how the grain charge influences the optical properties. Comparison of various rates for ionized and neutral grains is presented in Sect. 3.1.

We adopt optical properties for neutral dust grains presented in the works of Jones (2012a, b, c) with some modifications concerning the bands around $7-8\,\mu$m described in the work of Jones et al. (2013). The full list of mid-infrared bands, corresponding to both aliphatic and aromatic bonds, can be found in Table C.1 of the work of Jones et al. (2013). Modified properties for dehydrogenated HAC were calibrated to the properties used in the DustEm model (Compiègne et al., 2011).

To compare efficiencies of destructive processes for charged and neutral grains one needs to know optical properties for charged grains as well. Unfortunately, studies of optical properties of ionized HAC grains are apparently lacking, while the properties of PAH cations are well investigated. In this situation the only way to proceed is to assume that ionization alters optical properties of HAC grains in the same way as it changes optical properties of PAHs. Perhaps, this approach is justified as most small dust grains in the ISM should be aromatised (Jones et al., 2013), and PAHs can be considered as a special case of small aromatic HACs. Moreover, we only need to know the properties of charged grains, when we calculate rates of those processes in which the smallest grains are involved. Thus, we describe the properties of such ionized grains as if they were PAHs, while the optical properties of bigger grains are assumed to be independent on the charge state. In this work we deal with conditions where the dust grains are either neutral or positively charged, so we consider only cations. Also, we do not take into account multiple ionized states in calculations of optical properties because in this case we would have to make too many assumptions. So, talking about optical properties of charged dust grains, we only mean singly ionized cations.

It has been shown (DeFrees et al., 1993; Schutte et al., 1993; Langhoff, 1996; Boersma et al., 2011) that charge has a strong influence on aromatic infrared bands. Specifically, intensities of C–H bands (around $3.3\,\mu$m) are reduced (Langhoff, 1996), while intensities of C–C stretching and in-plane C–H bending modes ($6-9\,\mu$m) are enhanced for ionized PAHs (Szczepanski & Vala, 1993; Langhoff, 1996; Hudgins & Allamandola, 1999; Kim et al., 2001; Bauschlicher et al., 2009). The bands corresponding to C–H out-of-plane bending are also shifted due to the charge (Hudgins & Allamandola, 1999; Boersma et al., 2013; Shannon et al., 2016). Based on these and other works, the optical properties for both neutral and ionized PAHs were reproduced in the PAH model presented in the works of Li & Draine (2001); Draine & Li (2007). We adopt the band variations related to the ionized state analogously to these works. In particular, we make the following modifications:

• •

The integrated strength of C–H stretching bands (3.25, 3.28, 3.31, 3.32, 3.35, 3.38, 3.42, 3.45, 3.47, 3.48, $3.51\,\mu$m) for ionized grains is smaller by a factor of 4.4 compared to neutral grains.

• •

The integrated strength of aromatic C–C stretching bands (6.10, 6.25, 6.67, 7.53, 7.69, 7.85, $8.60\,\mu$m) for ionized grains is larger by a factor of 8 compared to neutral grains.

Other modes are left unchanged.

In Fig. 1 we demonstrate how the adopted optical properties differ from the original properties from the work of Jones (2012c) for the grain size of 5 Å and the band gap energies 0.1 eV (left panel) and 2.67 eV (right panel). These $E_{\rm g}$ values correspond to aromatic (dehydrogenated) and aliphatic (hydrogenated) states, respectively. We focus on the wavelength range from 2 to 15 $\mu$m assuming that properties are identical in other spectral bands. The original data from Jones (2012c) are shown with blue lines. Green lines correspond to a newer version of optical data from Jones et al. (2013), where some bands around $8\,\mu$m have been added. We utilise these data for neutral grains. The adopted $Q_{\rm abs}$ for charged grains are shown with red lines.

In the presented version of the model a grain charge is treated as an input parameter for each bin and should be specified from some separate considerations.

3 Results

Here we present some results produced by the Shiva tool. We show how the grain charge influences grain destruction rates in subsection 3.1. Then we apply the model to some specific astrophysical conditions. In the general ISM both the ISRF and relative motions of gas and dust only slightly affect the distributions of sizes and aromatisation states of dust grains. But there are more extreme objects, like photo-dissociation regions (PDRs), Hii regions, and supernova remnants, where dust is exposed to more powerful destructive factors. For clarity, we consider separately the extreme medium with the strong radiation field and the extreme medium with high velocities and temperatures. We follow the evolution of dust grains in these two environments in subsections 3.2 and 3.3.

3.1 Rate comparison for charged and neutral grains

It was mentioned above that charge plays a role in defining the cooling rate of a grain and its ionization probability. Both these quantities are needed to estimate the photo-destruction rate. Here we consider how they differ for neutral and charged grains and how this difference influences the photo-destruction rate.

Fig. 2 (left) illustrates the dependence of $k_{\rm IR}$ on the photon energy for ionized and neutral grains with radii of 3.3 and 5 Å. We further assume that absorption efficiency $Q_{\rm abs}$ for charged HAC grains is determined by the very presence of non-zero charge, but does not depend on its value or sign. So we compare the IR emission rates ($k_{\rm IR}$) only for HAC0 and HAC*+*. For both sizes $k_{\rm IR}$ is several times higher for ionized grains, and this difference increases with the photon energy. This implies that a neutral grain is less stable to photo-destruction than a charged grain, because its IR relaxation rate is lower.

In Fig. 2 (right) we show the dependence of the ionization yield $Y_{\rm ion}$ computed using an approach from the work of Weingartner & Draine (2001a) on the photon energy for a 3.3 Å grain with charge numbers of $Z=0,1,2,3$. Obviously, $Y_{\rm ion}$ decreases with increasing charge. In other words, probability to ionize the grain drops with each subsequent ionization. Therefore, more energy is available for destruction or emission. In this respect, an ionized grain is less resistant to destruction than a neutral one, contrary to what is said in the previous paragraph.

Relying upon calculations of photo-destruction rates for charged and neutral grains and taking into account the combined dependence of $Y_{\rm ion}$ and $k_{\rm IR}$ on the charge state, we conclude that the net effect makes charged grains less stable to photo-destruction than neutral ones, which is further demonstrated in Fig. 3. For a 3.3 Å grain the rate is increased by a factor of $\sim 1.2$, while for the 5 Å grains the factor can be as high as $\sim 2.5$, depending on the charge value. These ratios are independent on the RF intensity. The larger are dust particles, the greater is the difference between the photo-destruction rates for ionized and neutral grains, but the rates themselves are too small to influence the relevant evolutionary time-scales. The analogous conclusions were made by Allain et al. (1996b), where they investigated the stability of PAHs. We obtain a similar difference for the sputtering rate due to elastic interaction with fast moving ions/electrons.

3.2 Evolution of dust in the medium with enhanced UV radiation field

3.2.1 Time-scales of photo-processing

Before we proceed to further results, let us introduce two time-scales, which we consider as characteristic evolutionary times for various processes. We define a time-scale $\tau$ as a time needed to reach an $e$-fold change in the number density in an $i$th size bin summed over all aromatisation bins ($\sum\limits_{j}N_{ij}(t)$) relative to the initial density $\sum\limits_{j}N_{ij}(0)$. This time-scale is related to the overall evolution of the grain size distribution. Its value can be both positive and negative. If the number density in a bin decreases with time, $\tau$ is positive and represents the time needed to reduce the number density by a factor of $e$. Conversely, if the number density increases with time, then $\tau$ is negative. Another time-scale, $\tau_{\rm a}$, is the time during which the number density $N_{ij_{\rm min}}(t)$ in the $i$th mass bin and the $(j_{\min})$th aromatisation bin (the highest aromatisation state) becomes $\approx e$ times larger than the number density $N_{ij_{\rm max}}(t)$ in the same mass bin $i$, but in the aromatisation bin $j_{\rm max}$ (the lowest aromatisation state). In this study, we assume that initially all the aromatisation bins, except for bins with $E_{\rm g}=E_{\rm g}^{\max}$, are empty. This time-scale roughly characterises the aromatisation rate in a given mass bin.

When the radiation prevails over other factors, we can consider only the photo-destruction processes (hydrogen and carbon loss). In Fig. 4 the dependence of the photo-destruction and aromatisation time-scales on the scaled radiation field, $U$, is shown. The red line shows the time-scale of a 5 Å-grain destruction due to carbon loss under the influence of the UV field. We see that the destruction time-scale for such grains drops below $10^{6}$ yr only when $U$ exceeds $10^{3}$. For example, in the Orion Bar, where $U$ is estimated to be $10^{4}-10^{5}$ (Tielens & Hollenbach, 1985), PAH-like grain should be destroyed on a time-scale of about $10^{3}$ yr, which seems to be relevant for this PDR (Kassis et al., 2006). In another PDR, Horsehead Nebula, the value of $U$ is smaller, about $10^{2}$ (Compiègne et al., 2007). In this case the PAH-like grains (or, at least, their carbon skeletons) should survive for a longer time, up to several Myr. The rate of hydrogen loss is greater, and the full dehydrogenation of aromatic species can lead to formation of fullerenes (Berné & Tielens, 2012), but in this work the corresponding evolutionary path is not considered. Also, we do not consider dehydrogenation of PAHs and HAC grains with $E_{\rm g}<0.1$ eV as was pointed above.

Blue, green, and magenta lines in Fig. 4 illustrate the aromatisation time-scales of 5, 50 and 1000 Å-grains, correspondingly. It can be seen that the time-scale for each grain is inversely proportional to the intensity of the radiation field. Among the considered sizes, aromatisation is most effective for 50 Å-grains. Analogous results were obtained in the work of Jones et al. (2014), which we rely upon in our calculations. This result comes from Eq. 15. Absorption efficiency grows with size while other values are constant until the size reaches 200 Å, as small grains are in the Rayleigh regime, so the aromatisation rate increases with size. At $a=200$ Å, $Y_{\rm diss}^{\rm CH}$ becomes inversely proportional to the grain size, while the absorption efficiency is independent on the size, because large grains are in the geometrical regime, so the time-scale of aromatisation for large grains (in our case 1000 Å-grains) turns out to be smaller than that for small grains (5 Å).

3.2.2 Evolution of grains and infrared spectra in PDRs

In Fig. 5 we present results of calculations that include only photo-processes for HAC particles with the initial J13 size distribution (Jones et al., 2013), assuming that all grains are initially hydrogenated. We took the evolution time interval of 1 Myr which is a typical lifetime of PDRs. The considered intensities are varied from 1 to 105 MMP83. We take into account the charge of dust grains, which is calculated using approach from the work of Weingartner & Draine (2001a), so that the average charge of small grains is about 1 for $U>10^{3}$ and 0 for lower RF intensities. Evolution of big grains is insensitive to the charge as their optical properties are the same for charged and neutral states. In Fig. 5 (left) blue lines show the dependence of the fractional dust loss from a 5 Å size bin on the radiation field intensity. At $U<10^{2}$ grains of this size are not destroyed at all for 1 Myr. At higher intensities the fraction of destroyed dust, $f_{\rm destr}$, starts to grow, reaching 100% at $U=10^{4}$. Aromatisation process is more efficient. The dependencies of fractions of fully aromatised dust grains with radii of 5, 50 and 1000 Å on the radiation field intensity are shown by green, red, and cyan lines, correspondingly. The 50 Å-grains are aromatised independently on the value of $U$ (in the considered range), other grains are fully aromatised at $U>10$.

In Fig. 5 (right) the infrared dust emission spectra after 1 Myr of evolution are presented. Spectra were computed using DustEm555https://www.ias.u-psud.fr/DUSTEM/index.html tool (Compiègne et al., 2011). The optical properties for the HAC grains with different $E_{\rm g}$ as well as grain size distributions were taken from this work, while heat capacities are assumed to be the same as those for HAC grains in J13 model. The infrared bands are most sensitive to environmental conditions. At $U<10^{2}$, the band near 3 $\mu$m corresponds to vibrations of aliphatic bonds, at $U\geqslant 10^{2}$ the contribution of aromatic bonds starts to dominate, and in this case the band at 3.3 $\mu$m is more intense than that at 3.4 $\mu$m. At intensities $U\gtrsim 5\cdot 10^{3}$ the 3 $\mu$m band is barely seen, because the smallest grains responsible for this band are destroyed.

Bands around 6–8 $\mu$m attain features more typical for aromatic grains at lower intensities than bands around 3 $\mu$m. At $U=5$, the spectrum in the 6–8 $\mu$m range corresponds to the mix of aliphatic and aromatic grains but at higher intensities only bands characteristic of aromatic grains remain. More rapid aromatisation of bands in the 6–8 $\mu$m range is related to faster aromatisation of bigger grains, which are responsible for the emission at these wavelengths. Moreover, these aromatic bands do not completely disappear even at the highest considered intensities, because the grains are too big to be destroyed by the UV radiation in the current model.

We illustrate the evolution of infrared spectra at different RF intensities in Fig. 6. At $U=1$, the spectrum near 3 $\mu$m remains almost invariable during the considered time span, while at longer wavelengths (around 10–20 $\mu$m) the spectrum changes due to aromatisation of medium-sized grains. At higher irradiation intensities the spectra at 3 $\mu$m also change as the smallest grains become aromatic, too. At $U=10^{3}$, the transition from aliphatic-dominant spectrum to aromatic-dominant spectrum occurs between $10^{4}$ and $10^{5}$ yr, when the band at 3.3 $\mu$m becomes more intense than the band at 3.4 $\mu$m. At $U=2\cdot 10^{4}$ this transition occurs faster by approximately 2 orders of magnitude, and the bands near 3 $\mu$m disappear during further evolution due to destruction of corresponding grains. We chose the specific value of $U=2\cdot 10^{4}$ as this intensity is representative for the Orion Bar (Tielens & Hollenbach, 1985; Goicoechea et al., 2015). Bands at both 3.3 and 3.4 $\mu$m are observed in the Bar, although they are not prominent. On the other hand, the bands near 11.2 and 12.6 $\mu$m are strong in the Bar, while, according to our calculations, these bands should not be so intense at this RF value. Thus, our model indicates that properties of amorphous carbons alone are not sufficient to explain emission around 11 $\mu$m in this region, and some other processes are important in addition to those already included in the model. In fact, at such intense radiation field amorphous carbonaceous grains likely fragment to slices forming PAHs or small hydrocarbons. The inconsistency between our model and the Orion Bar observations can be relaxed if we replace or amend HAC optical properties with the properties of genuine PAHs.

This is further illustrated in Fig. 6 (bottom), where we show evolving infrared spectra at $U=2\cdot 10^{4}$ for an initial size distribution from the work of Weingartner & Draine (2001b) (WD01) and the optical properties from the work of Draine & Li (2007) (DL07). This distribution explicitly contains PAHs with their corresponding optical properties that we use to compute both photo-rates and infrared spectra. Indeed, in this case we do see all the PAH infrared bands. Obviously, they also become less intense with time, but still exist after 1 Myr of evolution. This confirms that there should be an additional step in the evolution of HAC particles, when fragments of HAC have optical properties similar to that of planar aromatic macromolecules like PAHs. In the current version of the model we do not consider this step and the corresponding evolution of optical properties.

3.2.3 Photo-processing in IR bubbles

Another astrophysical question that can be answered with the Shiva model is related to non-uniform spatial dust emission in the so-called IR bubbles (Anderson et al., 2012). Most of these objects are believed to be expanding Hii regions around massive stars surrounded by gas-dust clouds. Their name stems from a ring-like structure of the dust emission in these objects. While big grains can be relatively easy swept by radiation pressure (Draine, 2011; Akimkin et al., 2015), small nano-size grains cannot be easily blown-out from the inner part of the ionized region. Some other factor is needed to explain the lack of small grains there. The radiation field in the star vicinity is very strong, so one may suggest that small grains are efficiently destroyed by UV photons. The Shiva model provides an easy way to check this hypothesis and show how IR spectra change during the Hii region evolution.

As a reference for IR bubble parameters, we use results of numerical modelling, pre-computed with the MARION code and presented in the works of Pavlyuchenkov et al. (2013) and Akimkin et al. (2015, 2017). To track the changes in the object physical structure, we choose three representative locations within the modelled object. One of this location approximately corresponds to the inner wall of the dense shell that envelops the ionized region. Specifically, it is chosen to be the location where the density of neutral atomic hydrogen reaches its maximum ($n_{\rm H}^{\max}$). The second location corresponds to the density of neutral hydrogen, which is about half of its maximum value ($0.5n_{\rm H}^{\max}$). This is roughly a transition from neutral hydrogen to ionized hydrogen. Finally, the third location is selected inside the Hii region, at the point, where density of ionized hydrogen has its maximum ($n_{\rm H+}^{\max}$). As the Hii region expands, these locations also move away from the star. At each time moment from 0 to about 300 kyr, we extract physical parameters for these locations, including the radiation field, and use them to simulate dust evolution with the Shiva code.

Radiation field intensities for two time moments, 5 and 150 kyr, at the considered locations are illustrated in Fig. 7. For comparison, the MMP83 radiation field is also shown. At 5 kyr the radiation field intensity at all points is very strong, 4–5 orders higher than MMP83. Also, at the point inside the Hii region there are photons with $\lambda<912$ Å, which are absent at other locations. At 150 kyr the radiation field intensity remains high only inside the Hii region, while the field at other points becomes relatively weak and comparable to MMP83.

We adopt the WD01 initial size distribution and assume that initially dust grains are fully hydrogenated HACs. In Fig. 8 the evolution of number density of dust grains in the size bins corresponding to radii of 4, 6, and 8 Å at three selected locations is presented. Note that we only show the density of dehydrogenated grains. While initially their density is zero (because all grains are hydrogenated), but all grains become fully dehydrogenated during first 5 kyr, and afterwards they can lose only carbon atoms. The smallest grains experience significant destruction at all the considered locations, but only in the inner part of the Hii region destruction of 4 Å-grains occurs over the entire calculation time, 300 kyr. Outside the inner part of the ionized region destruction practically stops after $\sim 20$ kyr. Due to the region expansion, the radiation field becomes too weak for photo-destruction there.

Grains with a radius of 6 Å are destroyed in the Hii region, but stay mostly intact in more remote locations (with half-maximum and maximum neutral hydrogen density). The 8 Å grains are not destroyed anywhere. This means that in spite of very strong radiation field around massive stars only the smallest grains can be photo-destroyed within Hii regions, at least, according to the evolutionary model implemented in Shiva. We use the RRK approach to estimate the dissociation rate, and this rate strongly depends on the number of carbon atoms in a grain. This causes large decrease in the dissociation rate, when the size of a grain varies only a little.

Let us now consider how this evolution of small grains influences infrared spectra for the whole grain ensemble, shown in Fig. 9 for the three locations described above. As initially all dust grains are hydrogenated, the spectra at $t=0$ shown with blue lines in all three panels demonstrate IR features corresponding to fully hydrogenated amorphous carbon grains. After only 5 kyr of the evolution the spectra change: the features of aliphatic bonds disappear, while IR features corresponding to aromatic bonds emerge in the spectrum. First they are quite strong, but later gradually weaken, especially at 3 $\mu$m, because this band arises from the smallest grains, which are destroyed over the entire region. Features at 6–8 $\mu$m only slightly stand out above the continuum spectrum of hot large grains in the right panel, i.e. in the ionized hydrogen region, as 6 Å grains are mostly destroyed there. At 150 kyr the spectra only change due to the radiation field change, but the number density of dust grains at the selected locations does not vary significantly since then.

Summarising, in a typical IR bubble small grains with radius of about 4 Å cannot survive. Larger grains with radius 6 Å can survive only close to the dense shell, where neutral hydrogen dominates and there are no photons with energy higher than 13.6 eV. Grains with radius of 8 Å can survive even inside the bubble, where the ionized hydrogen dominates. As no emission at 8 $\mu$m is observed in the central parts of real bubbles, there should be some other mechanism, responsible for the absence of 8 Å grains in the centres of infrared bubbles.

3.3 Evolution of grains in the medium with high gas/dust velocities and temperature

To demonstrate the evolution of dust in the medium with high velocities and temperatures we have performed calculations on a grid of velocities $v_{\rm gas}$ in the range of $[30...200]$ km s*-1* and temperatures $T_{\rm gas}$ in the range of $[10^{3}...10^{6}]$ K with a fixed value for a product $n_{\rm H}T_{\rm gas}=10^{4}$ cm*-3* K, taking into account that temperature and density can be coupled and anti-correlated. Abundances of helium and carbon are assumed to be $n_{\rm He}=0.1n_{\rm H}$ and $n_{\rm C}=10^{-4}n_{\rm H}$, respectively. We deliberately assume zero RF in these computations to highlight the effects of sputtering and shattering. The gas-dust relative velocity ($v_{\rm gd}$) is taken to be $3/4$ of the gas velocity, which is appropriate for adiabatic shock waves. We estimate the grain-grain collision velocity between bins $i_{1}$ and $i_{2}$ as

[TABLE]

where $\sigma_{\rm d}$ is grain velocity standard deviation, which is equal to the gas-dust relative velocity. It is also assumed that dust differential drift is negligible.

In Fig. 10 (top row) we demonstrate maps of grain destruction time-scales $\tau$ for various grain sizes. We stop the calculations, when the ratio $\sum\limits_{j}N_{ij}(0)$/$\sum\limits_{j}N_{ij}(t)$ is around $e$. In the simulations an exact $e$ value is never reached, and this may introduce some minor uncertainties in determining the time-scales, but they do not alter the overall picture. In the top and middle rows of this diagram the results are presented for amorphous carbon grains with the J13 size distribution, which are assumed to be fully hydrogenated initially. The grains with radii of 5, 50 and 1000 Å are chosen as typical representatives of the smallest nano-size grains, very small grains, and big grains, correspondingly.

The time-scale for each grain type is a result of the balance of the considered processes. In Fig. 10 (top row) for 5Å-grain the destruction time-scale gradually increases from high velocities and low temperatures (high density) to low velocities and high temperatures (low density) reflecting that the density plays a key role in determining the destruction efficiency. The same trend is seen for grains with radii of about 50 Å and 1000 Å. This implies that, despite the sputtering process being very efficient at high temperature, grains have large time-scales due to low gas density. The decreasing time-scale can only be noted for highest values of considered temperatures. Presence of various heterogeneities on the maps clearly shows how subtle the balance of evolutionary processes is and how strongly it depends on the adopted conditions.

Smoother changes are seen for aromatisation time-scales in Fig. 10 (middle row) as the only process that influences this time-scale is sputtering. The most efficient aromatisation occurs at the highest velocities and again it slows down with the decreasing gas density.

As the outcome of the evolution strongly depends on the abundance of small grains, it is obviously sensitive to the adopted initial size distribution. In the bottom row of Fig. 10 we show the same map of time-scales as in the top row, but for graphites and the WD01 initial size distribution.

In the case of the WD01 initial size distribution, the number density in a specific mass bin can both increase and decrease during the evolution, which we denote with the sign of the corresponding time-scales. Negative values indicate increasing number density, and this may indeed happen with smaller grains at the considered conditions. Specifically, the replenishment of 50 Å-grains due to shattering of bigger grains dominates over their destruction via sputtering, when $v_{\rm gas}\lesssim 150$ km s*-1* and $T_{\rm gas}\lesssim 5\cdot 10^{5}$ K. Under different conditions the destruction dominates, and the time-scales are positive. The smallest and the largest grains have only positive time-scales. Destruction time-scales for them always decrease with velocity and increase with temperature as in the adopted setup growing temperature implies decreasing density.

While grain destruction by energetic particles has already been considered in the literature, the novel feature of our model is that we also consider simultaneous changes in the grain structure. We demonstrate the efficiency of aromatisation by energetic ions and electrons, using the results of supernova remnant (SNR) shock simulations performed in the work of Nozawa et al. (2006). The sputtering process was treated in their calculations, so we simulate only aromatisation process adopting the evolution of grain sizes from their work. Calculations by Nozawa et al. (2006) are intended to reproduce conditions in the early Universe, but this is not a critical issue for our study as we only use their results as a reference for a sample computation.

Initial shock velocity is about 3000 km s*-1*, but the starting time of dust evolution calculation is 2000 yr in their model. By this time the shock velocity falls to $\sim 1000$ km s*-1*, temperature of the shock front is $\sim 6\cdot 10^{7}$ K, and gas number density is $\sim 3$ cm*-3*. As dust grains of different radii move differently, the conditions for them at each time also differ. Distance, over which a grain moves inside the shock, varies depending on the time-scale of its destruction. The smallest grains are destroyed rapidly, while large grains are not completely destroyed and gradually turn around in the direction of shock wave motion. So the conditions (relative gas-dust velocity, density, temperature) cover different ranges for different grain sizes. A dust grain with a radius of $a=8.6\cdot 10^{-7}$ cm is exposed to the following conditions: $v_{\rm gd}$ in $[1000\ldots 15]$ km s*-1*, $T_{\rm gas}$ in $[6\cdot 10^{7}\ldots 4\cdot 10^{7}]$ K and $n_{\rm H}$ in $[3-1]$ cm*-3*, where the ranges are defined by the lifetime of the grains. Grains of larger radii ($a>1.7\cdot 10^{-6}$ cm) survive in the shock, so ranges of $v_{\rm gd}$, $T_{\rm gas}$ and $n_{\rm H}$ for them are wider: $[1000\ldots 0]$ km s*-1*, $[6\cdot 10^{7}\ldots 1\cdot 10^{7}]$ K, $[3-0.1]$ cm*-3*, correspondingly.

In Fig. 11 (left) we show how various grains change their radii (solid lines) and band gap energies (dashed lines) as they traverse the shock. A time-scale of aromatisation is always smaller than a grain destruction time-scale. Grains with radii of about $8.6\cdot 10^{-7}$ cm and smaller disappear after the supernova shock wave passage. Thus, the aromatisation level of these grains is irrelevant. Larger grains survive in the shock wave though their radii are reduced. Grains smaller than $\approx 1.4\cdot 10^{-5}$ cm are fully aromatised by ions and electrons, while larger grains remain aromatised only partially.

In Fig. 11 (right) we present the infrared spectra of the modelled SNR shock region, integrated spatially to include grains of all sizes, which have different velocities and positions in the considered model. Initially the spectrum is that of fully hydrogenated HACs (black line). By the time of $3.5\cdot 10^{3}$ yr small grains are completely aromatised, and the mid-IR spectrum has the shape typical for aromatic grains. Later, at $4.8\cdot 10^{3}$ yr, number density of small grains becomes smaller, so the emission intensity in the wavelength range of 2–50 $\mu$m decreases almost by the order of magnitude. At even later times small dust grains are destroyed, and the emission in the mid-IR range is nearly absent. Only large grains contribute to the spectrum at longer wavelengths. Such computations can be useful for interpreting real observations of the IR spectrum evolution across the SNR shock, like those reported by Tappe et al. (2006) and Arendt et al. (2010). Specifically, comparing observed spectra and Shiva predictions, we can infer the time-scale of the supernova remnant expansion similarly to what has been done by Tappe et al. (2012). Also, it may turn out that the fading brightness of the Spitzer 3.6 $\mu$m band observed in the remnant of SN1987A (Arendt et al., 2016) is at least partially caused by grain restructuring.

The Shiva tool is well suited to provide quick preliminary answers to these questions, which will be useful for interpreting SOFIA and forthcoming JWST spectral observations. Other possible scenarios of Shiva usage may include studies of grain evolution in the intergalactic medium (Gjergo et al., 2018), establishing the origin for various correlations in gas and dust properties in star-forming regions (Smirnova et al., 2017; Smirnova & Wiebe, 2019), etc.

4 Summary

Dust evolution in the ISM is a hot topic in the modern astrophysics, and it is important even in those studies, which are not directly concerned with the dust itself. We present a physical model Shiva for simulating the dust evolution in warm and hot ISM phases. The model is implemented as a numerical tool that can be used for estimating the evolutionary time-scales for various grains at different conditions in the ISM and for studying dust size evolution. The main destructive processes included in the presented version of the model are photo-processing, shattering, and sputtering. The input model parameters are the radiation flux, hydrogen density (and also densities of helium and carbon as well as fractions of ionized H, He, C), gas velocity and temperature, dust velocity and its dispersion. It is also possible to change the size distribution and type of dust: silicates, HAC grains or graphites+PAHs. The presented numerical tool is publicly available.

We have performed test calculations and presented the dependence of characteristic time-scales on radiation field, gas velocity and temperature. For the last two parameters we constructed the maps of destruction and aromatisation time-scales. These calculations show that the dust evolution is a complicated process and that photo-destruction, sputtering and shattering must be considered simultaneously. We also demonstrated evolutionary changes of the infrared dust emission spectra due to destruction and aromatisation processes in the media with enhanced radiation field and with high gas velocity. Beside the conditions in a PDR and a supernova shock, we also applied our model to the conditions in an IR bubble around an expanding Hii region and show how the model can be used to estimate time-scales of dust evolution in these objects.

We have considered the influence of dust charge on the evolutionary processes and concluded that the charge influences the photo-destruction process. Changes in ionization potential and cooling rate for charged dust grains leads to their faster destruction in comparison to neutral grains.

The present model does not account for hydrogenation of HACs and their subsequent aliphatisation, volatile accretion and coagulation processes that dominate in the cold neutral medium, so it cannot be used for modelling molecular clouds or protoplanetary discs, where these processes are crucial.

The authors are very grateful to Anthony Jones for interest to this work, fruitful discussions and comments. We also thank Takaya Nozawa for providing the shock profiles and Sergey Khoperskov for stimulating discussions. The CVODE package (SUNDIALS software) was used in the work and we express our gratitude to the SUNDIAL team for providing this tool. The work was supported by the RFBR grants 17-02-00521 and 18-52-52006.

Appendix A Sputtering rate coefficients

A rate coefficient of aromatisation due to inelastic collision with thermal and non-thermal ions (H, He, C) and electrons for grains with less than 1000 atoms is calculated as follows:

[TABLE]

A rate coefficient of sputtering of carbon atoms due to inelastic collision with thermal and non-thermal ions (H, He, C) and electrons for grains with less than 1000 atoms is calculated as follows:

[TABLE]

A rate coefficient of aromatisation due to elastic collision with thermal and non-thermal ions (H, He, C) for grains with less than 1000 atoms is calculated as follows:

[TABLE]

A rate coefficient of sputtering of carbon atoms due to elastic collision with thermal and non-thermal ions (H, He, C) for grains with less than 1000 atoms is calculated as follows:

[TABLE]

A rate coefficient of aromatisation due to elastic collision with thermal and non-thermal ions (H, He, C) for grains with more than 1000 atoms is calculated as follows:

[TABLE]

A rate coefficient of sputtering of carbon atoms due to elastic collision with thermal and non-thermal ions (H, He, C) for grains with more than 1000 atoms is calculated as follows:

[TABLE]

Expressions for ejection rates of C and H atoms were given in Paper I.

Appendix B Photo-destruction rate coefficients

A rate of change of band gap energy $\varepsilon^{(1)}_{ij}$ due to aromatisation by photons can be expressed as:

[TABLE]

where $\delta$ is a conversion coefficient between hydrogen atom fraction $X_{\rm H}$ and band gap energy, $E_{\rm g}=\delta X_{\rm H}$, $Y_{\rm diss}^{\rm CH}$ is the dissociation probability, $\sigma^{\rm CH}_{\rm loss}$ is the C–H dissociation cross section. We adopt $\delta=4.3$ eV (Tamor & Wu, 1990). Following Jones et al. (2014), we adopt $Y_{\rm diss}^{\rm CH}=1$ for grains with $a<20$ Å and $Y_{\rm diss}^{\rm CH}=20$Å$/a$ for larger grains, and $\sigma^{\rm CH}_{\rm loss}=10^{-19}$ cm2.

In the model of Allain et al. (1996a), which is used to calculate the second step of destruction, when the grain is dehydrogenated, a hydrogen atom, a hydrogen molecule, or an acetylene molecule are detached from a grain with the rates of $k_{\rm H}$, $k_{\rm H_{2}}$, $k_{\rm C_{2}H_{2}}$, s*-1*, correspondingly, which are calculated using the RRK (Rice-Ramsperger-Kassel) theory (Nic et al., 2012) based on the experimental data (Jochims et al., 1994). The sum of these rates is the total rate dissociation coefficient, $k_{\rm diss}=k_{\rm H}+k_{\rm H_{2}}+k_{\rm C_{2}H_{2}}$. The dissociation probability of C–C bond by a single photon (sp) can be expressed as:

[TABLE]

where $k_{\rm IR}$ is a rate coefficient of IR photon emission (see more details in Paper I), and factor 2 means that two C–C bonds are destroyed after acetylene detachment. We adopt $Y_{\rm ion}$ from the work of Weingartner & Draine (2001a) developed for PAHs and graphite-like grains supposing that ionization yields for dehydrogenated HACs are similar to those for PAHs. While there is no calculated dependence of $Y_{\rm ion}$ for HAC on their size, energy and charge, it was experimentally found that the work function of HAC films is about 4–5 eV depending on the sp${}^{2}/{\rm sp}^{3}$ ratio (Ilie et al., 2000). The analogous value for PAHs is 4.4 eV, which is somewhere between the values for the HAC material. In addition, this value almost does not depend on the film thickness if the HAC material includes sp2-nanosized clusters (Shakerzadeh et al., 2012). So, these studies give grounds to believe that our assumption on $Y_{\rm ion}$ of HACs is reasonable.