11institutetext:  Sterrenkundig Observatorium Department of Physics and
Astronomy Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium
22institutetext: Department of Physics, Graduate School of Science, The
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
33institutetext: Institute of Space and Astronautical Science, Japan Aerospace
Exploration Agency, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210,
Japan 44institutetext: Institute of Astronomy and Astrophysics, Academia
Sinica, Astronomy-Mathematics Building, No. 1, Section 4, Roosevelt Road,
Taipei 10617, Taiwan 55institutetext: Theoretical Astrophysics, Department of
Earth and Space Science, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka
560-0043, Japan 66institutetext: Kavli IPMU (WPI), The University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan 77institutetext: Department
of Physics and Astronomy, University of Nevada, Las Vegas, 4505 S. Maryland
Pkwy, Las Vegas, NV 89154-4002, USA 88institutetext: Dept. Física Teórica y
del Cosmos, Universidad de Granada, Spain 99institutetext: Universitäts-
Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München,
Scheinerstr. 1, D-81679 München, Germany 1010institutetext: Max-Planck-
Institut für extraterrestrische Physik, Giessenbachstr. 1, D-85741 Garching,
Germany 1111institutetext: Excellence Cluster ORIGINS, Boltzmannstr. 2,
D-85748 Garching, Germany 1212institutetext: Instituto Universitario Carlos I
de Física Teórica y Computacional, Universidad de Granada, 18071, Granada,
Spain 1313institutetext: STAR Institute, Université de Liège, Quartier Agora,
Allée du six Aout 19c, B-4000 Liege, Belgium 1414institutetext: Physics
Department, Ben-Gurion University of the Negev, Be’er-Sheva 84105, Israel

# Observational signatures of the dust size evolution in isolated galaxy
simulations

Kosei Matsumoto , 112233<EMAIL_ADDRESS>Hiroyuki Hirashita 4455
Kentaro Nagamine 556677 Stefan van der Giessen 1188 Leonard E. C. Romano
9910101111 Monica Relaño 881212 Ilse De Looze 11 Maarten Baes 11 Angelos
Nersesian 111313 Peter Camps 11
Kuan-chou Hou 1414 Yuri Oku 55

(Received )

###### Abstract

Aims. We aim at providing observational signatures of the dust size evolution
in the ISM, in particular exploring indicators of polycyclic aromatic
hydrocarbon (PAH) mass fraction ($q_{\mathrm{PAH}}$) defined as the mass
fraction of PAHs relative to total dust grains. Additionally, we validate our
dust evolution model by comparing the observational signatures from our
simulations to those from observations.

Methods. We model dust properties of Milky Way-like and NGC 628-like galaxies
representing star-forming galaxies with a hydrodynamic simulation code,
GADGET4-OSAKA, which incorporates the evolution of grain size distribution
driven by dust production and interstellar processing. Furthermore, we perform
post-processing dust radiative transfer calculations with SKIRT based on the
hydrodynamic simulations to predict the observational properties of the
simulations.

Results. We find that the intensity ratio between 8 $\rm\mu m$ and 24 $\rm\mu
m$ ($I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$) correlates with
$q_{\mathrm{PAH}}$ and can be used as an indicator of PAH mass fraction.
However, this ratio is influenced by the local radiation field. As another
indicator of PAH mass fraction, we suggest the 8 $\rm\mu m$-to-total infrared
intensity ratio ($\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$), which tightly correlates with the PAH mass fraction.
Furthermore, we explore the spatial evolution of the PAH mass fraction in the
simulated Milky Way-like galaxy using $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$. We find that the spatially resolved PAH mass fraction
increases with metallicity in a lower metallicity regime ($Z\lesssim 0.2\
Z_{\odot}$) due to the interplay between accretion and shattering while it
decreases with metallicity in a high metallicity regime ($Z\gtrsim 0.2\
Z_{\odot}$) due to coagulation process. Also, the coagulation decreases the
PAH mass fraction in regions with a high hydrogen surface density. Finally, we
compare the above indicators in the NGC 628-like simulation with those
observed in NGC 628 by Herschel, Spitzer, and JWST. Consequently, we find that
our simulation underestimates the PAH mass fraction throughout the entire
galaxy. This is probably because PAH is too efficiently lost by coagulation in
the interstellar medium in our model, which suggests that the inhibition of
coagulation of the PAHs in the ISM is a key to enhancing PAH formation in our
model.

###### Key Words.:

Radiative transfer – ISM: molecules – Dust extinction – Methods: numerical

##  1 Introduction

Dust in the interstellar medium (ISM) plays a key role in shaping the observed
properties of galaxies as well as indirectly affecting galaxy evolution. Dust
governs various radiative processes in galaxies: dust absorbs and scatters
ultraviolet (UV) and optical light from young stars, and reemits that light as
infrared (IR) radiation. These processes are reflected in the spectral energy
distributions (SEDs) of galaxies; thus, analysis of the UV-to-IR SEDs of
galaxies provides various clues to their physical properties such as star
formation rate (SFR) and dust mass (e.g., Conroy 2013). In addition, dust
grains provide favorable conditions for the formation of hydrogen molecules
through the grain-surface reaction, and shielding of UV radiation from stars
by dust grains suppresses the dissociation of molecular hydrogen (Hollenbach &
Salpeter 1971; Hollenbach & McKee 1979; Cazaux & Tielens 2004; Grieco et al.
2023). Such a molecular-rich environment is favorable for star formation, with
dust acting as a cooling channel in molecular clouds (Ostriker & Silk 1973;
Vogelsberger et al. 2019). On the contrary, the photoelectric heating by dust
grains is the main heating source of the photodissociation and X-ray
dissociation regions in the ISM, changing the ISM conditions and causing
chemical reactions (Tielens & Hollenbach 1985; Maloney et al. 1996; Wolfire et
al. 2022). Therefore, dust is a crucial component of galaxy formation and
evolution.

To derive dust properties in the nearby universe, many observational studies
have been performed by analyzing SEDs with IR photometric data after the
appearance of Herschel and Spitzer telescopes, which cover wavelengths of
3–670 $\rm\mu m$ (e.g., da Cunha et al. 2008; Leja et al. 2017; Nersesian et
al. 2019; Burgarella et al. 2020). Observed galaxy IR SEDs reflect the grain
temperatures dependent on the grain size (including the stochastic heating;
Draine & Anderson 1985) and the material-dependent emission features at
certain wavelengths. In particular, the most prominent features seen at
$\sim$3–20 $\rm\mu m$ are considered to originate from PAH bands (Li & Draine
2001; Draine & Li 2007; Li 2020). Thus, the SED fitting analysis is capable of
extracting the grain size distribution and the PAH abundance (Galliano et al.
2008b, 2021; Relaño et al. 2020, 2022). Rémy-Ruyer et al. (2015) reveal an
increasing trend of PAH mass fraction, which is defined as the mass fraction
of PAHs relative to total dust grains, with respect to metallicity (see also
Draine et al. 2007; Aniano et al. 2020; Shim et al. 2023). Seok et al. (2014)
suggest that PAHs are produced by the shattering of carbonaceous grains, which
form more efficiently by the accretion of metals on dust grains in regions
with higher metallicity, naturally explaining the increasing trend of PAH mass
fraction against metallicity. Chastenet et al. (2019) investigate the PAH mass
fraction with spatially resolved data of the Small and Large Magellanic Clouds
(SMC and LMC, respectively) and find a higher PAH mass fraction in the diffuse
medium in the LMC compared to that in the SMC. To explain this result, they
suggest shattering large grains to PAHs is preferable to other PAH formation
scenarios (e.g., PAH formation in the envelope of asymptotic giant branch
(AGB) stars; Cherchneff et al. 1992; Sloan et al. 2008, 2014; Matsuura et al.
2014) and that the shattering works more efficiently in the LMC with higher
metallicity. Therefore, these observations using PAH features provide us clues
to understanding dust processes in various ISM.

In the SED fitting methods, however, the spatial resolutions of the
photometric data are always reduced to the worst one, compromising the spatial
dust characteristics obtained. Photometric and spectroscopic observations with
only a few prominent PAH emission features are useful to derive more detailed
spatial PAH properties (mass, size, and ionization), offering a more
comprehensive understanding of dust processes in the ISM (Li 2020). Today, the
James Webb Space Telescope (JWST) starts to investigate spatially resolved PAH
features of several nearby galaxies. Egorov et al. (2023) suggest the
intensity ratio of the emission at 7.7 and 11.3 $\rm\mu m$ to the 21 $\rm\mu
m$ emission ($R_{\mathrm{PAH}}$) for an indicator of the PAH mass fraction
given by

$R_{\mathrm{PAH}}=\frac{I_{\nu}(\mathrm{7.7\,\rm\mu
m})+I_{\nu}(\mathrm{11.3\,\rm\mu m})}{I_{\nu}(\mathrm{21\,\rm\mu m})}.$ (1)

They find that $R_{\mathrm{PAH}}$ decreases with the ionization fraction and
implies that PAH destruction by UV radiation works in the star-forming
regions. Chastenet et al. (2023b) exhibit a decreasing trend of the 3.35 to
11.3 $\rm\mu m$ intensity ratio against molecular gas fraction, implying that
larger grains are more abundant in denser regions. Although those observations
imply the spatially resolved dust processes in the ISM, it is difficult to
quantitatively infer the actual dust properties such as mass and grain sizes.
Therefore, to interpret those observations and discuss the dust processes in
ISM, it is essential to theoretically model dust evolution and compare
observational signatures of the dust evolution of the model to those of
observations.

In the past decades, many studies have been conducted on dust evolution using
physically motivated one-zone models (e.g., Lisenfeld & Ferrara 1998; Dwek
1998; Hirashita 1999). Essentially, it calculates grain size distributions for
predictions of dust signatures in extinction curves and dust emission in SEDs.
Some models describe the evolution of grain size distribution (e.g., O’Donnell
& Mathis 1997; Asano et al. 2013; Hirashita & Aoyama 2019), which have indeed
been used to predict extinction curves. Hirashita et al. (2020); Nishida et
al. (2022) investigate the evolution of the dust emission SED based on the
model of grain size distribution.

Today, these evolution models of grain size distributions, including a
simplified two-size approximation (Hirashita 2015), have been applied to
hydrodynamical simulations (McKinnon et al. 2018; Aoyama et al. 2017, 2020;
Hou et al. 2019; Graziani et al. 2020; Granato et al. 2021; Romano et al.
2022a; Narayanan et al. 2023) to calculate the evolution of dust content in
galaxies in a manner consistent with the hydrodynamic evolution of the ISM.
Among them, Aoyama et al. (2020) solve the evolution of grain size
distribution represented by 32-grain radius bins in an isolated galaxy
simulation with a hydrodynamic simulation code, GADGET3-OSAKA, and find that
the grain size distribution depends on the physical condition (especially
density and metallicity) of ambient ISM. This result underlines the importance
of solving the grain size distribution in a manner consistent with the
hydrodynamical and chemical evolution of the ISM. Romano et al. (2022a) extend
their models with more sophisticated sub-grid models, in which the fraction on
a sub-grid scale is varied according to the ambient density, and dust
diffusion between SPH gas particles (hereafter, gas particles) is taken into
account. In particular, they show that the inhomogeneity in the grain size
distribution in the ISM is strongly affected by the treatment of diffusion.
These simulations broadly succeed in reproducing the Milky Way extinction
curve and predicting the evolution of grain size distribution, but they have
not studied the influences of the radiative processes on the observational
properties based on the stellar and dust distribution of the hydrodynamic
simulations.

