# A new galaxy spectral energy distribution model consistent with the
evolution of dust

Kazuki Y. Nishida,1 Tsutomu T. Takeuchi,1,2 Takuma Nagata,1 and Ryosuke S.
Asano1
1Division of Particle and Astrophysical Science, Nagoya University, Furo-cho,
Chikusa-ku, Nagoya, 464–8602, Japan
2The Research Center for Statistical Machine Learning, The Institute of
Statistical Mathematics, 10-3 Midori-cho, Tachikawa, Tokyo 190–8562, Japan
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The spectral energy distribution (SED) of galaxies provides fundamental
information on the related physical processes. However, the SED is
significantly affected by dust in its interstellar medium. Dust is mainly
produced by asymptotic giant branch stars and Type II supernovae. In addition,
the dust mass increases through the metal accretion, and the grain size
changes by the collisions between the grains. The contribution of each process
and the extinction depend on the size distribution. Therefore, the SED model
should treat the evolution of the dust mass and size distribution. In spite of
the importance of dust evolution, many previous SED models have not considered
the evolution of the total mass and size distribution in a physically
consistent manner. In this work, we constructed a new radiative transfer SED
model, based on our dust evolution model consistent with the chemical
evolution. To reduce the computational cost, we adopted the mega-grain and the
one-dimensional plane parallel galaxy approximation. As a fiducial case, we
calculated Milky Way-like galaxy SEDs at various ages under the closed-box
model. We found that a galaxy at the age of 100 Myr does not produce small
grains such as polycyclic aromatic hydrocarbons. After 1 Gyr, we observed a
drastic increase of infrared emission and attenuation caused by a rapid
increase of dust mass. This phenomenon can be treated appropriately for the
first time by our new model. This model can be used for the SED fitting to a
galaxy at any stage of evolution.

###### keywords:

dust, extinction – galaxies: evolution – ISM: evolution – radiative transfer –
galaxies: ISM – galaxies: disc

††pubyear: 2020††pagerange: A new galaxy spectral energy distribution model
consistent with the evolution of dust–A new galaxy spectral energy
distribution model consistent with the evolution of dust

##  1 Introduction

The spectral energy distribution (SED) fitting is a fundamental method to
extract the information of the physical processes in galaxies (e.g., star
formation rate: SFR, stellar mass, dust mass) from observational data. Stars
emit photons with wavelengths ranging from ultraviolet (UV) to near-infrared
(NIR). Dust grains absorb and scatter the photons emitted from stars, and re-
emit the absorbed energy at mid-infrared (MIR) to far-infrared (FIR). In
addition to the radiative aspect of a galaxy, dust grains promote the
formation of hydrogen molecules on the surface of the grains (e.g., Hollenbach
& McKee, 1979; Hirashita & Ferrara, 2002; Cazaux et al., 2005). Since hydrogen
molecules are one of the fundamental ingredients of the star formation, dust
grains directly activate the star formation in galaxies.

Dust is a solid grain consisting of the elements heavier than helium, and is
produced by stellar mass loss and supernovae (SNe). Outflows from low- and
intermediate-mass stars during the thermally pulsing asymptotic giants branch
(TP-AGB) phase and Type II SNe (SNe II) are considered to be the primary
sources of dust (e.g., Nozawa et al., 2007; Bianchi & Schneider, 2007;
Zhukovska et al., 2008). Dust grains are formed by condensation of heavy
elements in the atmosphere of massive stars, and only the grains that could
survive the reverse shock of the SN are finally expelled into the interstellar
medium (ISM). Then, the blast waves from SNe propagating in the ISM also
destroy dust grains (e.g., Jones et al., 1994, 1996; Nozawa et al., 2003;
Nozawa et al., 2006; Zhukovska et al., 2008; Yamasawa et al., 2011). This
destruction process has been confirmed by observations of several supernova
remnants (e.g., Borkowski et al., 2006; Arendt et al., 2010). Details of the
dust grain survival still remain controversial. It might depend on their
composition and size111Terminologies such as grain size, size distribution,
etc. are often used in articles related to dust. Throughout this paper, when
we mention “size” of dust grain, it always mean the dust grain radius under
the assumption of spherical shape. (e.g., Nozawa et al., 2007; Gall et al.,
2014; Slavin et al., 2020), and others claim that it depends on the clumpiness
of the ejecta (e.g., Biscaro & Cherchneff, 2016). In addition, Matsuura et al.
(2019) argue that the dust destruction by the SN are suppressed because the
atoms can stick to the surviving dust grain in the passage of the forward
shock region and it can reform or increase dust grains. Observations of SN
remnants (SNRs) do not give a final answer, since dust mass and composition
are significantly different among observed SNRs.

In addition to the dust production from stars, dust growth in the ISM is also
an important process, and necessary to explain the large amount of dust
present in galaxies (e.g., Asano et al., 2013a; Zhukovska, 2014; Michałowski,
2015; Lesniewska & Michałowski, 2019). In the cold phase of ISM, metal is
accreted onto the dust grain surface to increased the size and total mass of
the dust (e.g., Dwek, 1998; Zhukovska et al., 2008; Michałowski et al., 2010;
Hirashita & Kuo, 2011; Asano et al., 2013a). The size distribution of the dust
is also changed through the collisions between dust grains (e.g., Yan et al.,
2004; Jones et al., 1996; Hirashita & Yan, 2009; Kuo & Hirashita, 2012).

Which process controls the mass of dust varies greatly, depending on the age
and the environment of the galaxy, and it is still actively debated from
several different points of view. The SNe II dominates the dust production,
especially in very young galaxies, because SNe II have a shorter lifetime
($<30$ Myr) than AGB stars ($>150$ Mry) (e.g., Morgan & Edmunds, 2003;
Marchenko, 2006; Dwek et al., 2007; Valiante et al., 2009; Gall et al., 2011a,
b; Liu & Hirashita, 2019; De Looze et al., 2020; Burgarella et al., 2020;
Nanni et al., 2020). We should note, however, that the contribution of AGB
cannot be ignored even in galaxies with the age of 500 Myr, when the star
formation rate (SFR) is high (Valiante et al., 2009).

The debate of dust grain growth in high-$z$ galaxies has not been settled. In
high-$z$ galaxies, several studies claim that the process is not effective
because there is not enough time to growth, low gas density, and high
temperature (e.g., Ferrara et al., 2016; Ceccarelli et al., 2018). Since dust
growth in the ISM is strongly affected by the metallicity in a galaxy (e.g.,
Inoue, 2003; Asano et al., 2013a), they claim that it might not be very
important in young galaxies with low metallicity. However, other studies have
shown that the dust mass in distant galaxies cannot be explained without
considering metal accretion (e.g., Pipino et al., 2011; Valiante et al., 2011;
Zhukovska, 2014; Michałowski, 2015; Mancini et al., 2015; Lesniewska &
Michałowski, 2019; Rouillé et al., 2020).

Asano et al. (2013a) defined the critical metallicity, $Z_{\mathrm{cr}}$, as
the metallicity of the ISM, at which the production rate of dust from stars
(AGB and SNe II) become equal to the mass growth rate in the ISM. When the
metallicity reaches $Z_{\mathrm{cr}}$, the dust mass increases suddenly and
nonlinearly. This rapid increase in dust mass is caused by the following
process (e.g., Hirashita & Yan, 2009; Asano et al., 2013a, b). First, the
metal accretion depends on the metallicity and total surface area of the dust
grains. As the dust size increases through the metal accretion, yet another
process of dust evolution in the ISM, shattering, is more likely to occur.
This process is basically the collision of grains with each other and
redistribute the mass of dust into smaller-sized grains. When the shattering
becomes effective, the total surface area of the dust grain per mass
increases, and the metal accretion becomes more efficient. This cycle leads to
the sharp increase of the total dust mass along with metallicity. Therefore,
not only the total dust mass, but also it is of vital importance to take into
account the grain size distribution to discuss the evolution of dust in
galaxies.

The absorption and scattering coefficients of dust as a function of wavelength
depend on the size and composition of grains. When a dust grain absorbs light,
its temperature rises and the absorbed energy is re-emitted at longer
wavelengths (mainly IR). The wavelength of the re-emission depends on the
instantaneous temperature of the grain, and the temperature strongly depends
on the size of the grain. Thus, as already mentioned, the mass, size, and
composition of dust grains play a fundamental role in shaping the SED of a
galaxy. Fitting to the SED of distant galaxies by an empirical dust emission
model without evolution may lead to erroneous results. This happens, for
example, when we use a model in which the dust size distribution is constant.
Recently, some galaxies with a large amount of dust
($M_{\mathrm{dust}}>10^{6}~{}M_{\odot}$) have been observed at $z>6$ (e.g.,
Watson et al., 2015; Laporte et al., 2017; Tamura et al., 2019). Now it is a
proper moment to develop a new SED model based on the theory of dust
evolution, after the advent of the Atacama Large Millimeter/Submillimeter
Array (ALMA) and other large facilities at this wavelength range.

The dust evolution in the ISM has been considered by a number of previous
studies (e.g., Dwek, 1998; Calura et al., 2008; da Cunha et al., 2010; Asano
et al., 2013a, b, 2014; Mancini et al., 2015; Schneider et al., 2016; Ginolfi
et al., 2018; De Vis et al., 2017, 2019; Hirashita & Aoyama, 2019; De Looze et
al., 2020; Nanni et al., 2020; Burgarella et al., 2020). We put our basis on
the theoretical framework of dust evolution proposed by Asano et al. (2013a,
b, 2014) and Nozawa et al. (2015) (hereafter Asano model) to develop a new
radiative transfer SED model. The Asano model considers SNe II and AGB stars
as dust production sources, and not only the metal accretion but also
shattering and coagulation as the dust evolution process in the ISM, which
enables us to determine the dust mass and size distribution in all galaxy ages
from a first principle. Hirashita & Aoyama (2019) developed a dust evolution
model also based on the Asano model, but with a better computational
performance to apply to cosmological simulations. Recall that the Asano model
considers different physical quantities (e.g., ambient gas density of SN,
hydrogen gas density, and magnetic field strength) for various galaxies to
treat the dust destruction by SN and dust collisions. However,considering dust
on the cosmological scale, it is impossible to reach a galaxy scale
resolution. Therefore, Hirashita & Aoyama (2019) adopt many simplifications to
optimize their model for the simulations. In contrast, since we aim at
calculating SED of an individual galaxy, we make a maximal use of the Asano
model.

There have been several SED models that include the evolution of dust in the
ISM. For example, Schurer et al. (2009) is the SED model which is based on the
dust model of Calura et al. (2008). They calculate the chemical evolution in a
single gas phase and dust evolution including metal accretion. However, since
Schurer et al. (2009) do not consider shattering and coagulation, the rapid
increase of the total dust mass does not occur. Version 3 of Pégase(Fioc &
Rocca-Volmerange, 2019) considers the evolution of dust mass, and can
calculate not only the radiation from stars but also the extinction by dust
grain and the radiation of dust with a stochastic temperature distribution.
Pégasedoes not take into account the dust size distribution, and assumes that
the fraction and the size distribution of each grain species do not evolve. In
this paper, we construct a new galaxy SED model, including the dust evolution
theory proposed by Asano et al. (2013a, b, 2014). We adopt the mega-grain
approximation (MGA) with a one-dimensional plain parallel galaxy (Varosi &
Dwek, 1999; Inoue, 2005) to make the radiative transfer calculation faster.

This paper is organized as follows. In Section 2, we introduce how to
calculate the SED for each component. In Section 3, as an example of our SED
model, we show the SED of a Milky Way (MW)-like galaxy. In Section 4, we
discuss the effect of parameters on the model SEDs. Section 5 is devoted to
the conclusions.

##  2 Methods: Construction of SED model

To construct a galaxy SED model, we synthesis the stellar SED calculated by
version 2 of Pégase (Fioc & Rocca-Volmerange, 1999, Pégase.2), the dust
evolution model based on Asano et al. (2013a, b, 2014), dust attenuation
calculated by radiative transfer with MGA in a one-dimensional galaxy (Varosi
& Dwek, 1999; Inoue, 2005), and the dust emission by a Monte Carlo simulation.
In this section, we present how to calculate each component. §2.1 introduces
the equation of mass evolution in a galaxy. §2.2 shows the details of the dust
chemical evolution model. §2.3 is the overview of calculation of stellar SED
by Pégase.2. §2.4 and 2.5 show dust properties with mega-grain approximation
and radiative transfer in one-dimensional galaxy, respectively. In §2.6 and
2.7, we explain the calculation of the dust temperature distribution by Monte
Carlo simulation and the dust emission by the distribution.

###  2.1 Equations governing galaxy evolution

We consider stars, gases and dust grains as the components of a model galaxy.
For simplicity, we assume a one-zone galaxy model, where physical quantities
vary uniformly over the entire galaxy. The time evolution of the total stellar
mass $M_{\ast}(t)$, the ISM mass $M_{\mathrm{ISM}}(t)$, the metal mass
$M_{\mathrm{Z}}(t)$, and the dust mass $M_{\mathrm{d}}(t)$ at an age of galaxy
$t$ are represented as (Lisenfeld & Ferrara, 1998; Asano et al., 2013a),

$\displaystyle\frac{\mathrm{d}{M_{\ast}(t)}}{\mathrm{d}{t}}$
$\displaystyle=\mathrm{SFR}(t)-R(t),$ (1)
$\displaystyle\frac{\mathrm{d}{M_{\mathrm{ISM}}(t)}}{\mathrm{d}{t}}$
$\displaystyle=-\mathrm{SFR}(t)+R(t)+\frac{\mathrm{d}{M_{\mathrm{infall}}(t)}}{\mathrm{d}{t}},$
(2) $\displaystyle\frac{\mathrm{d}{M_{\mathrm{Z}}(t)}}{\mathrm{d}{t}}$
$\displaystyle=-Z(t)\mathrm{SFR}(t)+R_{\mathrm{Z}}(t)+Y_{\mathrm{Z}}(t),$ (3)
$\displaystyle\frac{\mathrm{d}{M_{\mathrm{d}}(t)}}{\mathrm{d}{t}}$
$\displaystyle=-D(t)\mathrm{SFR}(t)+Y_{\mathrm{d}}(t)$
$\displaystyle-\Big{(}\frac{\mathrm{d}M_{\mathrm{d}}(t)}{\mathrm{d}t}\Big{)}_{\mathrm{SN}}+\Big{(}\frac{\mathrm{d}M_{\mathrm{d}}(t)}{\mathrm{d}t}\Big{)}_{\mathrm{acc}},$
(4)

where $\mathrm{SFR}(t)$ is the star formation rate, and $R(t)$ and
$R_{\mathrm{Z}}(t)$ are the mass of the gas and metal taken into stars from
ISM and returned to ISM when stars die per unit time, respectively.
$\mathrm{d}M_{\mathrm{infall}}/\mathrm{d}t$ is the infall gas rate, which is
assumed to be zero in this paper except §4.5. $Z(t)\equiv
M_{\mathrm{Z}}/M_{\mathrm{ISM}}$ is the metallicity, and $D(t)\equiv
M_{\mathrm{d}}/M_{\mathrm{ISM}}$ is the mass fraction of dust with respect to
the total amount of metal. $Y_{\mathrm{Z}}(t)$ and $Y_{\mathrm{d}}(t)$ are the
metal and dust masses newly produced by stars per unit time, respectively.
$(\mathrm{d}M_{\mathrm{d}}(t)/\mathrm{d}t)_{\mathrm{SN}}$ and
$(\mathrm{d}M_{\mathrm{d}}(t)/\mathrm{d}t)_{\mathrm{acc}}$ is the change of
grain mass cased by SN shock and metal accretion. In the Asano model, three
phases are supposed in the ISM: warm neutral medium (WNM, with gas temperature
$T_{\mathrm{gas}}=6000$ K, and hydrogen number density
$n_{\mathrm{H}}=0.3~{}\mathrm{cm^{-3}}$), cold neutral medium (CNM, with
$T_{\mathrm{gas}}=100$ K, $n_{\mathrm{H}}=30~{}\mathrm{cm^{-3}}$), and
molecular cloud (MC, with $T_{\mathrm{gas}}={25}$ K,
$n_{\mathrm{H}}=300~{}\mathrm{cm^{-3}}$) (Nozawa et al., 2015). In the MC,
dust grains form icy mantle on the grain surface. (e.g., Kalvans, 2017;
Ceccarelli et al., 2018). However, since the properties of the icy mantle are
not well understood yet, we do not consider its effect in this work. As for
the dust growth in the ISM, metal acceleration occurs only in CNM and MC, and
shattering and coagulation occur in all three phases. Because metal accretion
occurs effectively in high density regions, the grain growth in the MC is more
prominent. In this paper, we fix the phase fraction of WNM, CNM, and MC to
$\eta_{\mathrm{WNM}}=0.5$, $\eta_{\mathrm{CNM}}=0.3$, and
$\eta_{\mathrm{MC}}=0.2$, respectively, the same values to that of Nozawa et
al. (2015) used to reproduce the MW extinction curve. These fractions are
constant throughout the calculation of the dust model. For each time step, the
dust grain is redistributed into each ISM phase so that the mass fraction of
each ISM phase is maintained. Thus, we consider the dust grain cycling between
the different ISM phases. Also, we do not consider the outflow effects in this
model. We assume that at the age of $t=0$, the galaxy contains no stars and
dust, and contains only zero-metallicity gas (i.e.,
$M_{\mathrm{star}}(0)=M_{\mathrm{Z}}(0)=M_{\mathrm{d}}(0)=0$, and
$M_{\mathrm{ISM}}(0)$ is total galaxy mass).

