© The Authors 2026

1 Introduction

Dust in the interstellar medium (ISM) plays a crucial role in shaping the observable characteristics of galaxies. It absorbs and scatters stellar light across UV to optical wavelengths and re-emits energy in the mid-to-far-infrared (MIR-FIR), significantly affecting a galaxy’s spectral energy distribution (SED). The galaxy attenuation curve, which quantifies the wavelength-dependent attenuation of stellar emission, is essential for deriving intrinsic galaxy spectra. As a result, accurate modeling of the attenuation curve is critical for reliably estimating key properties such as the stellar mass, star formation rate, and age via SED fittings (e.g., Walcher et al. 2011; Conroy 2013).

Although attenuation curves are often discussed alongside extinction curves (see Salim & Narayanan 2020, for a review), they are fundamentally different in nature. Extinction curves are derived along the line of sight by comparing the spectrum of a reddened star to that of an unreddened star of the same type (Cardelli et al. 1989; Fitzpatrick 1999). Extinction curves reflect only dust properties along the line of sight, such as grain size distribution and dust composition, and they are useful for understanding local dust properties. However, extinction curves require resolving individual stars and have thus only been obtained for nearby galaxies such as the MW, Small Magellanic Cloud (SMC), Large Magellanic Cloud (LMC), and M31 (e.g., Koornneef & Code 1981; Nandy et al. 1981; Prevot et al. 1984; Bouchet et al. 1985; Cardelli et al. 1989; Mathis & Cardelli 1992; Fitzpatrick 1999; Valencic et al. 2004; Clayton et al. 2015; Gordon et al. 2023).

In contrast, attenuation curves for galaxies depend not only on dust properties but also on the star-dust geometry and dust scattering, meaning they do not directly reflect the extinction curves of individual stars. The attenuation curve can be estimated with several methods. One method of measuring attenuation curves involves assessing the deviation of multiple Balmer emission line ratios from theoretically estimated intrinsic Balmer emission line ratios (Calzetti et al. 1994; Calzetti 1997; Reddy et al. 2015; Battisti et al. 2016, 2017b; Reddy et al. 2020; Sanders et al. 2025). Using this method, Calzetti (1997) derived the Calzetti attenuation law (A_λ ∝ λ^−0.75) based on the average attenuation of local starburst galaxies (Calzetti et al. 1994, 2000). In addition, numerous studies have fitted parameters of flexible attenuation laws via SED fitting (A_λ ∝ λ⁻ⁿ, Buat et al. 2012, 2018; Salim et al. 2018; Barišić et al. 2020; Nersesian et al. 2024, 2025; Boquien et al. 2022; Fisher et al. 2025; Battisti et al. 2022). This approach was applied to high-redshift galaxies by Markov et al. (2023, 2025a), who adopted a more flexible dust attenuation law with additional free parameters. Wild et al. (2011) employed the pair method, where they compared galaxies with the same type of intrinsic SED but different levels of dust attenuation to derive attenuation curves. Salim et al. (2018) reported a wide diversity in both the slopes and the strengths of the 2175 Å bump among local galaxies. They also identified a correlation between the slope and the V-band attenuation (A_V), a trend that has been confirmed by several other observational studies (e.g., Salmon et al. 2016; Leja et al. 2017; Decleir et al. 2019; Boquien et al. 2022). Battisti et al. (2017a) found that galaxies with lower inclination angles tend to exhibit steeper attenuation curves and weaker 2175 Å bumps (see also Wild et al. 2011; Barišić et al. 2020). While various studies have suggested that the attenuation curve depends on physical parameters such as metallicity, star formation rate, and stellar mass (e.g., Shivaei et al. 2020; Battisti et al. 2020; Salim et al. 2018), the fundamental drivers of these variations remain uncertain (Salim & Narayanan 2020).

The dust attenuation curves of galaxies have been extensively studied across a wide range of redshifts. Battisti et al. (2022) reported that attenuation curves become steeper with increasing redshift up to z < 3, while the strength of the 2175 Å bump remains largely unchanged (see also Zeimann et al. 2015; Scoville et al. 2015). In contrast, recent James Webb Space Telescope (JWST) observations (Markov et al. 2023, 2025a,b) indicate that attenuation curves flatten and tend to lack prominent 2175 Å bumps at higher redshifts, suggesting a dust population dominated by larger grains in the early universe. Similarly, Fisher et al. (2025) found flat attenuation curves in Lyman-break galaxies at z = 6.5–7.7. On the other hand, Witstok et al. (2023) and Ormerod et al. (2025) identified strong 2175 Å bumps in galaxies at z = 6–7, indicating that small grains were already present in some early systems. Collectively, these findings underscore the critical role of dust properties, particularly the grain size distribution, in shaping the attenuation curves of galaxies in the early universe.

To better interpret these observations, it is crucial to theoretically investigate physical mechanisms that shape attenuation curves. Although radiative transfer calculations based on toy-model geometries have provided valuable insights into the underlying physics of attenuation curves (e.g., Witt & Gordon 2000; Gordon et al. 2001; Chevallard et al. 2013; Seon & Draine 2016; Nersesian et al. 2020; Lin et al. 2021), those based on hydrodynamical simulations allow for the study of attenuation curves in the context of more realistic and self-consistent galaxy structures. Trayford et al. (2020) found a strong relation between A_V and the Near-UV slope based on attenuation curves in the EAGLE simulation (Schaye et al. 2015), which varies with galaxy inclination. Sommovigo et al. (2025) showed a correlation between A_V and the star formation rate surface density (Σ_SFR) in galaxies at z = 0.07 from the TNG50 and TNG100 simulations (Pillepich et al. 2018, 2019), suggesting that galaxies with higher star formation rate densities exhibit higher levels of dust attenuation and have flatter attenuation curves. Narayanan et al. (2018) demonstrated that star-dust geometry primarily drives the slopes of attenuation curves in galaxies at z = 0–6 from the MUFASA simulations (Davé et al. 2016): galaxies with less obscured young stars show flatter attenuation curves. Both Narayanan et al. (2018) and Sommovigo et al. (2025) observed a decreasing trend in the strength of the 2175 Å bump with increasing A_V. These studies highlight the importance of the star–dust geometry in shaping both the slope and the 2175 Å bump. However, they assume a fixed dust-to-metal mass ratio and a static dust model, where the grain size distribution and dust composition are not variable. This assumption may introduce uncertainties in the predicted dust distributions of galaxies.

Many studies have modeled physically motivated dust mass evolution alongside hydrodynamic simulations for galaxies (Li et al. 2019; Graziani et al. 2020; Choban et al. 2022, 2024, 2025), as well as the evolution of grain size distribution (McKinnon et al. 2018; Aoyama et al. 2017, 2020; Hou et al. 2019; Romano et al. 2022a; Narayanan et al. 2023; Li et al. 2021; Dubois et al. 2024). These models highlight the nonlinear evolution of the dust-to-gas and dust-to-metal mass ratios, which vary according to the conditions of the ISM, where dust grains grow and are destroyed, as also inferred from several observational studies (e.g., Rémy-Ruyer et al. 2014; De Vis et al. 2019; Galliano et al. 2021; Park et al. 2024; Algera et al. 2025). These studies emphasize the importance of considering dust mass evolution. Li et al. (2021) suggested that both grain size distribution and the graphite-to-silicate ratio influence the slope and 2175 Å bump in extinction curves. Similarly, Dubois et al. (2024) showed that SMC-like curves occur at low metallicity and MW-like curves occur at high metallicity. However, the influence of evolving dust mass and grain size distribution, in conjunction with the evolving stellar populations and their geometry, on galaxy attenuation curves remains largely unexplored.

In this study, we explore the impact of dust mass and grain size evolution on attenuation curves using a MW-like galaxy simulation. We further examine the roles of scattering, star-dust geometry, and grain size distribution based on hydrodynamic simulations with dust evolution (Matsumoto et al. 2024). This paper is structured as follows. Section 2 provides an overview of the physical models used in our hydrodynamic simulations with GADGET4-OSAKA and the methodology for post-processing radiative transfer with SKIRT. In Sect. 3, we examine the time evolution of the average extinction curve of the MW-like galaxy simulation, focusing on how the evolution of the grain size distribution impacts the extinction curves. Section 4 explores the time evolution of the global attenuation curve, identifying the primary factors that influence the attenuation curve, including scattering, star-dust geometry, and grain size distribution. In Sect. 5 we compare our simulations with previous theoretical and observational studies, and present physical insights into dust attenuation curve variations. Finally, in Sect. 6, we present our conclusions.

2 Method

2.1 A MW-like galaxy simulation with GADGET4-OSAKA

We modeled an isolated MW-like galaxy using the methodology outlined in (Matsumoto et al. 2024). We summarize the methodology here, referring to Romano et al. (2022a); Matsumoto et al. (2024) for more details. This simulation was performed using a smooth particle hydrodynamic (SPH) simulation code, GADGET4-OSAKA (Romano et al. 2022a,b; Oku et al. 2022; Oku & Nagamine 2024), a modified version of the GADGET4 code (Springel et al. 2021), incorporating the OSAKA feedback models (Shimizu et al. 2019) and an evolutionary model for grain size distributions in the ISM (Aoyama et al. 2020; Romano et al. 2022a). The simulation follows the evolution of grain size distributions across 30 bins ranging from 3.0 × 10⁻⁴ to 10 µm in a MW-like galaxy over 10 Gyr¹, taking into account interstellar dust processes and dust production from asymptotic giant branch (AGB) stars and supernovae (SNe). The initial condition for the MW-like galaxy simulation is taken from the AGORA project² (Kim et al. 2016). The components of the initial conditions of the MW-like galaxy simulation are summarized in Table 1. The gas and stellar metallicity are initially set at 0.0001 Z_⊙, where Z_⊙ = 0.013 (Asplund et al. 2009), and increase to 0.7 Z_⊙ by 10 Gyr. The gravitational softening length is set to 40 pc, and the minimum smoothing length is allowed to reach 0.1 ϵ_grav.

Star formation occurs stochastically, governed by a star formation efficiency of ϵ_∗ = 0.05, as per the Kennicutt–Schmidt law (Kennicutt 1998). Star formation occurs when the density of gas particles exceeds n_H > 20 cm⁻³ and the temperature of the gas is below T_gas < 10⁴ K. In the OSAKA feedback model, the feedback processes of type Ia SNe, type II SNe, and AGB stars affect the physical conditions of the ISM within the galaxy (Shimizu et al. 2019) via thermal and kinetic feedback processes. Simultaneously, SNe and AGB stars eject metal elements and dust.

The ejected metal mass is calculated using the CELib library (Saitoh 2017), and 10% of the ejected metals are assumed to condense into dust grains (Asano et al. 2013). We employed the GRACKLE-3 radiative cooling library (Smith et al. 2017)³, which solves the nonequilibrium primordial chemical network, taking into account metal cooling as well as photo-heating and photo-ionization due to the UV background (Haardt & Madau 2012).