Narayanan et al. (2023) employ hydrodynamical simulations and post-processing
radiative transfer calculations, which handle the grain size distributions and
dust composition including PAHs, and investigate characteristics of dust
emission in different types of galaxies. They find a significant influence of
the radiation field on the ratio between global PAH luminosity and mass,
inferring the importance of the radiative processes to derive actual dust
properties in observations. However, it is still uncertain how the effect of
radiative processes can be canceled out in observations, and observational
signatures for deriving dust processes, which are not affected by the
radiation field, are unknown. Therefore, it is crucial to investigate how
grain size distributions are reflected on IR spectral features in simulations
and infer observational signatures, which can be directly compared with those
in observations.

In this paper, we model the evolution of grain size distribution in isolated
galaxies with hydrodynamic simulations (Romano et al. 2022a) and perform
radiative transfer calculations by taking into account the grain size
distributions and different grain species such as silicate, carbonaceous dust
grains, and PAHs in post-processing. Here, we assume that the main driver of
PAH formation is the shattering of large carbonaceous grains, which is implied
by both observational and theoretical studies (Chastenet et al. 2019; Seok et
al. 2014; Rau et al. 2019; Hirashita & Murga 2020; Narayanan et al. 2023). We
focus on two isolated galaxy simulations representing star-forming galaxies
similar to the Milky Way and NGC 628. The Milky Way-like galaxy is well-
studied in different types of simulations as a typical star-forming galaxy in
explaining the impact of their subgrid models on the galaxy evolution (Kim et
al. 2016; Shimizu et al. 2019; Romano et al. 2022a, b). The Milky Way-like
galaxy encapsulates a broad range of gas densities and metallicities,
rendering it an apt experimental setup to probe the long-term dust evolution
and the impact of the evolution on observable properties. On the other hand,
NGC 628 is among the most extensively studied galaxies in observations of the
nearby universe, offering high spatial resolutions in the PHANGS surveys
(Kennicutt et al. 2011; Leroy et al. 2021; Lee et al. 2023). NGC 628 is one of
the most nearby galaxies and a nearly face-on galaxy, allowing us to
separately study dust processes in diffuse and dense ISM with minimal
interference from dust extinction. Furthermore, NGC 628 is a moderately star-
forming galaxy without the signature of the active galactic nucleus observed
in the JWST program 2107 (PI: J. Lee; Lee et al. 2023). Therefore, comparing
the observable properties of NGC 628 between the simulation and observations
serves as a robust test to validate the dust evolution model, particularly
within the context of star-forming galaxies. More comprehensive comparisons
including other nearby galaxies will be presented in our future work (van der
Giessen et al. in preparation).

The main goals of this paper are the following: (i) We explore what drives the
variations of grain size distributions and PAH abundances in different galaxy
environments. (ii) We spatially investigate how dust evolution is reflected in
SEDs and dust emission at various stages of galaxy evolution. (iii) We infer
and validate various indicators of PAH mass fractions in ISM by making spatial
comparisons with the actual PAH mass fraction in hydrodynamic simulations.
(iv) We validate our dust evolution model by spatially comparing the
indicators of the PAH mass fraction between the NGC 628-like galaxy simulation
and actual observations.

This paper is organized as follows. Section 2 describes the models in our
hydrodynamic simulations of GADGET4-OSAKA, post-processing radiative transfer
calculations by SKIRT, and the observational analysis. In Section 3, we
explore the differences of the dust evolution and PAH abundance between the
Milky Way- and NGC 628-like galaxies. In Section 4, we investigate the
evolution of the spectral features in SEDs of the Milky Way-like galaxy and
indicators of the PAH mass fraction using intensity ratios at various
wavelengths. In Section 5, we perform spatially resolved comparisons for the
indicators of the PAH mass fraction between observations and simulations, and
discuss additional physical processes required for our dust model. In Section
6, we give our conclusions.

##  2 Methods

###  2.1 Isolated galaxy simulations

Table 1: Properties of each component laid in the initial conditions of the Milky Way-like and NGC 628-like galaxy simulations | Milky Way-like |  | NGC 628-like |
---|---|---|---|---
Components | Mass | Particle number | Mass | Particle number
Gaseous disk | $8.6\times 10^{9}\ {\rm M_{\odot}}$ | $1.0\times 10^{6}$ | $3.5\times 10^{10}\ {\rm M_{\odot}}$ | $9.0\times 10^{5}$
Gaseous halo | $1.0\times 10^{9}\ {\rm M_{\odot}}$ | $4.0\times 10^{5}$ | $1.1\times 10^{10}\ {\rm M_{\odot}}$ | $3.0\times 10^{5}$
Dark matter halo | $1.0\times 10^{12}\ {\rm M_{\odot}}$ | $1.0\times 10^{6}$ | $1.0\times 10^{12}\ {\rm M_{\odot}}$ | $3.0\times 10^{5}$
Stellar disk | $3.4\times 10^{10}\ {\rm M_{\odot}}$ | $1.0\times 10^{6}$ | $2.5\times 10^{10}\ {\rm M_{\odot}}$ | $3.0\times 10^{5}$
Bulge stars | $4.3\times 10^{9}\ {\rm M_{\odot}}$ | $1.25\times 10^{5}$ | $9.2\times 10^{9}\ {\rm M_{\odot}}$ | $3.0\times 10^{4}$

To perform isolated galaxy simulations, we use a smooth particle hydrodynamic
(SPH) simulation code, GADGET4-OSAKA (Romano et al. 2022a, b), which is a
modified version of GADGET4 (Springel et al. 2021). GADGET4-OSAKA incorporates
the OSAKA feedback model along with a dust evolution model (Aoyama et al.
2020; Romano et al. 2022a, b). The integrated feedback and dust models allow
us to approach a self-consistent treatment of star formation, feedback
processes, and metal/dust formation within the context of hydrodynamic
evolution in galaxies. Here, we briefly summarize those models and refer the
reader to the references for more details.

The star formation and feedback processes are incorporated into our code as
subgrid prescriptions. Star formation occurs stochastically, guided by a star
formation efficiency of $\epsilon_{*}=0.05$ that follows the Kennicutt-Schmidt
law (Kennicutt 1998). This process is initiated when gas particles satisfy
certain conditions, specifically when $n_{\mathrm{H}}>20\ \mathrm{cm^{-3}}$
and $T_{\mathrm{gas}}<10^{4}$ K ($n_{\mathrm{H}}$ is the hydrogen number
density and $T_{\mathrm{gas}}$ is the gas temperature) are satisfied. The
Osaka feedback model distinctly addresses the feedback processes of Type Ia
supernovae (SNe), Type II SNe, asymptotic giant branch (AGB) stars, and young
stars. The model manages energy ejections from both thermal and kinetic
feedback modes (Shimizu et al. 2019). Additionally, SNe and AGB stars eject
various metal elements and the ejected metal mass is calculated using the
CELib library (Saitoh 2017). Moreover, we utilize the GRACKLE-3 radiative
cooling library (Smith et al. 2017)111https://grackle.readthedocs.org/, which
solves the non-equilibrium primordial chemical network (H, He, and D), taking
into account metal cooling as well as photo-heating and photo-ionization due
to the UV background (Haardt & Madau 2012).

The dust evolution model adopted by Romano et al. (2022a) is an enhanced
version of those used by Hirashita & Aoyama (2019) and Aoyama et al. (2020).
We calculate the grain size distribution in each gas particle. This implicitly
assumes a perfect dynamic coupling between dust and gas within each gas
particle. The grains are assumed to be compact spheres with material density
$s$, such that the grain mass $m$ and radius $a$ are related as
$m=\frac{4}{3}\pi a^{3}s$. The grain size distribution, denoted as $n(a)$, is
defined so that $n(a)\,\mathrm{d}a$ is the number density of grains whose
radii are between $a$ and $a+\mathrm{d}a$. We consider the grain radii ranging
from $3.0\times 10^{-4}$ to 10 $\mathrm{\mu}$m with 30 logarithmically spaced
discrete bins. In our simulations, AGB stars and SNe produce large dust
grains. Here, 10% of the ejected metals from AGB stars and SNe are assumed to
condense into dust grains with a log-normal grain size distribution
characterized by the mean grain radius of 0.1 $\rm\mu m$ and variance of 0.47
(Asano et al. 2013). Our model also takes into account a variety of
interstellar dust processes including shattering and thermal sputtering in the
diffuse ISM and dust size growth through accretion and coagulation in the
dense ISM (Aoyama et al. 2020). Dust destruction in SN shocks and metal/dust
diffusion between the diffuse and dense ISM are also considered.

Furthermore, We note that our simulations do not have sufficient spatial
resolutions for cold and dense gas ($T_{\mathrm{gas}}<50$ K and
$n_{\mathrm{H}}>10^{3}$ cm-3), which could cause inefficient coagulation and
accretion processes. Thus, Romano et al. (2022a) employ a two-phase ISM
subgrid model for gas particles, where it is assumed that dense and diffuse
gas media co-exist in each gas particle. The mass fraction of the dense gas
medium is given by

$f_{\mathrm{dense}}=\mathrm{min}\,\Bigg{(}\alpha\,\frac{n_{\mathrm{H}}}{1.0\
\mathrm{cm}^{-3}},1.0\Bigg{)},$ (2)

where we adopt $\alpha=0.12$. While the temperature and density of the dense
gas medium are assumed to be constant ($T_{\mathrm{dense}}=50$ K and
$n_{\mathrm{H,\ dense}}=10^{3}$ cm-3), we estimate the temperature and density
of the diffuse gas medium in each gas particle by considering the conservation
of the internal energy in each gas particle. This subgrid model allows us to
treat dust processes in the diffuse and dense ISM separately according to
their respective properties within each gas particle. Further details of the
separation treatment of the two-phase ISM are described in Romano et al.
(2022a), and the parameters operating the dust processes in the diffuse and
dense ISM are described in Aoyama et al. (2020).

Using these models, we simulate Milky Way-like and NGC 628-like galaxies.
While the initial condition for the Milky Way-like galaxy is taken from the
AGORA project222https://sites.google.com/site/santacruzcomparisonproject/ (Kim
et al. 2016), that for the NGC 628-like galaxy is created by
DICE333https://bitbucket.org/vperret/dice/src/master/ (Perret 2016). To create
the initial condition of the NGC 628-like galaxy, we need the morphological
parameters of stellar and gas distributions in NGC 628. Thus, we compute those
parameters based on the observationally obtained radial profiles of stars and
gas (for the details on creating initial conditions, see Appendix A).
Moreover, we set dark matter and gaseous halos, bulge stars, and stellar and
gaseous disks in the initial conditions. The properties of these geometries
are summarized in Table 1.

We run the Milky Way-like and NGC 628-like galaxy simulations until the
simulation time of 10 Gyr when the metallicities in both simulations are
consistent with the observed ones (Berg et al. 2015, for NGC 628). Here, we
adopt an initial metallicity of $0.0001\ Z_{\odot}$, where $Z_{\odot}$ is the
solar metallicity ($Z_{\odot}=0.013$; Asplund et al. 2009), to start our
simulations in the condition of the early universe. For the hydrodynamic
simulations, we permit the minimum smoothing length to reach $0.1\
\epsilon_{\mathrm{grav}}$, where $\epsilon_{\mathrm{grav}}$ represents the
gravitational softening length. The softening lengths of the Milky Way-like
and NGC 628-like galaxy simulations are 40 and 50 pc, respectively.