We adopt the Schmidt law (Schmidt, 1959), $\mathrm{SFR}(t)\propto
M^{n}_{\mathrm{ISM}}$ for SFR with $n=1$ for simplicity, as

$\mathrm{SFR}(t)=\frac{M_{\mathrm{ISM}}(t)}{\tau_{\mathrm{SF}}},$ (5)

where $\tau_{\mathrm{SF}}$ is the timescale of star formation. In this paper,
initial galaxy mass $M_{\mathrm{ISM}}(0)$ and $\tau_{\mathrm{SF}}$ are set to
be $10^{11}~{}\mathrm{M_{\odot}}$ and $3$ Gyr as fiducial values. $R(t)$,
$R_{\mathrm{Z}}(t)$, and $Y_{\mathrm{d}}(t)$ are represented as

$\displaystyle R(t)$
$\displaystyle=\int^{100~{}M_{\odot}}_{m_{\mathrm{min}}(t)}[m-\omega(m,Z(t-\tau_{m}))]\phi(m)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\mathrm{SFR}(t-\tau_{m})\,\mathrm{d}m,$
(6) $\displaystyle R_{\mathrm{Z}}(t)$
$\displaystyle=\int^{100~{}M_{\odot}}_{m_{\mathrm{min}}(t)}[m-\omega(m,Z(t-\tau_{m}))]\phi(m)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\mathrm{SFR}(t-\tau_{m})Z(t-\tau_{m})\,\mathrm{d}m,$
(7) $\displaystyle Y_{\mathrm{Z}}(t)$
$\displaystyle=\int^{100~{}M_{\odot}}_{m_{\mathrm{min}}(t)}m_{\mathrm{Z}}(m,Z(t-\tau_{m}))\phi(m)\mathrm{SFR}(t-\tau_{m})\,\mathrm{d}m,$
(8) $\displaystyle Y_{\mathrm{d}}(t)$
$\displaystyle=\int^{100~{}M_{\odot}}_{m_{\mathrm{min}}(t)}m_{\mathrm{d}}(m,Z(t-\tau_{m}))\phi(m)\mathrm{SFR}(t-\tau_{m})\,\mathrm{d}m,$
(9)

where $m_{\mathrm{min}}(t)$ is the lower limit mass of star which can explode
at time $t$, $\phi(t)$ is the initial mass function (IMF),
$\omega(m,Z(t-\tau_{m}))$, $m_{\mathrm{Z}}(m,Z(t-\tau_{m}))$, and
$m_{\mathrm{d}}(m,Z(t-\tau_{m}))$ are the remnant mass which remains after a
star explodes, and the metal mass and the dust mass newly produced by a star
of mass $m$ and metallicity $Z(t-\tau_{m})$. As for $\omega$ and
$m_{\mathrm{Z}}$, we adopt Ventura et al. (2013) for AGB stars with mass
$m=1$–$8~{}M_{\odot}$ and metallicity $Z=(0.015$, 0.4) $Z_{\odot}$, and
Kobayashi et al. (2006) for SNe II with progenitor mass
$m=13$–$40~{}M_{\odot}$ and metallicity $Z=(0.0$, 0.05, 0.3, 1.0) $Z_{\odot}$.
We interpolate and extrapolate all data tables over mass and metallicity in
this paper. $\tau_{m}$ is the lifetime of a star with mass $m$ and we use the
following equation by Raiteri et al. (1996),

$\log\tau_{m}=a_{0}(Z)+a_{1}(Z)\log m+a_{2}(Z)(\log m)^{2},$ (10)

with

$\displaystyle a_{0}(Z)$ $\displaystyle=10.13+0.07547\log Z-0.008084(\log
Z)^{2},$ (11) $\displaystyle a_{1}(Z)$ $\displaystyle=-4.424-0.7939\log
Z-0.1187(\log Z)^{2},$ (12) $\displaystyle a_{2}(Z)$
$\displaystyle=1.262+0.3385\log Z+0.05417(\log Z)^{2}.$ (13)

This equation was obtained by fitting the calculation result of stars with
stellar mass range 0.6–120 $M_{\odot}$ and metallicity range 0.0004–0.05 by
the Padova group (Alongi et al., 1993; Bressan et al., 1993; Bertelli et al.,
1994). We use the Salpeter IMF (Salpeter, 1955):

$\phi(m)\propto m^{-2.35}.$ (14)

The IMF is normalized as

$\int^{100~{}M_{\odot}}_{0.1~{}M_{\odot}}\phi(m)m~{}\mathrm{d}m=1~{}M_{\odot}.$
(15)

###  2.2 Dust evolution model

In this paper, we adopt the Asano model. The Asano model takes into account
ten species for dust; C, Si, SiO2, Fe, FeS, Al2O3, MgO, MgSiO3, Mg2SiO4, and
Fe3O4 (Nozawa et al., 2007; Zhukovska et al., 2008). Since the Asano model
calculates the evolution for each dust species, the dust composition evolves
with time. For simplicity in this work, we divide dust grains into two
representative families, silicate and carbonaceous grains, for the grain
growth and grain-grain collision, because the optical properties of other dust
species are not well understood yet. Further, among the carbonaceous grains,
the smaller ones are treated as polycyclic aromatic hydrocarbon (PAH) grains.
The fraction of graphite in carbonaceous dust grain is obtained by the
following formula (Draine & Li, 2007),

$f_{\mathrm{gra}}=\begin{cases}0.01&(a<50~{}\mathrm{\mathring{A}})\\\
0.01+0.99\left[1-\left(\frac{50~{}\mathrm{\mathring{A}}}{a}\right)^{3}\right]&(a>50~{}\mathrm{\mathring{A}})\end{cases}.$
(16)

Where $a$ is the dust grain radius and the fraction of PAHs is defined as
$f_{\mathrm{PAH}}=1-f_{\mathrm{gra}}$. Since the PAHs are divided into ionized
and neutral PAHs, which have different optical properties. The fraction of
ionized PAH is shown in Figure 7 of Li & Draine (2001). A part of carbonaceous
grains could consist of amorphous carbon. According to Nozawa et al. (2015),
graphite is accepted for reproducing the attenuation curve of nearby galaxy
like MW, and amorphous carbon is accepted for galaxies that do not have 2175 Å
bumps in the attenuation curve, such as high-$z$ quasars. Since the
differences between these compositions are not well understood yet, we only
consider graphite grains in this work. It is also possible to take into
account amorphous carbon in our SED model.

This model considers only AGB and SNe II as the sources of dust grains. For
simplicity, we do not consider SNe Ia contributions because they are
considered to be a very minor contributor to the total mass of dust (e.g.,
Calura et al., 2008).

The dust size distribution is represented by the dust number
$f_{\mathrm{X}}(a,t)$ and mass $\rho_{\mathrm{X}}(m_{\mathrm{d}},t)$
distribution. $f_{\mathrm{X}}(a,t)\,\mathrm{d}a$ and
$\rho_{\mathrm{X}}(m_{\mathrm{d}},t)\,\mathrm{d}m_{\mathrm{d}}$ are the number
and mass density of dust grains with radii [$a,a+\mathrm{d}a$] and mass
[$m_{\mathrm{d}},m_{\mathrm{d}}+\mathrm{d}m_{\mathrm{d}}$] at time $t$,
respectively. Where ’X’ represents the dust species (C: carbonaceous or Si:
silicate grain), $a$ is the dust radius, and $m_{\mathrm{d}}$ is the dust
grain mass. We assume that the dust grain has a constant density $s$ and is
spherical grain, so dust grain mass is

$m_{\mathrm{d}}=\frac{4}{3}\pi a^{3}s.$ (17)

The relation of dust number and mass density is expressed as

$\rho_{\mathrm{X}}(m_{\mathrm{d}},t)\,\mathrm{d}m_{\mathrm{d}}=m_{\mathrm{d}}f_{\mathrm{X}}(a,t)\,\mathrm{d}a.$
(18)

In the initial condition of the galaxy, the dust mass of all sizes is set to
be zero. Therefore, in the first time step of the computation, only the dust
grain is produced by the stars. In the next time step, the dust size
distribution produced by the stars evolves through the dust evolution process
in the ISM. In the following, we explain the details of the Asano model.

####  2.2.1 Dust production by AGB

AGB stars are the final phase of the evolution of low and intermediate mass
stars ($<8~{}M_{\odot}$). They have the carbon-oxygen core and are burning
hydrogen and helium surrounding the core. They release heavier elements in the
ISM, and dust grains are formed in the ejecta. The dust size distribution
produced by AGB stars depends on the progenitor mass and suggested that it is
represented by a log-normal distribution with a peak at $\sim
0.1~{}\mathrm{\mu m}$ by Winters et al. (1997). Further, Yasuda & Kozasa
(2012) calculates dust formation by AGB stars with the hydrodynamical
simulation including SiC production. They suggest that the mass distribution
per unit logarithmic bin $a^{4}f(a)$ is described by a log-normal distribution
with a peak at 0.2–0.3 $\mu$m. We assume the size distribution of dust grains
from AGB stars is represented by the log-normal with a peak at 0.1 $\mu$m with
a standard deviation of $\sigma=0.47$. This shape reproduces Figure 7 of
Yasuda & Kozasa (2012). We assume the same size distributions for all dust
species.

As for the dust mass produced by AGB stars, we adopt Ventura et al. (2012);
Ventura et al. (2013). They consider AGB stars with a mass range of 1–8
$M_{\odot}$ and metallicity range of $Z=(0.05$, 0.4) $Z_{\odot}$. The
condensation fraction of the key elements is $\sim 0.3$ for silicate with
progenitor stellar mass $M_{\mathrm{AGB}}=6~{}M_{\odot}$ and initial
metallicity $Z=0.05~{}Z_{\odot}$, and $\sim 0.05$ for carbon with
$M_{\mathrm{AGB}}=3~{}M_{\odot}$ (Ventura et al., 2012). We interpolate and
extrapolate their data to obtain the dust yield at the required stellar mass
and metallicity, as in §2.1.

####  2.2.2 Dust production by SNe II

Massive stars end their lives as supernovae (SNe), and dust grains are formed
in the ejecta of the SNe. The element synthesis determines the dust grains
composition in stars and the mechanism of explosion (Nozawa et al., 2003).
Furthermore, the reverse shock of SN destroys the dust grains by sputtering
(Nozawa et al., 2007; Bianchi & Schneider, 2007). Dust destruction by SN
reverse shock is still under debate and there is no common agreement (e.g.,
Gall et al., 2014; Biscaro & Cherchneff, 2016; Matsuura et al., 2019; Slavin
et al., 2020), but we use Nozawa et al. (2007) as a working hypothesis in this
paper. Nozawa et al. (2007) calculates the grain size distribution produced by
SNe II, and we adopt it for SNe in a progenitor mass range of 13–30
$M_{\odot}$. They calculate only production by SNe II from zero-metallicity
star, but the size distribution and composition of dust produced by the SNe II
are less dependent on the metallicity of the progenitor stars (e.g., Todini &
Ferrara, 2001; Kozasa et al., 2009). Therefore, we assume that the dust
production from the SN does not depend on metallicity and uses the stellar
mass interpolated and extrapolated.

Nozawa et al. (2007) discusses two extreme cases for the structure of the
cores of progenitor stars, mixed and unmixed. According to Hirashita et al.
(2005), the unmixed model gives a better fit to the observed high-$z$
extinction curve of SDSS J1048+4637 at $z=6.2$ (Maiolino et al., 2004). Thus,
we adopt the unmixed model in this paper. The condensation fraction is about
0.003–0.006 (Nozawa et al., 2007).

####  2.2.3 Dust destruction by supernova shock

Dust grains in the ISM are partially destroyed by SN shocks (e.g., Jones et
al., 1996; Nozawa et al., 2006). The SN shocks decrease the total dust mass,
and change the size distribution through the sputtering process (Nozawa et
al., 2006). The sputtering is separated into thermal and non-thermal. The
thermal sputtering is caused by the motion of hot gas, and nonthermal one is
caused by relative motion between gas and dust grain. The sputterings depend
on grain size, gas density, temperature (for thermal sputtering), and the
relative velocity between dust grain and gas (for nonthermal sputtering).

We adopt the result by Yamasawa et al. (2011) for the treatment of the SN
destruction. The grain number density after the destruction by SN shocks,
$f_{\mathrm{X}}^{\prime}(a,t)$, is formulated as

$f_{\mathrm{X}}^{\prime}(a,t)=\int^{a_{\mathrm{max}}}_{a}\eta_{\mathrm{X}}(a,a^{\prime})f_{\mathrm{X}}(a^{\prime},t)\,\mathrm{d}a^{\prime}.$
(19)

Where $\eta_{\mathrm{X}}(a,a^{\prime})$ is the conversion efficiency of SN
sputtering defined as the conversion rate of dust grains from radii
$[a,a+\mathrm{d}a]$ to $[a^{\prime},a^{\prime}+\mathrm{d}a^{\prime}]$.
$a_{\mathrm{max}}$ is the maximum radius of dust grain and we adopt
$a_{\mathrm{max}}=8~{}\mathrm{\mu m}$ (Asano et al., 2013b). This value is
large enough to represent the maximum size produced by shattering and
coagulation (Hirashita & Yan, 2009). Yamasawa et al. (2011) calculate
$\eta_{\mathrm{X}}$ by the method developed by Nozawa et al. (2006). In this
process, the size of dust grains is only reduced by destruction, if
$a>a^{\prime}$, $\eta_{\mathrm{X}}=0$. Equation (19) represents the increasing
amount of dust with radii $[a,a+\mathrm{d}a]$ by SN destruction of dust larger
than $a$. The actual upper limit of the integration corresponds to the maximum
dust size of the distribution before the shock passes. The change of grain
number density caused by SN shock is represented as

$\displaystyle\mathrm{d}f_{\mathrm{X}}(a,t)$
$\displaystyle=f_{\mathrm{X}}^{\prime}(a,t)-[1-\eta_{\mathrm{X}}(a,a)]f_{\mathrm{X}}(a,t)$
$\displaystyle=\int_{0}^{a_{\mathrm{max}}}\eta_{\mathrm{X}}(a,a^{\prime})f_{\mathrm{X}}(a^{\prime},t)\,\mathrm{d}a^{\prime}-f_{\mathrm{X}}(a,t).$
(20)

The change of grain mass density by SN shock at a grain radius $a$ and the
time $t$ is represented by Eq. (20) as,

$\displaystyle\left(\frac{\mathrm{d}{\rho_{\mathrm{X}}(m_{\mathrm{d}},t)}}{\mathrm{d}{t}}\right)_{SN}=$
$\displaystyle
m_{\mathrm{d}}\frac{\mathrm{d}{f_{\mathrm{X}}(a,t)}}{\mathrm{d}{t}}$
$\displaystyle=$
$\displaystyle-\tau_{\mathrm{SN,X}}^{-1}\Big{[}\rho_{\mathrm{X}}(m_{\mathrm{d}},t)$
$\displaystyle-
m_{\mathrm{d}}\int^{a_{\mathrm{max}}}_{0}\eta_{\mathrm{X}}(a,a^{\prime})f_{\mathrm{X}}(a^{\prime},t)\,\mathrm{d}a^{\prime}\Big{]}.$
(21)

If we integrate this equation with respect to $a$ and summing up the dust
species, it agrees with the third term of the right hand side of Equation (4).
The timescale of dust destruction $\tau_{\mathrm{SN}}(t)$ by SN is expressed
as,

