Evolution of the grain size distribution in Milky Way-like galaxies in post-processed IllustrisTNG simulations

This model treats the galaxy as a one-zone object, in which dust enrichment is calculated in a consistent manner with metal enrichment. We apply it to individual subhalos. The grain size distribution at time $t$ is denoted as $n(a,\,t)$ such that $n(a,\,t)\,\mathrm{d}a$ is the number density of dust grains per hydrogen atom in a grain radius range of $(a,\,a+\mathrm{d}a)$. We solve the evolution of $n(a,\,t)$ with an initial condition of $n(a,\,t=0)=0$ (no dust). The total dust abundance is traced by the dust-to-gas ratio, which is evaluated by

We consider the following five processes for dust evolution: stellar dust production (dust condensation in stellar ejecta), dust destruction in SN shocks in the ISM, dust growth by the accretion of gas-phase metals and by coagulation in the dense ISM, and dust disruption by shattering in the diffuse ISM. The evolution of the grain size distribution is calculated for each of the above processes as summarized below.

The stellar dust production is based on the metallicity evolution given by TNG. We estimate the increase of dust-to-gas ratio for a galaxy between two successive simulation snapshots from the increment of metallicity, by assuming the metal-to-dust condensation fraction, denoted as $f_{\mathrm{in}}$, to be 0.1. When the metallicity decreases because of e.g. infall of low-metallicity gas we manually dilute the dust in proportion to the decrease in the metallicity (further details are given in Section 2.3.2). We assume that the newly-formed dust mass arising from stellar evolution has a log-normal grain mass distribution with an average grain radius of 0.1 $\micron$ and a standard deviation of 0.47 (in the log). This assumes that the dust grains produced by stars have submicron sizes (see Hirashita & Aoyama 2019 and references therein for a discussion of arguments in support of this assumption). For dust destruction by SN shocks, we adopt a grain-size-dependent destruction efficiency, which is multiplied by the rate at which the interstellar dust is swept up by SN shocks.

The remaining three mechanisms only occur in specific phases of the ISM. Thus, we consider two ISM phases defined by the hydrogen number density $n_{\rm{H}}\,(\rm{cm^{-3}})$ and the gas temperature $T_{\rm{gas}}\,(\rm{K})$. One is the diffuse/warm ISM with $(n_{\rm{H}},\,T_{\rm{gas}})\sim(10^{-1},\,10^{4})$, and the other is the dense/cold ISM with $(n_{\rm{H}},\,T_{\rm{gas}})\sim(10^{2},\,10^{2})$. In principle, the mass of gas in each phase could be determined from the density and temperature distributions of all the gas elements associated with each subhalo in the original simulation. However, we instead treat the fraction of gas in each phase as a parameter of our model and adopt the following fixed values for the density and temperature in each phase: $(n_{\rm{H}},\,T_{\rm{gas}})=(0.3,\,10^{4})$ for the diffuse ISM and $(n_{\rm{H}},\,T_{\rm{gas}})=(300,\,100)$ for the dense ISM.

We apply this treatment for the following reasons. First, dense gas is not directly resolved in the simulation and instead relies on a sub-resolution two-phase model (Springel & Hernquist, 2003). Thus, coagulation and accretion need to be treated with sub-grid prescriptions (e.g. Aoyama et al., 2018). Second, the density and temperature are degenerate with the assumed turbulent velocities in coagulation and shattering, and uncertainty in the turbulent velocities can be absorbed in the coagulation and shattering efficiencies that are adjustable in this work. We will discuss more details in the next section. For these reasons, we simply adopt a fixed characteristic density and temperature for each ISM phase.