###  2.2 Post-processing radiative transfer with SKIRT

We employ the state-of-the-art three-dimensional (3D) Monte Carlo radiative
transfer code, SKIRT444www.skirt.ugent.be. (Baes et al. 2011; Camps & Baes
2015, 2020), to produce the SEDs and images from UV to millimeter wavelengths.
While SKIRT is evolving into a generic radiative transfer code with
capabilities to handle many radiative processes relevant to dust and gas media
(e.g., Camps et al. 2021; Gebek et al. 2023; Vander Meulen et al. 2023; Kapoor
et al. 2023; Matsumoto et al. 2023), here we consider extinction, emission,
scattering, and self-absorption only by dust. We include stochastic heating of
dust grains, which is particularly important for small grains and PAHs (Camps
et al. 2015). The code is set up to deal with dust mixtures with spatially
varying grain size distributions and optical properties555Configuring custom
dust mixes is described in
https://skirt.ugent.be/root/_tutorial_custom_dust.html. SKIRT has been
optimized to post-process snapshots from hydrodynamical galaxy simulations and
contains mechanisms to read in SPH, adaptive mesh refinement, and moving-mesh
simulations (Saftly et al. 2013, 2014; Camps et al. 2013; Baes & Camps 2015).
The methodology to set up SKIRT simulations based on various types of
hydrodynamical simulations is described in detail in Camps et al. (2018,
2022), Kapoor et al. (2021), Trčka et al. (2022), and Baes et al. (2024). In
the following, we describe the details of the post-processing models of the
stellar sources and dust composition used in our SKIRT simulations:

Stellar populations:

We utilize the Bruzual & Charlot (2003) templates for single stellar
populations with the Chabrier (2003) initial mass function. The intrinsic SED
of each star particle is determined based on its stellar age, metallicity, and
initial stellar mass. We note that, under the initial conditions of the
hydrodynamic simulations, the disk and bulge star particles are pre-assigned
to make the galaxy dynamically stable. This means that some stars already
exist at the simulation time of $t=0$ Gyr. Thus, we need to construct the star
formation history (SFH) before $t=0$ Gyr for the disk and bulge star
particles. Since we run the simulations until $t=10$ Gyr and assume the galaxy
age to be comparable to the cosmic age ($13.8$ Gyr), the disk and bulge stars
have formed between $t=-3.8$ and $t=0$ Gyr. Therefore, in post-processing, we
randomly assign their ages to obtain a constant SFH in the $3.8$ Gyr interval
before $t=0$ Gyr. The metallicity of the disk and bulge star particles is the
same as the initial metallicity of $0.0001\ Z_{\odot}$.

Dust models:

We apply the following dust composition model to each gas particle containing
dust grains in the post-processing. First, larger dust grains ($a>1.3$ nm) are
designated as silicate and carbonaceous dust grains. The grain size boundary
of $a=1.3$ nm between smaller and larger grains is taken from the definition
of the maximum size of PAHs in Draine & Li (2001) and Draine et al. (2007).
The mass fractions of silicate and carbonaceous dust (denoted as
$f_{\mathrm{silicate}}$ and $f_{\mathrm{carbon}}$, respectively) in each gas
particle are estimated based on the abundance of Si and C within the gas
particle following Hirashita & Murga (2020) as

$\displaystyle f_{\mathrm{silicate}}$ $\displaystyle=$
$\displaystyle\frac{6Z_{\mathrm{Si}}}{6Z_{\mathrm{Si}}+Z_{\mathrm{C}}},$ (3)
$\displaystyle f_{\mathrm{carbon}}$ $\displaystyle=$
$\displaystyle(1-f_{\mathrm{silicate}}),$ (4)

where $Z_{\mathrm{Si}}$ and $Z_{\mathrm{C}}$ represent the mass abundances of
Si and C within each gas particle, respectively. The factor of 6 is obtained
from the mass fraction of Si in silicate, as described in Hirashita & Kuo
(2011) and Hirashita & Murga (2020). As shown by Hirashita et al. (2020), our
model tends to underpredict PAH emission (see also Chang et al. 2022). To
improve our predictions in a maximally acceptable manner within our model
framework, we assume that all the smaller grains ($a\leq 1.3$ nm) are
carbonaceous. This is consistent with the grain model obtained by Weingartner
& Draine (2001) as a result of fitting to the Milky Way extinction curve: in
their model, the grain size distribution of silicate has a steeper decline
toward small grains while that of carbonaceous dust is enhanced at small
radii. As a consequence of this assumption, the grains with $a\leq 1.3$ nm are
composed of PAH0, PAH+, and non-aromatic carbonaceous grains. Since the
aromatization and its inverse reaction predominantly occur in the diffuse and
dense media, respectively, the mass fractions of aromatic (PAH0 and PAH+)
carbonaceous grains correspond to the diffuse gas fraction
($f_{\mathrm{aroma}}=1-f_{\mathrm{dense}}$; Hirashita & Murga 2020). In the
analysis of this paper, the PAH mass fraction is given by

$q_{\mathrm{PAH}}=\frac{M_{\mathrm{PAH\ (a\leq 1.3\
\mathrm{nm})}}}{M_{\mathrm{(total\ dust)}}}\times 100.$ (5)

Finally, to estimate the mass fraction for the PAH0 and PAH+ in the aromatic
carbonaceous dust grains of each gas particle, we use the ionization fraction
from the “standard” ionization model in Draine et al. (2021) and Hensley &
Draine (2023):

$\displaystyle f_{\mathrm{ion}}(a)=1-\frac{1}{1+a/10\,\AA}.$ (6)

The opacities and calorimetric properties of PAHs, silicate, and carbonaceous
grains (assumed to be graphite) follow the prescription of Draine & Li (2007).

Instrument setups of the Milly Way-like galaxy:

We observe the galaxy at a distance of 10 Mpc and inclination angles of 0 and
90 ∘. For generating the synthetic observational data, we set the fields of
views to ($40$ kpc)2 and adopt a pixel size of 50 pc, which is close to the
gravitational softening length of the Milky Way-like galaxy simulation. We
select pixels of each band data with an intensity ($I_{\nu}>0.01$
$\mathrm{MJy\,sr^{-1}}$) to avoid pixels with significant Monte Carlo noise in
SKIRT.

Instrument setups of the NGC 628-like galaxy:

NGC 628 is located at a distance of 9.77 Mpc, and the inclination angle is 8.9
∘ (Leroy et al. 2021). Thus, the simulated galaxy is observed under the same
conditions as the actual galaxy. For comparing this simulation with the actual
observations in Section 5, we use the fields of views of ($20$ kpc)2 and a
pixel size of 50 pc, which corresponds to the gravitational softening length
of the NGC 628-like galaxy simulation. The pixels of each band data with an
intensity ($I_{\nu}>0.01$ $\mathrm{MJy\,sr^{-1}}$) are selected. When we
compare the simulation data with the Herschel and Spitzer data, the
resolutions of the simulation data are reduced to 600 pc to make the
consistency of the beam size of the Herschel PACS $160$ $\rm\mu m$ for NGC 628
(see also Section 2.3).

###  2.3 Observational data and analysis

For the observational photometric data of NGC 628, we use the Dustpedia
database (Davies et al. 2017; Clark et al.
2018)666http://dustpedia.astro.noa.gr, the Spitzer Local Volume Legacy survey
(Dale et al. 2009), and the PHANGS (Physics at High Angular resolution in
Nearby GalaxieS) survey (Lee et al.
2023)777https://sites.google.com/view/phangs/home/data for Herschel, Spitzer,
and JWST data, respectively. More specifically, we use the Herschel PACS $70$
and $160\leavevmode\nobreak\ \rm\mu m$ bands, the Spitzer IRAC
$8\leavevmode\nobreak\ \rm\mu m$ and MIPS $24\leavevmode\nobreak\ \rm\mu m$
bands, and JWST MIRI $7.7$, $11.3$, and $21\leavevmode\nobreak\ \rm\mu m$
bands. After we obtain the data, we convolve each image with the kernel
function of the point spread function of each instrument. The kernel functions
are taken from the Astronomical Convolution-kernel repository (Aniano et al.
2011)888https://www.astro.princeton.edu/~draine/Kernels.html. To accomplish
the comparison between the observations and the NGC 628-like galaxy simulation
in Section 5, we averaged the observational data so that the spatial
resolutions of JWST images are reduced to 50 pc, which corresponds to the
softening length of the simulation. The resolutions of the Herschel and
Spitzer images are unified to 600 pc, which corresponds to the beam size of
the Herschel PACS $160$ $\rm\mu m$. After these processes, we select the
averaged pixels with a 1 sigma level of the original photometric data.

##  3 Overall evolution of the galaxy properties of the isolated galaxy
simulations

Figure 1: Time evolution of global physical properties of the Milky Way-like (
lines) and NGC 628-like galaxies (blue solid and orange dashed lines,
respectively). The evolution of (a) SFR, (b) dust mass, (c) dust-to-metal mass
ratio, and (d) PAH mass fraction are shown.

We provide overviews of the simulated galaxies and show the time evolution of
global quantities that are important to understand the dust evolution for the
Milky Way-like and NGC 628-like cases.

Figure 1(a) shows the SFH of the Milky Way-like and NGC 628-like galaxies. In
the Milky Way-like galaxy, a large amount of gas assumed in the galactic
center starts to contract, leading to explosive star formation in the early
stages of the simulation (hereafter, referred to as the starburst phase).
Subsequently, star formation and its feedback gradually deplete the gas in the
center of the galaxy, resulting in a decrease in the SFR. In contrast, star
formation occurs moderately in the NGC 628-like galaxy, since the gas radial
profile of the NGC 628-like galaxy is flat and not dense enough to cause
explosive star formation. Thus, the gas is slowly consumed by star formation
with $\mathrm{SFR}=1.0$ $\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$ over
10 Gyr.

The different SFHs between the Milky Way-like and NGC 628-like galaxies lead
to their different dust evolution. The total dust mass increases more rapidly
for the Milky Way-like galaxy at $t<1.0$ Gyr than for the NGC 628-like galaxy,
as shown in Fig. 1(b) since more dust is produced in the starburst phase of
the Milky Way-like galaxy.

The dust-to-metal mass ratio gives us insights into the efficiency of dust
accretion in the ISM since dust grains grow by acquiring gas-phase metals
(e.g., Hirashita 1999; Chiang et al. 2018). Figure 1(c) shows that the dust-
to-metal mass ratio of the Milky Way-like galaxy increases more steeply than
that of the NGC 628-like galaxy since the quicker metal enrichment in the
starburst phase of the Milky Way-like galaxy leads to more efficient accretion
of gas-phase metals onto dust grains.