$\tau_{\mathrm{SN}}(t)=\frac{M_{\mathrm{ISM}}(t)}{\epsilon
m_{\mathrm{swept}}\gamma_{\mathrm{SN}}(t)},$ (22)

where $\epsilon$ is the efficiency of the dust destruction by SN shocks, and
$\gamma_{\mathrm{SN}}(t)$ is the SN rate. The SN rate is expressed as follows
by (McKee, 1989; Nozawa et al., 2006),

$\gamma_{\mathrm{SN}}(t)=\int^{40~{}M_{\odot}}_{\max(m_{\mathrm{min}}(t),8~{}M_{\odot})}\phi(m)\mathrm{SFR}(t-\tau_{\mathrm{m}})\mathrm{d}m.$
(23)

The integration range is determined by when the SNe can occur (Heger et al.,
2003). When $t<\tau(40~{}M_{\odot})$, $\gamma_{\mathrm{SN}}(t)=0$. We assume
$\epsilon=0.1$ (McKee, 1989; Nozawa et al., 2006).

$m_{\mathrm{swept}}$ is the ISM mass swept by a SN shock. $m_{\mathrm{swept}}$
depends on the density and metallicity of the ISM (Nozawa et al., 2006;
Yamasawa et al., 2011). When the ISM density is high, $m_{\mathrm{swept}}$ is
small because there are many particles that slow down the SN shock. When the
metallicity is high, efficient line cooling with metal results in a faster
shock deceleration and smaller $m_{\mathrm{swept}}$. We use the following
formulae fitted by Yamasawa et al. (2011),

$m_{\mathrm{swept}}=1535n^{-0.202}_{\mathrm{SN}}\left[\left(Z/Z_{\odot}\right)+0.039\right]^{-0.289}~{}M_{\odot},$
(24)

where $n_{\mathrm{SN}}$ is the ISM density surrounding SNe. The fitting
accuracy is within 16% for $0.03~{}\mathrm{cm}^{-3}\leq n_{\mathrm{SN}}\leq
30~{}\mathrm{cm}^{-3}$ and for $10^{-4}\leq Z/Z_{\odot}\leq 1.0$ (Yamasawa et
al., 2011), and we use $n_{\mathrm{SN}}=1.0~{}\mathrm{cm^{-3}}$.

####  2.2.4 Grain growth by metal accretion

In the cold phase of the ISM, the metal in the gas phase is accreted onto the
pre-existing dust grain surface, which increases the radius of the dust grain
and the total mass of the dust, known as grain growth (e.g., Dwek & Scalo,
1980; Draine, 2009; Jones & Nuth, 2011). We assume that grain growth occurs in
the CNM ($T_{\mathrm{gas}}=100$ K and $n_{\mathrm{H}}=30~{}\mathrm{cm^{-3}}$)
and the MC ($T_{\mathrm{gas}}=25$ K and
$n_{\mathrm{H}}=300~{}\mathrm{cm^{-3}}$). The total mass fraction in these ISM
phases is assumed to be 0.5 (Nozawa et al., 2015). We treat only refractory
grains (silicate and carbonaceous dust) and do not consider the icy mantle. We
assume that a grain instantly becomes a sphere with a smooth surface, and we
adopt only the geometric cross-section, i.e., we do not consider the effect of
the Coulomb interaction. Here, we consider only two species of dust grain,
carbonaceous with key element C and silicate with Si. Jones & Nuth (2011)
indicate that in an H2, CO and H2O-rich environment, Si, Fe and Mg accretion
forms silicates through complex chemical reactions such as ice formation.
However, they also say that the spectrum of silicates formed in such a
scenario is inconsistent with actual observations. In this paper, we do not
know the chemical properties well, so it is assumed that only the same key
element X accretes onto the dust of key elements X. Since grain growth
requires a sufficient amount of metals and dust grains in the ISM, grain
growth is difficult to occur in a very young galaxy, but it becomes
efficiently processed at 1 Gyr in general (Asano et al., 2013a).

In the following, we introduce the formalism of size evolution by Hirashita &
Kuo (2011). The collision rate at which an atom of element X with radius $a$
with the surface of the dust grain is expressed as follows (Evans, 1994):

$\mathcal{R}=\pi a^{2}n_{\mathrm{X}}(t)v_{\mathrm{th}},$ (25)

where $n_{\mathrm{X}}(t)$ is the number density of the key element X in gas
phase and $v_{\mathrm{th}}$ is the thermal velocity

$v_{\mathrm{th}}=\left(\frac{8kT_{\mathrm{gas}}}{\pi
m_{\mathrm{X}}}\right)^{1/2},$ (26)

where $k$ is the Boltzmann constant, $T_{\mathrm{gas}}$ is the gas temperature
and $m_{\mathrm{X}}$ is the atomic mass of the key element X. In reality,
since metals other than the corresponding key element may accrete onto the
grain surface, Equation (25) represents the accretion rate associated with the
key element. The evolution of grain mass
$\mathrm{d}m_{\mathrm{d}}(a,t)/\mathrm{d}t$ is

$\frac{\mathrm{d}{m_{\mathrm{d}}(a,t)}}{\mathrm{d}{t}}=g_{\mathrm{X}}^{-1}m_{\mathrm{d}}\alpha_{\mathrm{acc}}\mathcal{R},$
(27)

where $g_{\mathrm{X}}$ is the mass fraction of the key element X in a specific
grain species (silicate: 0.166, graphite: 1.00). We assume Mg1.1Fe0.9SiO4 for
the composition of silicate (Draine & Lee, 1984). $\alpha_{\mathrm{acc}}$ is
the sticking probability of atoms that collide with grains. It is very
difficult to quantify whether the sticking atoms become a part of the grain
(e.g., Jones & Nuth, 2011). This value may be almost 1 in the low temperature
environment (e.g., Zhukovska et al., 2008), so this paper sets
$\alpha_{\mathrm{acc}}=1$ for simplicity.

$n_{\mathrm{X}}(t)$ is estimated as

$n_{\mathrm{X}}(t)=\frac{\rho_{\mathrm{ISM}}^{\mathrm{eff}}}{m_{\mathrm{d}}}\frac{M_{\mathrm{X}}(t)-g_{\mathrm{X}}M_{\mathrm{d,X}}(t)}{M_{\mathrm{ISM}}(t)},$
(28)

where $M_{\mathrm{X}}(t)$ is the total mass of element X (including gas and
dust), $M_{\mathrm{ISM}}(t)$ is the total mass of gas, $M_{\mathrm{d,X}}(t)$
is the dust mass associated with element X, and
$\rho^{\mathrm{eff}}_{\mathrm{ISM}}$ is the effective ISM mass density which
is averaged mass density of the cloud where accretion process occurs.
$\rho^{\mathrm{eff}}_{\mathrm{ISM}}$ is calculated as
$\rho_{\mathrm{ISM}}^{\mathrm{eff}}=\mu m_{\mathrm{H}}n_{\mathrm{H,acc}}$,
where $\mu=1.4$ is the mean atomic weight, $m_{\mathrm{H}}$ is the hydrogen
atom mass, and $n_{\mathrm{H,acc}}$ is the mean hydrogen number density in the
ISM where the accretion process takes place. When
$\eta=\eta_{\mathrm{CNM}}+\eta_{\mathrm{MC}}=0.5$, $n_{\mathrm{H,acc}}$ is 130
$\mathrm{cm^{-3}}$. The second term on the right hand side represents the gas
mass of element X to total ISM mass ratio. From (25)–(28), the grain mass
growth rate is

$\displaystyle\frac{\mathrm{d}{m_{\mathrm{d}}(a,t)}}{\mathrm{d}{t}}=$
$\displaystyle\frac{\pi
a^{2}m_{\mathrm{d}}\alpha_{\mathrm{acc}}v_{\mathrm{th}}}{g_{\mathrm{X}}}n_{\mathrm{X}}(t)$
$\displaystyle=$ $\displaystyle\frac{\pi
a^{2}\alpha_{\mathrm{acc}}\rho_{\mathrm{ISM}}^{\mathrm{eff}}v_{\mathrm{th}}}{g_{\mathrm{X}}}\frac{M_{\mathrm{X}}(t)-g_{\mathrm{X}}M_{\mathrm{d,X}}(t)}{M_{\mathrm{ISM}}(t)}.$
(29)

The total mass growth rate by accretion process
$\left(\mathrm{d}M_{\mathrm{d}}(t)/\mathrm{d}t\right)_{\mathrm{acc}}$ is
calculated by integrating the Equation (29) with respect to the grain radius
and summing up for all dust species in all ISM phases. From Equation (17) and
(29), grain radius growth rate is represented as

$\frac{\mathrm{d}{a}}{\mathrm{d}{t}}=\frac{\alpha_{\mathrm{acc}}\rho_{\mathrm{ISM}}^{\mathrm{eff}}v_{\mathrm{th}}}{4sg_{\mathrm{X}}}\frac{M_{\mathrm{X}}(t)-g_{\mathrm{X}}M_{\mathrm{d,X}}(t)}{M_{\mathrm{ISM}}(t)}.$
(30)

In computation, we solve this equation for each size bin at each time step. By
transferring all the dust that was in the radius bin before metal accretion to
the radius bin after accretion, the evolution of size distribution by metal
accretion can be calculated. The number of dust grains does not change in this
process.

####  2.2.5 Grain-grain collision

We consider two types of gain-grain collisions, shattering, and coagulation.
They only change dust size distribution and conserve total grain mass in the
ISM. Which processes occur is determined by the relative velocity between two
collisional grains. In the case that the relative velocity is fast, shattering
is easy to occur. On the contrary, when the relative velocity is small,
coagulation occurs. Relative velocities between dust grains can be caused by
ubiquitous ISM turbulence (e.g., Draine & Anderson, 1985; Ossenkopf, 1993;
Lazarian & Yan, 2002; McKee & Ostriker, 2007; Hirashita & Yan, 2009; Ormel et
al., 2009). Furthermore, because the dust is thought to be magnetized (Arons &
Max, 1975), it is necessary to consider the motion of grains due to
magnetohydrodynamics (MHD) turbulence. We consider dust collisions by applying
the relative velocity of grains in MHD turbulence calculated by Yan et al.
(2004). The velocity is calculated in consideration of gas drag (hydro drag)
and gyroresonance. When a grain with a relative velocity higher than the
threshold collides, it is shattered into small pieces. Since the larger grains
are affected by turbulence strongly, they have a large relative velocity and
it can occur shattering. Yan et al. (2004) indicate that grains with
$a>10^{-6}$ cm can accelerates to velocity (1–2 $\mathrm{km\,s^{-1}}$) close
to the shattering threshold in CNM. In WNM, gyroresonance accelerates grain
with $a>2$–$3\times 10^{-5}~{}\mathrm{cm}$ to high relative velocity ($\sim
20~{}\mathrm{kms^{-1}}$). On the contrary, small dust has a small relative
velocity and causes coagulation. When the relative velocity is small, the
collision cross section is small, so a high density environment (e.g., MC) is
required for coagulation. We consider silicate and graphite as grain species
while grain-grain collision, and they collide only with the same species.
Furthermore, We treat spherical grains which have a constant density.
Hirashita & Yan (2009) calculates in various ISM phases (including CNM, WNM,
and MC) by shattering and coagulation, and the Asano model applied the same
method.

We consider four types of grain-grain collisions, in other words, the relative
velocity is divided into four types. This treatment is the same as Jones et
al. (1994) and Hirashita & Yan (2009). Considering the shattered and
coagulated grains are radii $a_{1}$ and $a_{2}$, which are called grain 1 and
2, respectively. The mass of grain 1 and 2 is denoted as $m_{1}$ and $m_{2}$.
The relative collisional velocities between grain 1, and 2 are as follows,

  * •

front collision ($v_{1,2}=v_{1}+v_{2}$)

  * •

back-end collision ($v_{1,2}=|v_{1}-v_{2}|$)

  * •

side collision ($v_{1,2}=v_{1}$)

  * •

another side collision ($v_{1,2}=v_{2}$)

where $v_{1}$ and $v_{2}$ are the velocity of grain with radii $a_{1}$ and
$a_{2}$, respectively. We assume that collisions in all directions have the
same probability.

Jones et al. (1996) suggests that the shattering is significantly affecting
the grain size distribution in the ISM. The time evolution of grain mass
density by shattering process is

$\displaystyle\left[\frac{\mathrm{d}{\rho_{\mathrm{X}}(m_{\mathrm{d}},t)}}{\mathrm{d}{t}}\right]_{\mathrm{shat}}$
$\displaystyle=-m_{\mathrm{d}}\rho_{\mathrm{X}}(m_{\mathrm{d}},t)$
$\displaystyle\times\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\alpha\left[m_{\mathrm{d}},m_{1}\right]\rho_{\mathrm{X}}(m_{1},t)~{}\mathrm{d}m_{1}$
$\displaystyle+\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\alpha\left[m_{1},m_{2}\right]m^{1,2}_{\mathrm{shat}}(m_{\mathrm{d}})$
$\displaystyle\times\rho_{\mathrm{X}}(m_{1},t)\rho_{\mathrm{X}}(m_{2},t)~{}\mathrm{d}m_{1}\mathrm{d}m_{2},$
(31)

where $\alpha[m_{1},m_{2}]$ is the collision frequency normalized by two grain
masses and grain number density and expressed as

$\alpha[m_{1},m_{2}]=\begin{cases}0&(v_{1,2}<v_{\mathrm{shat}})\\\
\frac{\sigma_{1,2}v_{1,2}}{m_{1}m_{2}}&(v_{1,2}>v_{\mathrm{shat}})\end{cases},$
(32)

$m_{\mathrm{shat}}^{1,2}(m_{\mathrm{d}})$ represents the total mass of
fragments with masses between $m_{\mathrm{d}}$ and
$m_{\mathrm{d}}+\mathrm{d}m_{\mathrm{d}}$ as the result of collision between
grain 1 and 2. We assume that the distribution of shattered fragments is
proportional to $a^{-3.3}$ (Hellyer, 1970; Jones et al., 1996). $\sigma$ is
the collisional cross-section and represented as

$\sigma_{1,2}=\beta\pi(a_{1}+a_{2})^{2},$ (33)

$\beta$ is the coefficient connecting the cross-section and the geometric
cross-section, assumed $\beta=1$ for simplicity. $v_{\mathrm{shat}}$ is the
threshold of shattering, we assume $1.2~{}\mathrm{km\,s^{-1}}$ and
$2.7~{}\mathrm{km\,s^{-1}}$ for silicate and graphite grains, respectively
(Jones et al., 1996). $a_{\mathrm{min}}$ and $a_{\mathrm{max}}$ are minimum
and maximum radius and we adopt $a_{\mathrm{min}}=0.0003~{}\mathrm{\mu m}$ and
$a_{\mathrm{max}}=8~{}\mathrm{\mu m}$, respectively (Asano et al., 2013b). The
minimum grain radius in the ISM is less well understood, even if
$a_{\mathrm{min}}=0.001~{}\mathrm{\mu m}$, the dust size distribution does not
change significantly (Hirashita, 2012).

The first term on the right hand side of Equation (31) represents the decrease
of grain mass $m_{\mathrm{d}}$ due to destruction by the collisions with other
grains. The second term represents the grain mass $m_{\mathrm{d}}$ increase
due to the fragments resulting from the collision between grain 1 and 2.
Shattering does not produce larger fragments than the original grain, so it
only contributes if either grain 1 or 2 is heavier than $m_{\mathrm{d}}$.

The coagulation occurs when the relative velocity is low. The time evolution
for coagulation is expressed as a similar form to shattering,

$\displaystyle\left[\frac{\mathrm{d}{\rho_{\mathrm{X}}(m_{\mathrm{d}},t)}}{\mathrm{d}{t}}\right]_{\mathrm{coag}}$
$\displaystyle=-m_{\mathrm{d}}\rho_{\mathrm{X}}(m_{\mathrm{d}},t)$
$\displaystyle\times\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\alpha\left[m_{\mathrm{d}},m_{1}\right]\rho_{\mathrm{X}}(m_{1},t)~{}\mathrm{d}m_{1}$
$\displaystyle+\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\alpha[m_{1},m_{2}]m_{\mathrm{coag}}^{1,2}(m_{\mathrm{d}})$
$\displaystyle\times\rho_{\mathrm{X}}(m_{1},t)\rho_{\mathrm{X}}(m_{2},t)~{}\mathrm{d}m_{1}\mathrm{d}m_{2}.$
(34)