We denote the mass fraction of the dense ISM as $\eta_{\rm{dense}}$ and thus the mass fraction of the diffuse ISM is $1-\eta_{\rm{dense}}$. These fractions are included in the calculation in the form of an effective time-step. Let $\Delta t$ be the time-step in our calculations. We distribute this time-step to the dense and diffuse ISM by adopting $\eta_{\mathrm{dense}}\Delta t$ for the processes occurring in the dense ISM (accretion and coagulation) and $(1-\eta_{\mathrm{dense}})\Delta t$ for the process in the diffuse ISM (shattering). The evolution of the grain size distribution by accretion is solved by considering the grain growth rate dependent on the metallicity and the grain radius. The change of the grain size distribution by coagulation and shattering are treated by the Smoluchowski-type equations with kernel functions determined consistently with the grain–grain collision rates in a turbulent media. The turbulent velocity in each ISM phase is given by the model of Ormel et al. (2009) albeit with an adjusted velocity normalization.

We discretize the entire grain radius range ($a=0.3$ nm–10 $\micron$) into 64 logarithmic bins. We apply $n(a,\,t)=0$ at the upper and lower radii as boundary conditions. That is, we remove grains that evolve beyond the above radius range. The interval between TNG snapshots is often longer than the characteristic time-scale of interstellar dust processing, in which case we divide the time interval into many sub-steps, each of which is shorter than the timescale of the most rapid process.

Shattering and coagulation are both associated with grain-grain collision in the ISM. As mentioned above, the collision rate is regulated by the kernel function, which is a product of grain cross-section and grain velocity, in the Smoluchowski-type equations. HM20 simply used the geometrical cross-sections of compact spheres. The cross-sections depend on the grain shape. As mentioned above, the grain velocity also depends on the adopted turbulence model. Since the balance between coagulation and shattering determines the shape of grain size distribution in evolved (metal-enriched) galaxies, it is worth adjusting the coagulation and shattering efficiencies to clarify the robustness and uncertainty in the model. To regulate the shattering and coagulation efficiencies, we introduce dimensionless parameters $\omega_{\mathrm{shat}}$ and $\omega_{\mathrm{coag}}$ and multiply the respective kernel functions. These parameters are useful to modify the relative strength between coagulation and shattering. Setting $\omega_{\mathrm{shat}}=\omega_{\mathrm{coag}}=1$ gives the same rates as in HM20.

Determination of the dense gas fraction $\eta_{\mathrm{dense}}$ is critical since some processes occur only in one of the two ISM phases. Given the resolution limitations of cosmological simulations such as TNG, we do not measure this directly from the hydrodynamical outputs, but instead assume a relation between $\eta_{\mathrm{dense}}$ and the SFR:

where $M_{\mathrm{gas}}^{*}$ is the interstellar gas mass, $\epsilon_{*}$ is the star-formation efficiency, and $\tau_{\mathrm{ff}}=2.51(n_{\mathrm{H}}/300~{}\mathrm{cm}^{-3})^{-1/2}$ Myr is the free-fall time. We obtain $M_{\mathrm{gas}}^{*}$ and SFR for each subhalo from TNG300-1, and fix $\epsilon_{*}=0.01$ (Krumholz & McKee, 2005). For consistency with the assumed density of the dense ISM, we fix $n_{\mathrm{H}}=300$ cm^-3 in the estimation of the free-fall time.

As shown later, the resulting grain size distributions are sensitive to the choice of $\eta_{\mathrm{dense}}$. Thus, we also examine some models with fixed values of $\eta_{\mathrm{dense}}$ to understand the response of the results to the change of $\eta_{\mathrm{dense}}$ (Section 3.1).

Galaxies in TNG grow through hierarchical merging as well as via the smooth accretion of baryons and dark matter from the cosmic web. HM20 considered the dust evolution in an isolated field galaxy and did not account for either mode of growth. Thus, we need to extend HM20’s model to deal with subhalo growth through smooth accretion (we hereafter refer to this process as gas accretion to avoid confusion with dust growth by accretion) and merging.

To model gas accretion, we note that galaxies accrete gas from the circumgalactic medium (CGM). We assume that the CGM contributes no dust and only few metals. Under this assumption, gas accretion from the CGM simply dilutes the dust-to-gas ratio. We therefore regard decreasing metallicity as an indicator of significant gas accretion. Furthermore, we assume that neither the dust-to-metal ratio nor the functional shape of the grain size distribution changes during a gas accretion episode. Thus, if the gas accretion changes the metallicity from $Z(t)$ to $Z(t+\Delta t)[<Z(t)]$, we update the grain size distribution $n(a,\,t)$ according to the ratio between $Z(t)$ and $Z(t+\Delta t)$. That is, the grain size distribution at $t+\Delta t$ is evaluated as

