© The Authors 2024

1 Introduction

Dust in the interstellar medium (ISM) plays a key role in shaping the observed properties of galaxies, in addition to indirectly affecting galaxy evolution. Dust governs various radiative processes in galaxies: dust absorbs and scatters ultraviolet (UV) and optical light from young stars, then reemiting the absorbed energy in the infrared (IR). These processes are reflected in the spectral energy distributions (SEDs) of galaxies; thus, analyses of the UV-to-IR SEDs of galaxies provide various clues to their physical properties, such as the 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; also, the 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 molecularly rich environment is favorable to the progress of star formation, with dust acting as a cooling channel in molecular clouds (Ostriker & Silk 1973; Vogelsberger et al. 2019). On the contrary, photoelectric heating by dust grains is the main heating source of the photodissociation and X-ray dissociation regions in the ISM, thereby altering the ISM conditions and leading to chemical reactions (Tielens & Hollenbach 1985; Maloney et al. 1996; Wolfire et al. 2022). Therefore, dust certainly stands as 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 Spitɀer telescopes, which cover wavelengths of 3–500 μm (e.g., da Cunha et al. 2008; Leja et al. 2017; Nersesian et al. 2019; Burgarella et al. 2020). Observed 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 in ~3–20 μ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) revealed 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) suggested that PAHs are produced by the shattering of carbonaceous grains, which form more efficiently via the accretion of metals on dust grains (this process is simply referred to as accretion) in regions with higher metallicity, naturally explaining the increasing trend of PAH mass fraction against metallicity. Chastenet et al. (2019) investigated the PAH mass fraction with spatially resolved data on the Small and Large Magellanic Clouds (SMC and LMC, respectively), finding a higher PAH mass fraction in the diffuse medium in the LMC compared to that in the SMC. To explain this result, they suggested shattering large grains to PAHs is preferable over other PAH formation scenarios (e.g., PAH formation in the envelope of 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 at higher metallic-ities. Therefore, these observations using PAH features provide us with clues to understanding the dust processes in various ISM phases.

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. Meanwhile, photometric and spectroscopic observations that trace prominent PAH emission features are useful for deriving more detailed spatial PAH properties (mass, size, and ionization) and offering a more comprehensive understanding of dust processes in the ISM (Li 2020). Today, the James Webb Space Telescope (JWST) has started to investigate the spatially resolved PAH features of several nearby galaxies. Egorov et al. (2023) suggested the use of the intensity ratio of the emission at 7.7 and 11.3 μm to the 21 μm emission (R_PAH) for an indicator of the PAH mass fraction given by $$R_{\text{PAH}}=\frac{I_v(7.7\text{\hspace{0.17em}μm})+I_v(11.3\text{\hspace{0.17em}μm})}{I_v(21\text{\hspace{0.17em}μm})}.$$(1)

They found that R_PAH decreases with the ionization fraction and implies that PAH destruction by UV radiation works in the star-forming regions. Chastenet et al. (2023b) exhibited a decreasing trend of the 3.35–11.3 μm intensity ratio against molecular gas fraction, implying that larger or colder grains are more abundant in denser regions. Although those observations imply spatially varying dust processes in the ISM, it is difficult to quantitatively infer the actual dust properties such as mass and grain sizes. Theoretically modeling the observable signatures of dust evolution would be useful to interpret observed dust and PAH emission.

In recent 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, such models calculate the grain size distributions for predictions of dust signatures, such as extinction curves and dust emission SEDs. Some models describe the evolution of grain size distributions (e.g., O’Donnell & Mathis 1997; Asano et al. 2013; Hirashita & Aoyama 2019) and have indeed been used to predict extinction curves. Hirashita et al. (2020) and Nishida et al. (2022) investigated the evolution of the dust emission SED based on the calculated grain size distribution.

Today, these evolution models of grain size distributions have been applied to hydrodynamical simulations (McKinnon et al. 2018; Aoyama et al. 2017, 2020; Hou et al. 2019; 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) solved the evolution of grain size distributions represented by 32 grain radius bins in an isolated galaxy simulation with a hydrodynamic simulation code, GADGET3-OSAKA. These authors found that the grain size distribution depends on the physical condition (especially the density and metallicity) of the 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) extended their models to more sophisticated subgrid models, in which the fraction of dense gas on a sub-resolution scale is varied according to the ambient density, where 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 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 the radiative properties are yet to be clarified based on the simulated stellar and dust distributions.

Narayanan et al. (2023) employed hydrodynamical simulations and post-processing radiative transfer calculations, which handle the grain size distribution and dust composition including PAHs, to investigate characteristics of dust emission in different types of galaxies. They found a significant influence of the radiation field on the ratio between global PAH luminosity and mass, inferring the importance of the radiative processes in deriving actual observational properties of dust. However, it has not yet been clarified how the PAH abundance can be extracted from such a strong influence of the radiation field. Therefore, it is crucial to investigate how grain size distributions are reflected on IR spectral features in simulations and infer observable signatures.

In this work, we model the evolution of grain size distributions in isolated galaxies with hydrodynamic simulations (Romano et al. 2022a) and perform radiative transfer calculations taking into account 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 has been well studied in different types of simulations as a typical star-forming galaxy useful for explaining the impact of 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 metallicity, 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 (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 prep.).