Our simulations lacked the necessary spatial resolution to accurately resolve cold and dense gas (T_gas < 50 K and n_H > 10³ cm⁻³). Therefore, we implemented a two-phase ISM sub-grid model for gas particles (Romano et al. 2022a), assuming that each gas particle contains both dense and diffuse gas phases. The mass fraction of the dense gas phase is given by $$f_{\text{dense}}=\text{min}\left(\alpha \frac{n_{\text{H}}}{1.0{\text{cm}}^{-3}},1.0\right),$$(1) for which we adopted α = 0.12. This equation is based on the observation that regions with higher gas concentrations in the ISM tend to host a larger fraction of dense molecular gas clouds (Gnedin & Kravtsov 2011; Gnedin & Draine 2014). Romano et al. (2022a) calibrated the parameter α to reproduce the observed mean molecular gas fraction for the MW-mass galaxies with stellar masses around 10¹⁰ M_⊙ (Catinella et al. 2018). While the temperature and density of the dense gas phase were assumed to be constant (T_dense = 50 K and n_{H, dense} = 10³ cm⁻³), the temperature and density of the diffuse gas phase in each gas particle were determined by considering the conservation of internal energy. This subgrid model enabled us to separately treat dust processes in the diffuse and dense ISM, accounting for their distinct properties within each gas particle.

We took into account various interstellar dust processes in both the diffuse and dense gas phases. In the ISM, shattering breaks down larger grains into smaller ones, while coagulation leads to the growth of small grains into larger ones. Shattering and coagulation are assumed to occur through grain-grain collisions in the dense and diffuse gas phases, respectively. Both processes are influenced by the grain size and velocities, with the grain velocity determined by its coupling to turbulent motion, which follows a Kolmogorov spectrum (Kolmogorov 1941) and reflects the local physical conditions. Additionally, metal accretion onto dust grains predominantly takes place in the dense gas phase, causing an increase in the size and mass of the grains. The accretion rate is determined by the collision frequency between gas-phase metals and dust, and it is parameterized by the dust-to-metal ratio, grain size, temperature, and metallicity. On the other hand, sputtering, which occurs in the diffuse gas phase, reduces the size and mass of each grain. The sputtering rate in hot gas is based on Tsai & Mathews (1995), and it is also influenced by the sweeping rate of supernova shocks that contribute to the destruction of dust. The chemical composition of grains is not explicitly included in the dust evolution model of the GADGET4-OSAKA and, thus, we modeled the dust composition of each dust grain in post-processing (see Sect. 2.2). We refer to Aoyama et al. (2020) and Romano et al. (2022a) for more details on the treatment of dust processes in the simulation.

2.2 Post-processing radiative transfer with SKIRT

We performed radiative transfer in post-processing based on the snapshots of the MW-like galaxy simulation with the radiative transfer code, SKIRT⁴ (version 9.0; Baes et al. 2011; Camps & Baes 2015, Camps & Baes 2020), to generate attenuation curves. SKIRT has been widely used for dust radiative transfer calculations (Camps et al. 2018; Liang et al. 2018; Ma et al. 2019; Parsotan et al. 2021; Cochrane et al. 2019, 2022, 2023a,b, 2024; Camps et al. 2022; Kapoor et al. 2021, 2023, 2024; Trčka et al. 2022; Baes et al. 2024a,b; Gebek et al. 2024, 2025; Sommovigo et al. 2025), although it is 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; Matsumoto et al. 2023; Vander Meulen et al. 2023, 2024a,b).

First, we used the Bruzual & Charlot (2003) templates with the Chabrier (2003) initial mass function for single stellar populations. The intrinsic SED of each star particle is determined based on its stellar age, metallicity, and initial stellar mass. 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. Since the simulations run until t = 10 Gyr, and we assume the final galaxy age to be comparable to the cosmic age (13.8 Gyr), the disk and bulge stars are assumed to have formed between t = −3.8 and t = 0 Gyr. Therefore, in post-processing, we randomly assigned their ages to obtain a constant star formation history within the 3.8 Gyr interval prior to t = 0 Gyr. We note that the pre-assignment of stellar age did not significantly impact the evolution of the galaxy attenuation curves on Gyr timescales, as the stellar emission is primarily dominated by the stars formed during the simulation. Here, we assume that star particles are point sources with a fixed size of 0.1 pc. This assumption is reasonable because star particles are treated as collisionless particles in the GADGET4-Osaka simulation and their assumed size does not significantly affect the galaxy SED.

Second, for modeling the dust composition, we decomposed the dust grains in each gas particle into silicate, carbonaceous dust, and PAHs in post-processing. Here, we follow the method of decomposition of dust grains into the various chemical compositions described in Matsumoto et al. (2024). First, we assigned larger dust grains (with sizes a > 0.0013 µm) as silicate and carbonaceous dust grains according to the abundances of Si and C in each gas particle (Hirashita & Murga 2020). We distinguished among small and large grains using a grain size boundary of a = 0.0013 µm, derived from the definition of the maximum size of PAHs as described in Draine & Li (2001) and Draine et al. (2007). Small grains (with sizes a ≤ 0.0013 µm) were designated as PAHs and nonaromatic carbonaceous grains. The mass fraction of nonaromatic carbonaceous grains can be determined by the dense gas fraction in each gas particle (f_aroma = 1 − f_dense; Hirashita & Murga 2020). Based on this assumption, the majority of small grains with a ≤ 0.0013 µm were assigned to PAHs, as these small grains are primarily produced in the diffuse ISM. Finally, to decompose PAHs into neutral (PAH⁰) and ionized (PAH⁺) forms, we used the ionization fraction derived from the “standard” ionization model in Draine et al. (2021) and Hensley & Draine (2023). As a result, PAH⁺ are primarily assigned to the smallest grains (a ∼ 0.0003 µm).

The opacity and calorimetric properties of PAHs, silicate, and nonaromatic carbonaceous grains (assumed to be graphite) follow the prescription of Draine & Li (2007). The upper panels of Fig. 1 show the extinction curve for carbonaceous and silicate grains with a given grain size (upper left and right, respectively) taken from Draine & Li (2007). The extinction curve is derived from the wavelength dependence of the dust grain opacity, κ_λ, and is normalized at the V-band wavelength (λ = 5500 Å). As the grain size increases, the extinction curves for silicate and graphite grains become flatter. The 2175 Å bump strength of graphite grains develops when the grain size is lower than 0.05 µm. PAH shows a steep extinction curve with strong 2175 Å bumps compared to the extinction curve for graphite grains. We note that the extinction curves for graphite grains with a size of 10⁻⁴ µm are nearly identical to those for grains of 10⁻³ µm in size.

The lower panels of Fig. 1 show the wavelength dependence of the albedo, which is the fraction of scattering against total extinction, for carbonaceous and silicate grains with a given dust grain size (lower left and right, respectively). The albedo of graphite increases with grain size in the optical-to-UV wavelength range, while smaller grains exhibit lower albedo at longer wavelengths. Smaller silicate grains show higher albedo around 2000 Å, whereas larger silicate grains, up to a size of 0.2 µm, have higher albedo in the optical wavelengths. For even larger silicate grains, the albedo tends to decrease at optical wavelengths. Thus, larger grains are responsible for scattering at optical wavelengths. Additionally, silicate grains tend to cause more scattering compared to graphite. We note that graphite and silicate grains with sizes smaller than 10⁻⁴ µm, as well as PAHs, do not contribute to scattering, resulting in an albedo of 0.0 across all wavelengths.

Given the dust and stellar distributions from the GADGET4-OSAKA simulation and dust properties from Draine & Li (2007), we performed post-processing radiative transfer calculations with SKIRT, incorporating spatially varying grain size distributions and chemical composition of the grains. For the SKIRT simulations, the dust distribution of SPH particles, smoothed with the cubic kernel function, is resampled onto an octree grid with refinement based on a dust mass fraction threshold of 10⁻⁶ relative to the total dust mass. We show the time evolution of dust mass and star formation history of the MW-like galaxy simulation in Appendix A. For generating synthetic images with SKIRT, the field of view with (40 kpc)² is composed of spaxels with a pixel size of 50 pc and a spectral resolution of $R=\frac{\lambda }{\Delta \lambda }=500$ in the wavelength ranges from 1000 to 10 000 Å. We set the inclination angles from 0^◦ to 90^◦ in increments of 10^◦ (from face-on to edge-on views).

Fig. 1 Upper panels: extinction curves for carbonaceous and silicate grains with a given dust grain size (left and right, respectively). The extinction curve is normalized by the opacity at the V-band wavelength (λ = 5500 Å). Lower panels: wavelength dependence of albedo, which is the fraction of scattering against extinction, of carbonaceous and silicate grains (left and right, respectively).

2.3 Modeling global extinction and attenuation curves

We assume spherical grains with material density, ρ, such that the grain mass, m(a), and radius, a, are related as $m(a)=\frac{4}{3}\pi a^3\rho$. The grain size distribution, denoted as $\frac{dn(a)}{da}$, is defined so that n(a) is the number distribution of grains as a function of grain radius a. Using the global grain size distribution across the entire galaxy, we calculated the average extinction curve for each dust component, normalized by the total extinction at the V-band (λ = 5500 Å), as follows, $$\frac{{\tau }_{\lambda }^{\text{comp}}}{{\tau }_V^{\text{total}}}=\frac{{{\displaystyle \int }}^{\text{}}{\kappa }_{\lambda }^{\text{comp}}(a)m^{\text{comp}}(a){\left(\frac{dn(a)}{da}\right)}_{\text{comp}}da}{{{\displaystyle \sum }}_{\text{comp}}{{\displaystyle \int }}^{\text{}}{\kappa }_V^{\text{comp}}(a)m^{\text{comp}}(a){\left(\frac{dn(a)}{da}\right)}_{\text{comp}}da},$$(2) where ${\kappa }_{\lambda }^{\text{comp}}(a)$ represents the extinction opacity for each dust component as a function of grain size at a given wavelength. The total average extinction curve is given by ${\tau }_{\lambda }^{\text{total}}/{\tau }_V^{\text{total}}$ = $\Sigma {\tau }_{\lambda }^{\text{comp}}/{\tau }_V^{\text{total}}$ and influenced solely by the grain size distribution and the chemical composition of the dust.

On the other hand, the attenuation curve is estimated based on the observational properties given by $$A_{\lambda }=-2.5\text{log}\left(\frac{F_{\lambda }^{\text{observed}}}{F_{\lambda }^{\text{intrinsic}}}\right),$$(3) where $F_{\lambda }^{\text{observed}}$ and $F_{\lambda }^{\text{intrinsic}}$ are the observed and intrinsic flux densities of the galaxy, respectively. We note that the attenuation is not an additive property, as it is influenced not only by the dust properties but also by scattering and star-dust geometry.