and

$\alpha[m_{1},m_{2}]=\begin{cases}\frac{\sigma_{1,2}v_{1,2}}{m_{1}m_{2}}&(v_{1,2}<v_{\mathrm{coag}})\\\
0&(v_{1,2}>v_{\mathrm{coag}})\end{cases},$ (35)

where $m_{\mathrm{coag}}^{1,2}$ is the total mass of coagulated grains:

$m_{\mathrm{coag}}^{1,2}(m_{\mathrm{d}})=\begin{cases}m_{\mathrm{d}}&\mathrm{when}~{}m_{\mathrm{d}}\leq
m_{1}+m_{2}<m_{\mathrm{d}}+\mathrm{d}m_{\mathrm{d}})\\\
0&\mathrm{otherwise}\end{cases}.$ (36)

We use Equation (33) as the collisional cross-section for coagulation.
$v_{\mathrm{coag}}$ is the threshold velocity of coagulation, and grains with
higher relative velocity do not stick. Chokshi et al. (1993) calculate the
threshold velocity as $10^{-3}$–$10^{-1}~{}\mathrm{km\,s^{-1}}$ and it depends
on the grain size. Here we assume that the dust grain is a smooth sphere, but
in the real picture the grain is fluffy (Ossenkopf, 1993). It has been
suggested that the coagulation threshold relative velocity is higher because
fluffiness increases the cross-section of grain collisions (Ormel et al.,
2009; Hirashita & Kobayashi, 2013). In addition, Asano et al. (2014) indicate
that the coagulation threshold suppresses the production of large grain,
producing smaller grain ($a<0.01~{}\mathrm{\mu m}$) than the Mathis et al.
(1977). Therefore, in this paper, it is assumed that coagulation can occur at
all relative velocities without setting the coagulation threshold. The first
term of Equation (34) indicates the decrease of the grain with mass
$m_{\mathrm{d}}$ by coagulation with other grains. The second term indicates
the increase of the grain with mass $m_{\mathrm{d}}$ by coagulation between
grain 1 and 2. Since coagulation works effectively on small size dust, it
becomes effective after shattering becomes effective and small dust increases.
Coagulation shifts the dust size distribution to the larger one.

####  2.2.6 Result of dust evolution model

We show the dust grain (including carbon and silicate) size distribution
calculated with the star formation timescale $\tau_{SF}=3$ Gyr and total
galaxy mass $10^{11}~{}M_{\odot}$ in Figure 1. We note that the total galaxy
mass is merely a normalization for our dust model, and can be rescaled freely.

Figure 1: The time evolution of all species of dust grain size distribution.
Blue, orange, green, and purple curves indicate the age of 100 Myr, 1 Gyr, 5
Gyr, and 13 Gyr, respectively. Black line represents the slope of the MRN
distribution.

In this representation, the ratio of ISM phases are $\eta_{\mathrm{WNM}}=0.5$,
$\eta_{\mathrm{CNM}}=0.3$, $\eta_{\mathrm{MC}}=0.2$. Blue, orange, green, and
purple curves indicate the age of 100 Myr, 1 Gyr, 5 Gyr, and 13 Gyr,
respectively. Black curve indicates the slope of the grain size distribution
suggested by Mathis et al. (1977) which reproduces the MW extinction curve.
This is known as the MRN distribution, expressed by a single power law,

$f(a)\mathrm{d}a\propto a^{-3.5}\mathrm{d}a~{}(0.005~{}\mathrm{\mu
m}<a<0.25~{}\mathrm{\mu m}).$ (37)

The overview of the time evolution of dust size distribution is as follows.

  * •

$<100$ Myr

    * –

Dust production from SNe dominates, and the original size distribution of SN
dust is reflected in the overall dust size distribution.

  * •

100 Myr–1 Gyr

    * –

Metal accretion dominates the evolution of dust.

    * –

Shattering and coagulation become effective, and consequently dust mass
rapidly increases because of the increase of the total surface area per dust
mass.

    * –

Production of PAHs is dominated by shattering.

  * •

1–5 Gyr

    * –

Shattering and coagulation become more effective.

    * –

Metal accretion also becomes more effective thanks to the increased amount of
small dust continuously generated by shattering.

    * –

The increase of dust mass is the most rapid between 1 Gyr and 2 Gyr.

  * •

5–13 Gyr

    * –

The production of dust grain by stars decreases, and shattering and
coagulation dominate the evolution of the grain size distribution.

Details of each step in the evolution of dust size distribution is explained
as follows. At the age of 100 Myr, only a tiny amount of grains exist in the
ISM, and the slope of size distribution is completely different from the slope
of the MRN distribution. Particularly, PAHs are not produced in early
galaxies. Figure 2 shows how different processes (i.e., production by AGB and
SNe, and grain growth in the ISM) contribute to the increase of the total PAH
mass. We note that the contribution of destruction processes including the SN
shock and astration are not shown here. The dust destruction process only
depends on the grain size and species, then it works in the same way for dust
from any production source. Therefore, even if dust reduction is taken into
account, the ratio of dust mass for each source remains the same. In galaxies
younger than 100 Myr, short-lived SNe II (lifetime $\sim 10^{6}$–$10^{7}$ yr)
is the main source of dust supply (e.g., Maiolino et al., 2004; Hiraki &
Hirak, 2008), as the age of the galaxy is too young for stars to evolve into
AGB (lifetime $\sim 10^{8}$–$10^{9}$ yr) (e.g., Morgan & Edmunds, 2003;
Marchenko, 2006). However, smallest grains such as PAHs are supplied from SN
stars, though the amount is very small (e.g., Nozawa et al., 2007).

Figure 2: The evolution of PAH mass per each production source. The total
galaxy mass is $10^{11}~{}M_{\odot}$. Blue, orange, and green curves represent
PAH mass produced by SN, AGB, and evolution in the ISM, respectively.

As the chemical evolution proceeds in the galaxy, the amount of metal in the
ISM increases. Very young galaxies ($\mbox{age}\simeq 20$ Myr) have only a
small supply of dust from SNe. When the galaxy age reaches $\sim 100$ Myr, the
smallest grains (PAH) are gradually formed by shattering, and the mass of the
PAH increases. The galaxy must evolve to reach the critical metallicity for
dust growth to work effectively (Inoue, 2011; Asano et al., 2013a). Since AGB
provides larger size dust grains ($>0.1~{}\mathrm{\mu m}$), their contribution
to PAHs is not significant (Winters et al., 1997; Yasuda & Kozasa, 2012).

At 1 Gyr, the total dust mass continues to increase gradually, while the PAH
mass starts to increase significantly, because the shattering in the ISM
becomes effective. The bump in 10-3–10${}^{-2}~{}\mathrm{\mu m}$ in Figure 1
is the consequence of the activated shattering process. We show the evolution
of the total dust mass of the model galaxy in Figure 3. Solid and dashed lines
represent the dust grain evolution in the ISM (fiducial) and without evolution
(no evolution) case, respectively. Figure 3 clearly demonstrates that, if dust
grains evolve in the ISM, dust mass rapidly increases by metal accretion in
1–2 Gyr.

Figure 3: The evolution of dust to gas mass ratio calculated by the Asano
model.

It has been suggested that in the MW-like galaxy model, when the metallicity
exceeds 0.1 $Z_{\odot}$, the metal accretion process becomes effective and the
dust mass drastically increases (Asano et al., 2013a). When the metal
accretion becomes effective, dust collisions with each other in the ISM become
more likely to occur, and the shattering and coagulation also become
effective. The shattering process results in a significant increase in the
amount of small dust grains including PAH. Since the metal accretion depends
on the total surface area of dust grain (Equation (25)), the shattering
promotes the accretion. Such a dust growth cycle causes the dust mass to
increase nonlinearly. This cycle is effective between 1 Gyr and 2 Gyr in this
model galaxy.

After that, the mass of dust increases and peaks at $\sim 3$ Gyr. This peak
time depends on the timescale of star formation $\tau_{\mathrm{SF}}=3$ Gyr.
After 3 Gyr, the dust mass decreases due to the destruction by SN shocks. The
smaller the dust size, the more effectively the SN destruction works (Nozawa
et al., 2006). In addition, production of dust from stars also decreases due
to the decrease of the SFR. Thus, in total, the dust mass gradually decreases
by SN shock and astration. As the production of dust by stars decreases,
coagulation dominates the evolution of dust size distribution. Due to the SN
shock destruction and coagulation, the dust size distribution is biased toward
larger radius. For galaxies with fully grown dust after 5 Gyr, the dust grain
size distribution finally converges to a similar function to that obtained
from observations of nearby galaxies such as Schurer et al. (2009). A galaxy
with the age in $5\mbox{--}13$ Gyr has a dust distribution with a power-law
slope similar to the MRN.

In contrast, for the no evolution case, the total grain mass does not increase
rapidly and only increases by stellar production with a constant rate up to
$\tau_{\mathrm{SF}}=3$ Gyr. In the age of the $<1$ Gyr galaxy, no evolution
case has a larger grain mass than the mass of the fiducial case. This is
because the no evolution case does not consider the destruction of dust due to
SN shocks. After 3 Gyr, the dust mass decreases by astration. As described
above, if the evolution of dust in the ISM is not taken into account in the
calculation, a rapid increase of dust mass in 1–2 Gyr does not appear.

###  2.3 Stellar SED

We use the version 2 of Pégase (Fioc & Rocca-Volmerange, 1999, Pégase.2) to
produce stellar SEDs. Pégase calculates the stellar emission by stellar
population synthesis (SPS) method with simple stellar populations (SSPs). The
SSP represents the time variation of the SED of a single contemporaneous
stellar population with a single metallicity and abundance pattern. The
monochromatic luminosity per unit wavelength of SSP is expressed as

$L_{\lambda}^{\mathrm{SSP}}(t,Z)=\int^{m_{\mathrm{max}}}_{m_{\mathrm{min}}}L_{\lambda}^{\mathrm{star}}(T_{\mathrm{eff}}(t,m),\log
g(t,m),Z)\phi(m)~{}\mathrm{d}m,$ (38)

where $L^{\mathrm{star}}_{\lambda}$ is the monochromatic luminosity of a star
with the mass in the interval effective temperature $T_{\mathrm{eff}}$,
surface gravity of stellar $g$, metallicity $Z$, and an age of galaxy $t$
(e.g., Conroy, 2013). $m_{\mathrm{max}}$ and $m_{\mathrm{min}}$ are the upper
and lower limit of stellar mass, set to be 100 $M_{\odot}$ and 0.1
$M_{\odot}$, which is the same as the IMF integration range. The effective
temperature $T_{\mathrm{eff}}(t,m)$ and the surface gravity $\log g(t,m)$ are
from the stellar evolutionary track. Pégase.2 uses the evolutionary track
based on the Padova tracks (Bressan et al., 1993; Fagotto et al., 1994b, c, a;
Girardi et al., 1996). The metallicity of the ISM $Z$ evolves with galaxy age
$t$ and it is calculated from Woosley & Weaver (1995) SN II models. Since only
the evolutionary track table with metallicities
$Z=(0.005,0.02,0.2,0.4,1,0,2.5,5.0)~{}Z_{\odot}$ is prepared, they use the
interpolated value. Pégase assumes that a star releases metal into ISM only at
the end of its life, and the recycling model is not instantaneous. The library
of stellar spectra used by Pégase.2 is divided into two according to effective
temperature $T_{\mathrm{eff}}$. For $T_{\mathrm{eff}}<50000~{}\mathrm{K}$, the
library comes from Lejeune et al. (1997, 1998, corrected version (BaSeL-2.0)).

A monochromatic luminosity from total stars at time $t$ is calculated by
weighting $L_{\lambda}^{\mathrm{SSP}}$ at galaxy age $t^{\prime}$ with star
formation rate SFR,

$L_{\lambda}(t)=\int^{t^{\prime}=t}_{t^{\prime}=0}\int^{Z=Z_{\mathrm{max}}(t-t^{\prime})}_{Z=0}\mathrm{SFR}(t-t^{\prime})L_{\lambda}^{\mathrm{SSP}}(t^{\prime},Z(t-t^{\prime}))~{}\mathrm{d}Z\mathrm{d}t^{\prime},$
(39)

where $Z_{\mathrm{max}}(t-t^{\prime})$ is the maximum metallicity at time
$t-t^{\prime}$. In order to take into account the stars that were born at the
time of the galaxy’s birth to the stars that are just born, the value is
integrated over time. The time lag $t-t^{\prime}$ represents the time
difference between the formation of a star and the end of the evolution of the
star. We chose the Schmidt law (Schmidt, 1959, Equation (5)) for the SFR.

###  2.4 Dust properties

Radiative transfer is the method to calculate the propagation of energy in
systems of various sizes (from isolated gas clouds to galaxies). In the galaxy
ISM, radiation is mainly affected by absorption and scattering by dust grains.
One of the easiest ways to calculate radiative transfer is to assume that the
ISM has a homogeneous distribution. However, it has been observed that actual
galaxies have a more complex structure in general (e.g., Field et al., 1969;
McKee & Ostriker, 1977). If a homogeneous dust distribution is assumed, the
optical depth of the dust is larger than that in the case of inhomogeneous
distribution. In other words, if the dust mass is estimated with a homogeneous
distribution, the attenuation per dust grain is overestimated, and the dust
mass would be underestimated as a consequence. Therefore, in this paper, we
consider a clumpy dust distribution.

The calculation of three-dimensional radiative transfer with clumpy dust
usually requires substantial computational cost. Neufeld (1991) and Hobson &
Padman (1993) introduced the method that solves the radiative transfer with
MGA in a one-dimensional plain parallel galaxy (Varosi & Dwek, 1999; Inoue,
2005). The MGA treats the dusty region as a kiloparsec-size huge grain called
mega-grain, and regards absorption and scattering behave in the same way as
typical grains with effective optical properties. We approximate the complex
distribution of stars, dust grains, and gas in the model galaxy, to simplify
costly calculations in three-dimensional space. In this approximation, the
distribution of young stars is clumpy and the young stars are embedded by
mega-grain. In contrast, old stars are supposed to distribute smoothly in a
diffuse way. The light emitted by young stars is stronger attenuated than the
light emitted by older stars due to the surrounding mega-grains. Inoue (2005)
researched the effect of changing criterion of young star $t_{\mathrm{y}}$. He
conclude 10 Myr is best fit to the MW attenuation, and we apply it in this
paper. Assuming thermal and chemical equilibrium with temperature $T<10^{4}$
K, the ISM is represented by two phases, WNM and CNM (e.g., Wolfire et al.,
2003; Koyama & Inutsuka, 2002). The relation of the thermal pressure with
hydrogen density is expressed by fitting the phase diagram (Inoue, 2005) as,

$\displaystyle\frac{p/k}{10^{4}~{}\mathrm{K\,cm^{-3}}}$
$\displaystyle=\frac{n_{\mathrm{H,WNM}}}{1~{}\mathrm{cm^{-3}}},~{}$
$\displaystyle(\mathrm{WNM})$ (40)
$\displaystyle\frac{p/k}{10^{4.5}~{}\mathrm{K\,cm^{-3}}}$
$\displaystyle=\left(\frac{n_{\mathrm{H,CNM}}}{10^{3}~{}\mathrm{cm^{-3}}}\right)^{0.7},~{}$
$\displaystyle(\mathrm{CNM})$ (41)

where $p$ is the pressure. $n_{\mathrm{H,WNM}}$ and $n_{\mathrm{H,CNM}}$ are
the hydrogen density of WNM and CNM, respectively. We regard the WNM as a
homogeneous interclump medium and the CNM as a clump. The clump radius
$r_{\mathrm{cl}}$ is calculated by assuming it to be self-gravitating (Inoue,
2005),

$r_{\mathrm{cl}}=\frac{1}{\rho_{\mathrm{cl}}}\sqrt{\frac{15p}{4\pi
G}}=\frac{1}{\mu m_{\mathrm{p}}n_{\mathrm{H,CNM}}}\sqrt{\frac{15p}{4\pi
G}}\sim 10.4~{}\mathrm{pc},$ (42)

where $\rho_{\mathrm{cl}}=\mu m_{\mathrm{p}}n_{\mathrm{H,CNM}}$ is the clump
density, $\mu$ is the mean atomic density, $G$ is the gravitational constant
and $m_{\mathrm{p}}$ is the proton mass. Clumps exist in the interclump
medium. We assume all clumps are spherical and have a constant radius and
density. In this approximation, we use the mass absorption and scattering
coefficient and the scattering asymmetry parameter of dust grain,
$k_{\mathrm{abs}}$, $k_{\mathrm{scat}}$ and $g_{\mathrm{d}}$ respectively,
averaged by dust size distribution calculated by the Asano model:

$\displaystyle k_{\mathrm{abs}}$
$\displaystyle=\frac{\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\pi
a^{2}Q_{\mathrm{abs}}(a)f(a)\mathrm{d}a}{\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}m_{\mathrm{d}}(a)f(a)\mathrm{d}a},$
(43) $\displaystyle k_{\mathrm{scat}}$
$\displaystyle=\frac{\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\pi
a^{2}Q_{\mathrm{scat}}(a)f(a)\mathrm{d}a}{\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}m_{\mathrm{d}}(a)f(a)\mathrm{d}a},$
(44) $\displaystyle g_{\mathrm{d}}$
$\displaystyle=\frac{\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}g(a)\pi
a^{2}Q_{\mathrm{scat}}(a)f(a)\mathrm{d}a}{\int^{a_{\mathrm{max}}}_{a_{\mathrm{min}}}\pi
a^{2}Q_{\mathrm{scat}}(a)f(a)\mathrm{d}a},$ (45)

where $f(a)$ is the dust number distribution, $Q_{\mathrm{abs}}(a)$ and
$Q_{\mathrm{scat}}(a)$ are the absorption and scattering coefficient, and
$g(a)$ is the scattering asymmetry parameter of a grain, respectively. In this
model $Q_{\mathrm{abs}}(a)$, $Q_{\mathrm{scat}}(a)$, and $g(a)$ are calculated
by the Mie theory (Bohren & Huffman, 1983). We adopt Draine & Lee (1984) and
Laor & Draine (1993) for silicate and graphite, and Li & Draine (2001) for PAH
as optical parameters. The mass extinction coefficient and scattering albedo
are defined as

$k_{\mathrm{d}}=k_{\mathrm{abs}}+k_{\mathrm{scat}}.$ (46)

In the MGA, we replace the optical properties, namely the extinction
coefficient per unit length $\kappa$, the scattering albedo $\omega$, and the
scattering asymmetry parameter $g$ with the effective ones. The relative
optical depth of a clump with the interclump medium is

$\tau_{\mathrm{cl}}=(\rho_{\mathrm{cl}}-\rho_{\mathrm{icm}})k_{\mathrm{d}}Dr_{\mathrm{cl}}.$
(47)

where $D$ is the dust-to-gas mass ratio calculated by the Asano model. The
extinction coefficient per unit length of the medium by clump is

$\kappa_{\mathrm{mg}}=n_{\mathrm{cl}}\pi
r^{2}_{\mathrm{cl}}P_{\mathrm{int}}(\tau_{\mathrm{cl}})=\frac{3f_{\mathrm{cl}}}{4r_{\mathrm{cl}}}P_{\mathrm{int}}(\tau_{\mathrm{cl}}),$
(48)

where $n_{\mathrm{cl}}$ is the number density of clump and $f_{\mathrm{cl}}$
is the clump filling fraction,

$f_{\mathrm{cl}}=\frac{n_{\mathrm{H}}-n_{\mathrm{H,WNM}}}{n_{\mathrm{H,CNM}}-n_{\mathrm{H,WNM}}}.$
(49)

We assume that mean hydrogen number density in the galaxy $n_{\mathrm{H}}$ has
a constant value 1 cm3. $P_{\mathrm{int}}(\tau)$ is the interaction
probability against parallel light by a sphere with optical depth $\tau$, and
represented as

$P_{\mathrm{int}}(\tau)=1-\frac{1}{2\tau^{2}}+\left(\frac{1}{\tau}+\frac{1}{2\tau^{2}}\right)e^{-2\tau}.$
(50)

This equation is obtained by integrating the light incident on the sphere in
all directions and taking the ratio to the case where the optical depth of the
sphere is zero (see appendix C of Varosi & Dwek (1999) for details). The
extinction coefficient of interclump medium is

$\kappa_{\mathrm{icm}}=k_{\mathrm{d}}D\rho_{\mathrm{icm}}.$ (51)

Thus, the effective extinction coefficient is expressed as

$\kappa_{\mathrm{eff}}=\kappa_{\mathrm{mg}}+\kappa_{\mathrm{icm}}.$ (52)

The scattering albedo of clump is

$\omega_{\mathrm{cl}}=\omega_{\mathrm{d}}P_{\mathrm{esc}}(\tau_{\mathrm{cl}},\omega_{\mathrm{d}}),$
(53)

where $\omega_{\mathrm{d}}=k_{\mathrm{scat}}/k_{\mathrm{d}}$ is the scattering
albedo of normal grain averaged by grain size distribution and

$P_{\mathrm{esc}}(\tau,\omega)=\frac{\frac{3}{4\tau}P_{\mathrm{int}}(\tau)}{1-\omega\left[1-\frac{3}{4\tau}P_{\mathrm{int}}(\tau)\right]},$
(54)

is the photon escape probability from a sphere grain. The effective scattering
albedo is

$\omega_{\mathrm{eff}}=\frac{\omega_{\mathrm{cl}}\kappa_{\mathrm{mg}}+\omega_{\mathrm{d}}\kappa_{\mathrm{icm}}}{\kappa_{\mathrm{eff}}}.$
(55)

The light entering the clump is scattered by the dust in the clump and escapes
in various directions. An optical parameter that indicates in which direction
light escapes is called an asymmetry parameter of clump and it is defined as
$g_{\mathrm{cl}}=\langle\cos\theta_{\mathrm{esc}}\rangle$, where
$\theta_{\mathrm{esc}}$ is the angle between enter and escape directions.
$g_{\mathrm{cl}}$ is given by fitting the Monte Carlo calculation result in
Varosi & Dwek (1999), and represented by the following empirical formula,

$g_{\mathrm{cl}}(\tau_{\mathrm{cl}},\omega_{\mathrm{cl}},g_{\mathrm{d}})=g_{\mathrm{d}}-C\left(1-\frac{1+e^{-B/A}}{1+e^{(\tau_{\mathrm{cl}}-B)/A}}\right),$
(56)

where

$\displaystyle A$ $\displaystyle\equiv
1.5+4g_{\mathrm{d}}^{3}+2\omega_{\mathrm{d}}\sqrt{g_{\mathrm{d}}}\exp(-5g_{\mathrm{d}}),$
(57) $\displaystyle B$ $\displaystyle\equiv
2-g_{\mathrm{d}}(1-g_{\mathrm{d}})-2\omega_{\mathrm{d}}g_{\mathrm{d}},$ (58)
$\displaystyle C$
$\displaystyle\equiv\frac{1}{3-\sqrt{2g_{\mathrm{d}}}-2\omega_{\mathrm{d}}g_{\mathrm{d}}(1-g_{\mathrm{d}})}.$
(59)

The effective asymmetry parameter is

$g_{\mathrm{eff}}=\frac{g_{\mathrm{cl}}\kappa_{\mathrm{mg}}+g_{\mathrm{d}}\kappa_{\mathrm{icm}}}{\kappa_{\mathrm{eff}}}.$
(60)

###  2.5 Radiative transfer in a one-dimensional galaxy

Figure 4: The geometry of a one-dimensional plain parallel galaxy.

We assume a one-dimensional plane-parallel galaxy along the $z$-axis for
solving radiative transfer shown in Figure 4. We set two kinds of disks in the
model. One is a gas+dust disk containing young stars (disk 1). Stars are born
in cold, dense regions (CNM), so we assume young stars are surrounded by
clamps, and disk 1 is also full of interclump medium. Disk 1 has a constant
density of stars, gases, and grains. Optical depth $\tau$ is defined with
constant effective extinction coefficient $\kappa_{\mathrm{eff}}$ as

$\mathrm{d}\tau=-\kappa_{\mathrm{eff}}\,\mathrm{d}z,$ (61)

with $\tau=0$ at $z=h_{\mathrm{d}}$ and
$\tau=\kappa_{\mathrm{eff}}h_{\mathrm{d}}$ at $z=0$. Where $2h_{\mathrm{d}}$
is the thickness of disk 1. The other disk contains only exponentially and
smoothly distributed old stars (disk 2). The thickness of the disk 2 is
$4h_{\mathrm{d}}$, which is twice larger than that of disk 1. These two disks
are stacked so that their centers are aligned. In the above condition,
radiative transfer is formulated as,

$\mu\frac{\mathrm{d}{I(\tau,\mu)}}{\mathrm{d}{\tau}}=-I(\tau,\mu)+S(\tau,\mu),$
(62)

where $I(\tau,\mu)$ is the specific intensity at $\tau$ and $\mu\equiv
1/\cos\theta$. The $\theta$ is the angle between the ray and $z$-axis. Source
function $S$ is represented as

$S(\tau,\mu)=\frac{\eta_{\ast}(\tau)}{\kappa_{\mathrm{eff}}}+\omega_{\mathrm{eff}}\int^{1}_{-1}I(\tau,\mu^{\prime})\Phi(g_{\mathrm{eff}},\mu,\mu^{\prime})\mathrm{d}\mu^{\prime},$
(63)

where $\eta_{\ast}$ is the stellar emissivity and $\Phi$ is the scattering
phase function. Here we adopt the Henyey-Greenstein phase function (Henyey &
Greenstein, 1941)

The first term in the right hand side of Equation (63) represents the
intensity of light emitted from the star that escapes from the clump where the
star was born. The second term in the right hand side represents the integral
of the light scattered in the considering direction $\mu$ among the scattered
light by dust grain. The boundary conditions in this galaxy at $z=0$ and
$z=h_{\mathrm{d}}$ are

$\displaystyle I(\tau=\kappa_{\mathrm{eff}}h_{\mathrm{d}},\mu)$
$\displaystyle=I(\tau=\kappa_{\mathrm{eff}}h_{\mathrm{d}},-\mu),$ (64)
$\displaystyle I(\tau=0,\mu<0)$
$\displaystyle=-\frac{\int^{\infty}_{h_{\mathrm{d}}}\eta_{\ast}(z)\,\mathrm{d}z}{\mu}.$
(65)

The stellar emissivity is normalized by

$\int^{\infty}_{-\infty}\eta_{\ast}(z)\mathrm{d}z=1.$ (66)

The intrinsic emissivity from young stars in disk 1
$\eta^{\mathrm{young}}_{\mathrm{\ast}}$is $1/2h_{\mathrm{d}}$ for $|z|\leq
h_{\mathrm{d}}$ and zero for $|z|>h_{\mathrm{d}}$ because it is normalized by
Equation (66). The energy emitted by young star is absorbed by the dust in
clump surrounding the star. Since the escape probability from the clump is
Equation (54), the emissivity from young stars is represented as

$\eta^{\mathrm{young}}_{\mathrm{\ast}}\begin{cases}P_{\mathrm{est}}(\tau_{\mathrm{cl}},\omega_{\mathrm{cl}})/2h_{\mathrm{d}}&(|z|\leq
h_{\mathrm{d}})\\\ 0&(|z|>h_{\mathrm{d}})\end{cases}.$ (67)

The old stars in disk 2 are distributed with exponential diffusion along the
$z$-axis. From the normalization Equation
$(\ref{equ:emissivity_normalization})$, the emissivity from the old stars is

$\eta_{\mathrm{\ast}}^{\mathrm{old}}(z)=\frac{e^{-|z|/2h_{\mathrm{d}}}}{4h_{\mathrm{d}}}.$
(68)

The total stellar emissivity at $z$ is represented as

$\eta_{\ast}(z)=f_{\mathrm{y}}(t)\eta_{\mathrm{\ast}}^{\mathrm{young}}(z)+(1-f_{\mathrm{y}}(t))\eta_{\mathrm{\ast}}^{\mathrm{old}}(z),$
(69)

where $f_{\mathrm{y}}(t)$ is the luminosity fraction emitted by young stellar
at age $t$, it is calculated by

$f_{\mathrm{y}}(t)=\frac{\int^{\mathrm{min}[t_{\mathrm{y}},t]}_{0}\int^{\mathrm{Z}_{\mathrm{max}}(t-t^{\prime})}_{0}\mathrm{SFR}(t-t^{\prime})L_{\lambda}^{\mathrm{SSP}}(t^{\prime},\mathrm{Z}(t-t^{\prime}))\,\mathrm{dZ}\,\mathrm{d}t^{\prime}}{L_{\lambda}(t)}.$
(70)

In the radiative transfer calculation, we calculate Equation (62) and (63)
iteratively until the ratio to the previous loop of the source function in all
directions on the galaxy surface converged to 10-10.

###  2.6 Dust temperature distribution

The UV and optical photons emitted by stars heat dust grains. The heated dust
grains release the energy by emission from MIR to FIR wavelength photons.
Large grains have an equilibrium temperature determined by the stellar
radiation field. In contrast, very small grains cannot establish radiative
equilibrium, and it does not have a stable equilibrium temperature (Draine &
Anderson, 1985; Draine & Li, 2001; Li & Draine, 2001; Takeuchi et al., 2003,
2005; Horn et al., 2007). Since they have small heat capacities, they are
easily heated by photons and then rapidly cooled Thus, the (instantaneous)
temperature of very small grains is inevitably stochastic, we calculate
temperature distribution by Monte Carlo simulation.

####  2.6.1 Stochastic heating

The rate at which a dust grain absorbs photons in the energy range
$[E,E+\mathrm{d}E]$ and time interval $[t,t+\mathrm{d}t]$ is expressed as

$\mathrm{d}p(a,\lambda)=\pi
a^{2}Q_{\mathrm{abs}}(a,\lambda)\bar{u}_{\lambda}\frac{\lambda^{3}}{h_{\mathrm{p}}^{2}c}\mathrm{d}E\mathrm{d}t,$
(71)

where $\bar{u}_{\lambda}$ is the mean energy density per wavelength in a
galaxy, $h_{\mathrm{p}}$ is the Plank constant, and $c$ is the speed of light,
respectively (e.g., Draine & Anderson, 1985; Takeuchi et al., 2003, 2005). The
energy density actually varies depending strongly on its spatial position in a
galaxy. Therefore, considering the different energy densities for each
position, the calculation time becomes enormous. Thus, we use the mean energy
density $\bar{u}_{\lambda}$ which is calculated in the same way as Fioc &
Rocca-Volmerange (2019):

$L_{\lambda}^{0}-L_{\mathrm{obs}}=c\bar{u}_{\lambda}k_{\mathrm{abs}}M_{\mathrm{d}},$
(72)

where $L_{\lambda}^{0}$ and $L_{\mathrm{obs}}$ are the intrinsic stellar
luminosity and observed luminosity calculated by transfer radiation.
Therefore, the mean energy density is represented as

$\bar{u}_{\lambda}=\frac{L_{\lambda}^{0}-L_{\mathrm{obs}}}{ck_{\mathrm{abs}}M_{\mathrm{d}}}.$
(73)

Equation (71) can be regarded as the probability density distribution if the
time interval $\mathrm{d}t$ is appropriately small. For each dust size,
$\mathrm{d}t$ is determined so that the maximum collision probability among
all wavelengths is 0.01,

$\mathrm{d}t(a)=0.01\left[\pi
a^{2}Q_{\mathrm{abs}}(a,\lambda)\bar{u}_{\lambda}\frac{\lambda^{3}}{h_{\mathrm{p}}^{2}c}\mathrm{d}E\right]^{-1}.$
(74)

For simplicity, we assume that the energy of an absorbed photon is totally
used to heat the dust grains, represented as

$E(T+\Delta T)=E(T)+\frac{h_{\mathrm{p}}c}{\lambda},$ (75)

where $E(T)$ is the enthalpy of dust grains at temperature $T$ and $\Delta T$
is the increment of temperature. We adopt the Debye model for calculating the
enthalpy of dust grains (Li & Draine, 2001). The enthalpy of silicate and
graphite grains are

$\displaystyle E_{\mathrm{sil}}(T)$
$\displaystyle=(N_{\mathrm{atom}}-2)k\left[2f_{2}\left(\frac{T}{500~{}\mathrm{K}}\right)+f_{3}\left(\frac{T}{1500~{}\mathrm{K}}\right)\right],$
(76) $\displaystyle E_{\mathrm{gra}}(T)$
$\displaystyle=(N_{\mathrm{C}}-2)k\left[f_{2}\left(\frac{T}{863~{}\mathrm{K}}\right)+2f_{2}\left(\frac{T}{2504~{}\mathrm{K}}\right)\right],$
(77)

where