The main goals of this paper are as follows. First, we explore what drives the variations of grain size distributions and PAH abundances in different galaxy environments. Second, we investigate how dust evolution is reflected in spatially resolved dust emission and SEDs at various stages of galaxy evolution. Then, we infer and validate various indicators of PAH mass fraction through spatially resolved comparisons with the actual PAH mass fraction in hydrodynamic simulations. Finally, 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 Sect. 3, we explore the differences in the dust evolution and PAH abundance between Milky Way-like and NGC 628-like galaxies. In Sect. 4, we investigate the evolution of the spectral features in SEDs of the Milky Way-like galaxy and the indicators of PAH mass fraction using intensity ratios at various wavelengths. In Sect. 5, we perform spatially resolved comparisons for the indicators of PAH mass fraction between observations and simulations, and discuss additional physical processes required for our dust model. In Sect. 6, we give our conclusions.

2 Methods

2.1 Isolated galaxy simulations

To perform isolated galaxy simulations, we used a smooth particle hydrodynamic (SPH) simulation code, GADGET4-OSAKA (Romano et al. 2022a,b), which is a modified version of GAD-GET4 (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 and dust formation within the context of hydrodynamic evolution in galaxies. Here, we briefly summarize those models. We refer to the cited 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 ϵ_* = 0.05 that follows the Kennicutt-Schmidt law (Kennicutt 1998). This process is initiated when gas particles satisfy certain conditions, specifically when n_H > 20 cm⁻³ and T_gas < 10⁴ K are satisfied, where n_H and T_gas are the hydrogen number density and gas temperature, respectively. 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 utilized the GRACKLE-3 radiative cooling library (Smith et al. 2017)¹, 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 models used by Hirashita & Aoyama (2019) and Aoyama et al. (2020). Using this model, we calculated 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) da is the number density of grains whose radii are between a and a + da. We considered the grain radii ranging from 3.0 × 10⁻⁴ to 10 μm with 30 logarithmically spaced discrete bins. Although we assumed a single dust species in the calculation of the grain size distribution, we later decomposed the resulting grain size distribution in various species (as described in Sect. 2.2). This approach has the advantage of avoiding complexity in interstellar processing caused by inter-species collisions and we can still separately address the uncertainties by changing the adopted grain material properties (noted later in this work).

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 μm and variance of 0.47 (Asano et al. 2013). Simultaneously, dust is also destroyed by sputtering in SN shocks. Our model also takes into account a variety of interstellar dust processes including dust size growth through accretion and coagulation in the dense ISM and shattering and thermal sputtering in the diffuse ISM. The treatment of each process in the simulation is described by Aoyama et al. (2020) and Romano et al. (2022a). We refer to these papers for more details.

We summarize the treatment of each process in the following. Coagulation and shattering are assumed to occur by grain–grain collisions in the dense and diffuse ISM, respectively. Both of these processes are governed by a Smoluchowski equation (or its extended version) with a Kernel function determined by the radii of colliding grains and the grain velocities determined by the coupling with the turbulent motion, which has a Kolmogorov spectrum (Kolmogorov 1941) and a characteristic velocity reflecting the local physical condition. The total mass of fragments produced by shattering in each collision is determined following Kobayashi & Tanaka (2010), and is redistributed according to a power-law fragment size distribution with an index of −3.3 (Jones et al. 1996). For sputtering and accretion, the grain size distribution evolves according to an advection equation in the grain radius space with an appropriate rate of grain radius increase or decrease. The sputtering rate is taken from Tsai & Mathews (1995) in the hot gas, while it is regulated by the sweeping rate of SN shocks for SN destruction. The accretion rate is determined by the collision frequency between gas-phase metals and dust.

In our hydrodynamic simulations, the parameters for the interstellar dust processes are evaluated based on the properties of silicate. In reality, however, the efficiency of interstellar dust processing should vary for different dust compositions, and the accretion efficiency of dust also varies for each gas-phase metal element (Granato et al. 2021; Choban et al. 2022; Narayanan et al. 2023). Consequently, the timescale for grain size evolution will differ for each dust species. However, previous studies found that both silicate and carbonaceous dust trace the same evolutionary trends in the grain size distribution within a ~Gyr time lag (Hou et al. 2019; Hirashita & Murga 2020). Therefore, as far as the evolutionary path of grain size distribution on Gyr timescales is concerned, the treatment of the grain composition does not have a significant impact on the grain size and SED evolution over 10 Gyr time. In particular, the final grain size distribution achieved at t ~ 10 Gyr, suitable for nearby galaxies in comparison, is robustly determined by the balance between shattering and coagulation, which is not sensitive to the assumed grain composition.

We note that our simulations do not have sufficient spatial resolutions for cold and dense gas (T_gas < 50 K and n_H > 10³ cm⁻³). On the other hand, Romano et al. (2022a) employed a two-phase ISM subgrid model for gas particles, where it is assumed that dense and diffuse gas phases co-exist. The mass fraction of the dense gas phase is given by $$f_{\text{dense\hspace{0.17em}}}=\min \left(\alpha \frac{n_{\text{H}}}{1.0{\text{\hspace{0.17em}cm}}^{-3}},1.0\right),$$(2)

for which we adopted α = 0.12. Equation (2) is motivated by the fact that denser gas hosts more dense clouds with molecular hydrogen (Gnedin & Kravtsov 2011; Gnedin & Draine 2014). Romano et al. (2022a) calibrate α for the Milky Way-like galaxy to obtain a global dense gas fraction of 20%, which corresponds to a mean value of the molecular gas fraction for galaxies with stellar masses of around 10¹⁰ M_⊙ (Catinella et al. 2018). While the temperature and density of the dense gas medium were assumed to be constant (T_dense = 50 K and n_{H, dense} = 10³ cm⁻³), we estimated the temperature and density of the diffuse gas medium in each gas particle by considering the conservation of the internal energy. This subgrid model allowed 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).

