Evolution of dust extinction curves in galaxy simulation

TL;DR

This study uses galaxy simulations to model dust evolution and extinction curves, revealing how dust processes and galaxy age influence observed extinction features like the 2175 Å bump and FUV rise.

Contribution

It introduces a comprehensive dust evolution model in galaxy simulations that reproduces observed extinction curves across different galaxy ages and environments.

Findings

Extinction curves are flat at early galaxy ages due to stellar dust production.

Dust growth by accretion enhances the 2175 Å bump and FUV rise over time.

Simulated extinction curves at >3 Gyr match Milky Way observations.

Abstract

To understand the evolution of extinction curve, we calculate the dust evolution in a galaxy using smoothed particle hydrodynamics simulations incorporating stellar dust production, dust destruction in supernova shocks, grain growth by accretion and coagulation, and grain disruption by shattering. The dust species are separated into carbonaceous dust and silicate. The evolution of grain size distribution is considered by dividing grain population into large and small gains, which allows us to estimate extinction curves. We examine the dependence of extinction curves on the position, gas density, and metallicity in the galaxy, and find that extinction curves are flat at $t \lesssim 0.3$ Gyr because stellar dust production dominates the total dust abundance. The 2175 \AA\ bump and far-ultraviolet (FUV) rise become prominent after dust growth by accretion. At $t \gtrsim 3$ Gyr, shattering works efficiently in the outer disc and low density regions, so extinction curves show a very strong 2175 \AA\ bump and steep FUV rise. The extinction curves at $t\gtrsim 3$ Gyr are consistent with the Milky Way extinction curve, which implies that we successfully included the necessary dust processes in the model. The outer disc component caused by stellar feedback has an extinction curves with a weaker 2175 \AA\ bump and flatter FUV slope. The strong contribution of carbonaceous dust tends to underproduce the FUV rise in the Small Magellanic Cloud extinction curve, which supports selective loss of small carbonaceous dust in the galaxy. The snapshot at young ages also explain the extinction curves in high-redshift quasars.

Tables

Table 1. Table 1: Parameters adopted in equations ( 2 )–( 8 ). In this work, we set the solar elemental abundance pattern ( Z X ⁣ ⊙ subscript 𝑍 X direct-product Z_{\mathrm{X\odot}} ) for C and Si, and assume that Si occupies a mass fraction of 0.166 in silicate. ρ X subscript 𝜌 X \rho_{\mathrm{X}} is the material density, and f in , X subscript 𝑓 in X f_{\mathrm{in,X}} is condensation efficiency in stellar ejecta (see Section 2.1.1 )

Species	X	$Z_{X\odot }$	${\rho }_X\left(g{\mathrm{cm}}^{-3}\right)$	$f_{\mathrm{in},X}$
Carbonaceous dust	C	$2.47\times {10}^{-3}$	2.24	0.34
Silicate	Si	$8.17\times {10}^{-4}$	3.5	0.13

Equations

$${\mathcal{D}}_{i,\mathrm{X}(j)}\equiv m_{i,\mathrm{X}(j)}/m_{\mathrm{gas}(j)},$$

$${\mathcal{D}}_{\mathrm{L},\mathrm{X}}(t+\Delta t)$$

$$-\left(\frac{{\mathcal{D}}_{\mathrm{L},\mathrm{X}}(t)}{{\tau }_{\mathrm{sh}}}-\frac{{\mathcal{D}}_{\mathrm{S},\mathrm{X}}(t)}{{\tau }_{\mathrm{co}}}\right)\Delta t$$

$$+f_{\mathrm{in},\mathrm{X}}{\mathcal{Y}}_{\mathrm{X}}^{\prime }\frac{\Delta M_{\mathrm{return}}}{m_{\mathrm{gas}}}(1-\delta ),$$

$${\mathcal{D}}_{\mathrm{S},\mathrm{X}}(t+\Delta t)$$

$$+\left(\frac{{\mathcal{D}}_{\mathrm{L},\mathrm{X}}(t)}{{\tau }_{\mathrm{sh}}}-\frac{{\mathcal{D}}_{\mathrm{S},\mathrm{X}}(t)}{{\tau }_{\mathrm{co}}}+\frac{{\mathcal{D}}_{\mathrm{S},\mathrm{X}}(t)}{{\tau }_{\mathrm{acc}}}\right)\Delta t,$$

$$\Delta {\mathcal{D}}_{L/S,\mathrm{X}}^{\mathrm{SN}}(t)=\left[1-(1-\eta {)}^N\right]{\mathcal{D}}_{L/S,\mathrm{X}}^{\mathrm{SN}}(t),$$