$f_{n}(a)\equiv n\int^{1}_{0}\frac{y^{n}\mathrm{d}y}{\exp(y/x)-1}.$ (78)

The subscripts ’sil’ and ’gra’ represent the silicate and graphite grains,
respectively. Equation (78) is the $n$ dimensional Debye function.
$N_{\mathrm{atom}}$ and $N_{\mathrm{C}}$ are the number of atoms in a grain,
they are expressed as,

$N=\dfrac{\frac{4}{3}\pi a^{3}\rho N_{\mathrm{A}}}{M},$ (79)

where $\rho$ is the mass density, $M$ is the mass number and $N_{\mathrm{A}}$
is the Avogadro constant. For carbonaceous (graphite or PAH) grain,
$\rho=2.26~{}\mathrm{g/cm^{3}}$ and $M=12.0~{}\mathrm{g/mol}$(Draine & Lee,
1984). For silicate grain, $\rho=3.50~{}\mathrm{g/cm^{3}}$ and
$M=172.25~{}\mathrm{g/mol}$(Li & Draine, 2001). In polycyclic aromatic
hydrocarbon (PAH) grains, we consider C-C bond modes are same as graphite and
C-H bond modes component is added to Equation (77),

$E_{\mathrm{pah}}(T)=E_{\mathrm{gra}}+\frac{\mathrm{H}}{\mathrm{C}}N_{\mathrm{C}}\sum^{3}_{j=1}\left(\frac{h_{\mathrm{p}}\nu_{j}}{\exp(h_{\mathrm{p}}\nu_{j}/kT)-1}\right).$
(80)

The index $j$ represents the C-H out-of-plane bending modes
($\nu_{1}/c=886~{}\mathrm{cm^{-1}}$), in-plane bending modes
($\nu_{2}/c=1161~{}\mathrm{cm^{-1}}$), and stretching modes
($\nu_{3}/c=3030~{}\mathrm{cm^{-1}}$), respectively (Draine & Li, 2001).
$\frac{\mathrm{H}}{\mathrm{C}}$ is the hydrogen to carbon ratio. We adopt the
following empirical formula (Li & Draine, 2001),

$\frac{\mathrm{H}}{\mathrm{C}}=\begin{cases}0.5&(N_{\mathrm{C}}<25)\\\
\frac{0.5}{\sqrt{N_{\mathrm{C}}/25}}&(25<N_{\mathrm{C}}<100)\\\
0.25&(N_{\mathrm{C}}>100)\end{cases}.$ (81)

####  2.6.2 Dust cooling

The equation of emission of dust grains with radius $a$ is formulated as,

$4\pi\epsilon(T,a)=4\pi\left(\pi a^{2}\right)\int
Q_{\mathrm{abs}}(\lambda)\frac{2h_{\mathrm{p}}c^{2}}{\lambda^{5}}\frac{\mathrm{d}\lambda}{\exp\left(\frac{h_{\mathrm{p}}c}{\lambda
kT}\right)-1},$ (82)

where $T$ is the temperature of dust grain and $\epsilon(T,a)$ is the emission
power per unit time per unit solid angle. In the case that dust grains do not
absorb energy while cooling, since emission energy and changes in internal
energy are balanced, the following equation holds,

$\frac{\mathrm{d}E(T,a)}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}t}=-4\pi\epsilon(T,a).$
(83)

This equation can not be solved analytically, but we can calculate, but we can
get the temperature variation by numerical calculation.

####  2.6.3 Result of dust temperature distribution

Figure 5 is the result of dust temperature distribution of silicate.

Figure 5: Temperature distribution of the several grain sizes of silicate
calculated by Monte Carlo calculation. The blue, orange, green, red, and
purple curves indicate $3.98\times 10^{-8}$, $1.26\times 10^{-7}$, $3.98\times
10^{-7}$, $1.26\times 10^{-6}$, and $3.98\times 10^{-5}$ cm grains,
respectively.

The galaxy is the face on ($\mu=1$) MW-like galaxy model at the age of 13 Gyr.
The condition of the galaxy is the same as §2.2.6. The radius of galaxy
$R_{\mathrm{gal}}$ is 10 kpc and scale height of dust $h_{\mathrm{d}}$ is 150
pc which is the typical scale height of cold dust in the MW (e.g., Binney &
Merrifield, 1998).

The temperature of small grains is very widely distributed from 1 to 4,000 K.
When grain size becomes larger, the temperature range becomes narrower and
approaches the equilibrium temperature. The equilibrium temperature is
represented by Draine & Lee (1984); Takeuchi et al. (2003) as

$T_{\mathrm{eq}}\simeq\left(\frac{h_{\mathrm{p}}c}{\pi
k}\right)\left[\frac{945u}{960\pi(2\pi Aa)h_{\mathrm{p}}c}\right]^{1/6},$ (84)

and

$u\equiv\int^{\infty}_{0}u_{\lambda}\mathrm{d}\lambda,$ (85)

where $h_{\mathrm{p}}$ is the Planck constant, and we adopt
$A_{\mathrm{sil}}=1.34\times 10^{-3}$ cm for silicate grains (Drapatz &
Michel, 1977) and $A_{\mathrm{C}}=3.20\times 10^{-3}$ cm for carbonaceous
grains (Draine & Lee, 1984). When grain size is $3.98\times 10^{-5}$ cm, the
equilibrium temperature of that grain is about 19 K by Equation (84) and it is
equal to the result calculated by the Monte Carlo calculation.

Figure 6 is the temperature distribution of graphite. Graphite grains have
broad temperature compared with silicate grains in Figure 5. The difference
comes from differences of internal energy between silicate and graphite.

Figure 6: Temperature distribution of the several grain sizes of graphite
grain. Calculation parameters and color coordinates are the same as Figure 5.

We show the temperature distribution of PAH in Figure 7. Comparing PAH
temperature distribution with graphite, PAHs stay in a narrower temperature
range because PAHs have an additional term in the equation of internal energy
(Equation 80). Almost the same behavior is seen in the result of Draine & Li
(2007). From the result of temperature distribution, some dust grains might
exceed the sublimation temperature, which is 1500 K or higher (e.g., Baskin &
Laor, 2018). If we assume that the grain above 1500 K has sublimated and
calculated its effect by removing it from the galaxy, the mass of the
sublimated grains is only a few percent of the total dust mass. Thus, the
effect of the sublimation for the result is negligible, and we do not consider
the effect of sublimation temperature in our model for simplicity.

Figure 7: Temperature distribution of the several grain sizes of PAHs.
Calculation parameters are the same as Figure 5. The blue, orange, green, and
red curves indicate $3.98\times 10^{-8}$, $1.26\times 10^{-7}$, $3.98\times
10^{-7}$, $1.00\times 10^{-6}$ cm grains, respectively.

###  2.7 Dust radiation

The dust radiation depends on the temperature distribution
$\frac{\mathrm{d}{P_{i}(a)}}{\mathrm{d}{T}}$ calculated by the method of the
above sections. The monochromatic luminosity of a dust grain of species $i$
(silicate, graphite, neutral PAH, or ionized PAH) is expressed as

$L_{i}^{\mathrm{grain}}(a,\lambda)=4\pi a^{2}\pi\int
Q_{\mathrm{abs}}^{i}(\lambda)B_{\lambda}(T)\frac{\mathrm{d}{P_{i}(a)}}{\mathrm{d}{T}}~{}\mathrm{d}T,$
(86)

where $B_{\lambda}$ is the blackbody radiation and $Q_{\mathrm{abs}}^{i}$ is
the absorption coefficient of dust species $i$. Total luminosity at wavelength
$\lambda$ is represented as,

$L(\lambda)=\sum_{i}Q_{\mathrm{abs}}^{i}(\lambda)\int
L_{i}^{\mathrm{grain}}(a,\lambda)f_{i}(a)~{}\mathrm{d}a.$ (87)

$f_{i}(a)$ is the dust number distribution of dust species $i$ from the dust
evolution model.

##  3 Result: Milky Way-like galaxy model SED

In Figure 8, we show the model result with a face-on ($\mu=1$) MW-like galaxy
model (§2.6.3 and §2.2.6) at the age of 13 Gyr (the same setting as §2.6.3).

Figure 8: The result of our SED model with MW-like galaxy parameters at an age
of 13 Gyr. The black curve represents the overall emission of the galaxy.
Other color curves express each dust species (blue: ionized PAH, orange:
neutral PAH, green: silicate, red: graphite).

Each curve in Figure 8 represents the corresponding emission species. At the
912 $\mathrm{\mathring{A}}$ wavelength, we see the cutoff of the Lyman break.
The UV to IR wavelength region is dominated by stellar emission. Numerous PAH
lines are prominent in the mid-IR, and the far-IR range is dominated by the
continuum emission from large graphite grains. The emission of silicate is
weaker than that of graphite in the wavelength range of $200~{}\mathrm{\mu m}$
or less, and is effective only in the longer wavelength range. The temperature
of graphite and silicate fitted by a gray body are 28 K and 26 K,
respectively. The difference in emission and temperature between graphite and
silicate is caused by the number of grains.

Figure 9 is the time evolution of our galaxy SED model with the MW-like galaxy
model parameters.

Figure 9: The evolution of the SED of a MW-like galaxy. Parameters are the
same as Figure 8. Blue, orange, green, red, and purple curves indicate age of
100 Myr, 1, 5, 10, 13 Gyr, respectively.

The purple curve represents the SED of a MW-like galaxy at the age of 13 Gyr,
as Figure 8. Blue, orange, green, and red curves represent the age of 100 Myr,
1, 5, and 10 Gyr, respectively. Figure 9 shows that the UV region emitted by
stars monotonically decreases with the evolution. Since we assume a closed
box, the gas mass decreases monotonically as it is consumed by star formation.
The SFR is proportional to the gas mass, hence SFR also decreases
monotonically, and the UV radiation also decreases. The overview of the time
evolution of the SED is as follows.

  * •

<100 Myr

    * –

Stellar emission dominates the SED and PAH emission does not exist yet.

  * •

100 Myr–1 Gyr

    * –

Stellar emission still dominates the SED, but the dust emission including PAH
gradually becomes prominent.

  * •

1–5 Gyr

    * –

Dust emission dominates the SED and dust emission becomes strongest at this
age.

  * •

5–13 Gyr

    * –

The emission both from stars and dust gradually decreases, along with the
decline of the star formation rate.

The details of the SED evolution are explained below.

At 100 Myr, since only a very small amount of dust has been produced, stellar
radiation is not attenuated, and dust radiation is weak. Particularly, PAHs
are not produced in young galaxies, hence the MIR radiation is very faint. For
an SED model that supposes a constant size distribution without considering
the evolution of the dust size distribution (e.g., Schurer et al., 2009), many
PAHs are observed even in such young galaxies, and different conclusions are
deduced.

The evolution of metallicity, dust mass and bolometric luminosity for each
component are shown in Figure 10.

Figure 10: The evolution of metallicity, dust mass, and bolometric luminosity
of each component of the galaxy. Parameters are set to be the same as Figure
8. Metallicity, dust mass, and bolometric luminosities are normalized by solar
metallicity $Z_{\odot}=0.02$ (Anders & Grevesse, 1989), maximum value of it,
and overall bolometric luminosity, respectively. The calculation was performed
with the age of the logarithmic scale bin.

The dust mass is normalized with respect to its maximum value. Dust mass and
luminosity are tightly correlated. Here we adopt $Z_{\odot}=0.02$ (Anders &
Grevesse, 1989). At 1 Gyr, the dust mass is gradually increasing, and along
with this, the IR radiation from dust becomes prominent. Because the PAH mass
increases by the evolution in the ISM, their characteristic mid-infrared line
emission can be seen.

The dust emission becomes strongest at 3 Gyr if we adopt the star formation
timescale $\tau_{\mathrm{SF}}=3$ Gyr. As predicted, in the MW-like model, when
the metallicity exceeds 0.1 $Z_{\odot}$, the metal accretion process becomes
effective and the dust mass increases (Asano et al., 2013a). Star formation is
active, but the UV continuum from young stars is strongly attenuated due to
the increase of dust mass.

After 3 Gyr, the dust mass decreases due to the destruction by SN shocks and
astration. Dust radiation also decreases with the age of the galaxy. This is
not only due to the reduction of dust mass, but also due to the decline of the
UV light from young stars to heat dust grains.

Focusing on the metallicity, it reaches 1.6 $Z_{\odot}$ at 13 Gyr, which is
larger than the solar metallicity. This is due to the assumption of the closed
box model. In the closed box model, there is no inflow of gas and outflow of
ISM and the metallicity increases monotonically. However, considering the
infall model, the metallicity is reduced because the ISM is diluted by the
inflow of gas (Erb et al., 2006).

We should note the difference in the evolution of each species of dust grains
in Figure 10. The increase of graphite emission is more gradual than the
increase of other components. The metallicity of a galaxy exceeds the critical
metallicity in accretion onto the dust surface, leading to the sharp rise of
the dust mass and emission (Asano et al., 2013a). Shattering produces small
dust grains and makes the surface area of dust larger, and consequently leads
to the boost of the accretion efficiency. However, since we regard almost all
small carbonaceous grains as PAHs, the mass of graphite grains does not have
the discontinuous increase. The bolometric luminosity of dust emission is
dominated by graphite in all epochs.

##  4 Discussion

###  4.1 Effect of star formation timescale

The effect of star formation (SF) timescale $\tau_{\mathrm{SF}}$ is shown in
Figure 11.

Figure 11: The effect of star formation timescale $\tau_{\mathrm{SF}}$ for
galaxy SED. The geometrical model parameters are the same as Figure 11. The
galaxy age is increasing from left to right, SF timescale increases from top
to bottom. The age and SF timescale are written on each plot. Black, orange,
green, red thick curves represent overall, graphite, silicate, and PAHs
luminosity, respectively. Blue thin curve is an intrinsic (unattenuated)
stellar emission.

In Figure 11, the SEDs are calculated with different SF timescale and age from
Figure 8, while geometrical parameters are kept the same. The galaxy age
increase from left to right the panels (100 Myr, 1, and 10 Gyr), and SF
timescale increase from top to bottom the panels ($\tau_{\mathrm{SF}}=$ 0.5,
1, 5 Gyr). Three characteristic trends are observed in Figure 11. First, the
UV light emitted by stars and the IR light emitted by graphite and silicate
grains at the age of 100 Myr decrease with increasing $\tau_{\mathrm{SF}}$.
This is because the age of the galaxy is sufficiently small with respect to
the SFR, and the larger the time scale of the star formation rate, the fewer
stars form. In these young galaxies, PAHs are not produced and the PAHs
emission is not observed in the model with any $\tau_{\mathrm{SF}}$ yet. This
indicates that the dust mass in early galaxies is dominated by production from
stars instead of accretion processes.

Second, the overall bolometric luminosities tend to be stronger when the age
of the galaxy is equal to the SF timescale. In this evolutionary phase, the
SFR is still large and a large amount of dust exists in the galaxy.

Lastly, when the age of the galaxy is older than the SF timescale, the galaxy
has very weak stellar emission due to consumption of most of the gas in the
ISM which is an ingredient of star formation. The dust emission is also very
weak in the galaxy because both, decreasing dust mass and UV light which is
the source of heating dust grains.

###  4.2 Effect of geometrical parameters

Figure 12 shows the effect of changing the dust scale height of galaxy
$h_{\mathrm{d}}$ for our galaxy SED model at an age of 13 Gyr.

Figure 12: Galaxy SEDs with various dust scale heights $h_{\mathrm{d}}$ at the
age of 13 Gyr. The star formation history is the same as Figure 8. Blue,
orange, and green curves represent 75, 150 (fiducial), and 300 pc,
respectively.

The SFH is the same as Figure 8. Blue, red, and green curves represent
$h_{\mathrm{d}}=$ 75, 150 (fiducial), and 300 pc, respectively. Intrinsic
stellar radiation does not depend on $h_{\mathrm{d}}$. Since the optical depth
of the galaxy is defined as $\tau=\kappa_{\mathrm{eff}}h_{\mathrm{d}}$, $\tau$
increases as $h_{\mathrm{d}}$ increases. Then, the absorption by the dust
grain becomes stronger, and the observed UV radiation becomes weaker. Since
the energy absorbed by the dust grain increases, the radiation in the IR
region becomes stronger.

