© The Authors 2024

1. Introduction

The spectral energy distribution (SED) of galaxies characterises the distance and properties of galaxies (e.g., Fioc & Rocca-Volmerange 1997, 2019; Devriendt et al. 1999; Noll et al. 2009; Chevallard & Charlot 2016; Boquien et al. 2019). As dust efficiently absorbs light in the optical or in the ultraviolet (UV), it re-emits this light in the infrared (IR), which is subsequently absorbed and re-emitted in the IR (producing IR multi-scattering). The exact SEDs of galaxies across the entire spectral range are, thus, tied in part to the properties of the dust, and in particular to the size and chemical composition of grains.

Unlike extinction, corresponding to a line-of-sight emitted by a single point source whose light is absorbed and scattered, the attenuation is the result of the collection of multiple sources with different extinctions into a single line-of-sight. While galaxy attenuation curves offer insights into the dust properties (Calzetti et al. 2000), it is crucial to note that they result from the interplay of complex radiative transfer effects within the intricate geometry of the interstellar medium (ISM) and stellar distribution (e.g., Gordon et al. 1997; Witt & Gordon 2000; Granato et al. 2000; Inoue 2005; Narayanan et al. 2018). A key source of knowledge about dust properties comes from the extinction curves of multiple lines-of-sights of individual stars (for which the complexity of dust absorption is better characterised) in the Milky Way or from the Magellanic Clouds. The Milky Way extinction curve is characterised by a moderate extinction in the ultraviolet, as opposed to Magellanic Clouds, and by a significant feature in the optical range, called the 2175 Å bump, very faint in the Large Magellanic Cloud and almost indistinguishable in the Small Magellanic Cloud (Pei 1992). Thanks to the description with Mie theory of light absorption and scattering by solid spherical bodies, and optical properties of grains, it is possible to interpret the shape of the Milky Way extinction curve by a number density distribution skewed towards small grains (a few nm in size) but dominated by large grains (0.1 μm) in mass, as well as a mixed composition of carbonaceous grains and of silicate grains (Weingartner & Draine 2001a). Although, the exact nature of the carbonaceous grains (amorphous, aromatic, aliphatic) and silicate grains (olivine, pyroxene, iron inclusions) and how they mix up in a single grain body is complex (Jones et al. 2017), it is possible to make predictions about the grain properties as a function of the galaxy properties and test them against observables (e.g., Mathis et al. 1977; Draine & Lee 1984; Weingartner & Draine 2001a; Zubko et al. 2004; Draine & Li 2007; Compiègne et al. 2011; Galliano et al. 2018; Hensley & Draine 2023), such as the extinction or attenuation laws, or the depletion of the various elements in the gas phase.

Not only is dust key to reconstruct the properties of galaxies from their SED, but it also intervenes directly or indirectly in the thermodynamics of the ISM. Dust is a catalyst for molecular H₂ formation, that is an important tracer of star formation with the so-called Kennicutt–Schmidt relation (Kennicutt 1998), and which is a key component for gas cooling by rotational transitions or for heating through H₂ formation (Hollenbach & McKee 1979). Also dust directly affects the heating of the gas via photo-electric heating, or the cooling of the gas by depletion of elements in the gas phase, by dust cooling at low temperatures in the dark cores of molecular clouds, or by energy exchange through collisions with electrons in X-ray bright gas. As a contributor to the reprocessing of light in the ISM, dust can have a role in radiative feedback at galactic scales. In particular, the role of dust is certainly exacerbated in massive starburst galaxies which appear as ultra-luminous IR galaxies (or sub-millimetre galaxies at high redshift), as they can obscure up to 99% of their UV-visible emission (Casey et al. 2014).

Observations show that the dust-to-gas ratio (DTG) has a strong correlation with galaxy metallicity (e.g., Lisenfeld & Ferrara 1998; Draine et al. 2007; Galametz et al. 2011; Sandstrom et al. 2013; Rémy-Ruyer et al. 2014), and that the dust-to-metal ratio (DTM) shows a weak but positive correlation with the gas phase metallicity of galaxies (De Cia et al. 2013, 2016; Wiseman et al. 2017). Evolved galaxies, including the Milky Way, have large values of DTM ≃ 0.4, while younger galaxies (with larger gas fraction and lower metallicity) have lower DTM. Large samples of galaxies show that there is also a break in both the DTG- and the DTM–metallicity relation at around 0.1 Z_⊙ (Rémy-Ruyer et al. 2014; De Vis et al. 2019). This break predicted by numerical models (Popping et al. 2017; Hou et al. 2019; Li et al. 2019; Parente et al. 2022) is explained by a transition in the leading mechanism of dust growth in the ISM: from a stellar ejecta-dominated phase at low metallicity to an accretion-dominated phase at high metallicity (Asano et al. 2013a; Zhukovska 2014; Popping et al. 2017; Graziani et al. 2020; Li et al. 2021; Lewis et al. 2023).

Several dust evolution models in galaxies have been developed which include the dust released by stars, the growth of dust by accretion of elements from the gas phase, and the destruction of dust by astration, supernovae and thermal sputtering. In the seminal work of Dwek (1998), they developed a one zone analytical model to follow the evolution of carbonaceous and silicate grains separately, and they were able to reproduce the observed amount of dust in Milky Way-like galaxies. Such one zone models were further refined to track the distribution of grains in sizes including the coagulation and shattering of grains (Asano et al. 2013b). Evolution of dust, followed as a separate component to that of metals from the gas phase, is also modelled in complete hydrodynamical simulations of galaxies (e.g., Bekki 2013; McKinnon et al. 2016; Aoyama et al. 2017; Gjergo et al. 2018; Li et al. 2019, 2021; Graziani et al. 2020; Granato et al. 2021; Trebitsch et al. 2021; Choban et al. 2022). Some of these hydrodynamical simulations make various assumptions about the dust physics, either neglecting the chemical differentiation of grains (as in e.g., Aoyama et al. 2017) or the grain sizes (as in e.g., Graziani et al. 2020; Choban et al. 2022), or both (McKinnon et al. 2016; Trebitsch et al. 2021), or using a unique size distribution for different dust grain species (Li et al. 2021), or using too low resolution to separate ab initio the dense gas structures, in which accretion and coagulation efficiently proceed, from the diffuse ISM (Gjergo et al. 2018; Granato et al. 2021; Li et al. 2021).

Our goal in this paper is to provide a theoretical understanding of the physical origin of the dust properties and, hence, of the different shapes of extinction curves observed as in the Milky Way, or the Magellanic Clouds. Therefore, our approach is to implement and test against various observables a model of dust evolution in the adaptive mesh refinement code RAMSES (Teyssier 2002) that evolves the sizes and the chemical composition of dust. This model assumes a two-size grain approximation, at 5 nm and at 0.1 μm, following separately the evolution of carbonaceous and silicate grains, and is tested in numerical simulations with sufficient resolution to capture the structure of a turbulent multiphase ISM.

The paper is structured as follows. In Sect. 2, we provide an overview of the galactic physics (cooling, star formation, supernova feedback, chemical enrichment) used in RAMSES, and the setup of initial conditions to simulate isolated disc galaxies. In Sect. 3, we describe the details of the model of dust evolution. In Sect. 4, we test our model in a typical Milky Way-like galaxy including stringent tests on extinction curves. In Sect. 5, we extend our predictions to various galaxy masses and explain how Magellanic Cloud-like extinction curves can be obtained. We discuss the limitations of this work in Sect. 6 and, in Sect. 7, we outline our conclusions.

2. Numerics

2.1. Initial conditions and resolution

Our simulations were performed with the adaptive mesh refinement code RAMSES (Teyssier 2002). It solves the hydrodynamics on a cartesian grid with refinement with the Godunov method using the unsplit MUSCL-Hancock scheme with the minmod slope limiter to preserve the total variation diminishing property of the scheme, and the Harten-Lax-van Leer-Contact approximate Riemann solver to obtain the flux at a cell interface. Dark matter and stars are represented by particles that are evolved with a leap-frog scheme. The gravitational potential is obtained by a particle-mesh method where particles are projected on the grid with a cloud-in-cell interpolation.

The spatial resolution was varied from Δx = 72 pc to 9 pc by multiples of 2 with a resolution of 18 pc for our fiducial simulations. The halo mass was either M_h = 10¹⁰, 10¹¹, or 10¹² M_⊙, with a total galaxy mass (gas and stars) that represented f_bar = 3.8% of the halo mass, hence that was respectively M_gal = 3.8 × 10⁸, 3.8 × 10⁹, or 3.8 × 10¹⁰ M_⊙, for simulations that we named respectively G8, G9 and G10. The halo, stellar disc and bulge were initialised with respectively 10⁶, 10⁶ and 10⁵ particles. We allowed refinement where the total mass exceeded 8 × 10³ M_⊙ starting at level 6. In addition to pseudo-Lagrangian mass refinement, we added a refinement criterion on the Jeans length so that it is refined with at least 4 cells down to the maximum level of refinement. The box sizes of respectively 150, 300, and 600 pc for G8, G9, and G10 were evolved with outflowing boundary conditions.

The dark matter (DM), star and gas distributions were initialised with MAKEDISK (Springel et al. 2005), following the strategy introduced in Rosdahl et al. (2017) for RAMSES. The DM mass distribution followed a Navarro et al. (1997) profile with concentration c = 10 and a spin parameter of λ_DM = 0.04 (Macciò et al. 2008). The disc distribution was decomposed into a stellar and a gaseous component (with mass M_⋆,d and M_gas resp.), with density profiles of ρ ∝ e^−R/R_dsech²(z/(2H)), where R is the cylindrical radius, R_d is the disc scale length, z is the disc height, and H is the disc characteristic scale height that we took as H = 0.1R_d. The stellar bulge component (with mass M_⋆,b) followed a Hernquist density profile with bulge scale length equals to H. The gas disc density profile with a temperature of 10⁴ K was used up to a cut radius and a height radius of respectively R_c = 5R_d and H_c = 5H after which, a uniform density of 10⁻⁶ H cm⁻³ and a uniform temperature of 10⁷ K were imposed. For the G8, G9, and G10 galaxies, we used respectively a disc scale radius of 0.7, 1.5, and 3.2 kpc. We varied the gas fraction f_gas =M_gas/(M_gas + M_⋆,d + M_⋆,b) between two extreme values. The initial properties of the simulated galaxies and their halos are shown in Table 1.

The gas metallicity was initialised with a decreasing radial profile, following

$$\begin{array}{cc}\hfill Z(r_{\mathrm{cyl}})& =Z_0{10}^{0.5-R/R_c},\mathrm{if}R<R_{\mathrm{c}}\mathrm{and}z<H_{\mathrm{c}},\hfill \\ \hfill & =0,\mathrm{elsewhere},\hfill \end{array}$$(1)

where R_c = 7.42 kpc, and Z₀ is the scaling metallicity of the profile, which took values of Z₀ = 0.045, 0.15, 0.45, 0.90 Z_⊙, with Z_⊙ = 0.01345, which gave respectively a mean metallicity of the gas of 0.1, 0.3, 1, 3 Z_⊙. The various elements followed in this work (H, He, C, N, O, Mg, Fe, Si, and S) were initialised in the gas phase with solar mass fractions (Asplund et al. 2009). The exact mass fraction of each individual element is a function of the metallicity and of star formation history of the galaxy, as can be observed for low metallicity stars in the Milky Way, but we neglected this aspect. We assumed that 10% of the initial refractory¹ elements were trapped in dust grains for the G10LG galaxy and 0.1 per cent for the lower mass galaxies G8 and G9, whose mass were equally distributed in sizes, which corresponded to a DTM = 3.6 × 10⁻² and 3.6 × 10⁻⁴ respectively.

For galaxies simulated with the large gas fractions (i.e., galaxy G8HG, G9HG, and G10HG), a special procedure to smoothly relax the initial conditions was required. Due to the very idealised nature of the disc initial conditions, during the first 10s of Myr, the gas loses pressure support due to rapid cooling and produces a strong burst of star formation rate (SFR). In turn this burst of SFR creates a massive outflow that can remove very significant amounts of gas from the galaxy and turn a gas-rich galaxy into a gas-poor galaxy. In order to circumvent this shortcoming, and after some tests, we introduced a perturbation in the initial vertical velocity field with an amplitude of 30 km s⁻¹ modulated by two sinusoids in the x and y axis of the cartesian box (the z axis is oriented along the spin axis of the disc) with wavelength of 200 pc. In addition to this initial perturbation of the velocity field, we also imposed that the star formation efficiency (see Sect. 2.3) linearly ramped up with time with a timescale of 100 Myr.

2.2. Gas radiative cooling and heating

Gas was allowed to cool down to 10 K through H and He collisions with a contribution from metals in the gas phase using rates tabulated by Sutherland & Dopita (1993) above 10⁴ K and those from Dalgarno & McCray (1972) below 10⁴ K. Dust did not directly contribute to the gas cooling rates at low temperature, however, these rates can be significantly reduced in regions where the dust fraction over total metallicity is large, by decreasing the amount of metals in the gas phase available for gas cooling. The contribution of metals to the overall cooling rates was reduced by the corresponding amount of metals locked into the dust phase. We note that this approach is not consistent since the cooling curves, as tabulated in this version of RAMSES, are for a total metallicity with a solar composition of elements, while the amount of depletion is not uniform across all elements, and, hence, this should affect differently the cooling curves as a function of temperature. A more consistent approach will be explored using the PRISM ISM model for RAMSES from Katz et al. (2022) that computes cooling rates from the amount of each individual element and their ionisation state in the gas phase (Rodríguez Montero et al., in prep.).

At high temperatures, fast moving electrons frequently collide with grains, leading to an exchange of internal energy from the gas phase to the dust phase. Since dust radiates the stored energy into the IR, this leads to a net cooling for the gas phase. We used the corresponding cooling rates from Dwek & Werner (1981):

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{D}}(a,T)=\mathcal{D}\frac{3m_{\mathrm{H}}}{X_{\mathrm{H}}4\pi s{a}^3}{n}_{\mathrm{H}}{n}_{\mathrm{e}}\mathcal{H}(a,T),\end{array}$$(2)

where a is the grain radius, s is the grain material density, 𝒟 is the dust-to-gas ratio, n_H is the hydrogen number density, n_e is the (free) electron number density (for which we assume a gas fully ionised at these temperatures), X_H = 0.76 is the hydrogen mass fraction, m_H is the hydrogen mass, and