We parameterized the extinction and attenuation curves using the UV-to-optical slope and 2175 Å bump strength (hereafter referred to as the slope and bump, respectively). The slope is defined as $$S=\frac{A_{FUV}}{A_V},$$(4) where A_FUV and A_V represent the attenuation at FUV- and V-bands (λ = 1535 and 5500 Å), respectively. To derive the bump strength, we fitted the extinction and attenuation curves by the following function from Li et al. (2008), expressed as $$\begin{array}{lll}\frac{A_{\lambda }^{\text{fit}}}{A_V}\hfill & =\hfill & \frac{c_1}{{(\lambda /0.08)}^{c_2}+{(0.08/\lambda )}^{c_2}+c_3}\hfill \\ \hfill & \hfill & +\frac{233[1-c_1/({0.145}^{-c_2}+{0.145}^{-c_2}+c_3)-c_4/4.60]}{{(\lambda /0.046)}^2+{(0.046/\lambda )}^2+90}\hfill \\ \hfill & \hfill & +\frac{c_4}{{(\lambda /0.2175)}^2+{(0.2175/\lambda )}^2-1.95},\hfill \end{array}$$(5)

where c₁, c₂, c₃, and c₄ are dimensionless parameters and λ is the wavelength in µm units (Markov et al. 2023, 2025a; Sommovigo et al. 2025). The way these parameters are reflected in the attenuation curves is discussed in Sommovigo et al. (2025). Consequently, the bump strength is defined as $$B=\frac{A_{2175Å}-A_{2175Å,0}}{A_{2175Å}},$$(6) where A_{2175Å, 0} is the baseline of the attenuation at 2175 Å, calculated by $A_{2175Å}^{\text{fit}}\left(c_1,c_2,c_3,c_4=0\right)$. For the extinction curves, the A_λ in Eqs. (4), (5), and (6) is replaced by τ_λ.

Fig. 2 Time evolution of the global grain size distribution in our MW-like galaxy simulation. The solid and dashed gray lines correspond to the grain size distributions of carbonaceous and silicate grains in the static dust model from Draine & Li (2007).

3 Time evolution of the average extinction curve

First, we briefly introduce the global evolution of the grain size distribution in the MW-like galaxy simulation shown in Fig. 2. Full details of the dust evolution scenario are described in Romano et al. (2022a) and Matsumoto et al. (2024). The solid lines in Fig. 2 represent the grain size distribution of the MW-like galaxy at various simulation times (t = 0.1, 0.25, 0.5, 3.0, and 10.0 Gyr). We note that in our model, the evolution of the grain size distributions for silicate and carbonaceous grains does not differ significantly. At t = 0.1 Gyr, only large dust grains (a ∼ 0.2 µm) are produced by type II SNe in the galaxy. Then, shattering converts the large grains to small ones in the diffuse ISM, and accretion helps dust grains to increase their sizes. Thus, small grains (a < 0.1 µm) gradually form in the galaxy up to t = 0.5 Gyr. Later, small grains start being coagulated to large grains (a > 0.1 µm) in the dense ISM, while shattering continues producing small grains in the diffuse ISM.

Figure 3 shows the evolution of the average extinction curve of the MW-like galaxy simulation, derived using Eq. (2). The extinction curve purely reflects the global grain size distribution and chemical composition of dust grains within the MW-like galaxy simulation. Although it is impossible to observationally obtain the average extinction curve, it is essential to theoretically explore how the evolution of grain size distributions is reflected in the average extinction curve. The average extinction curve becomes progressively steeper up to t = 0.5 Gyr along with a development of the 2175 Å bump, driven by the production of smaller grains through shattering and accretion. After that, the extinction curve becomes gradually flatter up to t = 3 Gyr because small grains are coagulated. Later, the extinction slightly steepens until t = 10 Gyr because of the formation of very small grains including PAH through shattering. At all times, the evolution of the 2175 Å bump strength is coherent with that of the slope. The detailed evolution of the slope and the 2175 Å bump strength is presented in Sects. 4.1 and 4.2.

Figure 4 shows the contribution of each dust component to the average extinction curve of the MW-like galaxy simulation at various times. At t = 0.1 Gyr, the extinction curve is mainly driven by silicate at all wavelengths. At this time, gas-phase silicon and carbon are produced by SNe II. After t = 0.5 Gyr, the extinction curve at shorter wavelengths is determined by graphite, since gas-phase carbon is produced by SNe Ia (Saitoh 2017). Also, the 2175 Å bump becomes prominent due to the formation of small graphite grains (a < 0.05 µm). At t = 3 Gyr, the extinction curve flattens because of the coagulation of small grains. Later, the bump becomes stronger, and the slope steepens thanks to the production of very small grains, especially PAHs. Overall, silicate grains are responsible for the extinction curve at FUV wavelengths at all times, while carbonaceous grains become the dominant contributors to the extinction curve at wavelengths longer than 2000 Å, similar to previous studies (Weingartner & Draine 2001; Jones et al. 2013; Dubois et al. 2024).

Fig. 3 Time evolution of the average extinction curve of the MW-like galaxy simulation indicated by colored solid lines. The extinction curve is estimated by Eq. (2) and normalized by the opacity at the V-band wavelength (λ = 5500 Å). The dashed gray, dotted orange, and dash-dotted red lines represent the extinction curves of the MW, SMC, and LMC, respectively (Pei 1992; Li et al. 2008).

4 Time evolution of the global attenuation curve

The attenuation curve is influenced by various factors such as (i) scattering (Chevallard et al. 2013), (ii) star-dust geometries (Seon & Draine 2016; Narayanan et al. 2018), and (iii) dust properties including grain size distribution and chemical composition. Ultimately, these factors can be reflected in the relation between the attenuation curve and global galaxy properties and inclination effects (Trayford et al. 2020; Sommovigo et al. 2025).

However, the question of which particular factor shapes the attenuation curve is not fully understood. To examine and isolate these competing factors, we generated attenuation curves using two dust models: the first is a dynamic dust model utilizing dust grain size distribution and chemical composition from our GADGET4-OSAKA simulation (the details of the decomposition are described in Sect. 2.2); the second is a static dust model with a fixed grain size distribution and composition based on the MW model from Draine & Li (2007). Although the second model accounts for dust mass evolution, it does not use the grain size distribution from the GADGET4-OSAKA simulation. A comparison of the grain size distributions of the two dust models is shown in Fig. 2. Moreover, we generated attenuation curves considering scattering explicitly and handling scattering as absorption (i.e., without scattering) to probe the impact of scattering on attenuation curves. We note that, even when distinguishing between these models, the scattering processes still depend on the star–dust geometry and dust properties. Consequently, four types of attenuation curves are created and summarized below:

1. Static-no-scattering: the attenuation curves are derived from the static dust model, excluding scattering effects. The attenuation curves are affected only by the star-dust geometry.

2. Static-with-scattering: the attenuation curves are derived using the static dust model, and scattering is explicitly taken into account. The attenuation curves are affected by the stardust geometry and scattering.

3. Dynamic-no-scattering: the attenuation curves are derived using the dynamic dust model, excluding scattering effects. The attenuation curves are driven by the star-dust geometry and dust properties including grain size distribution and dust composition.

4. Dynamic-with-scattering: the attenuation curves are derived using the dynamic dust model, and scattering is explicitly taken into account. The attenuation curves are determined by a combination of star-dust geometry, scattering, and evolving dust grain properties. This is the most physically comprehensive model of the four considered here.

First, we present an overview of the time evolution of the attenuation curve and its dependence on different inclination angles. Figure 5 shows the evolution of the attenuation curves of our MW-like galaxy simulation with the dynamic-with-scattering and static-with-scattering models. In both models, the attenuation curves for an inclination angle of 90^◦ (i.e., edge-on view) are flatter than those at 0^◦ (i.e., face-on view), and the 2175 Å bumps become less pronounced. The attenuation curves for the static-with-scattering model become progressively steeper over time in the face-on view, despite assuming a constant grain size distribution and dust composition, while the slope remains nearly constant in the edge-on view. The dynamic-with-scattering model shows a similar evolutionary trend. We caution that spikes in the attenuation curves are caused by star–dust geometry effects. Lyα and Lyβ absorption lines from older stellar populations, as modeled in the Bruzual & Charlot (2003) SSPs, appear as convex-shaped features in the attenuation curves. In contrast, absorption lines originating exclusively from younger stars present as concave-shaped features. We quantify the time evolution of the slope and bump strength more accurately in the following sections where we discuss the following factors: (i) scattering, (ii) star-dust geometry, and (iii) dust properties by examining the difference in the slope and bump strength of the attenuation curves from the four different models (dynamic-with-scattering, dynamic-no-scattering, static-with-scattering, and static-no-scattering).

Fig. 4 Contribution of each dust component to the average extinction curve of the MW-like galaxy simulation at various times (t = 0.1, 0.5, 3.0, and 10.0 Gyr). The black thick line is the total average extinction curve, while the orange solid, green dashed, blue dotted, and light blue dash-dotted thin lines represent the contributions of silicate, graphite, PAH0, and PAH+, respectively. At t = 0.1 Gyr, the extinction curve is dominated by silicate grains across all wavelengths. However, at later times, carbonaceous grains become the dominant contributors to extinction at wavelengths longer than 2000 Å.

Fig. 5 Time evolution of the global attenuation curve for the MW-like galaxy simulation using the static-with-scattering and dynamic-with-scattering models shown in left and right panels, respectively. The solid and dashed lines represent inclination angles of 0◦ and 90◦ (i.e., face-on and edge-on view, respectively). The attenuation curve is estimated using Eq. (3) and normalized by the attenuation at the V-band wavelength (λ = 5500 Å).

4.1 Time evolution of the slope at various inclination angles

We examine the time evolution of the slope of the global attenuation curve of the MW-like galaxy simulation by comparing that of the extinction curve. Figure 6 presents the time evolution of the slope of the global attenuation curve for the four different models of the MW-like galaxy simulation at various inclination angles. Firstly, in the static-no-scattering model, the slope decreases until 250 Myr and then slightly rises. At all times, most of the intrinsic stellar emission at UV and optical wavelengths comes from young and old stars (stellar age < 10 and >10 Myr, respectively; Narayanan et al. 2018). At t = 100 Myr, dust is produced only by SNe and hence preferentially obscures only young stars in the star-forming region. Later, dust starts growing in the ISM due to metal accretion on dust grains and attenuates emission from old stars, increasing A_V and flattening the slope until t = 250 Myr. At higher inclination angles, older stars are likely to be obscured more by the dust along the line of sight, and hence, the slope further flattens until t = 500 Myr. At later times, dust becomes optically thick to part of the UV emission from young stars, and only a small amount of UV emission from young stars leaks from the galaxy. Consequently, the UV observed spectrum is more contributed to by the emission from old stars. The intrinsic spectra of old stars become redder as the simulation time increases. As a result, the difference between the intrinsic UV spectra from young stars and the observed UV spectra from older stars leads to higher A_FUV and a steeper attenuation curve (see also Sect. 4.3 in Narayanan et al. 2018). After t = 7 Gyr, the intrinsic UV spectra are also increasingly dominated by older stars, resulting in lower A_FUV and slightly flattening the slope.