Since our SED model assumes an axisymmetric one-dimensional disk, the galaxy
has no structure in the radial direction and does not determine the radius.
However, in reality, when the radius of the galaxy changes and the volume
changes, the density of the dust clump changes and the optical depth also
changes. In our model, the optical depth depends on clump filling fraction
(Equation (49)), and we assume that $n_{\mathrm{H}}$ is constant. Therefore,
if we consider a galaxy with a volume in which $n_{\mathrm{H}}$ changes
significantly from 1 cm-3, it is necessary to consider it, but it is not
implemented in our model and is a future work.

###  4.3 Effect of the ISM phase fraction

In the current model, we consider three phases in the ISM: WNM, CNM and MC.
Figure 13 and 14 are the effects of the ISM phase fraction for dust size
distribution. The parameters except the ISM fraction are the same as the MW-
like galaxy. The value of the fraction of cold region $\eta$ is changed while
keeping the ratio in the cold region constant to
$\eta_{\mathrm{CNM}}:\eta_{\mathrm{MC}}=3:2$ in this section.

Figure 13: The dust size distribution with cold ISM region fraction $\eta=0.5$
(fiducial, solid), 1 (dashed), and 0 (dot-dashed). Blue, orange, and purple
curves are the age of 100 Myr, 1, and 13 Gyr galaxies, respectively. Note, the
100 Myr galaxy has three overlapping curves.

Solid, dashed, and dot-dashed curves represent fiducial, $\eta=1$
($\eta_{\mathrm{WNM}}=0.0$, $\eta_{\mathrm{CNM}}=0.6$, and
$\eta_{\mathrm{MC}}=0.4$), and $\eta=0$ ($\eta_{\mathrm{WNM}}=1.0$,
$\eta_{\mathrm{CNM}}=0.0$, and $\eta_{\mathrm{MC}}=0.0$) case, respectively.

At 100 Myr galaxy, there is no difference in three cases, since the stellar
production dominates dust size distribution and is not affected by the ISM
fractions.

At 1 Gyr galaxy, the size distribution of small $\eta$ cases have small amount
of grain in radius of $>2\times 10^{-1}~{}\mathrm{\mu m}$ region because
shattering is more likely to occur thanks to collisions between larger grains.
The bump in 10-3–10-2 $\mathrm{\mu m}$ is generated by accretion on the
$>10^{-3}$ size of the grain surface. The bump is not observed in the $\eta=0$
case, because in this case the accretion process on grains in the cold regions
is not included. The $\eta=1$ result has a larger bump than the fiducial case.
This is because shattering is not effective yet, and the larger the fraction
of cold regions is, the more effective the metal accretion. On the contrary,
large amounts of intermediate size grains ($2\times 10^{-2}$–$2\times
10^{-1}~{}\mathrm{\mu m}$) in a small $\eta$ case. This results from the
coagulation process in WNM.

At the 13 Gyr galaxy, a small amount of grain is observed in the $\eta=0$
case. When grain evolution occurs in only WNM, a strong shattering process
generates the large amount of small grains. Small dust grains are largely
destroyed by SN shock (Nozawa et al., 2006), hence the mass of dust grains
effectively decreases. The large bump in $10^{-1}~{}\mathrm{\mu m}$ is caused
by the balance between strong shattering and coagulation. At the grain sizes
of $<1~{}\mathrm{\mu m}$, the size distribution of the $\eta=1$ case has a
smaller amount of dust. This is because the shattering is weak in the $\eta=1$
case and the metal accretion does not occur as effectively as the fiducial
case, because the WNM is not considered in the calculation of the dust
evolution. Further, maximum grain radius in the $\eta=1$ case reaches
$>1~{}\mathrm{\mu m}$, as weak shattering efficiency in CNM and MC.

The effect of $\eta$ for total dust mass of the MW-like galaxy model is shown
in Figure 14.

Figure 14: The evolution of total dust mass with $\eta=0.5$ (fiducial, orange
solid), 1 (blue dashed), and 0 (green dot-dashed).

Very small amount of total dust mass is observed in the $\eta=0$ case because
the case does not consider the mass-increasing process other than the
production from the stars. Around 1 Gyr, the $\eta=1$ case has a larger amount
of total dust mass than the fiducial case. This is because shattering is still
less effective in this age, and metal accretion dominates the increase of
total dust mass. After 1 Gyr, the increase of total dust mass in the $\eta=1$
case is slower than that in the fiducial case, since the rapid increase cycle
is less effective in the $\eta=1$ case. The total mass at 13 Gyr is determined
by the balance between destruction by SN shocks and the dust growth by metal
accretion.

The galaxy SED at 13 Gyr with the cold ISM region fraction $\eta=0.5$
(fiducial, solid orange), 1 (dashed blue), and 0 (dot-dashed green).

Figure 15: The galaxy SED at 13 Gyr with $\eta=0.5$ (fiducial, orange solid),
1 (blue dashed), and 0 (green dot-dashed). The parameters are the same as the
MW-like galaxy model except $\eta$.

In $\eta=0$ case, the SED has weaker dust attenuation and dust radiation than
that of the fiducial case (Figure 14), since dust mass in all radius is
smaller than the fiducial case. Though the difference of total dust mass
between the fiducial and the $\eta=1$ case is small, the size distribution of
the two cases has a large difference. The $\eta=1$ case has a lot of large
dust grains and a few small dust grains. The attenuation in the
$0.1~{}\mathrm{\mu m}$ wavelength is mainly dominated by the
$a<1~{}\mathrm{\mu m}$ radius of grain. As the large grain has large heat
capacity, the radiation from the large grain is weaker than that from the
small grain. Therefore, $\eta=1$ case has weaker attenuation in the UV region
and weaker radiation in the IR region than the fiducial case.

###  4.4 Effect of the coagulation threshold

Coagulation can occur when the relative velocity of grains $v_{\mathrm{coag}}$
is slower than the threshold velocity. However, we do not adopt the threshold
of coagulation to reproduce the dust size distribution of the MW (Asano et
al., 2013a; Nozawa et al., 2015). If $v_{\mathrm{coag}}$ is adapted to the
dust model, since the small radius grains have lower relative velocity, the
small grains are more likely to coagulate. Conversely, the large radius grains
have large relative velocity and the grains cannot occur the coagulation.
Hirashita & Yan (2009) show that the grain radius increases only up to
0.01–0.1 $\mathrm{\mu m}$, because the radii $a>0.1~{}\mathrm{\mu m}$ grains
have the relative velocity larger than $v_{\mathrm{coag}}$. Therefore, a lower
coagulation threshold velocity suppresses the effect of coagulation. First, We
show the effect of $v_{\mathrm{coag}}$ for total dust grain mass in Figure 16.
The dashed and solid lines represent the fiducial case (no
$v_{\mathrm{coag}}$) and the adopted $v_{\mathrm{coag}}$ case (we call it a
suppressed coagulation case). Galaxy parameters are the same as §2.2.6 except
$v_{\mathrm{coag}}$. We calculate the coagulation velocity threshold in the
same formula as Hirashita & Yan (2009). $v_{\mathrm{coag}}$ between grain 1
and 2 is represented as

$v_{\mathrm{coag}}=21.4\left[\frac{a^{3}_{1}+a^{3}_{2}}{(a_{1}+a_{2})^{3}}\right]^{1/2}\frac{\gamma^{5/6}}{E^{1/3}R_{1,2}^{5/6}s^{1/2}},$
(88)

where suffix 1 and 2 represents the each value of grain 1 and 2,
$R_{1,2}\equiv a_{1}a_{2}/(a_{1}+a_{2})$ is the reduced radius of the grains,
$\gamma$ is the surface energy per unit area, and $E$ is related to the
Poisson ratios ($\nu_{1}$ and $\nu_{2}$) and the Young modulus ($E_{1}$ and
$E_{2}$) by $1/E\equiv(1-\nu_{1})^{2}/E_{1}+(1-\nu_{2})^{2}/E_{2}$. The value
of $\gamma$, $\nu$, and $E$ are 25 $\mathrm{erg/cm^{2}}$, 0.17
$\mathrm{erg/cm^{2}}$ and $5.4\times 10^{11}~{}\mathrm{dyn/cm^{2}}$ for
silicate, and 12 $\mathrm{erg/cm^{2}}$, 0.5 and $3.4\times
10^{10}~{}\mathrm{dyn/cm^{2}}$ for graphite from Chokshi et al. (1993).

Figure 16: The effect of coagulation threshold velocity for total dust mass.
Solid and dashed curves represent suppressed coagulation and fiducial case,
respectively. The galaxy properties are the same as the MW-like galaxy model
except $v_{\mathrm{coag}}$.

From Figure 16, the effect of $v_{\mathrm{coag}}$ for total dust mass is very
small. Coagulation itself decrease total surface area of grains and suppress
the cross section of the metal accretion. On the other hand, as the grain size
increases, shattering is more likely to occur, and smaller radius grain
increases. Since these effects are balanced, coagulation only suppresses the
increase in total dust mass and has a small effect. Coagulation is difficult
to occur in young galaxies (age > 1 Gyr), and becomes effective after the
shattering process becomes effective. Therefore, when the coagulation becomes
effective, the rapid increase in the total dust grain mass has already
finished, and the coagulation does not significantly affect the total mass,
but only changes the size distribution of the dust grain.

Second, we show the effect of $v_{\mathrm{coag}}$ for dust size distribution
in Figure 17.

Figure 17: The effect of $v_{\mathrm{coag}}$ for dust size distribution. Blue,
orange, green, and purple curves represent the dust grain size distribution
with the age of 100 Myr, 1 Gyr, 5 Gyr, and 13 Gyr, respectively. Solid and
dashed lines represent suppressed coagulation and fiducial case, respectively.

In the age of the 100 Myr and 1 Gyr galaxies, since effective dust evolution
in the ISM has not started yet, there is no difference in the dust size
distribution between the two cases. On the other hand, $v_{\mathrm{coag}}$
strongly affects the dust distribution after 1 Gyr galaxy. $v_{\mathrm{coag}}$
suppresses the coagulation between larger grains and determines the maximum
radius of dust grain. Coagulation shifts the size distribution to larger
sizes, thus, in the suppressed coagulation case, the dust size distribution is
biased toward the smaller one, and the slope is also different from MRN.
Therefore, adopting the low $v_{\mathrm{coag}}$ in the MW-like galaxy model
leads to the dust size distribution to be different from the MRN distribution,
and thus we do not adopt $v_{\mathrm{coag}}$ in our model.

###  4.5 Comparison between closed-box and infall model

The comparison between the result of the closed-box and infall model is shown
in Figure 18. We adopt the following equation for the infall rate (Inoue,
2011):

$\frac{\mathrm{d}{M_{\mathrm{infall}}}}{\mathrm{d}{t}}=\frac{M_{\mathrm{infall}}}{\tau_{\mathrm{infall}}}\exp\left(-\frac{t}{\tau_{\mathrm{infall}}}\right),$
(89)

where $\tau_{\mathrm{infall}}$ is the timescale of infall, and
$M_{\mathrm{infall}}$ is the total mass that flows into the galaxy by infall
as $t\rightarrow\infty$. For the infall model, the initial mass of a galaxy is
set to zero, and primordial gas (zero-metallicity) fall onto the galaxy with
$M_{\mathrm{infall}}=10^{11}~{}M_{\odot}$ and $\tau_{\mathrm{infall}}=6$ Gyr.

Figure 18: The comparison of MW-like model galaxy SED at age
$t_{\mathrm{gal}}=13$ Gyr with infall and closed-box model. The orange curve
represents closed-box model (same as Figure 8), the blue curve represents
infall model with infall time scale $\tau_{\mathrm{infall}}=6$ Gyr.

The time evolution of the SFR and dust mass are plotted in Figure 19. The star
formation history is very different between the two models.

Figure 19: The time evolution of dust mass and SFR of the galaxy with closed-
box and infall model (same galaxies as Figure 18). Colors represent the
difference of quantities: the ratio of $\mathrm{SFR}(t_{\mathrm{gal}})$ and
maximum value of it (red), and the ratio of dust mass $M_{\mathrm{dust}}$ and
maximum value of it (blue). The solid and dashed curves represent closed-box
and infall models, respectively.

While the SFR of the closed-box model case monotonically decrease, the SFR of
the infall model gradually increase and it reaches to a peak at
$t_{\mathrm{gal}}=5$ Gyr (close to infall timescale $t_{\mathrm{infall}}=6$
Gyr), after that the SFR decrease gradually. Since the SFR at 13 Gyr of the
infall model is higher than that of the closed-box model, the ratio of younger
stars is increased in the infall model, and then the luminosity of the UV
region becomes stronger. On the other hand, the continuum at near IR
wavelengths emitted from old stars becomes slightly weaker due to the smaller
amount of old stars. The metallicity is 1.6 $Z_{\odot}$ in the closed box
model, while it is 0.86 $Z_{\odot}$ in the infall model which is closer to the
solar metallicity.

In the infall model case, the peak of the dust mass comes later because of the
different star formation history. It leads to an increase of the IR emission
which is emitted by dust grains. In general, the infall model tends to delay
the evolution of the galaxy.

###  4.6 Radio emission

Our model does not include the radio emission. Because, in a normal galaxy,
the luminosity of the radio region is only $<10^{-4}$ of the overall
bolometric luminosity of a galaxy (Condon, 1992). Above $\sim 1$ mm, radio
emission is swamped in dust emission for normal galaxies. The radio is mainly
emitted by synchrotron radiation from relativistic electrons accelerated in
supernova remnants and free-free emission by Hii region, which is a ionized by
the radiation of heavy and young stars (Klein et al., 1988; Carlstrom &
Kronberg, 1991). Since both radio sources are associated with SN explosions,
their radiation is considered to depend on the SN rate (SNR) (Condon, 1992).
In particular, many galaxies with strong synchrotron radiation by a jet from
active galactic nuclei have been observed (e.g., Carilli et al., 1991; Laing &
Bridle, 2002), and we will take it into account in our future work.

##  5 Conclusions

In this paper, we construct a new galaxy SED model including the dust
evolution in galaxies consistent with the chemical evolution (Asano et al.,
2013a, b, 2014). The dust model considers several evolutionary processes of
the dust production by AGB stars and SNe II, the destruction by SN shocks in
the ISM, the grain growth by metal accretion to grain surface, and the two
types of grain-grain collision, the shattering and coagulation. The stellar
radiation is calculated by PÉGASE.2 (Fioc & Rocca-Volmerange, 1999). Based on
this, we constructed a radiative transfer model with a one-dimensional plain
parallel geometry equipped with the mega-grain approximation for fast
computation (Varosi & Dwek, 1999; Inoue, 2005). For the radiation from dust,
we take into account the stochastic heating of dust grains by Monte Carlo
simulation. As a fiducial model, we assumed the Schmidt law with star
formation time scale $\tau_{\mathrm{SF}}=3~{}\mathrm{Gyr}$, the Salpeter IMF
(Salpeter, 1955), and the closed box model. The ISM phase fractions were set
as $\eta_{\mathrm{WNM}}=0.5$, $\eta_{\mathrm{CNM}}=0.3$, and
$\eta_{\mathrm{MC}}=0.2$, scale height of dust is $h_{\mathrm{d}}=150$ pc, and
the threshold of coagulation velocity is removed. Our model indicates that
early galaxies ($\sim 100$ Myr) produce a small amount of dust. The PAHs,
which dominate the MIR wavelength region, have not been produced yet, in
particular. The SED at the age of 100 Myr is dominated by stellar emission.
Then the amount of dust mass and emission explosively increases at the age of
about 3 Gyrs. Subsequently, the dust mass and emission from both the stars and
dust decreases, along with the decline of the star formation rate. Since this
model treats the evolution of dust appropriately, we can apply it to any age
of a galaxy as far as the model assumptions are valid.

## Acknowledgements