$$\begin{array}{c}\hfill \frac{\mathcal{H}(a,T)}{{\mathrm{erg}\mathrm{s}}^{-1}c{m}^3}=\{\begin{array}{c}0,x\ge x_{\max },\hfill \\ 5.38\times {10}^{-18}{\left(\frac{a}{\mathrm{\mu }\mathrm{m}}\right)}^2{\left(\frac{T}{\mathrm{K}}\right)}^{1.5},4.5\le x<x_{\max },\hfill \\ 3.37\times {10}^{-13}{\left(\frac{a}{\mathrm{\mu }\mathrm{m}}\right)}^{2.41}{\left(\frac{T}{\mathrm{K}}\right)}^{0.88},1.5\le x<4.5,\hfill \\ 6.48\times {10}^{-6}{\left(\frac{a}{\mathrm{\mu }\mathrm{m}}\right)}^3,x<1.5,\hfill \end{array}\end{array}$$(3)

where x = 2.71 × 10⁸(a/μm)^2/3(T/1 K)⁻¹, and x_max = 14 000(a/μm)^2/3, which was introduced to impose a sharp cut-off to the dust cooling function at temperatures below 2 × 10⁴ K in order to prevent the contribution to this cooling from dust in the cold phase. While the dust contribution to the gas cooling rate can largely exceed that of the fiducial gas cooling at T > 10⁶ K, thermal sputtering from ion collisions efficiently erodes dust grains (see Sect. 3.3), such that grains are destroyed on a timescale shorter than that of the gas cooling from dust (Guillard et al. 2009), which reduces the effective contribution to cooling at high temperatures (Montier & Giard 2004; Pointecouteau et al. 2009; Vogelsberger et al. 2019).

2.3. Star formation

Star formation followed a Schmidt law: ${\dot{\rho }}_{\star }={\epsilon }_{\star }{\rho }_{\mathrm{g}}/t_{\mathrm{ff}}$, where ${\dot{\rho }}_{\star }$ is the star formation rate mass density, ρ_g the gas mass density, t_ff the local free-fall time of the gas, and ε_⋆ is the star formation efficiency. Star formation occurred in regions with gas number density above n₀ which varied with resolution (n₀ = 0.6, 2.5, 10, 40 H cm⁻³ for a minimum spatial resolution of respectively Δx = 72, 36, 18, 9 pc). This corresponds to a stellar mass resolution for newly formed stars of respectively = 7.4 × 10³, 3.7 × 10³, 1.9 × 10³, 0.9 × 10³ M_⊙.

As in Dubois et al. (2021) (see for details of the implementation and the underlying parameters), the star formation efficiency varied with the virial parameter α_vir = 2E_turb/|E_g|< 1, with E_turb and E_g the turbulent kinetic energy and the gravitational energy of the gas respectively, and with the turbulent Mach number ℳ = σ_t/c_s, where σ_t is the gas turbulent velocity and c_s is the sound speed. In practice, σ_t was computed by taking for each cell its own gas velocity and that of its neighbouring cells into account, removing the local bulk (mass-weighted) mean velocity, and constructing σ through σ² = sum(∇ ⊗ udx)². In this theory (e.g., Krumholz & McKee 2005; Padoan & Nordlund 2011; Hennebelle & Chabrier 2011; Federrath & Klessen 2012), star formation (SF) is driven by how much gas mass passes a given density threshold. This amount is controlled by the level of turbulence and how much the molecular cloud is gravitationally bound. The unresolved density distribution can be described by a log-normal probability density function (PDF), with standard deviation given by ℳ and α_vir: the more turbulent and bound the gas is, the larger the SF efficiency, reaching values as high as ε_⋆ ∼ 1. The mean SFR between 200 Myr and 400 Myr for the different simulated galaxies are 0.01 M_⊙ yr⁻¹, 0.03 M_⊙ yr⁻¹, 0.5 M_⊙ yr⁻¹, and 0.7 M_⊙ yr⁻¹ for respectively G8HG, G9LG, G9HG, and G10LG.

2.4. Stellar yields and feedback

We limited our simulations to individually track the amount of C, N, O, Mg, Si, Fe, S (which represented ∼85% of the mass of solar heavy elements), of total metallicity Z, and H, but these can easily be extended to other individual elements. We assumed a Chabrier (2005) initial mass function for the distribution of zero age star (ZAS) masses. Stellar gas plus dust yields for intermediate stars were those of Karakas (2010) for ZAS masses M_ZAS < 8 M_⊙ (asymptotic giant branch, AGB, stars), while yields from Limongi & Chieffi (2018) for massive stars M_ZAS ≥ 8 M_⊙ were employed. The yields for massive stars were parametrised by the rotation velocity of the star (either 0, 150 or 300 km s⁻¹). We followed the initial distribution of rotational velocities as a function of the initial metallicity of the star from Prantzos et al. (2018; obtained from multiple Milky Way chemical constraints), with the corresponding grid² of velocity-to-metallicity of ZAS: V_ZAS = 150, 100, 50, 50, 50, 50, 50 km s⁻¹ for Z_ZAS = 10⁻³, 10⁻², 10⁻¹, 10^−0.6, 10^−0.3, 1, 10^0.3 Z_⊙ (where Z_⊙ = 0.01345 is the solar value of metallicity from Asplund et al. 2009). Since the Limongi & Chieffi (2018) stellar yields largely underproduce the amount of Mg, with respect to Si and solar abundances (Prantzos et al. 2018), we artificially boosted up by a factor of 2 the amount of Mg released by massive stars. The mass distribution of stars was assumed to cover the [0.1, 100] M_⊙ mass range, assuming that massive stars only successfully explode in SNII and release their corresponding SN-processed yields for M_ZAS ≤ 30 M_⊙, while they also release mass through winds over their entire mass range ([8, 100] M_⊙). The death of each individual intermediate mass star followed the Padova stellar tracks (Girardi et al. 2000) with thermally pulsating AGB (Vassiliadis & Wood 1993). We further assumed that the mass release of individual AGB stars happens instantaneously at the end of their evolution (pulsating evolution is not individually time-resolved).

For type Ia SN (SNIa), we assumed the delay time distribution of Maoz & Mannucci (2012) with a time-integrated number of SNIa of $N_{\mathrm{Ia}}(<t)=2.35\times {10}^{-4}(log(t/\mathrm{Gyr})+3)M_{\odot }^{-1}$ (limiting the range of SNIa between 50 Myr and 13.7 Gyr), together with the stellar yields of Iwamoto et al. (1999; using their W70 carbon-deflagration model). The resulting mass release for each individual elements are shown in Appendix A.

We simplified further the feedback from SNII by assuming that all mass and energy are released in one single explosive event after 5 Myr to maximise the impact of SNII feedback. We assumed a Chabrier (2005) initial mass function, with a SNII explodability range of M_ZAS = [8, 30] M_⊙ of canonical individual kinetic energy of 10⁵¹ erg, hence, corresponding to a SNII specific energy of ${10}^{49}\text{erg}{M}_{\odot }^{-1}$. We included the SN feedback model from Kimm & Cen (2014), where the model follows either the energy conserving or momentum conserving phase of the explosion. If the SN explosion is still in the energy-conserving phase, only internal energy is given to the gas since the code will be able to handle the Sedov explosion and get the correct amount of momentum at the end of the Sedov phase. However, out of the energy-conserving phase, the correct amount of momentum is given to the gas. We also included a minimal model for UV radiation from OB young stars, by considering the larger amount of momentum SNII can impart to the gas thanks to the pre-heating by the UV-radiation (Geen et al. 2015).

3. Dust

The two bin size decomposition of Hirashita (2015; see also the implementations of Aoyama et al. 2017; Gjergo et al. 2018, or Granato et al. 2021) was adopted for our subgrid model of dust evolution. These bins correspond to a small bin and a large bin of dust sizes (radii; a) of respectively a_S = 5 nm and a_L = 0.1 μm. This decomposition allows to capture the bulk of the grain size distribution (Hirashita 2015), while having a negligible impact on the computational footprint (and on the memory load) of hydrodynamical simulations of galaxy formation.

We separated the composition of dust between carbonaceous grains and silicate grains. For silicates, we assumed a fixed amorphous olivine composition Mg_2 − xFe_xSiO₄ with iron inclusions balancing the amount of Mg (x = 1), i.e., MgFeSiO₄, as it is supposed to make the bulk of silicates (Kemper et al. 2004; Min et al. 2007), although there are still large modelling uncertainties plaguing the reconstructed compositions of silicates, particularly on the depletion of oxygen (Jenkins 2009). Assuming pyroxene with MgFeSi₂O₆ instead of olivine typically increased the amount of dust mass released as silicates by 20% for our stellar production models. The respective grain material density of carbonaceous and silicate grains were respectively s_C = 2.2 g cm⁻³ (i.e. we assumed a solid structure close to graphite) and s_Sil = 3.3 g cm⁻³.

The equations of evolution of the dust mass content D_i, j of each grain size bin (with subscript j = S and L for respectively small and large grains) with chemical composition from the ‘key’ element i, which is C for carbonaceous grains, and Si for silicate grains, are:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{D}_{i,\mathrm{S}}}{\mathrm{d}t}& =\frac{\mathrm{d}{D}_{\mathrm{acc},i,\mathrm{S}}}{\mathrm{d}t}+\frac{\mathrm{d}{D}_{\mathrm{ej},i,\mathrm{S}}}{\mathrm{d}t}+\frac{\mathrm{d}{D}_{\mathrm{sha},i,\mathrm{L}}}{\mathrm{d}t}\hfill \\ \hfill & -\frac{\mathrm{d}{D}_{\mathrm{coa},i,\mathrm{S}}}{\mathrm{d}t}-\frac{\mathrm{d}{D}_{\mathrm{spu},i,\mathrm{S}}}{\mathrm{d}t}-\frac{\mathrm{d}{D}_{\mathrm{SN},i,\mathrm{S}}}{\mathrm{d}t}-\frac{\mathrm{d}{D}_{\star ,i,\mathrm{S}}}{\mathrm{d}t},\hfill \\ \hfill \frac{\mathrm{d}{D}_{i,\mathrm{L}}}{\mathrm{d}t}& =\frac{\mathrm{d}{D}_{\mathrm{acc},i,\mathrm{L}}}{\mathrm{d}t}+\frac{\mathrm{d}{D}_{\mathrm{ej},i,\mathrm{L}}}{\mathrm{d}t}+\frac{\mathrm{d}{D}_{\mathrm{coa},i,\mathrm{S}}}{\mathrm{d}t}\hfill \\ \hfill & -\frac{\mathrm{d}{D}_{\mathrm{sha},i,\mathrm{L}}}{\mathrm{d}t}-\frac{\mathrm{d}{D}_{\mathrm{spu},i,\mathrm{L}}}{\mathrm{d}t}-\frac{\mathrm{d}{D}_{\mathrm{SN},i,\mathrm{L}}}{\mathrm{d}t}-\frac{\mathrm{d}{D}_{\star ,i,\mathrm{L}}}{\mathrm{d}t},\hfill \end{array}$$(4)

where for each equation the first three terms on the right-hand side stand for the increase in dust mass and the last four terms for the decrease in dust mass. The various terms stand for: the accretion of elements from the gas phase to the dust phase (see Sect. 3.4); the ejecta mass release from SSP evolution (see Sect. 3.1); the shattering and coagulation of dust grains which transfer dust from, respectively, large to small sizes (Sect. 3.5), and small to large sizes (Sect. 3.6); the thermal sputtering that returns dust elements to the gas phase (see Sect. 3.3); and the destruction of dust by shocks from type II SNe (see Sect. 3.2).

In this work we neglected entirely the differential dynamical motion of dust with respect to gas dynamics, therefore, the dust content of each bin size and chemical composition can be treated as a passive variable that is transported along with the flow of gas in the exact same vein as for the transport of metals in the gas phase. This is possible due to the relatively short stopping time associated to Epstein aero-dynamical drag in dense gas t_st ≃ 0.5 a_0.1(n/1 cm⁻³)⁻¹(c_s/10 km s⁻¹)⁻¹ Myr compared to dynamical time t_dyn ≃ 10 (L/100 pc)(σ_gr/10 km s⁻¹)⁻¹ Myr, where σ_gr is the grain velocity dispersion. This gives a diffusion length scale of L ≃ 5 a_0.1(n/1 cm⁻³)⁻¹(σ_gr/c_s) pc, which is in general below the resolution scale of this work. We defer the investigation of the diffusion of dust grains on cloud scales to future work.

3.1. Stellar production

The amount of silicates m_key, traced by its key element (here Si), released by a single star is given by the least abundant element entering its composition

$$\begin{array}{c}\hfill m_{\mathrm{key}}=\underset{\ell =\mathrm{Mg},\mathrm{Fe},\mathrm{Si},\mathrm{O}}{min}\left(\frac{M_{\ell }}{A_{\ell }{N}_{\ell }}\right),\end{array}$$(5)

with M_ℓ, A_ℓ, N_ℓ, respectively, the amount of mass of the ℓth element released by the star, the atomic weight of the ℓth element, and the number of atoms of the ℓth element entering the composition of the silicate. This available number of atoms of the key element condensates into dust within the ejecta at a given efficiency δ_Sil:

$$\begin{array}{c}\hfill D_{\mathrm{ej},\mathrm{sil},\ell }^k={\delta }_{\mathrm{Sil}}^k{m}_{\mathrm{key}}{A}_{\ell }{N}_{\ell },\end{array}$$(6)

which guarantees a constant mass ratio (A_ℓN_ℓ/∑(A_ℓN_ℓ) = 0.14, 0.33, 0.16, 0.37 for resp. ℓ = Mg, Fe, Si, and O in x = 1 olivine³) between the various elements entering the composition of silicate grains. The superscript k stands for the different type of mass release by stars with k = AGB, SNII or SNIa, The corresponding efficiencies were equal to ${\delta }_{\mathrm{Sil}}^{\mathrm{SNII}}=0.8$ for SNII following Dwek (1998; and ${\delta }_{\mathrm{Sil}}^{\mathrm{SNIa}}=1$ for SNIa⁴). For AGB stars, the condensation efficiency depended on the ratio of C/O: if C/O ≥ 1 all the oxygen was associated to carbon in CO molecules, and, hence, silicates did not condensate anymore as a result of unavailable O, while with C/O < 1 oxygen was still available to condensate into silicates. We adopted ${\delta }_{\mathrm{Sil}}^{\mathrm{AGB},\mathrm{C}/\mathrm{O}<1}=0.8$ and ${\delta }_{\mathrm{Sil}}^{\mathrm{AGB},\mathrm{C}/\mathrm{O}\ge 1}=0$ as in Dwek (1998). For the same reasons, accretion of elements onto silicates from the ISM was limited by the least abundant of the elements in the ISM entering the composition of the grain, replacing the mass M_ℓ by the mass of gas of the given element ℓ in the cell.

Similarly to silicates, the amount of carbonaceous dust and the condensation efficiency δ_C varied with the nature of the stellar ejecta. It was assumed to be ${\delta }_{\mathrm{C}}^{\mathrm{SNII}}=0.5$ and ${\delta }_{\mathrm{C}}^{\mathrm{SNIa}}=0.5$, while for AGB stars this efficiency depended on the ratio of C/O in the ejecta. If the ratio of C/O < 1, then all the C elements formed into CO molecules that were not available for dust ${\delta }_{\mathrm{C}}^{\mathrm{AGB},\mathrm{C}/\mathrm{O}<1}=0$, while for C/O ≥ 1 the efficiency was 1 with a mass released into C dust:

$$\begin{array}{c}\hfill D_{\mathrm{ej},\mathrm{C}}^{\mathrm{AGB}}={\delta }_{\mathrm{C}}^{\mathrm{AGB},\mathrm{C}/\mathrm{O}\ge 1}(M_{\mathrm{C}}^{\mathrm{AGB}}-\frac{A_{\mathrm{C}}}{A_{\mathrm{O}}}{M}_{\mathrm{O}}^{\mathrm{AGB}}).\end{array}$$(7)

The production of dust in stellar ejecta is balanced between a population of small and large grains, i.e., ${\dot{D}}_{\mathrm{ej},i,\mathrm{S}}=f_{\mathrm{ej},i,\mathrm{S}}{\dot{D}}_{\mathrm{ej},i}$ and ${\dot{D}}_{\mathrm{ej},i,\mathrm{L}}^k=(1-f_{\mathrm{ej},i,\mathrm{S}}){\dot{D}}_{\mathrm{ej},i}$, where f_ej, i, S is the fraction of small grains released for a given chemical type i (carbonaceous or silicate grains). For sake of simplicity, we assumed that all the dust released in the ejecta is condensed into the large grain size population only, thus, f_ej, i, S = 0 (see Appendix B for the predictions of extinction curves with different values of f_ej, i, S, and the discussion in Sect. 6). Finally, it is important to note that the predicted and observed dust condensation efficiencies are very uncertain and can vary by an order of magnitude (see Schneider & Maiolino 2023).

3.2. Supernova-driven destruction

SNe (type II and Ia) destroy dust already present in the ISM, where the dust destruction is produced by inertial (non-thermal) sputtering and grain collisions (Kirchschlager et al. 2022). The amount of SN-destroyed dust is a fraction of the mass of gas M₁₀₀ shocked at above 100 km s⁻¹ defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }{M}_{\mathrm{SN},i,j}={\epsilon }_{\mathrm{SN},i}(a_j)\min (\frac{M_{100}}{M_{\mathrm{g}}},1)M_{\mathrm{D},i,j},\end{array}$$(8)

with the size-dependent destruction efficiency of Aoyama et al. (2020)

$$\begin{array}{c}\hfill {\epsilon }_{\mathrm{SN},i}(a_j)=1-exp(-\frac{{\delta }_{\mathrm{SN},i}}{a_{0.1,j}}),\end{array}$$(9)

where M_g is the local gas mass, M_d, i, j is the local dust mass, a_0.1, j is the dust grain size in units of 0.1 μm, and δ_SN, i is the destruction efficiency. Small grains are more efficiently decelerated by drag forces, trapping them near the shock region where thermal sputtering can quickly destroy them (Nozawa et al. 2006). This functional form of the destruction efficiency has been constructed to capture this behaviour with grain size a. The mass of shocked gas is provided by the Sedov solution in a medium of homogeneous gas density, i.e. M₁₀₀ = 6800E_SN, 51 M_⊙, where E_SN, 51 is the energy of the SN explosion in units of 10⁵¹ erg. Multiple SNe were released at once over one time step Δt (each individual SN engulfing the gas swept up by the previous explosion), following (see Hou et al. 2017):

$$\begin{array}{c}\hfill \mathrm{\Delta }{M}_{\mathrm{SNe},i,j}=[1-{(1-\frac{\mathrm{\Delta }{M}_{\mathrm{SN},i,j}}{M_{\mathrm{D},i,j}})}^{N_{\mathrm{SN}}}]M_{\mathrm{D},i,j},\end{array}$$(10)

with N_SN the number of SNe. We stress that this destruction step by the SN blast solution does explicitly destroy the returned ejecta from the stellar production, and destroys only the background dust material, hence, the dust condensation efficiency implicitly takes this term into account. We used a different value of the dust destruction efficiency for carbonaceous and silicate grains with respectively δ_SN, C = 0.10 and δ_SN, Sil = 0.15 following the qualitative behaviour of Hu et al. (2019), where silicate grains are approximately destroyed 50 per cent more than carbonaceous grains.

3.3. Thermal sputtering

Dust is also destroyed by thermal sputtering, the ejection of atoms from grains by the transfer of kinetic energy from gas particles at temperatures high enough to overcome the energy barrier of the binding energy. We adopted the fitting form of Hu et al. (2019) to the thermal sputtering yields Y_th calculated by Nozawa et al. (2006), which relates to the thermal sputtering timescale as (recall that t_spu,i,j = m_gr,i,j/|ṁ_gr,i,j| = a_i,j/(3|ȧ_i,j|) with m_gr, i, j the mass of a single dust grain):

$$\begin{array}{c}\hfill t_{\mathrm{spu},i,j}=\frac{a_i}{3n_{\mathrm{H}}{Y}_{\mathrm{th},j}(T)}.\end{array}$$(11)