In the dynamic-no-scattering model, the slope increases up to t = 0.5 Gyr, following the trend of the extinction slope. This is due to the formation of small grains (a < 0.1 µm) via the interplay between shattering and accretion in the ISM. Indeed, small grains absorb UV emission more significantly than larger dust grains, increasing A_FUV and the slope. After t = 0.5 Gyr, the small grains coagulate, forming larger grains, and as a consequence, the slope slightly decreases. Here, silicate grains are responsible for the attenuation at UV wavelengths and the change in the slope, as discussed in the evolution of the extinction curve in Sect. 3. Before t = 1 Gyr, the evolution of the slope follows the trend of the average extinction curve, as shown by the gray thick line in Fig. 6. Later, the grain size distribution becomes relatively stable over time, and the post-processed assigned ratio of silicate to carbonaceous grains shows little variation (see the bottom panel of Fig. A.1). Consequently, the time evolution of the dynamic-no-scattering model closely resembles that of the static-no-scattering. As in the static-no-scattering model, the slope in the dynamic-no-scattering model decreases with inclination angle.

Once scattering is included (lower panels of Fig. 6), the time evolution of the slope deviates significantly from that of the extinction curve, with the extent of deviation depending on the inclination angle. The slopes for both dynamic-with-scattering and static-with-scattering models are significantly higher compared to the no-scattering models, especially at lower inclination angles. When the dust column density along the line of sight is low, optical photons can escape more easily from the galactic plane along the line of sight thanks to scattering, effectively reducing attenuation at optical wavelengths (A_V). However, the dust opacity at UV wavelengths is much higher than that at optical wavelengths, meaning that UV photons are less likely to escape even for the same dust column density independently from accounting for scattering or not. Thus, the attenuation at optical wavelengths is effectively reduced by scattering, resulting in an overall steepening of the attenuation curves. The effect of scattering is maximal when line-of-sight optical depth at optical wavelengths is low (A_V < 1) and, hence, the slope becomes steeper at lower inclination angles. In fact, the attenuation curves with scattering are almost identical to those without scattering at the inclination angle of 90^◦. In addition, we also note that, as UV photons have stronger forward scattering by dust, they are unlikely to change their path and escape from the galactic disk (Draine 2003).

In the static-with-scattering model, the slope becomes flatter up to t = 1 Gyr, since the effectiveness of scattering at optical wavelengths is reduced as the optical depth along the line of sight increases due to dust formation via metal accretion. After t = 1 Gyr, the slope becomes significantly steeper over time. Dust formation shifts from the center to the outskirts of the galaxy over time (see Fig. D.1) as the total dust mass monotonically increases. Consequently, the enhancement of dust obscuration of stars within the galaxy increases the possibility of scattering in the galactic plane, steepening the slope over time. This result suggests that the effect of scattering on attenuation curves is not solely determined by the dust column density along the line of sight. Rather, the dust column density perpendicular to the line of sight also plays a significant role in shaping the attenuation curve (we discuss the three-dimensional (3D) effect of scattering in Sects. 5.1 and 5.2 thoroughly).

The time evolutionary trend of the dynamic-with-scattering model is similar to that of the static-with-scattering model. However, the slope for the dynamic-with-scattering model starts with lower values at higher inclination angles due to the lack of small grains at t < 0.5 Gyr and the suppression of scattering processes. Furthermore, overall, the slope for the dynamic-with-scattering model is flatter compared to that for the static-with-scattering model. This is because large silicate dust grains (a > 0.1 µm) are more populated in the static-with-scattering model (Draine & Li 2007). Large silicate grains have a high albedo and are likely to scatter at optical wavelengths (V-band), resulting in the steeper attenuation curve for the static-with-scattering model. These effects demonstrate the importance of the dust grain properties in shaping the attenuation curves.

In brief, in our simulations, the evolutionary trend of the slope for the dynamic-with-scattering model differs from that of the underlying extinction curve due to the interplay of multiple effects. At low inclination angles, the attenuation curves are strongly shaped by scattering, which is sensitive to changes in the star-dust geometry. At higher inclination angles, the effect of the scattering is suppressed, and thus, the attenuation curves are dominated by the variations in the star-dust geometry and grain size distribution. Importantly, changes in grain size distribution affect not only the underlying extinction curve, but also scattering processes, both of which influence the slope of the resulting attenuation curves. Appendix C further discusses the impacts of star-dust geometry, scattering, and dust properties on A_V and A_FUV.

Fig. 6 Time evolution of the slope for the global attenuation curves from the four different models (dynamic, dynamic-no-scattering, static, and static-no-scattering) of the MW-like galaxy simulation. The inclination angle is color-coded in each panel. Here, all models with Av < 0.01 are excluded. The gray thick lines represent the time evolution of the slope of the average extinction curves for the static and dynamic dust models in the left and right panels, respectively. The value of slope for the static dust model corresponds that of the MW (Pei 1992; Li et al. 2008).

Fig. 7 Same as Fig. 6, but considering the time evolution of the 2175 Å bump strength.

4.2 Time evolution of the bump at various inclination angles

We examine the main factors driving the time evolution of the 2175 Å bump of the global attenuation curve in the MW-like galaxy simulation by comparing it with that of the extinction curve. Figure 7 shows the time evolution of the bump of the global attenuation curve of the MW-like galaxy simulation at various inclination angles.

In the static-no-scattering model, the 2175 Å bump strength is initially as high as in the extinction curve from t = 0.05 to 0.1 Gyr. Thereafter, it gradually decreases until t = 3 Gyr and then remains nearly constant up to t = 10 Gyr. Initially, dust remains optically thin to stellar emission at λ = 2175 Å, meaning that most stars contribute to the observed emission with minimal differential attenuation across different lines of sight. As a result, the shape of the attenuation curve closely follows that of the intrinsic extinction curve, including the 2175 Å bump. However, after t = 0.1 Gyr, as dust column density increases, the UV optical depth exceeds unity (A_2175Å > 1.0) along some sight lines. This leads to a stronger bias toward emission from u006Eobscured stars in the observed UV light (Narayanan et al. 2018), causing the observed attenuation curve to deviate from the intrinsic extinction curve and reducing the apparent strength of the 2175 Å bump. Consequently, the bump strength becomes weaker over time. At higher inclination angles, the area where dust is optically thick to the stellar emission at 2175 Å increases, weakening the bump strength more.

In the dynamic-no-scattering model, the 2175 Å bump strength gradually increases until t = 0.5 Gyr due to the formation of small carbonaceous grains in the ISM via shattering and accretion. Here, small carbonaceous dust grains, including PAH and graphite, are responsible for the change in the bump strength as discussed in Sect. 3. After t = 1.0 Gyr, the bump strength becomes weaker because small grains are coagulated, and larger grains form. After 3 Gyr, very small carbonaceous grains such as PAH (a < 0.0012 µm) begin to form through shattering, strengthening the 2175 Å bump again. Similarly to the static-no-scattering model, the bump strength weakens at higher inclination angles, but the bump strength increases over time after t = 3.0 Gyr at higher inclination angles. This is because the observed stellar emission at higher inclination angles is more affected by the dust grains in the outer disk of the galaxy, where small grains are more produced by shattering (Matsumoto et al. 2024).

When scattering is included (lower panels of Fig. 7), the 2175 Å bump strength increases in both the dynamic-with-scattering and static-with-scattering models. The 2175 Å attenuation feature is primarily caused by small grains with near-zero albedo, including PAHs, whereas the surrounding wavelengths are attenuated by somewhat larger grains with nonzero albedo (see also Appendix B). As a result, scattering is more efficient at wavelengths near 2175 Å than at the feature itself. This leads to enhanced scattering out of the line of sight in the adjacent continuum, suppressing the surrounding attenuation and increasing the contrast of the 2175 Å bump in the attenuation curve, making it appear more pronounced.

The evolutionary trend of the 2175 Å bump strength for the dynamic-with-scattering model follows that of the extinction curve, even when scattering is taken into account. Therefore, we conclude that, while the 2175 Å bump strength is influenced by scattering and the star-dust geometry, the evolutionary trend of the 2175 Å bump strength still traces the change in the grain size distribution. Remarkably, the attenuation curves in the dynamic-with-scattering dust model show a shallower or even absent 2175 Å bump in the galaxy at the early evolutionary phase of the galaxy (t < 250 Myr), as the interplay of shattering and accretion has not yet fully been activated in the ISM. Therefore, the 2175 Å bump strength can serve as a good proxy for the mass fraction of small carbonaceous grains relative to the total dust mass in galaxies, at a given inclination angle.

Fig. 8 Upper panels: relationships between the slope and AV for the static-with-scattering and dynamic-with-scattering models (left and right, respectively) at different inclination angles. Diamonds indicate simulation snapshots at t = 0.05, 0.25, 0.5, 1, 3, and 10 Gyr, with sizes increasing with time. Solid lines connect diamonds at each inclination angle, representing evolutionary tracks. Gray and purple contours illustrate 25th, 50th, 75th, 95th, and 99th percentiles for the galaxies from the GSWLC observational catalog (Salim et al. 2018) and the Large Early Galaxy Astrophysics Census survey (Barišić et al. 2020), respectively. The blue symbols indicates the observational measurements made for high-redshift galaxies from the REBELS ALMA large program (Fisher et al. 2025). Lower panels: relationships between the bump and AV for static-with-scattering and dynamic-with-scattering models (left and right panels, respectively).

4.3 A_V versus slope and bump

The upper panels of Fig. 8 show the relationship between the slope and A_V at different inclination angles for the static-with-scattering and dynamic-with-scattering models. The gray and purple contours in each panel of Fig. 8 show the relationships between slope and A_V obtained from observations of local galaxies with stellar masses of M_∗ = 10⁹–10¹¹ M_⊙ from the GALEX SDSS-WISE Legacy Catalog (GSWLC; Salim et al. 2018) and star-forming galaxies at higher redshifts of z ∼ 0.8 from the Large Early Galaxy Astrophysics Census survey (LEGA-C survey; Barišić et al. 2020), respectively.