First of all, we offer our sincere thanks to the anonymous referee for her/his
enormous effort to read through the article and invaluably important comments
and suggestions that improved the quality of the paper very much. We are
grateful to the colleagues in the Lab for fruitful discussions and comments.
We thank H. Kobayashi and A.K. Inoue for helpful comments on the coding of
dust evolution model. This work has been supported by JSPS Grants-in-Aid for
Scientific Research (17H01110, 19H05076, and 21H01128). This work has also
been supported in part by the Sumitomo Foundation Fiscal 2018 Grant for Basic
Science Research Projects (180923), and the Collaboration Funding of the
Institute of Statistical Mathematics “New Development of the Studies on Galaxy
Evolution with a Method of Data Science”.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Alongi et al. (1993) Alongi M., Bertelli G., Bressan A., Chiosi C., Fagotto F., Greggio L., Nasi E., 1993, A&AS, 97, 851
  * Anders & Grevesse (1989) Anders E., Grevesse N., 1989, Geochim. Cosmochim. Acta, 53, 197
  * Arendt et al. (2010) Arendt R. G., et al., 2010, ApJ, 725, 585
  * Arons & Max (1975) Arons J., Max C. E., 1975, Astrophys. J., 196, L77
  * Asano et al. (2013a) Asano R. S., Takeuchi T. T., Hirashita H., Inoue A. K., 2013a, Earth, Planets Sp., 65, 213
  * Asano et al. (2013b) Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2013b, MNRAS, 432, 637
  * Asano et al. (2014) Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2014, Mon. Not. R. Astron. Soc., 440, 134
  * Baskin & Laor (2018) Baskin A., Laor A., 2018, Mon. Not. R. Astron. Soc., 474, 1970
  * Bertelli et al. (1994) Bertelli G., Bressan A., Chiosi C., Fagotto F., Nasi E., 1994, A&AS, 106, 275
  * Bianchi & Schneider (2007) Bianchi S., Schneider R., 2007, Mon. Not. R. Astron. Soc., 378, 973
  * Binney & Merrifield (1998) Binney J., Merrifield M., 1998, Galactic Astronomy
  * Biscaro & Cherchneff (2016) Biscaro C., Cherchneff I., 2016, Astron. Astrophys., 589, 1
  * Bohren & Huffman (1983) Bohren C. F., Huffman D. R., 1983, Absorption and scattering of light by small particles. New York: Wiley, https://ui.adsabs.harvard.edu/abs/1983asls.book.....B
  * Borkowski et al. (2006) Borkowski K. J., et al., 2006, ApJ, 642, L141
  * Bressan et al. (1993) Bressan A., Fagotto F., Bertelli G., Chiosi C., 1993, A&AS, 100, 647
  * Burgarella et al. (2020) Burgarella D., Nanni A., Hirashita H., Theulé P., Inoue A. K., Takeuchi T. T., 2020, Astron. Astrophys., 637, A32
  * Calura et al. (2008) Calura F., Pipino A., Matteucci F., 2008, A&A, 479, 669
  * Carilli et al. (1991) Carilli C. L., Perley R. A., Dreher J. W., Leahy J. P., 1991, ApJ, 383, 554
  * Carlstrom & Kronberg (1991) Carlstrom J. E., Kronberg P. P., 1991, ApJ, 366, 422
  * Cazaux et al. (2005) Cazaux S. M., Caselli P., Walmsley M., Tielens A. G., 2005, Proc. Int. Astron. Union, 1, 325
  * Ceccarelli et al. (2018) Ceccarelli C., Viti S., Balucani N., Taquet V., 2018, Mon. Not. R. Astron. Soc., 476, 1371
  * Chokshi et al. (1993) Chokshi A., Tielens A. G. G. M., Hollenbach D., 1993, ApJ, 407, 806
  * Condon (1992) Condon J. J., 1992, ARA&A, 30, 575
  * Conroy (2013) Conroy C., 2013, ARA&A, 51, 393
  * De Looze et al. (2020) De Looze I., et al., 2020, Mon. Not. R. Astron. Soc., 496, 3668
  * De Vis et al. (2017) De Vis P., et al., 2017, Mon. Not. R. Astron. Soc., 471, 1743
  * De Vis et al. (2019) De Vis P., et al., 2019, Astron. Astrophys., 623, A5
  * Draine (2009) Draine B. T., 2009, Space Sci. Rev., 143, 333
  * Draine & Anderson (1985) Draine B. T., Anderson N., 1985, ApJ, 292, 494
  * Draine & Lee (1984) Draine B. T., Lee H. M., 1984, ApJ, 285, 89
  * Draine & Li (2001) Draine B. T., Li A., 2001, ApJ, 551, 807
  * Draine & Li (2007) Draine B. T., Li A., 2007, ApJ, 657, 810
  * Drapatz & Michel (1977) Drapatz S., Michel K., 1977, A&A, 56, 353
  * Dwek (1998) Dwek E., 1998, ApJ, 501, 643
  * Dwek & Scalo (1980) Dwek E., Scalo J. M., 1980, ApJ, 239, 193
  * Dwek et al. (2007) Dwek E., Galliano F., Jones A. P., 2007, ApJ, 662, 927
  * Erb et al. (2006) Erb D. K., Shapley A. E., Pettini M., Steidel C. C., Reddy N. A., Adelberger K. L., 2006, Astrophys. J., 644, 813
  * Evans (1994) Evans A., 1994, The dusty universe. Wiley, Chichester, https://ui.adsabs.harvard.edu/abs/1994duun.book.....E
  * Fagotto et al. (1994a) Fagotto F., Bressan A., Bertelli G., Chiosi C., 1994a, A&AS, 104, 365
  * Fagotto et al. (1994b) Fagotto F., Bressan A., Bertelli G., Chiosi C., 1994b, A&AS, 105, 29
  * Fagotto et al. (1994c) Fagotto F., Bressan A., Bertelli G., Chiosi C., 1994c, A&AS, 105, 39
  * Ferrara et al. (2016) Ferrara A., Viti S., Ceccarelli C., 2016, Mon. Not. R. Astron. Soc. Lett., 463, L112
  * Field et al. (1969) Field G. B., Goldsmith D. W., Habing H. J., 1969, ApJ, 155, L149
  * Fioc & Rocca-Volmerange (1999) Fioc M., Rocca-Volmerange B., 1999, arXiv e-prints
  * Fioc & Rocca-Volmerange (2019) Fioc M., Rocca-Volmerange B., 2019, A&A, 623, A143
  * Gall et al. (2011a) Gall C., Andersen A. C., Hjorth J., 2011a, Astron. Astrophys., 528, A13
  * Gall et al. (2011b) Gall C., Andersen A. C., Hjorth J., 2011b, Astron. Astrophys., 528, A14
  * Gall et al. (2014) Gall C., et al., 2014, Nature, 511, 326
  * Ginolfi et al. (2018) Ginolfi M., Graziani L., Schneider R., Marassi S., Valiante R., Dell’Agli F., Ventura P., Hunt L. K., 2018, Mon. Not. R. Astron. Soc., 473, 4538
  * Girardi et al. (1996) Girardi L., Bressan A., Chiosi C., Bertelli G., Nasi E., 1996, A&AS, 117, 113
  * Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003, ApJ, 591, 288
  * Hellyer (1970) Hellyer B., 1970, Obs., 90, 55
  * Henyey & Greenstein (1941) Henyey L. C., Greenstein J. L., 1941, Nature, 147, 613
  * Hiraki & Hirak (2008) Hiraki A., Hirak H., 2008, Rev. Mex. Física, 54, 44
  * Hirashita (2012) Hirashita H., 2012, MNRAS, 422, 1263
  * Hirashita & Aoyama (2019) Hirashita H., Aoyama S., 2019, MNRAS, 482, 2555
  * Hirashita & Ferrara (2002) Hirashita H., Ferrara A., 2002, Mon. Not. R. Astron. Soc., 337, 921
  * Hirashita & Kobayashi (2013) Hirashita H., Kobayashi H., 2013, Earth, Planets Sp., 65, 1083
  * Hirashita & Kuo (2011) Hirashita H., Kuo T. M., 2011, MNRAS, 416, 1340
  * Hirashita & Yan (2009) Hirashita H., Yan H., 2009, MNRAS, 394, 1061
  * Hirashita et al. (2005) Hirashita H., Nozawa T., Kozasa T., Ishii T. T., Takeuchi T. T., 2005, Mon. Not. R. Astron. Soc., 357, 1077
  * Hobson & Padman (1993) Hobson M. P., Padman R., 1993, MNRAS, 264, 161
  * Hollenbach & McKee (1979) Hollenbach D., McKee C. F., 1979, ApJS, 41, 555
  * Horn et al. (2007) Horn K., Perets H. B., Biham O., 2007, arXiv e-prints
  * Inoue (2003) Inoue A. K., 2003, Publ. Astron. Soc. Japan, 55, 901
  * Inoue (2005) Inoue A. K., 2005, Mon. Not. R. Astron. Soc., 359, 171
  * Inoue (2011) Inoue A. K., 2011, Earth, Planets Sp., 63, 1027
  * Jones & Nuth (2011) Jones A. P., Nuth J. A., 2011, A&A, 530, 1
  * Jones et al. (1994) Jones A. P., Tielens A. G. G. M., Hollenbach D. J., McKee C. F., 1994, ApJ, 433, 797
  * Jones et al. (1996) Jones A. P., Tielens A. G. G. M., Hollenbach D. J., 1996, ApJ, 469, 740
  * Kalvans (2017) Kalvans J., 2017, Proc. Int. Astron. Union, 13, 374
  * Klein et al. (1988) Klein U., Wielebinski R., Morsi H. W., 1988, A&A, 190, 41
  * Kobayashi et al. (2006) Kobayashi C., Umeda H., Nomoto K., Tominaga N., Ohkubo T., 2006, ApJ, 653, 1145
  * Koyama & Inutsuka (2002) Koyama H., Inutsuka S., 2002, in 8th Asian-Pacific Reg. Meet. Vol. II. pp 159–160
  * Kozasa et al. (2009) Kozasa T., Nozawa T., Tominaga N., Umeda H., Maeda K., Nomoto K., 2009, in Cosm. Dust - Near Far. p. 43
  * Kuo & Hirashita (2012) Kuo T. M., Hirashita H., 2012, Mon. Not. R. Astron. Soc. Lett., 424, 34
  * Laing & Bridle (2002) Laing R. A., Bridle A. H., 2002, MNRAS, 336, 1161
  * Laor & Draine (1993) Laor A., Draine B. T., 1993, ApJ, 402, 441
  * Laporte et al. (2017) Laporte N., et al., 2017, ApJ, 837, L21
  * Lazarian & Yan (2002) Lazarian A., Yan H., 2002, Astrophys. J., 566, L105
  * Lejeune et al. (1997) Lejeune T., Cuisinier F., Buser R., 1997, A&AS, 125, 229
  * Lejeune et al. (1998) Lejeune T., Cuisinier F., Buser R., 1998, Astron. Astrophys. Suppl. Ser., 130, 65
  * Lesniewska & Michałowski (2019) Lesniewska A., Michałowski M. J., 2019, Astron. Astrophys., 624, 4
  * Li & Draine (2001) Li A., Draine B. T., 2001, Astrophys. J., 554, 778
  * Lisenfeld & Ferrara (1998) Lisenfeld U., Ferrara A., 1998, ApJ, 496, 145
  * Liu & Hirashita (2019) Liu H.-M., Hirashita H., 2019, Mon. Not. R. Astron. Soc., 490, 540
  * Maiolino et al. (2004) Maiolino R., Schneider R., Oliva E., Bianchi S., Ferrara A., Mannucci F., Pedani M., Roca Sogorb M., 2004, Nature, 431, 533
  * Mancini et al. (2015) Mancini M., Schneider R., Graziani L., Valiante R., Dayal P., Maio U., Ciardi B., Hunt L. K., 2015, Mon. Not. R. Astron. Soc., 451, L70
  * Marchenko (2006) Marchenko S. V., 2006, in Stellar Evol. Low Met. Mass Loss, Explos. Cosmol.. p. 299, https://ui.adsabs.harvard.edu/abs/2006ASPC..353..299M
  * Mathis et al. (1977) Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425
  * Matsuura et al. (2019) Matsuura M., et al., 2019, Mon. Not. R. Astron. Soc., 482, 1715
  * McKee (1989) McKee C. F., 1989, ApJ, 345, 782
  * McKee & Ostriker (1977) McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, Annu. Rev. Astron. Astrophys., 45, 565
  * Michałowski (2015) Michałowski M. J., 2015, Astron. Astrophys., 577, 1
  * Michałowski et al. (2010) Michałowski M., Watson D., Hjorth J., 2010, ApJ, 712, 942
  * Morgan & Edmunds (2003) Morgan H. L., Edmunds M. G., 2003, MNRAS, 343, 427
  * Nanni et al. (2020) Nanni A., Burgarella D., Theulé P., Côté B., Hirashita H., 2020, Astron. Astrophys., 641, A168
  * Neufeld (1991) Neufeld D. A., 1991, ApJ, 370, L85
  * Nozawa et al. (2003) Nozawa T., Kozasa T., Umeda H., Maeda K., Nomoto K., 2003, Astrophys. J., 598, 785
  * Nozawa et al. (2006) Nozawa T., Kozasa T., Habe A., 2006, Astrophys. J., 648, 435
  * Nozawa et al. (2007) Nozawa T., Kozasa T., Habe A., Dwek E., Umeda H., Tominaga N., Maeda K., Nomoto K., 2007, Astrophys. J., 666, 955
  * Nozawa et al. (2015) Nozawa T., Asano R. S., Hirashita H., Takeuchi T. T., 2015, Mon. Not. R. Astron. Soc. Lett., 447, L16
  * Ormel et al. (2009) Ormel C. W., Paszun D., Dominik C., Tielens A. G. G. M., 2009, Astron. Astrophys., 502, 845
  * Ossenkopf (1993) Ossenkopf V., 1993, Astron. Astrophys., 280, 617
  * Pipino et al. (2011) Pipino A., Fan X. L., Matteucci F., Calura F., Silva L., Granato G., Maiolino R., 2011, A&A, 525, A61
  * Raiteri et al. (1996) Raiteri C. M., Villata M., Navarro J. F., 1996, A&A, 315, 105
  * Rouillé et al. (2020) Rouillé G., Jäger C., Henning T., 2020, Astrophys. J., 892, 96
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  * Schneider et al. (2016) Schneider R., Hunt L., Valiante R., 2016, Mon. Not. R. Astron. Soc., 457, 1842
  * Schurer et al. (2009) Schurer A., Calura F., Silva L., Pipino A., Granato G. L., Matteucci F., Maiolino R., 2009, MNRAS, 394, 2001
  * Slavin et al. (2020) Slavin J. D., Dwek E., Mac Low M.-M., Hill A. S., 2020, Astrophys. J., 902, 135
  * Takeuchi et al. (2003) Takeuchi T. T., Hirashita H., Ishii T. T., Hunt L. K., Ferrara A., 2003, Mon. Not. R. Astron. Soc., 343, 839
  * Takeuchi et al. (2005) Takeuchi T. T., Ishii T. T., Nozawa T., Kozasa T., Hirashita H., 2005, MNRAS, 362, 592
  * Tamura et al. (2019) Tamura Y., et al., 2019, ApJ, 874, 27
  * Todini & Ferrara (2001) Todini P., Ferrara A., 2001, MNRAS, 325, 726
  * Valiante et al. (2009) Valiante R., Schneider R., Bianchi S., Andersen A. C., 2009, MNRAS, 397, 1661
  * Valiante et al. (2011) Valiante R., Schneider R., Salvadori S., Bianchi S., 2011, MNRAS, 416, 1916
  * Varosi & Dwek (1999) Varosi F., Dwek E., 1999, ApJ, 523, 265
  * Ventura et al. (2012) Ventura P., et al., 2012, MNRAS, 2357, 2345
  * Ventura et al. (2013) Ventura P., Criscienzo M. D., Carini R., Antona F. D., 2013, MNRAS, 3653, 3642
  * Watson et al. (2015) Watson D., Christensen L., Knudsen K. K., Richard J., Gallazzi A., Michałowski M., 2015, Nature, 519, 327
  * Winters et al. (1997) Winters J. M., Fleischer A. J., Le Bertre T., Sedlmayr E., 1997, A&A, 326, 305
  * Wolfire et al. (2003) Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, Astrophys. J., 587, 278
  * Woosley & Weaver (1995) Woosley S. E., Weaver T. A., 1995, Astrophys. J. Suppl. Ser., 101, 181
  * Yamasawa et al. (2011) Yamasawa D., Habe A., Kozasa T., Nozawa T., Hirashita H., Umeda H., Nomoto K., 2011, ApJ, 735
  * Yan et al. (2004) Yan H., Lazarian A., Draine B. T., 2004, Astrophys. J., 616, 895
  * Yasuda & Kozasa (2012) Yasuda Y., Kozasa T., 2012, ApJ, 745, 159
  * Zhukovska (2014) Zhukovska S., 2014, Astron. Astrophys., 562, 1
  * Zhukovska et al. (2008) Zhukovska S., Gail H. P., Trieloff M., 2008, A&A, 479, 453
  * da Cunha et al. (2010) da Cunha E., Eminian C., Charlot S., Blaizot J., 2010, Mon. Not. R. Astron. Soc., 403, 1894