We note that sputtering yields are approximately 3 times lower at T ≳ 10⁶ K for carbonaceous grains than for silicate grains (see also Draine & Salpeter 1979; Tielens et al. 1994). Hence, they differ from the widely adopted yields of Tsai & Mathews (1995) giving a unique sputtering timescale based on the intermediate values for silicates and carbonaceous grains of Draine & Salpeter (1979) and Tielens et al. (1994; see Appendix C). Using cosmological hydrodynamical simulations of galaxy clusters, Gjergo et al. (2018) and Vogelsberger et al. (2019) show that the dust mass content in observed galaxy clusters can be reproduced only for ten times lower sputtering yields than the canonical values. We defer this investigation to future work, but we immediately stress that thermal sputtering is already a sub-dominant destruction mechanism in our set of simulations compared to direct destruction by SNe. We performed a G10LG simulation with sputtering time artificially increased by a factor of 10: there is no noticeable difference in any of the dust properties investigated in Sect. 4.

3.4. Gas accretion

The dust mass content can finally grow by accretion through the metals in the gas phase, following Dwek (1998)

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc},i,j}=(1-\frac{M_{\mathrm{D},i,j}}{M_{Z,i}})\frac{M_{\mathrm{D},i,j}}{t_{\mathrm{acc},i,j}},\end{array}$$(12)

where M_Z, i is the mass of metals (gas + dust) of the key element, t_{acc, i, j} is the accretion timescale

$$\begin{array}{c}\hfill t_{\mathrm{acc},i,j}^{-1}=\alpha \frac{\pi {\overline{a}}_j^2}{f_X{m}_{\mathrm{gr},i,j}}{\rho }_X{u}_{\mathrm{th},X},\end{array}$$(13)

where $m_{\text{gr},i,j}=s_i{a}_j^{3}4\pi /3$ is the grain mass, s_i is the grain material density, f_X is the mass fraction of the limiting element in the chemical composition of the grain, u_th, X = (8k_BT/(πm_X))^1/2 is the gas thermal velocity, m_X = A_Xm_H is the atomic mass, and ρ_X is the total (gas plus dust phase) mass density of the limiting element. $\pi {\overline{a}}_{i,j}^2$ is the surface of the grain with Coulomb enhancement factor E_i(a_j) due to electrostatic effects caused by ionised gas interacting with charged grains (Weingartner & Draine 1999), which is integrated over the whole distribution of the grain sizes. In our case, where the size distribution was discretised over two bins and assuming a Dirac distribution, we simply got ${\overline{a}}_{i,j}=E_i(a_j)a_j$. For carbonaceous grains, X = i = C, however, for silicate grains, the limiting element depends on the grain composition, i.e.,

$$\begin{array}{c}\hfill \frac{{\rho }_X{u}_{\mathrm{th},X}}{f_X}=\min {\left[\frac{{\rho }_{\ell }{u}_{\mathrm{th},\ell }}{f_{\ell }}\right]}_{\ell =\mathrm{Mg},\mathrm{Fe},\mathrm{Si},\mathrm{O}}.\end{array}$$(14)

The accretion timescale can be rewritten as