Note that we only evolve the dust properties in this way when no mergers occur in a given time-step.

Although TNG provides information of progenitors and descendants for each subhalo, we do not know exactly when a merger happens because TNG stores snapshot data with a time interval of roughly 0.1 Gyr. To overcome this limitation we assume that (i) the merger happens at the end of the interval between two snapshots, (ii) the amount of accreting gas mass right before the merger is proportional to the gas mass in the progenitor, and (iii) the total gas mass is conserved during the merging. In practice, we insert a virtual node right before the descendant subhalo, evolve grain size distributions from the progenitor to this virtual node, and combine the grain size distributions from two virtual nodes.

Our method is illustrated in Fig. 2. $S_{n}$ denotes the $n$th snapshot, with $n$ increasing with cosmic time. Blue circles are the original subhalos from TNG and yellow circles show the virtual nodes. The unprimed quantities, $M_{\mathrm{gas}}$ and $n(a)$, stand for the total gas mass directly obtained from TNG and the grain size distribution calculated by our model, respectively. We use the primed quantities $M_{\mathrm{gas}}^{\prime}$ and $n^{\prime}(a)$ for virtual nodes (pre-merging stage), and estimate them as follows. The progenitors (blue circles) in $S_{n}$ eventually evolve into the pre-merging stage (yellow circle) in $S_{n+1}^{\prime}$. Under assumption (ii) we have $M_{\mathrm{gas1}}^{\prime}/M_{\mathrm{gas2}}^{\prime}=M_{\mathrm{gas1}}/M_{\mathrm{gas2}}$. Additionally, under assumption (iii), we have $M_{\mathrm{gas1}}^{\prime}+M_{\mathrm{gas2}}^{\prime}=M_{\mathrm{gas}}$. Hence, the grain size distribution in the merged system $n(a)$ is estimated by the following average:

The evolution from $n_{i}(a)$ to $n_{i}^{\prime}(a)$ is calculated using our dust evolution model, as described in Section 2.2. In this way the information available at snapshot $S_{n}$ allows us to obtain the grain size distribution at $S_{n+1}$.

As explained in Section 2.1, we consider two definitions of gas mass associated with the dust evolution: one is the the total gas mass within the subhalo ($M_{\mathrm{gas,tot}}$), and the other is the gas mass within twice the stellar half mass radius ($M_{\mathrm{gas,2R*}}$). We compare these two alternatives to investigate the uncertainty caused by the definition of gas mass associated with star formation and dust processing in galaxies.

In Fig. 3 we show the grain size distributions obtained using $M_{\mathrm{gas,tot}}$ and $M_{\mathrm{gas,2R*}}$ in Panels (a) and (b), respectively. The case with $M_{\mathrm{gas,tot}}$ tends to have more intermediate-sized ($a\sim 0.03~{}\micron$) grains than the case with $M_{\mathrm{gas,2R*}}$ at $z\lesssim 1$. Although the grain size distribution has a large dispersion for small grains at early epochs, the grain size distributions are not significantly affected by the choice of gas mass for large grains or at later epochs. The sensitive behaviour at early epochs is due to dust growth by accretion, which causes a drastic increase of small grains on a short time-scale. However, the overall median behaviour is not strongly affected by the choice of gas mass, which motivates our fiducial choice of $M_{\mathrm{gas,2R*}}$ for our modeling below.

