The impact of binary stars on the dust and metal evolution of galaxies

We present detailed implementations of (a) binary stellar evolution (using binary_c) and (b) dust production and destruction into the cosmological semi-analytic galaxy evolution simulation, L-Galaxies. This new version of L-Galaxies is compared to a version assuming only single stars and to global and spatially-resolved observational data across a range of redshifts ($z$). We find that binaries have a negligible impact on the stellar masses, gas masses, and star formation rates of galaxies only if the total mass ejected by massive stars is unchanged. This is because massive stars determine the strength of supernova (SN) feedback, which in turn regulates galaxy growth. Binary effects, such as common envelope ejection and novae, affect carbon and nitrogen enrichment in galaxies, however heavier alpha elements are more affected by the choice of SN and wind yields. Unlike many other simulations, the new L-Galaxies reproduces observed dust-to-metal (DTM) and dust-to-gas (DTG) ratios at $z\sim{}0-4$. This is mainly due to shorter dust accretion timescales in dust-rich environments. However, dust masses are under-predicted at $z>4$, highlighting the need for enhanced dust production at early times in simulations, possibly accompanied by increased star formation. On sub-galactic scales, there is very good agreement between L-Galaxies and observed dust and metal radial profiles at $z=0$. A drop in DTM ratio is also found in diffuse, low-metallicity regions, contradicting the assumption of a universal value. We hope that this work serves as a useful template for binary stellar evolution implementations in other cosmological simulations in future.