Figure 1(d) shows the evolution of the PAH mass fractions, which indicates the
efficiency of shattering since PAHs represent small grains. The PAH mass
fraction in the Milky Way-like galaxy less steeply increases since the Milky
Way-like galaxy contains more amounts of dense gas than the NGC 628-like
galaxy, and hence, shattering is less efficient. The difference in the ISM
conditions between these galaxies also reflects their grain size distributions
at $t=10$ Gyr as shown in Fig. 2. The mean grain size distributions have
different slopes between the two simulations in such a way that the Milky Way-
like galaxy (blue solid line) is more dominated by large grains than the NGC
628-like galaxy (orange dashed line). At $t=10$ Gyr, the shape is determined
by the balance between coagulation and shattering. Coagulation turns more
small grains into large grains in the Milky Way-like galaxy, since the galaxy
contains more amounts of dense gas than the NGC 628-like galaxy. These results
suggest that different ISM conditions and SFHs cause different dust evolution.

Figure 2: Grain size distribution in the Milky Way-like (blue solid line) and
NGC 628-like (orange dashed line) galaxies at $t=10$ Gyr. The solid lines and
the shades show the mean values and the range of 10th–90th percentiles for the
sample of gas particles, respectively.

##  4 Observational signatures of the grain size evolution in the Milky Way-
like galaxy

To investigate observational signatures of the dust size evolution, we focus
on the results from the Milky Way-like galaxy simulations. We note that the
following discussion in this section does not significantly change even if we
involve the results of the NGC 628-like galaxy simulation. Also, to scrutinize
the ISM processes across the whole galaxy with minimal interference from dust
extinction, we use the synthetic observation data of a face-on view of the
Milky Way-like galaxy in Section 4.2.

###  4.1 Influences of the dust size evolution on dust emission

Figure 3: Evolution of the global SED for the Milky Way-like galaxy
simulation. The different colored lines correspond to the total SEDs at
$t=0.05,\ 0.5,\ 1,\ 3,$ and $10$ Gyr (cyan, orange, green, red, and purple
lines, respectively). The observed flux density is estimated at a distance of
10 Mpc. Figure 4: Evolution of stellar (first and second rows) and dust
emission maps (third and fourth rows) of the Milky Way-like galaxy simulation.
The panels show the emission maps at $t=0.05,\ 0.5,\ 1,\ 3,$ and $10$ Gyr from
left to right. The images on the first and third (second and fourth) rows are
face-on (edge-on) views. The stellar emission maps consist of the SDSS $i$
(red), $r$ (green), and $g$ (blue) fluxes with an additionally enhanced blue
level for the GALEX FUV flux. The dust emission maps are synthesized by the
Spitzer IRAC $8\leavevmode\nobreak\ \rm\mu m$ (blue), Herschel PACS
$70\leavevmode\nobreak\ \rm\mu m$ (green), and Herschel PACS
$160\leavevmode\nobreak\ \rm\mu m$ (red) fluxes.

We investigate how dust evolution is reflected in SEDs at various stages of
galaxy evolution. Figure 3 shows the total SEDs of the Milky Way-like galaxy
simulation at $t=0.05,\ 0.5,\ 1,\ 3$, and $10$ Gyr. At $t=0.05$ Gyr, only
emission of large dust grains can be observed around $\lambda=70$ $\rm\mu m$,
as dust production is dominated by stellar sources, which predominantly
produce large grains. At $t\geq 0.5$ Gyr, emission from small dust grains and
PAHs become noticeable around $\lambda=10\leavevmode\nobreak\ \rm\mu m$. This
is due to shattering, which acts as a source of small grains, and accretion,
which increases the abundance of small grains (see also Asano et al. 2013;
Hirashita & Aoyama 2019). Given the monotonic decrease in the SFR in the Milky
Way-like galaxy at later epochs, the intrinsic radiation from young stars at
UV to optical wavelengths is weakened with time. Consequently, the peak of
dust emission in the SED shifts from $\lambda=70$ to $150$ $\rm\mu m$ over
time. This shift corresponds to a decrease in dust temperature from
$T_{\mathrm{dust}}=40$ to $20$ K. On the other hand, the dust mass
monotonically increases with time in the Milky Way-like galaxy through mainly
the accretion of gas-phase metals. Thus, at $t=1$ Gyr, the total IR luminosity
achieves its maximum.

Figure 4 shows stellar and dust emission maps from the Milky Way-like galaxy
simulation at the simulation times of $t=0.05,\ 0.5,\ 1.0,\ 3.0,$ and $10$
Gyr. The first two rows show stellar emission with colors synthesized using
the SDSS $i$ (red), $r$ (green), and $g$ (blue) fluxes and the GALEX FUV flux
(additionally enhanced blue color) as seen from face-on and edge-on. As viewed
face-on, old stars radiate throughout the entire galaxy (appearing in white),
while UV radiation from young stars (blue) originates from star-forming
regions, which tend to be concentrated in dense clumps along the spiral arms.
As viewed face-on and edge-on, the effect of dust attenuation is prominent at
later stages of the galaxy evolution when dust formation has proceeded enough
to shield the stellar radiation. Spiral-shaped and horizontally extended dark
lanes in the stellar emission develop from face-on and edge-on views,
respectively.

The bottom two rows of Fig. 4 show the dust emission maps consisting of
Spitzer IRAC $8\leavevmode\nobreak\ \rm\mu m$ (red), Herschel PACS
$70\leavevmode\nobreak\ \rm\mu m$ (green), and Herschel PACS
$160\leavevmode\nobreak\ \rm\mu m$ (blue) fluxes as seen from face-on and
edge-on. The IRAC ${8\leavevmode\nobreak\ \rm\mu m}$ band mainly traces the
emission from PAHs heated by radiation from young stars, while the PACS
${70\leavevmode\nobreak\ \rm\mu m}$ and PACS ${160\leavevmode\nobreak\ \rm\mu
m}$ bands trace emission from warm and cold large dust grains, respectively.
At $t=0.05$ Gyr, large grains are produced by stars in the central region of
the galaxy, and the emission from the warm large grains heated by young stars
is prominent. At $t\geq 0.5$ Gyr, large dust grains are produced by stars in
the spiral arms, and simultaneously accretion and shattering become prevalent
in the dense and diffuse ISM, respectively. These processes increase small
grains, leading to significant emission from small dust grains and PAHs in the
MIR (i.e., the color becomes bluer in the images). At $t\geq 1$ Gyr,
coagulation also becomes efficient in the dense ISM, enhancing emission from
large dust grains in the central regions and spiral arms of the galaxy. The
emission from small dust grains and PAHs is prominent in diffuse regions along
the spiral arms. At later times, as star formation activity declines, the
overall dust emission becomes weaker, and the color becomes redder because of
lower dust temperature. Finally, at $t=10$ Gyr, the functional shape of the
grain size distribution at each region converges to the one determined by the
balance between the shattering and coagulation rates (see also Fig. 2). Since
the overall density is higher in the central regions (thus, coagulation is
more enhanced), the emission from large and small grains is relatively
prominent in the inner and outer regions of the galaxy, respectively.

###  4.2 Indicator of the PAH mass fraction

We investigate promising indicators of PAH mass fractions in the ISM using
multi-wavelength intensity ratios in the SKIRT simulation of the Milky Way-
like galaxy. The validity of the adopted ratios as a good indicator is tested
against actual PAH mass fractions of the Milky Way-like galaxy simulation.
Since the PAH emission is the most prominent in the 8 $\rm\mu m$ band among
the bands available for Spitzer and JWST, we construct the indicators by
combining the 8 $\rm\mu m$ band and other wavelengths. We focus on the
following two: one is the intensity ratio using representative MIR bands, that
is, the 8 to 24 $\rm\mu m$ intensity ratio,
$I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu
m})$, and the other is the $8$ $\rm\mu m$ to the total infrared intensity
ratio, $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$. Here we
denote the intensity per frequency at wavelength $\lambda$ as
$I_{\nu}(\lambda)$. Also, we note that we investigate those indicators with
the spatially resolved information of the 50-pc pixel data constructed in
Section 2.2.

####  4.2.1 $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$

Figure 5: Pixel-based relation between the 8-to-24 $\rm\mu m$ intensity ratio
and the PAH mass fraction of the Milky Way-like galaxy simulation at $t=0.5$,
1.0, 3.0, and 10 Gyr from top left to bottom right. The color represents
density-weighted mean intensity at 1000 Å relative to the Milky Way value as
shown by the color bar on the top of each panel. Blue dashed lines indicate
the best-fit relation given by Eq. (7), and the legends show the values of
coefficients $k$ and $c$. Figure 6: Same as Fig. 5 but for $R_{\mathrm{PAH}}$.

Figure 5 shows the pixel-based relation between PAH mass fraction of
$q_{\mathrm{PAH}}$ and $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$ at $t=\ 0.5,\ 1.0,\ 3.0,$ and
$10$ Gyr. To investigate the effect of the radiation field (or dust heating),
we color-code the data points with the density-weighted mean intensity at 1000
Å relative to the Milky Way value of $10^{6}$ photons
$\mathrm{cm^{-2}s^{-1}eV^{-1}}$, $U_{\mathrm{MW}}$. We fit the data at each
time by the following equation:

$\mathrm{log}_{10}(\mathrm{Intensity\
Ratio})=k\,\mathrm{log}_{10}(q_{\mathrm{PAH}})+c,$ (7)

where $k$ and $c$ are coefficients derived from the fitting, and $k$ is
associated with the radiative processes of dust grains. The best-fit relation
at each time is indicated as blue dashed lines in Fig. 5. We find that
$I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu
m})$ correlates with $q_{\mathrm{PAH}}$ at each time. The 8 $\rm\mu m$ band is
almost purely dominated by the PAH emission, while the $24$ $\rm\mu m$ band is
contaminated by emission from larger dust grains with $a\gtrsim 20$ Å (Draine
& Li 2007), so that $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$ increases with
$q_{\mathrm{PAH}}$. However, we also find that stronger radiation fields
decrease the ratio at a fixed $q_{\mathrm{PAH}}$, as color indicates. This is
because stronger radiation fields cause higher dust temperatures, shifting the
large grain emission toward shorter wavelengths and raising emission at 24
$\rm\mu m$ (Draine & Li 2007; Schreiber et al. 2018). Consequently, given the
monotonic decrease in SFR with time in the Milky Way-like galaxy, the best-
fitting slope becomes shallower with time. The above dependence on the
radiation field, however, is not strong enough to break the correlation
between $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$ and $q_{\mathrm{PAH}}$ at each
time. Therefore, $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$ can be an indicator of the PAH
mass fraction with a tendency of underestimating the PAH mass fraction in
regions of strong radiation fields such as star-forming regions. In
particular, $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$ has the advantage of using only
two bands in the MIR, which is useful if the available wavelength coverage is
limited. We note that the above results are also applicable to an indicator of
PAH mass fraction, $R_{\mathrm{PAH}}$, for JWST as shown in Fig. 6.

####  4.2.2 $\nu I_{\nu}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$

Here, we use the total infrared intensity estimated by the following equation:

$I_{\mathrm{TIR}}=0.95\,\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})+1.15\,\nu
I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})+\nu
I_{\mathrm{\nu}}(\mathrm{70\,\rm\mu m})+\nu
I_{\mathrm{\nu}}(\mathrm{160\,\rm\mu m}),$ (8)