Using these models, we simulated Milky Way-like and NGC 628-like galaxies. While the initial condition for the Milky Way-like galaxy was taken from the AGORA project² (Kim et al. 2016), the one for the NGC 628-like galaxy was created by DICE³ (Perret 2016). To create the initial condition of the NGC 628-like galaxy, we needed the morphological parameters of stellar and gas distributions in NGC 628. Thus, we computed those parameters based on the observationally obtained radial profiles of stars and gas (see Appendix A for details). Moreover, we set dark matter and gaseous halos, bulge stars, and stellar and gaseous disks in the initial conditions as summarized in Table 1.

We ran both simulations until the simulation time of 10 Gyr, where the metallicity in both simulations is consistent with the observed ones (Berg et al. 2015, for NGC 628). Here, we adopted an initial metallicity of 0.0001 Z_⊙, where Z_⊙ is the solar metal-licity (Z_⊙ = 0.013; Asplund et al. 2009), to start our simulations in the condition of the early universe. For the hydrodynamic simulations, we permited the minimum smoothing length to reach 0.1 ϵ_grav, where ϵ_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 employed the state-of-the-art three-dimensional (3D) Monte Carlo radiative transfer code, SKIRT⁴ (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 (e.g., Camps et al. 2021; Gebek et al. 2023; Vander Meulen et al. 2023; Matsumoto et al. 2023), here we consider absorption, scattering, emission, and self-absorption only by dust. We include stochastic heating of dust grains in SKIRT (Camps et al. 2015), which is particularly important for small grains and PAHs (Duley 1973; Purcell 1976; Draine & Anderson 1985). The code is set up to deal with dust mixtures with spatially varying grain size distributions and optical properties⁵. 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 utilized 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 needed to construct the star formation history (SFH) before t = 0 Gyr for the disk and bulge star particles. Since we ran the simulations until t = 10 Gyr and assumed 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 assigned 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_⊙.

Dust models. We applied the following dust composition model to each gas particle containing dust grains in the postprocessing. 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_silicate and f_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 $$f_{\text{silicate\hspace{0.17em}}}=\frac{6Z_{\text{Si}}}{6Z_{\text{Si}}+Z_{\text{C}}},$$(3) $$f_{\text{carbon\hspace{0.17em}}}=\left(1-f_{\text{silicate\hspace{0.17em}}}\right),$$(4)

where Z_Si and Z_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 ≤ 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 ≤ 1 .3 nm are composed of PAH⁰, 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 (PAH⁰ and PAH⁺) carbonaceous grains correspond to the diffuse gas fraction (ƒ_aroma = 1 − ƒ_dense; Hirashita & Murga 2020). In the analysis of this paper, the PAH mass fraction is given by $$q_{\text{PAH}}=\frac{M_{\text{PAH}(\text{a}\le 1.3\text{\hspace{0.17em}nm})}}{M_{\text{(total\hspace{0.17em}dust)\hspace{0.17em}}}}\times 100.$$(5)

Finally, to estimate the mass fraction for the PAH⁰ 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): $$f_{\text{ion\hspace{0.17em}}}(a)=1-\frac{1}{1+a/10Å}.$$(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 observed 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 view to (40 kpc)² and adopt a pixel size of 50 pc, which is close to the gravitational softening length of the Milky Way-like galaxy simulation. We selected pixels from each band data with an intensity (I_ν > 0.01 MJy sr⁻¹) to exclude those 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. To compare this simulation with the actual observations in Sect. 5, we used the field of view of (20 kpc)² and a pixel size of 50 pc, which corresponds to the gravitational softening length of the NGC 628-like galaxy simulation. The pixels of the data in each band with an intensity (I_ν > 0.01 MJy sr⁻¹) are selected. When we compared the simulation data with the Herschel and Spitɀer data, the resolution was reduced to 600 pc to make it consistent with the beam size of the Herschel PACS 160 μm for NGC 628 (see also Sect. 2.3).

2.3 Observational data and analysis

For the observational photometric data of NGC 628, we used the Dustpedia database (Davies et al. 2017; Clark et al. 2018)⁶, Spitɀer Local Volume Legacy survey (Dale et al. 2009), and PHANGS (Physics at High Angular resolution in Nearby GalaxieS) survey (Lee et al. 2023)⁷ for Herschel, Spitɀer, and JWST, respectively. More specifically, we used the Herschel PACS 70 and 160 μm bands, the Spitɀer IRAC 8 μm and MIPS 24 μm bands, and JWST MIRI 7.7, 11.3, and 21 μm bands. After we obtained the data, we convolved each image with the kernel function of the point spread function of each instrument. The kernel functions were taken from the Astronomical Convolution-kernel repository (Aniano et al. 2011)⁸. To accomplish the comparison between the observations and the NGC 628-like galaxy simulation in Sect. 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 Spitɀer images were unified to 600 pc, corresponding to the beam size of the Herschel PACS 160 μm. After these processes, we only adopted the pixels whose intensity is above the 1 sigma level of the background noise in the original photometric data.

Fig. 1 Time evolution of global physical properties of the Milky Waylike 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.

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

Here, 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 1a shows the SFH of the Milky Way-like and NGC 628-like galaxies. In the Milky Way-like galaxy, a large amount of gas assigned 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. The radial gas profile of the NGC 628-like galaxy is flatter than that of the Milky Way-like galaxy (see Appendix A) and the gas in the entire NGC 628-like galaxy is not dense enough to cause explosive star formation. Thus, the gas is slowly consumed by star formation with SFR ~ 1.0 M_⊙ yr⁻¹ over 10 Gyr.

The different SFHs between the two 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. 1b since more dust is produced in the starburst phase of the former 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; Draine 2009). Figure 1c 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 1d 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 increases less steeply than the NGC 628-like one, since the former contains more amounts of dense gas, and hence, shattering is less efficient. The difference in the ISM conditions between these galaxies also affects 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. These results suggest that different ISM conditions and SFHs cause different dust evolution.