$${\tau }_{\mathrm{sh}}=5.408\times 10^8\mathrm{yr}\left(\frac{{\rho }_{\mathrm{X}}}{3\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}\right){\left(\frac{n_{\mathrm{gas}}}{1{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{{\mathcal{D}}_{\mathrm{L},\mathrm{X}}}{0.01}\right)}^{-1},$$

$${\tau }_{\mathrm{co}}=2.704\times 10^5\mathrm{yr}\left(\frac{{\rho }_{\mathrm{X}}}{3\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}\right){\left(\frac{{\mathcal{D}}_{\mathrm{S},\mathrm{X}}}{0.01}\right)}^{-1}/f_{\mathrm{dense}},$$

$${\tau }_{\mathrm{acc},\mathrm{C}}=0.95\times 10^6\mathrm{yr}{\left(\frac{Z_{\mathrm{C}}}{Z_{\mathrm{C}\odot }}\right)}^{-1}{\left(1-\frac{{\mathcal{D}}_{\mathrm{tot},\mathrm{C}}}{Z_{\mathrm{C}}}\right)}^{-1}/f_{\mathrm{dense}},$$

$${\tau }_{\mathrm{acc},\mathrm{Si}}=1.05\times 10^6\mathrm{yr}{\left(\frac{Z_{\mathrm{Si}}}{Z_{\mathrm{Si}\odot }}\right)}^{-1}{\left(1-\frac{{\mathcal{D}}_{\mathrm{tot},\mathrm{Si}}}{Z_{\mathrm{Si}}}\right)}^{-1}/f_{\mathrm{dense}},$$

$$n_{i,\mathrm{X}}(a)=\frac{C_{i,\mathrm{X}}}{a^4}\exp \left\{-\frac{{\left[\ln (a/a_{0,i})\right]}^2}{2{\sigma }^2}\right\},$$

$$\mu m_{\mathrm{H}}{\mathcal{D}}_{i,\mathrm{X}}={\int }_0^{\infty }\frac{4}{3}\pi a^3{\rho }_{\mathrm{X}}{n}_{i,\mathrm{X}}(a)\mathrm{d}a,$$

$$\frac{A_{\lambda ,\mathrm{X}}}{N_{\mathrm{H}}}=2.5{\log }_{10}\mathrm{e}\sum _i{\int }_0^{\infty }{n}_{i,\mathrm{X}}(a)\pi a^2{Q}_{\mathrm{ext}}(a,\lambda ,\mathrm{X})\mathrm{d}a,$$

Full text

††pubyear: 2016††pagerange: Evolution of dust extinction curves in galaxy simulation–Evolution of dust extinction curves in galaxy simulation

1 Introduction

Studying dust is necessary for the understanding of galaxy evolution. Dust surfaces are the main sites for the efficient production of molecular hydrogen, which is the main constituent of molecular clouds and an important coolant in low-metallicity environments (e.g. Hirashita & Ferrara, 2002; Cazaux & Spaans, 2004). Moreover, dust cooling determines the mass of final fragments of star-forming clouds (Omukai et al., 2005; Schneider et al., 2006). Therefore, star formation in galaxies is strongly affected by dust grains. In terms of the observational properties of galaxies, dust shapes the spectral energy distribution (SED) by absorbing stellar light and re-emitting it in the far infrared (e.g. Yajima et al., 2014; Schaerer et al., 2015, for recent modelling); thus, correcting for dust extinction is required to derive the intrinsic stellar SED.

Extinction curves, which represent the wavelength dependence of dust extinction (absorption + scattering), provide useful information on the grain size distribution and chemical composition (e.g. Weingartner & Draine, 2001), both of which are important in determining the aforementioned processes, i.e. grain surface reactions and dust cooling (Yamasawa et al., 2011). The Milky Way extinction curve has been observed in detail (Pei, 1992; Fitzpatrick & Massa, 2007), and by fitting that, Mathis et al. (1977) derived a grain size distribution of a power law $\propto a^{-3.5}$ in the range of $a\simeq 0.005$–0.25 $\micron$ ($a$ is the grain radius) with a mixture of silicate and graphite. Furthermore, Weingartner & Draine (2001) suggested detailed functional forms of grain size distribution for carbonaceous dust and silicate to explain the extinction curves of the Milky Way and Large/Small Magellanic Clouds (LMC/SMC). Although these models are successful in explaining the extinction curves, it is still necessary to clarify what processes govern the evolution of grain size distribution in galaxies.

There have been a lot of efforts not only to model the evolution of dust abundance in the solar neighbourhood (Dwek & Scalo, 1980; Zhukovska et al., 2008), and in both local galaxies (Dwek, 1998; Lisenfeld & Ferrara, 1998; Hirashita, 1999; Inoue, 2003) and high redshift galaxies (Kuo & Hirashita, 2012; de Bennassuti et al., 2014; Mancini et al., 2016; Popping et al., 2016; Wang et al., 2017), but also to model the evolution of grain size distribution in galaxies (Liffman & Clayton, 1989; O’Donnell & Mathis, 1997; Yamasawa et al., 2011). Asano et al. (2013b), treating a galaxy as a single zone, constructed a framework for calculating the evolution of grain size distribution in a consistent manner with galaxy evolution over the entire galaxy history. Furthermore, Asano et al. (2014) used their model to calculate the evolution of extinction curve and found the following evolutionary features. The extinction curve is flat at the earliest evolutionary stage because the dust content is dominated by stellar sources. The subsequent contribution from shattering and accretion increases the abundance of small grains drastically; as a consequence, the extinction curve becomes steep, and eventually even steeper than the Milky Way extinction curve. Using the same framework, Nozawa et al. (2015) successfully reproduced the Milky Way extinction curve by newly considering dense clouds that increases the efficiency of coagulation. The above studies basically treat a galaxy as a single-zone object, and this treatment inevitably and implicitly adopt some strong assumptions such as instantaneous mixing and homogeneity (i.e., the spatial structure in a galaxy is neglected). Fitzpatrick & Massa (2007) showed the variation of the Milky Way extinction curves toward different lines of sight, which indicates the inhomogeneity in dust properties (see also Nozawa & Fukugita, 2013). This inhomogeneity might cause uncertainty in interpreting the extinction properties of unresolved extragalactic objects.

To calculate the evolution of dust in realistic conditions and in a spatially resolved manner, we use a hydrodynamical simulation of a galaxy. Hydrodynamic simulations have been recognized as a useful tool in studying galaxy formation and evolution. There have been some hydrodynamic simulations incorporating the effects of dust in recent years. Yajima et al. (2015) estimated the dust abundance in high-resolution zoom-in cosmological simulations of galaxies with a fixed dust-to-metal ratio. The ultraviolet (UV) and infrared luminosities of individual galaxies are estimated by post-processing using a radiative transfer code. Bekki (2015) performed a smoothed particle hydrodynamics (SPH) simulation with a dust evolution model, which includes dust formation in stellar ejecta, dust growth by accretion, and dust destruction in SN shocks. In particular, they treated dust as new particles besides dark matter, gas and stellar particles. McKinnon et al. (2016b) treated dust as an attached property of gas. They implemented dust formation and destruction in a moving-mesh code and performed cosmological zoom-in simulations. Their result showed that dust growth by accretion is important and that an accurate treatment of stellar and active galactic nuclei feedback is needed to reproduce the observation. Moreover, they also ran cosmological simulations (McKinnon et al., 2016a) to investigate the dust mass function and radial profile of dust in the circum-galactic environment. The simulated dust mass function is consistent with observations in the local Universe but has a tendency of underestimating the number of dust-rich galaxies at high redshift. Zhukovska et al. (2016) focused on the dominant processes, dust growth and destruction in the interstellar medium (ISM), to examine the temperature-dependent sticking coefficient in the accretion of gas-phase metals onto dust. However, the grain size distribution is not considered in the above simulations.

Our previous work (Aoyama et al., 2017, hereafter Paper I) simulated dust abundance in an isolated galaxy using SPH simulations with the dust enrichment model developed by Hirashita (2015). The model includes the dust production in stellar ejecta, destruction in supernova (SN) shocks, dust growth by accretion, grain growth by coagulation, and grain disruption by shattering. Shattering and coagulation, which are caused by grain–grain collisions, are important in determining the grain size distribution. To treat the evolution of grain size distribution within a reasonable computational cost, Paper I adopted the two-size approximation proposed by Hirashita (2015) in which the whole range of grain size is represented by two sizes divided at around $0.03\mu$m. In Paper I, we found the following evolutionary features of dust. Stellar dust production dominates the dust content at young ages ($t\lesssim$ 0.2 Gyr). After that, dust processing in the ISM is much more important than stellar dust production. In particular, dust growth by the accretion of gas-phase metals is triggered by small grains produced as a result of shattering, and drastically increases the dust-to-gas ratio and the dust-to-metal ratio. Coagulation becomes efficient at later times of $t\gtrsim$ 1 Gyr. We also computed the radial profile of dust-to-gas mass ratio and dust-to-metal mass ratio (i.e. depletion), both of which are broadly consistent with the observational data of a sample of spatially resolved nearby galaxies in Mattsson & Andersen (2012).

In this paper, we further compare our simulation results with observed extinction curves in the nearby Universe to obtain more constraints on the grain size distribution. To calculate extinction curves, we separate the dust species into carbonaceous dust and silicate. For the first time, extinction curves can be predicted by hydrodynamical galaxy evolution simulations with a physically plausible model including dust production and destruction. Taking advantage of the simulations, we examine how the evolution of extinction curve depends on the physical conditions such as the galactic radius, gas density and metallicity. We compare the calculated extinction curves with observational extinction curves in the Milky Way, LMC and SMC, in which a precise determination of extinction curves are possible because individual bright stars can be spatially resolved (e.g. Pei, 1992; Gordon et al., 2003). We also discuss the evolution of extinction curves and its implication for the extinction curves at high redshifts.

This paper is organized as follows. In Section 2, we summarise the simulation, the dust evolution model and the extinction curve calculation. The evolution of extinction curves and the dependence on the radial distance, gas density, and metallicity are shown in Section 3. In Section 4, we make a comparison with observational data and discuss the results. Finally, we provide the conclusions in Section 5.

2 Model

2.1 Galaxy simulation with dust enrichment

We briefly summarize the simulation with the dust enrichment model. Our model is based on Paper I, and the main difference in this work is that we separate the dust species into carbonaceous dust and silicate. Separate treatment of grain species is fundamental in calculating extinction curves (e.g. Draine & Lee, 1984). We focus on the difference from Paper I, but also provide a brief summary of Paper I. We refer the interested reader to Paper I for details.

We adopt the modified version of GADGET-3 $N$-body/SPH code for this study (Springel, 2005), and the Grackle111https://grackle.readthedocs.org/ chemistry and cooling library (Bryan et al., 2014; Kim et al., 2014; Smith et al., 2017) to solve non-equilibrium primordial chemistry network. Initial condition is the low-resolution model of AGORA simulations (Kim et al., 2014) with the minimum gravitational softening length of $\epsilon_{\rm grav}=80$ pc. The star formation is assumed to occur on a local free-fall time with an efficiency of 0.01. Stellar feedback from SNe and stellar winds is also taken into account. All the feedback and metal enrichment from stars are assumed to occur in 4 Myr after the star formation; thus, in our framework, the injection of metals and dust occurs almost instantaneously after the star formation, and it is difficult to include delayed dust and metal input from asymptotic giant branch (AGB) stars. Although we discuss the effect of dust production by AGB stars later (Section 2.1.1), we leave the treatment of delayed dust production by AGB stars to the future work.

For dust enrichment, we adopt the two-size dust enrichment model developed by Hirashita (2015). In this model, we divide the whole grain size range into two at around $a\sim 0.03~{}\micron$ ($a$ is the grain radius) and consider the evolution of large and small grains. This approximate treatment is introduced because solving the full grain size distribution in individual SPH particles is computationally too expensive. Hirashita (2015) has already clarified that the two-size approach reproduces the evolution of extinction curves calculated by a full grain size treatment in Asano et al. (2014) if a functional form of grain size distribution is properly chosen (see Section 2.2).

For dust evolution processes, we consider stellar dust production, SN destruction, grain disruption by shattering in the diffuse ISM, and grain growth by coagulation and accretion in the dense ISM. On each gas (SPH) particle, we consider the mass of small and large grains, $m_{i,\mathrm{X}(j)}$, where subscripts $i$, X, and $j$ denote the grain size ($i=$ L and S for large and small grains, respectively), grain species (X=Si and C for silicate and carbonaceous dust, respectively), and the label of the gas particle, respectively. We formulate the dust evolution using the dust-to-gas mass ratio:

[TABLE]

where $m_{\mathrm{gas}(j)}$ is the gas mass of the $j$-th particle. We calculate the evolution of $\mathcal{D}_{i,\mathrm{X}(j)}$ for each gas particle. We follow Paper I for the basic equations of dust evolution except that we consider two species. The time evolution of the large and small grain abundances in the $j$-th particle from time $t$ to the next time step $t+\Delta t$ is formulated as (for brevity, we omit the subscript $j$, and we always solve the equations for each gas particle):

[TABLE]

where $\mathcal{D}_{{\rm L,X}}$ ($\mathcal{D}_{{\rm S,X}}$) is large (small) grain dust-to-gas mass ratio defined in Eq. (1), $\Delta\mathcal{D}^{\rm SN}_{i,\mathrm{X}}(t)$ is the abundance of the pre-existing dust destroyed by SN shock, $f_{\rm in,X}$ is dust condensation efficiency in stellar ejecta, $\Delta M_{{\rm return}}$ is returned gas mass, $\delta$ is the fraction of newly formed dust that is destroyed by SNe, and $\mathcal{Y}^{\prime}_{\rm X}$ is the metal yield, and $\tau_{\rm sh}$, $\tau_{\rm co}$ and $\tau_{\rm acc}$ are the time-scales of shattering, coagulation and accretion respectively. For the metal yield, we adopt $\mathcal{Y}^{\prime}_{\rm Z}=0.02896$ (Paper I) and multiply $Z_{\mathrm{X\odot}}/Z_{\odot}$ (the mass ratio of the element X to the total metal content in the solar abundance) for the Si and C yields by assuming the solar metallicity pattern (Lodders, 2003), where $Z_{\mathrm{X\odot}}$ is the solar abundance (X = Si and C). Following Paper I, the solar metallicity $Z_{\odot}=0.02$ is adopted. We notice that the precise solar metallicity is less than 0.02 (Lodders, 2003; Asplund et al., 2009). However, in this work, the solar metallicity is just for normalizing metallicity in the figures, and the results are not affected by which value we adopted. We assume that Si occupies a mass fraction of 0.166 in silicate while C is the only constituent of carbonaceous dust. The parameters adopted in this work are summarised in Table 1. For simplicity, we assume that the two dust species evolve independently in this work.

We briefly explain each dust evolution process in what follows.

2.1.1 Stellar dust formation

Stars produce dust at their final stage of evolution. The formed dust is distributed to the surrounding SPH particles in the same way as metals by assuming the dust condensation efficiency $f_{\mathrm{in,X}}$. The dust condensation efficiencies in stellar ejecta for carbonaceous dust and silicate (X = C and Si) are calculated based on Hou et al. (2016) except that we adopt a Chabrier (2003) initial mass function with a stellar mass range of 0.1–100 $M_{\odot}$. The values are listed in Table 1. In our simulation, metals are assumed to be produced only by SNe. Thus, to derive the dust condensation efficiency in stellar ejecta, we integrate the metal and dust yields of SNe for the progenitors in the mass range $m=8$–40 $M_{\odot}$. The metal and dust yields taken from Kobayashi et al. (2006) and Nozawa et al. (2007, with an ambient hydrogen number density of 1 cm*-3* and the unmixed helium core), respectively. We assume that neither metals nor dust is ejected from stars whose mass is above 40 $M_{\odot}$ (Heger et al., 2003).

In this study, dust production by AGB stars is not included for the following two reasons. The first reason is that the stellar dust production is more important in the early stage ($\lesssim 0.3$ Gyr) than in the later stage, at least in our simulation: at $\lesssim 0.3$ Gyr, the progenitors of AGB stars have no time to evolve to a dust production phase (Valiante et al., 2009, 2011). The second reason is related to a technical aspect: our simulation assumes almost instantaneous metal production after star formation. Thus, the metal and dust production by AGB stars, whose progenitor lifetimes have a large variety, cannot be included in the current simulation framework. We could effectively include the dust production by AGB stars by adopting the dust condensation efficiency ($f_{\mathrm{in,X}}$) that takes into account the contribution from AGB stars. Paper I already clarified the effects of changing $f_{\mathrm{in,X}}$. Extinction curves are governed by the small-to-large grain abundance ratio ($\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$) in our model, but the value of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ is affected not by stellar dust production but by interstellar dust processing, especially accretion, at the epoch when AGB stars contribute to the dust production ($t\gtrsim 1$ Gyr). At this stage, the dust evolution is insensitive to stellar dust production or $f_{\mathrm{in,X}}$. Moreover, as shown by Asano et al. (2013a), the epoch at which interstellar processing starts to dominate the dust evolution over stellar dust production is determined by the metallicity and is insensitive to the choice of $f_{\mathrm{in,X}}$. Therefore, we do not attempt a detailed modelling of $f_{\mathrm{in,X}}$ by including the contribution from AGB stars, although we plan to modify the framework to include the delayed dust input from AGB stars in the future. The parameter $f_{\mathrm{in,X}}$ is important for determining the $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ and C/Si ratio at early evolutionary stage, which will be discussed in Sections 4.2 and 4.3.

2.1.2 SN destruction

Each SN destroys the dust in its sweeping radius. Because the simulations do not resolve individual SNe, a sub-grid model is necessary to treat dust destruction by SNe. We separately treat the destruction of the pre-existing dust and the newly formed dust by the same SNe. The terms for the destruction of the pre-existing dust in equations (2) and (3) can be written as

[TABLE]

where $\eta\equiv\min[\epsilon_{\mathrm{SN}}(m_{\mathrm{sw}}/m_{\mathrm{gas}}),\,\epsilon_{\mathrm{SN}}]$ ($m_{\mathrm{sw}}$ is the gas mass swept by a single SN; see Paper I for its estimation, and $\epsilon_{\mathrm{SN}}$ is the efficiency of dust destruction in a single SN blast), and $N$ is the number of SN explosions in the gas particle. We adopt $\epsilon_{\mathrm{SN}}=0.4$ based on Nozawa et al. (2006). The SN destruction of newly formed dust is absorbed in the stellar ejecta term by giving a destruction fraction $\delta$ in Eq. (2). The derivation of this quantity was shown in Paper I.

2.1.3 Shattering

Shattering is a process in which large grains collide with each other and shattered or fragmented into small grains. Shattering occurs in the diffuse ISM (Hirashita & Yan, 2009) because it needs high relative velocity between grains; therefore, we allow for shattering only in gas particles with $n_{\mathrm{gas}}<1~{}\mathrm{cm}^{-3}$.

The shattering time-scale is estimated as

[TABLE]

where $\rho_{\mathrm{X}}$ is the material density of the dust species listed in Table 1. For the particles with $n_{\mathrm{gas}}\geq 1~{}\mathrm{cm}^{-3}$, we turn off shattering, i.e. $\tau_{\mathrm{sh}}=\infty$.

2.1.4 Coagulation and accretion

Coagulation is a process in which small grains collide with each other and turn into large grains. It occurs under the condition that grains have low relative velocity. Accretion is a process in which small grains gain their mass by accreting gas-phase metals. Coagulation and accretion are expected to occur in dense clouds, where relative velocities of grains and gas are low but collision still happens because of high number density; however, our simulations are not capable of resolving them since the highest density is determined by the minimum gravitational softening length of $\epsilon_{\rm grav}=80$ pc, which is larger than the typical size of a dense cloud. We apply a sub-grid model for gas particles with $n_{\mathrm{gas}}\geq 10~{}\mathrm{cm}^{-3}$ and $T_{\mathrm{gas}}<10^{3}$ K by assuming that a fraction, $f_{\rm dense}$, of gas content is in the dense cloud phase with fixed density and temperature as $n_{\mathrm{gas}}=10^{3}~{}\mathrm{cm}^{-3}$ and $T_{\mathrm{gas}}=50$ K. The time-scale of coagulation is formulated as

[TABLE]

and the time-scales of accretion for carbonaceous dust and silicate are formulated as

[TABLE]

where $\mathcal{D}_{\mathrm{tot,X}}=\mathcal{D}_{\mathrm{S,X}}+\mathcal{D}_{\mathrm{L,X}}$ (X = C and Si). We assume $f_{\rm dense}=0.5$ as a fiducial value. We set $\tau_{\mathrm{acc,X}}=\infty$ (X = C and Si) and $\tau_{\mathrm{co}}=\infty$ if the particles do not satisfy the criteria ($n_{\mathrm{gas}}\geq 10~{}\mathrm{cm}^{-3}$ and $T_{\mathrm{gas}}<10^{3}$ K).

2.2 Extinction curve

We calculate the extinction curve for each gas particle in the following way. To calculate the extinction curve, grain size distribution is required; however, the whole range of grain size is represented by the large and small grain populations in our simulation. Thus, we need to assume a specific functional form for the grain size distribution of each population. Following Hirashita (2015), we adopt a modified-lognormal function for the grain size distribution. In the $j$th particle, the grain size distribution is (we omit the subscript $j$):

[TABLE]

where subscript $i$ indicates the small ($i=\mathrm{S}$) or large ($i=\mathrm{L}$) grain component, subscript X represents dust species [carbonaceous dust (X = C) or silicate (X = Si)], $C_{i,\mathrm{X}}$ is the normalization constant for the $j$th particle and $a_{0,i}$, and $\sigma$ are the central grain radius and the standard deviation of the lognormal part, respectively. We adopt $a_{\mathrm{0,s}}=0.005~{}\micron$, $a_{\mathrm{0,l}}=0.1~{}\micron$ and $\sigma=0.75$ since these values reproduce the Milky Way extinction curve when the small-to-large grain abundance ratio is the same as the Mathis et al. (1977, MRN) size distribution (Hirashita, 2015). The normalization $C_{i,\mathrm{X}}$ is determined by

[TABLE]

where $\mu=1.4$ is the gas mass per hydrogen nucleus, and $m_{\mathrm{H}}$ is the mass of hydrogen atom.

The extinction $A_{\lambda,\mathrm{X}}$ (in units of magnitude) normalized to the column density of hydrogen nuclei ($N_{\mathrm{H}}$) is written as

[TABLE]

where $Q_{\mathrm{ext}}(a,\lambda,\mathrm{X})$ is the extinction coefficient (extinction cross section normalized to the geometric cross section) as a function of grain size, wavelength and dust species. $Q_{\mathrm{ext}}(a,\lambda,\mathrm{X})$ is calculated by using the Mie theory (Bohren & Huffman, 1983) based on the same optical constants for silicate and graphite as in Weingartner & Draine (2001) unless otherwise stated, on the assumption that all the grains are spherical with a uniform composition.

3 Results

The new features calculated by our models are the time evolution and the spatially resolved information of extinction curves. In this section, we describe how the calculated extinction curves depend on various quantities that characterize the environmental or evolutionary properties. We focus on the silicate–graphite modelling as a standard (Draine & Lee, 1984; Pei, 1992), but the qualitative behaviour of overall steepness (other than the carbon bump) does not change even if we adopt other dust species (we mention other possible dust species in Section 4.5).

Following Paper I, we adopt the ages 0.3, 1, 3, and 10 Gyr to represent four galaxy evolution stages that cover a wide range of specific star formation rate (sSFR), 3–0.1 Gyr*-1*. To discuss the extinction curves, we mainly focus on two features, 2175 Å bump and far-UV (FUV) slope. Since the extinction curves reflect the variation of grain size distribution, the small-to-large grain abundance ratio ($\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$) determines the extinction curve in our model.

3.1 Radial dependence

We investigate the radial dependence here to represent the spatial variation. For extinction curves, we usually obtain an average over a line of sight, which is often through a galactic disc for a disc galaxy. Thus, for the radial dependence, we consider a scale larger than the disc scale height (but still smaller than the disc scale length) by adopting a radius bin width of 2 kpc and take the mass-weighted average for all gas particles in each radius bin. We analyse extinction curves up to $r$ = 8 kpc ($r$ is the radius in cylindrical coordinate), which is close to $R_{25}$ of our simulated galaxy. The centre of the galaxy is determined by the centre of mass for the stellar particles. Fig. 1 shows the extinction curves in the four radial bins in the four evolution epochs.

We also show the small-to-large grain abundance ratio as a function of radial distance for silicate and carbonaceous dust (note that we plot the distribution of all gas particles on this diagram), since this ratio determines the shape of the extinction curve. In general, a higher small-to-large grain abundance ratio ($\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$) shows an extinction curve with a more prominent 2175 Å bump and a steeper UV slope. The evolution of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ has already been investigated in Paper I, but we repeat some discussions on $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$, since this is the central quantity in interpreting the calculated extinction curves. We refer the interested reader to Paper I for detailed discussions on $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$.

We observe in Fig. 1 that the features of extinction curves (the 2175 Å bump and the UV slope) become stronger with age, which is caused by the increase of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$. The small-to-large grain abundance ratio increases in the centre first since dust growth by accretion becomes the most efficient in the most metal-enriched region. The regions with high small-to-large grain abundance ratio extend to the outer radii with age.

At 0.3 Gyr, only the extinction curve of the central region shows a slight FUV rise, and the other regions have flat extinction curves. A portion of gas particles have high small-to-large grain abundance ratio in the central region, where some particles are dense and already metal-rich so that accretion already enhances the $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$. The small-to-large grain abundance ratios at outer radii are smaller than $10^{-2}$ and $10^{-3}$ for carbonaceous dust and silicate, respectively. The low $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ is due to the dominance of stellar dust production, and the existence of small amounts of small grains is explained by shattering.

At 1 Gyr, the extinction curves at all radii show a prominent 2175 Å  bump and FUV rise, and these extinction curve features are the strongest in the central part of the galaxy. The small-to-large grain abundance ratio is greater than $10^{-1}$ at $r\lesssim 6$ kpc. There is a transition of the small-to-large grain abundance ratio at $6\lesssim r\lesssim 8$ kpc, and silicate shows a stronger transition than carbonaceous dust (see Section 4.2 for further discussions). This clear rise of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at $r\sim$ 6–8 kpc is caused by a sensitive dependence of accretion on metallicity. Therefore, the strong radial dependence of extinction curves is driven by accretion in this stage.

The radial variation of extinction curves at 3 Gyr shows the opposite trend to that at 1 Gyr, with outer regions having steeper extinction curves. The small-to-large grain abundance ratio shows a decrease toward the centre. Because accretion dominates the dust evolution in the entire galaxy regions, the small-to-large grain abundance ratio is as high as $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ $\gtrsim 0.1$. The decrease of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ in the central region is caused by coagulation since the ISM is relatively dense in the central region.

At 10 Gyr, the features in the extinction curves become more prominent than at 3 Gyr except in the centre. At this stage, the small grain abundance is almost comparable to the large grain abundance ($\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ $\sim 1$) in the outer disc, and $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ is smaller in the centre than at 3 Gyr. The highest value of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ gradually increases at $t\gtrsim$ 3 Gyr. We notice that shattering continuously converts large grains to small grains, while coagulation becomes less efficient because of the general trend that the gas becomes less dense in the late evolutionary stage mainly as a result of gas consumption by star formation.

3.2 Density dependence

In our model, the efficiencies of shattering, coagulation and accretion have gas density dependence, so it is worth investigating how the extinction curve varies in different density environments. To this aim, we separate the gas particles according to the logarithmic density in $n_{\mathrm{gas}}=0.1$ to $10^{3}$ cm*-3* and derived the mass-weighted average in each density bin. The resulting extinction curves are shown in Fig. 2.

In the early evolutionary stage ($t\sim 0.3$ Gyr), all the density ranges show a weak 2175 Å bump and a flat FUV curve. There is a slight trend of higher densities showing relatively steep FUV slopes. Although $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ increases with density, most of the particles have low $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$, which is determined by stellar dust production. The effect of shattering and accretion can be seen in a small fraction of particles with high density, causing the upward spread of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at $n_{\mathrm{gas}}\gtrsim 1$ cm*-3*. Since the steepening of FUV slope is seen in the dense medium, it is caused by accretion rather than shattering, although shattering is also necessary to produce the ‘seed’ small grains for accretion.

At 1 Gyr, all the density ranges show a prominent 2175 Å bump and a steep FUV rise. These features become clearer as the density becomes higher. Therefore, at this stage, the extinction curves are steeper in the dense ISM than in the diffuse ISM. This is explained by the $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ ratio as follows. The small-to-large grain abundance ratio rapidly increases in $n_{\mathrm{gas}}\gtrsim 0.1$ $\mathrm{cm}^{-3}$ from 0.3 to 1 Gyr. Gas particles with $n_{\mathrm{gas}}\gtrsim 10$ $\mathrm{cm}^{-3}$ all reach $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}\gtrsim 0.1$ by active accretion. When such dense regions become diffuse by stellar feedback, etc., the small grains formed in the dense regions are dispersed into the diffuse ISM. Thus, the high-$\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ regions extend down to low densities. In the diffuse regions with $n_{\mathrm{gas}}<10$ cm*-3*, there is another concentration of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ around $\sim 10^{-2}$ in carbonaceous dust ($\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}\sim 10^{-3}$ in silicate), which is the branch of inefficient accretion (i.e. the main source of small grains is shattering). Thus, there is a large dispersion in $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at $n_{\mathrm{gas}}<10$ cm*-3*. We expect that, although the extinction curves in the diffuse ISM are on average flatter than in the dense ISM, the dispersion in the curves is larger than the difference between the mean extinction curves in the dense and diffuse ISM because of the large dispersion in $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$.