The relationships between the slope and A_V for both static-with-scattering and dynamic-with-scattering models show a similar trend to those of the observations. Overall, the slope becomes steeper at lower A_V. Indeed, as the line-of-sight optical depth at optical wavelengths decreases, optical photons easily escape along the line of sight from the galactic plane, decreasing A_V and steepening the slope, especially at lower inclination angles. These results confirm that the decreasing trend of the slope against A_V is explained by the effect of scattering (Chevallard et al. 2013, we further discuss in Sect. 5.2) and is consistent with the observational results from Salim et al. (2018) as well as the ones from simulations in a static dust model framework and at lower resolution presented by Trayford et al. (2020) and Sommovigo et al. (2025). This relationship is seen not only in nearby galaxies (Salim et al. 2018; Leja et al. 2017) but also in higher redshift galaxies (Barišić et al. 2020; Boquien et al. 2022). Moreover, compared to the static-with-scattering model, the dynamic-with-scattering model produces flatter attenuation curves at lower A_V, as observed in high-redshift galaxies (Fisher et al. 2025), since small grains are absent in the early evolutionary stages of galaxies.

The lower panels of Fig. 8 show the relationship between the 2175 Å bump strength and A_V at different inclination angles for the static-with-scattering and dynamic-with-scattering models.

A larger fraction of the dust becomes optically thick to stellar emission at UV wavelengths at higher A_V, resulting in a stronger bias toward emission from unobscured stars in the observed UV light (Narayanan et al. 2018), thereby reducing the apparent bump strength. In the dynamic-with-scattering model, the bump strength is initially low and increases over time as small dust grains form within the galaxy. However, its peak strength diminishes with increasing A_V or inclination angle, similar to the behavior seen in the static-with-scattering model. While the static-with-scattering model exhibits a tight correlation between the bump strength and A_V, the dynamic-with-scattering model predicts a broader distribution of bump strengths and A_V, varying with time and inclination angle, and encompassing the observational results of Barišić et al. (2020). In contrast, no clear trend is observed in the data from Salim & Boquien (2019), likely due to the 2715 Å bump falling between the GALEX/FUV and NUV bands in photometric observations of local galaxies, which complicates its accurate estimation (see Appendix A.1 of Salim et al. 2018). In the dynamic-with-scattering model, young galaxies with low A_V exhibit shallower or absent bumps in their attenuation curves, regardless of the inclination angle. This finding is consistent with recent observations of galaxies at redshifts z > 6 (Fisher et al. 2025), which also show the absence or shallowness of these bumps in galaxies with low A_V. This could suggest a lack of small grains or a small amount of them in the early evolutionary stages of galaxies.

Fig. 9 (a) Schematic picture of a toy model with a single star embedded in a dust disk. τLOS and τ⊥LOS are optical depths at V-band (λ = 5500 Å) along the line of sight and perpendicular to the line of sight. The two parameters that serve to vary the dust geometry. The resulting attenuation curves from this model are influenced solely by scattering processes, which vary in response to changes in the star–dust geometry. (b) Relationship between AV and slope based on this simple toy model. Different line styles represent different τ⊥LOS.

5 Discussion

5.1 Understanding the 3D effect of scattering with a toy model

In Sect. 4.3 we show that the attenuation curve becomes steeper over time in our MW-like galaxy simulation due to the increasing possibility of scattering as the dust mass increases. This behavior reflects the 3D radiative effect of scattering, which depends not only on the optical depth along the line of sight but also on the optical depth perpendicular to it. To demonstrate the effect of scattering on the slope of the attenuation curve, we consider a toy model with a single star embedded in a dust disk as shown in Fig. 9a. Here, dust properties are taken from Draine & Li (2007), and the geometry of the dust disk is parameterized by the optical depths at the V-band (5500 Å) along the line of sight and perpendicular to the line of sight (τ_LOS and τ_⊥LOS, respectively). We then conduct 3D radiative transfer calculations using SKIRT for this toy model to calculate the attenuation curve for the star. In this model, only scattering influences the attenuation curves.

Figure 9b shows the relationship between A_V and the slope of the attenuation curves based on the radiative transfer calculations for the simple toy model. The slope decreases with increasing A_V because, at lower τ_LOS, scattered photons from the central star can more easily escape along the line of sight. This causes the observed A_V to be lower than the actual optical depth τ_LOS, resulting in a steeper attenuation curve. This effect is further enhanced at higher optical depths perpendicular to the line of sight (τ_⊥LOS), where scattering is more likely to happen in the dusty disk. As a result, the observed A_V becomes lower than the actual τ_LOS, further steepening the attenuation curve.

Therefore, the relation between A_V and slope and the associated scatter arises from the configuration of dust optical depths along and perpendicular to the line of sight. These two optical depths are ultimately governed by the inclination angle of galaxies and the star-dust geometry.

5.2 Physical mechanisms shaping the attenuation curve

Figure 10a shows the relation between the slope and A_V for the dynamic-with-scattering and dynamic-no-scattering models. The relation for the dynamic-no-scattering exhibits a flat trend, while the dynamic-with-scattering model shows a decreasing trend of the slope with increasing A_V due to scattering, as discussed in Sect. 4.3. For further discussion for the dynamic-with-scattering model, the key processes influencing the slope are schematically illustrated in Fig. 10b. At lower inclination angles for the dynamic-with-scattering model, optical photons escape along the line of sight via scattering, which reduces A_V and steepens the attenuation curve slope. Up to t = 1 Gyr the slope decreases with increasing A_V,driven by the growing optical depth along the line of sight (τ_LOS). After t = 1 Gyr, the trend reverses: the slope increases with decreasing A_V. This is because the optical depth perpendicular to the line of sight increases (τ_⊥LOS), enhancing the probability of scattering within the galactic plane. We confirm that the decrease in A_V during this later phase is due to scattering effects rather than a decline in dust column density along the line of sight over time (see Appendix C). Thus, changes in the star–dust geometry modulate the efficiency and directionality of scattering, producing time-dependent variations in both the slope and A_V, as also demonstrated in Sect. 5.1 with a toy model. As a result, at lower inclination angles, galaxies with higher dust masses at a given A_V tend to exhibit steeper slopes in our simulation. We also note that during this period large grains (a > 0.1 µm) begin to form through coagulation, enhancing the probability of scattering in the V-band. However, this trend can be significantly degenerate with that produced by large variations in the star-dust geometry (Narayanan et al. 2018; Trayford et al. 2020; Vijayan et al. 2024; Sommovigo et al. 2025).

At higher inclination angles, scattering becomes less effective, and the attenuation curve slope is instead shaped by changes in the grain size distribution and star-dust geometry. Consequently, in Fig. 10a, the evolutionary trends of the slopes for the dynamic-with-scattering model at higher inclination angles align with those for the dynamic-no-scattering model. Up to t = 3 Gyr, the slope steepens slightly with increasing A_V, as young stars are preferentially obscured by dust and small grains form through shattering within the galaxies. During this phase, the observed V-band flux is dominated by older stars, while the intrinsic flux includes contributions from both young and old stars. Also, at earlier times (t ≤ 0.5 Gyr), the slope remains relatively flat due to the lack of small grains in the galaxy, which reduces A_FUV. After t = 3 Gyr, the star formation rate declines significantly, and the intrinsic V band emission becomes increasingly dominated by older stars. As a result, the intrinsic and observed light originate from similar stellar populations, reducing the difference between them, which lowers A_V and steepens the slope. After t = 7 Gyr, the intrinsic FUV band emission is also increasingly dominated by older stars, resulting in lower A_FUV and slightly flattening the slope. Overall, the relative importance of scattering, star–dust geometry, and grain size distribution varies nontrivially depending on A_V and inclination angles.

Figure 10c shows the relationship between the 2175 Å bump strength and A_FUV for the dynamic-with-scattering and dynamic-no-scattering models. We use A_FUV instead of A_V because it is less affected by scattering. For further discussion, the key processes affecting the bump strength are schematically illustrated in Fig. 10d. In Fig. 10c, the 2175 Å bump strength in the dynamic-with-scattering model is slightly higher than in the dynamic-no-scattering model, particularly at lower inclination angles. This is because scattering is more efficient at wavelengths adjacent to 2175 Å than at the bump itself (see Appendix B). As this effect is more pronounced at lower A_FUV or lower inclination angles, the contrast of the bump in the attenuation curve is enhanced. At low inclination angles for both dynamic-with-scattering and dynamic-no-scattering models, the bump strengths dramatically increase over time up to t = 0.5 Gyr, as small grains (a < 0.05 µm) form in the galaxy through the interplay between accretion and shattering. After t = 0.5 Gyr, small grain formation becomes less efficient, since some small grains start being coagulated into larger grains, causing the bump strength to nearly plateau. After t = 1.0 Gyr, coagulation slightly weakens the bump strength with time. At the same time, as dust column density increases with time, the UV optical depth exceeds unity (A_2175Å > 1.0) along some sight lines. This leads to a stronger bias toward emission from unobscured stars in the observed UV light, which further reduces the apparent bump strength (Narayanan et al. 2018). After t = 5 Gyr, very small carbonaceous grains such as PAH (a < 0.0012 µm) begin to form through shattering, leading to a renewed increase in bump strength.

At higher inclination angles, the bias toward emission from unobscured stars becomes more pronounced, systematically reducing the observed bump strength. Early on (up to t = 0.5 Gyr), the bump strength shows a slight increase due to the formation of small carbonaceous grains. Afterward, as the dust column density grows, this leads to higher UV optical depths and an enhanced bias toward unobscured stellar emission. After t = 3 Gyr, the bump strength significantly increases again. This is because small carbonaceous grains preferentially form in the diffuse ISM, which is more likely to be along the line of sight at higher inclination angles. After t = 7 Gyr, the intrinsic FUV emission becomes increasingly dominated by older stars, leading to a reduction in A_FUV. Together, these effects highlight the critical roles of grain size evolution and star–dust geometry in shaping the strength and visibility of the 2175 Å bump over time.

Fig. 10 (a) Relation between the slope and AV. Colored diamonds represent the dynamic-with-scattering model, while gray diamonds shows the dynamic-no-scattering model. The color indicates the total dust mass within the galaxy for the dynamic dust model, which depends on simulation time as indicated by the legend. Here, the diamonds correspond to simulation snapshots at t = 0.05, 0.25, 0.5, 1, 3, and 10 Gyr. For the dynamic-with-scattering and no-scattering models, solid and dashed lines connect the diamonds, respectively, representing the evolutionary tracks at each inclination angle. (b) Physical explanation of the relation between the slope and AV. Solid arrows indicate the effects shaping the evolution of the attenuation curve slope, while the dashed arrow highlights the dependence of the slope on the inclination angle. The dark gray-shaded region corresponds to the early phase of our MW-like galaxy (t < 1 Gyr), whereas the light gray-shaded region marks the later phase (t ≥ 1 Gyr). (c) Relation between the 2175 Å bump strength and AFUV. The color indicates the mass fraction of small grains (a < 0.05 µm) relative to the total dust mass within the galaxy for the dynamic dust model, which depends on simulation time as indicated by the legend. Here, the diamonds correspond to simulation snapshots at t = 0.05, 0.25, 0.5, 1, 3, and 10 Gyr. (d) Physical explanation of the relation between the 2175 Å bump strength and AFUV. Thin solid arrows indicate the effects shaping the evolution of the 2175 Å bump, while thick solid arrows represent effects that systematically change the trend across inclination angles or AFUV. The dashed arrow highlights the dependence of the slope on the inclination angle.