Fig. 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 shades in the same color show the range of 10th–90th percentiles for the sample of gas particles.

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

To investigate observational signatures of the grain size evolution, we focus on the results from the Milky Way-like galaxy simulations. The Milky Way-like galaxy simulation from the AGORA project captures various ISM environments with a broad range of gas densities and metallicities and can track the dust evolution over long periods of 10 Gyr. 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 Sect. 4.2.

4.1 Influences of the dust size evolution on dust emission

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 λ = 70 μm, as dust production is dominated by stellar sources, which predominantly produce large grains. At t ≥ 0.5 Gyr, the emission of small dust grains and PAHs becomes noticeable around λ = 10 μ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 mono-tonic 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 λ = 70 to 150 μm over time. This shift corresponds to a decrease in the dust temperature from T_dust = 40 to 20 K. On the other hand, the dust mass monotonically increases with time mainly through the accretion of gas-phase metals. Thus, at t = 1 Gyr, the total IR luminosity reaches its maximum.

Figure 4 shows the stellar and dust emission maps from the Milky Way-like galaxy simulation at simulation times of t = 0.05, 0.5, 1.0, 3.0, and 10 Gyr. The first two rows show the stellar emission seen from face-on and edge-on with colors synthesized using the SDSS i (red), r (green), and 𝑔 (blue) fluxes and the GALEX FUV flux (additionally enhanced blue color). 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 viewed face-on and edge-on consisting of Spitɀer IRAC 8 μm (red), Herschel PACS 70 μm (green), and Herschel PACS 160 μm (blue) fluxes. The IRAC 8 μm band mainly traces the emission from PAHs heated by radiation from young stars, while the PACS 70 μm and 160 μ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 ≥ 0.5 Gyr, stars produce large dust grains in the spiral arms, and simultaneously accretion and shattering become prevalent in the dense and diffuse ISM, respectively. These processes increase small grains and PAHs, leading to significant emission in the MIR (i.e., the color becomes bluer in the images). At t ≥ 1 Gyr, coagulation also becomes efficient in the dense ISM, enhancing the emission from large grains in the central regions and spiral arms of the galaxy. The emission of small 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 in 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.

Fig. 3 Evolution of the global SED for the Milky Way-like galaxy simulation. The lines in different colors 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.

4.2 Indicators 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 μm band among the bands available for Spitɀer and JWST, we construct the indicators by combining the 8 μ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 μm intensity ratio, I_v(8 μm)/I_v(24 μm), and the other is the 8 μm to the total infrared intensity ratio, vI_v(8 μm)/I_TIR. Here, we denote the intensity per frequency at wavelength λ as I_v(λ). Also, we note that we investigate those indicators with the spatially resolved information of the 50-pc pixel data constructed in Sect. 2.2. The following discussion focuses on the Milky Way-like galaxy but does not significantly change for the NGC 628-like galaxy simulation. We discuss the difference of the PAH indicator, vI_v(8 μm)/I_TIR, between the Milky Way-like and NGC 628-like galaxy simulations in Appendix B.

4.2.1 I_v(8 μm)/I_v(24 μm)