At 3 Gyr, the trend of extinction curves for the density is reversed: the diffuse regions have extinction curves with a high 2175 Å bump and a steep FUV rise. There are few particles that have $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}<10^{-2}$ over all the density range, because the increase of small grains by accretion has been prevalent in the entire galaxy after the increase of metallicity. We observe a slight trend of decreasing $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at $n_{\mathrm{gas}}\gtrsim 10$ $\mathrm{cm}^{-3}$ because of coagulation. This effect of coagulation makes the 2175 Å bump and FUV slope less prominent. In contrast, shattering efficiently occurs in $0.1\lesssim n_{\mathrm{gas}}\lesssim 1$ $\mathrm{cm}^{-3}$, creating an enhanced small grain abundance in this density range. This as well as the enhanced accretion efficiency is the reason for a prominent 2175 Å bump and a steep FUV rise in the extinction curve in the lowest density range.

At 10 Gyr, the trend seen at 3 Gyr is enhanced. The extinction curve features are very strong in $0.1\lesssim n_{\mathrm{gas}}\lesssim 1$ $\mathrm{cm}^{-3}$ and less prominent in denser regions. This is due to the small grain formation by shattering as seen in the high values of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ between $n_{\mathrm{gas}}\sim 0.1$ and 1 $\mathrm{cm}^{-3}$. At the highest densities, $n_{\mathrm{gas}}\gtrsim 10$ $\mathrm{cm}^{-3}$, the extinction curve is as flat as that seen at 0.3 Gyr, but this flatness at 10 Gyr is due to coagulation, not the predominant stellar dust production. $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ decreases at $n_{\mathrm{gas}}\gtrsim 1$ $\mathrm{cm}^{-3}$: accretion is saturated because a large fraction of gas-phase metals are locked into dust, while coagulation is efficient.