The dense gas fraction $\eta_{\mathrm{dense}}$ affects dust evolution processes, especially accretion and coagulation, which are more significant in the dense ISM. To demonstrate the impact of dense gas fraction, we first fix the dense gas fraction at different values ($\eta_{\mathrm{dense}}=0.1$, 0.3, and 0.5) and show the results in Fig. 4. It is clear that different dense gas fractions not only change the dust abundance but also alter the grain size distribution. A larger value of $\eta_{\mathrm{dense}}$ leads to an enhanced small-grain abundance in the early phases because accretion is more efficient. At later epochs, we observe a higher abundance of large grains for higher $\eta_{\mathrm{dense}}$ because coagulation is more efficient. With $\eta_{\mathrm{dense}}=0.5$, the overall slope of grain size distribution approaches the MRN distribution, which is consistent with the results of one-zone models in HM20.

The shape of the grain size distribution is affected by the efficiency of shattering as well as coagulation, especially at later stages. As introduced in Section 2.2.2, we adjust the strength of coagulation ($\omega_{\mathrm{coag}}$) and shattering ($\omega_{\mathrm{shat}}$) to investigate if we obtain grain size distributions more similar to the MRN distribution at $z=0$. We show the results in Figs. 5 and 6 for coagulation and shattering, respectively. We find that the enhancement of the coagulation efficiency increases the large grain abundance and tends to produce a MRN-like slope at $z\lesssim 1$. On the other hand, although a reduction of shattering strength increases the large-grain abundance, this is still not enough to reproduce the MRN slope (Fig. 6b). Efficient coagulation is essential to produce a slope similar to the MRN distribution as noted in Hirashita & Aoyama (2019). Note that neither $\omega_{\mathrm{coag}}$ nor $\omega_{\mathrm{shat}}$ affects the grain size distributions at $z\gtrsim 3$ because these processes only play a minor role at that early epoch.

In light of the results in Figs. 5 and 6, the enhancement of coagulation might be the more appropriate way to reproduce the MRN slope. Although our purpose is not to reproduce the MRN grain size distribution exactly, it is still worth exploring which solutions yield grain size distributions similar to that of the Milky Way. Considering that grain cross-sections and grain velocities are uncertain to an order of magnitude, we use $\omega_{\mathrm{coag}}=10$ as the fiducial model for the full sample. To minimize fine-tuning we fix $\omega_{\mathrm{shat}}=1$. Despite this tuning, we still expect that the physical interpretations given below will hold for any simulations in the future, at least qualitatively. Since there have been no calculations of full grain size distributions in cosmological simulations, our fiducial model serves as a robust basis for investigating which process dominate the evolution of the grain size distribution in the MW-like galaxies at various epochs.

Based on the above tests, we fix parameters for the full-sample calculations described in the next subsection. As shown in Section 3.1.1, there are two definitions of galactic gas mass available in TNG, but the choice between these does not significantly affect the grain size distributions predicted by our model. Since most of the dust formation and destruction processes occur around stars, rather than in the dark matter halo, we use the gas mass within twice the stellar half-mass radius in our dust model. Hereafter, we simply denote $M_{\mathrm{gas,2R}}$ as $M_{\mathrm{gas}}$.

We assumed various constant values for the dense-gas fraction ($\eta_{\mathrm{dense}}$) in Section 3.1.2 for the purpose of testing the effect of this parameter. However, we do not fix $\eta_{\mathrm{dense}}$ in our fiducial model, but instead use the recipe described in Section 2.3.1 (equation 2) to determine its value in each subhalo individually. As discussed in Section 3.1.3, we adopt $\omega_{\mathrm{coag}}=10$, but leave the shattering efficiency unchanged ($\omega_{\mathrm{shat}}=1$).

The relation between dust-to-gas ratio and metallicity has often been used to test dust enrichment models (e.g. Lisenfeld & Ferrara, 1998; Dwek, 1998). To examine if our model reproduces the dust enrichment reasonably, Fig. 7 shows the evolution of dust-to-gas ratio for all galaxies in our sample, together with their progenitors at earlier times, over the redshift range $z=0$–3. We compare the calculated $D_{\mathrm{tot}}$–$Z$ relation to observational results for local galaxies (Rémy-Ruyer et al., 2014), for which the dust mass is estimated from IR-to-submm photometry, and the gas mass is estimated from H i and CO data with a metallicity-dependent CO-to-H₂ conversion factor.