INTRODUCTION
Stellar evolution plays a critical role in overall galaxy evolution.In particular, it determines the timing, energetics, and chemical composition of the material ejected by stellar winds and supernovae (SNe).This material drives the stellar feedback and chemical enrichment of galaxies, and consequently many other key evolutionary processes, such as gas cooling and star formation.It is therefore crucial for galaxy evolution simulations to take proper account of stellar evolution in their 'sub-grid modelling'.
However, a major aspect of stellar evolution that is rarely considered in cosmological-scale simulations is stellar multiplicity.Observations now suggest that over 50 per cent of massive stars (i.e.those above ∼ 8 M ⊙ ) evolve in binary or higher-multiple systems (e.g.Sana et al. 2012).This binary stellar evolution (BSE) can have a significant impact on the lifetime, nucleosynthesis, and end state of these stars, due to the complex interactions with their lower-mass counterparts (e.g.Podsiadlowski et al. 1992;De Marco & Izzard 2017;Moe & Di Stefano 2017).Given that these massive stars are also the major ★ E-mail: r.yates3@herts.ac.uk contributors to ionising photons, stellar feedback, and chemical enrichment in galaxies (e.g.Nomoto et al. 2013;Götberg et al. 2019;Kobayashi et al. 2020) it is important to assess how their binarity impacts galaxy evolution as a whole.
Although BSE has been included in some stand-alone galaxy chemical enrichment (GCE) models (e.g.De Donder & Vanbeveren 2004;Mennekens & Vanbeveren 2016;Côté et al. 2018), and individual aspects of massive binaries have been studied on larger scales, such as transient event rates, compact binary mergers, and r-process element production e.g.Van de Voort et al. 2015;Schneider et al. 2017;Stanway et al. 2018;Chruslinska et al. 2018;Côté et al. 2019;Artale et al. 2020;Kobayashi et al. 2023;Peron et al. 2023), most cosmological-scale galaxy evolution simulations to date only consider single-star stellar populations in their sub-grid modelling (e.g.Santa Cruz, Arrigoni et al. 2010; L-Galaxies, Yates et al. 2013;Illustris, Vogelsberger et al. 2013;Gaea, De Lucia et al. 2014;Eagle, Schaye et al. 2015;Fire-2, Hopkins et al. 2018;Sag, Collacchioni et al. 2018;Illustris-TNG, Pillepich et al. 2018;Simba, Davé et al. 2019;Simba-C, Hough et al. 2023).Consequently, important binary phenomena such as type Ia supernovae (SNe-Ia) are also only treated approximately in these simulations.Typically, an analytic delay-time distribution (DTD) is assumed and the rate of SN-Ia events per stellar population scaled to match the iron abundances observed in early-type galaxies (ETGs) or the solar neighbourhood (e.g.Greggio 2005;Matteucci et al. 2006Matteucci et al. , 2009;;Yates et al. 2013).While this approach returns a good match to galaxy-scale observations, it does so essentially by construction.Other binary phenomena, such as common envelope (CE) ejection and novae, are often ignored altogether.
Similarly, dust can play a critical role in galaxy evolution.Observationally, intervening dust attenuates the UV and optical light emitted by galaxies, re-emitting it at longer wavelengths and therefore changing their observed colours or obscuring them completely (for a recent review, see Galliano et al. 2018).It can also provide a useful independent tracer of the metal and gas content in galaxies (e.g.Eales et al. 2012;Brinchmann et al. 2013).Physically, dust is a source of depletion for refractory elements, thus changing the gas-phase abundances in the interstellar medium (ISM) (e.g.Savage & Sembach 1996;Jenkins 2009) and affecting cooling rates (e.g.Sutherland & Dopita 1993).
Therefore, in this work we implement new models for (a) binary stellar populations and (b) dust into the semi-analytic galaxy evolution simulation, L-Galaxies 2020 (Henriques et al. 2020;Yates et al. 2021a), in order to self-consistently constrain the impact of BSE on the dust and metal evolution of galaxies.In particular, we present global and radially-resolved measurements of gas and stellar masses, chemical abundances, DTM and DTG ratios, SN rates, and cosmic densities, and compare these to the latest multi-redshift observations.
In Section 2, we describe the L-Galaxies base simulation used in this work.Section 3 presents the binary_c BSE model and its implementation into L-Galaxies, while Section 4 presents our new dust production and destruction model.In Section 5, we provide results from two versions of L-Galaxies, one with BSE included and one without, and make detailed comparisons to observational data.Section 6 provides our conclusions and future prospects.Throughout, we take a Hubble parameter of ℎ ≡  0 /100 = 0.673, log to the base 10, and distances in co-moving units.
Semi-analytic simulations like L-Galaxies treat all baryonic processes analytically, in a similar way to sub-grid models in hydrodynamical simulations.These models allow the transfer of mass and energy between 'homogeneous spatial components' inside dark matter haloes and galaxies, rather than between computational particles or cells.This method is highly computationally efficient, allowing many millions of galaxies to be modelled from high  to the present day at a computational cost (in terms of CPU time) many thousands of times smaller than for hydrodynamical simulations of comparable size.L-Galaxies is able to estimate global properties such as stellar masses and specific star formation rates (sSFR = SFR/ * ) to within ∼ 0.3 dex of those predicted by fully hydrodynamical simulations run on the same dark matter haloes, for systems where active galactic nucleus (AGN) feedback is not significant (Ayromlou et al. 2021).
There is a wide range of baryonic physics included in L-Galaxies, as described in Henriques et al. (2020) and the L-Galaxies 2020 Model Description2 .Most recently, models for (a) the partition of HI and H 2 within the ISM (Fu et al. 2010), (b) delayed chemical enrichment from SNe and stellar winds (Yates et al. 2013), and (c) radially-resolved ISM and stellar discs (Fu et al. 2013) were incorporated.Resolved discs are represented by twelve 'radial rings' with outer radii given by log(  ) = log(0.44)+  log(1.5)kpc/ℎ, where  = 0 − 11 (Fu et al. 2013).Gas (and dust) is then allowed to flow inwards over time with a velocity given by  inflow =   , where   = 0.6 km s −1 kpc −1 .L-Galaxies can therefore resolve the gas and stellar components of galaxies down to sub-kpc scales in their inner ∼ 1 kpc, as well as at lower spatial resolution in their outskirts and surrounding subhalo.The following spatial components are considered: central supermassive black hole (SMBH), stellar bulge, radially-resolved stellar disc, radially-resolved cold-gas disc (representing the interstellar medium, ISM), hot-gas reservoir (representing the circumgalactic medium, CGM), stellar halo, and an ejecta reservoir (representing gas ejected out of a galaxy's dark matter subhalo via feedback).
Star formation histories (SFHs) are also stored for each stellar component (Yates et al. 2013;Shamshiri et al. 2015).These comprise 20 bins, structured such that there is higher temporal resolution in the recent past, enabling properties such as rest-frame magnitudes and colours to be optimally calculated on-the-fly within the simulation.These SFHs also allow the timing of chemical element ejection by stars and SNe to be accurately calculated (see Section 3).
Finally, following Yates et al. (2021a), we base our current study on a version of L-Galaxies known as the 'modified model' (hereafter, MM).This version allows the efficient ejection of newly-synthesised metals out of galaxies.90 per cent of the material ejected by SNe-II and 70 per cent of the material ejected by SNe-Ia is directly deposited into the CGM surrounding galaxies, without first mixing with the ISM.This material can then cool back onto the disc at later times.The remaining ejecta material is allowed to fully mix with the local ISM, before being available for star formation or later removal through entrainment in SN-driven winds.This version of L-Galaxies was designed to emulate the metal-rich galactic outflows observed in many real (e.g.Martin et al. 2002;Strickland et al. 2004;Tumlinson et al. 2011) and simulated (e.g.Emerick et al. 2020;Gutcke et al. 2021) star-forming galaxies, and produces a gradual enrichment of galaxies over cosmic time which is more in line with recent obser-vational inferences from ALMA, Keck, and JWST (e.g.Jones et al. 2020;Sanders et al. 2021;Schaerer et al. 2022;Arellano-Córdova et al. 2022;Trump et al. 2023;Sanders et al. 2023;Curti et al. 2023;Brinchmann 2023).

BINARY STELLAR EVOLUTION (BSE) MODEL
In this work, we implement the mass-and metallicity-dependent chemical yields, stellar lifetimes, and event rates from synthetic stellar populations generated by the binary_c3 code (Izzard et al. 2006(Izzard et al. , 2009(Izzard et al. , 2018;;Izzard & Jermyn 2022) into the L-Galaxies simulation described above.Binary_c is a BSE framework designed to model the evolution of single and binary stars using semi-analytic methods.Included in this modelling is the calculation of chemical yields from a wide range of ejection processes, including stellar and Wolf-Rayet (WR) winds, Roche-lobe overflow (RLOF), CE ejection, various types of (super)nova explosion, and others.We make use of a recent extension to this framework, binary_c-python4 (Hendriks & Izzard 2023), to calculate yields and event rates on a "per population basis" as a function of time after their formation.Such yields are perfect for implementation into large-scale galaxy evolution simulations, which do not have the mass or time resolution required to resolve individual stars.binary_c is preferred in this work over other BSE models in the literature (e.g.SeBa, Portegies Zwart & Verbunt 1996;BPASS, Eldridge et al. 2017;and MOBSE, Giacobbo et al. 2018) because it is specifically designed to rapidly model the impact of (binary) stellar evolution on a very wide range of nucleosynthetic yields on-the-fly.
Below, we describe the key aspects of binary_c of relevance here.A list of the key binary_c parameters used to generate our stellar populations is also available online in the Supplementary Material.5

Binary stellar populations
In this work, we utilise stellar populations generated by binary_c v2.2.1 using binary_c-python v0.9.4 for six distinct metallicities ( ≡  Z / b = 0.0001, 0.001, 0.004, 0.008, 0.01, 0.03).The standard binary_c set-up is used for this work, which includes a Kroupa (2001) initial mass function (IMF) and a maximum star mass of 80 M ⊙ .The fraction of stars that are in binary systems, as well as the initial orbital period and mass-ratio distribution, are determined following Moe & Di Stefano (2017).This mass-dependent prescription predicts that ∼ 75 per cent of 8 M ⊙ stars are in higher-order (binary, triple, or quadruple) systems, with this fraction reaching ≳ 90 per cent at ≳ 25 M ⊙ (see fig. 39 from Moe & Di Stefano 2017).All multiple systems are treated as binaries here, leaving the higherorder effects caused by triple and quadruple interactions for future work.We also note that in the case of wide binaries, the component stars essentially evolve as single stars even though they are strictly in binary systems.Below, we outline the key prescriptions included in binary_c v2.2.1 which are of most relevance here.The impact of modifying these standard assumptions will be the focus of future work.

Stellar yields
binary_c provides yields for a total of 117 chemical elements (487 isotopes), opening-up a whole new window of investigation into the chemical evolution of galaxies in cosmological-scale simulations.The implementation of binary_c into L-Galaxies has been designed to allow the user to freely choose which elements/isotopes are tracked.For the purposes of this work, we only consider the 11 chemical elements commonly used in current cosmological-scale simulations, including previous versions of L-Galaxies: H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe.Follow-up papers in this series will consider additional elements.These 11 elements account for ∼ 98% of the total baryonic mass in the solar photosphere, so also allow for a robust study of the total metallicity in most astrophysical regions.
For AGB winds from both single and binary stars, the stellar yields calculated by Izzard et al. (2004Izzard et al. ( , 2006, based on those from Karakas et al. 2002) are used.For SNe-Ia, we assume the yields from the DD2 model (i.e.assuming a delayed detonation at a density of 2.2×10 7 g cm 3 ) from Iwamoto et al. (1999) for Chandrasekhar-mass detonations and from Livne & Arnett (1995) for sub-Chandrasekharmass detonations.For SNe-II, the metallicity-dependent yields of Chieffi & Limongi (2004) are used.For details on the yield sets for less-common transient events, along with information on the full range of ejection processes modelled in binary_c, see the binary_c online documentation.
Fig. 1 provides an illustration of the normalised metal ejection rate (i.e. the metal ejection rate per solar mass of stars formed) for a stellar population with birth metallicity  = 0.004, as a function of time after star formation.This is obtained by summing all nine chemical elements heavier than He considered in this work.The full complement of ejection processes modelled in binary_c is shown (see legend).Stellar winds dominate enrichment at very early times ( ≲ 6 Myr), after which prompt ejection events such as WR winds, SNe-Ibc, and SNe-II dominate until  ∼ 100 Myr.At later times, AGB winds and finally novae dominate.Enrichment from stellar winds is a strong function of metallicity in binary_c.For example, the normalised metal ejection rate from winds at  = 1 Myr increases by ∼ 4.5 dex from  = 0.0001 to 0.03.This is because metals increase the opacity of the surface layers, thus increasing the radiation pressure on them.
Yield sets such as these represent a significant improvement for GCE models in simulations, which typically just assume single-star yields for SNe-II and AGB winds, along with a simplified prescription for SN-Ia enrichment.

Common envelope (CE) ejection
During their evolution, stars in binary systems expand and may overflow their Roche lobe, i.e. the volume in which material is strictly bound to only that star.This leads to mass transfer which is either stable or unstable.We determine the stability of mass transfer based on the criteria from Claeys et al. (2014).When mass transfer is unstable, the system undergoes a CE evolution where the transferred material engulfs both stars, and orbital energy and angular momentum is lost through dynamical friction with this envelope.This leads to shrinkage of the orbit, and possibly the merging of the two stars.To determine the final orbital separation we use the  CE −  CE prescription of Webbink (1984) and Hurley et al. (2002) to calculate the energy necessary to expel the envelope, where  CE is the efficiency with which orbital energy is transferred to the CE and  CE scales the binding energy of the donor's envelope.We take  CE = 1, i.e. 100 per cent of the transferred orbital energy is transformed into energy to eject the envelope, and  CE = 0.5.

Remnants of massive stars
The binary_c code calculates the remnant mass of a massive star following the prescription of Hurley et al. (2000Hurley et al. ( , 2002)), which depends on the CO core mass at the end of the star's life,  rem = 1.17 + 0.09 CO . (1) Stars which undergo SNe-II typically have  CO ≲ 35 M ⊙ (see e.g.Limongi & Chieffi 2018), thus obtaining remnant masses of  rem ≲ 4.3 M ⊙ , meaning that the majority of the star's initial mass is ejected during the SN or in the pre-SN winds.We note that other prescriptions which account for dependencies on metallicity and explosion energy (e.g.Fryer et al. 2012;Ertl et al. 2016) tend to predict higher remnant masses, implying lower ejected yields of newlysynthesised material from SNe-II.However, the detailed effects of binary evolution on pre-SN mass loss could counteract these effects to some extent (e.g.Schneider et al. 2021;Farmer et al. 2023).For example, Schneider et al. (2021) find that the stripping of a massive star by its binary companion can lead to larger carbon abundances and/or lower  CO at core-helium exhaustion, and thus higher explosion energies, lower remnant masses, and greater ejection of heavy elements at core collapse.A detailed investigation into the impact of different remnant and fallback prescriptions on the chemical evolution of galaxies will be the focus of future work.

SN-Ia progenitors
The evolution of SN-Ia progenitor systems is modelled in bi-nary_c following Claeys et al. (2014).In this formalism, three SN-Ia pathways are considered; the single-degenerate scenario (SNIa-CHAND), double-degenerate scenario (SNIa-CHAND-Coal), and sub-Chandrasekhar-mass scenario (SNIa-ELD).The relative frequency of the latter two scenarios in binary_c depends on the assumed mass that must be accreted before a white dwarf can ignite as an edge-lit detonation (ELD).Here, this accretion mass is set to 0.15 M ⊙ .The double-degenerate scenario is the dominant pathway for latetime (i.e.≳ 200 Myr) SNe-Ia in the Claeys et al. (2014) formalism, whereas the single-degenerate scenario dominates at earlier times (see their fig.7).For the canonical CE efficiency of  ce = 1, a SN-Ia DTD that roughly traces the power-law shape inferred from observations of SN-Ia rates in local galaxies is obtained (e.g.Maoz et al. 2011Maoz et al. , 2014)).However, the normalisation of this DTD in binary_c is significantly lower than observed, reflecting the well-known difficulty stellar evolution models have in generating enough SNe-Ia per population to match observed SN-Ia rates and iron abundances (see e.g.Kobayashi et al. 2023).The impact of this discrepancy on the properties of galaxies in L-Galaxies is discussed in Section 5. Recent attempts to improve the modelling of mass transfer during binary evolution could improve the SN-Ia rate in models like bi-nary_c (Temmink et al. 2023).The impact of such approaches will be investigated using L-Galaxies in future work.

Implementation into galaxy evolution simulations
Traditionally, galaxy evolution simulations utilise stellar yields calculated for individual stars of various initial masses and metallicities.The overall mass ejection rate from a whole stellar population at time  is then calculated using the GCE equation (e.g.Tinsley 1980;Matteucci 2012), i.e. the integral over the death rate at time  of stars of mass  multiplied by the mass ejected by each of those stars (see e.g.section 4 of Yates et al. 2013).However, in this work we utilise the extension to binary_c called binary_c-python, which generates yields for entire stellar populations.In essence, binary_c-python calculates the integral of the yield ejected by each star multiplied by the IMF.Therefore, the GCE equation does not need to be solved again within L-Galaxies.Instead, only the synchronisation of these yields with the SFHs of model galaxies (see Section 2) is required.This process is described below, using the example of a single radial ring within the stellar disc of a model galaxy.
At each simulation timestep, the upper and lower edges of the th SFH bin ( upper,i and  lower,i ) are determined.Then, the normalised yield of element  ejected at time  by stars born in SFH bin  is calculated for each of the six metallicities () considered in binary_c, by numerically integrating the binary_c mass ejection rates between these two times, The true value of  x,i (, ) is then obtained for radial ring  by interpolating between the six values using the actual metallicity of the stars at their formation in L-Galaxies,  r,i .This true normalised yield is then multiplied by the SFR in the th SFH bin,  r,i , to obtain the true mass ejection rate.The mass ejection rates from all SFH bins are then added together to obtain the total mass ejection rate of element  at time  for ring , and the total mass of element  ejected in that timestep by stars in ring  is then where Δ is the simulation timestep width.
We note that carrying-out all these steps for every radial ring of every galaxy within the simulation would be computationally prohibitive.Therefore, we adopt the approach used for the GCE model in previous versions of L-Galaxies (Yates et al. 2013), by carrying-out the first part (i.e. up to Eqn. 2) in pre-processing.This takes advantage of the fact that the timestep spacings and SFH bin spacings are identical for all galaxies (and all radial rings) in L-Galaxies.Such a uniform time structure allows us to loop through all timesteps just once before the main simulation code is executed, to calculate  x,i (, )/ M ⊙ and store it in a look-up table with dimensions of timestep, SFH bin, initial metallicity, and chemical element.This look-up table is then used within L-Galaxies to solve Eqns. 3 and 4.

Comparing single and binary stellar population yields
Before assessing the impact of binaries in a complex galaxy evolution simulation, it is informative to first compare the basic ejection rates predicted by binary_c with those from single-star-only stellar populations.To do this, we compare three set-ups: (a) the mixed (i.e.binary+single star) stellar population from binary_c described in Section 3.1, (b) an equivalent single-star-only stellar population from binary_c, and (c) a single-star-only stellar population including an analytic prescription for SNe-Ia from the previous L-Galaxies GCE model, first presented by Yates et al. (2013, hereafter Y13).
This third set-up assumes a Chabrier ( 2003) IMF and includes Marigo (2001) yields for AGB stars, Thielemann et al. (2003) yields for SNe-Ia, and Portinari et al. (1998) yields for SNe-II, similar to many other cosmological-scale simulations.Metal enrichment from stars of mass 0.85 to 120 M ⊙ is assumed, along with a power-law SN-Ia delay-time distribution (DTD) of slope -1.12 and a SN-Ia factor of  SNIa = 0.035 (see Yates et al. 2021a for details).This  SNIa parameter determines the fraction of stellar systems in the mass range 3−16 M ⊙ that are SNe-Ia progenitor systems, and essentially corrects for any shortfall in the predicted yield of iron-peak elements by amplifying the expected SN-Ia fraction to compensate.The ejection rates for AGB winds, SNe-II, and SNe-Ia are then calculated using the standard GCE equation (eqn.5 in Y13).These are comparable to ejection rates from the binary_c set-ups using Eqn. 3 above.
We note that the different IMF upper-mass limits between the Y13 prescription (120 M ⊙ ) and that from binary_c (80 M ⊙ ) does not have a significant impact on the evolution of model galaxies.When changing this limit from 120 M ⊙ to 80 M ⊙ in the Y13 set-up, we find no significant changes to stellar and gas masses of galaxies at  = 0, and average ISM oxygen abundance only drops by ∼ 0.04 dex (due to having slightly fewer core-collapse SNe).
For the purposes of this comparison, the yields from binary_c have been grouped into the following three broad categories, to roughly match the three enrichment channels previously implemented into L-Galaxies: • Wind group: AGB, RLOF, Wind, Thorne-Zytkow, GB, CE • SN-II group: SNII, SNIIa, SNIbc, WR, SNeCAP, SNAIC • SN-Ia group: SNIaELD, SNIaChandCoal, SNIaChand, Novae These set-ups are then run through an arbitrary SFH (here a constant SFR of 1 M ⊙ /yr), to obtain the total mass ejection rate as a function of time.
Fig. 2 shows the total mass ejection rates (top panel) and metal mass ejection rates (bottom panel) for the mixed set-up from bi-nary_c (solid lines), the single-star-only set-up from binary_c (dashed lines), and the Y13 set-up (dotted lines), at an example metallicity of  = 0.004.
We can see that the total mass ejected by massive stars (i.e. from the SN-II group, orange) is quite similar between these three set-ups.Such a result was also found by De Donder & Vanbeveren (2004) when comparing single-star-only vs. binary-only stellar populations.This indicates that the total energy available for SN feedback in L-Galaxies will be relatively unchanged when switching from a single-star-only to mixed prescription (see Section 5.1).However, the bottom panel of Fig. 2 shows that the metal mass ejected by massive stars is slightly greater when using binary_c.This is partly due to the more efficient mass loss that massive stars undergo when in binary systems (comparing the dashed and solid orange lines), but is predominantly due to the different SN-II yield tables used, i.e.Chieffi &Limongi 2004 in binary_c andPortinari et al. 1998  binary_c.This is due to the under-estimation of the SN-Ia rate, as discussed in Section 3.1.4,and leads to a deficit in the production of elements such as Fe in L-Galaxies (see Section 5.4).
Ejection from SNe-II is largely independent of metallicity in bi-nary_c.Conversely, ejection rates from the Wind group increase with metallicity, predominantly due to increasing yields for elements heavier than N as opacity increases.There is also a large drop in the early-time ejection rate from SNe-Ia at both very low ( = 0.0001) and very high ( = 0.03) metallicity in binary_c, predominantly due to a decrease in prompt (∼ 100 Myr) SNe-Ia leading to a drop in the yield of heavy  elements and Fe.
We also note that binaries contribute up to 90 per cent of the elemental ejection from the 'SN-II group' in the binary_c mixed populations, mainly via Wolf-Rayet winds and traditional SNe-II.In the single-star-only scenario, most of this material is simply ejected by isolated massive stars instead.
Finally, Fig. 3 shows the individual chemical element ejection rates from a single 1 M ⊙ /yr burst of star formation with  = 0.004.There is a reasonable agreement between the AGB yields from Izzard et al. (2004Izzard et al. ( , 2006, used by binary_c) , used by binary_c) and Marigo (2001, used by Y13) at late times, although the Marigo (2001) yields tend to be higher (especially at lower metallicities) and more sensitive to star mass and therefore time.The inclusion of early stellar winds in binary_c extends its 'Wind group' ejection to much earlier times, and CE ejection provides a boost to the yields of all elements (particularly carbon) in binary_c when binaries are included.However, we note that the overall carbon ejection from the binary_c set-ups is ∼ 38 per cent lower than in the Y13 set-up at low metallicities (i.e. ≤ 0.001).This is predominantly caused by lower carbon yields from AGB stars.The minimum core mass and efficiency assumed for the third dredge up in the Karakas et al. (2002) models used in binary_c are low compared to those from Marigo (2001) used by Y13.This is because they are not explicitly calibrated to match carbon star luminosity functions (Z.Osborn, priv.comm.).Conversely, nitrogen is boosted by up to ∼ 76 per cent in the mixed binary_c model compared to Y13 at low metallicities.This is due to enhanced 14 N ejection from CE, novae, AGB stars, and SNe-Ia.However, this is no longer the case once metallicities reach the value of  ∼ 0.004 seen in Fig. 3.
There is little difference in the ejection rates from SNe-II at this metallicity, with the exception of a dip in the  element ejection rate at ∼ 5 Myr in the Portinari et al. (1998) single-star yields (used by Y13) compared to binary_c.This can be seen for the heavier  elements at  = 0.004 in Fig. 3 but is also present for lighter  elements at lower metallicities.The impact of this is discussed in Section 5.1.The late-time 'SN-II group' ejection from the binary_c mixed population is predominantly from WR winds, along with more minor contributions from accretion-induced collapse SNe (SNe-AIC) and the SNe-IIa of binary systems (see Zapartas et al. 2017).
However, there is a clear difference in the predicted SN-Ia group ejection rates for many elements.For example, N (mostly 14 N) is ejected much more significantly after ∼ 2 Gyr in the binary_c yields than those used by Y13 due to the inclusion of novae.Conversely, the total amount of Fe ejected by SNe-Ia is lower in binary_c by ∼ 1.0 dex at all times, as discussed above.The H and He ejection seen for the 'SN-Ia group' in binary_c is predominantly from novae, with an additional contribution of accreted He from SNe-ELD.
Versions of Figs. 1, 2, and 3 for other metallicities considered in binary_c are available online in the Supplementary Material.

DUST MODEL
The dust model implemented here is based on that developed for an earlier version of L-Galaxies by Vĳayan et al. (2019, hereafter V19).This includes the production of four main dust types (silicates, carbon, iron oxides, and silicon carbides), their distribution across the diffuse and molecular components within the ISM, and their destruction by various mechanisms.
In this work, we extend that model to include the survival and gradual destruction of dust in the hot CGM and Ejecta components surrounding galaxies, as well as the spatial distribution of dust within the radially-resolved ISM.The partitioning of dust into diffuse and molecular gas components in the ISM is now also fully integrated with our HI-H 2 partitioning model.
We model the transport of dust into the CGM in the same way as for metals; (a) partial ejection of newly-formed dust from stars and SNe in the stellar disc, (b) entrainment in subsequent galactic outflows, and (c) full ejection of newly-formed dust from stars and SNe in the bulge and stellar halo.We find that allowing dust to escape galaxies, particularly via direct enrichment from disc SNe, has an important impact on the dust masses and therefore grain growth rates inside galaxies at high  (see Sections 5.2 & 5.3).The key dust production and destruction processes modelled in L-Galaxies are outlined below.

Dust production and growth
In this work, dust is assumed to be produced by three sources: winds from asymptotic giant branch (AGB) stars, type Ia supernovae (SNe-Ia), and type II supernovae (SNe-II).In addition, the dust mass can increase through in-situ grain growth in the ISM.These processes are described in detail below.

Dust production by AGB winds
For AGB stars, we use the mass-and metallicity-dependent dust yields presented by Ferrarotti & Gail (2006).The total dust production rate is then given by, where  d,AGB (, ) is the dust mass ejected by one AGB star of mass  and birth metallicity , ( −  M ) is the star formation rate (SFR) at the birth of a star with lifetime  M , and () is the IMF which is assumed to be constant across all time and space.
The mass-and metallicity-dependent lifetimes from Portinari et al. (1998) are used for  M here.We note that dust is assumed to be formed at the end of the AGB star's life, and that we ensure the dust mass produced does not exceed the mass of the constituent chemical elements initially ejected by the AGB stars.

Dust production by SNe
For dust production by SNe-Ia and SNe-II, we follow the method recommended by Zhukovska et al. (2008), whereby the mass of dust formed is assumed to be proportional to the mass of the "key element" ejected by the SNe (Si for silicates and silicon carbides, C for carbon dust, and Fe for iron dust).The total dust production rate from SNe-Ia + SNe-II is therefore given by, where  is each of the four dust types considered,  is the corresponding key element,  ,SN is the ejection rate of key element  (which directly depends on the yields calculated by binary_c, see Section 3),  ,SN is the dust condensation efficiency parameter for dust type  (see table 1 of V19) which accounts for the destruction of newly-formed dust by the SN reverse shock, and   and   are the atomic weights of dust type  and its corresponding key element , respectively.For AGB stars and SNe, the fraction of dust ejected that goes into the molecular or diffuse components of the ISM is set such that their dust mass ratio,  dust,clouds / dust,diff , is unchanged.

Grain growth in the ISM
Growth of dust grains through accretion in the ISM is assumed to only take place inside molecular clouds.In L-Galaxies, the H 2 fraction in the ISM is calculated at each timestep (for each radial ring) following the model developed by Fu et al. (2010Fu et al. ( , 2013)).This model uses the metallicity-and density-dependent H 2 partitioning formalism of Krumholz et al. (2009), as discussed in section 2.2.3 of Henriques et al. 2020.We note that a metallicity floor of  = 0.01 is assumed when calculating the gas clumping factor in this formalism, in order to stimulate H2 formation at very early times.
At the start of each timestep, we assume that the fraction of dust in molecular clouds is the same as the H 2 fraction, thereby ensuring that the chemical composition of diffuse and molecular components within each radial ring is initially the same.This is reasonable given that the expected lifetime of a molecular cloud (∼ 10 Myr) is comparable to the timestep widths in L-Galaxies (< 20 Myr).
Once the total mass of molecular material is known, the maximum amount of each element in molecular clouds available for accretion onto dust grains is calculated.This is known as the maximum condensation fraction,  max .For the refractory elements Mg, Si, Ca, and Fe, we assume  max = 1.0.For the volatile elements N, Ne, and S, we assume  max = 0.0.For C, we assume an  max equal to the fraction of C in the molecular gas phase that is not in CO molecules.In essentially all cases, the CO fraction for carbon is set to 0.3 (leading to  max,C = 0.7), in line with observations of molecular clouds in the Milky Way (Irvine et al. 1987;Van Dishoeck & Blake 1998).However, it can be lower in environments where the oxygen abundance is particularly low.
The maximum condensation fraction for O,  max,O , can then be determined at each timestep from the abundance of other refractory elements left available for grain growth (excluding oxygen atoms in CO), assuming that O forms only into silicate and iron oxide dust.
To calculate this, we follow Zhukovska et al. (2008) such that where  O,dust,clouds is the mass of O already in dust in molecular clouds and  O,clouds is the total mass of O in molecular clouds.
Once the maximum condensation fractions are known, the actual amount of each refractory element that accretes onto dust grains in that timestep can be calculated.To do this, we need to define (a) the timescale for the exchange of material between molecular and diffuse gas (i.e. the effective molecular cloud lifetime),  ex , and (b) the timescale for accretion of gas-phase elements onto existing dust grains,  acc .For the exchange timescale, we assume  ex = 10 Myr (see Zhukovska et al. 2008).For the accretion timescale, we assume, This accretion timescale formalism differs from that presented by Asano et al. (2013) and commonly used in other cosmological simulations (e.g.De Bennassuti et al. 2014;Popping et al. 2017;Triani et al. 2020;Graziani et al. 2020), in that the efficiency of grain growth here depends on the mass of dust present in molecular clouds, rather than on gas-phase properties such as the ISM metallicity and density.This means that our efficiency parameter,  acc,0 , takes on a slightly different meaning to other works.Nonetheless, as in those works (but not V19), we tune our  acc,0 to obtain log( acc /yr) ∼ 6.08 at  = Z ⊙ = 0.0134 and  mol ∼ 250 cm −3 , as expected from the Asano et al. (2013) formalism.This returns  acc,0 = 5 × 10 3 yr, which we assume to be fixed across space and time.This calibration leads to a wider dynamic range of accretion timescales than found in most other simulations.For example, at fixed H 2 density, it produces shorter  acc in dust-rich environments and longer  acc in dust-poor environments than seen in Popping et al. (2017).This has an important impact on the dust-to-metal ratios predicted for galaxies (see Section 5.2).
The first panel of Fig. 4 shows the relation between ISM oxygen abundance and  acc at  = 0 for the two versions of L-Galaxies we consider in this work.These versions are defined in detail in Section 5.
The fraction of the available mass of a given chemical element that accretes onto dust grains within one timestep is then given by, where Δ is the timestep width and a constant  acc is implicitly assumed across the timestep.We also simultaneously solve for the fraction of material in molecular clouds that isn't yet exchanged with diffuse gas, given by where  =  H2 /( HI +  H2 ) is the molecular gas fraction.In the limit where  acc ≫  ex , the dust fractions in the molecular and diffuse phases are identical after grain growth.However, for cases where  acc ≲  ex (i.e. in dustier environments), the dust fraction in molecular clouds will exceed that present in the surrounding diffuse ISM after grain growth (see also section 3.2 of V19).This model returns values for  C,dust,clouds of typically 0.6 − 0.7 at  = 0.There is a scatter down to lower values at high metallicity in early-type galaxies (ETGs), which have had most of their interstellar dust removed or destroyed by  = 0.For oxygen,  O,dust,clouds is always ≲ 0.2, which can has a noticeable impact on measured gasphase oxygen abundances in galaxies (see Section 5.1).
We note that the grain growth efficiency could also be affected by other properties, such as the grain size distribution (Aoyama et al. 2020) and charge distribution (Glatzle et al. 2022).Such dependencies will be investigated in L-Galaxies in future works.

Dust destruction by SN shocks and astration
Following McKee (1989), dust in the ISM is assumed to be destroyed via sputtering in SN-induced shocks.The destruction timescale for this process is where  cold is the mass of cold gas within each ISM radial ring,  SN = 0.36 is the expected fraction of all SNe that interact with the ISM (see section 4 of McKee 1989), and  SN is the combined rate of SNe that explode within each radial ring at each timestep.
cleared represents the mass of cold gas assumed to be cleared of dust by destruction from a single SN.In V19, this value was fixed at 1200 M ⊙ .In this work, we instead follow the findings of Hu et al.
(2019) from their high-resolution hydrodynamical simulations by allowing  cleared to vary between 1180 and 1660 M ⊙ depending on the composition of the dust, with the lower limit corresponding to  C,dust /( C,dust +  Si,dust ) = 1 and the upper limit corresponding to  C,dust /( C,dust +  Si,dust ) = 0.For typical dust-to-gas ratios of ∼ 0.001 at  = 0 (De Vis et al. 2019), this leads to a typical dust-destruction rate of 1.18 − 1.66 M ⊙ per SN, which is consistent with the upper limit of 3 M ⊙ estimated by Ferrara & Peroux (2021).
We then compute the mass of dust in the ISM that is destroyed per timestep as The second panel of Fig. 4 shows  shock as a function of overall gas-phase oxygen abundance for regions of the ISM within model galaxies at  = 0. Most regions exhibit destruction timescales around 1 Gyr.This is in reasonable agreement with the values predicted for the solar neighbourhood by the hydrodynamical simulations of Slavin et al. (2015) and Hu et al. (2019).The scatter up to very large values of  shock is predominantly caused by regions with low current SFR and therefore low SN-II rates.At fixed metallicity,  shock is longer in MM+binary_c than in MM+singleStars due to lower overall SN rates (see Section 5.3).
In addition to SN shocks, dust is also destroyed in the ISM through astration.In L-Galaxies, all dust present in the molecular gas used to form stars is assumed to be destroyed, with it's constituent elements becoming incorporated into the newly-formed stars.

Dust destruction by thermal sputtering
We also consider the survival of dust in galactic outflows and its subsequent gradual destruction via thermal sputtering in the CGM and Ejecta phase surrounding galaxies.Dust in extragalactic hot gas has been observed out to large distances, with a spatial distribution similar to the CGM itself (e.g.Ménard et al. 2010;Peek et al. 2015).Our model for thermal sputtering is based on that of Tsai & Mathews (1995) The thermal-sputtering timescale depends on the hot gas density ( hot ) and temperature ( hot ) by  sput = 0.17 Gyr where the initial grain radius at each timestep is assumed to be  0 = 0.1m, the factor controlling the sputtering timescale in the low-temperature regime is fixed at  = 2.5, the critical temperature above which the sputtering rate becomes roughly constant is  crit = 2 × 10 6 K, and the hot gas is assumed to be at the virial temperature given by  hot ≡  vir = 0.5 H  2 vir / B .The constants  H and  B are the mean atomic mass and the Boltzmann constant, respectively.The dust-grain radius after sputtering, , is calculated as where  sput = 3.2 × 10 −18 cm 4 s −1 , and  p is the mass of a proton.The mass of dust destroyed through sputtering per timestep is then given by where the factor of three arises from the grain mass ( g ) being related to the grain radius and density ( g ) by  g = (4/3) 3  g , as pointed out by McKinnon et al. (2017).
When assessing how to measure  hot for use in Eqn. 13, we note that assuming a uniform density of  hot =  hot /(4/3) 3 vir (i.e. the single isothermal sphere, or SIS, approximation) would likely underestimate the true average in the CGM, as both gas and dust are observed to exhibit higher densities in the core than the outskirts (e.g.Vikhlinin et al. 2006;Ménard et al. 2010).Therefore, we approximate the gas density profile for Eqn. 13 with a  model given by () = [1+ (/ c ) 2 ] −3 s /2 (Cavaliere & Fusco-Femiano 1976).Such a profile is commonly used to approximate the hot-gas density distribution in galaxy groups and clusters in observations (e.g.Reiprich & Böhringer 2002) and semi-analytic simulations (Yates et al. 2017).Here, we assume the canonical shape parameter of  s = 2/3 and that the core radius,  c , equals the scale length of the dark-matter halo assuming an NFW profile (Navarro et al. 1997).The core density is then obtained from  hot,0 =  hot / ∫  vir 0 4 2 ()d, and the density in Eqn. 13 as ρhot =  hot,0 We find that this assumption leads to shorter dust sputtering timescales in the CGM surrounding lower-mass galaxies compared to the simpler SIS approximation.However, we caution that the degree of this difference varies depending on the assumed gas density profile.Therefore, we take our preferred model as a reasonable intermediate case, which allows additional system-dependent flexibility in the CGM sputtering timescale.
For the Ejecta component surrounding haloes, we assume the simpler SIS density when calculating  sput , as this is likely to better represent the lower ambient densities found in the intergalactic medium (IGM).
The third and fourth panels of Fig. 4 show  sput for the CGM and Ejecta components, respectively, as a function of stellar mass for galaxies of all types at  = 0. Sputtering timescales are much shorter in high-mass systems because of their higher core densities and virial temperatures.At log( * / M ⊙ ) ≳ 11.0,  sput tends to be shorter than the galaxy-averaged  shock within the ISM, meaning that thermal sputtering outside galaxies dominates the dust destruction in such systems (see Section 5.3).

RESULTS
To assess how binary stars affect galaxy evolution, we compare two versions of the L-Galaxies MM simulation described in Section 2: (a) a version using the mixed stellar populations from binary_c presented in Section 3.1 (hereafter, MM+binary_c), and (b) a version using the single-star stellar populations plus SNe-Ia formalism presented in Section 3.3 (hereafter, MM+singleStars).These two versions are chosen to best compare the new and old L-Galaxies simulations, and our new dust model is included in both.Unless stated otherwise, in this section we combine model galaxies from runs on both Millennium and Millennium-II, with each galaxy weighted by the inverse of the effective volume of the underlying N-body simulation, considering only "resolved" systems above log( * / M ⊙ ) = 8.0 from Millennium and log( * / M ⊙ ) = 7.0 from Millennium-II.For plotting, we use representative sub-volumes of the Millennium and Millennium-II boxes, containing a total of ∼ 495, 000 model galaxies by  = 0.

General galaxy properties
We first present the key galaxy scaling relations used to assess the success of cosmological-scale simulations.Fig. 5 shows the stellar mass function (SMF), HI mass function (HIMF),  * -sSFR relation, SN rates, and mass -metallicity relations (MZ g R and MZ * R) for galaxies at  = 0. Observational data (orange) is shown alongside the MM+binary_c (blue) and MM+singleStars (red) versions of L-Galaxies.
The main conclusion from this comparison is that the use of binary stellar populations does not appear to significantly affect the general properties of galaxies at low .For example, panels (a), (b), and (c) of Fig. 5 show that the SMF, HIMF, and galaxy sSFRs are essentially unchanged when switching from the Y13 set-up to the binary_c set-up in L-Galaxies (see Section 3.3).This is chiefly due to the similar mass ejection rates from massive stars in these two set-ups (as shown in Fig. 2), which determines the total SN feedback energy in the simulation.A similar conclusion was also drawn by De Donder & Vanbeveren (2004), who found that binaries have a negligible impact on the star formation history of Milky-Way-type galaxies in one-zone GCE models with no outflows.
Panel (d) shows that SN-II rates are also in reasonable agreement between both the L-Galaxies versions and observations.However, SN-Ia rates are severely under-estimated in MM+binary_c, as expected from the yield analysis in Section 3.3.This has an impact on the iron abundances of galaxies, which is discussed in Section 5.4.Conversely, MM+singleStars is in good agreement with the observed SN-Ia rate from the LOSS survey Graur et al. (2017), primarily because the fraction of SN-Ia systems per stellar population (i.e. the  SNIa parameter) is scaled to match global chemical properties at  = 0.
Panel (e) shows the MZ g R for star-forming disc galaxies at  = 0, where gas-phase metallicity is calculated as the SFR-weighted average ISM oxygen abundance (excluding oxygen in dust) in units of 12+log(O/H) for L-Galaxies, where O/H =  O / H = ( O / H ) • (  H / O ) and  H and  O are the atomic weights of H and O, respectively.This method best matches the emission-line methods used for the observations to which we compare.
There is a slight but systematic difference in gas-phase metallicity between the two versions of L-Galaxies, especially for lowmass galaxies.MM+binary_c returns higher 12+log(O/H) than MM+singleStars by ∼ 0.1 dex on average, and is in marginally better agreement with the latest metallicity observations (although see Cameron et al. 2023).This difference is predominantly caused by differences in the SN-II yield tables assumed, rather than binary effects.In the Portinari et al. (1998) yields used in MM+singleStars, there is a drop in the O ejection rate from low-metallicity stars of mass ∼ 50 − 60 M ⊙ due to the build-up of large CO cores, which leads to larger remnant masses and therefore a lower ejected yield (see Portinari et al. 1998, section 5).This effect can also be seen as a dip in the single-star burst yields in Fig. 3 for most heavy elements ∼ 5 Myr after star formation.In the Chieffi & Limongi (2004) yields used in MM+binary_c, this effect is weaker, leading to higher O enrichment in low-metallicity environments, and consequently higher 12+log(O/H) in the ISM by  = 0 in low-mass galaxies.
It is interesting to note that this boost in 12+log(O/H) at low mass in MM+binary_c is almost exactly offset by the drop in 12+log(O/H) we find when accounting for the depletion of oxygen onto dust grains.This means that MM+binary_c maintains the good agreement with the observed MZ g R at  = 0 that was found for previous versions of L-Galaxies that did not account for dust depletion.
We also note that metallicities are not converged between Millennium and Millennium-II.For example, gas-phase metallicities at  = 0 are, on average, ∼ 0.1 − 0.3 dex higher at fixed stellar mass when L-Galaxies is run on the higher-resolution Millennium-II.This is due to the known effect of higher SFRs and gas masses in lower-mass galaxies for Millennium-II (Henriques et al. 2020;Yates et al. 2021b), which in turn leads to higher metal and dust abundances at larger radii in the ISM.
Finally, panel (f) shows the MZ * R, where  * is the total massweighted stellar metallicity measured within a 3 arcsec aperture centred on each galaxy for both observations and L-Galaxies.These metallicities are normalised to the bulk metallicity of the Sun as measured by Asplund et al. (2009).Interestingly, higher stellar metallicities are not seen in MM+binary_c within 3 arcsec.This illustrates that the higher oxygen abundances discussed above are predominantly at larger galactocentric radii in  = 0 galaxies, where the metallicity of the star-forming gas is lower.

Dust scaling relations
We now consider the impact of binary stellar evolution on the dust properties of galaxies at low and high .Fig. 6 shows the dust-to-metal ratio (DTM =  dust / metal,tot , top panel), dust-to-gas ratio (DTG =  dust / HI+H2 , middle panel), and dust mass ( dust , bottom panel) in the ISM of galaxies at  = 0 as a function of their SFR-weighted average gas-phase oxygen abundance for the MM+singleStars (red) and MM+binary_c (blue) versions of L-Galaxies.
There is a very good correspondence between both versions of L-Galaxies and observational data (orange points) in Fig. 6.At lower metallicities, galaxies undergo an increase in their DTM, DTG, and  dust over time in L-Galaxies, due to the increasing efficiency of grain growth as the dust mass in molecular clouds increases.At higher metallicities, a typical value of log(DTM) ∼ −0.6 is found for star-forming systems at  = 0, along with a significant scatter down to low DTM which is also seen in the DustPedia observations (De Vis et al. 2019, diamonds).In L-Galaxies, this scatter at high metallicity contains mostly quenched galaxies, which have undergone net dust destruction at late times.
When comparing the two versions of L-Galaxies, we find that .In all cases, the MM+binary_c (blue) and MM+singleStars (red) versions of L-Galaxies are shown.Overall, stellar masses, gas masses, and star formation rates are largely unaffected by switching from our single-star to binary_c set-up, however gas-phase metallicities are higher at  = 0 by ∼ 0.1 dex on average.
MM+binary_c returns slightly higher 12+log(O/H) at fixed DTM, DTG, and  dust by  = 0.This is due to low-mass galaxies being more oxygen rich in MM+binary_c (see Section 5.1), which has a greater effect on the gas than the dust in our model due to the lower maximum condensation factor for oxygen than other refractory elements (see Section 4.1.3).This effect shifts the low-mass galaxy population right and down in the top panel and right in the middle and bottom panels.
In Fig. 7, we show the evolution back to  ∼ 4 of the median DTM and DTG ratios as a function of gas-phase oxygen abundance (top two rows) and median  dust as a function of galaxy  * (bottom row).A range of observational data is included for comparison (see figure caption).A similar observational trend is seen at high  when only considering H 2 mass in the DTG ratio (Saintonge et al. 2013;Popping et al. 2023).We note that oxygen abundances for the QSO and GRB DLA data from Péroux & Howk (2020, squares) and Heintz et al. (2023, pentagons) were derived from their measured metalto-hydrogen mass ratios as follows: 12 + log(O/H) = 12 + [M/H] + log( O,⊙ / H,⊙ )+log(  H / O ).This approach assumes the fraction of all metal mass that is oxygen in these systems is the same as in the Sun, i.e. ( O / Z ) ⊙ = 0.43 from Asplund et al. (2009), which may under-estimate the true 12+log(O/H) at high  by ∼ 0.1 dex.In order to better compare to these absorption-line-based measurements, we plot mass-weighted average gas-phase oxygen abundances for L-Galaxies in the top two rows of Fig. 7.
We find that both L-Galaxies versions match observations well back to  ∼ 4, although the slightly lower DTM and DTG ratios found in MM+binary_c (blue) at low  appear to match the data better.An evolution in both the DTM and DTG ratios for metal-poor galaxies is also found, similar to that seen in GRB samples (e.g.Wiseman et al. 2017;Heintz et al. 2023).This is a significant success of the new L-Galaxies simulation, as previous cosmological-scale simulations have struggled to simultaneously match these dust properties across this redshift regime (Ginolfi et al. 2018;Popping & Péroux 2022).The key reason for this success is the shorter accretion timescales returned by our grain-growth formalism in dustier environments, which allow more efficient dust production at later times.Nonetheless, DTM ratios in metal-poor galaxies at  ≳ 4 appear under-estimated in L-Galaxies compared to absorption-line measurements from QSO and GRB DLAs (squares and pentagons).This is due to the inefficiency of grain growth in dust-poor environments in our model, which suggests either the need for enhanced dust production from other sources or a more detailed modelling of the dependence of the dust production rate on environmental properties (Graziani et al. 2020).A similar result has also been found for the CROC simulations (Esmerian & Gnedin 2022, 2023) using a simpler dust model in post-processing.However, we note that the difficulty in detecting galaxies in this regime leaves open the possibility that current observational samples are biased towards dustier systems.The sharp increase in DTM ratio seen at low metallicity in L-Galaxies is driven by the transition between SN-dominated and grain-growth-dominated dust production.At low metallicities, dust is predominantly produced by SNe and AGB stars.However, once the metallicity, and thus in-situ dust mass, begins to increase, grain growth becomes increasingly efficient, causing a rapid increase in the DTM ratio until a rough equilibrium between grain growth and dust destruction is reached.This transition can also be seen in the grain growth accretion timescales of galaxies in Fig. 4. A similar effect is found in the DustyGadget hydrodynamical simulation developed by Graziani et al. (2020), but is not present in most other cosmologicalscale simulations (Popping & Péroux 2022).Other simulations tend to assume longer accretion timescales at low H 2 densities, leading to less efficient grain growth at late times (see Section 4.1.3).Such a sharp transition in DTM ratio with metallicity is unconfirmed in observations, however the upper limits obtained via absorption-line measurements (Péroux & Howk 2020, grey triangles) suggest that such a drop is not ruled out.
Finally, the bottom row of Fig. 7 illustrates that both versions of L-Galaxies return very similar  * - dust relations at all redshifts.This indicates that net dust production from AGB stars and SNe is largely independent of the effects of binary stellar evolution.Both versions of L-Galaxies are also in marginal agreement with the dust masses measured by Santini et al. (2014) andDa Cunha et al. (2015) back to  ∼ 4. We would expect these observational data to lie in the upper envelope of the model relation, as they include the most star-forming, dust-rich systems.

Dust production and destruction rates
The top panel in Fig. 8 shows the evolution of various cosmic dust growth and destruction rate densities for MM+binary_c (solid lines) and MM+singleStars (dashed lines) run on Millennium-II.This higher-resolution N-body simulation is preferred here, as it resolves more of the low-mass galaxies which dominate cosmic densities at high  (see Yates et al. 2021b).For this reason, there is also no stellar mass cut applied for this analysis.
We can see significant differences in the evolution of some rates between the two versions of L-Galaxies.For example, the evolution of the cosmic destruction rate in the ISM (cyan), which is dominated by SN shock destruction (  shock ), is flatter in MM+binary_c than in MM+singleStars.This is because, at high ,  shock is driven by its dependence on  dust (see Eqn. 12), which is higher at this epoch in MM+binary_c due to higher dust and metal yields from lowmetallicity SNe-II (see Section 3.3).Whereas at low ,  shock is limited by the dust destruction timescale (i.e. shock ), which is always longer in MM+binary_c due to lower SN rates (see Fig. 4).
A similar trend is seen for the cosmic grain growth rate (  acc , red).At high , there are more galaxies with high metal and dust masses in MM+binary_c, leading to shorter accretion timescales (see Eqn. 8) and thus higher  acc .At low , there is less mass available for accretion onto dust grains in MM+binary_c, due to the lower dust destruction rates described above.In other words, grain growth saturates sooner in MM+binary_c than MM+singleStars.Indeed, when turning-off dust destruction in L-Galaxies, we find that  acc is very similar between MM+binary_c and MM+singleStars, indicating that grain growth is limited by  shock at low .In addition, accretion timescales are slightly longer in MM+binary_c by  = 0, due to lower DTG ratios (see Fig. 4).This coupled evolution between the main dust production and destruction mechanisms  The stellar mass -dust mass relation.In all cases, the median relation is shown (solid lines) along with the 1-and 2- spread (shaded regions), with metallicity measured as the mass-weighted oxygen abundance.Observational data is from Rémy-Ruyer et al. (2015, circles), the DustPedia data set (De Vis et al. 2019, diamonds), the collection of Péroux & Howk (2020, squares, with upper limits as grey triangles), Heintz et al. (2023, pentagons), Shapley et al. (2020, stars), Santini et al. (2014, pluses), andDa Cunha et al. (2015, crosses, considering only those galaxies with photometric redshift errors < 0.5), binned in redshift to match the L-Galaxies outputs ±0.5 (see labels).Both versions of L-Galaxies match observations well back to  ∼ 4, although dust masses may be under-estimated in metal-poor galaxies at ≳ 3. nicely demonstrates the rapid cycling of dust in the ISM occurring in L-Galaxies.
The cosmic SN-II dust production rate (orange) is higher at all epochs in MM+binary_c due to its dependence on the mass of metals ejected by SNe-II, which is higher than in MM+singleStars.Conversely, the SN-Ia dust production rate (green) is significantly lower in MM+binary_c due to the under-estimation of the SN-Ia rate in binary_c (see Section 3.3).
We also find that the total dust production and destruction rates in MM+binary_c also exceed those in MM+singleStars above  ∼ 4.4 and ∼ 5.4, respectively.
In both versions of L-Galaxies, grain growth overtakes SNe-II as the dominant dust production mechanism at  ∼ 6.4.This "crossover redshift" is later than predicted by most other cosmological simulations (e.g.Popping et al. 2017;Vĳayan et al. 2019;Graziani et al. 2020;Dayal et al. 2022;Parente et al. 2023), likely due to our lessefficient early grain growth.A notable exception is DustySAGE (Triani et al. 2020), in which a crossover redshift of  ∼ 2.2 was found, due to lower grain growth efficiencies and a metallicity-dependent dust destruction efficiently (see Popping & Péroux 2022).
The bottom panel of Fig. 8 compares L-Galaxies to observational measurements of the cosmic dust mass density,  dust , from SED fitting (orange circles, Dunne et al. 2011;Beeston et al. 2018;Driver et al. 2018;Pozzi et al. 2020) and neutral gas absorption lines in QSO-DLAs (orange squares, Péroux & Howk 2020).At low , the total cosmic dust mass begins to level-off in both versions of L-Galaxies, but does not decrease.This is because dust continues to build-up in the CGM and ejecta reservoirs surrounding galaxies, with the latter actually becoming the dominant dust reservoir by  = 0.However, there is a net decrease in the ISM dust mass at late times, In both panels, values at discreet output redshifts have been fit with 5th order polynomials.Dust evolution in the two versions of L-Galaxies differs significantly at high , predominantly due to differences in ISM metallicity.
due to the combination of decreased cosmic star formation (thus decreased stellar dust production) and decreased ISM gas reservoirs (thus decreased grain growth), as discussed by Yates et al. (2021b).This leads to a drop in  dust within molecular clouds (dark blue) by a factor of ∼ 1.6 from  ∼ 1 to 0. However, this factor is less than the factor of 3 − 4 observed from the attenuation measurements using SED fitting.
To help reconcile this difference, Parente et al. (2023) have proposed allowing more efficient SMBH growth, and thus stronger AGN feedback, in simulations, to reduce star formation at late times and thus cause a stronger drop in  dust .However, in our case, this could also lead to an under-estimation of dust masses at  = 0, which match observations well (Fig. 6).Alternatively, an increase in star-formation at  ≳ 1 would cause a similar drop in  dust by boosting dust production at earlier times.The increased gas-cooling rates this requires would also stimulate dust production while leaving unaffected the good match L-Galaxies has to observed DTG ratios back to  ∼ 4 (M.Parente, priv. comm.).Indeed, larger gas reservoirs in resolved subhaloes at high  are already required in semi-analytic simulations to better match the observed cosmic Hi density (see Yates et al. 2021b), as well as the cosmic SFR density from sub-mm and radio observations.Therefore, further investigation into the viability of such a modification in L-Galaxies will be the focus of future work.
Finally, at high , L-Galaxies returns significantly lower  dust than is implied from DLA observations (orange squares).This discrepancy could be even larger if very-dust-rich DLAs, which are typically missed from optically-selected QSO samples, are also considered (Krogager et al. 2019).Part of this discrepancy could be mitigated by the increased gas cooling and star formation discussed above.Increased dust production, including from additional binary sources such as CE ejection and novae, or superluminous SNe (Chen et al. 2021), could also help.However, a residual offset in the cosmic metal mass density between L-Galaxies and QSO-DLAs was also found by Yates et al. (2021b), even when accounting for the simulation's gas deficit in resolved subhaloes.That work determined that this remaining difference could be reconciled if current high- DLA metallicity samples are biased towards galaxies with log( * / M ⊙ ) ≳ 7.5.It is possible that such a mass bias could also play a role in an over-estimation of  dust here.

Radial profiles
We now turn to studying the spatial distributions of dust and metals in the discs of massive star-forming galaxies at  = 0.For this, we select model galaxies with log( * / M ⊙ ) ≥ 9.0,  bulge / * ,tot < 0.3, and log(sSFR/yr −1 ) ≥ − 10.9 from the MM+binary_c and MM+singleStars versions of L-Galaxies run on Millennium.This larger N-body simulation is preferred here, as it provides better statistics while also being able to comfortably resolve higher-mass galaxies.The effective radius,  e , is measured as the half light radius in the r-band for L-Galaxies, to best match that measured in the observations to which we compare.We also note that "local" processes, such as star formation, dust production and destruction, and SN feedback, are calculated independently in each radial ring in L-Galaxies.This means that radial profiles are allowed to evolve in an unparameterised way (i.e.we do not impose any particular radial distribution function).However, we do assume an exponential radial profile for the accretion of gas and dust onto galaxy discs (see Henriques et al. 2020, section 2.2.1).
The top row of Fig. 9 shows the median ISM gas-phase oxygen abundance profiles (excluding oxygen in dust) for star-forming disc galaxies at  = 0 in three bins of stellar mass.The observational data shown (orange points) are from the MaNGA sample of Yates et al. (2021a) and the MUSE sample of (Erroz-Ferrer et al. 2019) with spaxel-based metallicities re-calculated by Easeman et al. (in prep.).In both cases, the Dopita et al. (2016) strong-line metallicity diagnostic is used, as this best reproduces electron-temperature metallicities in low-mass galaxies, super-solar metallicities in very high mass galaxies, and has a negligible secondary dependence on ionisation parameter (see Yates et al. 2021a, section 3.2).
We find very good agreement between both L-Galaxies models and the high-resolution MUSE data, in all mass bins.In Yates et al. (2021a), we reported that the L-Galaxies profiles are too flat at large  Asplund et al. (2009).Third row: The DTM ratio in the ISM for molecular clouds (solid lines) and diffuse gas (dashed lines).Fourth row: The DTG ratio.Fifth row: Dust mass surface density, with observational data of six galaxies from the DustPedia dataset of Casasola et al. (2017).Both versions of L-Galaxies match observational data well, and predict a significant drop in DTM and DTG ratios at large radii.
radii in low-mass galaxies compared to this MUSE sample.However, this discrepancy is now greatly reduced by the careful re-analysis and removal of diffuse ionised gas (DIG) dominated spaxels from these data (Easeman et al. in prep.).This provides a purer sample of HII region measurements, which return higher metallicities at large radii in low-mass galaxies.
At high mass, the two versions of L-Galaxies also exhibit similar 12+log(O/H) profiles to each other.However, at low mass, MM+binary_c returns higher 12+log(O/H) (by up to ∼ 0.15 dex), particularly at large galactocentric radii.Similar trends are seen at large radii in the stellar metallicity profiles (not shown).This extra oxygen is a consequence of the smaller remnant masses predicted by the yield sets used in binary_c at low metallicity (see Section 5.1).We note that the strongest suppression of oxygen from the Portinari et al. (1998) yields used in MM+singleStars occurs below  ∼ 0.004, which is roughly equivalent to 12+log(O/H) < 8.4.
Conversely, the gas-phase [Fe/H] profiles (not shown) exhibit an under-abundance of iron in MM+binary_c compared to MM+singleStars, particularly at low radii.This is due to the lower SN-Ia rates in MM+binary_c.The combination of these two factors (namely, more oxygen and less iron in the MM+binary_c model) is most clearly seen in the second row of Fig. 9, which shows the median azimuthally-averaged [O/Fe] profiles for bulge plus (thick + thin) disc stars of all ages.These profiles are flatter in MM+binary_c, with a normalisation that is higher than expected in, for example, starforming disc galaxies of Milky Way mass (see e.g.Hayden et al. 2015;Rojas-Arriagada et al. 2019).This result suggests that measuring the  enhancements in the centres of disc galaxies could provide useful constraints on stellar yields.
The third row of Fig. 9 shows the DTM profiles for molecular clouds (solid lines) and diffuse gas (dashed lines) in the ISM.In both L-Galaxies versions, the diffuse gas DTM ratio drops at larger radii.This indicates that the common assumption of a universal value (see e.g.Rémy-Ruyer et al. 2015) does not hold on sub-galactic scales.We find that the degree of this drop depends on (a) the efficiency of grain growth, and (b) the amount of dust that is able to re-accrete onto the ISM after being ejected by galactic outflows.A similar conclusion was made by De Vis et al. (2019) when interpreting DustPedia data.In the extreme case where no dust is allowed to survive in the CGM, values for the DTM and DTG in the outer disc can become very low in L-Galaxies, as re-accretion cannot supplement the low in-situ dust formation rates at these radii.Central DTM ratios are also lower in the MM+binary_c model than the MM+singleStars model, due to lower dust masses.This dust deficit is what drives the lower global DTM ratios in the MM+binary_c model seen in Figs. 6 and 7.
We also note that the apparent difference between the average molecular and diffuse DTM profiles at large radii is due to averaging effects when stacking profiles.In L-Galaxies, the few galaxies which have molecular gas at large radii are relatively evolved and therefore dust rich, whereas galaxies with only diffuse gas at large radii are relatively dust poor.
The fourth row of Fig. 9 shows DTG profiles for the same model disc galaxies.There is a clear decrease in DTG ratio with radius for all galaxies, including intermediate-mass systems which only exhibit the commonly-assumed universal DTG ratio of ∼ 0.01 in their innermost regions.There are also slightly lower DTG ratios at low radii in MM+binary_c than MM+singleStars, due to lower central dust masses.The radial profiles of general properties such as Hi, H 2 , SFR, and  * densities are very similar in both versions of L-Galaxies, reflecting their similar global values shown in Fig. 5.
Finally, the fifth row of Fig. 9 shows dust mass surface density, Σ dust , profiles for L-Galaxies.These are compared to six DustPedia  galaxies from the sample of Casasola et al. (2017) for which effective radii could be obtained from the literature, namely NGC 2403 (circles), NGC 5457 (squares), NGC 5194 (diamonds), NGC 628 (pentagons), NGC 4736 (pluses), and HGC 5055 (crosses).Spatially resolved properties were obtained by the DustPedia team for these galaxies using multi-wavelength photometric maps.
There is a very good qualitative agreement between the observed Σ dust profiles and those from L-Galaxies beyond / e ∼ 1.The flatter inner Σ dust profiles found in most DustPedia galaxies appear in conflict with the steeper central oxygen abundance profiles seen in MaNGA and MUSE data (top row of Fig. 9).This could suggest either enhanced dust destruction (possibly by AGN feedback or older stellar populations in the bulge) or suppressed star formation and thus dust production in galaxy centres.
A further comparison to spatially resolved observations is shown in Fig. 10 where we compare the MM+binary_c model at  = 0 (lines) to stacked radial profiles from the ten HERACLES/THINGS galaxies presented by Abdurro'uf, et al. (2022, points).Similar results are found for the MM+singleStars model.There is remarkable agreement in both slope and normalisation for the Σ * , Σ H2 , and Σ dust profiles between the observations and L-Galaxies.However, the central SFR and Hi densities are slightly over-predicted in L-Galaxies compared to these observations.As for the comparison to DustPedia galaxies above, this could be due to ∼ 50% of the HERACLES sample hosting AGN, noting that most of their AGN hosts have flattened SFR profiles.Direct suppression of star formation via ISM gas removal from AGN feedback is not currently modelled in L-Galaxies, which instead only allows AGN to indirectly suppress star formation by reducing further cooling onto the disc.

CONCLUSIONS
In this work, we implement models for (a) binary stellar evolution (BSE) using the binary_c code and (b) dust production and destruction into the cosmological semi-analytic galaxy evolution simulation, L-Galaxies.This opens up a whole new range of studies for cosmological-scale simulations, including comparisons to the growing plethora of observations of chemical abundances and dust properties in galaxies at low and high .We hope that this work will serve as a template for implementing BSE sub-grid models into other galaxy evolution simulations in future.
When comparing stellar populations produced by binary_c to those from a simpler single-star plus SNe-Ia set-up (Section 3.3), we find that the total mass ejected by massive stars is relatively unchanged.This is due to the similar assumptions made about explodability in both set-ups.However, the metal mass ejected by massive stars is slightly larger in binary_c.This is partly due to the more efficient mass loss in binary systems, but mostly due to the different SN-II yields assumed between the set-ups.Binary_c also predicts significant ejection of light alpha elements (especially C) during the common envelope phase -a process which is not possible when considering only single stars.However, the ejection of heavier elements such as iron from SNe-Ia is significantly under-predicted by binary_c, due to an under-estimation of the SN-Ia rate.This issue is common to stellar evolution codes which attempt to model this process from first principles.
Our model for dust production and destruction is based on that developed by V19, but with the following key improvements: the inclusion of dust survival in the hot gas surrounding galaxies and its gradual destruction via thermal sputtering, the radial resolution of dust within the ISM, improvements to the SN shock destruction modelling, and full synchronisation with the L-Galaxies HI-H 2 partitioning model.The new dust model also has important differences to those commonly used in other cosmological simulations, for example by allowing a wider range of accretion timescales during grain growth (Section 4).
With these models in place, we obtain the following results from two versions of L-Galaxies, one including binary stars (MM+binary_c) and one including only single stars (MM+singleStars): • The inclusion of binary stars into L-Galaxies does not have a significant impact on key galaxy properties such as stellar mass, gas mass, and star formation rate.This is because the amount of mass and energy ejected by massive stars remains relatively unchanged compared to the single-star-only approximation.The mass ejected by massive stars determines the strength of the SN feedback in L-Galaxies, which in turn drives gas ejection rates and consequently star formation and stellar mass growth (Section 5.1).
• MM+binary_c exhibits higher gas-phase oxygen abundances in low-mass galaxies than MM+singleStars.This is predominantly due to differences in the low-metallicity SN-II yields assumed, rather than binary effects.This additional oxygen is typically found in the outer disc of star-forming galaxies by  = 0 (see Section 5.1).
• Binary effects do impact the ejection of some lighter elements.
For example, CE ejection, novae, and SNe-Ia help boost nitrogen yields at very early or very late times at low metallicities.
• L-Galaxies matches observed dust properties back to  ∼ 4, with slightly better agreement for MM+binary_c.An evolution in the DTM and DTG ratios of metal-poor galaxies is also present, in qualitative agreement with that seen in GRB-DLA samples.These successes are largely due to the shorter accretion timescales in dustier environments predicted by our grain growth formalism, thus allowing more efficient dust production at late times (Section 5.2).
• At high , dust production and destruction rates are significantly higher in MM+binary_c than MM+singleStars.This is due to higher dust masses and shorter accretion timescales at early times.These rates then drop below those found in MM+singleStars at low , due to longer destruction timescales.In general, there is a coupled evolution between the main dust production and destruction mechanisms in L-Galaxies as the key timescales and dust reservoirs change over time (see Section 5.3).
• Grain growth overtakes SNe-II as the dominant dust production mechanism at  ∼ 6.4 in both versions of L-Galaxies.This suggests that the balance between the main sources of dust production and destruction in galaxies is relatively independent of the details of stellar evolution (Section 5.3).
• On sub-galactic scales, there is very good overall agreement between L-Galaxies and observed dust and metal radial profiles at  = 0.However, dust observations suggest slightly flatter Σ dust profiles in the centres of galaxies (i.e.≲ 1 e ), in contrast to both L-Galaxies and observed O/H profiles.This could suggest either enhanced dust destruction or suppressed dust production at low radii, perhaps associated with the presence of AGN or stellar bulge (Section 5.4).
• The DTG ratio in the diffuse ISM drops by ∼ 0.8 dex at large radii (i.e.≳ 3 e ) on average in L-Galaxies, particularly in low-mass (low-metallicity) systems and when binary stars are included.This demonstrates that the assumption of a constant DTM ratio does not hold on sub-galactic scales.The degree of this drop in DTM depends on the grain growth efficiency in molecular clouds and the amount of dust that is able to survive in the CGM and re-accrete onto galaxy discs at large radii (Section 5.4).
There are a number of avenues of further investigation made possible by this work.For example, the effect of modifying the standard parameters in binary_c away from their default values will be studied, in particular the effect of more complex prescriptions for the fate of massive stars, SN-Ia progenitor evolution, and higher-order multiple systems.The evolution of more exotic chemical species, such as s-and r-process elements, individual isotopes, and rare elements of interest to Milky Way studies, will also be presented in follow-up papers.The potential for enhanced gas cooling at higher  to reconcile a number of inconsistencies between L-Galaxies and observations will also be investigated, as well as additional dust production and destruction mechanisms and the dependence of H 2 formation on the dust mass.Finally, we highlight that the choice of stellar yields always has a significant impact on the chemical evolution of galaxies seen in simulations.In this work, we have chosen yield sets that are commonly used in the literature, to facilitate a more meaningful comparison between old and new versions of L-Galaxies, and to other simulations.However, the inclusion of more detailed chemical yields will be the main the focus of future work.
The authors would like to thank the referee for their helpful and insightful comments, as well as Kasper Heintz, Chiaki Kobayashi, Zara Osborn, Massimiliano Parente, and Patricia Schady for valuable discussions during the undertaking of this work.We also thank Abdurro'uf, Viviana Casasola, Bethan Easeman, and Kasper Heintz for sharing their observational data.APV acknowledges support from the Carlsberg Foundation (grant no CF20-0534).The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140.

DATA AVAILABILITY
The full source code for the L-Galaxies versions used in this work are publicly available via GitHub (see the L-Galaxies website for more details, including installation and running instructions).Associated L-Galaxies output catalogues can be directly obtained via Zenodo.The binary_c ensembles used (in JSON format), including the python script used to generate L-Galaxies yield tables from them, are also available via the L-Galaxies GitHub and Zenodo repositories.Additional simulation and observation data products derived by the authors and presented here can be obtained upon request.

Figure 1 .
Figure 1.The normalised metal ejection rate (  Z,norm,SP = d Z /d 1/ M ⊙ ), for a binary_c stellar population including single and binary stars with an initial metallicity of  = 0.004.The 16 different ejection processes modelled in binary_c are shown (see binary_c documentation for details).

Figure 2 .
Figure 2. Top panel:Total mass ejection rate from stellar winds (blue), SNe-II type enrichment channels (orange), and SNe-Ia + novae (green), assuming a constant SFR of  = 1 M ⊙ /yr and a fixed metallicity of  = 0.004.Dotted lines denote the single-star-only set-up used in previous versions of L-Galaxies, dashed lines denote the single-star-only set-up from binary_c, and solid lines denote single+binary stellar populations from binary_c.Bottom panel: Same as top panel, but for the metal mass ejection rate, considering C, N, O, Ne, Mg, Si, S, Ca, and Fe.

Figure 3 .
Figure 3. Mass ejection rates from a single 1 M ⊙ /yr burst of star formation for the 11 chemical elements considered in this work.Dotted lines denote the singlestar-only set-up used in previous versions of L-Galaxies, dashed lines denote the single-star-only set-up from binary_c, and solid lines denote single+binary stellar populations from binary_c.

Figure 4 .
Figure 4. First panel: Dust accretion timescale for grain growth in molecular clouds as a function of ISM oxygen abundance for individual radial rings.Second panel: Dust destruction timescale for SN shocks in the ISM versus ISM oxygen abundance for individual rings.Third panel: Dust destruction timescale for thermal sputtering in the CGM as a function of galaxy stellar mass.Fourth panel: Dust destruction timescale for thermal sputtering in the Ejecta phase as a function of galaxy stellar mass.In all cases, the MM+binary_c (blue) and MM+singleStars (red) versions of L-Galaxies are shown, run on Millennium and Millennium-II combined.Contours represent the 1-4  spread in the distributions.

Figure 6 .
Figure 6.Relation between SFR-weighted metallicity and the DTM ratio (top panel), DTG ratio (middle panel), and  dust (bottom panel).In all cases, the MM+binary_c (blue) and MM+singleStars (red) versions of L-Galaxies are shown, run on Millennium and Millennium-II combined, with metallicity measured as the SFR-weighted oxygen abundance.Orange points represent observational data from Rémy-Ruyer et al. (2015) (circles), DustPedia (De Vis et al. 2019; Casasola et al. 2022) (diamonds), and the collection of Péroux & Howk (2020) (squares).There is a good match between both L-Galaxies versions and observations at  = 0.

Figure 8 .
Figure 8. Top panel: Evolution of various cosmic dust production and destruction rates from  = 8 to 0 for the MM+binary_c (solid lines) and MM+singleStars (dashed lines) versions of L-Galaxies run on Millennium-II.Bottom panel: Evolution of the dust mass density in various galaxy and halo components for the same L-Galaxies versions.Observational data is taken from Dunne et al. (2011); Beeston et al. (2018); Driver et al. (2018); Pozzi et al. (2020) (orange circles) and Péroux & Howk (2020) (orange squares).In both panels, values at discreet output redshifts have been fit with 5th order polynomials.Dust evolution in the two versions of L-Galaxies differs significantly at high , predominantly due to differences in ISM metallicity.

Figure 9 .
Figure 9. Stacked radial profiles for star-forming galaxies at  = 0 from the MM+binary_c (blue) and MM+singleStars (red) versions of L-Galaxies run on Millennium, separated into three mass bins.First row: Gas-phase oxygen abundance (excluding oxygen depleted onto dust).Observational data is from Yates et al. (2021a) (circles) and Erroz-Ferrer et al. (2019); Easeman et al. (in prep.)(squares).Second row: Oxygen enhancement in disc + bulge stars normalised to the solar photospheric value of (  O / Fe ) ⊙ = 4.44 fromAsplund et al. (2009).Third row: The DTM ratio in the ISM for molecular clouds (solid lines) and diffuse gas (dashed lines).Fourth row: The DTG ratio.Fifth row: Dust mass surface density, with observational data of six galaxies from the DustPedia dataset ofCasasola et al. (2017).Both versions of L-Galaxies match observational data well, and predict a significant drop in DTM and DTG ratios at large radii.

Figure 10 .
Figure 10.Radial profiles from the HERACLES sample (Abdurro'uf, et al. 2022, points) and star-forming disc galaxies within the same mass range [9.7 ≤ log(  * / M ⊙ ) < 11.2] from the MM+binary_c version of L-Galaxies (lines) run on Millennium.MM+binary_c matches the observed profiles very well, albeit with a slight excess of SFR at low radii.