3.3 Metallicity dependence

Metallicity could also be an important driver of dust evolution because of the coevolution of dust and metals (e.g. Rémy-Ruyer et al., 2014). To investigate the effect of metal enrichment, we show the metallicity dependence of extinction curves and small-to-large grain abundance ratios at various ages in Fig. 3. We derived mass-weighted extinction curves for gas particles that have the corresponding metallicity in each metallicity bin. The whole metallicity range considered is from $0.1Z_{\odot}$ to $2.5Z_{\odot}$.

At 0.3 Gyr, there are few particles with $Z$ > 1 $Z_{\odot}$ and the intensity of 2175 Å bump and the steepness of FUV rise increase with metallicity, although both of these features are weak and the extinction curves are overall flat. The $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ value has a rapid increase at $Z\gtrsim 0.1Z_{\odot}$; thus the features of extinction curves become prominent. The value of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at $Z\lesssim 0.01Z_{\odot}$ is determined by stellar dust production and shattering. Although there is no observational data at such low metallicity, we calculate the extinction curve for a metallicity range of $0.005\lesssim Z\lesssim 0.05Z_{\odot}$ to show prediction for an extremely metal-poor phase, which could be important for future high-redshift observations. The extinction curve in the lowest metallicity range is much flatter than the other metallicity ranges because of the lowest $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$. At $Z\gtrsim 0.01Z_{\odot}$, the efficiencies of shattering and accretion increase with metallicity. As a result, the small grain abundance increases drastically, especially at $Z\gtrsim 0.1Z_{\odot}$.