which is suggested by Draine & Li (2007). They theoretically indicate that the
estimated total infrared intensity differs from the actual one within the
relative difference of 10 $\%$ for the radiation fields of
$0.1<U_{\mathrm{MW}}<10^{2}$. Here, we expect to maintain higher spatial
resolutions of the observed data by utilizing the estimation with bands at up
to 160 $\rm\mu m$ although there are other estimation methods with longer
wavelength bands for the total infrared intensity (Boquien et al. 2011;
Galametz et al. 2013). Figure 7 shows the pixel-based relation between PAH
mass fraction and $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$ at $t=0.5,\ 1.0,\ 3.0,$ and $10$ Gyr. The color
represents the density-weighted mean intensity at 1000 Å relative to the Milky
Way value, $U_{\mathrm{MW}}$.We find that the ratio tightly correlates with
$q_{\mathrm{PAH}}$ and is not affected by the local radiation field of
$U_{\mathrm{MW}}$ unlike $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$. Remarkably, the parameters of
the best-fit relation (blue) are consistent among the different times within
the relative difference of 5%. This is because both total and 8 $\rm\mu m$
intensities are proportional to the strength of the local radiation field
(Draine & Li 2007), and hence, the effect of the radiation field is canceled
out.

Figure 7: Pixel-based relation between the 8 $\rm\mu m$-to-total IR intensity
ratio and the PAH mass fraction of the Milky Way-like galaxy simulation at
$t=0.5,\ 1.0,\ 3.0,$ and $10$ Gyr from top left to bottom right. The color
represents the density-weighted mean intensity at 1000 Å relative to the Milky
Way value as indicated by the color bar on the top of each panel. Blue dashed
lines indicate the best-fit relation of Eq. (7).

Additionally, we perform linear regression for all of the data including
snapshots at different times. We find that the PAH mass fraction of
$q_{\mathrm{PAH}}$ can be formulated by the following equation,

$q_{\mathrm{PAH}}=66\,\left(\frac{\mathrm{\nu}I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})}{I_{\mathrm{TIR}}}\right)^{1.2}.$ (9)

Therefore, we suggest that $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$ is a good indicator of the PAH mass fraction in the ISM
of galaxies in various stages of galaxy evolution.

####  4.2.3 Dependence of the PAH mass fraction on the metallicity and
hydrogen surface density

To demonstrate the usefulness of $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$ as an indicator of the PAH mass fraction, we investigate
the spatially resolved relation among the PAH mass fraction, gas-phase
metallicity, and hydrogen surface density. Appendix B shows the projected maps
of PAH mass fraction, metallicity, and hydrogen surface density.

Figure 8: Time evolution of pixel-based relation among the 8 $\rm\mu m$-to-
total IR intensity ratio, gas-phase metallicity, and hydrogen surface density
for the Milky Way-like galaxy simulation. The panels show the relation at
$t=0.5,\ 1.0,\ 3.0,$ and $10$ Gyr from top left to bottom right. The color bar
indicates the hydrogen surface density of each pixel. In each panel, the time
evolution of the entire distribution over the four epochs is also indicated by
the purple dashed lines. The purple star symbols represent the global
intensity ratio and metallicity of the galaxy at $t=0.5,\ 1.0,\ 3.0,$ and $10$
Gyr from the lower left to the upper right in each panel (the symbol
transparency decreases) and are connected by the dashed line. We note that the
intensity ratio and metallicity are the intensity- and mass-weighted values,
respectively. The grey line shows the evolutionary path in the one-zone dust
evolution model by Hirashita & Murga (2020) under the condition of
$f_{\mathrm{dense}}=0.5$, $\tau_{\mathrm{SF}}=5$ Gyr, and $U_{\mathrm{MW}}=1$.

Figure 8 shows the pixel-based relation between $\nu
I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ and metallicity at
$t=0.5,\ 1.0,\ 3.0,$ and $10$ Gyr. The color represents the hydrogen surface
density of each pixel. The purple star symbols represent the global intensity
ratio and metallicity of the galaxy at four epochs, and the dashed lines
indicate the time evolution of global properties over the epochs (from the
lower left to the upper right in each panel). We find that the global $\nu
I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ (PAH mass fraction)
increases with time and metallicity. This global trend is also suggested by
observations using many galaxy samples (Galliano et al. 2008a; Rémy-Ruyer et
al. 2015; Shim et al. 2023). Shattering and accretion happen in separate ISM
phases, and their interplay increases small grains (and PAHs) on a galactic
scale in our simulations (Hirashita 2012). Hence, the global PAH mass fraction
evolves with metallicity.

In the pixel-based analysis, the local PAH mass fraction is associated with
the local ISM conditions such as metallicity and hydrogen surface density. At
$t=0.5$ Gyr, $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$
increases with metallicity in a lower metallicity regime ($Z\lesssim 0.2\
Z_{\odot}$) while it decreases with metallicity in a higher metallicity regime
($Z\gtrsim 0.2\ Z_{\odot}$), as suggested by a one-zone dust evolution model
(grey line; Hirashita & Murga 2020). We note that the number of pixels in the
lower metallicity regime is significantly reduced by the pixel selection
because those with significant noises are excluded in our analysis in Section
2.3. At low metallicity, the interplay between accretion and shattering
increases the PAH mass fraction, producing a trend of increasing the PAH mass
fraction with metallicity. At higher metallicity, dust grows sufficiently
through metal accretion, and the effect of coagulation becomes prominent. The
coagulation turns small grains, including PAHs, into larger grains and plays
an important role in decreasing the PAH mass fraction at a high-metallicity
regime. At $\geq 1.0$ Gyr, $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$ anticorrelates with metallicity since most of the
regions reach high metallicity. We also find that $\nu
I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ depends on the
hydrogen surface density since coagulation converts small grains, including
PAHs, to larger grains in the dense ISM.

These results are roughly consistent with observational studies. Egorov et al.
(2023) report a weak anti-correlation between $R_{\mathrm{PAH}}$ and oxygen
abundance based on the spatially resolved observations of four galaxies with
high metallicity (12+log(O/H)$\gtrsim$8.4). Chastenet et al. (2023a) show a
tendency to increase $R_{\mathrm{PAH}}$ with oxygen abundance at low oxygen
abundance (12+log(O/H)$\lesssim$8.5) but flatten out at high oxygen abundance
(12+log(O/H)$\gtrsim$8.5). However, those observations cover only four
galaxies, and future observations are necessary to verify whether there are
general trends for the spatially resolved PAH mass fractions to increase and
decrease with metallicity at a lower and higher metallicity, respectively.

Therefore, the spatially resolved relation of the PAH mass fraction against
metallicity is related to the efficiency of coagulation, shattering, and
accretion across diverse ISM conditions, and the analysis of the spatially
resolved relation leads us to interpret dust evolution in observations.

##  5 NGC 628-like galaxy simulation

NGC 628 is an ideal galaxy to validate our dust model based on spatially
resolved data since it is well-observed with a high spatial resolution in the
PHANGS surveys (Kennicutt et al. 2011; Leroy et al. 2021; Lee et al. 2023).
Here, we model the NGC 628-like galaxy simulation according to the
morphological parameters of the NGC 628 (see Appendix A) to make a fair
comparison. We use the simulation snapshot at $t=10$ Gyr as a representative
age of present-day galaxies. More comprehensive comparisons including other
nearby galaxies will be presented in our future work (van der Giessen et al.
in preparation).

###  5.1 SED comparisons

Figure 9: Comparison of the global SED between the NGC 628-like galaxy
simulation and actual observations. Upper panel: The black solid line shows
the global SED of the NGC 628-like galaxy simulation at $t=10$ Gyr. Blue
dashed and red dash-dotted lines represent the intrinsic stellar and dust
emission, respectively. The orange and blue diamonds correspond to the
simulated and observed photomeries of the total SEDs (Relaño et al. 2020),
respectively. Lower panel: Difference in the photometries between the
simulation and observation data.

Figure 9 shows the total SED comparison between the simulation and the actual
observations of NGC 628 (Relaño et al. 2020, blue and orange diamonds,
respectively). The total, intrinsic stellar, and dust emission SEDs of the NGC
628-like simulation are also indicated by black solid, blue dashed, and red
dash-dotted lines, respectively. The lower panel indicates the difference in
the photometries between the simulation and the observations. Stellar emission
at optical and near-IR wavelengths of the simulated galaxy is consistent with
that of the observations while the emission at UV wavelengths in the
simulation is weaker than that of the observations. This is because the SFR of
the simulated galaxy at $t=10$ Gyr ($\dot{M}_{*}=0.75$
$\mathrm{M_{\odot}\leavevmode\nobreak\ yr^{-1}}$) is lower than that reported
by the observation ($\dot{M}_{*}=1.7$ $\mathrm{M_{\odot}\leavevmode\nobreak\
yr^{-1}}$) (Leroy et al. 2021). On the other hand, dust emission at longer
wavelengths of $\lambda\geq 20$ $\rm\mu m$ in the simulation is consistent
with that in observations while the PAH emission at MIR wavelengths in the
simulation is lower than that in observations by about 0.5 dex. This result
implies that our dust model underestimates the PAH mass, as reported in
Hirashita et al. (2020) and Hirashita (2023). To explore where the PAH mass
fraction is underestimated in our simulation, we perform the pixel-based
comparison of the PAH mass fraction in the next section.

###  5.2 Spatial comparison of PAH mass fraction

We compare the spatially resolved data of $\nu
I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ between our
simulation and Herschel/Spitzer observations. The observations cover large
scales up to the galactic radius of $r\sim 10$ kpc, but the spatial resolution
is about 600 pc. Figure 10(a) shows the comparison of the relation between
$\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ and the
intensity at $\mathrm{160\,\rm\mu m}$. The line and shaded area correspond to
the median value and the range between the 10th and 90th percentiles,
respectively. Here, the $\mathrm{160\,\rm\mu m}$ band traces large grains in
radiative equilibrium and represents the dust mass distribution We find that
$\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ from the
simulation is lower than that of the observation in the entire range of
$\mathrm{160\,\rm\mu m}$. This indicates an underestimation of the PAH mass
fraction in our model, regardless of location. To estimate how much additional
PAH mass is required in the NGC 628-like galaxy simulation (hereafter,
referred to as fiducial model), we assume eight times the PAH mass of the
fiducial model and find $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$ can be consistent with that of the observation.
Therefore, we conclude that we underestimate the PAH mass abundance by a
factor of eight on average in our simulated area.

Figure 10: Pixel-based comparison of the indicators of the PAH mass fraction
between the NGC 628-like galaxy simulation and actual observations. Panel (a):
$\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ versus the
$\mathrm{160\,\rm\mu m}$ intensity for the NGC 628-like galaxy simulation at
$t=10$ Gyr (orange) and based on the Herschel and Spitzer observations for NGC
628 (blue). Panel (b): $R_{\mathrm{PAH}}$ versus the $\mathrm{21\,\rm\mu m}$
intensity for the NGC 628-like galaxy simulation at $t=10$ Gyr (orange) and
based on the JWST observations for NGC 628 (blue). In both panels, the green
color represents the model, in which the PAH mass is enhanced by a factor of
eight. The solid and dashed lines represent the mean value over all pixels of
the photometric data of the models and observations, and the shaded area
corresponds to the 10th–90th percentile range, respectively.  Figure 11:
Comparison of $21$ $\rm\mu m$ intensity and $R_{\mathrm{PAH}}$ maps between
the NGC 628-like galaxy simulation and JWST observations. Panels (a) and (b):
$21$ $\rm\mu m$ intensity maps of the NGC 628-like galaxy simulation at $10$
Gyr on 100 and 20 kpc scales, respectively. Panels (c) and (d):
$R_{\mathrm{PAH}}$ maps of the simulation on 100 and 20 kpc scales,
respectively. Panels (e) and (f): maps of the $21$ $\rm\mu m$ intensity and
$R_{\mathrm{PAH}}$ based on the JWST/MIRI observations for NGC 628 (Lee et al.
2023), where the spatial resolutions are reduced to 50 pc.