Figure 5 shows the pixel-based relation between PAH mass fraction of q_PAH and I_v(8 μm)/I_v(24μ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⁶ photons cm⁻² s⁻¹ eV⁻¹, U_мw. We fit the data at each time by the following equation: $${\log }_{10}(\text{\hspace{0.17em}Intensity\hspace{0.17em}Ratio\hspace{0.17em}})=k{\log }_{10}\left(q_{\text{PAH}}\right)+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_v(8 μm)/I_v(24 μm) correlates with q_PAH at each time. The 8 μm band is almost purely dominated by the PAH emission, while the 24 μm band is contaminated by emission from larger dust grains with a ≳ 20 Å (Draine & Li 2007), so that I_v(8 μm)/I_v(24 μm) increases with q_PAH. However, we also find that stronger radiation fields decrease the ratio at a fixed q_PAH. This is because stronger radiation fields cause higher dust temperatures, shifting the large grain emission toward shorter wavelengths and raising emission at 24 μm (Draine & Li 2007; Schreiber et al. 2018). Consequently, given the mono-tonic 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_v(8 μm)/I_v(24 μm) and q_PAH at each time. Therefore, I_v(8 μm)/I_v(24 μ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_v(8 μm)/I_v(24 μ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_PAH, for JWST as shown in Fig. 6.

Fig. 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 𝑔 (blue) fluxes with an additionally enhanced blue level for the GALEX FUV flux. The dust emission maps are synthesized by the Spitɀer IRAC 8 μm (blue), Herschel PACS 70 μm (green), and Herschel PACS 160 μm (red) fluxes.

Fig. 5 Pixel-based relation between the 8-to-24 μ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 the density-weighted mean intensity at 1000 Å relative to the Milky Way value as shown by the color bar on the top. Blue dashed lines indicate the best-fit relation given by Eq. (7), and the legends show the values of the coefficients k and c.

4.2.2 vI_v(8 μm)/І_TIR

Here, we use the total infrared intensity estimated by the following equation: $$\begin{array}{c}{I}_{\text{TIR}}=0.95v{I}_v(8\text{\hspace{0.17em}μm})+1.15v{I}_v(24\text{\hspace{0.17em}μm})+v{I}_v(70\text{\hspace{0.17em}μm})\\ +v{I}_v(160\text{\hspace{0.17em}μm}),\end{array}$$(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_MW < 10². Here, we expect to maintain higher spatial resolutions of the observed data by utilizing the estimation with bands at up to 160 μ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 νI_ν(8 μm)/I_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_MW .We find that the ratio tightly correlates with q_PAH and is not affected by the local radiation field of U_MW unlike I_ν(8 μm)/I_ν(24 μm). Remarkably, the parameters of the best-fit relation (blue) are consistent among the different times within a relative difference of 5%. This is because both the total and 8 μ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. Additionally, we perform linear regression analysis for all of the data including the snapshots at all the ages shown above. We find that the PAH mass fraction of qPAH can be estimated using the following equation: $$q_{\text{PAH}}=66{\left(\frac{\nu I_v(8\text{\hspace{0.17em}μm})}{I_{\text{TIR}}}\right)}^{1.2}.$$(9)

In Appendix B, we obtain a similar relation for the NGC 628-like galaxy with a deviation of less than a factor of one at 10⁻¹ < q_PAH < 50. The deviation is greater at lower q_PAH due to contamination of the continuum of non-aromatic small grains at 8 μm. Therefore, we suggest that vI_v(8 μm)/I_TIR is a good indicator of the PAH mass fraction in the ISM of galaxies in various stages of galaxy evolution as long as the PAH emission is prominent enough.

Fig. 6 Same as Fig. 5 but for RPAH.

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

To demonstrate the usefulness of vI_v(8 μm)/I_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. Figure 8 shows the pixel-based relation between vI_v(8 μm)/I_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 vI_v(8 μm)/I_TIR (PAH mass fraction) increases with time and metallicity. This global trend is also supported by larger 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 (the projected maps of the hydrogen surface density, metallicity, and PAH mass fraction are shown in Appendix C). At t = 0.5 Gyr, vI_v(8 μm)/I_TIR increases with metallicity in a lower metallicity regime (Z ≲ 0.2 Z_⊙) while it decreases with metallicity in a higher metallicity regime (Z ≳ 0.2 Z_⊙), 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 (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 high metallicity, dust grows sufficiently through metal accretion, and the effect of coagulation becomes prominent. 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 ≥ 1.0 Gyr, vI_v(8 μm)/I_TIR anticorrelates with metallicity since most of the regions reach high metallicity. We also find that νI_ν(8 μm)/I_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_PAH and oxygen abundance based on the spatially resolved observations of four galaxies with high metallicity (12 + log(O/H) ≳ 8.4). Chastenet et al. (2023a) show a tendency to increase R_PAH with oxygen abundance at low oxygen abundance (12 + log (O/H) ≲ 8.5) but flatten out at high oxygen abundance (12 + log (O/H) ≳ 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 efficiencies of coagulation, shattering, and accretion across diverse ISM conditions, and the analysis of the spatially resolved relation leads us to the understanding of dust evolution.

Fig. 7 Pixel-based relation between the 8 μ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. The blue dashed line at each age indicates the best-fit relation of Eq. (7).

Fig. 8 Time evolution of pixel-based relation among the 8 μ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 ƒdense = 0.5, τSF = 5 Gyr, and UMW = 1.

Fig. 9 Comparison of the global SED between the NGC 628-like galaxy simulation and actual observations. Upper panel: global SED of the NGC 628-like galaxy simulation at t = 10 Gyr shown with the black solid line. 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 (van der Giessen et al., in prep.), respectively. Lower panel: difference in the photometries between the simulation and observation data.

5 NGC 628-like galaxy simulation

NGC 628 is an ideal galaxy for validating our dust model based on spatially resolved data since it is well observed with a high spatial resolution in the PHANGS surveys (Leroy et al. 2021; Lee et al. 2023). Here, we model the NGC 628-like galaxy simulation according to the morphological parameters of 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 prep.).

5.1 SED comparisons

Figure 9 shows the total comparison of SED 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 (Ṁ_* = 0.75 M_⊙ yr⁻¹) is lower than that reported by the observation (Ṁ_* = 1.7 M_⊙ yr⁻¹; Leroy et al. 2021). On the other hand, dust emission at longer wavelengths of λ ≥ 20 μ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.

Fig. 10 Pixel-based comparison of the indicators of the PAH mass fraction between the NGC 628-like galaxy simulation and actual observations. Panel a: vIv(8 μm)/ITIR versus the 160 μm intensity for the NGC 628-like galaxy simulation at t = 10 Gyr (orange) and based on the Herschel and Spitɀer observations for NGC 628 (blue). Panel b: RPAH versus the 21 μ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, respectively, and the shaded area corresponds to the 10th–90th percentile range.

5.2 Spatial comparison of PAH mass fraction

We compare the spatially resolved data of vI_v(8 μm)/I_TIR between our simulation and Herschel/Spitɀer observations. The observations cover large scales up to the galactic radius of r ~ 10 kpc, but the spatial resolution is about 600 pc. Figure 10a shows the comparison of the relation between vI_v(8 μm)/I_TIR and the intensity at 160 μm. The line and shaded area correspond to the median value and the 10th-90th percentile range, respectively. Here, the 160 μm band traces large grains in radiative equilibrium and represents the dust mass distribution. We find that vI_v(8 μm)/I_TIR from the simulation is lower than that of the observation in the entire range of 160 μ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 the fiducial model), we assume eight times the PAH mass in the fiducial model and find that vI_v(8 μm)/I_TIR can be consistent with that of the observation as shown in Fig. 10a. Therefore, we conclude that the PAH mass abundance is underestimated by a factor of eight on average in our simulated area.

Furthermore, we compare the spatially resolved data from our simulation and the JWST observations (Lee et al. 2023), which have 12 times higher spatial resolution than the Herschel observations (a pixel size ~50 pc). Here, instead of I_v(8 μm)/I_v(24 μm), we use the intensity ratio of R_PAH (Eq. (1)). Figure 10b shows the comparison of R_PAH to the intensity at 21 μm between our simulation (orange) and the observation (blue). We find that, in both simulation and observation, R_PAH decreases as the 21 μm intensity increases. This trend can be explained by our model in which coagulation is more efficient in denser regions with higher 21 μm intensities. Figures 11a/c and b/d show the intensity/R_PAH maps at 21 μm of the entire NGC 628-like galaxy at t = 10 Gyr and the zoomed central 6 kpc region, respectively. In Fig. 11c, R_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. 11d), R_PAH is lower in the spiral arms and clumps, where the intensity at 21 μ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. Figures 11e and f show the observed maps of 21 μm intensity and R_PAH maps, respectively (Lee et al. 2023). These observations also indicate that R_PAH is lower in the spiral arms and clumps, where the 21 μ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 21 μ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. 2024). We also find that the value of R_PAH from the simulation is lower than that from the observation, as discussed above: R_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 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, the shape of the grain size distribution for the NGC 628 galaxy at t = 10 Gyr is determined by the balance between coagulation and shattering, and hence, PAHs might be lost too efficiently by coagulation in the 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, because of the higher metallicity in NGC 628 than in the Magellanic Clouds, dust is likely produced more efficiently by accretion than by stellar dust production. Thus, it is not probable that stellar PAH production is dominant over the PAH formation by the shattering of large grains (predominantly formed by accretion) in NGC 628.

In our model, coagulation in the dense ISM efficiently depletes PAHs. Hirashita (2023) constructs a hypothetical model in which small carbonaceous grains are not involved in some interstellar processing, especially coagulation, and succeeds in enhancing the PAH abundance to a level consistent with the observed intensity of PAH emission in the Milky Way. There are some possible mechanisms that suppress the coagulation of small grains. Akimkin et al. (2023) suggest that collisions between charged small grains are inhibited by Coulomb repulsion, whereas collisions between large grains are energetic enough to overcome the Coulomb barrier. This study indicates that a key to the PAH deficiency problem in our model lies in the treatment of grain charges. Alternatively, the interaction of charged grains with magnetic fields may lead to grain velocities that could cause shattering even in the dense ISM (Yan et al. 2004). If PAH formation by shattering occurs in the dense star-forming regions, PAHs could be efficiently illuminated by UV.

PAH formation may also be activated through rotational disruption caused by radiative torques in such illuminated regions (Hoang 2019; Hoang et al. 2019). These additional PAH formation mechanisms in irradiated regions could also contribute to resolving the underestimate of PAH emission. These arguments imply that some processes not involved in our simulations may provide keys to resolving the underestimate of the PAH abundance.

Fig. 11 Comparison of 21 μm intensity and RPAH maps between the NGC 628-like galaxy simulation and the JWST observations. Panels a and b: 21 μm intensity maps of the NGC 628-like galaxy simulation at 10 Gyr on 100 and 20 kpc scales, respectively. Panels c and d: RPAH maps of the simulation on 100 and 20 kpc scales, respectively. Panels e and ſ: maps of the 21 μm intensity and RPAH based on the JWST/MIRI observations for NGC 628 (Lee et al. 2023), where the spatial resolutions are reduced to 50 pc.

6 Conclusion

In this paper, we investigate how dust evolution is reflected on IR spectral features and we provide observational signatures of the grain size evolution in ISM, in particular, by exploring indicators of the PAH mass fraction. We modeled 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 simulated the evolution of grain size distributions in a manner consistent with the physical conditions of the ISM (see Sect. 2.1) by taking into account interstellar dust processes, including stellar dust production, SN destruction, shattering, coagulation, and metal accretion.

Furthermore, we separated the dust compositions into silicate and carbonaceous dust grains, and PAHs according to the abundance of Si and C and the mass fraction of the dense gas-phase in each gas particle in post-processing. Then, we performed radiative transfer calculations with SKIRT under the dust optical properties consistent with the grain size distribution and grain composition of the hydrodynamic simulations (see Sect. 2.2) to investigate the evolution ofSED and emission from the galaxies and indicators of PAH mass fraction. Finally, we compared the indicators of the PAH mass fraction between the NGC 628-like simulation and the observations.

We summarize our results as follows:

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 it contains a larger 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. We investigate the detailed time evolution of observable properties for the Milky Way-like simulation. In the early phase, only the emission of large dust grains can be observed around λ = 70 μm in the SED, as dust production is dominated by stellar sources (see Fig. 3). Afterward, the emission of small grains and PAHs becomes noticeable around λ = 10 μm, since shattering acts as a source of small grains in diffuse regions. In addition, accretion also increases the abundance of small grains (see panels with t = 0.5 Gyr in Fig. 4). Finally, the functional shape of the grain size distribution in each region of the ISM converges to that determined by the balance between shattering and coagulation rate according to the local physical condition. Consequently, the emission from large and small grains is relatively prominent in the inner and outer regions of the galaxy, respectively (see panels with t = 10 Gyr in Fig. 4).

3. We find that I_v(8 μm)/I_v(24 μm) is correlated with the PAH mass fraction, but the local radiation field influences it (see Fig. 5). The 8 μm band is almost purely dominated by PAH band emission, while the 24 μm band is contaminated by the emission of large grains; as a result, I_v(8 μm)/I_v(24μm) increases 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 μm emission. Consequently, I_v(8 μm)/I_v(24μm) decreases in regions with strong radiation. Hence, I_v(8 μm)/I_v(24 μ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. We find that vI_v(8 μm)/I_TIR is tightly correlated with the PAH mass fraction and is not affected by the local radiation field (see Fig. 7 and Eq. (9)). This is because both the total and 8 μm intensities are proportional to the strength of the local radiation field; hence, the effect of the radiation field is canceled out. Therefore, we suggest that vI_v(8 μm)/I_TIR is a good indicator of the PAH mass fraction in various stages of galaxy evolution.

5. We find that the average PAH mass fraction or vI_v(8 μm)/I_TIR in the entire galaxy evolves with metal-licity (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 values of vI_v(8 μm)/ITIR decrease with metallicity in the ISM at higher metallicity (Z > 0.2 Z_⊙), while it increases at lower metallicity (Z ≤ 0.2 Z_⊙). At low metallicity, metal enrichment causes more efficient shattering, which tends to increase the PAHs. At high metallicity, coagulation is rather important in depleting PAHs. We also find that vI_v(8 μm)/I_TIR depends on the gas surface density since coagulation converts PAHs to large grains more efficiently in denser environments. Therefore, the spatially resolved PAH mass fraction helps us interpret the efficiency of accretion, coagulation, and shattering in various ISM phases.

6. We compared the total SED in the NGC 628-like simulation with the observed one (see Fig. 9). Although dust emission at longer wavelengths of λ ≥ 20 μ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. We compared the spatially resolved data of vI_v(8 μm)/I_TIR in the NGC 628-like galaxy simulation and the Herschel/ Spitɀer data (see Fig. 10a). We find that the values of vI_v(8 μm)/I_TIR from the simulation are lower than those from the observations, implying that eight times more PAH mass is required in our simulation.

8. We also compared the spatially resolved data from our simulation with the JWST observations with high spatial resolution (see Fig. 10b). We find that R_PAH, a PAH indicator composed of JWST bands (Eq. (1)), decreases with the 21 μm intensity in both the simulation and observations. This trend can be explained by the fact that coagulation is more efficient in a denser environment. This comparison also shows that the values of R_PAH from the simulation are lower than those from the observations. 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 spatially resolved comparison between observations and theoretical models using observable signatures of grain size evolution (e.g., vI_v(8 μm)/I_TIR, I_v(8 μm)/I_v(24 μm), R_PAH) leads us to a quantitative investigation of dust evolution in various ISM phases. In particular, future JWST observations for a large sample of galaxies with high spatial resolution will systematically reveal the efficiency of dust processes in the diffuse and dense ISM.

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 performed on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, and on the SQUID at the Cyberme-dia Center, Osaka University as part of the HPCI system Research Project (hp220044,hp230089). This work is supported in part by the MEXT/JSPS KAK-ENHI grant numbers 20H00180 and 22K21349 (KN), and 21J20930, 22KJ2072 (YO). KN acknowledges the support from the Kavli IPMU, World Premier Research Center Initiative (WPI), where a part of this work was conducted. H.H. thanks the National Science and Technology Council for support through the grant 111-2112-M-001-038-MY3 and the Academia Sinica for Investigator Award AS-IA-109-M02. M.R. acknowledges the support from the project PID2020-114414GB-100, financed by MCIN/AEI/10.13039/501100011033. S.v.d.G. and I.D.L. acknowledge financial support from the European Research Council (ERC) under de European Union’s Horizon 2020 research and innovation program 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.

Appendix A Initial conditions of the NGC 628-like galaxy

We constructed the initial condition for the NGC 628-like simulation using the DICE code (Perret 2016), which computes the gravitational potential on a multilevel Cartesian mesh and considers the dynamical equilibrium of multiple material components. The process is detailed in the following. 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 × 10¹² M_⊙, a virial velocity of v₂₀₀ = 150 km s⁻¹, and the spin parameter λ = 0.04. These parameters are determined so that the rotation 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 × 10¹² M_⊙ follows the Navarro–Frenk–White (Navarro et al. 1997) profile with a concentration parameter of c = 25 and an effective radius of 100 kpc. The gaseous halo follows the same profile but has 1 % of the total mass, 1.2 × 10¹⁰ M_⊙. 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 densities (S. van der Giessen et al. in prep.) with the Sersic profile given by, $$\text{Σ}(r)={\text{Σ}}_0\exp -b_n\left[{\left(\frac{r}{R_e}\right)}^{1/n}-1\right],$$(A.1)

where r, R_e, ∑₀, 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 describing the steepness of the exponential decline, respectively. The parameters obtained in the fitting are summarized in Table A.1. Furthermore, we assume that the vertical density profile of the galaxy is described by $$\rho (z)\propto \exp \left(-z/h_z\right),$$(A.2)

where ɀ and h_ɀ 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 × 10¹⁰ and 3.5 × 10¹⁰ M_⊙, respectively. Finally, the stellar bulge follows the Einasto profile (Einasto 1965) with a structural parameter of α = 1.0 and a mass of 9.26 × 10⁹ M_⊙.

Figure A.1 shows the comparison of the radial profiles of the hydrogen gas surface densities between NGC 628-like and Milky Way-like galaxies at t = 0 Gyr (blue and orange lines, respectively). For reference, the radial profile from the actual observations for NGC 628 (van der Giessen et al. in prep.) is also indicated by the green circle. Under the initial conditions, the gas radial profile of the NGC 628-like galaxy is flatter than that of the Milky Way-like galaxy, and the gas density in the central region of the NGC 628-like galaxy is not as high. Thus, the star formation of the NGC 628-like galaxy occurs slowly compared to that of the Milky Way-like galaxy (see Section 3).

Fig. A.1 Radial profiles of the hydrogen gas surface densities of NGC 628-like and Milky Way-like galaxies at t = 0 Gyr (orange and blue lines, respectively), compared with that from the actual observations of NGC 628 (green circle; van der Giessen et al. in prep.).

Appendix B Comparison of νI_v(8 μm)/I_TIR between the Milky Way-like and NGC 628-like galaxies

To discuss the robustness of the indicator of the PAH mass fraction, vI_v(8 μm)/I_TIR, in different galaxies, we compare the pixel-based relation for the Milky Way-like and NGC 628-like galaxies, as shown in Fig. B.1. The blue dashed lines in each panel show the linear regression for the relation between vI_v(8 μm)/I_TIR and the PAH mass fraction, which is formulated by the following equation: $$q_{\text{PAH}}=s{\left(\frac{v{I}_{\nu }(8\text{\hspace{0.17em}μm})}{I_{\text{TIR}}}\right)}^w.$$(B.1)

We find that the fitting relation between the Milky Way-like and NGC 628-like galaxies is slightly different Relative to the relation of the Milky Way-like galaxy, that of the NGC 628-like galaxy has a 70, 50, and 30 % difference at q_PAH = 0.1, 1, and 10, respectively This is because the 8 μm band is more contaminated by the continuum of non-aromatic small grains at lower q_PAH, causing the difference of significance of the PAH band emission against the continuum level at 8 μm between the two galaxies Therefore, the relation between vI_v (8 μm)/I_TIR and q_PAH can vary between galaxies, but the variance might be small.

Fig. B.1 Pixel-based relation between the 8 μm-to-total IR intensity ratio and the PAH mass fraction of the Milky Way-like and NGC 628-like galaxies using all snapshots at t = 0.5, 1.0, 3.0, and 10 Gyr (left and right panels, respectively). The color represents the density-weighted mean intensity at 1000 Å relative to the Milky Way value as indicated by the color bar. Blue dashed lines indicate the best-fit relation of Eq. (B.1), and the legend indicates the best-fit parameter.

Appendix C Various maps of the Milky Way-like galaxy

Figure C.1 shows the maps of the hydrogen surface density, metallicity, local radiation field, and PAH mass fraction, q_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_ν > 0.01 MJy sr⁻¹ at all the Spitɀer 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 solar 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_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.

Fig. C.1 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).

References

All Tables

All Figures

	Fig. 1 Time evolution of global physical properties of the Milky Waylike 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.
In the text

	Fig. 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 shades in the same color show the range of 10th–90th percentiles for the sample of gas particles.
In the text

	Fig. 3 Evolution of the global SED for the Milky Way-like galaxy simulation. The lines in different colors 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.
In the text

Fig. 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 𝑔 (blue) fluxes with an additionally enhanced blue level for the GALEX FUV flux. The dust emission maps are synthesized by the Spitɀer IRAC 8 μm (blue), Herschel PACS 70 μm (green), and Herschel PACS 160 μm (red) fluxes.

In the text

Fig. 5 Pixel-based relation between the 8-to-24 μ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 the density-weighted mean intensity at 1000 Å relative to the Milky Way value as shown by the color bar on the top. Blue dashed lines indicate the best-fit relation given by Eq. (7), and the legends show the values of the coefficients k and c.

In the text

	Fig. 6 Same as Fig. 5 but for RPAH.
In the text

Fig. 7 Pixel-based relation between the 8 μ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. The blue dashed line at each age indicates the best-fit relation of Eq. (7).

In the text

Fig. 8 Time evolution of pixel-based relation among the 8 μ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 ƒdense = 0.5, τSF = 5 Gyr, and UMW = 1.

In the text

Fig. 9 Comparison of the global SED between the NGC 628-like galaxy simulation and actual observations. Upper panel: global SED of the NGC 628-like galaxy simulation at t = 10 Gyr shown with the black solid line. 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 (van der Giessen et al., in prep.), respectively. Lower panel: difference in the photometries between the simulation and observation data.

In the text

Fig. 10 Pixel-based comparison of the indicators of the PAH mass fraction between the NGC 628-like galaxy simulation and actual observations. Panel a: vIv(8 μm)/ITIR versus the 160 μm intensity for the NGC 628-like galaxy simulation at t = 10 Gyr (orange) and based on the Herschel and Spitɀer observations for NGC 628 (blue). Panel b: RPAH versus the 21 μ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, respectively, and the shaded area corresponds to the 10th–90th percentile range.

In the text

Fig. 11 Comparison of 21 μm intensity and RPAH maps between the NGC 628-like galaxy simulation and the JWST observations. Panels a and b: 21 μm intensity maps of the NGC 628-like galaxy simulation at 10 Gyr on 100 and 20 kpc scales, respectively. Panels c and d: RPAH maps of the simulation on 100 and 20 kpc scales, respectively. Panels e and ſ: maps of the 21 μm intensity and RPAH based on the JWST/MIRI observations for NGC 628 (Lee et al. 2023), where the spatial resolutions are reduced to 50 pc.

In the text

	Fig. A.1 Radial profiles of the hydrogen gas surface densities of NGC 628-like and Milky Way-like galaxies at t = 0 Gyr (orange and blue lines, respectively), compared with that from the actual observations of NGC 628 (green circle; van der Giessen et al. in prep.).
In the text

Fig. B.1 Pixel-based relation between the 8 μm-to-total IR intensity ratio and the PAH mass fraction of the Milky Way-like and NGC 628-like galaxies using all snapshots at t = 0.5, 1.0, 3.0, and 10 Gyr (left and right panels, respectively). The color represents the density-weighted mean intensity at 1000 Å relative to the Milky Way value as indicated by the color bar. Blue dashed lines indicate the best-fit relation of Eq. (B.1), and the legend indicates the best-fit parameter.

In the text

	Fig. C.1 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).
In the text