At 1 Gyr, the extinction curves of all the metallicity range show a clear 2175 Å bump and FUV rise because of the overall increase of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at all metallicities. The FUV slope is relatively flat at the lowest metallicity range, $0.1\lesssim Z\lesssim 0.2Z_{\odot}$. The $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ value does not increase any further when particles reach $Z\sim 0.5Z_{\odot}$, because the coagulation efficiency is comparable to the accretion efficiency at such high metallicities.

At 3 Gyr, the features of extinction curves are stronger in $0.1\lesssim Z\lesssim 0.5Z_{\odot}$ than in $0.5\lesssim Z\lesssim 2.5Z_{\odot}$. $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ shows a clear decline at $Z\gtrsim 0.5Z_{\odot}$ since coagulation is efficient enough to convert small grains to large grains. At $Z\sim 0.1Z_{\odot}$, $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ rises to the highest value in the entire metallicity range, because shattering dominates the grain size distribution around that metallicity. Thus, low-metallicity ($Z\sim 0.1$ $Z_{\sun}$) gas has steeper extinction curves than solar-metallicity gas. This trend at 3 Gyr is opposite to that seen above at 1 Gyr.

At 10 Gyr, particles with $0.2\lesssim Z\lesssim 1Z_{\odot}$ show a very strong 2175 Å  bump and steep FUV slope. Although the extinction curves of $Z\lesssim 0.2~{}\mathrm{Z}_{\odot}$ have relatively weak features, the 2175 Å  bump and FUV rise are still prominent. In this stage, shattering completely dominates the grain size distribution at $Z\lesssim 1Z_{\odot}$, raising $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ to the highest value seen at $Z\sim 1Z_{\odot}$. We note that $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ has a bimodal distribution especially for carbonaceous dust. As mentioned in Paper I, the higher-$\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ component, which makes the extinction curve features prominent, is associated with the diffuse medium in the galactic disc, while the lower component is associated with low density ($n_{\mathrm{gas}}<10^{-3}$ cm*-3*) and hot ($T_{\mathrm{gas}}>10^{4}$ K) gas more extended in the halo. This clear separation is due to stellar feedback and is further discussed in Section 4.4. Therefore, shattering in the diffuse ISM in the disc is important in the steepening of extinction curves in the latest ($\sim 10$ Gyr) evolutionary phase of the galactic disc.