Furthermore, we compare the spatially resolved data from our simulation and
the JWST observations (Lee et al. 2023), which has 12 times higher spatial
resolution than the Herschel observations (a pixel size $\sim$ 50 pc). Here,
instead of $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$, we use the intensity ratio of
$R_{\mathrm{PAH}}$ (Eq. (1)). Fig. 10(b) shows the comparison of
$R_{\mathrm{PAH}}$ to the intensity at $\mathrm{21\,\rm\mu m}$ between our
simulation (orange) and the observation (blue). We find that, in both
simulation and observation, $R_{\mathrm{PAH}}$ decreases as the
$\mathrm{21\,\rm\mu m}$ intensity increases. This trend can be explained by
our model in which coagulation is more efficient in denser regions with higher
$\mathrm{21\,\rm\mu m}$ intensities. Figures 11(a) and (b) show the intensity
maps at $21$ $\rm\mu m$ of the entire NGC 628-like galaxy at $t=10$ Gyr and
the central 6 kpc regions of the galaxy, respectively. Figs. 11(c) and (d)
show the $R_{\mathrm{PAH}}$ maps of the entire NGC 628-like galaxy at $t=10$
Gyr and the central 6 kpc regions of the galaxy, respectively. In Fig. 11(c),
$R_{\mathrm{PAH}}$ gradually increases toward the outer regions of the galaxy,
since shattering is more dominant in the diffuse ISM than coagulation. In the
central region of the galaxy (Fig. 11(d)), $R_{\mathrm{PAH}}$ is lower in the
spiral arms and clumps, where the intensity at $21$ $\rm\mu m$ is higher. This
is because coagulation is dominant in the dense ISM. These trends are
consistent with the results revealed by recent JWST observations. Figs. 11(e)
and (f) show the observed maps of $21$ $\rm\mu m$ intensity and
$R_{\mathrm{PAH}}$ maps, respectively (Lee et al. 2023). These observations
also indicate that $R_{\mathrm{PAH}}$ is lower in the spiral arms and clumps,
where the $21$ $\rm\mu m$ intensity is higher. Alternatively, it is also
possible that PAHs are destroyed by UV radiation from young stars in star-
forming regions in the dense ISM, as discussed in the literature (e.g., Egorov
et al. 2023; Chastenet et al. 2023b) since $\mathrm{21\,\rm\mu m}$ intensity
is expected to increase as UV radiation field becomes stronger. Thus, to
discriminate these scenarios, we have to implement PAH destruction due to UV
radiation from young stars in future simulations (e.g., Nanni et al. 2023). We
also find that the value of $R_{\mathrm{PAH}}$ from the simulation is lower
than that from the observation, as discussed above: $R_{\mathrm{PAH}}$ of the
model with the PAH mass enhanced by a factor of eight is consistent with that
of the observations, implying that further enhancement of small dust grains in
our model is required.

###  5.3 Insights into additional PAH formation mechanisms

As we discussed above, our dust evolution model underestimates the PAH mass
fraction in the entire galaxy, and thus, additional mechanisms to enhance the
PAH production in our model are required. In our model, we assume that the
main driver of PAH formation is the shattering of large carbonaceous grains,
but PAH formation in the envelopes of carbon-rich AGB stars may also be
considered. Furthermore, in our model, PAHs might be too destroyed by
coagulation in dense ISM. Here, we discuss the possible scenarios to enhance
PAHs in the ISM of our simulations.

Cherchneff et al. (1992) suggest that the stellar envelopes of carbon-rich AGB
stars can be a site of the PAH formation. Indeed, many observational studies
reveal evidence of the PAH formation in the stellar envelopes of AGB stars by
identifying the PAH features in the spectra of the individual stars in the
Magellanic clouds and our Galaxy (Sloan et al. 2007, 2014; Matsuura et al.
2014). However, some studies based on the Magellanic clouds report that the
presence of PAHs in the ISM can not be explained entirely only by carbon-rich
AGB stars (Matsuura et al. 2013; Sandstrom et al. 2010; Chastenet et al.
2019). Moreover, since NGC 628 has higher metallicity than the Magellanic
Clouds, stellar dust production containing PAHs should be even less efficient
than dust growth by accretion in the ISM in NGC 628. Thus, it is not probable
that stellar PAH production is dominant over PAH formation in the ISM in NGC
628.

In our model, coagulation in the dense ISM efficiently depletes PAHs by
converting small grains to larger grains. Hirashita (2023) model a process
that prevents small carbonaceous grains from being lost through coagulation by
hypothesizing that small carbonaceous grains are not involved in interstellar
processing. They succeed in enhancing small grains, especially PAHs, and
reproduce the observed global dust SED with the PAH emission at the MIR
wavelengths. Although the chemical or physical processes underlying their
hypothesis have yet to be justified, their success indicates that a key to the
PAH deficiency problem in our model lies in the treatment of interstellar
processing of PAHs.

##  6 Conclusion

In this paper, we aim to investigate how dust evolution is reflected on IR
spectral features and provide observational signatures of the dust size
evolution in ISM, in particular exploring indicators of PAH mass fraction. We
model the dust evolution for two simulations representing star-forming
galaxies similar to the Milky Way and NGC 628 using GADGET4-OSAKA (Romano et
al. 2022a). We simulate the evolution of grain size distribution in a manner
consistent with the physical conditions of the ISM (see Section 2.1) by taking
into account interstellar dust processes including stellar dust production, SN
destruction, shattering, coagulation, and metal accretion. Furthermore, we
separate the dust compositions into silicate and carbonaceous dust grains, and
PAHs in post-processing. Then, we perform radiative transfer calculations with
SKIRT under the dust optical properties consistent with the grain size
distribution and grain composition of the hydrodynamic simulations (see
Section 2.2).

We summarize the results as follows:

  1. 1.

We show the history of dust evolution and star formation for the Milky Way-
like and NGC 628-like galaxy simulations (see Fig. 1). The Milky Way-like
galaxy shows more rapid dust formation because of the earlier starburst and
metal production than the NGC 628-like galaxy. However, the Milky Way-like
galaxy shows less PAH formation in the ISM since the Milky Way-like galaxy
contains a large amount of dense gas, where shattering is inefficient. These
results suggest that different ISM conditions and SFHs cause different dust
evolution and PAH mass fractions among galaxies.

  2. 2.

We investigate the detailed time evolution of the observable properties for
the Milky Way-like simulation. In the early phase, only the emission of large
dust grains can be observed around $\lambda=70$ $\rm\mu m$ in the SED, as dust
production is dominated by stellar sources (see Fig. 3). After that, emissions
from small dust grains and PAHs become noticeable around $\lambda=10$ $\rm\mu
m$, since shattering acts as a source of small grains in diffuse regions, and
accretion also increases the abundance of small grains (see panels of $t=0.5$
Gyr in Fig. 4). Finally, the functional shape of the grain size distribution
at each region of the ISM converges to the one determined by the balance
between shattering and coagulation rate according to the condition of the ISM.
Consequently, the emission from large and small grains is relatively prominent
in the inner and outer regions of the galaxy, respectively (see panels of
$t=10$ Gyr in Fig. 4).

  3. 3.

We find that $I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu m})$ correlates with the PAH mass
fraction, but the local radiation field influences it (see Fig. 5). The 8
$\rm\mu m$ band is almost purely dominated by PAH band emission, while the 24
$\rm\mu m$ band is contaminated by the emission of large grains, resulting in
$I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu
m})$ increasing with the PAH mass fraction. However, stronger radiation fields
cause higher dust temperatures, shifting the large grain emission toward
shorter wavelengths and raising the 24 $\rm\mu m$ emission. consequently,
$I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu
m})$ decrease in regions with strong radiation. Hence,
$I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu
m})$ can be used as an indicator of the PAH mass fraction with a tendency of
underestimating PAH mass fraction in regions with strong radiation.

  4. 4.

We find that $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$
tightly correlates with the PAH mass fraction and is not affected by the local
radiation field (see Fig. 7 and Eq. 9). This is because both total and 8
$\rm\mu m$ intensities are proportional to the strength of the local radiation
field, and hence, the effect of the radiation field is canceled out.
Therefore, we suggest that $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$ is a good indicator of the PAH mass fraction in various
stages of galaxy evolution.

  5. 5.

We find that the global PAH mass fraction or $\nu
I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ in the entire galaxy
evolves with metallicity (see Fig 8). This is because the interplay between
shattering and accretion produces PAHs in the entire galaxy. On the contrary,
based on the pixel-based analysis, local $\nu
I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ decrease with
metallicity in the ISM at higher metallicity ($Z>0.2\ Z_{\odot}$), while it
increases at lower metallicity ($Z\leq 0.2\ Z_{\odot}$). At low metallicity,
the metal enrichment causes more efficient shattering, which tends to increase
PAHs. At high metallicity, coagulation is rather important in depleting PAHs.
We also find that $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu
m})/I_{\mathrm{TIR}}$ depends on the gas surface density since the coagulation
converts PAHs to large grains in dense ISM. Therefore, the spatially resolved
relation of the PAH mass fraction helps us interpret the efficiency of
accretion, coagulation, and shattering in various ISM phases.

  6. 6.

We compare the total SED in the NGC 628-like simulation with the observed one
(see Fig. 9). Although dust emission at longer wavelengths of $\lambda\geq 20$
$\rm\mu m$ in the simulation is consistent with that in the observations, the
PAH emission in the MIR wavelengths in the simulation is lower than that in
observations by about 0.5 dex. This result implies that our dust model
underestimates PAH masses.

  7. 7.

We compare the spatially resolved data of $\nu
I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$ in the NGC 628-like
galaxy simulation and the Herschel/Spitzer data (see Fig. 10(a)). We find that
the values of $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$
from the simulation are lower than those from the observations, implying that
eight times more PAH mass is required in our simulation.

  8. 8.

We also compare the spatially resolved data from our simulation to the JWST
observations with high spatial resolutions (see Fig. 10(b)). We find that
$R_{\mathrm{PAH}}$, a PAH indicator composed of JWST bands (Eq. 1), decreases
with the $\mathrm{21\,\rm\mu m}$ intensity in both the simulation and
observations. This trend can be explained by the fact that coagulation is more
efficient in denser ISM. This comparison also shows that the values of
$R_{\mathrm{PAH}}$ from the simulation are lower than those from the
observation. In our model, coagulation may deplete PAHs too efficiently in
dense ISM. Therefore, it may be required to prevent small carbonaceous grains
from being lost through coagulation in the ISM.

These results suggest that the spatial comparison between observations and
theoretical models using observable signatures of the dust size evolution
(e.g., $\nu I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{TIR}}$,
$I_{\mathrm{\nu}}(\mathrm{8\,\rm\mu m})/I_{\mathrm{\nu}}(\mathrm{24\,\rm\mu
m})$, $R_{\mathrm{PAH}}$) leads us to a quantitative investigation of dust
evolution in various ISM while maintaining the original spatial resolution of
the observational data. In particular, the comparison with future JWST
observations for many galaxies with high spatial resolutions will reveal the
efficiency of dust processes in the diffuse and dense ISM in various galaxies.