5.3 Comparison with theoretical studies based on large-volume cosmological simulations

In this study, we demonstrated a strong relation between A_V and the slope of the attenuation curves. The scatter in this relation is attributed to variations in inclination angles and in the configuration of dust optical depths around stars. (see Figs. 8 and 10a). However, we handle only a MW-like galaxy simulation and lack the variation of star-dust geometry. Previous studies based on large-volume cosmological simulations offer an advantage in capturing the variation of star-dust geometry, although they typically assume a fixed dust-to-metal ratio and constant dust properties based on dust models (e.g., Weingartner & Draine 2001; Zubko et al. 2004; Draine & Li 2007), and are limited by lower spatial resolution. Trayford et al. (2020) identified a correlation between A_V and the UV slope in EAGLE galaxies across redshifts z = 0–2, with variations attributed to galaxy inclination. Specifically, galaxies viewed at higher inclination angles tend to exhibit flatter attenuation curves. This trend is consistent with our findings, as well as with results from local galaxies in TNG simulations (Sommovigo et al. 2025) and high-redshift galaxies in the FirstLight and MUFASA simulations (Mushtaq et al. 2023; Narayanan et al. 2018, respectively).

Furthermore, the importance of the star-geometry is discussed by examining the factors to change A_V. Trayford et al. (2020) found that A_V is predominantly influenced by dust surface densities of galaxies. Mushtaq et al. (2023) found that galaxies with higher stellar masses contain a larger amount of dust. Concurrently, Sommovigo et al. (2025) revealed that galaxies with higher Σ_SFR (i.e., compact star-forming galaxies) tend to show higher A_V. These results suggest that A_V is mainly driven by the star-dust geometries of galaxies. As a result, the star-dust geometry changes the slope of the attenuation curves according to the relation between A_V and slope. Narayanan et al. (2018) provided physical insights into variations of the slope of attenuation curves. The relative fraction of dust-obscured young stars compared to nonobscured old stars changes the shape of observed spectra of galaxies at UV wavelengths: galaxies with heavily obscured young stars exhibit steeper curves, as the observed UV light is dominated by older stellar populations. In our simulation, we confirm that changes in the dust obscuration fraction of young stars affect the slope of the attenuation curves in the static-no-scattering model (see the upper left panel in Fig. 6). This effect, however, becomes less pronounced when scattering is included, as scattering significantly changes the slope of the attenuation curve. We note that this lower responsiveness of the pure star–dust geometry effect might be due to the limited variation in star–dust geometry in our simulations. Therefore, a more comprehensive investigation using a large suite of galaxy simulations is required to fully assess the impact of star–dust geometry on attenuation curves along with dust evolution.

Moreover, we found that a decreasing trend of the 2175 Å bump strength with A_V in the static-with-scattering model (see the lower left panel in Fig. 8) is consistent with Narayanan et al. (2018) and Sommovigo et al. (2025). Narayanan et al. (2018) found a correlation between the fraction of unobscured young stars in galaxies and the bump strengths, suggesting the importance of the effect of the star-dust geometry. Specifically, Seon & Draine (2016) argued that an inhomogeneous distribution of optically thick clumps obscuring stars weakens the 2175 Å bump. Additionally, Narayanan et al. (2018) indicated that the bump strength slightly decreases when scattering is less effective in galaxies and concluded that scattering is a minor effect that influences the bump strength. These findings agree with our results.

5.4 Comparison with observational studies

The relation between A_V and the slope of the attenuation curve is also reported in observational studies (e.g., Salmon et al. 2016; Leja et al. 2017; Salim et al. 2018; Boquien et al. 2022). We argued that the scatter in the relation arises from variations in galaxy inclination: the attenuation curves at higher inclination angles tend to be flatter and exhibit less pronounced 2175 Å bumps. Observationally, Battisti et al. (2017a) found flatter curves and weaker bumps at higher inclination angles in local star-forming galaxies, consistent with our results (see also Wild et al. 2011). This trend is also confirmed at higher redshifts (z ∼ 0.8; Barišić et al. 2020).

In the lower right panel of Fig. 8, the comparison of the relation between A_V and 2175 Å bump with observations (Battisti et al. 2017a) shows that the dynamic-with-scattering model reproduces the observed decreasing trend of the maximum bump strength with A_V, which is also found in high redshift observed galaxies (Markov et al. 2025b). We attributed this trend to a stronger bias toward unobscured stellar emission in the observed UV light at higher UV optical depths (Narayanan et al. 2018). We also argued that the bump strength increases according to the formation of small grains within galaxies via the interplay between shattering and accretion, reproducing the observed distribution of the bump strengths and A_V (Barišić et al. 2020).

Recent JWST observations have revealed the characteristics of the attenuation curves in galaxies at high redshifts (z > 4). Markov et al. (2025a) show the average attenuation curves steepening with cosmic time, indicating that galaxies at higher redshifts of z = 4.5-11.5 exhibit flatter attenuation curves with less pronounced 2175 Å bumps. Additionally, Fisher et al. (2025) exhibits flat attenuation curves from 12 Lyman-break galaxies at z = 6.5–7.7 with weak 2175 Å bump. Notably, high-redshift galaxies tend to exhibit flatter curves even at low A_V compared to nearby galaxies (see also Markov et al. 2025b). Markov et al. (2025a) suggested that the flat attenuation curves observed at high redshifts are driven by large grains produced by supernovae (McKinney et al. 2025), which have low UV opacity. This scenario also provides an explanation for the weaker 2175 Å bump seen in the average attenuation curves at high redshifts. In the dynamic-with-scattering model, we find that the attenuation curves exhibit shallower slopes and weaker 2175 Å bumps during the early stages of galaxy evolution (see the lower right panels in Figs. 6 and 7). In our simulations, attenuation curves at high inclination angles tend to be flatter during the early epoch (t < 1 Gyr), due to the lack of small grains. At low inclination angles, attenuation curves are more strongly influenced by scattering and thus appear steeper. In the early stages of galaxy evolution, however, the effect of scattering is reduced. This is because the dust distribution remains compact in three dimensions at early times (t < 1 Gyr) – unlike the extended dusty disk that forms later (see Fig. D.1) – resulting in weaker scattering and correspondingly flatter attenuation curves.

On the other hand, our dynamic-with-scattering model initially shows a weak 2175 Å bump, which becomes more prominent on a timescale of ∼250 Myr as small grains (<0.05 µm) are produced through the interplay between accretion and shattering, consistent with the scenario proposed by Markov et al. (2025a). Ormerod et al. (2025) reported the detection of a 2175 Å bump in a Lyman-break galaxy at z = 7.1, indicating the rapid formation of small grains. They proposed the presence of an older stellar population (∼252 Myr), suggesting that small grains could form through shattering processes following initial star formation. In contrast, Witstok et al. (2023) identified a strong bump in a galaxy at z = 6.7 without signs of an old stellar population, implying that small grains may form shortly after the onset of star formation, potentially through mechanisms such as carbonaceous grain production in Wolf–Rayet stars (Lau et al. 2023, 2024).

5.5 Limitations of our simulation

First, we use an isolated galaxy simulation with an initial stellar disk and bulge, and a dust-free gas disk. Without cosmological effects like mergers and with steady IGM inflow, the galaxy follows a moderate star formation history and weak stellar feedback. Thus, our model does not fully capture the diversity of star-dust geometry seen across galaxy populations (see Sect. 5.3).

Second, our simulation does not have sufficient spatial resolution to resolve the star formation processes embedded in individual molecular clouds (MCs), which require subparsec resolution. Di Mascia et al. (2025) discussed the significant impact of dust evolution in MCs on the attenuation curves of these regions. Additionally, many observations report that the optical reddening, E(B−V), for nebular emission lines from star-forming regions is larger than that for optical stellar continuum emission (e.g., Calzetti et al. 1994, 2000; Buat et al. 2018; Reddy et al. 2015, 2020; Shivaei et al. 2020), implying the importance of modeling the steep attenuation curve in birth clouds (Charlot & Fall 2000). On the other hand, Trayford et al. (2020) demonstrated that modeling the steep attenuation curve in birth clouds does not significantly affect galaxy attenuation curves in EAGLE simulations.

Finally, our dust evolution model also contains some uncertainties. In our previous studies, Matsumoto et al. (2024) compared the dust properties from our simulations with observations of nearby galaxies and found that the model underestimates the abundance of very small grains such as PAHs (a < 0.0013 µm). In contrast, van der Giessen et al. (2024) reported that our model overestimates the amount of small grains with sizes a < 0.015 µm. These findings suggest that the Draine & Li (2007)-like grain size distribution for carbonaceous grains characterized by distinct peaks at both small and large sizes (a = 0.0006 and 0.4 µm, respectively) may be more appropriate. Furthermore, our post-processing dust decomposition recipe does not account for variations in dust accretion efficiency depending on the type of metal elements, nor for the depletion of key elements such as O, Si, Fe, and C within galaxies (see Choban et al. 2022; Dubois et al. 2024). Li et al. (2021) found that the graphite-to-silicate ratio of dust grains influences the strength of the 2175 Å bump in the extinction curve of a MW-like galaxy simulation. Dubois et al. (2024) reported that the accretion efficiency of carbonaceous grains is lower than that of silicate grains, implying that the emergence of the prominent 2175 Å bump, which appears on a timescale of ∼250 Myr in our model, could be delayed in their model. In addition, the extinction curve in our simulation at t = 10 Gyr is steeper than that for the MW observation (Pei 1992, see Fig. 3). This is because our model does not account for the depletion of key metal elements, potentially leading to an overestimation of the amount of silicate grains. Therefore, future simulations should incorporate dust composition evolution for a more accurate treatment.

6 Conclusions

In this paper, we investigate the average extinction curve and global attenuation curves of a MW-like galaxy simulation and examined the physical mechanisms that shape the attenuation curves. To accomplish this, we performed post-processing dust radiative transfer calculations using the SKIRT code, based on the MW-like galaxy simulation (Matsumoto et al. 2024). Our model incorporates the evolution of grain size distributions across 30 sizes from 3.0 × 10⁻⁴ to 10 µm and considers the hydrodynamic evolution of the galaxy, including star formation and the OSAKA feedback model. For modeling the dust composition, we decomposed the dust grains into silicate and carbonaceous components, as well as PAHs, based on the abundance of Si and C and the mass fraction of dense gas in each gas particle in post-processing. To disentangle the impacts of star-dust geometry, scattering, and grain size distribution on the attenuation curves, we generated four different types of attenuation curves (static-no-scattering, static-with-scattering, dynamic-no-scattering, and dynamic-with-scattering). The details of the model prescription are explained in Sect. 4. We summarize our results as follows:

1. The extinction curve evolves from flat to steep by t = 0.5 Gyr, gradually developing the 2175 Å bump, driven by the formation of small grains via shattering and accretion. It then flattens due to coagulation up to t = 3 Gyr. Later, the extinction curve slightly steepens again up to t = 10 Gyr, while the 2175 Å bump grows as PAHs form in the ISM (see Fig. 3);

2. At t = 0.1 Gyr, the extinction curve is dominated by silicate grains across all wavelengths. However, at later times, carbonaceous grains become the dominant contributors to extinction at wavelengths longer than 2000 Å (see Fig. 4);

3. At low inclination angles, the attenuation curve for the dynamic-with-scattering model evolves differently from the extinction curve because of scattering. The attenuation curves slightly flatten over time up to t = 1 Gyr due to the increasing dust column density along the line of sight. Later, they become progressively steeper due to the interplay between scattering and extending dust disk geometry (see Fig. 6);

4. At high inclination angles, the attenuation curves starts with lower values at higher inclination angles due to the lack of small grains at t < 0.5 Gyr. At later times, the attenuation curves become slightly steeper due to differences in the contributions of young and old stars to the observed and intrinsic spectra and the formation of small grains;

5. Comparing the attenuation curves between the dynamic-with-scattering and static-with-scattering models shows that the underlying extinction curve influences attenuation curves at high inclination angles and during the early epoch (t < 1 Gyr). The attenuation curves for the dynamic-with-scattering is significantly flatter than static-with-scattering due to the lack of small grains. At low inclination angles, variations in the grain size distribution affect the albedo, thereby influencing the slope primarily through scattering (see also Appendices B and C). These results highlight the importance of dust properties, both in terms of composition and grain size distribution, in determining the shape of attenuation curves;

6. The interplay between scattering and star-dust geometry plays a key role in shaping the relation between the slope and A_V, as seen in observations (Fig. 8; Salim et al. 2018; Barišić et al. 2020). The scatter of this relation arises due to the dust optical depth along and perpendicular to the line of sight; thus, it reflects the inclination angle and star-dust geometry (see Figs. 9, 10a, and 10b). Changes in the grain size distribution and the star-dust geometry (Narayanan et al. 2018) also play a key role in shaping slope of the attenuation curves at high inclination angles or high A_V (> 0.5);

7. At lower inclination angles, the 2175 Å bump becomes stronger on a timescale of ∼250 Myr due to the formation of small carbonaceous grains (a < 0.05 µm) through the interplay between shattering and accretion (see Fig. 7).

   After t = 1.0 Gyr, the bump strength become lower as small grains are coagulated. At the same time, as the dust column density along the line of sight increases with time, the bump strength becomes weaker due to the stronger bias toward emission from unobscured stars; this effect is governed by star-dust geometry. After t = 5 Gyr, very small carbonaceous grains such as PAHs (a < 0.0012 µm) begin to form through shattering, enhancing bump strength. At higher inclination angles, the bias toward emission from unobscured stars becomes more pronounced, further reducing the observed bump strength;

8. Variations in the bump strength are mainly caused by the variations in the grain size distributions and the star-dust geometry (see Figs. 10c and 10d). The maximum bump strength become systematically weaker at higher A_FUV or higher inclination angles because of the bias toward emission from unobscured stars.

These results provide physical insights into recent JWST observations, which reveal flatter attenuation curves and less prominent 2175 Å bumps in galaxies at high redshifts (z > 4; Markov et al. 2023, 2025a,b; Fisher et al. 2025). However, the systematic evolution of the attenuation curves could vary depending on the diversity of galaxies, especially star-dust geometry. Therefore, to better understand how dust evolution impacts the systematic evolution of attenuation curves, future studies using cosmological simulations that incorporate dust evolution across a wide range of galaxy types will be crucial.

Data availability

The data used in this paper are publicly available⁵.

Acknowledgements

We thank Dr. Hirashita and Prof. Nakagawa for their valuable discussions. We gratefully acknowledge Samir Salim for sharing data from his study. We are grateful to Lucy Reading-Ikkanda for her invaluable help in creating and refining the schematic figures for this paper. KM acknowledges the financial support from the Flemish Fund for Scientific Research (FWO-Vlaanderen) provided through Grant number 1169822N and 12AJQ26N. Numerical computations were performed on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, on the SQUID at the Cybermedia Center, Osaka University, as part of the HPCI system Research Project (hp220044, hp230089, hp240141), and on the Flatiron Institute’s research computing facilities (Rusty). This work is supported in part by the MEXT/JSPS KAKENHI grant numbers 20H00180, 22K21349, 24H00002, 24H00241 (KN). KN acknowledges the travel support from the Kavli IPMU, World Premier Research Center Initiative (WPI), where a part of this work was conducted.

References

Appendix A Global properties of the MW-like galaxy simulation

Figure A.1 shows the global evolution of the galactic properties of the MW-like galaxy simulation. The simulation shows a high initial star formation rate, which gradually declines over time, as shown in the top panel. The middle panel presents the evolution of the total dust mass and the dust masses across different size ranges. The total dust mass increases over time mainly due to the accretion on dust grains. Initially, dust is dominated by large dust grains (a ≥ 0.05 µm) due to the dust production of SNe II, while the small grains (a < 0.05 µm) are gradually produced by the interplay of shattering and accretion over time. The very small dust grains (a < 0.002 µm) including PAHs, form later compared to the small dust grains. The bottom panel shows the mass evolution of silicate, graphite, and PAHs within the MW-like galaxy simulation. Initially, silicate grains produced by SNe II dominate the dust mass in the galaxy. At later times, SNe Ia contribute additional gas-phase carbon, increasing the abundance of carbonaceous grains (Saitoh 2017). After t = 0.5 Gyr, PAHs are produced through shattering, resulting in a dramatic increase in their mass fraction. We note that our simulation does not take into account condensation efficiencies or accretion efficiencies that depend on the metal species. These factors could potentially affect the mass evolution of the dust compositions (Li et al. 2021; Choban et al. 2022; Dubois et al. 2024).

Fig. A.1 Top panel: Star formation history of the MW-like galaxy simulation. Middle panel: Evolution of the total dust mass over time, along with the dust masses of small grains (a < 0.05 µm) and very small grains (a < 0.002 µm). Bottom panel: Time evolution of the dust masses for silicate, graphite, and PAHs.

Appendix B Evolution of the average albedo of the MW-like galaxy simulation

Figure B.1 shows the evolution of the wavelength dependency of the average albedo of the MW-like galaxy simulation. Overall, the average albedo at UV wavelengths is lower compared to optical wavelengths, since small grains with lower albedo are responsible for the opacity at UV wavelengths. At t = 0.10 Gyr, large grains dominate in the galaxy, resulting in a high albedo across all wavelengths. Later, small grains form through shattering and accretion, and albedo gradually becomes lower up to t = 0.5 Gyr. At the same time, the albedo at 2175 Å becomes lower compared to other wavelengths. After t = 0.5 Gyr, the albedo at optical wavelengths becomes higher due to the formation of larger grains via coagulation, while the albedo at UV wavelengths continuously goes down due to the formation of very small grains including PAH through shattering. At t = 10 Gyr, the albedo at all wavelengths is lower than that of Draine & Li (2007) model.

Fig. B.1 Evolution of the wavelength dependency of the average albedo of our MW-like galaxy simulation. The gray dashed line represent the wavelength dependency of the albedo of the Draine & Li (2007) model.

Appendix C Evolution of the global dust attenuation of the MW-like galaxy simulation

The upper left panel of Fig. C.1 shows the time evolution of the global A_V for the static dust model, comparing results that include scattering (with-scattering) and those that neglect it (no-scattering). At all inclination angles, A_V increases over time as the total dust mass in the galaxy grows. At an inclination angle of 90^◦, however, A_V shows a slight decline after t = 3 Gyr for both the static-with-scattering and no-scattering models. This decline is attributed to the decreasing fraction of young stars over time (Narayanan et al. 2018). In contrast, at an inclination angle of 0^◦, A_V for the static-with-scattering model begins to decrease after t = 2 Gyr, whereas A_V for the static-no-scattering model continues to increase. This different behavior is driven by the effect of scattering. The lower left panel of Fig. C.1 presents the time evolution of the ratio between A_V in the static-with-scattering and static-no-scattering models, providing a quantitative measure of scattering effects on A_V. At all inclination angles, this ratio increases over time, indicating that the impact of scattering is gradually reduced as the optical depth along the line of sight becomes larger. At an inclination angle of 90^◦, the ratio is close to unity, suggesting that scattering has a minimal effect. In contrast, at lower inclination angles, the ratios are smaller, meaning that scattering reduces the observed A_V. Moreover, after t = 1 Gyr, the ratio declines further over time, implying that the contribution from scattering becomes increasingly significant because the dust disk geometry becomes more extended and scattering events become more likely.

Fig. C.1 Time evolution of the global AV for the static and dynamic dust models in our Milky Way-like galaxy simulation is shown in the left and right panels, respectively. In each panel, solid lines indicate attenuation including scattering (with-scattering), while dashed lines represent attenuation neglecting scattering (no-scattering). The inclination angle is color-coded. The bottom panels show the ratio between the with-scattering and no-scattering attenuation. This ratio quantifies the impact of scattering on the total attenuation, with lower values indicating a stronger contribution from scattering.

The upper right panel of Fig. C.1 shows the time evolution of the global A_V for the dynamic dust model. The overall trends of A_V for both the dynamic-with-scattering and no-scattering models are similar to those for the static models. However, at lower inclination angles, A_V for the dynamic-with-scattering model is higher than that of the static-with-scattering model. This is because large silicate grains (a > 0.05 µm) are more abundant in the static model (Draine & Li 2007), leading to enhanced scattering and consequently lower A_V. The lower right panel of Fig. C.1 shows the time evolution of the ratio between A_V in the dynamic-with-scattering and no-scattering models. At all inclination angles, the ratios for the dynamic model increase more significantly than those for the static model up to t = 0.5 Gyr, driven by the formation of small grains with low albedo. At lower inclination angles, the dynamic model exhibits higher ratios than the static model, indicating that scattering is less effective in reducing A_V. While the ratio remains nearly constant at an inclination angle of 90^◦ after t = 1 Gyr, it declines at lower inclination angles. This trend is driven by the growth of large grains through coagulation in the ISM, which enhances scattering and further reduces A_V.