4 discussion

4.1 Comparison with the Milky Way extinction curve

In Fig. 4, we compare our results at 3 and 10 Gyr with the observed mean Milky Way extinction curve and the 1$\sigma$ dispersion of the extinction curve in various lines of sight. The Milky Way extinction curves are determined by the stars at distances $\sim$0.1–2 kpc from the Sun (Fitzpatrick & Massa, 2007). Because spatially resolved information on a single line of sight (or on a scale below the path length) cannot be obtained, we need to smooth the spatially resolved extinction curves obtained in our simulation over a length comparable to the path length in the observations. However, because every observational line of sight (for every extinction curve) has different path length, it is difficult to apply a single smoothing length for comparison with observational data. Instead, we estimate a range of extinction curves predicted for various densities and radii as possible sources of the variation in extinction curves, and compare them with the observed curves and their variance.

We select 3 and 10 Gyr as two appropriate epochs because the representative metallicities of these epochs are close to that of the Milky Way. Our results are consistent with both the mean value and the dispersion except that the extinction curves at low density and in the outer disc region at 10 Gyr show a stronger 2175 Å  bump than the observation. Cardelli et al. (1989) and Fitzpatrick & Massa (2007) also show that the 2175 Å bump strength in the Milky Way extinction curves has a positive correlation with the FUV slope. This correlation can be interpreted as a sequence of small-to-large grain abundance ratio and is consistent with our results above. Reproducing the Milky Way extinction curve implies that we implemented all the processes that drive dust evolution in Milky-Way-like galaxies. We emphasize that we basically adopted the parameters adopted in Paper I without any fine-tuning specific for the extinction curve modelling.

Schlafly et al. (2016) recently showed that the Milky Way extinction curves are rather uniform across the entire disc, but that a significant variation in extinction curves is indeed observed. This variation occurs on scales much larger than the typical molecular-cloud size. As shown in Fig. 2, the large dispersion in $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at a given density shows that the grain size distribution (i.e. $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$) is not governed uniquely by the density. This is firstly because the current grain size distribution is determined by the past environments in which the dust has resided. The second reason is that the ISM density structure also changes especially by stellar feedback, which also disperses dust in a wide area. Therefore, the large-scale variation in extinction curves shown in the above study indicates the importance of treating dust processing in a consistent manner with the variation of ISM density structures.

We mention the following two caveats. First, most of the Milky Way extinction curves are taken from solar neighbourhood stars ($\lesssim 2$ kpc) that only reflects the local properties of dust around the Sun. Extinction curves in galactic centre and outer disc still remain uncertain because of lacking observation. Although there are some observational studies of extinction curves toward the centre of the Milky Way at near-infrared wavelengths (Nishiyama et al., 2006; Nishiyama et al., 2008; Nataf et al., 2013), it is hard to compare our results with them because our model is not able to predict precise extinctions at infrared wavelengths.

Second, we focus on dust evolution with a commonly used dust species model and succeed in reproducing observations. However, some previous efforts of modelling various dust species and grain size distributions suggest that there are some different dust species that fit the Milky Way extinction curves (Zubko et al., 2004; Jones et al., 2013). Although the trend that small grains produce steep and strong-feature extinction curves does not change against the change of dust materials, the effect of varying dust material properties is worth investigating in the future.

4.2 $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ in silicate and carbonaceous dust

Overall, silicate and carbonaceous dust share the same dust enrichment mechanisms. The differences between these two species lie in the condensation efficiency ($f_{\mathrm{{in,X}}}$), mass density ($\rho_{\mathrm{X}}$) and elemental abundance ($Z_{\mathrm{X\odot}}$; these quantities are listed in Table 1). Carbonaceous dust has a higher $f_{\mathrm{{in,X}}}$ than silicate so that stars produce more carbonaceous dust than silicate in our models. With a fixed dust mass, carbonaceous dust has higher shattering and coagulation efficiencies because a lower mass density means a larger number of dust grains.

At the beginning of galaxy evolution, the important processes are stellar dust production and shattering: Stars produce only large grains so that $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ is determined by shattering. From the above argument, carbonaceous dust has, on average, higher $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ than silicate at $t\lesssim 0.3$ Gyr (Figs. 1–3).

Equations (7) and (8) show that the elemental abundance ($Z_{\mathrm{X}}$) determines the efficiency of accretion: because we set the Milky Way elemental abundance pattern for the metal yields, accretion of silicate is more efficient than that of carbonaceous dust (Table 1; recall that Si occupies a mass fraction of 0.166 in silicate). Thus, when accretion dominates the dust abundance (at $0.3\lesssim t\lesssim 3$ Gyr), silicate has, on average, higher $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ than carbonaceous dust (Figs. 1–3).

Carbonaceous dust has higher shattering and coagulation efficiencies than silicate has as mentioned above, so in the late evolutionary stage ($t\gtrsim 3$ Gyr), carbonaceous dust shows a clearer decrease of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ by coagulation in the dense medium or in metal-rich regions while it presents higher $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ in the diffuse medium or sub-solar metallicity regions because of shattering (Figs. 1–3).

4.3 Dust species abundance ratio

The mass ratio of carbonaceous dust to silicate (C/Si) is important in producing various 2175 Å bump strengths under a similar FUV slope or vice versa. We adopted the classical graphite–silicate dust species model to estimate the extinction curve. In this model, the 2175 Å bump is contributed from small graphite grains while the FUV rise is mainly dominated by small silicates. Although the behaviour of extinction curve is broadly understood by the small-to-large grain abundance ratio, we here inspect C/Si to analyse further details of the extinction curve shapes.

In Fig. 5, we find that C/Si decreases with increasing age. In the early phase of galaxy evolution ($t\lesssim$ 0.3 Gyr), C/Si is determined by the dust condensation efficiency in stellar ejecta ($f_{\rm in,X}$) multiplied by the available elements composing the dust (C/Si $\sim$ 1.2). After accretion becomes efficient, C/Si is rather governed by the available gas-phase elements. Thus, in the late evolution stage ($t\sim$ 10 Gyr), C/Si approaches a value ($\sim$0.5) determined by the elemental abundance pattern assumed in the model (i.e. the solar abundance pattern). We note that the solar abundance pattern we adopted should be achieved as a result of contribution from various stellar populations including AGB stars (note that we did not include the contribution from AGB stars in the metal yield adopted in the simulation). Therefore, at ages younger than 1 Gyr, when less stars evolve into the AGB phase than at the age of the Milky Way ($\sim 10$ Gyr), adopting the solar metallicity pattern would overproduce the carbon abundance, which is raised significantly by AGB stars (Ferrarotti & Gail, 2006). Thus, the 2175 Å bump may be overestimated at the early evolutionary stage. The dust species abundance ratio is also important when we attempt to explain the SMC/LMC extinction curves in Section 4.5, where we mention a possible variation of C/Si again from a different point of view.