## Acknowledgements

KM is a Ph.D. fellow of the Flemish Fund for Scientific Research (FWO-
Vlaanderen) and acknowledges the financial support provided through Grant
number 1169822N. Numerical computations were carried out on the Cray XC50 at
the Center for Computational Astrophysics, National Astronomical Observatory
of Japan, and the SQUID at the Cybermedia Center, Osaka University as part of
the HPCI system Research Project (hp220044,hp230089). This work is supported
in part by the MEXT/JSPS KAKENHI grant numbers 20H00180 and 22K21349 (K.N.),
and 21J20930, 22KJ2072 (YO). K.N. acknowledges the support from the Kavli
IPMU, World Premier Research Center Initiative (WPI), where a part of this
work was conducted. HH thanks the National Science and Technology Council for
support through grant 111-2112-M-001-038-MY3 and the Academia Sinica for
Investigator Award AS-IA-109-M02. MR acknowledge support from project
PID2020-114414GB-100, financed by MCIN/AEI/10.13039/501100011033. SvdG and IDL
acknowledge financial support by the European Research Council (ERC) under de
European Union’s Horizon 2020 research and innovation programme DustOrigin
(ERC-2019-StG-851622) and by the Flemish Fund for Scientific Research (FWO-
Vlaanderen) through the research project G023821N. We thank Prof. Nakagawa for
many fruitful discussions.

## References

  * Aniano et al. (2011) Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218
  * Aniano et al. (2020) Aniano, G., Draine, B. T., Hunt, L. K., et al. 2020, ApJ, 889, 150
  * Aniyan et al. (2018) Aniyan, S., Freeman, K. C., Arnaboldi, M., et al. 2018, MNRAS, 476, 1909
  * Aoyama et al. (2020) Aoyama, S., Hirashita, H., & Nagamine, K. 2020, MNRAS, 491, 3844
  * Aoyama et al. (2017) Aoyama, S., Hou, K.-C., Shimizu, I., et al. 2017, MNRAS, 466, 105
  * Asano et al. (2013) Asano, R. S., Takeuchi, T. T., Hirashita, H., & Nozawa, T. 2013, MNRAS, 432, 637
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Baes & Camps (2015) Baes, M. & Camps, P. 2015, Astronomy and Computing, 12, 33
  * Baes et al. (2024) Baes, M., Gebek, A., Trcka, A., et al. 2024, arXiv e-prints, arXiv:2401.04224
  * Baes et al. (2011) Baes, M., Verstappen, J., De Looze, I., et al. 2011, ApJS, 196, 22
  * Berg et al. (2015) Berg, D. A., Skillman, E. D., Croxall, K. V., et al. 2015, ApJ, 806, 16
  * Boquien et al. (2011) Boquien, M., Calzetti, D., Combes, F., et al. 2011, AJ, 142, 111
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Burgarella et al. (2020) Burgarella, D., Nanni, A., Hirashita, H., et al. 2020, A&A, 637, A32
  * Camps & Baes (2015) Camps, P. & Baes, M. 2015, Astronomy and Computing, 9, 20
  * Camps & Baes (2020) Camps, P. & Baes, M. 2020, Astronomy and Computing, 31, 100381
  * Camps et al. (2013) Camps, P., Baes, M., & Saftly, W. 2013, A&A, 560, A35
  * Camps et al. (2021) Camps, P., Behrens, C., Baes, M., Kapoor, A. U., & Grand, R. 2021, ApJ, 916, 39
  * Camps et al. (2022) Camps, P., Kapoor, A. U., Trcka, A., et al. 2022, MNRAS, 512, 2728
  * Camps et al. (2015) Camps, P., Misselt, K., Bianchi, S., et al. 2015, A&A, 580, A87
  * Camps et al. (2018) Camps, P., Trčka, A., Trayford, J., et al. 2018, ApJS, 234, 20
  * Cazaux & Tielens (2004) Cazaux, S. & Tielens, A. G. G. M. 2004, ApJ, 604, 222
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chang et al. (2022) Chang, C.-Y., Huang, Y.-H., Hirashita, H., & Cooper, A. P. 2022, MNRAS, 513, 2158
  * Chastenet et al. (2019) Chastenet, J., Sandstrom, K., Chiang, I.-D., et al. 2019, ApJ, 876, 62
  * Chastenet et al. (2023a) Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023a, ApJ, 944, L11
  * Chastenet et al. (2023b) Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023b, ApJ, 944, L12
  * Cherchneff et al. (1992) Cherchneff, I., Barker, J. R., & Tielens, A. G. G. M. 1992, ApJ, 401, 269
  * Chiang et al. (2018) Chiang, I.-D., Sandstrom, K. M., Chastenet, J., et al. 2018, ApJ, 865, 117
  * Clark et al. (2018) Clark, C. J. R., Verstocken, S., Bianchi, S., et al. 2018, A&A, 609, A37
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393
  * da Cunha et al. (2008) da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595
  * Dale et al. (2009) Dale, D. A., Cohen, S. A., Johnson, L. C., et al. 2009, ApJ, 703, 517
  * Davies et al. (2017) Davies, J. I., Baes, M., Bianchi, S., et al. 2017, PASP, 129, 044102
  * Draine & Anderson (1985) Draine, B. T. & Anderson, N. 1985, ApJ, 292, 494
  * Draine et al. (2007) Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866
  * Draine & Li (2001) Draine, B. T. & Li, A. 2001, ApJ, 551, 807
  * Draine & Li (2007) Draine, B. T. & Li, A. 2007, ApJ, 657, 810
  * Draine et al. (2021) Draine, B. T., Li, A., Hensley, B. S., et al. 2021, ApJ, 917, 3
  * Dwek (1998) Dwek, E. 1998, ApJ, 501, 643
  * Egorov et al. (2023) Egorov, O. V., Kreckel, K., Sandstrom, K. M., et al. 2023, ApJ, 944, L16
  * Einasto (1965) Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Galametz et al. (2013) Galametz, M., Kennicutt, R. C., Calzetti, D., et al. 2013, MNRAS, 431, 1956
  * Galliano et al. (2008a) Galliano, F., Dwek, E., & Chanial, P. 2008a, ApJ, 672, 214
  * Galliano et al. (2008b) Galliano, F., Madden, S. C., Tielens, A. G. G. M., Peeters, E., & Jones, A. P. 2008b, ApJ, 679, 310
  * Galliano et al. (2021) Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18
  * Gebek et al. (2023) Gebek, A., Baes, M., Diemer, B., et al. 2023, MNRAS, 521, 5645
  * Granato et al. (2021) Granato, G. L., Ragone-Figueroa, C., Taverna, A., et al. 2021, MNRAS, 503, 511
  * Graziani et al. (2020) Graziani, L., Schneider, R., Ginolfi, M., et al. 2020, MNRAS, 494, 1071
  * Grieco et al. (2023) Grieco, F., Orsolini, Y., & Pérot, K. 2023, Journal of Atmospheric and Solar-Terrestrial Physics, 252, 106147
  * Haardt & Madau (2012) Haardt, F. & Madau, P. 2012, ApJ, 746, 125
  * Hensley & Draine (2023) Hensley, B. S. & Draine, B. T. 2023, ApJ, 948, 55
  * Hirashita (1999) Hirashita, H. 1999, ApJ, 510, L99
  * Hirashita (2012) Hirashita, H. 2012, MNRAS, 422, 1263
  * Hirashita (2015) Hirashita, H. 2015, MNRAS, 447, 2937
  * Hirashita (2023) Hirashita, H. 2023, MNRAS, 518, 3827
  * Hirashita & Aoyama (2019) Hirashita, H. & Aoyama, S. 2019, MNRAS, 482, 2555
  * Hirashita et al. (2020) Hirashita, H., Deng, W., & Murga, M. S. 2020, MNRAS, 499, 3046
  * Hirashita & Kuo (2011) Hirashita, H. & Kuo, T.-M. 2011, MNRAS, 416, 1340
  * Hirashita & Murga (2020) Hirashita, H. & Murga, M. S. 2020, MNRAS, 492, 3779
  * Hollenbach & McKee (1979) Hollenbach, D. & McKee, C. F. 1979, ApJS, 41, 555
  * Hollenbach & Salpeter (1971) Hollenbach, D. & Salpeter, E. E. 1971, ApJ, 163, 155
  * Hou et al. (2019) Hou, K.-C., Aoyama, S., Hirashita, H., Nagamine, K., & Shimizu, I. 2019, MNRAS, 485, 1727
  * Kapoor et al. (2023) Kapoor, A. U., Baes, M., van der Wel, A., et al. 2023, MNRAS, 526, 3871
  * Kapoor et al. (2021) Kapoor, A. U., Camps, P., Baes, M., et al. 2021, MNRAS, 506, 5703
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541
  * Kennicutt et al. (2011) Kennicutt, R. C., Calzetti, D., Aniano, G., et al. 2011, PASP, 123, 1347
  * Kim et al. (2016) Kim, J.-h., Agertz, O., Teyssier, R., et al. 2016, ApJ, 833, 202
  * Lee et al. (2023) Lee, J. C., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJ, 944, L17
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170
  * Leroy et al. (2021) Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, ApJS, 257, 43
  * Li (2020) Li, A. 2020, Nature Astronomy, 4, 339
  * Li & Draine (2001) Li, A. & Draine, B. T. 2001, ApJ, 554, 778
  * Lisenfeld & Ferrara (1998) Lisenfeld, U. & Ferrara, A. 1998, ApJ, 496, 145
  * Maloney et al. (1996) Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 466, 561
  * Matsumoto et al. (2023) Matsumoto, K., Camps, P., Baes, M., et al. 2023, A&A, 678, A175
  * Matsuura et al. (2014) Matsuura, M., Bernard-Salas, J., Lloyd Evans, T., et al. 2014, MNRAS, 439, 1472
  * Matsuura et al. (2013) Matsuura, M., Woods, P. M., & Owen, P. J. 2013, MNRAS, 429, 2527
  * McKinnon et al. (2018) McKinnon, R., Vogelsberger, M., Torrey, P., Marinacci, F., & Kannan, R. 2018, MNRAS, 478, 2851
  * Nanni et al. (2023) Nanni, A., Cristallo, S., Donevski, D., et al. 2023, arXiv e-prints, arXiv:2312.15998
  * Narayanan et al. (2023) Narayanan, D., Smith, J. D., Hensley, B., et al. 2023, arXiv e-prints, arXiv:2301.07136
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Nersesian et al. (2019) Nersesian, A., Xilouris, E. M., Bianchi, S., et al. 2019, A&A, 624, A80
  * Nishida et al. (2022) Nishida, K. Y., Takeuchi, T. T., Nagata, T., & Asano, R. S. 2022, MNRAS, 514, 2098
  * O’Donnell & Mathis (1997) O’Donnell, J. E. & Mathis, J. S. 1997, ApJ, 479, 806
  * Ostriker & Silk (1973) Ostriker, J. & Silk, J. 1973, ApJ, 184, L113
  * Perret (2016) Perret, V. 2016, DICE: Disk Initial Conditions Environment, Astrophysics Source Code Library, record ascl:1607.002
  * Rau et al. (2019) Rau, S.-J., Hirashita, H., & Murga, M. 2019, MNRAS, 489, 5218
  * Relaño et al. (2022) Relaño, M., De Looze, I., Saintonge, A., et al. 2022, MNRAS, 515, 5306
  * Relaño et al. (2020) Relaño, M., Lisenfeld, U., Hou, K. C., et al. 2020, A&A, 636, A18
  * Rémy-Ruyer et al. (2015) Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2015, A&A, 582, A121
  * Romano et al. (2022a) Romano, L. E. C., Nagamine, K., & Hirashita, H. 2022a, MNRAS, 514, 1441
  * Romano et al. (2022b) Romano, L. E. C., Nagamine, K., & Hirashita, H. 2022b, MNRAS, 514, 1461
  * Saftly et al. (2014) Saftly, W., Baes, M., & Camps, P. 2014, A&A, 561, A77
  * Saftly et al. (2013) Saftly, W., Camps, P., Baes, M., et al. 2013, A&A, 554, A10
  * Saitoh (2017) Saitoh, T. R. 2017, AJ, 153, 85
  * Sandstrom et al. (2010) Sandstrom, K. M., Bolatto, A. D., Draine, B. T., Bot, C., & Stanimirović, S. 2010, ApJ, 715, 701
  * Schreiber et al. (2018) Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30
  * Seok et al. (2014) Seok, J. Y., Hirashita, H., & Asano, R. S. 2014, MNRAS, 439, 2186
  * Shim et al. (2023) Shim, H., Hwang, H. S., Jeong, W.-S., et al. 2023, AJ, 165, 31
  * Shimizu et al. (2019) Shimizu, I., Todoroki, K., Yajima, H., & Nagamine, K. 2019, MNRAS, 484, 2632
  * Sloan et al. (2007) Sloan, G. C., Jura, M., Duley, W. W., et al. 2007, ApJ, 664, 1144
  * Sloan et al. (2008) Sloan, G. C., Kraemer, K. E., Wood, P. R., et al. 2008, ApJ, 686, 1056
  * Sloan et al. (2014) Sloan, G. C., Lagadec, E., Zijlstra, A. A., et al. 2014, ApJ, 791, 28
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217
  * Springel et al. (2021) Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871
  * Tielens & Hollenbach (1985) Tielens, A. G. G. M. & Hollenbach, D. 1985, ApJ, 291, 722
  * Trčka et al. (2022) Trčka, A., Baes, M., Camps, P., et al. 2022, MNRAS, 516, 3728
  * Vander Meulen et al. (2023) Vander Meulen, B., Camps, P., Stalevski, M., & Baes, M. 2023, A&A, 674, A123
  * Vogelsberger et al. (2019) Vogelsberger, M., McKinnon, R., O’Neil, S., et al. 2019, MNRAS, 487, 4870
  * Weingartner & Draine (2001) Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296
  * Wolfire et al. (2022) Wolfire, M. G., Vallini, L., & Chevance, M. 2022, ARA&A, 60, 247