The upper left panel of Fig. C.2 shows the time evolution of the global A_FUV for the static-with-scattering and no-scattering models. At all inclination angles, A_FUV for both the static-with-scattering and no-scattering models increases over time as the total dust mass in the galaxy grows up to t = 7 Gyr. At later times, A_FUV declines because the intrinsic UV spectrum becomes dominated by old stars, reducing the contrast between the intrinsic and observed spectra. The lower left panel of Fig. C.2 presents the time evolution of the ratio between A_FUV in the static-with-scattering and static-no-scattering models. Overall, the ratio at the FUV wavelength is larger than that at the V-band wavelength, indicating that scattering has a less significant effect at the FUV wavelength. This is primarily because the opacity at the FUV wavelength is generally higher than at the V-band wavelength, causing most FUV photons to be eventually absorbed. The ratio at FUV wavelengths increases over time until t = 1.0 Gyr, since the effect of scattering is reduced as the optical depth along the line of sight becomes larger. After t = 1 Gyr, the ratio at the FUV wavelength remains constant.

The upper right panel of Fig. C.2 shows the time evolution of the global A_FUV for the dynamic dust model. The overall trends of A_FUV for both the dynamic-with-scattering and no-scattering models are similar to those for the static models. However, A_FUV for the dynamic model starts with a lower value compared to the static models, since only large grains, which have low opacity at FUV wavelengths, are present at the beginning of the simulation. As small grains form through accretion and shattering over time, A_FUV for the dynamic model becomes comparable to that for the static model. The lower right panel of Fig. C.2 shows the time evolution of the ratio between A_FUV in the dynamic-with-scattering and no-scattering models. At all inclination angles, the ratios for the dynamic model increase more significantly than those for the static model up to t = 0.5 Gyr, driven by the formation of small grains with low albedo at the FUV wavelength. After t = 1.0 Gyr, the ratio at the FUV wavelength remains constant.

Fig. C.2 Same as Fig. C.1, but time evolution of the AFUV.

In summary, A_V is significantly influenced by scattering, particularly at lower inclination angles. The scattering processes are sensitive to changes in the star–dust geometry, such as inclination angle and the configuration of dust optical depths, as well as to variations in albedo driven by the evolving grain size distribution. In contrast, A_FUV is less influenced by scattering due to its generally higher optical depths. Instead, it is more sensitive to the abundance of small grains and the star–dust geometry.

Appendix D Maps of attenuation at various wavelengths

The top panels of Fig. D.1 show the time evolution of the A_V map viewed face-on. Initially, dust is produced by SNe II primarily in the central region of the galaxy. As time progresses, dust forms via accretion, extending outward along the spiral arms from the center to the outskirts. Large grains are responsible for A_V, and thus, the A_V maps exhibit thin spiral arms tracing the dense gas in the ISM. The second panels from the top in Fig. D.1 show the time evolution of the attenuation maps at λ = 2175 Å (A_2175Å), where the attenuation is primarily caused by small grains. Small grains form via shattering, especially in the diffuse regions surrounding the spiral arms. As a result, the A_2175Å maps exhibit thicker spiral arms. The third panels from the top in Fig. D.1 present the time evolution of the A_FUV map. Overall, its morphology is similar to that of the A_2175Å maps. At FUV wavelengths, scattering can be significant in regions with low dust surface density, allowing scattered photons to escape from the regions. As a result, A_FUV becomes negative in such regions, particularly where supernova feedback sweeps out gas and dust. The bottom panels in Fig. D.1 show the time evolution of the slope maps. Overall, the slopes are flatter in the central regions of the galaxy compared to the outskirts. This is because dust in the center is optically thick to stellar emission at optical-to-UV wavelengths, suppressing scattering. In contrast, dust is optically thin at optical wavelengths in the outer regions, allowing scattered optical photons to escape and resulting in steeper slopes.

Fig. D.1 Time evolution of attenuation maps at V-band, 2175 Å, and FUV-band wavelengths as well as the map of slope (top to bottom panels). The time evolution is exhibited from the left to right panels.

All Tables

All Figures

	Fig. 1 Upper panels: extinction curves for carbonaceous and silicate grains with a given dust grain size (left and right, respectively). The extinction curve is normalized by the opacity at the V-band wavelength (λ = 5500 Å). Lower panels: wavelength dependence of albedo, which is the fraction of scattering against extinction, of carbonaceous and silicate grains (left and right, respectively).
In the text

	Fig. 2 Time evolution of the global grain size distribution in our MW-like galaxy simulation. The solid and dashed gray lines correspond to the grain size distributions of carbonaceous and silicate grains in the static dust model from Draine & Li (2007).
In the text

	Fig. 3 Time evolution of the average extinction curve of the MW-like galaxy simulation indicated by colored solid lines. The extinction curve is estimated by Eq. (2) and normalized by the opacity at the V-band wavelength (λ = 5500 Å). The dashed gray, dotted orange, and dash-dotted red lines represent the extinction curves of the MW, SMC, and LMC, respectively (Pei 1992; Li et al. 2008).
In the text

Fig. 4 Contribution of each dust component to the average extinction curve of the MW-like galaxy simulation at various times (t = 0.1, 0.5, 3.0, and 10.0 Gyr). The black thick line is the total average extinction curve, while the orange solid, green dashed, blue dotted, and light blue dash-dotted thin lines represent the contributions of silicate, graphite, PAH0, and PAH+, respectively. At t = 0.1 Gyr, the extinction curve is dominated by silicate grains across all wavelengths. However, at later times, carbonaceous grains become the dominant contributors to extinction at wavelengths longer than 2000 Å.

In the text

Fig. 5 Time evolution of the global attenuation curve for the MW-like galaxy simulation using the static-with-scattering and dynamic-with-scattering models shown in left and right panels, respectively. The solid and dashed lines represent inclination angles of 0◦ and 90◦ (i.e., face-on and edge-on view, respectively). The attenuation curve is estimated using Eq. (3) and normalized by the attenuation at the V-band wavelength (λ = 5500 Å).

In the text

Fig. 6 Time evolution of the slope for the global attenuation curves from the four different models (dynamic, dynamic-no-scattering, static, and static-no-scattering) of the MW-like galaxy simulation. The inclination angle is color-coded in each panel. Here, all models with Av < 0.01 are excluded. The gray thick lines represent the time evolution of the slope of the average extinction curves for the static and dynamic dust models in the left and right panels, respectively. The value of slope for the static dust model corresponds that of the MW (Pei 1992; Li et al. 2008).

In the text

	Fig. 7 Same as Fig. 6, but considering the time evolution of the 2175 Å bump strength.
In the text

Fig. 8 Upper panels: relationships between the slope and AV for the static-with-scattering and dynamic-with-scattering models (left and right, respectively) at different inclination angles. Diamonds indicate simulation snapshots at t = 0.05, 0.25, 0.5, 1, 3, and 10 Gyr, with sizes increasing with time. Solid lines connect diamonds at each inclination angle, representing evolutionary tracks. Gray and purple contours illustrate 25th, 50th, 75th, 95th, and 99th percentiles for the galaxies from the GSWLC observational catalog (Salim et al. 2018) and the Large Early Galaxy Astrophysics Census survey (Barišić et al. 2020), respectively. The blue symbols indicates the observational measurements made for high-redshift galaxies from the REBELS ALMA large program (Fisher et al. 2025). Lower panels: relationships between the bump and AV for static-with-scattering and dynamic-with-scattering models (left and right panels, respectively).

In the text

Fig. 9 (a) Schematic picture of a toy model with a single star embedded in a dust disk. τLOS and τ⊥LOS are optical depths at V-band (λ = 5500 Å) along the line of sight and perpendicular to the line of sight. The two parameters that serve to vary the dust geometry. The resulting attenuation curves from this model are influenced solely by scattering processes, which vary in response to changes in the star–dust geometry. (b) Relationship between AV and slope based on this simple toy model. Different line styles represent different τ⊥LOS.

In the text

Fig. 10 (a) Relation between the slope and AV. Colored diamonds represent the dynamic-with-scattering model, while gray diamonds shows the dynamic-no-scattering model. The color indicates the total dust mass within the galaxy for the dynamic dust model, which depends on simulation time as indicated by the legend. Here, the diamonds correspond to simulation snapshots at t = 0.05, 0.25, 0.5, 1, 3, and 10 Gyr. For the dynamic-with-scattering and no-scattering models, solid and dashed lines connect the diamonds, respectively, representing the evolutionary tracks at each inclination angle. (b) Physical explanation of the relation between the slope and AV. Solid arrows indicate the effects shaping the evolution of the attenuation curve slope, while the dashed arrow highlights the dependence of the slope on the inclination angle. The dark gray-shaded region corresponds to the early phase of our MW-like galaxy (t < 1 Gyr), whereas the light gray-shaded region marks the later phase (t ≥ 1 Gyr). (c) Relation between the 2175 Å bump strength and AFUV. The color indicates the mass fraction of small grains (a < 0.05 µm) relative to the total dust mass within the galaxy for the dynamic dust model, which depends on simulation time as indicated by the legend. Here, the diamonds correspond to simulation snapshots at t = 0.05, 0.25, 0.5, 1, 3, and 10 Gyr. (d) Physical explanation of the relation between the 2175 Å bump strength and AFUV. Thin solid arrows indicate the effects shaping the evolution of the 2175 Å bump, while thick solid arrows represent effects that systematically change the trend across inclination angles or AFUV. The dashed arrow highlights the dependence of the slope on the inclination angle.

In the text

	Fig. A.1 Top panel: Star formation history of the MW-like galaxy simulation. Middle panel: Evolution of the total dust mass over time, along with the dust masses of small grains (a < 0.05 µm) and very small grains (a < 0.002 µm). Bottom panel: Time evolution of the dust masses for silicate, graphite, and PAHs.
In the text

	Fig. B.1 Evolution of the wavelength dependency of the average albedo of our MW-like galaxy simulation. The gray dashed line represent the wavelength dependency of the albedo of the Draine & Li (2007) model.
In the text

Fig. C.1 Time evolution of the global AV for the static and dynamic dust models in our Milky Way-like galaxy simulation is shown in the left and right panels, respectively. In each panel, solid lines indicate attenuation including scattering (with-scattering), while dashed lines represent attenuation neglecting scattering (no-scattering). The inclination angle is color-coded. The bottom panels show the ratio between the with-scattering and no-scattering attenuation. This ratio quantifies the impact of scattering on the total attenuation, with lower values indicating a stronger contribution from scattering.

In the text

	Fig. C.2 Same as Fig. C.1, but time evolution of the AFUV.
In the text

	Fig. D.1 Time evolution of attenuation maps at V-band, 2175 Å, and FUV-band wavelengths as well as the map of slope (top to bottom panels). The time evolution is exhibited from the left to right panels.
In the text