The dust-to-gas ratios our model predicts at early times are slightly below the linear relation expected from the dust condensation efficiency in stellar ejecta ($=0.1$; dotted line in Fig. 7). Dust destruction by SNe is responsible for the shift towards lower dust-to-gas ratios. The dust-to-gas ratio and metallicity of galaxies increase with decreasing redshift, but the relation stays on a well defined narrow track, except for a stripe of low dust-to-gas ratios around $Z\sim 0.2$ Z_☉. This relation between dust-to-gas ratio and metallicity is broadly consistent with observations, including the nonlinear (steep) increase of dust-to-gas ratio around $Z\sim 0.3$ Z_☉ (see also Aniano et al., 2020, for a recent analysis). The fraction of points in the band toward low $D_{\mathrm{tot}}$ around $Z\sim 0.2$ Z_☉ is only about 1 percent of all the points. We will discuss the reason for these outliers in Section 4.

To illustrate the redshift evolution of the $D_{\mathrm{tot}}$–$Z$ relation more clearly, we plot the $D_{\mathrm{tot}}$–$Z$ diagram at $z=0$, 1, 2, and 3 separately in Fig. 8 (colored circles) overlaid on the combined distribution across all redshifts (gray background cloud). At $z=3$, most galaxies are near the line predicted from the stellar yield ($Z\sim f_{\mathrm{in}}D_{\mathrm{tot}}$), but some systems at $Z\gtrsim 0.3$ Z_☉ have started to increase their dust abundance by accretion. At $z=2$ the distribution moves upwards and to the right, and most galaxies are located in the regime where the increase in dust-to-gas ratio becomes non-linear. During this epoch, dust production is no longer dominated by stellar sources but instead by dust growth through accretion.

After this epoch, the distribution gradually approaches the upper-right corner of the diagram, near the saturation limit of $D_{\mathrm{tot}}=Z$ at $z=1$. It is clear from Fig. 8c that the outliers start to appear at this epoch. Progressively more outliers appear at lower redshift, such that at $z=0$ the fraction reaches roughly 4 percent. Meanwhile, the tight correlation gradually broadens with larger scatter at $z=0$. As we will show later, the average metallicity of our sample reaches its peak at $z\sim 0.5$ and then drops (Fig. 11a). This decrease of metallicity affects the dust-to-gas ratio, because we disable accretion when the gas-phase metallicity is decreasing (although SN destruction still takes place). The large dispersion of $\eta_{\mathrm{dense}}$ at $z\sim 0$ (as shown later in Fig. 11b) could also increase the scatter in Fig. 8 because the dense gas fraction also affects the accretion and the relative strength of coagulation and shattering, which will have a mild effect on the total dust abundance.

In summary, our model is largely consistent with the observed relation between dust-to-gas ratio and metallicity. This indicates that we successfully reproduce the dust enrichment history of the MW-like galaxy sample. We further predict that the relation between dust-to-gas ratio and metallicity does not depend strongly on redshift. There is an increasing trend of dust-to-gas ratio with decreasing redshift, but the evolutionary track on the $D_{\mathrm{tot}}$–$Z$ is almost invariant. This is consistent with other theoretical models of dust enrichment in cosmological simulations by Hou et al. (2019) (at $z\lesssim 2$) and Li et al. (2019), and the semi-analytic model of Popping et al. (2017) and Triani et al. (2020) (although Vijayan et al. 2019 showed an additional dependence on the stellar age). The tight $D_{\mathrm{tot}}$–$Z$ relation is also compatible with the analytical result that metallicity is the main driver of dust enrichment (Inoue, 2011).

We present the grain size distributions at different redshifts ($z=0$–$3$) for the progenitors of our MW-like galaxy sample in Fig. 9. At $z>3$, the size distributions are dominated by stellar dust production, and their shapes are roughly described by a log-normal function centred at $a\sim 0.1~{}\micron$, as assumed for the grain size distribution of stellar dust (Section 2.2.1). To quantify the distribution we show the median with 25th and 75th percentiles at each grain radius at each redshift. We first discuss the evolving median behavior and then focus on the scatter in the sample.