##  Appendix A Initial conditions of the NGC 628-like galaxy

We construct the initial condition for the NGC 628-like simulation using the
DICE code (Perret 2016), which computes the gravitational potential on a
multi-level Cartesian mesh and considers the dynamical equilibrium of multiple
material components. We build up the structural properties of the five
components: the dark matter halo, stellar disk and bulge, and gaseous disk and
halo. First, the total halo is characterized by a viral mass of $1.2\times
10^{12}$ ${\rm M_{\odot}}$, a virial velocity of $v_{200}$=150
$\mathrm{km\leavevmode\nobreak\ s^{-1}}$, and spin parameter $\lambda=0.04$.
These parameters are determined so that the velocity curve of the simulation
is consistent with that of the NGC 628 observations (Aniyan et al. 2018). The
dark matter halo with a mass of $1.0\times 10^{12}$ ${\rm M_{\odot}}$ follows
the Navarro–Frenk–White (Navarro et al. 1997) profile with the concentration
parameter of $c=25$ and the effective radius of 100 kpc. The gaseous halo
follows the same profile but has 1 % of the total mass, $1.2\times 10^{10}$
${\rm M_{\odot}}$. Second, we obtain the morphological parameters of the
gaseous and stellar disk by fitting the observationally estimated radial
profiles of the gas and stellar surface density of NGC 628 (S. van der Giessen
et al. in prep.) with the Sersic profile given by,

$\Sigma(r)=\Sigma_{0}\exp{-b_{n}\Bigg{[}\Big{(}\frac{r}{R_{e}}\Big{)}^{1/n}-1\Bigg{]}},$
(10)

where $r$, $R_{e}$, $\Sigma_{0}$, $n$, and $b_{n}$ are the radius in the
cylindrical coordinate, the half-mass radius, the surface density at the half-
light radius, Sérsic index, and a coefficient, respectively. The parameters
obtained in the fitting are summarized in Table 2. Furthermore, we assume the
vertical density profile of the galaxy to be the exponential profile defined
by

$\rho(z)\propto\exp(-z/h_{z}),$ (11)

where $z$ and $h_{z}$ are the height in the cylindrical coordinate and the
scale height, respectively. The scale height was obtained from Aniyan et al.
(2018). The gaseous and stellar disks have masses of $2.6\times 10^{10}$ and
$3.5\times 10^{10}$ ${\rm M_{\odot}}$, respectively. Finally, the stellar
bulge follows the Einasto profile (Einasto 1965) with a structural parameter
of $\alpha=1.0$ and a mass of $9.26\times 10^{9}$ ${\rm M_{\odot}}$.

Table 2: Morphological parameters of the gas and stellar components in NGC 628-like galaxy simulations Type |  $\Sigma_{0}$ $(\mathrm{{\rm M_{\odot}}/pc^{-2}})$ |  $R_{e}$ $(\mathrm{kpc})$ | $n$ | $b$ |  $h_{z}$ $(\mathrm{kpc})$
---|---|---|---|---|---
Star | 12 | 12 | 1.5 | 3.9 | 0.4
Gas | 15 | 6.5 | 0.5 | 0.26 | 0.4

##  Appendix B Various maps of the Milky Way-like galaxy

Figure 12 shows the maps of the hydrogen surface density, metallicity, local
radiation field, and PAH mass fraction, $q_{\mathrm{PAH}}$, (left to right
columns) of the Milky Way-like galaxy simulation at $t=0.5,\ 1.0,\ 3.0,$ and
$10$ Gyr (top to bottom panels). The metallicity and local radiation field are
normalized by the Milky Way values. We only show the pixels that satisfy
$I_{\nu}>0.01$ $\mathrm{MJy\,sr^{-1}}$ at all the Spitzer and Herschel bands;
that is, the pixels used for the analysis (Section 2.2). The maps of the
hydrogen surface density show spiral arms and clumps inside the arms, which
reflect the resulting complex structures in the ISM containing diffuse and
dense regions. In contrast, the metallicity maps show a smoother axisymmetric
gradient in the radial direction. The metallicity reaches the solar one at
$t=10$ Gyr in the entire disk. In the maps of the local radiation field, the
star-formation regions associated with stellar clusters along the spiral arms
are prominent (redder colors). Given the monotonic decrease in the SFR of the
Milky Way-like galaxy, the star-forming regions shrink with time in terms of
spatial distributions and sizes. The PAH mass fraction, $q_{\mathrm{PAH}}$, in
the entire galaxy increases with time. PAHs are relatively deficient in dense
regions, especially in spiral arms and the galaxy center, because of
inefficient PAH formation by shattering and efficient PAH depletion by
coagulation.

Figure 12: Maps of the hydrogen surface density, metallicity, local radiation
field, and PAH mass fraction (left to right) of the Milky Way-like galaxy
simulation at $t=0.5,\ 1.0,\ 3.0,$ and $10$ Gyr (top to bottom).


* [23] S. Kim and K.-A. Lee.  Parabolic Harnack inequality of viscosity solutions on Riemannian manifolds.  J. Funct. Anal., 267(7):2152–2198, 2014.
  * [24] Y.-C. Kim and K.-A. Lee.  Regularity results for fully nonlinear parabolic integro-differential operators.  Math. Ann., 357(4):1541–1576, 2013.
  * [25] Y.-C. Kim and K.-A. Lee.  Cordes-Nirenberg type estimates for nonlocal parabolic equations.  Nonlinear Anal., 130:330–360, 2016.
  * [26] Y.-C. Kim and K.-A. Lee.  The Evans-Krylov theorem for nonlocal parabolic fully nonlinear equations.  Nonlinear Anal., 160:79–107, 2017.
  * [27] A. Laforgia.  Bounds for modified Bessel functions.  J. Comput. Appl. Math., 34(3):263–267, 1991.
  * [28] X.-N. Ma, N. S. Trudinger, and X.-J. Wang.  Regularity of potential functions of the optimal transportation problem.  Arch. Ration. Mech. Anal., 177(2):151–183, 2005.
  * [29] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.  NIST handbook of mathematical functions.  U.S. Department of Commerce, National Institute of Standards and Technology, Washington, DC; Cambridge University Press, Cambridge, 2010.  With 1 CD-ROM (Windows, Macintosh and UNIX).
  * [30] J. G. Ratcliffe.  Foundations of hyperbolic manifolds, volume 149 of Graduate Texts in Mathematics.  Springer, New York, second edition, 2006.
  * [31] L. Saloff-Coste.  Uniformly elliptic operators on Riemannian manifolds.  J. Differential Geom., 36(2):417–450, 1992.
  * [32] J. Segura.  Bounds for ratios of modified Bessel functions and associated Turán-type inequalities.  J. Math. Anal. Appl., 374(2):516–528, 2011.
  * [33] A. Terras.  Harmonic analysis on symmetric spaces and applications. I.  Springer-Verlag, New York, 1985.
  * [34] W. P. Thurston.  Three-dimensional geometry and topology. Vol. 1, volume 35 of Princeton Mathematical Series.  Princeton University Press, Princeton, NJ, 1997.  Edited by Silvio Levy.
  * [35] A. A. Ungar.  Thomas rotation and the parametrization of the Lorentz transformation group.  Found. Phys. Lett., 1(1):57–89, 1988.
  * [36] A. A. Ungar.  Thomas precession: its underlying gyrogroup axioms and their use in hyperbolic geometry and relativistic physics.  Found. Phys., 27(6):881–951, 1997.
  * [37] A. A. Ungar.  Analytic hyperbolic geometry.  World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2005.  Mathematical foundations and applications.
  * [38] A. A. Ungar.  Analytic hyperbolic geometry and Albert Einstein’s special theory of relativity.  World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2008.
  * [39] A. A. Ungar.  From Möbius to gyrogroups.  Amer. Math. Monthly, 115(2):138–144, 2008.
  * [40] C. Villani.  Optimal transport, volume 338 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, 2009.  Old and new.
  * [41] H. Wang, X. Bian, and X. Bian.  Möbius transformation and a version of Schwarz lemma in octonionic analysis.  Math. Methods Appl. Sci., 44(1):27–42, 2021.
  * [42] Y. Wang and X. Zhang.  An Alexandroff-Bakelman-Pucci estimate on Riemannian manifolds.  Adv. Math., 232:499–512, 2013.
  * [43] S. T. Yau.  Harmonic functions on complete Riemannian manifolds.  Comm. Pure Appl. Math., 28:201–228, 1975.