4.4 Extinction curves inside and outside the disc

To examine the impact of stellar and SN feedback on the extinction curve, we examine the distribution of various quantities by separating the gas spatially inside and outside of galactic disc by the following criteria: $r\leq 15$ kpc and $|z|\leq 2$ kpc for inside, and vice versa for outside, where $|z|$ is the distance perpendicular to the disc plane.

Before calculating the extinction curves, we produce a phase diagram at 10 Gyr (Fig. 6) to illustrate how stellar and SN feedback kick gas particles out of the disc and heat them up. Particles outside the disc occupy a low density ($n_{\mathrm{gas}}<10^{-2}$ cm*-3*) and high temperature ($T_{\mathrm{gas}}>10^{4}$ K) regions, as shown by the red contour. Gas inside the disc (blue contours) have higher densities and lower temperature, but there is also a plume of gas at $n_{\mathrm{gas}}\sim 1$ cm*-3* which is heated up to 106 K in the disc by feedback. The heated gas will try to escape from the disc and either rain down onto the disc to become cold gas again, or they will expand and adiabatically cool and join the red contour by moving to the left.

We examine the metallicity dependence of extinction curves again while this time separating inside and outside the disc components at 1 and 10 Gyr (Fig. 7). In general, the outer disc extinction curves have weaker features than the inner disc curves. At 1 Gyr, the particles inside the disc dominate the total particle number; thus, the extinction curves inside the disc are almost the same as those of all particles. Extinction curves of the outer disc show the weak 2175 Å bump and a flat FUV slope, and there is no particle outside the disc with $Z>1Z_{\odot}$. At 10 Gyr, except $Z>1Z_{\odot}$, the features of extinction curves inside the disc become more prominent than those in the entire galaxy, because the $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ values of carbonaceous dust concentrate on a higher value owing to accretion.

The $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ of silicate does not show a clear bimodal distribution at 10 Gyr for the following reasons: (i) Silicate has a higher accretion efficiency than carbonaceous dust (Table 1). Thus, silicate can reach higher $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ than carbonaceous dust before the particle is ejected out of the disc. (ii) Shattering efficiency is lower for silicate than for carbonaceous dust, so small silicate grains are not efficiently produced from large silicate grains in the disc. We also notice that $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ outside the disc evolves more slowly than that inside the disc. Once the particles are kicked out of the disc, only shattering occurs because of low density, and shattering is not efficient in such low density environment (see Eq. 5). Thus, the extinction curves outside the disc show weaker 2175 Å bumps and flatter FUV slopes compared to those inside the disc.

In Fig. 2, we calculated the extinction curves in the density range from 0.1 to $10^{3}$ cm*-3*, which automatically chose the particles inside the disc. The feedback effect is not apparent in the radial dependence of extinction curve shown in Fig. 1, because we considered only up to $r=8$ kpc, where the contribution from the particles outside the disc is negligible in the resulting extinction curves.

Because we are mainly interested in extinction curves in the galactic disc, we only considered supernova shocks as sources of dust destruction, but neglected sputtering in the diffuse hot gas which exists out of the disc. For more precise treatment of the extinction curves out of the disc, we need to include thermal sputtering. We expect that, if we include thermal sputtering, the extinction curves in the gas out of the disc would be further flattened because smaller grains are more efficiently destroyed by thermal sputtering (Hirashita et al., 2015). Thus, the tendency of flat extinction curves out of the galactic disc would be strengthened if we included thermal sputtering.

4.5 SMC/LMC extinction curves

In the above analysis, we have adopted graphite for the carbonaceous species. Small graphite grains have a strong 2175 Å bump, while some extinction curves in extragalactic objects show a much weaker bump than the Milky Way curve. We here examine a possibility of explaining the extinction curves in the SMC as a representative of bumpless extinction curves.

Nozawa et al. (2015) and Hou et al. (2016) adopted amorphous carbon instead of graphite to explain bumpless extinction curves in the nearby and distant Universe. Although the evolutionary link between amorphous carbon and graphite has not been clarified yet, it is worth examining if the bumpless SMC extinction curve can be reproduced with amorphous carbon.

There are other possible carbonaceous materials (Jones et al., 2013), but our two-size approach is not suitable for distinguishing detailed dust properties. Therefore, we just compare graphite and amorphous carbon in this paper to examine the two extremes in terms of the bump strength. In Fig. 8, we show the extinction curve at 1 Gyr, when the typical metallicity is $Z\sim 0.2Z_{\odot}$ (Fig. 3), with adopting amorphous carbon instead of graphite. Here, we adopt Zubko et al. (1996) for the optical constants for amorphous carbon (we adopt their ACAR). To keep consistency with the simulation, we still adopt the material density of graphite ($\rho_{\mathrm{C}}=2.24$ g cm*-3*) instead of that of amorphous carbon ($\rho_{\mathrm{C}}=1.81$ g cm*-3*) in Zubko et al. (2004) for carbonaceous dust. Because the difference in the material density between those two carbonaceous species is not large, the following results are hardly affected by the choice of the material density. Fig. 8 shows that, although the 2175 Å bump is eliminated successfully, there is a large discrepancy toward FUV wavelengths. This discrepancy is resolved only if the relative contribution of small silicate grains is increased in our model. Weingartner & Draine (2001) suggested a smaller graphite-to-silicate mass ratio than that of the Milky Way to explain the SMC extinction curve, and Pei (1992) fitted the SMC extinction curve exclusively with silicate. However, even the values of C/Si at 10 Gyr, when C/Si is the lowest (Fig. 5), cannot reproduce the SMC extinction curve (light grey lines in Fig. 8). This means that our simulation has difficulty in reproducing such a steep FUV rise as observed in the SMC extinction curve.

There have been some ideas proposed to reproduce the SMC extinction curve. Bekki et al. (2015) suggested that small carbonaceous dust grains are selectively transported outside the SMC by the most recent starburst. In their model, decoupling of grain motion from gas motion is important; however, our simulation assumes a complete dynamical coupling between dust and gas. This decoupling could be included in a scheme developed by Bekki (2015), who proposed an SPH modelling in which dust and gas are treated by different particle species. Thus, we will leave the possibility of selective transport of small carbonaceous grains to the future work. Hou et al. (2016) proposed more efficient SN destruction of small carbonaceous grains than that of silicate and showed that the extinction curve could become as steep as the SMC curve because of a resulting high fraction of small silicate grains.

The above previous ideas of reproducing the SMC extinction curve are based on selective elimination of small carbonaceous dust. Thus, we also examine an extreme case by assuming that the extinction curve is only contributed from silicate, i.e. we eliminate carbonaceous dust (red thin lines in Fig. 8). The results show too steep FUV slopes, which implies that a small fraction of carbonaceous dust is still necessary to reproduce the SMC extinction curve. The fact that the observed SMC extinction curve is between the pure silicate curve and our predictions supports the idea that a certain fraction of carbonaceous dust is lost by certain mechanisms such as outflow and SN destruction.

In addition, we also attempt to predict the LMC extinction curve by changing the graphite-to-amorphous carbon ratio and selecting snapshots with a higher metallicity than the SMC. We adopt the snapshot at 3 Gyr and assume that the graphite-to-amorphous carbon ratio is 1:1. The resulting extinction curves are shown in Fig. 9. In the figure, the 2175 Å bump strength is roughly reproduced in $n_{\mathrm{gas}}$ = 0.1–1 cm*-3* and $r=6$–8 kpc. Thus, the weak 2175 Å bump in the LMC compared to the Milky Way can be represented by the mixture of graphite and amorphous carbon. However, a discrepancy remains at FUV wavelengths as was also the case for the SMC. This implies that reducing C/Si is also necessary to reproduce the FUV slope of the LMC extinction curve.