$$\begin{array}{c}\hfill t_{\mathrm{acc},i,j}=0.28{\alpha }^{-1}{a}_{0.005,j}{s}_{3,i}\frac{A_X^{\frac{1}{2}}{f}_X}{E_i(a_j)Z_X}{\left(\frac{n}{1\mathrm{H}{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{T}{50\mathrm{K}}\right)}^{-\frac{1}{2}}\mathrm{Myr},\end{array}$$(15)

where a_0.005, j = a_j/5 nm, s_3, i = s_i/3 g cm⁻³, α is the sticking coefficient of gas particles onto dust, and Z_X = ρ_X/ρ_g is the mass fraction of the limiting element. Weingartner & Draine (1999) gave values of E(a) for a range of sizes of carbonaceous and silicate grains for three typical phases of the ISM. To mimic their complex dependencies, we simply adopted a value of E = 1 for all kind of grains everywhere, except for large carbonaceous and small silicate grains in the cold neutral medium with T < 2 × 10⁴ K and n > 10 H cm⁻³, for which we adopted E_C(0.1 μm) = 0 and E_Sil(5 nm) = 10.

Sub-grid accretion can be estimated assuming that the gas density with mean value $\overline{n}$ (the value of the gas density in a given resolution element) follows a log-normal PDF at unresolved scales, which is shaped by the turbulence:

$$\begin{array}{c}\hfill p(s)=\frac{1}{\sqrt{2\pi }{\sigma }_s}exp(-\frac{{(s-s_0)}^2}{2{\sigma }_s^2}),\end{array}$$(16)

where $s=log(n/\overline{n})$, $s_0=-{\sigma }_s^2/2$ and σ_s = log(1 + b²ℳ²), with b the compression ratio (we take b = 0.4 for mixed turbulence which is the same value used in the gravo-turbulent model for SF efficiency) and ℳ the turbulent Mach number. The effective timescale t_acc, eff was obtained by integrating the mass accretion rate over the log-normal PDF for a given set of mean density, temperature, and metallicity ($\overline{n},\overline{T},\overline{Z}$), hence:

$$\begin{array}{c}\hfill \frac{t_{\mathrm{acc}}(\overline{n},\overline{T},\overline{Z})}{t_{\mathrm{acc},\mathrm{eff}}}={\displaystyle {\int }_{-\infty }^{s_{\max }}}\frac{t_{\mathrm{acc}}(\overline{n},\overline{T},\overline{Z})}{t_{\mathrm{acc}}(n,T,Z)}\frac{n}{\overline{n}}p(s)\mathrm{d}s,\end{array}$$(17)

up to a maximum density $s_{\max }=log(n_{\max }/\overline{n})$, where the maximum density is n_max = 10⁴ H cm⁻³ where grains starts to be significantly coated with water ice mantle (n ≥ 10³ H cm⁻³, Cuppen & Herbst 2007; Hollenbach et al. 2009) suppressing the accretion of refractory material onto the grain surface. We used the value of the gas density and metallicity of a given cell for $\overline{n}$ and $\overline{Z}$, and we assumed that the gas temperature is always $\overline{T}=100\mathrm{K}$ for gas that fulfils the conditions to trigger this unresolved log-normal PDF of density (see below) for simplicity. In addition, we assumed that $T=\overline{T}$ and $Z=\overline{Z}$ for the sampled values of n by the log-normal PDF. This led to

$$\begin{array}{cc}\hfill \frac{t_{\mathrm{acc}}(\overline{n},\overline{T},\overline{Z})}{t_{\mathrm{acc},\mathrm{eff}}}& ={\displaystyle {\int }_{-\infty }^{s_{\max }}}exp(2s)p(s)\mathrm{d}s\hfill \\ \hfill & =\frac{e^{{\sigma }_s^2}}{2}\mathrm{Erfc}\left(\frac{3{\sigma }_s^2/2-s_{\max }}{\sqrt{2}{\sigma }_s}\right),\hfill \end{array}$$(18)

where Erfc is the complementary error function. This unresolved effective accretion timescale was estimated only for cells with gas density larger than 0.1 H cm⁻³, temperature lower than 10⁴ K, and Jeans length smaller than 4 times the local cell size, assuming a sticking coefficient of α_eff = 1/3, which value is representative of the sticking coefficient obtained at temperatures of T = 100 − 1000 K (see Fig. D.1). Otherwise the accretion time was that obtained for $\overline{n}$ (the local gas density) with the Le Bourlot et al. (2012) sticking coefficient (see Appendix D for a comparison of sticking coefficients).

For carbonaceous grains, we further assumed that carbon in the gas phase is fully molecular (CO) above a gas density of n > 10³ H cm⁻³. Indeed, simulations (e.g., Safranek-Shrader et al. 2017; Clark et al. 2019; Hu et al. 2021; Katz et al. 2022) give a transition to full CO gas at n = 2 × 10³ − 10⁴ H cm⁻³ depending on the cosmic ray ionisation rate. We note that while our value might seem fairly low (or corresponding to a very weak cosmic ray ionisation rate), our subgrid gas distribution neglects the fact that gas collapses further at smaller unresolved scales (higher density) and extends the log-normal PDF with a power-law tail (Vallini et al. 2018).

We will see in Sect. 4.2 that this turbulence-driven subgrid accretion model allows for better converged amounts of dust with respect to resolution (as opposed to ignoring the unresolved density structure of the gas). Other work in the literature had to make similar, though ad hoc, assumptions on how much gas concentrates into typical molecular cloud densities (10³ H cm⁻³, Aoyama et al. 2017; Granato et al. 2021), or using a Sobolev-like shielding length to infer the amount of unresolved dense gas (Choban et al. 2022).

3.5. Shattering

Shattering is the process by which grains with sufficiently large velocities collide and fragment into smaller size grains. Therefore, it is a transfer of mass from the large to the small grain population that conserves the total dust mass. The corresponding timescale is obtained by considering the timescale for grain collisions ${\tau }_{\mathrm{col}}={({\sigma }_{\mathrm{D}}\stackrel{\sim }{\sigma }{n}_{\mathrm{D}})}^{-1}$, where $\stackrel{\sim }{\sigma }$ is the grain cross section, σ_D the dust grain velocity dispersion, and n_D the dust number density, which boils down to (Aoyama et al. 2017):

$$\begin{array}{c}\hfill t_{\mathrm{sha},\mathrm{i}}=54a_{0.1,\mathrm{L}}{s}_{3,i}{\left(\frac{{\mathcal{D}}_{\mathrm{L}}}{0.01}\right)}^{-1}{\left(\frac{n}{1{\mathrm{H}\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{{\sigma }_{\mathrm{D},\mathrm{L}}}{10{\mathrm{km}\mathrm{s}}^{-1}}\right)}^{-1}\mathrm{Myr},\end{array}$$(19)

where 𝒟_L and σ_D, L are, respectively, the dust-to-gas ratio and velocity dispersion of the large grains population. Typical velocities required to shatter grains are above a few km s⁻¹ (Jones et al. 1996), which limit the effect of shattering to the diffuse phase of the ISM. Yan et al. (2004) provided grain-gas relative velocities for various phases of the ISM and showed that large grains (a_L = 0.1 μm) reach the turbulent velocity of the forcing scale, i.e. σ_D,L ≃ 10 km s⁻¹ in the warm ionised medium, ≃1 km s⁻¹ for the warm or cold neutral medium, and ≃0.1 − 1 km s⁻¹ for molecular clouds. We followed Granato et al. (2021) and adopted the following functional form:

$$\begin{array}{c}\hfill t_{\mathrm{sha}}=54a_{0.1}{s}_3{\left(\frac{{\mathcal{D}}_{\mathrm{L}}}{0.01}\right)}^{-1}{\left(\frac{n}{1\mathrm{H}{\mathrm{cm}}^{-3}}\right)}^{-p_{\mathrm{sh}}}\mathrm{Myr},\end{array}$$(20)

with p_sh = 1 for n < 1 H cm⁻³ and p_sh = 1/3 for 1 < n < 10³ H cm⁻³, which is equivalent to having a dispersion velocity varying from 10 km s⁻¹ to 1 km s⁻¹ from the warm ionised phase (1 H cm⁻³) to the molecular cloud density⁵ (10³ H cm⁻³).

3.6. Coagulation

When dust grains are embedded in gas that is dense and cold, they have small velocity dispersions. This leads to the coagulation of grains by their direct collisions, thereby, transferring mass from small grains to large grains (Yan et al. 2004). We set the timescale for grain coagulation to that obtained by considering the typical timescales obtained previously for shattering but now for the coagulation of small grains into large grains, i.e.,

$$\begin{array}{c}\hfill t_{\mathrm{coa}}=0.27\frac{a_{0.005}{s}_3}{F}{\left(\frac{{\mathcal{D}}_{\mathrm{S}}}{0.01}\right)}^{-1}{\left(\frac{n}{{10}^3\mathrm{H}{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{{\sigma }_{\mathrm{D},\mathrm{S}}}{0.1\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-1}\mathrm{Myr},\end{array}$$(21)

where F is a fudge factor. The density to achieve grain velocity dispersions that are low enough to pass below the coagulation threshold velocity of ≲0.1 − 1 km s⁻¹ for grains of size a = 0.005 μm (Chokshi et al. 1993; Poppe & Blum 1997) is of the order n ∼ 10² − 10³ H cm⁻³, which is barely resolved at resolutions of Δx ≃ 20 pc as in this work.

Since the velocity dispersion of small grains has complex dependencies with the gas turbulent velocity dispersion at the forcing scale, the magnetic field strength, the grain charge, and the ionisation state of the gas (Yan et al. 2004), we greatly simplified the problem (following Aoyama et al. 2017) by assuming that half (F = 0.5) the gas mass, which Jeans length is unresolved and with gas density above 0.1 H cm⁻³ and temperature below 10⁴ K, has an actually larger gas density of 10³ H cm⁻³ with a small grain turbulent velocity dispersion of 0.1 km s⁻¹ (Yan et al. 2004). We also performed two simulations with a coagulation model where F = 1 and where n is sampled by the log-normal PDF shaped by turbulence as for the subgrid model for dust growth by accretion. In this model, we tried two different density cuts, n_max, coa = 10⁵ or 10⁷ H cm⁻³, and their results show very negligible differences with respect to the fiducial coagulation model for G10LG (they produce respectively ∼0.2 and 3% less small grains than the fiducial G10LG simulation at final time).

4. Results for the Milky Way-like galaxy

The entire set of simulations with variations in resolution, metallicity, or dust physics are summarised in Table 2. We will first focus on the fiducial model of our Milky Way-like galaxy G10LG to validate the model of dust physics against available observations for the Milky Way.

4.1. Visual inspection

In Fig. 1, we show a mock image of the fiducial G10LG simulation at t = 400 Myr obtained with the radiative transfer code SKIRT (Camps & Baes 2015) in the g − r − i filter bands. For the purpose of this image, the stars set up from initial conditions have an age with a cylindrical radial dependency of age = (9 − R_cyl/2 kpc) Gyr ± 3 Gyr and have solar metallicities. One can see that the obtained light distribution shows a diffuse old stellar component structured spiral arms, together with more clumpy blue regions of star formation, with absorption features from dusty gas.

Fig. 1. Mock image of the fiducial G10LG simulation seen edge-on (top) or face-on (bottom) in SDSS g − r − i filter bands at t = 400 Myr.

Figure 2 gives the visual representation of various gas quantities in the galaxy of the fiducial run, as mass-weighted edge-on (30 kpc in depth) and face-on (0.6 kpc in depth) projections at time t = 400 Myr. The gas is clearly multiphase and strongly clustered into regions of neutral gas with high densities (n > 10 − 100 H cm⁻³) and low temperatures (T ≃ 100 K) embedded in a more diffuse ionised medium at intermediate density (n ≃ 0.1 − 1 H cm⁻³) and at warm temperatures (T ≃ 10⁴ K), and SN-driven hot pockets of ultra-diffuse gas (n < 0.1 H cm⁻³ and T ≳ 10⁶ K). The galactic wind exhibits only a two-phase structure with the warm and the hot phases from the ejected ISM.

Fig. 2. Mass-weighted projections of the gas density (top two left panels), temperature (top two middle panels), dust-to-metal mass ratio (top two right panels), and dust-to-gas mass ratio (bottom two left panels), small grain fraction (bottom two middle panels), and carbonaceous grain fraction (bottom two right panels), viewed edge-on (top panel with a depth of 30 kpc) or face-on (bottom panel with a depth of 600 pc) for the high resolution version of G10LG galaxy with the fiducial physics at t = 400 Myr. In black are shown isocontours of density with n = 1 H cm−3.

The dust-to-metal mass ratio DTM = ρ_D/ρ_Z and the dust-to-gas mass ratio DTG = ρ_D/ρ, where ρ_D is the total dust mass density and ρ_Z is the total metal (dust+gas) density, show significant spatial variations in the ISM. The DTM is higher in regions of higher densities due to the increase in accretion rates of refractory elements on dust with gas density, and to the higher destruction rates by thermal sputtering at higher temperatures. Dense gas also corresponds to regions where the mass fraction of small grains f_S is lower and the mass fraction of carbonaceous grains f_C is higher compared to the more diffuse ISM.

We will ignore this aspect in this work, but it can already be noted that the circum-galactic medium has an appreciable level of dust, with the DTM reaching at 0.1 − 0.2 composed mostly of large 0.1 μm grains. Due to the high values of temperature (T ≳ 10⁶ K) reached in the galactic outflow, the smallest populations of dust grains are faster destroyed by thermal sputtering.

4.2. Dust-to-metal ratio

The evolution over time of the DTM for the fiducial G10LG galaxy is shown in Fig. 3. This quantity is computed for all the gas and dust contained in a cylinder of r = 4 kpc cylindrical radius centred on the box and h/2 = 200 pc above and below the disc plane. The DTM rises steeply in less than 50 Myr to its approximately steady-state value at nearly DTM = 0.4 comparable to the canonical ratio in present-day galaxies with solar-like metallicities (e.g., Rémy-Ruyer et al. 2014; De Vis et al. 2019) or in the Milky Way (e.g., Jenkins 2009). The DTM is decomposed into the four simulated dust bins, i.e., into small (dotted) and large (dashed) grains, and into carbonaceous (blue) and silicate (red) grains. The bulk of the DTM is from large grains (80%) and from silicates (66%), and in particular from large silicate grains (53%).

Fig. 3. Dust-to-metal mass ratio (DTM) as a function of time for the fiducial G10LG galaxy. The DTM is further decomposed into the carbonaceous (’C’ in blue) and silicate (’Sil’ in red) grains, and into small 5 nm (’S’ in dotted) and large 0.1 μm (’L’ in dashed). The DTM is close to the canonical value of 0.4 in the Milky Way, where the bulk of the mass is composed of large silicate grains.

Since the accretion time scales inversely with the density of refractory elements in the gas phase, the value of the DTM is sensitive to the ability of the simulation to resolve the dense gas where most of the accretion and growth of dust proceeds. This motivates the use of a subgrid model (as introduced in Sect. 3.4) for accretion if the typical molecular cloud densities are not captured. This is illustrated in Fig. 4 where we compare the G10LG galaxy with the fiducial dust physics run at different spatial resolutions (solid lines) to the same simulated galaxies where the turbulence-driven subgrid dust accretion is turned off (dashed lines for the series of G10LG_LB simulated galaxies). The numerical convergence to DTM = 0.4 is obtained at lower resolution for the fiducial model for unresolved molecular cloud densities than in the naive model, although both give similar values of DTM when sufficiently resolved (few per cent relative difference at a resolution of 9 pc).

Fig. 4. Dust-to-metal mass ratio (DTM) as a function of time for the different resolutions of the G10LG galaxy using either the fiducial model for dust accretion with turbulence-driven for unresolved densities (solid lines) or the model without (dashed lines). From bottom to top are the simulations with Δx = 72, 36, 18, and 9 pc resolution (the 9 pc resolution simulation is for the dust fiducial model only). The subgrid turbulent model for accretion significantly improves the numerical convergence of the DTM.

4.3. Dust size distribution

The size distribution of grains in the Milky Way, n^D(a) the number density of grains at a given size a, can be well represented by a power law with slope n^D ∝ a^−3.5. This is the so-called MRN size distribution from Mathis et al. (1977).

The top panel of Fig. 5 shows the reconstructed size distribution of grains (in black solid line) from our four populations with size a_j (j standing for S or L) and chemical composition i (for carbonaceous or silicate grains), assuming that each population is sampled with a modified log-normal PDF (as in Hirashita 2015 and Hou et al. 2017) of

$$\begin{array}{c}\hfill n_{i,j}^{\mathrm{D}}(a)=\frac{C_{i,j}}{a^4}exp(-\frac{{[ln(a/a_j)]}^2}{2{\sigma }_{i,j}^2}),\end{array}$$(22)

Fig. 5. Respectively top and bottom panels: dust size distribution and extinction curve obtained from the fiducial G10LG simulation at t = 400 Myr for the total dust content (black solid), and the contribution from carbonaceous (blue, ’C’), silicate (red, ’Sil’) grains, and small (dotted, ’S’) and large grains (dashed, ’L’). The extinction curves at different times t = 100, 200, 300, and 400 Myr are shown from light to dark grey scales. The scatter of the extinction curves from the simulation are the brown solid lines. The dust size distribution is compared to that of the Mathis et al. (1977; MRN) size distribution in the Milky Way (in orange). The extinction from the Milky Way (MW), Large Magellanic Cloud (LMC), and Small Magellanic Cloud (SMC) from Pei (1992) are shown as labelled in the corresponding panel, with the typical scatter estimated from the data of Fitzpatrick & Massa (2007). Our simulated MW analogue is in remarkably good agreement with the MRN size distribution and the MW extinction curve.

where σ_i, j = 0.75 is the standard deviation of the log-normal part and C_i, j is the normalisation of the distribution. The normalisation is given by:

$$\begin{array}{c}\hfill \mu m_{\mathrm{p}}{\mathcal{D}}_{i,j}={\displaystyle {\int }_0^{\infty }}\frac{4}{3}\pi a^3{s}_i{n}_{i,j}^{\mathrm{D}}(a)\mathrm{d}a,\end{array}$$(23)

where μ = 1.22 is the mean molecular weight of the gas (assuming, for simplicity, that the dusty gas mainly contributing to the extinction curve is fully neutral). The size distribution of grains measured by 𝒟_i, j = ∑_ΩD_i, j/∑_ΩM_g, over a control volume Ω that is a cylinder of r = 4 kpc radius and h/2 = 0.2 kpc semi-height around the centre of the galaxy, is in good agreement with the expected MRN distribution of grains in the Milky Way, also shown in Fig. 5. It shows that the size distribution of carbonaceous and silicate grains are extremely similar, and that the amount of carbonaceous grains is lower than that of silicate grains (as expected from the decomposition of DTM in Fig. 3).

4.4. Extinction curve

Extinction curve is a key quantity for nearby galaxies including that of the Milky Way (Pei 1992). It is well known (Mathis et al. 1977; Draine & Lee 1984; Weingartner & Draine 2001a) that extinction curves are direct probes of the grain size distribution and of the grain chemical composition. In the particular case of the Milky Way extinction curve, it exhibits in the optical range a pronounced bump at λ = 2175 Å attributed to small carbonaceous grains⁶, and a moderate slope in the ultraviolet, compared to that of the Magellanic Clouds, for which their slopes are steeper and the bumps are erased. Extinction curves, or grain properties, are also central to understanding the emission of distant galaxies, shaped by their attenuation laws. Attenuation curve (not to be confused with the extinction curve) is the result of the combined radiative transfer effects with inhomogeneous stellar emission and dust distribution (Inoue 2005; Narayanan et al. 2018). However, we will keep the investigation of attenuation curves for future work.

The extinction curve A_λ is obtained from the contribution A_λ, i, j from each grain population, which is given by

$$\begin{array}{c}\hfill \frac{A_{\lambda ,i,j}}{N_{\mathrm{H}}}=2.5loge{\displaystyle {\int }_0^{\infty }}\pi a^2{n}_{i,j}^{\mathrm{D}}(a)Q_{\mathrm{ext},i}(a,\lambda )\mathrm{d}a,\end{array}$$(24)

where N_H is the hydrogen column number density and Q_ext, i(a, λ) is the extinction coefficient that is a function of grain size a, wavelength λ and grain composition i. The size distribution is obtained from the log-normal distribution⁷ of Eq. (21). Q_ext, i(a, λ) is obtained by using the Mie theory (Bohren & Huffman 1983) with the optical constants for carbonaceous and silicate grains from Weingartner & Draine (2001a). Bottom panel of Fig. 5 shows the extinction curve A_λ normalised by the extinction A_V in the V-band at λ_V = 0.55 μm, obtained for the fiducial G10LG simulation at time t = 100, 200, 300, and 400 Myr for all the dust in the same region (within a cylinder of r = 4 kpc and h/2 = 0.2 kpc) as used previously for measuring the dust size distribution (top panel of Fig. 5). We also show the scatter of the extinction curve at t = 400 Myr by sampling the corresponding extinction curves for each individual cells within the same region. The extinction curve from the Milky Way analogue G10LG shows an excellent agreement with that observed for the Milky Way (Pei 1992): it shows a similar slope at short wavelengths and the characteristic bump feature at λ = 2175 Å with similar strength. This shape is clearly distinct from the extinction curves from Magellanic Clouds. As expected, the 2175 Å bump is obtained from the small carbonaceous grains (blue dotted line) due to their optical properties, while the UV-to-optical slope is the result of the small silicate grains (red dotted line). It follows immediately that to obtain steeper slopes and a shallower 2175 Å bump (as in the Magellanic Clouds extinction curves), the fraction of small carbonaceous grains must be reduced and the fraction of small silicate grains must be enhanced. This will be investigated further in Sect. 5.

Diversity in the extinction curves is also a significant feature of Milky Way extinction curves with typical relative variations across lines-of-sights of A_λ/A_V of ∼20% at λ⁻¹ = 8 μm⁻¹ (Cardelli et al. 1989; Fitzpatrick & Massa 2007), as can be seen in Fig. 5. Although, extreme care must be taken when comparing directly with observational data (line-of-sight effects), we investigate what are the different extinction curves as a function of the various gas phases of the ISM to get a better indication of what are variations of the dust properties across them. In Fig. 6, we show the extinction curves at t = 400 Myr in G10LG, as relative variations of δ(A_λ/A_V)|_n = (A_λ/A_V)|_n/(A_λ/A_V)−1 at a given n with respect to the global extinction curve A_λ/A_V (as shown in Fig. 5), for different mean gas densities n ∼ 10⁻², 10⁻¹, 1, 10, and 10² H cm⁻³ (sampled in narrow bins of ±0.1 dex around the mean) collected within the same cylinder of r = 4 kpc and h/2 = 0.2 kpc. Gas at densities of n ∼ 10 and 10² H cm⁻³ exhibits extinction curves close to the global value. The variation of extinction curves is not monotonic with gas density. The extremely low value of gas density n ∼ 10⁻² H cm⁻³ corresponds to the larger decrement in the extinction curves: lower bump strength and shallower UV-to-optical slope. This is understood as the consequence of small grain destruction by SN explosions and thermal sputtering. At intermediate densities n ∼ 10⁻¹ and 1 H cm⁻³, the extinction curves approach their global values and stand above the global expectation at 1 H cm⁻³ as a signature of shattering. Shattering transfers mass from large grains to small grains preferentially in the diffuse ISM, reinforcing both the bump strength and the UV slope. The highest gas densities are the closest to the global value as this is where most of the dust mass is concentrated (50% of the total dust mass at n ∼ 10 H cm⁻³, and 30% at n ∼ 10² H cm⁻³). They have lower extinction curves relative to the lower ISM gas density of n ∼ 1 H cm⁻³ due to the efficient coagulation of small grains into large grains. Finally the different sampled densities all have large scatter of the order of the mean values, which is indicative of the intrinsic diversity of extinction curves at a given gas density.

Fig. 6. Relative variations of the mean extinction curve δ(Aλ/AV)|n (solid line) at a given gas density (color-coded as indicated in the panel) with respect to the galactic extinction curve in the G10LG simulation at t = 400 Myr. The dashed lines indicate the 1-σ scatter for each sampled density. The error bars stand for the scatter of the observed extinction curves of the Milky Way from Fitzpatrick & Massa (2007) relative to the mean extinction curve from Pei (1992).

4.5. Depletion factors

The presence of dust in the ISM can be characterised by the depletion of the corresponding atomic elements in the gas phase f_dep(X) = 1 − M_D(X)/M_Z(X), where M_D(X) and M_Z(X) are respectively the dust mass and the total (dust+gas) mass of the element X. We show in Fig. 7 the depletion factor of the various elements as a function of the neutral gas density n_neutral (assuming that gas with temperature below T ≤ 10⁴ K is neutral and fully ionised otherwise) compared to observations. The observational relations are shown with the raw value of Jenkins (2009) as the black dashed lines, and with the recommended renormalisation of their data in black solid line by Zhukovska et al. 2016 (a correction also adopted in Choban et al. 2022) that compensates for the fact that the observed densities are under-estimated since they are averaged along the line-of-sight. We also show the fit provided by Richings et al. (2022) of the data from De Cia et al. (2016), which extends the Jenkins (2009) data to low density lines-of-sight (with or without the rescaling in solid and dashed red lines respectively). For C depletion, we followed Choban et al. (2022), and we reduced the depletion values of Jenkins (2009) by a factor 2 (in that case, we did not rescale in density because the relation is sufficiently shallow) to account for the apparent over-estimate of the C gas-phase abundance from weak compared to strong CII lines (Sofia et al. 2011; Parvathi et al. 2012). We also show the data points from Jenkins (2009; rescaled by a factor of 2) and the data points from Parvathi et al. (2012) to underline the large uncertainties associated to the estimate of the depletion of C.

Fig. 7. Depletion factors as a function of gas density for the G10LG simulation at 400 Myr (black diamonds with standard deviation) compared to the results from Jenkins (2009; black dashed lines) and from De Cia et al. (2016; red dashed lines). The solid lines are the rescaled observational fits following Zhukovska et al. (2016). The purple and orange diamonds are the depletion factors for the G10LG simulation using a pyroxene compound (MgFeSi2O6), and an iron-poor olivine compound (Mg1.5Fe0.5SiO4), respectively, instead of the default olivine compound (MgFeSiO4) for silicates. The blue diamonds stand for the G10LG simulation without CO formation. For C depletion, we also show the data points from Jenkins (2009; triangles) and Parvathi et al. (2012; squares).

As expected, there is more depletion of individual elements in the gas phase (f_dep decreases) with increasing density as a result of larger accretion in the dense gas and of destroyed grains in the diffuse phase. Given the large uncertainties in the reconstruction of the observed gas densities, our depletion values are in broad agreement with the data, except for the Mg which clearly stands out of the data, while the other compounds of silicate grains (Fe, Si, and O) are in the observational ballpark. Iron is the most depleted of the four elements entering the composition of silicates, and, thus, is the limiting element in the accretion for the formation of olivine with that type of iron inclusions (i.e., MgFeSiO₄). Depletion factors of the silicate-bearing elements are naturally extremely sensitive to the exact composition of the silicate. We performed additional simulations where we have assumed a different silicate compound (MgFeSi₂O₆ or Mg_1.5Fe_0.5SiO₄). In the pyroxene compound MgFeSi₂O₆, the mass fraction of Si is 50% larger than in the olivine and Fe is 35% lower, and, thus, Si becomes the limiting element. Indeed, Fig. 7 shows that f_dep(Si) is lower for a G10LG simulation that is performed assuming pyroxene compound instead of olivine, and that Si, then, becomes the lowest of the silicate-bearing elements. Similarly, decreasing the amount of iron inclusions in olivine (or pyroxene) changes the picture: a G10LG simulation assuming iron-poor olivine Mg_1.5Fe_0.5SiO₄ instead of iron-rich olivine MgFeSiO₄ leads to lower values of f_dep(Mg) that makes Mg the new limiting element of silicate grains.

For C depletion, our resulting values are well within the range of observations although the observational data exhibits a large scatter. We also highlight the role of CO formation at high gas densities in limiting the depletion of C elements by grains (see also Choban et al. 2022). If we do not account for CO formation in the limitation of carbonaceous dust growth from ISM accretion, i.e. as for silicate grains, carbonaceous grain growth is only limited by ice mantle coating starting above n = 10⁴ H cm⁻³ (G10LG_NCO simulation), then carbonaceous grains accrete much more efficiently from the ISM, leading to too much depletion of C. Therefore, CO formation indirectly allows for more moderate values of f_dep(C) in better agreement with the data.

4.6. Which process dominates the dust growth?

The evolution of the dust mass in the simulated galaxies obeys different growth and destruction processes. Figure 8 shows the different dust mass variation rates for different processes involved in the dust evolution of the fiducial G10LG run: accretion, thermal sputtering, SN shock destruction, dust released by stars through stellar winds in SN ejecta, astration, coagulation and shattering. All rates are obtained within the whole simulation box, and curves are smoothed over a 5 Myr timescale for sake of readability. Apart from the first initial peak at a few 10 Myr after the start of the simulation, all curves show little variation with time due to the low SFR and no significant evolution in gas-phase metallicity (which could affect the accretion time since t_acc ∝ Z⁻¹). Accretion of the available refractory elements from the ISM onto the small grains is by far the dominant process of dust mass growth compared to dust released in the ejecta of SNe or of stellar winds (as expected for sufficiently metal-rich galaxies, e.g. Popping et al. 2017; Aoyama et al. 2017; Vijayan et al. 2019; Choban et al. 2022), or by accretion onto large grains. Destruction by SNe is a factor of 3 below the mass growth by accretion and is the dominant mechanism for dust mass removal. Thermal sputtering, which is the mechanism by which the dust is destroyed by pockets of hot (T > 10⁶ K) and diffuse (n < 0.01 H cm⁻³) gas, is a factor 2 below the destruction rate by SNe. Finally destruction of the dust by astration is lower by a factor 4 compared to destruction by SNe. Interestingly, the thermal sputtering removes more mass on large grains than on small grains, and SNII destruction removes an equal mass of small and large grains, while their physical timescale for dust destruction is an increasing function of the grain size. Hence, one would naively expect the sputtering to remove more dust mass from small grains. However, since the dust mass is 4 times larger for large grains than small grains (see Fig. 3), this effect is compensated by the available reservoir of dust mass in each bin of grain size.

Fig. 8. Dust mass variation rates as a function of time in the fiducial G10LG simulation for various dust processes: dust growth by gas accretion (‘acc’, black), dust destruction by thermal sputtering (‘spu’, red), by supernovae (‘SNII−’, pink), or by astration (‘ast’, purple), dust released by supernovae (‘SNII+’, orange) or stellar winds (‘SW’, green), and dust mass transfer between small and large grains by coagulation (‘coa’, blue) and vice versa by shattering (‘sha’, cyan).

The coagulation rate is slightly larger than the shattering rate and they both transfer more mass between their two respective dust grain sizes than accretion from the gas phase to the dust phase. Therefore, coagulation and shattering, and their exact balance, play the critical role in the dust size distribution rather than accretion balanced by destruction effects.

At a solar abundance, for which Fe is the limiting element to the growth by accretion of the silicate with olivine compound, the typical accretion times for silicate and carbonaceous grains for gas with n = 30 H cm⁻³, ℳ = 3 (corresponding to the average mass-weighted gas density and turbulent Mach number of the gas in G10LG) are t_acc, Sil ≃ 1.8 Myr and t_acc, C ≃ 9.5 Myr. Compared with the typical free-fall time at this density, that is t_ff ≃ 8.2 Myr as well as the time for the first SNe to explode (5 Myr), it shows that both silicate and carbonaceous grains can grow efficiently by accretion, and why there are larger amount of silicate grains with respect to carbonaceous grains (t_acc, Sil < t_acc, C).

4.7. Dust grain distribution in the multiphase interstellar medium

As shown in Fig. 2, the dust distribution shows a significant structure in the ISM. We investigate further how dust distributes in the ISM by measuring the DTM, the fraction of small grains and the fraction of carbonaceous grains as a function of the gas density and temperature. This is shown in Fig. 9 together with the distribution of gas mass (black contours).

Fig. 9. Dust-to-metal ratio (DTM, top left), dust-to-gas ratio (DTG, top right), fraction of small grains fS (bottom left), and fraction of carbonaceous grains fC (bottom right) as a function of density n and temperature T for the G10LG galaxy at time t = 400 Myr. Black contours are for the mass distribution of gas.

As could naively be expected by the scaling of the dust accretion time with n and T^−1/2, the DTM and the DTGs increase with increasing density and decreasing temperature, with the largest values corresponding to densities of n ≃ 100 cm⁻³ and temperatures of 100 K. We note that the quoted densities are the raw gas density values extracted from the simulation although dust accretion proceeds from unresolved gas densities at larger values driven by turbulence. It is interesting to see that appreciable levels of dust even in the hot component T > 10⁴ K and in particular in the galactic wind component at T ∼ 10⁶ K stays at DTM of ∼0.1.

The size distribution of grains shows a more complex structure. The fraction of small grains f_S is the largest in the cold neutral medium at n ≃ 1 cm⁻³ and T≃ a few 100 K. This corresponds to densities and temperatures where shattering is efficient at reducing the size of large grains, originated by coagulation in the densest star-forming regions. The fraction of small grains decreases in denser and colder regions as the consequence of efficient coagulation of grains in dense gas. Furthermore, the fraction of small grains diminishes as temperature increases from approximately 1000 K. This reduction occurs because of the enhanced efficiency of thermal sputtering of grains and the fact that high temperatures are remnants of regions propelled by SNe, which obliterate small grains.

Finally, the fraction of carbon f_C in Fig. 9 shows very little variation in the gas–temperature diagram as it is typically of ∼10%, except at the largest temperatures (T ≳ 10⁷ K), where f_C increases more significantly. This reflects the more efficient sputtering yields of silicate grains with respect to carbonaceous grains.

4.8. Effect of the metallicity

We now explore how changing the metallicity of the gas affects the dust mass content in the galaxy. Figure 10 (top panel) shows the evolution of the DTM in the G10LG galaxy with different initial gas metallicities, sampling from 0.1, 0.3, 1 and 2 Z_⊙. There is no significant evolution of the gas metallicity for the two highest metallicity bins, while there is an increase to, respectively, 0.3 and 0.5 Z_⊙ for the two lowest metallicity bins. The DTM shows a similar behaviour in all cases: a steep rise in 50–100 Myr that saturates to a steady state value that is larger for larger values of gas metallicity, which is in good agreement with observations (Rémy-Ruyer et al. 2014; De Vis et al. 2019) as shown in Sect. 5.1. This behaviour is naturally obtained by the scaling of dust accretion time with the inverse of the metallicity. However, by comparing the final value of the DTM with the final value of the gas phase metallicity, it is immediate to see that the relation of DTM with Z is not linear. The reason stems in the saturation of the depletion of refractory elements when metallicity is increasing. In particular, silicates are closer to saturation earlier on with respect to carbonaceous grains: there is extreme depletion of iron elements at solar metallicities, which is the limiting element of silicates for olivine, while there are still significant amounts of C in the gas phase (see Fig. 7). Consequently, the fraction of carbonaceous grains (bottom panel of Fig. 10) increases with metallicity.

Fig. 10. Top panel: dust-to-metal ratio (solid lines, left axis) and gas metallicity (dashed lines, right axis) as a function of time for the 4 different initial metallicities of the G10LG galaxy: 0.1, 0.3, 1, and 2 Z⊙ (from bottom to top). Bottom panel: fraction of carbonaceous grains (solid lines) and of small grains (dot-dashed lines). The DTM increases with metallicity, the fraction of carbonaceous grains increases with metallicity and the fraction of small grains decreases with metallicity.

We compare in Fig. 11 the resulting extinction curves for the different metallicities at t = 400 Myr. Simulations with lower metallicities have extinction curves with a steeper UV-to-optical slope. With lower metallicity of the gas, accretion onto grains is less efficient given that t_acc ∝ Z⁻¹. As a result, there is less dust in the dense phase, which also reduces the dust coagulation efficiency as it is directly proportional to the dust density of small grains. Therefore, the dust size distribution becomes further biased towards small grains (see bottom panel of Fig. 10). Combined to the lower fraction of carbonaceous grains, low metallicity leads to an enhancement of the UV part of the extinction curve since there is an increasing fraction of small silicate grains (see in Fig. 5 how small silicate grains contribute to the extinction curve in the UV part).

Fig. 11. Comparison of extinction curves in the G10LG galaxy at time t = 400 Myr for different initial metallicities (Zg, 0 = 0.1, 0.3, 1, and 2 Z⊙ for G10LG_VLZ, G10LG_LZ, G10LG, and G10LG_HZ respectively). The top panel shows the extinction curves and the bottom panel shows the relative variation of the extinction curve with respect to the fiducial simulation G10LG. Galaxies with lower metallicities produce steeper UV-to-optical slopes due lower accretion efficiencies in the dense gas where coagulation proceeds, hence, the fraction of small grains is larger.

4.9. Changing the coagulation and shattering rates

Our model for grain size evolution relies on several aspects that we explore further: coagulation in dense star forming gas and shattering in diffuse medium, which control the shape of the extinction curve by establishing the balance between small and large grains. To explore the effect of coagulation and shattering, we vary their baseline velocity dispersion – or timescale – by a factor of 3. Figure 12 shows the relative variation of the DTM and of the fraction of small grains with respect to the values from the fiducial G10LG simulation from their values averaged between time t = 350 and 400 Myr. As expected, faster (resp. slower) coagulation rates lead to less (resp. more) small grains, and vice versa for the shattering rates. The relative variations in DTM are less sensitive to variation in coagulation of shattering rates compared to variations in the fraction of small grains.

Fig. 12. Relative variations of the DTM (orange) and of the fraction of small grains fs (magenta) for the given run as labelled on the x-axis with respect to the G10LG simulation, decreasing or increasing the coagulation rate (resp. G10LG_cS and G10LG_cF), decreasing or increasing the shattering rate (resp. G10LG_sS and G10LG_sF), increasing both the coagulation and the shattering rates (G10LG_csF), and reducing the SN dust destruction rates (G10LG_sn). Quantities are averaged from 300 to 400 Myr. Variations of the rates of these processes have little impact on the DTM, however, variations of coagulation and shattering rates can significantly affect the fraction of small grains. Relative variations of the DTM (orange) and of the fraction of small grains fs (magenta) for the given run as labelled on the x-axis with respect to the G10LG simulation, decreasing or increasing the coagulation rate (resp. G10LG_cS and G10LG_cF), decreasing or increasing the shattering rate (resp. G10LG_sS and G10LG_sF), increasing both the coagulation and the shattering rates (G10LG_csF), and reducing the SN dust destruction rates (G10LG_sn). Quantities are averaged from 300 to 400 Myr. Variations of the rates of these processes have little impact on the DTM, however, variations of coagulation and shattering rates can significantly affect the fraction of small grains.

Extinction curves could eventually be used to discuss the validity of the different coagulation and shattering models. Since the coagulation and shattering rates change the balance between the amount of small and large grains, the extinction curve should change its UV-to-optical slope and its bump strength. Figure 13 shows the extinction curves for these simulations compared to the fiducial model for coagulation and shattering. As expected the simulations that produce larger fractions of small grains also exhibit a reinforced bump and a steeper UV-to-optical slope (simulations G10LG_cF and G10LG_sS), that seems to be significantly out of the characteristic Milky Way extinction curve. Simulations with larger fractions of large grains (simulations G10LG_cS and G10LG_sF) have a shallower UV-to-optical slope and a weaker bump feature. Although, for these two cases, the bump signature is acceptable compared to the Milky Way signal, the extinction in the UV seems slightly too low.

Fig. 13. Comparison of extinction curves in the G10LG galaxy at time t = 400 Myr for variations in coagulation (red), shattering rates (blue), and in SN dust destruction efficiency for small grains (cyan). An increase (decrease) in the corresponding rate is depicted by a dashed (respectively solid) line. The triple-dot-dashed orange line stands for an increased in both the shattering and the coagulation rates. Top panel shows the extinction curves and the bottom panel shows the relative variation of the extinction curve with respect to the fiducial simulation.

Finally, we investigate the run with both faster coagulation and shattering rates (with an increase by a factor of 3 in their rates). Figures 12 and 13 show that the DTM, size distributions and extinction curves remain largely unaffected by this simultaneous change in coagulation and shattering rates. As can be seen in Fig. 14, the increased coagulation and shattering rates produce larger variations in the extinction curves as a function of density compared to our fiducial model (compare with Fig. 6), which is reflective of the larger variety in the grain size distribution in the ISM. In particular, it is in the diffuse ISM (n ≲ 1 H cm⁻³) that the difference in the fraction of small grains is the largest. Indeed, at n = 1 H cm⁻³, the free fall time is t_ff ≃ 60 Myr, which is comparable to the canonical value of the shattering time (see Eq. (19)) adopted for the fiducial run t_sha ≃ 50 Myr. Increasing the shattering rates by a factor of 3, therefore, allows the diffuse ISM to efficiently shatter the big grains produced in dense clouds before this diffuse phase collapses into dense star-forming clouds.

Fig. 14. As Fig. 6 for the G10LG_csf simulation that has faster coagulation and shattering rates with respect to the fiducial G10LG simulation. There is more diversity in extinction curves across gas density compared to the fiducial run. The error bars stand for the scatter of the observed extinction curves of the Milky Way from Fitzpatrick & Massa (2007) relative to the mean extinction curve from Pei (1992).

4.10. Changing the dust destruction efficiency by supernovae

Since SN explosions are the leading destruction mechanism for dust (compared to thermal sputtering or astration), it is valuable to investigate how changing the SN dust destruction efficiency ε_SN affects the predicted dust distribution in the simulated galaxy. We recall that the SN dust destruction efficiency is stronger for smaller grains (see Sect. 3.2), and is ε_SN, Si(5 nm)≃0.95 and ε_SN, C(5 nm)≃0.86 for respectively small silicate and carbonaceous grains and ε_SN, Si(0.1 μm)≃0.14 and ε_SN, C(0.1 μm)≃0.095 for respectively large silicate and carbonaceous grains. We now impose the same value for the destruction efficiency of small and large grains by reducing the SN dust destruction efficiency of small grains to that of the large grains (simulation ‘G10LG_sn’).

Reducing the SN destruction efficiency for small grains by an order of magnitude still has a limited impact (% increase compared to G10LG) on both the final DTM and the fraction of small grains as can be seen in Fig. 12. This is the result of having an efficient accretion of refractory elements on grains that is mostly limited by the depletion of gas. Consequently, the corresponding extinction curve is not markedly different to that from the fiducial G10LG simulation, with relative variations of less than 5% (see Fig. 13).

5. Low mass galaxies

We now investigate the dust content and the properties of extinctions curves in ten and hundred times lower mass galaxies (respectively called G8 and G9) than the fiducial Milky Way-like galaxy G10.

5.1. Dust mass content

Low mass galaxies are generally younger and more metal-poor than more massive galaxies, typically decreasing by 0.3 dex in gas metallicity per 1 dex in stellar mass, albeit with very large dispersion in metallicities from one galaxy to another at any given mass (e.g., Sánchez et al. 2019; Curti et al. 2020). Therefore, for the lower mass simulations, we focus on lower gas phase metallicities. Figure 15 gives the visual representation of the gas density and the dust-to-metal (DTM) ratio for the low mass gas-rich galaxy, G8HG_VLZ, and for the intermediate mass galaxy, either gas poor, G9LG_LZ, or gas rich, G9HG_LZ. For the gas-rich cases, the gas is clearly multi-phase and strongly clustered into regions of neutral gas with high densities (n > 10 − 100 H cm⁻³) corresponding to low temperatures (T ≃ 100 K) surrounded by a more diffuse ionised medium with intermediate density (n ≃ 0.1 − 1 H cm⁻³) that is warm (T ≃ 10⁴ K), and SN-driven hot pockets of ultra-diffuse gas (n < 0.1 H cm⁻³ and T > 10⁶ K). Compared to the gas density distribution in the gas poor galaxies G9LG_LZ and G10LG, the gas and the corresponding DTM distributions in the gas rich cases appear to be more flocculent and spiral arms pattern are less well ordered. Also it immediately appears that for the gas rich G9HG_LZ case the DTM is lower than for the gas poor G9LG_LZ case.

Fig. 15. Density-weighted projections of the gas density (top two rows), and of the dust-to-metal ratio (bottom two rows) viewed edge-on (first and third rows) or face-on (second and fourth rows) for the G8HG, G9HG and G10LG galaxies from left to right columns. Overplotted black contours encapsulate regions of gas density larger than 1 H cm−3.

In Fig. 16, we show the DTG and the DTM as a function of the gas-phase metallicity of the galaxy measured in a cylinder of radius 2 kpc for G8 and G9, and 4 kpc for G10, with semi-height of 0.2 kpc. This is shown for the different simulated galaxy masses, metallicities, and gas fractions, at 4 different time (t = 100, 200, 300, and 400 Myr). It is compared with the observational relation obtained by Rémy-Ruyer et al. (2014) where they fit the data with a broken-power (we note that the data from De Vis et al. 2017 and also De Vis et al. 2019 show a similar break in DTM with metallicity although with a larger scatter at low metallicities). The values of DTG increase with the gas phase metallicity across all mass ranges with a significant departure from the linear relation corresponding to constant DTM, in agreement with predictions from cosmological simulations (e.g., Popping et al. 2017; Hou et al. 2019; Li et al. 2019; Graziani et al. 2020; Parente et al. 2022; Lewis et al. 2023). The relation is reasonably well reproduced with a nearly linear relation (i.e., constant DTM) around solar metallicities and with a departure from linear around a gas phase metallicity of a few tenths of solar. Galaxies with different initial gas fractions also show different behaviour in the DTM–metallicity relation. For G9HG galaxies (purple squares in Fig. 16), their DTM is nearly constant after 100 Myr, except for the lowest metallicity case, and their gas is continuously enriched as a result of important levels of SFR (≃0.5 M_⊙ yr⁻¹). Gas-poor intermediate mass galaxies G9LG have little evolution in metallicity (low SFR ≃ 0.03 − 0.04 M_⊙ yr⁻¹) – except during the first 100 Myr for the two lower metallicity cases 0.01 Z_⊙ and 0.03 Z_⊙ due to a short initial burst in SFR –, but they strongly increase their DTM. For the lowest mass galaxies G8HG (gas-rich), the two lowest metallicity cases with initial 0.01 Z_⊙ and 0.03 Z_⊙ show both a a strong evolution in metallicity and in DTM, and the higher metallicity cases with initial value of 0.1 and 0.3 Z_⊙ has only an important evolution in DTM. Therefore the various evolution histories due to varying initial metallicities, galaxy mass, or gas fractions contribute to the break in the DTM/DTG–metallicity relation at 0.1 Z_⊙ and to its large dispersion at around this value.

Fig. 16. Dust-to-gas ratio (DTG, top panel) and dust-to-metal (DTM, bottom panel) as a function of the gas phase metallicity for different simulated galaxies (varying galaxy mass and initial gas metallicity) using the fiducial model for dust evolution. Each symbol with connected lines represent a simulation at different times (100, 200, 300 and 400 Myr with increasing metallicities). Squares are for galaxies with high gas (HG) fractions, and circles for galaxies with low gas (LG) fractions. The dashed black line is the linear relation DTG ∝ [O/H], and the solid black lines are the broken power-law from observations of Rémy-Ruyer et al. (2014) for the mean (black thick solid) and the first and third quartiles (thin solid lines). The grey crosses are the observational data from De Vis et al. (2017) with error bars.

The G9 galaxy cases show a different DTM(DTG)–metallicity relation (normalisation) for the gas-poor and the gas-rich galaxies, but also G9HG galaxies can have different DTM at a given metallicity depending on their initial metallicity value. The reason is that for gas-rich G9HG galaxies, there is a strong evolution in metallicity. We assumed Milky Way chemical compositions, while due to the significant release of fresh metals over a few 100 Myr, there is an increase of α elements (Mg, Si, O which are the other silicate-bearing elements) over Fe (see the first and second panel of Fig. 17 for resp. [O/Fe] and [C/Fe]). We recall that Fe was the limiting element of silicate growth for the Milky Way analogue G10LG (see Sect. 4.5), which is also the case here, as can be seen in Fig. 17 (compare the third and fourth panels) due to the increase of [α/Fe]. Therefore, there is a larger reservoir of α elements that are not accreted due to Fe being strongly depleted in the gas phase. The mass of elements locked into dust 1 − f_dep (see third, fourth and fifth panel of Fig. 17) increases with metallicity (in qualitative agreement with the trends observed in SMC, LMC, and MW, see Roman-Duval et al. 2022), as a result of increased accretion rates of refractory elements with metallicity of the ISM. A crucial insight is that the break in the depletion of Fe as a function of metallicity happens at a different value than for the depletion of C (see also Choban et al. 2024), i.e., f_dep = 0.5 is passed at Z ≃ 0.1 Z_⊙ and Z ≃ Z_⊙ for Fe and C, respectively. This disparity might be key to explain variations in extinction curves (Li et al. 2021) from one simulation to another as we will see in the following section.

Fig. 17. Depletion factors for O (top panel), Fe (middle panel) and C (bottom panel) for the different simulations at different times. Similar symbol and color coding as in Fig. 16. Values of abundances and metallicities in SMC and LMC are taken from the compilation of data in Roman-Duval et al. (2022).

5.2. Extinction curves

We now characterise the extinction curves by measuring the UV-to-optical slope and the bump strength at 2175 Å. We follow Salim & Narayanan (2020) and Li et al. (2021), where the slope is given by the ratio of extinction at 1500 Å(A₁₅₀₀) over A_V, and the bump strength by $\mathrm{\Delta }{A}_{2175}/{\stackrel{\sim }{A}}_{2175}$, where ${\stackrel{\sim }{A}}_{2175}=A_{1500}/3+2A_{3000}/3$, and $\mathrm{\Delta }{A}_{2175}=A_{2175}-{\stackrel{\sim }{A}}_{2175}$. The bump versus slope relation is shown in Fig. 18 with the mean values for the global extinction curves of MW, LMC and SMC inferred from the Pei (1992) extinction curves, and some sampled lines-of-sights MW from Fitzpatrick & Massa (2007), and of the LMC and SMC from Gordon et al. (2003). As already shown in Sect. 4, the Milky Way analogue G10LG exhibits excellent agreement with the extinction features of the Milk Way. This is true at all times, in particular for the G10LG simulations starting with metallicities at solar values or above (i.e., for Z_g, 0 = 1 and 2 Z_⊙), and only at the latest times (t = 400 Myr) for simulations starting with metallicities below solar value (i.e., for Z_g, 0 = 0.1 and 0.3 Z_⊙), indicating that the gas phase metallicity is a key property to obtain particular features in the extinction curves of a galaxy.

Fig. 18. Bump strength versus UV-to-optical slope for the different simulations. Arrows give the direction of time sampled at 100, 200, 300 and 400 Myr. The big black points are the average values of the observed MW, LMC and SMC from Pei (1992) and the smaller grey points are values for sampled sight-lines from Fitzpatrick & Massa (2007) for MW, and from Gordon et al. (2003) for SMC and LMC. The left panel is coded by simulation mass (by their color as indicated in the panel), gas fraction (by their symbol as indicated in the panel), and initial metallicity (by their color shade, from dark to light for increasing initial metallicity Zg, 0, which corresponding metallicity can be seen in Figs. 16 or 17), while the right panel is color-coded by the gas metallicity at that given time. Milky Way-like galaxies (G10LG) are naturally attracted to extinction curve features that closely resemble that of the Milky Way, with the closest they get to the true metallicity of the Milky Way the closer they are to the extinction properties of the Milky Way. Lower mass galaxies (G8 and G9) explore a much larger span of the extinction curve features’ space, and can be attributed to their lower metallicity that correlates with the features of the extinction curves.

It is tempting to compare the results of the simulated G9 low mass galaxies to the extinction curves of the Magellanic Clouds. Indeed, the LMC stellar mass and gas mass are estimated at respectively 2 × 10⁹ M_⊙ and 5 × 10⁸ M_⊙ (Kim et al. 1998), similar to the SMC stellar and gas mass of 2 × 10⁹ M_⊙ and 5 × 10⁸ M_⊙ (Stanimirović et al. 2004). Hence, the stellar mass of G9HG is close to those of the Magellanic Clouds, while our adopted gas mass is larger. SMC and LMC share extended streams of gas of ∼5 × 10⁸ M_⊙ (Brüns et al. 2005), which indicates that they have endured important processing of their baryonic content in an interaction. Finally the stellar metallicities in SMC and LMC are respectively [Fe/H]≃ − 1 and ≃ − 0.45 (Cole et al. 2005; Choudhury et al. 2018, 2021). Additionally, it is important to note that there is a large scatter (and uncertainty) in the UV-to-optical slope of sampled lines-of-sights for SMC and LMC, but very little uncertainty in the bump strength (compare with the distribution of points for the MW in Fig. 18). In conclusion, the comparison with the extinction curves of SMC and LMC should be taken with extreme care given the dynamical interaction between the two systems, and the large dispersion in the UV-to-optical slope.

With this important warning in mind, we show the extinction curve features (bump and slope) of the simulated lower mass galaxies, G9HG and G8HG for different initial gas metallicities in Fig. 18 together with the values for the actual Magellanic Clouds (and Milky Way). Those two galaxies show a clear departure from the Milky Way extinction curve with (i) a lower bump strength, and (ii) a steeper slope. For G9HG, the extinction curve features are attracted to values in the broad range of SMC and LMC, in particular for initially low metallicities. Nonetheless, they strongly evolve in time mostly due to the strong evolution in metallicity they have (see the right panel of Fig. 18).

This metallicity effect has some important consequences in terms of the grain size distribution of dust. We recall that SN ejecta only releases large grains, while dust growth by accretion strongly favours the growth of small grains (t_acc ∝ a⁻¹). Consequently, the extinction curve should be flat in the UV-to-optical regime and the bump should be erased in galaxies whenever the growth by stellar ejecta becomes important with respect to the growth by accretion. This effect is illustrated in Fig. 19: the slope and bump for the low mass galaxies are color-coded by $f_{\mathrm{acc},i}={\dot{M}}_{\mathrm{acc},i}/({\dot{M}}_{\mathrm{acc},i}+{\dot{M}}_{\mathrm{ej}}+{\dot{M}}_{\mathrm{coa},i^{\ast }})$ (i^* = C or Sil, when i = Sil or C respectively) and are size-coded either by X_SNII,C = Ṁ_SNII,C/Ṁ_SNII,C + Ṁ_acc,C) or by X_coa,Sil = Ṁ_coa,Sil/(Ṁ_coa,Sil + Ṁ_acc,Sil). The corresponding growth rates are computed by taking the average value from the previous 30 Myr of evolution of the given time to smooth out the high-frequency variations (in particular for the ejection rates). The reason we employ M_coa,i^* and not M_acc,i^* in the calculation f_acc,i is because the damping of the bump (produced by small carbonaceous grains) or of the UV-to-optical slope (produced by small silicate grains) is due mostly to large grains.

Fig. 19. Similar to Fig. 18 where values of the extinction curve features are size-coded either by the SN ejecta XSNII, C growth strength for the case of carbonaceous grains, or by the coagulation growth strength Xcoa, Sil for the case of silicate grains, and are color-coded by the accretion fractions facc, C in blue or facc, Sil in red (see text for details). Values are shown every 25 Myr from 100 Myr to 400 Myr. Extinction curves features for different galaxy type (G8HG, G9HG, G9LG) are shown in different panels (respectively left, middle and right), and each panel contains simulations with different initial metallicities Zg, 0, from low to high initial metallicity moving from the bottom left of each panel to the top/top-right of the panel in anti-clock wise manner.

It appears that simulated galaxies have a low bump strength when both X_SNII, C are high and f_acc, C are low (dark colours), as a natural consequence of the larger fraction of large carbonaceous grains released in stellar ejecta over accretion in the ISM. When the accretion of C atoms onto carbonaceous dust grains is large with respect to the release by ejecta the bump strength becomes larger, and the extinction curves also evolve towards shallower UV-to-optical slopes. This increase in bump strength with X_SNII, C also corresponds to an increase in f_acc, C, i.e. in the fraction of accretion to the total growth of dust (carbonaceous and silicate grains). The increase in the UV-to-optical slope also follows an increase in accretion (f_acc, Sil), which also correspond to an increase in metallicity (compare with Fig. 18) until they get shallower at the highest metallicities. Importantly, one can see that these trajectories in extinction curve features bifurcate for some of these runs (e.g., G8HG_VVVLZ, G8HG_VVLZ, G9HG_VVVLZ, G9HG_VVLZ, G9LG_VVVLZ or G9LG_VVLZ) which corresponds to a significant increase in X_coa, S. For the largest metallicity runs (e.g., G8HG_LZ, G9HG_LZ, or G9LG_LZ) f_acc, Sil is lower compared to lower metallicities, i.e., the relative growth of carbonaceous grains increases (see top panels of Fig. 19), and they also present large relative contribution from coagulation X_coa, Sil, both conspiring to reduce the steepness of the UV-to-optical slope.

The consequences of the relative contribution of ejecta, accretion rates, and coagulation to the dust growth of the different size components can be seen in the fraction of small grains that constitute the dust population. Figure 20 shows the fraction of small silicate f_Sil, S and small carbonaceous f_C, S (relative to the total mass of dust) as a function of the gas metallicity. At very low metallicities (Z ≲ 0.1 Z_⊙), there is an increase of f_Sil, S, that stems from the increased contribution in silicate accretion over ejecta. It is followed by a decrease of f_Sil, S (Z ≳ 0.1 Z_⊙) caused by the increased coagulation rates. Carbonaceous grains have a significantly different behaviour. In qualitative agreement with the increased contribution from accretion (over mass released by ejecta) to the growth of carbonaceous grains (Fig. 19, top panels), the fraction of small carbonaceous grains f_C, S increases by ∼2 dex from ∼0.1% up to 8% from the lowest metallicities ∼0.1 Z_⊙ up to ∼0.5 Z_⊙. Above Z ≳ 0.5 Z_⊙f_C, S starts to saturate with the increasing contribution from coagulation. We also show the mass fraction of polycyclic aromatic hydrocarbons (PAHs) from the observations of nearby galaxies (Rémy-Ruyer et al. 2015), along with the results from our simulations for f_C, S, which show excellent agreement. Although this is reassuring in the ability of our model to capture a population of small carbonaceous structures in agreement with the data, we must stress that PAHs have additional growth and destruction channels to that of small carbonaceous grains (although they might well co-evolve). We will discuss further this aspect in Sect. 6.

Fig. 20. Fraction of small silicate (top) and of small carbonaceous (bottom) grains in the overall dust mass. Similar symbol and color coding as in Fig. 16. The grey crosses are the observational data from Rémy-Ruyer et al. (2015) with error bars for the fraction of dust mass in PAHs.

In summary for low metallicity galaxies, ejecta dominates the growth of dust until they reach sufficiently high metallicity and that the accretion of silicates dominates its growth, increasing the UV-to-optical slope. The behaviour is similar for carbonaceous grains but it happens at higher metallicities (and this can be traced back to the depletion of Fe and C in the gas phase, see Fig. 17), which explains the increase in the bump strength. Finally, the UV-to-optical slope gets shallower at the largest metallicities when both the amount of small carbonaceous grains (grown by accretion) and the coagulation of silicate grains increase.

5.3. The role of CO in the 2175 Å bump

In the dust accretion rate we introduced a limit to the growth of carbonaceous grains argued on the basis of rapid formation of CO molecules at high gas densities (for n > 10³ H cm⁻³, see Sect. 3.4). To explore the impact of CO, we also tested a range of models where we neglect the role of the removal of C atoms into CO molecules: i.e. both carbonaceous and silicate grains stoped their growth by accretion of refractory material at densities n > 10⁴ H cm⁻³ due to icy mantles coating the grain’s surface.

The effect on the depletion of C and on the extinction curve is shown in Fig. 21, resp. top and bottom panels, for several galaxy masses and for different gas fractions. Removing the consideration on CO formation for the subgrid model of carbonaceous grain growth results in a larger depletion of C, i.e., more carbonaceous grains are grown from the gas phase. Consequently, the 2175 Å bump strength of the extinction curve is more prominent as a result of a large amount of small carbonaceous grains due to their stronger growth by accretion, which favours the growth of small grains. In addition, without considering CO formation, the extinction curves are also shallower in their UV-to-optical slope due to the larger ratio of carbonaceous to silicate grains, since there is, proportionally, a lower amount of small silicate grains that are responsible for the steepness of the slope. A striking behaviour with CO modelling and its impact on carbonaceous dust growth is that galaxies with the low gas fractions exhibit weak variation of the amount of (small) carbonaceous grains with respect to the CO modelling, as opposed to simulated galaxies with a high fraction of gas. Indeed galaxies with high gas fractions have much more CO formed compared to galaxies with low gas fractions as they reach high gas surface densities, are more turbulent and, hence, get more gas that can reach typical threshold of CO formation at 10³ H cm⁻³. We can also see that for the simulations with the largest difference in bump strength due to the absence of CO corresponds also a significant shift in the UV-to-optical slope: there is overall more carbonaceous grains with respect to silicate grains, which produces a shallower slope of the extinction curve.

Fig. 21. Bump strength versus UV-to-optical slope of the extinction curve for some selected simulations considering or not the formation of CO as a limiter to carbonaceous grain growth (as indicated in the panel). Assuming capture of C atoms in molecules which are not available for grain growth has the effect of reducing the amount of small carbonaceous grains in the ISM and, hence, the bump strength of the extinction curve.

Figure 22 illustrates this last point by showing the estimated CO gas surface densities Σ_CO for two simulated galaxies of the same mass and with the same initial gas metallicity but which differ in initial gas fraction (G9LG_VLZ versus G9HG_VLZ). Due its larger gas fraction, the G9HG_VLZ run produces much more CO than the G9LG_VLZ run, with CO concentrated in regions of high gas densities. G9HG_VLZ contains ≃60 times more CO than G9LG_VLZ, while the ratio of their total gas mass is only ≃3. Since gas accretion primarily proceeds in dense gas, the effect of C capture in CO molecular formation is a strong limiter to growth of carbonaceous grains, and is, hence, the stronger in galaxies with large gas densities. The estimated CO masses are naive as we only assume instantaneous processing of the entire C gas mass into molecules at a threshold density: it does not take into account the time it takes to grow the CO or destroy it, nor the effect of destruction by UV radiation field and by cosmic rays. We defer such investigation and its impact on carbonaceous grain growth to future studies.

Fig. 22. Estimated CO gas surface densities ΣCO for the gas-poor G9LG_VLZ galaxy (left panels) and the gas-rich G9HG_VLZ galaxy (right panels) at t = 400 Myr. Red are for the isocontours of projected gas density of 10 H cm−3. The gas-rich galaxy produce much more molecular gas mass than the gas-poor galaxy.

6. Caveats

Despite the apparent success in reproducing several key observables of dust properties in our simulated galaxies, our model for dust evolution has made a few important assumptions that we will discuss.

6.1. Polycyclic aromatic hydrocarbons

The presence of multiple emission bands in the mid-IR is well attested and the vibrational excitation of small organic molecules, made of PAHs, is the favourite picture to explain these emission features (Leger & Puget 1984; Desert et al. 1990; Draine & Li 2007; Tielens 2013). Observations of the 3.3/11.3 μm band ratio near photo-dissociation regions suggests that, within the Milky Way ISM, PAHs are primarily composed from molecules with ∼50 − 100 C atoms (Croiset et al. 2016), well below the centroid of our small grain bin. Despite their small size, PAHs can lock up to ∼10% of C in the ISM. They are hypothesised to play a key role in the thermo-chemistry of the ISM, by heating HI gas and the warm phase via the photo-ejection of electrons (e.g., Bakes & Tielens 1994; Weingartner & Draine 2001b), catalysing the production of molecular hydrogen, or acting as the formation sites of small hydrocarbons (see the review by Kaiser et al. 2015, and references therein).

Theoretical modelling of PAH survival (subject to the diffuse ISM radiation field, SN shocks and thermal sputtering) implies that efficient mechanisms for their replenishment are necessary to account for their abundance (Micelotta 2009), may this be either through a ‘top-bottom’ picture via the shattering of small and/or large carbonaceous grains (Seok et al. 2014; Hirashita & Murga 2020; Narayanan et al. 2023), or from combustion-like reaction in the ISM giving rise to the most fundamental aromatic species (e.g., benzene, naphthalene, Reizer et al. 2022). Observations imply a significantly more linear relation of PAH mass with galaxy metallicity (e.g., Madden et al. 2006; Galliano et al. 2021; Shivaei et al. 2024), and a strong correlation of PAH band emission with the interstellar radiation field (e.g., Pilleri et al. 2012; Egorov et al. 2023), supporting that there are ISM processes that alter differently the PAH and the dust grain abundances. Of fundamental relevance for the present work, PAHs are also candidate sources of the 2175 Å feature. However, in our extinction curve modelling, the origin of this feature comes from the graphitic properties assumed for our small carbonaceous grains (5 nm) with the Weingartner & Draine (2001a) optical constants. Therefore, a model in which this feature is originated instead by smaller, less resilient particles (e.g., PAHs), could result in a larger variation (compared to the one presented here) of the 2175 Å bump strength with ISM conditions and galaxy metallicity.

Although the extremely high radiation fields (i.e., intensity in the Habing band of G₀ ≳ 10⁸) necessary for the thermal sublimation of dust grains are only attained near extreme sources like active galactic nuclei or γ-ray bursts (Guhathakurta & Draine 1989), radiation can have a direct effect on the content and properties of the dust grains considered in this work. In particular, if amorphous grains deviate from pure spherical geometry, they present differential scattering and absorption of photons. This means that under an anisotropic radiation field, they will experience an effective radiation torque, capable of accelerating them to the point that their cohesive forces cannot prevent the fragmentation of the grain by centrifugal forces (see Hoang 2020, and references therein). Given that small grains are more easily slowed down due to IR emission and drag from the gas atoms, this radiative-torque disruption (RATD) preferentially fragments large grains, thus altering the shape of the grain size distribution. This may play an important role within low-metallicity, highly star-forming systems, competing with SN shattering/destruction, allowing for steeper UV-slopes than the ones achieved in this work (Fudamoto et al. 2020). We will further explore the effect of RATD in the dust extinction curve in future work (Rodríguez Montero et al., in prep.).

6.2. Dust in stellar ejecta

We have seen that for low metallicity galaxies (≳0.1 Z_⊙), their extinction curves are reflective of their inability to grow carbonaceous dust by gas accretion, hence, their dust grain composition is skewed towards silicates, and carbonaceous dust is only composed of big grains shaped by the SNII ejecta. For extremely low metallicities galaxies (≲0.1 Z_⊙), the gas accretion onto silicate grains is also severely quenched, and SNII ejecta become the only source of both carbonaceous and silicate dust growth in the ISM. Dust condensation efficiencies for SNe are highly uncertain (Schneider & Maiolino 2023). Dust in SNe condenses in the shell of the young ejecta, with dust masses ranging [0.03 − 1] M_⊙ for 20 M_⊙ progenitor (e.g., Sarangi & Cherchneff 2015; Marassi et al. 2019; Brooker et al. 2022), before the destruction of dust takes on by the ignition of the SN reverse shock through sputtering and grain-grain collisions (see e.g., Kirchschlager et al. 2019), which settles the final amount of dust that is going to be mixed with the ISM. How much of the dust survives in this final phase depends on the magnetic field, initial dust size distribution, and clumpiness of the ejecta. Therefore, final dust masses remain largely uncertain, with dust survival rates ranging from a few % to a few 10%, and with large uncertainties in the final dust size distributions (e.g., Silvia et al. 2010; Micelotta et al. 2016; Kirchschlager et al. 2019, 2023; Slavin et al. 2020). On the observational side, dust mass estimates are also affected by strong uncertainties related to dust composition, size distribution and temperature, with best estimates for Cassiopeia A and Crab Nebula that are in the range [0.1 − 1] M_⊙ and [0.025 − 0.05] M_⊙, respectively (see the compilation of data and the discussion in Schneider & Maiolino 2023).

Given all those uncertainties in dust yields, it is difficult to predict accurately the dust composition and, hence, the extinction curves from low mass and low metallicity galaxies. As we have seen in our simulations, the contribution from ejecta is of critical importance in Z ≳ 0.1 Z_⊙, as opposed that of silicate grains in that same metallicity range. To gauge the impact of the uncertainties on dust ejecta, we ran some of the G8HG low mass galaxies with different fraction of small grain size fraction in stellar ejecta f_ej, S in Appendix B. Values of the extinction curve UV-to-optical slope and bump strength are significantly increased with higher f_ej, S as long as the grain growth of a given species is dominated by the contribution from stellar ejecta. If carbonaceous grains are the direct product of SN ejecta in low mass galaxies as suggested by these simulations, the extinction curves of the Magellanic Clouds, and young galaxies of the early Universe (Witstok et al. 2023; Markov et al. 2024) through their attenuation curves, could serve as indirect constraints on mean values of condensation efficiencies in SNII ejecta.

6.3. Approximations in the grain properties

An important foundation of the model at the heart of the obtained grain size distribution, especially for the MRN-like size distribution in the Milky Way analogue, is the velocity dispersion of grains. We assumed, as in other work (e.g., Aoyama et al. 2017; Gjergo et al. 2018; Granato et al. 2021), a unique mean grain velocity dispersion for coagulation and a unique σ_gr − n scaling relation for shattering, although different in magnitude for the small and big grains (for respectively coagulation and shattering rates). However, it is clear that the grain velocity dispersion should vary with the grain properties, such as its charge, and with the ISM properties, such as its density, ionisation, turbulence, and magnetisation (Yan et al. 2004; Moseley et al. 2023). Additionally, the velocity dispersion of grains should imply a grain-gas decoupling for the largest grain sizes at the molecular cloud scale, with a typical diffusion length of L ≃ 5 a_0.1(n/1 cm⁻³)⁻¹(σ_gr/c_s) pc, with the consequence that global shattering rates should increase. We defer such investigations with, using e.g., a diffusion-like approach (Lebreuilly et al. 2019), to future work.

We adopted a coarse two-size grain decomposition, as opposed to decompose fully the size spectrum of grains based on the model of Hirashita (2015; a similar approach that is employed in Hou et al. 2017; Aoyama et al. 2017; Gjergo et al. 2018; Granato et al. 2021). This approach significantly reduces the memory and computational requirements of the model, as it only necessitates 4 bins (2 sizes, 2 chemical compositions) instead of the 2N_bin, where N_bin ≳ 10 used in practice for full size spectrum decomposition (McKinnon et al. 2018; Aoyama et al. 2019; Li et al. 2021). Furthermore, the two-size bin approximation greatly simplifies the treatment of coagulation and shattering that otherwise would have needed cross-coupling between size bins, which makes the computations of coagulation and shattering rates to scale like $N_{\text{bin}}^2$. While a detailed description of fine-sampled grain size distribution associated with accurate models of grain size dispersion could be considered a gold standard, this might come at the expense of spatial resolution, which is crucial for capturing the multiphase ISM responsible for phase transitions between accretion/coagulation-dominated regions (in dense neutral or molecular gas) and shattering-dominated regions (in diffuse ionised gas).

7. Conclusion

We introduced a novel model of dust growth for galaxy simulations in the RAMSES code that was coupled to the chemical enrichment of the gas. This model decomposed dust into carbonaceous and silicate grains with two grain sizes of 5 nm and 0.1 μm. It included the modelling of dust growth by accretion, stellar ejecta, dust destruction by SNe, thermal sputtering and astration, and also included coagulation and sputtering to allow dust to cycle between small and large grain sizes. Using hydrodynamical simulations of idealised disc galaxies with a multiphase ISM structure captured at a resolution of 20 pc, we explored various galaxy masses, metallicities, gas fractions, and their effect on the galactic dust properties and extinction curves.

These simulations are able to reproduce the main dust properties of the Milky Way, and they produce dust–metallicity scaling relations in good agreement with observational constraints. They draw a scenario for the physical origin of the observed differences between the extinction curve of the Milky Way and those of the Magellanic Clouds, that stem in the emergence of different dominant processes for dust growth in the ISM (ejecta, accretion and coagulation) driven by the metallicity of the gas. In more details, our main results are the following:

• The numerical analogue of the Milky Way shows excellent agreement with the DTM, the MRN size distribution, the extinction curve, and the depletion factors of various elements of the Milky Way.

• Accretion is the leading growth mechanism in Milky Way-like of galaxies, and the DTM is sensitive to the gas metallicity which controls the amount of accretion onto dust grains.

• In Milky Way-like galaxies, coagulation and shattering play an important role to shape the most prominent features of the extinction curve: its bump and UV-to-optical slope.

• DTM is sensitive to spatial resolution effects but, as introduced here, a subgrid model that captures the unresolved density contrasts due to turbulence greatly improves the convergence properties of the accretion scheme.

• Varying the galaxy mass, gas fraction and metallicity has some impact on the DTM/DTG, and the DTM/DTG-Z relation is naturally reproduced in our model with a rapid decrease in DTM below 0.1 Z_⊙ as suggested by observations. This transition is due to a dominant growth mode by ejecta, at low metallicities, to growth by accretion, at high metallicities.

• A crucial insight is that the break in the depletion of silicate-bearing elements as a function of metallicity happens at a different value than that of C: at Z ≃ 0.1 Z_⊙ and Z ≃ Z_⊙ respectively. This offset, driven by the different accretion efficiencies of silicate and carbonaceous grains, is key to explain variations in extinction curves for low mass and low metallicity galaxies, such as for the Magellanic Clouds.

• By exploring different galaxy masses, gas fractions and metallicities, we obtain diversity in the bump strength and UV-to-optical slope of extinction curves. In particular, these features are strongly correlated with the metallicity of the gas, as it controls the relative amount of growth by accretion (small grains) relative to ejecta (large grains) and to coagulation.

• Due to their lower metallicity, we can interpret the steeper slopes and the weaker bumps of the extinction curves of the Magellanic Clouds with respect to that of the Milky Way as the result of a lower accretion of metals onto dust relative to ejecta (especially for C), of a more efficient accretion of silicate elements with respect to C, as well as of a less efficient coagulation of grains.

• Because of the capture of C atoms in CO molecules, there is a strong suppression of small carbonaceous grains in gas-rich galaxies, which guarantees the good agreement with the extinction curves of the Magellanic Clouds.

Those simulations were carried out in an idealised context without considering the effect of fresh gas replenishment by halo-scale gas accretion. A natural direct extension of this work is to test the validity of the dust model in fully cosmological hydrodynamical simulations to understand the influence of cosmic infall and mergers on the properties of dust. In particular, such cosmological simulations could be used to understand the cosmological build up of dust in galaxies, and test the fraction of galaxies with obscured star formation (Magnelli et al. 2020; Pozzi et al. 2021; Fudamoto et al. 2021; Inami et al. 2022). In this work, we limited ourselves to testing extinction curves of simulated galaxies, which gives a direct insight on dust properties. Although attenuation curves are the result of the complex interplay between the dust properties, and the ISM and stellar distributions (Narayanan et al. 2018), the flurry of observational evidence might serve as an important testbed for the cosmic evolution of dust properties.

The dust model introduced in this paper is currently extended to track PAHs, as well as to include a self-consistent coupling of the dust and ISM non-equilibrium thermo-chemistry in RAMSES-RTZ (Katz 2022) with radiation-hydrodynamics (Rodríguez Montero, et al., in prep.). This ‘Dusty-PRISM’ model would be capable of exploring the influence of radiation, chemistry, gas properties and dust grains on the properties of PAHs, shedding light on the origin of the correlations between mid-IR bands and the properties of the ISM across cosmic time.

Acknowledgments

We thank Hakim Atek, Benoît Commerçon, Pierre Cox, Frédéric Galliano, Pierre Guillard, Patrick Hennebelle, Harley Katz, Guillaume Laibe, Ugo Lebreuilly, Guillaume Pineau des Forêts, and Nicolas Prantzos for useful discussions. This work has made use of the Infinity cluster on which the simulation was run and post-processed, hosted by the Institut d’Astrophysique de Paris. We warmly thank S. Rouberol for running it smoothly. FRM is supported by the Wolfson Harrison UK Research Council Physics Scholarship. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 693024). MT acknowledges support from the NWO grant 016.VIDI.189.162 (“ODIN”). SKY acknowledges support from the Korean National Research Foundation (2020R1A2C3003769, 2022R1A6A1A03053472). Part of the interpretation was based on the calculations using KIST’s Nurion (KSC-2022-CRE-0088, KSC-2022-CRE-0344, KSC-2023-CRE-0343).

References

Appendix A: Released metallicities and dust

The time evolution of the resulting chemical and dust yields for a simple stellar population of different initial metallicities Z_ZAS = 10⁻³,10⁻²,10⁻¹,1 Z_⊙ are shown in Fig. A.1 with assumptions for stellar tracks, IMF, stellar rotations, and dust condensation efficiencies as described in sections 2.4 and 3.1. It shows that the total mass return, and that of O, Mg, Si, S and Ne are overall largely dominated by SNII. N mass release is largely dominated by AGB winds. C has an important contribution from AGB while still dominated by SNII, and, finally, Fe has an equal mass contribution from SNIa and SNII. For dust release, we see as much carbonaceous dust released by AGN as by SNII, while silicate dust is mostly produced by SNII.

Fig. A.1. Evolution of the SSP yield adopted in this work for total metallicity and different individual elements. The blue lines are for the total elemental release (gas plus dust). Red lines are the contribution from intermediate stars. Purple lines are the contribution from SNIa. Green lines are for dust. Lines goes from dark to light for an initial SSP metallicity of ZZAS = 10−3,10−2,10−1,1 Z⊙.

Appendix B: Ejecta sizes in low mass galaxies

We show in Fig. B.1, the changes in the UV-to-optical slope and bump strength of the extinction curves of the G8HG_VVLZ and G8HG_VLZ for different values of small grain size fraction in stellar ejecta f_ej, S = 0 (default value), 0.5, and 1. G8HG_VVLZ has been selected for its important contribution of SNII ejecta component (as opposed to accretion) for both carbonaceous and silicate grains in the fiducial model, and G8HG_VLZ because of its significant SNII ejecta component for carbonaceous grains (only). As expected, these two low mass and low metallicity galaxies are very sensitive to the grain size mixture of their stellar ejecta for both the UV-to-optical slope and bump strength for G8_VVLZ, but only significantly for the bump strength for G8_VLZ since G8_VLZ has silicate grain (mainly responsible for the UV-to-optical slope) growth dominated by ISM accretion.

Fig. B.1. Bump strength and UV-to-optical slope of the extinction curve for the low mass galaxies with initial gas metallicity 0.03 Z⊙ (G8_VVLZ) and 0.1 Z⊙ (G8_VLZ) for different fractions of small grain mass in stellar ejecta fej, S as indicated in the panel. Values of the extinction curve features for the Milky Way (‘MW’), the Large Magellanic Cloud (‘LMC’), and the Small Magellanic Cloud (‘SMC’) are also shown.

Appendix C: Sputtering times

Figure C.1 compares the sputtering times as a function of the gas temperature for a grain size of a = 0.1 μm and a gas density of n = 1 H cm⁻³ from Tielens et al. (1994), Tsai & Mathews (1995), and Hu et al. (2019). Hu et al. (2019) provides a convenient fifth-order polynomial fit to the predictions of the sputtering rates Y(T) from Nozawa et al. (2006). The values of the Hu et al. (2019) fit with sufficient digits are reported in table C.1 (private communication with C.-Y. Hu). The formula from Tsai & Mathews (1995) is intermediate to the sputtering times of carbonaceous and silicate grains from Tielens et al. (1994) and Hu et al. (2019), which both show a general good agreement up to a temperature of T ∼ 10⁸ K. We note that these calculations are all for a semi-infinite target, however, the finite size of the grain modifies the sputtering rates for grains with size smaller than the penetration depth of the impinging ion (e.g. Kirchschlager et al. 2019).

Fig. C.1. Sputtering time as a function of temperature for a gas density of n = 1 H cm−3 and a grain size of a = 0.1 μm. In red and blue solid and dotted are the results of Tielens et al. (1994) for respectively Si and C grains using their fitting functions valid up to T ≃ 2 × 108 K. The black line is the fitting formula from Tsai & Mathews (1995). The red and blue dashed lines are the fitting formulas for respectively Si and C grains from Hu et al. (2019).

Appendix D: Sticking coefficient

There are various predictions for the sticking coefficient α, which value depends on the temperature of the gas (and grains), e.g.:

1. α_C(T) = 0.95(1+βT/T₀)/(1+T/T₀)^β, β = 2.22 and T₀ = 56 K measured from laboratory experiments for chemisorption of H₂ molecules on silicate surfaces (Chaabouni et al. 2012).

2. α_LB(T) = (1+10⁻⁴T^1.5)⁻¹ from theoretical considerations from Le Bourlot et al. (2012).

3. α_L(T) = 1.9 × 10⁻²T(1.7 × 10⁻³T_D+0.4)exp(−7 × 10⁻³T) taken from the fit to the data from quantum calculations of Leitch-Devlin & Williams (1985) in Grassi et al. (2014), and where T_D is the dust temperature. The dust temperature for big grains is at local thermodynamical equilibrium around 20 K for local galaxies (e.g. Schreiber et al. 2018), however very small grains can be stochastically heated (Draine 2003) to large temperatures of about 1000 K.

We present in Fig. D.1 the behaviour of each of these predictions of the behaviour of sticking coefficient with T. They all show an increase of α with a decrease in temperature down to T ∼ 200 K, where α ranges from 0.3 to 0.9 at this temperature. Below this temperature, the values from Leitch-Devlin & Williams (1985) predicts a decrease of α with a decrease of T, while the values from Le Bourlot et al. (2012) and Chaabouni et al. (2012) show an increase of α with a decrease of T up to α = 1.

Fig. D.1. Sticking coefficient for gas accretion onto dust as a function of the gas temperature T for different models: Le Bourlot et al. (2012) in black, Chaabouni et al. (2012) in blue, and Leitch-Devlin & Williams (1985) in red for different dust temperatures TD = 10,20, 40, 80, 300 K for the dotted, solid, dashed, dot-dashed, and triple dot-dashed red lines respectively. The magenta line at α = 1/3 is our preferred value for unresolved gas densities (see Section 3.4).

Since the exact value and rate of change with temperature of the sticking coefficient α remains highly uncertain (see the discussion in Zhukovska et al. 2016), we adopted a constant value of α = 1/3 for our unresolved high-gas density regime driven by turbulence (see Section 3.4). For comparison, the turbulence-driven subgrid model for dust accretion is activated for λ_J < 4Δx and the effective temperature-density relation in our simulations is T ≃ 10⁴(n/0.1 H cm⁻³)^−/3 (see Fig. 9) for n > 0.1 H cm⁻³ and T < 10⁴ K, hence, the temperature at which the subgrid model is triggered for a fiducial resolution of Δx = 20 pc is T ≃ 750 K. This temperature gives a value of α from Le Bourlot et al. (2012) that is close to our adopted constant value of α in the subgrid model. We ran a G10LG simulation with an increased value of α = 1 in that regime: the DTM increases by 5%, mostly driven by an increase in the fraction of carbonaceous grains, which goes from a fiducial value of f_C = 0.34 to f_C = 0.38, and produces a significant decrease of the depletion factor of C by a factor of 2 at n ≃ 10 H cm⁻³.

The turbulence-driven subgrid model for gas accretion onto dust grains has a maximum density above which the accretion of new refractory elements is impossible (either n_max = 10³ or 10⁴ H cm⁻³ (see section 3.4). This value modifies the behaviour of the effective accretion timescale as a function of the turbulent Mach number and is shown in Fig. D.2 for the case of accretion onto carbonaceous dust grains. The effective timescale for accretion t_acc, eff exhibits a decrease with turbulent Mach number above ℳ > 1 until it reaches a minimum for strong ℳ ≃ 5 − 10 after which the effective timescale increases. The exact behaviour is driven by the maximum value n_max. To guide the eye, we show the corresponding free-fall time of the gas at the given mean gas densities $\overline{n}$, and the (single) time delay for SNII adopted in this work.

Fig. D.2. Effective accretion rate of refractory materials onto grains (here for carbonaceous grains) as a function of the turbulent Mach number for two different mean density $\overline{n}=1$ and 100 H cm−3 (respectively in solid blue and in dashed blue). The various colors stand for different density cut-offs nmax = 103 and 1004 H cm−3 (respectively in dark and light blue). For comparison we show the free-fall time corresponding to the quoted gas density $\overline{n}$ in black, and the SN time delay in red dot-dashed.

Appendix E: Sizes and extinction curves for low mass galaxies

In this appendix we show in Figs. E.1 and E.2 the dust grain size distribution and the extinction curves of the G9HG, G9LG, and G8HG listed in table 2.

Fig. E.1. Size distribution of grains and extinction curves at different times for different simulations as indicated in the panel. The colors (red and blue) indicate the contribution from different chemical compositions of the dust (resp. silicate and carbonaceous grains), and we only show this decomposition in the extinction curve for the final time of the simulation.

Fig. E.2. Continued from Fig. E.1.

All Tables

All Figures

	Fig. 1. Mock image of the fiducial G10LG simulation seen edge-on (top) or face-on (bottom) in SDSS g − r − i filter bands at t = 400 Myr.
In the text

Fig. 2. Mass-weighted projections of the gas density (top two left panels), temperature (top two middle panels), dust-to-metal mass ratio (top two right panels), and dust-to-gas mass ratio (bottom two left panels), small grain fraction (bottom two middle panels), and carbonaceous grain fraction (bottom two right panels), viewed edge-on (top panel with a depth of 30 kpc) or face-on (bottom panel with a depth of 600 pc) for the high resolution version of G10LG galaxy with the fiducial physics at t = 400 Myr. In black are shown isocontours of density with n = 1 H cm−3.

In the text

	Fig. 3. Dust-to-metal mass ratio (DTM) as a function of time for the fiducial G10LG galaxy. The DTM is further decomposed into the carbonaceous (’C’ in blue) and silicate (’Sil’ in red) grains, and into small 5 nm (’S’ in dotted) and large 0.1 μm (’L’ in dashed). The DTM is close to the canonical value of 0.4 in the Milky Way, where the bulk of the mass is composed of large silicate grains.
In the text

Fig. 4. Dust-to-metal mass ratio (DTM) as a function of time for the different resolutions of the G10LG galaxy using either the fiducial model for dust accretion with turbulence-driven for unresolved densities (solid lines) or the model without (dashed lines). From bottom to top are the simulations with Δx = 72, 36, 18, and 9 pc resolution (the 9 pc resolution simulation is for the dust fiducial model only). The subgrid turbulent model for accretion significantly improves the numerical convergence of the DTM.

In the text

Fig. 5. Respectively top and bottom panels: dust size distribution and extinction curve obtained from the fiducial G10LG simulation at t = 400 Myr for the total dust content (black solid), and the contribution from carbonaceous (blue, ’C’), silicate (red, ’Sil’) grains, and small (dotted, ’S’) and large grains (dashed, ’L’). The extinction curves at different times t = 100, 200, 300, and 400 Myr are shown from light to dark grey scales. The scatter of the extinction curves from the simulation are the brown solid lines. The dust size distribution is compared to that of the Mathis et al. (1977; MRN) size distribution in the Milky Way (in orange). The extinction from the Milky Way (MW), Large Magellanic Cloud (LMC), and Small Magellanic Cloud (SMC) from Pei (1992) are shown as labelled in the corresponding panel, with the typical scatter estimated from the data of Fitzpatrick & Massa (2007). Our simulated MW analogue is in remarkably good agreement with the MRN size distribution and the MW extinction curve.

In the text

Fig. 6. Relative variations of the mean extinction curve δ(Aλ/AV)|n (solid line) at a given gas density (color-coded as indicated in the panel) with respect to the galactic extinction curve in the G10LG simulation at t = 400 Myr. The dashed lines indicate the 1-σ scatter for each sampled density. The error bars stand for the scatter of the observed extinction curves of the Milky Way from Fitzpatrick & Massa (2007) relative to the mean extinction curve from Pei (1992).

In the text

Fig. 7. Depletion factors as a function of gas density for the G10LG simulation at 400 Myr (black diamonds with standard deviation) compared to the results from Jenkins (2009; black dashed lines) and from De Cia et al. (2016; red dashed lines). The solid lines are the rescaled observational fits following Zhukovska et al. (2016). The purple and orange diamonds are the depletion factors for the G10LG simulation using a pyroxene compound (MgFeSi2O6), and an iron-poor olivine compound (Mg1.5Fe0.5SiO4), respectively, instead of the default olivine compound (MgFeSiO4) for silicates. The blue diamonds stand for the G10LG simulation without CO formation. For C depletion, we also show the data points from Jenkins (2009; triangles) and Parvathi et al. (2012; squares).

In the text

Fig. 8. Dust mass variation rates as a function of time in the fiducial G10LG simulation for various dust processes: dust growth by gas accretion (‘acc’, black), dust destruction by thermal sputtering (‘spu’, red), by supernovae (‘SNII−’, pink), or by astration (‘ast’, purple), dust released by supernovae (‘SNII+’, orange) or stellar winds (‘SW’, green), and dust mass transfer between small and large grains by coagulation (‘coa’, blue) and vice versa by shattering (‘sha’, cyan).

In the text

	Fig. 9. Dust-to-metal ratio (DTM, top left), dust-to-gas ratio (DTG, top right), fraction of small grains fS (bottom left), and fraction of carbonaceous grains fC (bottom right) as a function of density n and temperature T for the G10LG galaxy at time t = 400 Myr. Black contours are for the mass distribution of gas.
In the text

Fig. 10. Top panel: dust-to-metal ratio (solid lines, left axis) and gas metallicity (dashed lines, right axis) as a function of time for the 4 different initial metallicities of the G10LG galaxy: 0.1, 0.3, 1, and 2 Z⊙ (from bottom to top). Bottom panel: fraction of carbonaceous grains (solid lines) and of small grains (dot-dashed lines). The DTM increases with metallicity, the fraction of carbonaceous grains increases with metallicity and the fraction of small grains decreases with metallicity.

In the text

Fig. 11. Comparison of extinction curves in the G10LG galaxy at time t = 400 Myr for different initial metallicities (Zg, 0 = 0.1, 0.3, 1, and 2 Z⊙ for G10LG_VLZ, G10LG_LZ, G10LG, and G10LG_HZ respectively). The top panel shows the extinction curves and the bottom panel shows the relative variation of the extinction curve with respect to the fiducial simulation G10LG. Galaxies with lower metallicities produce steeper UV-to-optical slopes due lower accretion efficiencies in the dense gas where coagulation proceeds, hence, the fraction of small grains is larger.

In the text

Fig. 12. Relative variations of the DTM (orange) and of the fraction of small grains fs (magenta) for the given run as labelled on the x-axis with respect to the G10LG simulation, decreasing or increasing the coagulation rate (resp. G10LG_cS and G10LG_cF), decreasing or increasing the shattering rate (resp. G10LG_sS and G10LG_sF), increasing both the coagulation and the shattering rates (G10LG_csF), and reducing the SN dust destruction rates (G10LG_sn). Quantities are averaged from 300 to 400 Myr. Variations of the rates of these processes have little impact on the DTM, however, variations of coagulation and shattering rates can significantly affect the fraction of small grains. Relative variations of the DTM (orange) and of the fraction of small grains fs (magenta) for the given run as labelled on the x-axis with respect to the G10LG simulation, decreasing or increasing the coagulation rate (resp. G10LG_cS and G10LG_cF), decreasing or increasing the shattering rate (resp. G10LG_sS and G10LG_sF), increasing both the coagulation and the shattering rates (G10LG_csF), and reducing the SN dust destruction rates (G10LG_sn). Quantities are averaged from 300 to 400 Myr. Variations of the rates of these processes have little impact on the DTM, however, variations of coagulation and shattering rates can significantly affect the fraction of small grains.

In the text

Fig. 13. Comparison of extinction curves in the G10LG galaxy at time t = 400 Myr for variations in coagulation (red), shattering rates (blue), and in SN dust destruction efficiency for small grains (cyan). An increase (decrease) in the corresponding rate is depicted by a dashed (respectively solid) line. The triple-dot-dashed orange line stands for an increased in both the shattering and the coagulation rates. Top panel shows the extinction curves and the bottom panel shows the relative variation of the extinction curve with respect to the fiducial simulation.

In the text

Fig. 14. As Fig. 6 for the G10LG_csf simulation that has faster coagulation and shattering rates with respect to the fiducial G10LG simulation. There is more diversity in extinction curves across gas density compared to the fiducial run. The error bars stand for the scatter of the observed extinction curves of the Milky Way from Fitzpatrick & Massa (2007) relative to the mean extinction curve from Pei (1992).

In the text

	Fig. 15. Density-weighted projections of the gas density (top two rows), and of the dust-to-metal ratio (bottom two rows) viewed edge-on (first and third rows) or face-on (second and fourth rows) for the G8HG, G9HG and G10LG galaxies from left to right columns. Overplotted black contours encapsulate regions of gas density larger than 1 H cm−3.
In the text

Fig. 16. Dust-to-gas ratio (DTG, top panel) and dust-to-metal (DTM, bottom panel) as a function of the gas phase metallicity for different simulated galaxies (varying galaxy mass and initial gas metallicity) using the fiducial model for dust evolution. Each symbol with connected lines represent a simulation at different times (100, 200, 300 and 400 Myr with increasing metallicities). Squares are for galaxies with high gas (HG) fractions, and circles for galaxies with low gas (LG) fractions. The dashed black line is the linear relation DTG ∝ [O/H], and the solid black lines are the broken power-law from observations of Rémy-Ruyer et al. (2014) for the mean (black thick solid) and the first and third quartiles (thin solid lines). The grey crosses are the observational data from De Vis et al. (2017) with error bars.

In the text

	Fig. 17. Depletion factors for O (top panel), Fe (middle panel) and C (bottom panel) for the different simulations at different times. Similar symbol and color coding as in Fig. 16. Values of abundances and metallicities in SMC and LMC are taken from the compilation of data in Roman-Duval et al. (2022).
In the text

Fig. 18. Bump strength versus UV-to-optical slope for the different simulations. Arrows give the direction of time sampled at 100, 200, 300 and 400 Myr. The big black points are the average values of the observed MW, LMC and SMC from Pei (1992) and the smaller grey points are values for sampled sight-lines from Fitzpatrick & Massa (2007) for MW, and from Gordon et al. (2003) for SMC and LMC. The left panel is coded by simulation mass (by their color as indicated in the panel), gas fraction (by their symbol as indicated in the panel), and initial metallicity (by their color shade, from dark to light for increasing initial metallicity Zg, 0, which corresponding metallicity can be seen in Figs. 16 or 17), while the right panel is color-coded by the gas metallicity at that given time. Milky Way-like galaxies (G10LG) are naturally attracted to extinction curve features that closely resemble that of the Milky Way, with the closest they get to the true metallicity of the Milky Way the closer they are to the extinction properties of the Milky Way. Lower mass galaxies (G8 and G9) explore a much larger span of the extinction curve features’ space, and can be attributed to their lower metallicity that correlates with the features of the extinction curves.

In the text

Fig. 19. Similar to Fig. 18 where values of the extinction curve features are size-coded either by the SN ejecta XSNII, C growth strength for the case of carbonaceous grains, or by the coagulation growth strength Xcoa, Sil for the case of silicate grains, and are color-coded by the accretion fractions facc, C in blue or facc, Sil in red (see text for details). Values are shown every 25 Myr from 100 Myr to 400 Myr. Extinction curves features for different galaxy type (G8HG, G9HG, G9LG) are shown in different panels (respectively left, middle and right), and each panel contains simulations with different initial metallicities Zg, 0, from low to high initial metallicity moving from the bottom left of each panel to the top/top-right of the panel in anti-clock wise manner.

In the text

	Fig. 20. Fraction of small silicate (top) and of small carbonaceous (bottom) grains in the overall dust mass. Similar symbol and color coding as in Fig. 16. The grey crosses are the observational data from Rémy-Ruyer et al. (2015) with error bars for the fraction of dust mass in PAHs.
In the text

Fig. 21. Bump strength versus UV-to-optical slope of the extinction curve for some selected simulations considering or not the formation of CO as a limiter to carbonaceous grain growth (as indicated in the panel). Assuming capture of C atoms in molecules which are not available for grain growth has the effect of reducing the amount of small carbonaceous grains in the ISM and, hence, the bump strength of the extinction curve.

In the text

	Fig. 22. Estimated CO gas surface densities ΣCO for the gas-poor G9LG_VLZ galaxy (left panels) and the gas-rich G9HG_VLZ galaxy (right panels) at t = 400 Myr. Red are for the isocontours of projected gas density of 10 H cm−3. The gas-rich galaxy produce much more molecular gas mass than the gas-poor galaxy.
In the text

Fig. A.1. Evolution of the SSP yield adopted in this work for total metallicity and different individual elements. The blue lines are for the total elemental release (gas plus dust). Red lines are the contribution from intermediate stars. Purple lines are the contribution from SNIa. Green lines are for dust. Lines goes from dark to light for an initial SSP metallicity of ZZAS = 10−3,10−2,10−1,1 Z⊙.

In the text

Fig. B.1. Bump strength and UV-to-optical slope of the extinction curve for the low mass galaxies with initial gas metallicity 0.03 Z⊙ (G8_VVLZ) and 0.1 Z⊙ (G8_VLZ) for different fractions of small grain mass in stellar ejecta fej, S as indicated in the panel. Values of the extinction curve features for the Milky Way (‘MW’), the Large Magellanic Cloud (‘LMC’), and the Small Magellanic Cloud (‘SMC’) are also shown.

In the text

Fig. C.1. Sputtering time as a function of temperature for a gas density of n = 1 H cm−3 and a grain size of a = 0.1 μm. In red and blue solid and dotted are the results of Tielens et al. (1994) for respectively Si and C grains using their fitting functions valid up to T ≃ 2 × 108 K. The black line is the fitting formula from Tsai & Mathews (1995). The red and blue dashed lines are the fitting formulas for respectively Si and C grains from Hu et al. (2019).

In the text

Fig. D.1. Sticking coefficient for gas accretion onto dust as a function of the gas temperature T for different models: Le Bourlot et al. (2012) in black, Chaabouni et al. (2012) in blue, and Leitch-Devlin & Williams (1985) in red for different dust temperatures TD = 10,20, 40, 80, 300 K for the dotted, solid, dashed, dot-dashed, and triple dot-dashed red lines respectively. The magenta line at α = 1/3 is our preferred value for unresolved gas densities (see Section 3.4).

In the text

Fig. D.2. Effective accretion rate of refractory materials onto grains (here for carbonaceous grains) as a function of the turbulent Mach number for two different mean density $\overline{n}=1$ and 100 H cm−3 (respectively in solid blue and in dashed blue). The various colors stand for different density cut-offs nmax = 103 and 1004 H cm−3 (respectively in dark and light blue). For comparison we show the free-fall time corresponding to the quoted gas density $\overline{n}$ in black, and the SN time delay in red dot-dashed.

In the text

	Fig. E.1. Size distribution of grains and extinction curves at different times for different simulations as indicated in the panel. The colors (red and blue) indicate the contribution from different chemical compositions of the dust (resp. silicate and carbonaceous grains), and we only show this decomposition in the extinction curve for the final time of the simulation.
In the text

	Fig. E.2. Continued from Fig. E.1.
In the text