At $z=3$, the median distribution has two peaks, and the larger-grain maximum ($a\sim 0.2~{}\micron$)²²2The centre of the lognormal distribution is at $a=0.1~{}\micron$, but since we multiply the distribution by $a^{4}$, the peak of the distribution in the figure is located around $a\sim 0.2~{}\micron$. is higher than the peak in the small-size regime. The larger-grain peak is created by stellar dust production (Section 2.2.1), while the smaller-grain peak arises from dust growth by accretion, which is more efficient for small grains (Kuo & Hirashita, 2012; Hirashita & Aoyama, 2019). The grain size distribution increases from $z=3$ to $z=2$ at almost all grain radii, especially for small sizes, reaching a nearly flat distribution at $a\lesssim 0.1~{}\micron$.

From $z=2$ to $z=1$ a further increase in dust abundance makes the integrated size distribution the largest among the four epochs. At $z=1$, because of efficient coagulation, the grain size distribution tends to a smooth power-law-like shape, approaching an MRN slope ($n\propto a^{-3.5}$). From $z=1$ to $z=0$ grains with $a\gtrsim 0.1~{}\micron$ are removed by shattering while the small-size end of the distribution remains nearly constant. As a result the grain size distributions we predict are more dominated by small grains than in the MRN distribution. Thus, shattering plays an important role in the grain size distribution at $z\sim 0$, as we discuss further below.

The scatter of the grain size distribution, shown by 25th–75th percentiles in Fig. 9, can arise from the different ISM properties in each galaxy. At $z=3$, the scatter in the grain size distribution is small at large radii, where the grain abundance is dominated by stellar dust production. The dispersion at large grain sizes is roughly the same as in the metallicity at $z=3$. On the other hand, the dispersion at small grain radii is large. Since the peak in the grain size distribution at small radii is caused by accretion (dust growth), the large scatter reflects variation in the fraction of dense gas in which that process occurs.

At $z=2$, the scatter becomes larger but then decreases by $z=1$. Thus, the typical epoch at which the small-grain abundance increases is around $z\sim 2$, and this is predominantly driven by accretion, which also induces a scatter at large grain radii through coagulation. The grain size distributions converge to a smooth MRN-like shape at $z\sim 1$, with extremely small scatter. This indicates the universality of the grain size distribution achieved as a result of efficient interstellar processing. The scatter tends to be larger at large grain radii if coagulation is weaker (Section 3.1). This means that coagulation plays an important role in realizing this ‘universality’ of the grain size distribution. At $z=0$, the scatter becomes larger again, especially at large grain radii. This is driven by the decrease of large grains by shattering and inefficient coagulation, due to a decrease in the dense gas fractions of the simulated galaxies, as discussed further below.

Based on the grain size distributions for the full sample above, we calculate extinction curves. Fig. 10 shows the median extinction curves at different epochs, together with their variation. Overall, the features in extinction curve shape can be interpreted in the context of the corresponding grain size distributions.

At $z=3$, the extinction curves are flat because of the lack of small grains. When the abundance of small grains grows at later times, the extinction curves become steeper at short wavelengths. At $z=1$, because the grain size distributions are similar to MRN, the resulting extinction curves are similar to the MW extinction curve. The 2175 Å bump is smaller in our calculation than in the observed MW curve. At $z=0$, because the large-grain abundance is reduced by shattering, our model predicts extinction curves that are steeper than the MW extinction curve. A typical Milky Way galaxy in our model has a $z=0$ extinction curve whose steepness is more similar to that of the Small Magellanic Cloud (SMC) extinction curve, although our predictions show a 2175 Å carbon bump, which is not present in the SMC. The scatter of predicted extinction curves is large at $z=0$, and marginally covers the MW curve in the near UV and optical.

It is clear that the assembly history of a galaxy can have a significant effect on its present-day dust properties and observables, including the extinction curve. Although we cannot reject the possibility that the MW extinction curve may not be representative for our sample, selected on stellar mass and $z=0$ star formation rate, the difference between the average prediction of our sample and the MW extinction curve at $z=0$ is worth further consideration.