4.6 Extinction curves in high redshift galaxies

At high redshift, very bright point sources such as quasars and gamma-ray bursts provide opportunities to derive the extinction curves in their host galaxies, taking advantage of their simple power-law-like SEDs at rest optical and UV wavelengths (Maiolino et al., 2004; Stratta et al., 2007; Li et al., 2008; Elíasdóttir et al., 2009; Gallerani et al., 2010; Perley et al., 2010; Zafar et al., 2010; Zafar et al., 2011; Zafar et al., 2015; Schady et al., 2012). Here we examine whether our simulation could also produce extinction curves in the distant Universe. In Fig. 10, we show the extinction curves at 1 Gyr (cosmic age at $z\sim 5$) using amorphous carbon for the carbonaceous dust species. The extinction curves are normalized to the value at rest 0.3 $\micron$ to compare with the extinction curve of a quasar, SDSS J1048+4637, at $z=6.2$ (Maiolino et al., 2004) as a representative extinction curve. As shown by Gallerani et al. (2010), the extinction curves in quasars at $z>5$ are similar in the sense that they show much flatter bumpless extinction curves than the SMC curve (but see Hjorth et al., 2013, for observational uncertainties). Our results are in good agreement with the observational data and are almost the same as the prediction by Nozawa et al. (2015), who successfully fitted the same extinction curve by using their single-zone dust evolution model but by taking account of the full grain size distribution (i.e. without the two-size approximation).

Since we have reproduced both the high-redshift extinction curve and the Milky Way extinction curve, there may exist an evolutionary link between these two extinction curves (or extinction curves in various redshifts in general). Li et al. (2016) examined extinction curves of star-forming disc galaxies at $0.5<z<2.5$ and found that extinction curves tend to be flatter at higher redshift. However, we note that the extinction curves (more precisely attenuation curves) derived by them represent effective extinction as a result of radiative transfer (absorption and scattering) effects of light emitted from multiple stellar populations (see also Calzetti et al., 1994; Inoue, 2005; Mancini et al., 2016). Radiative transfer effects, especially scattering, skew the shape of the extinction curve in highly extinguished objects (e.g. Scicluna & Siebenmorgen, 2015). Kriek & Conroy (2013) analysed SEDs to derive attenuation curves for galaxies at $0.3<z<2$. They found that there is a positive correlation between 2175 Å  bump strength and FUV slope. This correlation can be interpreted as a sequence of small-to-large grain abundance ratio. They also showed that galaxies with higher sSFR tend to have flatter attenuation curve. In this work, the radiative transfer calculation is not included yet, but simulation results tend to predict flat extinction curves at the early evolutionary stage (or in the phase of high sSFR; see Sec 5.1 of Paper I), which is similar to the above observational tendency of attenuation curves with sSFR.

5 conclusion

In this work, we investigated the evolution of extinction curve using the SPH simulation in which dust evolution is implemented. As discussed in Paper I, our simulation adopts the two-size model of dust by Hirashita (2015), and solves the formation and destruction of dust as a function of time together with gas dynamics. Additionally, dust species are separated into carbonaceous dust and silicate, and the two species are assumed to evolve independently. They share the same dust formation and destruction mechanisms but have different condensation efficiencies, material properties and elemental abundances.

We compare our results with observed galaxy extinction curves in the local Universe. To estimate the extinction curves, we adopt the optical constants of graphite and silicate based on Draine & Lee (1984) and Weingartner & Draine (2001). The simulation enables us to examine the dependence of extinction curves on the position, gas density, and metallicity in the galaxy as well as the time evolution.

Our results show that extinction curves are flat at the earliest evolutionary stage $t\lesssim$ 0.3 Gyr, because stellar sources dominate the dust production at that epoch (stars only produce large grains in our model). Dust growth by accretion increases the small-to-large grain abundance ratio $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ at $t\gtrsim$ 0.3 Gyr, making the 2175 Å bump and FUV rise more prominent. Coagulation becomes efficient after small grains become abundant, especially in the central region where the gas is dense and metal-rich. Because of coagulation, the 2175 Å bump and FUV rise are suppressed at $t\gtrsim$ 3 Gyr in dense regions. At $t\sim$ 10 Gyr, shattering dominates the dust evolution in low-density, low-metallicity, or outer regions of the galactic disc. Because shattering continuously converts large grains to small grains, the 2175 Å bump and FUV rise are strong in those regions.

To further constrain the grain size distribution, we make a comparison with observed extinction curves and find that the predicted extinction curves at $t\gtrsim 3$ Gyr are consistent with the mean Milky Way extinction curve and the dispersion of extinction curves toward different lines of sight. This implies that all dust processes that drive dust evolution in Milky-Way-like galaxies are implemented successfully in our simulation.

We also examine the difference in $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ between carbonaceous dust and silicate. The two dust species have different condensation efficiencies, mass densities and elemental abundances. At the stellar production dominated epoch ($t\lesssim 0.3$ Gyr), carbonaceous dust has higher $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ than silicate because the shattering efficiency of carbonaceous dust is higher than that of silicate. In the period in which accretion dominates the dust abundance ($0.3\lesssim t\lesssim 3$ Gyr), $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ of silicate is, on average, higher than that of carbonaceous dust because under the solar elemental abundance pattern we assumed, silicate has more gas-phase elements available for accretion than carbonaceous dust. After 3 Gyr, carbonaceous dust presents higher $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ in the diffuse medium or sub-solar metallicity regions and lower $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ in the dense medium or solar metallicity regions than silicate since the efficiencies of shattering and coagulation of carbonaceous dust are higher than those of silicate. The abundance ratio of carbonaceous dust to silicate, C/Si, is also affected by different efficiencies of enrichment mechanisms. Carbonaceous dust is more abundant than silicate in the stellar production dominated epoch ($t\lesssim 0.3$ Gyr) because carbonaceous dust has higher dust condensation efficiency in stellar ejecta than silicate has; after the accretion dominated epoch ($t\gtrsim 1$ Gyr), silicate becomes more abundant than carbonaceous dust because of the higher abundance of available gas-phase elements for accretion.

There are two sequences in the metallicity dependence of $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ corresponding to the inside and outside disc components. The outside disc component is transported there mainly by stellar and SN feedback. The extinction curves of the outside disc component show shallower FUV slopes and weaker 2175 Å bumps than those of the inside disc component because $\mathcal{D}_{\mathrm{S}}/\mathcal{D}_{\mathrm{L}}$ evolves more slowly outside the disc than inside. This indicates that the extinction curves in the circum-galactic medium are flatter than those in galactic discs.

We also attempt to reproduce the SMC extinction curve by adopting amorphous carbon instead of graphite. Although 2175 Å bump is successfully eliminated, there is a large discrepancy at FUV wavelengths. The steep FUV rise in the observed SMC extinction curve is reproduced only by a higher silicate fraction than calculated in this paper, which indicates that carbonaceous dust is selectively lost by, e.g., outflow and SN destruction. In addition, reproducing the LMC extinction curve by changing graphite-amorphous carbon ratio is not fully successful because we underpredict the FUV rise also for the LMC.

Finally, we also reproduced a representative extinction curve of high-redshift quasars using our simulation results at 1 Gyr. The observational trend that extinction curves are flatter for higher specific star formation rates is also interpreted by our model, in the sense that our model also predicts the change from flat to steep extinction curves as the specific star formation rate drops.

Acknowledgements

We thank Y.-H. Chu, T.-H Chiueh, and W.-H. Wang for useful comments. HH is supported by the Ministry of Science and Technology grant MOST 105-2112-M-001-027-MY3. KN, SA and IS acknowledges the support from JSPS KAKENHI Grant Number 26247022. Numerical simulations were in part carried out on the XC30 at the Centre for Computational Astrophysics, National Astronomical Observatory of Japan.