BEDE: Bayesian Estimates of Dust Evolution For Nearby Galaxies

We build a rigorous statistical framework to provide constraints on the chemical and dust evolution parameters for nearby late-type galaxies with a wide range of gas fractions ($3\%<f_g<94\%$). A Bayesian Monte Carlo Markov Chain framework provides statistical constraints on the parameters used in chemical evolution models. Nearly a million one-zone chemical and dust evolution models were compared to 340 galaxies. Relative probabilities were calculated from the $\chi^2$ between data and models, marginalised over the different time steps, galaxy masses and star formation histories. We applied this method to find `best fitting' model parameters related to metallicity, and subsequently fix these metal parameters to study the dust parameters. For the metal parameters, a degeneracy was found between the choice of initial mass function, supernova metal yield tables and outflow prescription. For the dust parameters, the uncertainties on the best fit values are often large except for the fraction of metals available for grain growth, which is well constrained. We find a number of degeneracies between the dust parameters, limiting our ability to discriminate between chemical models using observations only. For example, we show that the low dust content of low-metallicity galaxies can be resolved by either reducing the supernova dust yields and/or including photo-fragmentation. We also show that supernova dust dominates the dust mass for low metallicity galaxies and grain growth dominates for high metallicity galaxies. The transition occurs around $12+\log({\rm O/H})=7.75$, which is lower than found in most studies in the literature.


I N T RO D U C T I O N
Galaxies are complex systems consisting of stars, dust, heavy elements, and multiple gas phases.These complex systems evolve from relatively simple gas clouds under the influence of ongoing star formation (SF) and associated processes.Metals are expelled into the interstellar medium (ISM) when these stars end their lives either in the wind of asymptotic giant branch (AGB)-stars or when they explode as supernovae (SNe).Some of the metals that are expelled at the end of a star's life condense as dust grains, which are also mixed into the ISM.These dust grains absorb and scatter about one quarter of the stellar radiation (Bianchi et al. 2018 ) and re-emit most of the absorbed energy at f ar-infrared-submm w av elengths.The y thus strongly influence the way we observe galaxies.Dust grains also act as a catalyst for the formation of molecules, and metals continue to accrete on to the dust grains in the dense phases of the ISM (see Galliano, Galametz & Jones 2018 for a re vie w).In the local Universe, dust contains roughly one quarter (De Vis et al. 2019 ) to one half (Maiolino & Mannucci 2019 ) of the heavy elements in the ISM.It is thus important to study both dust and metals when trying to understand the chemical evolution of galaxies.
Chemical evolution models are a tool that can be used to interpret observed metal abundances and dust masses together with changes in other galaxy properties such as, e.g.gas and stellar masses and star formation histories (SFH).This is done by numerically solving first-order differential equations assuming some initial conditions, an initial mass function (IMF), stellar lifetimes, theoretical nucleosynthesis yields for various elements, a SF prescription and a prescription for the flow of gas in and out of the galaxy.Chemical evolution models focusing on the metal content of galaxies have been around for a long time, since the pioneering work of Schmidt ( 1959 ), and still form a rele v ant and acti v e field of study (e.g.Tinsle y 1980 ; Carigi, Hernandez & Gilmore 2002 ;Andrews et al. 2017 ;Zhang et al. 2018 ;Romano et al. 2019 ;Kang et al. 2021 ).Dust can also be included if the proper dust formation and destruction prescriptions are accounted for.Chemical evolution models can be divided into three categories.One-zone models study the integrated properties of galaxies without spatial resolution, and assume instant mixing of dust, gas, and metals (Dwek 1998 ;Lisenfeld & Ferrara 1998 ;Hirashita & Ferrara 2002 ;Inoue 2003 ;Morgan & Edmunds 2003 ;Valiante et al. 2009 ;Asano et al. 2013 ;R émy-Ruyer et al. 2014 ;Rowlands et al. 2014 ;Zhukovska 2014 ;Clark et al. 2015 ;Feldmann 2015 ;Popping, Somerville & Galametz 2017 ;De Vis et al. 2017b ;Zhukovska, Henning & Dobbs 2018 ).These models study the changing balance between the dust sources and sinks, and include processes such as the formation of dust in stellar winds and SNe, dust growth and destruction in the ISM, and radiation field effects.Most recent studies have found grain growth to be the dominant dust source for evolved galaxies (e.g.Dunne et al. 2011 ;Zhukovska 2014 ;De Vis et al. 2017b ;Zhukovska et al. 2018 ).
Multizone chemical evolution models are similar to one-zone models, b ut ha v e sev eral re gions in which the evolution is tracked separately.These models are often used to study radial gradients of the different properties (e.g.Boissier & Prantzos 1999 ;Spitoni et al. 2021 ).Finally, hydrodynamical simulations provide the most realistic framework and track the production of metals and dust and how these flow through galaxies (Bekki 2013(Bekki , 2015 ; ;Aoyama et al. 2016 ;McKinnon, Torrey & Vogelsberger 2016 ;McKinnon et al. 2018 ;Aoyama, Hirashita & Nagamine 2019 ).Ho we ver, these hydrodynamical simulations are much more computationally e xpensiv e, and it is thus not possible to explore the parameter space in the same way as for one-zone models.
In order to model the build-up of metals and dust within the galaxy, it is important to have a realistic prescription of the amount of gas flowing in and out of the galaxy.Pristine gas continuously accretes on to galaxies from the cosmic web as a result of gravity.At the same time, gas is blown out of the galaxy driven by the energy released from stars and SNe (Che v alier & Clegg 1985 ), as well as from supermassive black holes/active galactic nuclei (AGNs; Begelman, de Kool & Sikora 1991 ).Not all of the gas in these outflows will have sufficient energy to leave the gravitational potential well of the galaxy, and will thus eventually fall back on to the galaxy (e.g.Nelson et al. 2019 ).This recycling of the outflows is often called the galactic fountain.
There have been numerous studies that have used chemical evolution models including inflows and outflows (e.g.Calura, Pipino & Matteucci 2008 ;Feldmann 2015 ;De Vis et al. 2017b ;Zhukovska et al. 2018 ) as well as hydrodynamical simulations that intrinsically include prescriptions for the flow of gas throughout the galaxy (e.g.McKinnon et al. 2018 ;Aoyama et al. 2019 ).More recently, Mill án-Irigoyen, Moll á & Ascasibar ( 2020 ) included a framework to follow the evolution of dust in the atomic, ionized, and molecular gas phases.Typically, in all these studies, there are multiple parameters that can have a similar effect on the observed build-up of metals (e.g.changing the IMF or the SN metal yield tables can both affect the metallicity in the models in a similar way).To reveal the de generac y between such parameters, each combination of parameters has to be explored and e v aluated.Additionally, without a statistical frame work, e ven though it is possible to find a model that fits the observations well, it is impossible to know whether this is the only viable option.Without such a framework, one cannot put robust constraints on the input parameters of chemical models.On the other hand, using a formal statistical (e.g.Bayesian) framework, one can determine the most likely values for each parameter and even uncertainties (in case of numerical continuous variables).Bayesian frameworks have been used in previous chemical models to follow the evolution of elements (e.g.Rybizki, Just & Rix 2017 ;Belfiore et al. 2019 ;Spitoni et al. 2021 and references therein).
De Looze et al. ( 2020 ) used a Bayesian framework to put statistical constraints on a number of dust and chemical evolution parameters.This work provided an impro v ement on previous works by accounting for a wider range of model parameters yet did not sample the entire parameter space sufficiently.F or e xample, it can be seen from their results that the build-up of metals with decreasing gas fraction is not modelled well.Additionally, they combine their galaxy observations in six gas-fraction bins before fitting them and focus most of the discussion on two of these galaxy bins only.In this work, we use an impro v ed statistical framework that allows us to fit all individual galaxies simultaneously to a combined set of models.Additionally, we show that by sampling a wider parameter space, a good match to the observations can be found for all latetype galaxies.Including realistic models for galaxies spanning a wider range in gas fraction provides stronger constraints as it is often the most unevolved and evolved galaxies (i.e. at the extremes) that provide the most strain on the models.
In this work, we impro v e on the model of De Vis et al. ( 2017b ) where we found that low metallicity dwarfs required different dust production parameters to late-type spirals (see also Zhukovska 2014 ;Feldmann 2015 ).We include physically moti v ated dust models from laboratory data and more physical outflows to follow dust destruction, formation, and recycling in and out of clouds and it tracks both the interstellar and intergalactic medium (IGM).We use a grid of ∼40 000 models and a Monte Carlo Markov Chain (MCMC) framework to constrain the chemical evolution parameters related to the build-up of metals, and to rev eal an y de generac y between these parameters.Subsequently, we use the same approach to constrain the model parameters related to dust, using the metal-related model constraints from this work as input and comparing models with different dust parameters.In Section 2, we present the observational samples that are used to constrain our models.Section 3 details our chemical evolution model and the grid of models used.Our statistical framework is explained in Section 4. We illlustrate some of the main parameter dependencies of our models in Section 5 and our results from the statistical framework are presented in Section 6.Finally, we list some caveats in Section 7 and our conclusions are given in Section 8.

DustPedia
To best constrain our chemical evolution models, it is key to have a sizable sample of galaxies for which we have reliable gas, stellar, and dust masses, as well as star formation rates (SFRs), metallicities, and Nitrogen-to-Oxygen ratios.Few samples have all these data available, and the largest that does is the DustPedia (Davies et al. 2017 ) sample.DustPedia is a collaborative focused research project working towards a definitive understanding of dust in the local Universe, by capitalizing on the le gac y of Herschel .The full DustPedia sample consists of 875 nearby ( v < 3000 km s −1 ), extended ( D 25 > 1 arcmin) galaxies that have been observed by Herschel and have a near-infrared detected stellar component.These galaxies hav e e xcellent multiwav elength aperture-matched photometry available (typically 25 bands; Clark et al. 2018 ).The spectral energy distributions (SED) are fitted using CIGALE (Noll et al. 2009 ) and the resulting galaxy properties are presented in Nersesian et al. ( 2019 ).DustPedia uses the physically moti v ated (based on laboratory data) THEMIS dust model (Jones 2013 ;K öhler, Jones & Ysard 2014 ;K öhler, Ysard & Jones 2015 ;Ysard et al. 2015 ) as reference dust model, which was incorporated into CIGALE .
The gas masses, metallicities, and Nitrogen-to-Oxygen ratios for DustPedia are taken from De Vis et al. ( 2019 ).The authors compiled H I fluxes from the literature and combined them with H2 measurements from Casasola et al. ( 2020 ).The total gas masses used include both H I and H 2 (either measured or estimated from Casasola et al. 2020 ) as well as elements heavier than Hydrogen. 1For the metallicities, a literature compilation of emission line fluxes was done and combined with MUSE spectrophotometry.Characteristic metallicities were determined for each galaxy by fitting radial profiles to the available H II regions using a Bayesian framework.Various metallicity calibrations are available, and the Pilyugin & Grebel ( 2016 , hereafter PG16 ) calibration was chosen as the reference calibration in this work, following De Vis et al. ( 2019 ).We also give results for the IZI calibration in Appendix A for comparison.Nitrogen-to-Oxygen ratios are also available from PG16 .In total, 317 DustPedia sources have all the necessary observations and uncertainties (in the literature compilations, not all sources had uncertainties).Galaxies without uncertainties are discarded as they cannot be used as reliably in our statistical framework.In line with De Vis et al. ( 2019 ), the DustPedia sample has been divided into late-type galaxies (LTG) and early-type galaxies (ETGs).Only LTG galaxies are used in the statistical framework in Section 6 as our chemical evolution model is not representative for ETGs.Both LTGs and ETGs are shown throughout the plots in this work for completeness.

H I and dust-selected samples
Since DustPedia requires a 5 σ detection in the WISE W1 band and a diameter ( D 25 > 1 arcmin) as the selection criteria, it is somewhat biased against dwarf galaxies.In De Vis et al. ( 2017a , b ), we have shown the importance of including unevolved dwarf galaxies when studying dust and gas scaling relations.Given that our aim is to study dust in chemical evolution models, it is crucial that we have observations that span as wide a range of evolutionary states as possible.Therefore, we add the HAPLESS (Clark et al. 2015 ) and H I GH (De Vis et al. 2017a ) samples to increase our sample size at the high gas-fraction end.HAPLESS is a blind dust-selected, volume-limited sample of 42 local ( z < 0.01) galaxies detected at 250 μm from the H-ATLAS Phase 1 Version-3 internal data release, co v ering 160 de g 2 of the s (Bourne et al. 2016 ;Valiante et al. 2016 ).The H I -selected H I GH sample is extracted from the same H-ATLAS area and includes 40 unconfused H I sources identified in the H I Parkes All Sky Survey (Barnes & Hernquist 1992 ;Meyer et al. 2004 ) and the Arecibo Le gac y F ast ALFA Surv e y (Gio vanelli et al. 2005 ;Haynes et al. 2011 ;Haynes et al., pri v ate communication); 24 of these sources o v erlap with the HAPLESS sample.De Vis et al. ( 2017a ) compiled far -ultra violet-submm photometry for each of these samples, and subsequently derived dust masses, stellar masses, and SFRs consistently using MAGPHYS (da Cunha, Charlot & Elbaz 2008 ).For consistency with the THEMIS dust masses for DustPedia, the H I GH and HAPLESS dust masses are scaled down by a factor of 1.075 (De Vis et al. 2019 ).PG16 Metallicities for H I GH and HAPLESS are taken from De Vis et al. ( 2017b ).Nitrogen-to-Oxygen ratios were calculated for this work following exactly the method used for DustPedia (De Vis et al. 2019 ).We hav e remo v ed all the H I GH and HAPLESS sources that are already present in the DustPedia sample and all sources that do not have all the required data available.The remaining 18 H I GH sources and 5 HAPLESS sources are all LTGs and are added to our observational sample.

Build-up of metals
We use a chemical evolution model to build a consistent picture of how the metal, stellar, and gas content change as galaxies evolve.We use a one-zone model where only the integrated properties of the galaxies are modelled.We do separate the ISM into clouds and the diffuse ISM.Within their phases, the gas, dust, stars, and metals are assumed to be perfectly mixed.This model is directly based on that of De Vis et al. ( 2017b ), though with considerable changes (especially to the outflow and SF prescriptions).All models start as pristine clouds of gas, which are converted into stars as a result of the ongoing SF and a given IMF.The stars are then tracked throughout their lifetimes and when they end their life either as AGB stars or SN, the y e xpel dust and metals into the ISM.Inflows and outflows also alter the ISM content of the galaxies in our model.
To determine the total stellar, metal, dust, and gas content of the model galaxy, it is necessary to integrate over time, as well as o v er stellar mass.We use a numerical integration for this, with discrete stellar mass and time-steps.For the integration over stellar masses, 500 steps are used, logarithmically spaced between the upper and lower limit of the IMF (see equation 1).The size of the time-steps is set to 30 million years.This value is chosen to correspond to the time for the dissociation of a molecular cloud (McKee 1989 ), which simplifies the treatment of the cloud dissociation, as will be discussed in Section 3.4.Throughout this work, we show continuous integrals in the equations, corresponding to the theoretical dependencies.In practice, these are all implemented using a numerical integration o v er the discrete steps (i.e.d t = 30 Myr ).
The evolution of the gas mass can be described by where M g is the gas mass, ψ( t ) is the SFR, φ( m ) is the stellar IMF (normalized so that , and m R is the remnant mass of a star of mass m (Ferreras & Silk 2000 ).m U is the upper mass limit of the stars (which is set to 120 M ), m L is the lower mass limit of the stars (which is set to 0.8 M ), and m t is the lowest mass for which a star could have reached the end of its life by time t .The lifetime τ m of stars with initial mass m is derived using the model in Schaller et al. ( 1992 ).The first term in equation (1) accounts for gas consumed during SF and the second term for how much gas is returned at time t by stars of all masses combined.The third and fourth term, I ( t ) and O ( t ) are simple parametrizations of the rate at which gas is contributed or remo v ed via pristine inflows and outflows, respectively .Finally , the fifth term R ( t ) gives the rate at which the outflowing gas is recycled (see Section 3.2).Similarly, the evolution of the mass of metals ( M Z ) is given by Here, Z is the metal mass fraction defined as Z = M Z / M g .The first term accounts for metals expelled by stars and SNe.This includes metals re-released by stars after they die, and newly synthesized metals ejected in winds and SNe.mp Z are the metal yields that are taken separately from SN or AGB yield tables taken from the literature.We will explore a range of different yield tables in this work.The inflows in this work have a pristine metallicity of Z I = 0 (Coc et al. 2012 ) and the outflows use the current metallicity of the galaxy.The last term in equation (2) again corresponds to the recycled outflows and will be discussed in the next section.The factor 1.36 is to account for elements other than Hydrogen in the gas and 16 and 14 are the atomic weights of Oxygen and Nitrogen, respectively.Before calculating the 12 + log (O/H) and log (N/O), the oxygen masses are corrected for the amount of oxygen locked up in dust.Following the THEMIS model, we assume the average oxygen content of dust by mass is 23.8 per cent (Jones et al., in preparation).

Inflows and outflows
Galaxies continuously accrete gas from the surrounding IGM.This inflow of pristine gas is particularly strong at early evolutionary stages.We use the prescription of Zhukovska ( 2014 ): where τ inf is the infall time-scale, M inf is the amount of gas falling into the galaxy, and t G is the total amount of time o v er which this gas is accreted.The infall time-scale is set to τ inf = 2 Gyr; we hav e e xperimented with changing this value, and found v ery little difference to the results after a few Gyr, as long as τ inf τ G .M inf is set to half of the total mass of the galaxy M tot (which is a free parameter).This means the galaxy will start out as a gas cloud with mass of 0.5 M tot and the same amount of gas will be accreted by inflows.We have also explored models where M tot was left the same, but was di vided dif ferently between the primordial cloud and the inflowing material, without much change to the models after a few Gyr.The exact prescription for the inflows makes little difference as long as the majority happens at early evolutionary stages.We can even start with just a cloud of pristine gas, with little effect on the chemical evolution of a galaxy.Throughout the rest of this work, we will refer to the total galaxy mass M tot as the sum of the mass of the pristine cloud and the total infalling material (each set to 50 per cent of M tot ).
Galactic winds driven by the ongoing SF have a more significant ef fect.These outflo ws dri ve metals and gas from the ISM and slow down the build-up of metallicity.We express the rate of outflowing gas relative to the rate of SF as the mass loading factor: Due to their shallower gravitational potential wells, lower mass galaxies have higher mass loading factors.Yet at the same time, AGNs in massive galaxies efficiently blow out mass from the galaxy and have high mass loading factors too.We base our outflow prescription on mass loading factors taken from the Illustris TNG50 simulation (Marinacci et al. 2018 ;Naiman et al. 2018 ;Pillepich et al. 2018 ;Springel et al. 2018 ).For each halo in the simulation, they measure the outflow rate passing various radii with a range of velocities, as a function of redshift.We were provided with a set of tables giving the median mass loading factor for a radius of 10 kpc, as a function of outflow velocity, stellar mass ( M * ), and redshift (Nelson, pri v ate communicaon).Then, we use a simple bi-linear interpolation to estimate the outflow for any given M * and redshift,2 and sum over the velocity components.We also impose a strong limit so that no more than 50 per cent of the gas mass can be blown out of the galaxy in a single (30 Myr) time-step.
To determine what happens to the outflows after they leave the galaxy, we use the outflow velocity distribution to determine the median velocity of each component in each bin.The mean velocity for each bin is then used together with the average total (baryon + dark matter) halo mass profile for each bin from the Illustris TNG100 simulation (extracted using the ILLUSTRIS PYTHON package; Nelson et al. 2019 ) to determine the (friction-less) ballistic trajectory for each outflow component using a simple simulation.If the mean outflow velocity is larger than the escape velocity, the outflowing gas is lost into the IGM.Ho we ver, a large fraction of the outflows fall back on to the galaxy after a time τ rec determined by its ballistic trajectory.We refer to this returned gas as recycled gas and to the combined processes of outflow and recycling as the 'galactic fountain'.Before the gas is recycled, a fraction is lost to the IGM.The fraction of gas lost scales with the time the outflow spends in the IGM (i.e. the time before it is recycled) as where t i is the time at which the outflows are expelled and t is the time for which we are calculating the infall rate of recycled gas. is a scaling factor for how much of the gas that would otherwise fall back on to the galaxy is lost due to interactions with the IGM.In our work, is set to 0.2 Gyr −1 , corresponding to a loss of 18 per cent of the outflow mass per Gyr spent in the IGM.ρ( t − t i ) gives the fraction of the outflows that are falling back on to the galaxy after a time t − t i , and is determined by the velocity distribution of the outflows, and how long it takes for them to be recycled.Since we have discretized our outflows into three components with a single outflow velocity for each (we use the mean velocity of that component as described abo v e), ρ will be equal to 1 when t − t i = τ rec ,v out and 0 otherwise.
Here, τ rec ,v out is the time it takes the outflow to fall back on to the galaxy given the outflow velocity v out and the galaxy's mass profile.Equation ( 7) becomes Here, we have introduced the scaling factor f recy as an additional free parameter in our models.f recy scales the recycling time τ rec up or down to account for the uncertainty in our calculation of τ rec .
The sum in equation ( 8) is o v er the different components v out for which ρ = 1.Usually, this means the three outflow components are summed, but occasionally multiple components from different redshift or mass bins contribute to the current recycling rate.We note that if the outflow velocity is larger than the escape velocity of the halo, τ rec ,v out = ∞ and the contribution to R ( t ) will be zero.This gas is automatically lost to the IGM.The metal and dust content of the outflows is determined by the g alaxy's metal-to-g as ( Z ) and dust-to-g as ( δ) ratios at the time the gas is blown out ( t − f recy τ rec ,v out ).For the outflows with v out > 150 km s −1 , it is assumed that all dust is destroyed by shocks (Jones, Tielens & Hollenbach 1996 ) and transformed to gas-phase metals.By the time the outflows are recycled, the galaxy will have typically evolved to higher Z and dust-to-gas ratios, which are then diluted by the recycled gas.The galactic fountain thus slows down the build-up of metals and dust.

Star formation rates
In De Vis et al. ( 2017b ), we used a few template SFHs (one exponentially declining, one delayed, and one bursty SFH) for our models.The delayed SFH was found to be the best prescription for most sources.Ho we ver, for that study, the mass loading factors of the outflows were smaller and not mass dependent.When the mass loading factors are mass dependent, as is the case in our work, it is not possible anymore to define one SFH for all models.Galaxies with higher mass loading factors will run out of gas at different times, and it is not possible to account for this using a single SFH.Instead, a better approach is to define our SFR based on a gi ven SF ef ficiency ( SFE = ψ/ M g ).When using a fixed SFE, if a galaxy runs out of gas (due to the combination of ongoing SF and outflows), the SFR will automatically decrease, following the reduction in available gas mass.As a result, the outflow rate will also decrease (mass loading factor stays the same).The resulting SFH will have the shape of an exponentially declining SFH.
Although a constant SFE is an impro v ement compared to using one single SFH, it is still not ideal as not all galaxies have the same SFE.Indeed, SFE is likely correlated to the stellar mass of a galaxy as the higher mass surface density produces a higher hydrostatic pressure in the ISM (Elmegreen 1989 ;Wong & Blitz 2002 ), it may also be redshift dependent since a higher turbulence at a fixed stellar mass would result in lo wer ef ficiencies (Hayward & Hopkins 2017 ).Finally, it also well known that very low gas fraction sources are usually quenched.
The SF in this work is made up of (i) a continuous component with a slowly varying SF efficiency (SFE = ψ/ M g ) and (ii) bursts superimposed on the continuous component.Accounting for the dependencies described abo v e, we create an SFE prescription that varies with mass, redshift, and gas fraction, and empirically produces SFRs that match the observations (multiple power-law exponents were trialled).The prescription is given by where M * is the stellar mass, z is the redshift, 3 and SFE 0 is the reference SFE, which is a free parameter.Three options of SFE 0 are explored in this w ork: SFE 0 = 10 −8 . 5yr −1 (f ast), SFE 0 = 10 −9 yr −1 (average), and SFE 0 = 10 −9 . 5yr −1 (slow).We note that this prescription is not physically moti v ated, but does provide a sensible frame work, where lo w g as fraction g alaxies are quenched, yet immature low-mass galaxies also have rather low SFE.The resulting SFH that follows from equation ( 9) also match the delayed SFHs observed in various works (Fig. 1 ), and are naturally consistent with any ongoing changes to the gas mass.
For the bursty SFH, we follow Bruzual & Charlot ( 2003 ) and generate random bursts so that there is a 50 per cent probability that a galaxy had a burst within any 2 Gyr period of its evolution.The duration of the burst is randomly drawn from a uniform distribution between 30 and 300 Myr.The SFR during the burst is set so that the stellar mass formed during the burst is a given fraction of the stellar mass at the time.The amount of stars formed for each individual burst is drawn from a loguniform distribution between 0.004 and 0.4 times the stellar mass of the galaxy at the start of the burst.Some examples of the SFH used in this work are shown in Fig. 1 .After a burst, there is often a period of reduced SFR compared to the continuous SFH.This results from gas being used up (both by 3 The same redshifts were used as for the outflows (see previous footnote).).The colour of the model indicates its total galaxy mass, and the line type the reference SFE ( SFE 0 = 10 −8 . 5yr −1 : solid, SFE 0 = 10 −9 yr −1 : dashed and SFE 0 = 10 −9 . 5yr −1 : dotted).SFHs including b ursts ha ve been included using lighter shaded lines.The SFHs have the shape of a delayed SFH, which peaks at different times for different mass galaxies.
SF and outflows) and since the SFE only varies slowly in our model, the drop in gas mass after the burst results in a drop in SFR.The reduced SFR (and associated outflows), combined with recycling of the outflows expelled during the burst, means that the SFR in a bursty SFH converges back to where it would be for the continuous SFH.

Dust parameters
Interstellar dust forms in a range of environments, such as the winds of e volved lo w-to-intermediate mass stars (LIMS; Sargent et al. 2010 ), core-collapse SNe ejecta (Dunne et al. 2003 ;Gomez et al. 2012 ;De Looze et al. 2017 ) and grain growth and accretion in the ISM (Mattsson & Andersen 2012 ).Some of the produced dust and metals will stay in the molecular clouds they were formed in, and some will dissipate into the diffuse ISM.As the molecular clouds collapse to form the next generation of stars, the newly formed dust will be consumed together with the gas as fuel for the stars.SNe shocks also destroy dust as high-energy ions 'sputter' atoms from the surface of dust grains, and collisions between dust grains also break them up (Jones et al. 1996 ;Jones & Nuth 2011 ).Additional processes such as thermal sputtering, and ionizing destruction by cosmic rays, high-energy photons, and free electrons further reduce the dust mass (see Jones 2004 , for a re vie w).
To model the build-up and decline of dust, we include dust formation by stars (both LIMS and SN), dust grain growth in the diffuse and dense environments of the galaxy as well as dust destruction by SN shocks and photofragmentation of large grains.In this work, we essentially only track the evolution of large grains, as these make up the vast majority of the total dust mass.Small grains are easily destroyed (Bocchio, Jones & Slavin 2014 ) and thus will not be able to build-up to a significant amount of dust mass, even though they are essential to explain the MIR range of the SED.
The different dust processing mechanisms typically affect either only the diffuse ISM or only the cloud ISM.We therefore consider these two phases separately.Note that any reported dust-to-gas ratios will be the total dust-to-gas ratio, i.e. total dust mass divided by total gas mass, unless clearly indicated otherwise.f c gives the mass fraction of the ISM that is in cold dense molecular clouds.This fraction is kept constant throughout the ev olution, b ut as clouds are constantly dissociated and reformed on time-scales of ∼30 million years (e.g.McKee 1989 ), we reinitialize the ISM phases after every 30 million year time-step (the total gas and dust mass is redivided between clouds and diffuse ISM according to f c ).We do this by first mixing the cloud ISM into the diffuse ISM.During this cloud dissociation, the newly accreted dust mantled due to cloud grain growth will be exposed to a much harsher environment.As a result, the newly formed dust mass is reduced by 90 per cent to account for the e v aporation of ices in the dust mantles and to account for the processing of a-C:H mantle material into refractory a-C dust (Jones et al., in preparation).After the dissociation of these clouds, new clouds are formed using the updated dust-to-gas ratios from the diffuse ISM.The evolution of the total large grain dust mass is given by This expression is similar to equation ( 2), with the dust-to-gas ratio ( δ = M d / M g ) in place of Z .Noticeably, the astration, inflow, outflow, and recycling terms are essentially the same, where δ I is the dust-to-gas ratio of the inflowing material (set to zero) and R M d ( t) is the dust mass recycled by outflows.At the end of Section 3.2, we saw that only the low-velocity outflow component ( v out < 150 km s −1 ) contained any dust.Any time dust is blown out of the galaxy in this low-velocity component, we use the ballistic trajectory to determine if and when this dust will re-enter the galaxy as R M d ( t).
The first term in equation ( 10) gives the dust that is expelled by stars, which is determined from the amount of metals that is expelled, multiplied by a mass-dependent condensation efficienc y y d ( m ).F or LIMS, this efficiency is fixed at a value of 0.15 and will not be varied throughout this work. 4For SN, we use y d ( m ) based on the results of Todini & Ferrara ( 2001 ), and add a free parameter to scale the amount of dust produced up or down.
The dust grain growth, destruction by SN shocks and fragmentation rates apply to only one of the two ISM phases (as indicated by f c or 1 − f c ), and the various τ give the relevant formation and destruction time-scales for which the prescriptions are given below.As can be seen from the (1 − f c ) terms in equation ( 10), our dust destruction by SN shocks and fragmentation only remo v e dust in the diffuse ISM, as the SN shock velocities get too low to destroy dust in molecular clouds, and self-shielding in clouds makes the photofragmentation inefficient.f dis is a factor to account for how much cloud grain growth dust survives the dissociation of the clouds.f Si is the fraction of the dust that is made up of silicate cores, which are too robust to be affected by photofragmentation.
For the grain growth and destruction terms, we base our prescriptions on the THEMIS dust model (Jones 2013 ;K öhler et al. 2015 ;Ysard et al. 2015 ;Jones et al. 2017 ), which was also used in the determination of our observed dust masses (Nersesian et al. 2019 ). 4 After experimentation we found that changing the LIMS condensation efficiency has barely any affect on the final models.We thus do not include it as one of our free parameters in this work.
The THEMIS model consists of a mix of different grains with carbon and silicate core-mantle structures, for which the properties have been determined from laboratory-measured properties of physically reasonable interstellar dust analogue materials.THEMIS includes carbonaceous and silicate dust and grains of a range of different sizes.In this work we focus on the large grains, which make up the bulk of the mass.The ability of dust grains to evolve in response to the local physical conditions is one of the key concepts in THEMIS and different environments process grains in different ways.The cores of large grains typically consist of silicates or amorphous carbon (a-C).These cores usually accrete or form photoprocessed mantles of amorphous carbon (a-C).In transition to denser cloud environments, secondary a-C:H mantles form (a-C:H stands for hydrogenated amorphous carbon, i.e. the mantles were formed with more hydrogen in their molecular structure).The mantle formation increases both the mass of the individual grains as well as the total dust mass.In these dense clouds, the grains also begin to coagulate into aggregates (the mean grain size increases, but not the total dust mass).Then, in the the densest cloud environments, ice mantles will form on the aggregate grains (further increasing in the dust mass).We include grain growth terms for both the diffuse ( n H ∼ 10 2 cm −3 ) and dense ISM ( n H > 10 4 cm −3 ).Although we expect that the latter will dominate o v er the former due to the dependence on the accretion time-scales with n H , we include both phases here in order to not a priori make any assumptions.
The following prescriptions for the grain growth and destruction time-scales were taken from Jones et al. (in preparation).We refer to this work for further details on how these formulae were derived.Note that we express our prescriptions in terms of the associated timescales.Each of these enter equation (10 The grain growth time-scales are given by where Z 0 = Z / Z MW , i.e. the current metallicity relative to the Milky Way metallicity ( Z MW = 0.0134), k gg,dif is the diffuse grain growth scaling factor in Gyr −1 , and k gg,cloud is the dimensionless cloud grain growth scaling factor.The (1 − M d /( M Z × f )) factor accounts for the the depletion of the rele v ant elements.f dif and f cloud give the fraction of metals that are available for accretion in the diffuse ISM and molecular clouds.In the dense environments of clouds, the formation of ices and strong accretion of oxygen and carbon becomes possible and f cloud is 2.45 times higher than f dif .The SFR is included in the prescription for τ gg,cloud as stars form in dense ISM regions and this is where the bulk of the available matter that can form dust will accrete into mantles.When the molecular clouds are dissociated, the ice mantles e v aporate and only ≈10 per cent of the cloud dust mass accreted in dense clouds survives as a refractory material in the transition into the diffuse ISM (Jones & Ysard 2019 ).We set f dis = 0.1 (equation 10) to account for this.
Our prescription for the dust destruction by SN shocks is very similar to that from De Vis et al. ( 2017b ).Ho we ver, to be consistent with the THEMIS framework, it is now expressed in terms of the mass of dust destroyed per SN, M destr .The destruction time-scale is then given by where R SN is the SN rate.We convert the dust destruction rates for SN shocks from  7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5 Pre-existing gas mass and inflow mass combined SFH average.sfe,fast.sfe,slow.sfeSFE 0 = 10 −9 yr −1 , SFE 0 = 10 −8 . 5yr −1 , SFE 0 = 10 −9 . 5yr  ( 1989 ), where the dust mass shocked by an SN is proportional to the SN energy, the dust-to-gas ratio, and the velocity of the SN shock.
Assuming the THEMIS gas-to-dust ratio of 135 and shock speeds of 100 ( 200) km s −1 , McKee ( 1989 ) implies 27 (9) M SN −1 of dust is destroyed. 5In this work, we treat M destr as a free parameter, with range of values guided by Bocchio et al. ( 2014 ) and McKee ( 1989 ). Finally, we also include the photofragmentation of large a-C:H/a-C grains.Silicate grains are too robust to be affected by photofragmentation, though the Carbon mantles around them are.Photofragmentation is included as a destruction term as this mass is remo v ed from the large grains.The a-C nano-particle dust grains that are formed in this process will be rapidly destroyed and will never amount to a significant fraction of the dust mass (though they often account for a significant fraction of the dust luminosity).The photofragmentation time-scale is given by Here, we have assumed the diffuse UV radiation field G 0 is proportional to the specific star formation rate.This assumption can be made since a higher SFR per unit stellar mass would result in a larger contribution of UV photons from (hot) young stars, which in turn determines G 0 .The corresponding scaling factor between these has been folded in when going from k frag to k frag .We set f c = 0.5, f dis = 0.1, and f Si = 0.1, consistent with average ratios in the THEMIS model.

Grid of models
In order to better understand the parameter space, we build grids of models and compare to the observed properties of our galaxy samples.We vary the key parameters that determine the chemical and dust evolution of our models.The values and symbols we have used are listed in T able 1 .W e have made models with the following: 5 Assuming the relative amount of carbon and silicon is 1:2.
(i) Ten different total galaxy masses (sum of pre-existing initial gas mass and integrated mass of pristine inflows).These masses are logarithmically spaced between M tot = 10 7 M and M tot = 10 11 . 5M .
(iii) F our dif ferent IMFs φ( m ): We run models with the Chabrier ( 2003 ) IMF, the Salpeter ( 1955 ) IMF, the 'Galactic-field' (Kroupa & Weidner 2003 ) IMF [as opposed to the standard Kroupa ( 2001 ) IMF, which is more similar to the Chabrier ( 2003) IMF] and a top-heavy Chabrier IMF described by where m is the stellar mass in solar masses.
(iv) Fi v e metal yield tables for SNe ( mp Z, SN ): We use various tables from the literature to implement the amount of metals expelled by SNe.Our first three tables are the yield tables from Limongi & Chieffi ( 2018 , hereafter LC18 ) using their recommended set of yields (which is based on the mixing and fallback scheme and produces black holes for m > 25 M ) and rotation velocities of v rot = 0 km s −1 , v rot = 150 km s −1 , and v rot = 300 km s −1 , respectively.We have also implemented the yields from Maeder ( 1992 ) and Meynet & Maeder ( 2002 ).
(vi) Fi v e scaling f actors for the outflow recycling time-scale ( f recy ).
(vii) Four SN dust yields tables y d ( m ), which are taken from Todini & Ferrara ( 2001 ) and scaled down by , respectively , a factor SN red of 1, 5, 20, and 80.
(viii) F our v alues for the amount of dust mass destroyed per SN for a MW-like galaxy M destr .
(ix) Fi v e values for the photofragmentation parameter k frag .(x) Fi v e values for the cloud dust grain growth scaling factor k gg,cloud .
(xi) Three values for the diffuse dust grain growth scaling factor k gg,dif .
(xii) Three values for the fraction of metals that are available for grain growth in the diffuse ISM f dif .From these values, we have also determined the values for the fraction of metals that are available for grain growth in clouds f cloud as f cloud = f dif × 2.45 following the standard ISM THEMIS dust prescription (Jones et al., in preparation).
There are a large number of free parameters in our chemical evolution model.If we were to vary all of these parameters in one large grid, it would be necessary to run more than 97 million models.Unfortunately, this is not computationally feasible, even with the relatively simple and computationally light chemical evolution model we emplo y.Therefore, we tak e another approach and split the free parameters in two groups.The first group affects metal and SFrelated parameters.The second group of free parameters affects only the dust content of the simulated galaxy and can thus be decoupled from the metal and SF-related parameters.
We thus first vary only the free parameters that affect the metal properties (first six bullet points abo v e) and ignore dust for now, which results in a grid of 36 000 models (Table 1 ).We put statistical constrains on which models are most likely by comparing them to the observed properties (excluding dust) of our nearby galaxy samples.The results are presented in Sections 5.1 and 6.1.Next, we use the best-fitting metal parameters that come out of this analysis, and vary the model parameters that affect the galaxy dust mass, while keeping the metal parameters constant.The total galaxy masses and SFH are again varied here, as even with a fixed prescription for building metals, galaxies of different masses and at different evolutionary states are still necessary.The dust parameter grid consists of 162 000 models and is presented in Table 1 .This grid of dust parameters is illustrative (see Section 5.2 for the parameter dependencies).For the statistical results in Section 6.2, we use a direct (continuous) parameter search for the dust parameters.

S TAT I S T I C A L F R A M E W O R K
We have developed a statistical framework to compare our chemical evolution models described in the previous section to the observed galaxy properties in Section 2. In order to get constraints on the model parameters, we use a Bayesian MCMC approach.Here, we start with the metal grid; we simultaneously compare the observed gas mass, stellar mass, SFR, 12 + log(O/H), and log(N/O) to the model predictions for those parameters for the chemical evolution models in the metal grid described in the previous section.For the dust grid, we add dust masses to the abo v e observations and compare to models where the dust parameters are varied continuously, as described at the end of this section.
The likelihood function for our models is obtained by combining the likelihood of the individual galaxies: where x obs are the combined observations of all our galaxies and x obs ,i is a vector with all observed properties for each galaxy i .The individual likelihoods are not trivial to determine however for multiple reasons.In most cases, these kind of likelihoods are determined to compare an observed and a predicted value of an observ able.Ho we ver, in this case we compare an evolutionary track to a single observation.It is a priori not clear which point on the evolutionary track the observed values should be compared to.In order to get around this, we marginalize the likelihood o v er the different time-steps in the model in order to obtain the total likelihood of at any time.
In doing this, we account for the fact that the individual observed galaxies did not form at the same time and thus can be at a different point throughout their evolution.In addition, we also know that the individual galaxies did not all evolve from precursors of the same total galaxy mass and that they could have had quite different SFHs.
To account for this, we marginalize o v er all the M tot and SFH options in the grid and determine the total likelihood for any combination of IMF, mp Z ,SN , mp Z ,AGB , and f recy .In other words, we study 'families' of models (represented by the vector Models ) with same chemical evolution parameters, yet where the total galaxy mass, SFH, and formation time are allowed to vary.
Here, it has to be kept in mind that we are simultaneously fitting all observables [gas mass, stellar mass, SFR, 12 + log(O/H), log(N/O) and dust mass].The likelihood for a given time-step in one of the models is thus where x i, obs are the observables for each observation i , and x model the observables for a given model (at a given time t).σ x i and σ model are the uncertainties on the observations and the model uncertainty.Here, the model uncertainty is included to account for systematic errors in model conversion factors as well as the limited sampling we can do due to the discrete nature of our models (especially for bursty models, where the bursts themselves have a very limited sampling).For the stellar mass, gas mass, dust mass, and log(N/O) we use σ model = 0 .3 dex , which roughly corresponds to the intrinsic scatter in the observ ed sample.F or 12 + log(O/H), the observ ed scatter is smaller and we use σ model = 0 .15 dex .For the SFR, we use a larger amount of scatter as there is both larger intrinsic scatter in the observations and there is an inherent variation due to the b ursts.The interb urst periods and burst periods are randomly generated, so if the same model is generated multiple times there will be variation in the resulting SFR.Therefore, we use σ model = 0 .9 dex for the SFR.The final likelihood for comparing all available observations to a given family of models is obtained by combining equations (15-18): Using this likelihood function, it is relatively straightforward to fit the chemical evolution models to the data.We determined the poste-MNRAS 505, 3228-3246 ( 2021)   6 The priors are set uniformly (every model is a priori equally likely) and we use 250 walkers and a burn-in phase of 50 ( ×250) steps followed by 200 ( ×250) used steps.
For the dust parameters, we have also run a grid of models as described in Section 3.5, which will be used in the following section to illustrate the parameter dependencies.Ho we ver, for constraining the dust parameters, more realistic results can be obtained by varying the dust parameters continuously instead of on a discrete grid as for the metal parameters.Unlike the metal parameters, the dust parameters are all numerical values, and can thus be varied continuously.Therefore, we generate an MCMC sample where instead of using discrete indices on a grid of models, we generate the parameter values directly and run a chemical evolution model for each proposal in the MCMC chain.The priors on the parameters are uniform for M destr , k gg,dif , f dif , and log-uniform for SN red , k frag , and k gg,cloud .This division has been used because it mirrors the distribution of grid values for these parameters.The grid values were in turn distributed that way because they lead to similar spacing between models in Fig. 6 (Section 5.2).Table 2 shows the limits of the priors.
The probability for this model is then calculated in exactly the same way as described abo v e.This method is more computationally e xpensiv e as large numbers of models need to be run.Therefore, we use 100 w alk ers and a burn-in phase of 10 ( ×100) steps followed by 100 ( ×100) used steps.For each of these proposals, we run 60 models in order to marginalize our probabilities o v er the 6 SFH and 10 initial masses (for a total of 600 000 models analysed).We call this framework BEDE -Bayesian Estimates of Dust Evolution.

Metal parameters
Before we look at the statistical results, it is beneficial to visualize how our chemical evolution models typically depend on their input parameters and how they compare to the nearby galaxy observations.Since the parameter space is multidimensional and there are a lot of models, it is impossible to compare our models simultaneously in one plot.Therefore, we here select one model as a reference model, and vary the chemical evolution parameters one by one to illustrate the main effects the y hav e on the observables.As reference, we will use a model from the family of models with the highest probability from the MCMC fit (detailed in Section 6.1).Many of the chemical evolution parameters actually affect most or all of the observables, ho we ver, illustrating all these dependencies would take up too much space in this paper.Here, we look at some key plots to illustrate the main variation, with additional focus on the build-up of metals.
Using the MCMC trace described in Section 4, we determine the family of models with the highest probability as the one with a Salpeter IMF, standard recycling of outflows ( f recy = 1), LC18 ( v rot = 150 km s −1 ) SN metal yields, and Karakas et al. ( 2018 ) high mass-loss AGB metal yields (see Section 6.1).Here, we remind the reader that families of models have the same chemical evolution parameters, but the total galaxy mass, SFH, and formation time are allo wed to v ary.The observed data are compared to all models in the family, and thus the general location and spread of the whole family of models is more important than the location of an individual model.In Fig. 2 ( top ), we illustrate this spread and show the build-up of metals (12 + log(O/H)) as gas is converted into stars for models with different M tot that are part of this family.The gas fraction of the galaxy is defined as f g = M g /( M g + M s ).Each point along the plotted In Fig. 3 , we vary the remaining chemical evolution parameters one by one, and again show the build-up of metals with decreasing gas fraction.Here, each line belongs to a different family of models.These are the metal parameters we want to constrain in this work (Section 6.1).The top left-hand panel sho ws ho w the IMF affects the build-up of metals.We find considerable differences in how much metals are released by each IMF.Massive stars e xpel relativ ely more metals and evolve faster.IMFs with a larger fraction of massive stars will thus produce more metals.We indeed see that the topheavy Chabrier IMF described in Section 3.5 produces most metals, followed by the Chabrier IMF, the Salpeter IMF, and finally the Kroupa & Weidner ( 2003 ) IMF.The top right-hand panel of Fig. 3 shows the how the build-up of metals is affected by the SN yield tables used.Theoretical calculations of the amount of metals expelled in a SN event differs among authors, with the yield tables of Meynet & Maeder ( 2002 ) having the largest amount of expelled metals, followed by the rotating (300 and 150 km s −1 , respectively) SN models of LC18 .The yields from Maeder ( 1992 ) are next and finally the non-rotating models from LC18 .Note that here we are not varying multiple chemical evolution parameter simultaneously and are not showing all observables.We thus cannot make much inference about which IMF or yields are best by comparing the models to the data in only this one plot.Even though the non-rotating LC18 yields result to a poor fit to the data in this plot, they result in much a better fit if they are, e.g.combined with a top-heavy IMF.We explore the full metal parameter space in Section 6.1.
The bottom left-hand panel sho ws ho w the recyling time scaling factor f recy (see Section 3.2) affects the slope of the build-up of metals.When f recy is small, outflows are rec ycled v ery fast.As a result, there are not as many metals suspended in the IGM, and the build-up of metals within the galaxy is faster.Vice versa, a large f recy leads to a slow build-up of the galaxy's metallicity as a lot of metals are in the IGM.Changing f recy gives us a way to scale the amplitude of the effect the outflows have on the evolution of galaxies.
In the bottom right-hand panel panel of Fig. 3 , we see that the buildup of oxygen (or the total mass of metals) is little affected by which AGB metal yields are used.The vast majority of metals is expelled by  1 ).Changing the input AGB metal yields changes the N/O ratio in spite of not changing the metallicity (see also Fig. 3 ).Symbols for the observed data points are the same as in Fig. 2 .SN.This does not mean, ho we ver , that A GB yields are unimportant.In Fig. 4 , we sho w ho w the Nitrogen-to-Oxygen ratio changes with metallicity for models with different AGB yields.Here, it is clear that the predicted Nitrogen-to-Oxygen ratio is significantly affected by choice of AGB yields from the literature.Older AGB yield tables such as the ones from van den Hoek & Groenewegen ( 1997 ) or Karakas ( 2010 ) produce too high Nitrogen-to-Oxygen ratios at low metallicities.The other yields used provide a better match, with the Karakas et al. ( 2018) yields with high mass-loss rates providing the best match to the build-up of log(N/O) with increasing metallicity.
Given that our SFR prescription is empirical, we show that the resulting values are sensible in Fig. 5 .Here, we plot the SFR/ M baryon as this normalization allows us to more easily compare massive and dwarf galaxies in terms of the SFR relative to their size.The different SFH together span the parameter space quite well.7Note that the SF bursts are introduced randomly, so each bursty model will have its burst at different gas fractions.The various models thus together span the SFR parameter space e xcellently.F or each combination of chemical evolution parameters, we see there are six different SFHs based on three different reference SFE 0 (See Section 3.5).For each of these, there is one SFH where bursts have been included, and one SFH without bursts.During the bursts, the SFR are increased for a given number of time-steps as discussed abo v e.Because the mass loading factors do not change during the bursts, the increased SFR also infer significant outflows.The SFR and outflows together significantly reduce the gas fraction during the burst.Ho we ver, after the bursts, a significant fraction of the outflows will be recycled in a short time span, especially for more massive galaxies (with lower gas fractions).This results in the increase in gas fraction that can be seen after the bursts at lower gas fractions.

Dust parameters
We illustrate the variation introduced by the different dust parameters by plotting how the dust-to-gas mass ratio evolves with increasing metallicity for both the models and observed data.We choose the model with SN red = 5, M destr = 15 M , k frag = 0.05, k gg,cloud = 4000, k gg,dif = 5, and f dif = 0.2 as our reference model.In Fig. 6 , we vary the dust chemical evolution parameters one by one, compared to this reference model and the data.
The first panel of Fig. 6 shows how changing the SN dust yield affects the build-up of dust as galaxies evolve.When the SN dust yield is reduced by a larger factor (SN red ), the dust-to-gas ratio is reduced at low metallicities.At higher metallicities the difference is much smaller.This is because at early evolutionary stages, SN dust is the dominant contributor, whereas for more evolved galaxies (once the 'critical metallicity' is reached), dust grain growth is dominant.Reducing the SN dust content is a possible avenue to explain the low dust content of galaxies at early evolutionary stages (R émy-Ruyer et al. 2014 ; De Vis et al. 2017aVis et al. , 2019 ) ).
In the next two panels of Fig. 6 , we study the effects of varying the dust grain growth parameters in the diffuse ISM and in clouds.We find the cloud grain growth parameter has a much stronger influence than the diffuse grain growth parameter.For both environments, most variation is introduced at intermediate metallicities.At low metallicities, the grain growth is too inefficient since too few metals are available for accretion on to the dust grains.At large metallicities, grain growth is dominant and all of the available metals have accreted on to the dust grains.For each of the illustrated models, the maximum dust-to-metal ratio has been reached for the highest metallicities.The dust grain growth scaling factors k gg,cloud and k gg,dif thus affect when the grain growth becomes dominant (i.e. when the 'critical metallicity' is reached) and how fast the metal reservoir accretes on to the dust grains.
The fourth panel shows how the models change when we change the fraction of metals that are available for grain growth.This mainly determines the maximum dust-to-metal ratio (and thus dust-to-gas ratio at a given metallicity) that is reached at high metallicities.Given the majority of our models reach this maximum dust-to-metal ratio at high Z , the f dif parameter is the one that most strongly affects the dust-to-gas ratio at high metallicities.
The two bottom panels of Fig. 6 sho w ho w the dif ferent dust destruction mechanisms change the build-up of dust.Destruction of grains by SN shocks is not very efficient, the effect is minor have quite poorly determined SFR and molecular gas masses.The outlying sources also have quite large uncertainties and can thus safely be ignored.even if a single SN destroys 30 M of dust.At low metallicity, the SN shocks are inefficient because the dust-to-gas ratio is very low and so not much dust is present in the gas the shock propagates through.At high metallicities, significant amounts of dust are destroyed by SN, but by this point grain growth has become so efficient that the destroyed dust is rapidly re-accreted on to other dust grains.The photofragmentation of large a-C:H/a-C grains is poorly understood, but the last panel of Fig. 6 demonstrates a large value of k frag results in a reduction in the dust content of lowmetallicity galaxies.This could provide an alternative destruction mechanism for explaining the low dust content of low-metallicity galaxies.

Statistical constraints on metal parameters
By comparing our grid of chemical evolution models to all available observations simultaneously, we obtain probability functions for each of the four input parameters (IMF, mp Z ,SN , mp Z ,AGB , and f recy ).Additionally, it is also possible to make 2D probability distributions between any two of the input parameters.The corner plot showing these probability distributions is shown in Fig. 7 .From the histograms, it is immediately clear that some input values in our grid do not result in realistic models.Decreasing the recycling time scaling factor f recy by a factor of 10, or increasing it by a factor of 2 results in statistically poor fits.Similarly, a Chabrier IMF or top-heavy Chabrier IMF have low probability, as well as SN yields of Meynet & Maeder ( 2002 ) or the non-rotating models of LC18 .The AGByields of van den Hoek & Groenewegen ( 1997 ) and Karakas ( 2010 ) also have low probability.
We now study the 2D probability distributions of the MCMC trace in further detail.We find that the most likely combination of parameters is the model with f recy = 1, Salpeter IMF, LC18 SN yields with v rot = 150 km s −1 and AGB yields from Karakas et al. ( 2018 ) with high mass-loss rates.Another particularly good combination of parameter values is found for the model with f recy = 0.25/0.1,Salpeter IMF, Maeder ( 1992 ) SN yields and AGB yields from Karakas et al. ( 2018 ) with high mass-loss rates.There is a third option that also has non-negligible probability, with a Kroupa IMF, f recy = 0.5, LC18 SN yields with v rot = 150 km s −1 and AGB yields from Karakas et al. ( 2018 ) with high mass-loss rates.In addition to these three most likely options, there are additional variations that give nonneglibiible contributions.
There are thus three groups of models that each have entirely different combinations of IMF, f recy , and SN yield tables yet produce probable matches to the observations in spite of their different numerical values.This demonstrates a de generac y between the IMF, SN yields, and f recy .Using the 2D histograms is paramount for identifying the most likely parameter combinations and revealing degeneracies.

Statistical constraints on dust parameters
Fig. 8 shows the corner plot with the 1D and 2D probability distributions for the continuous MCMC run constraining the dust parameters.Allowing the dust parameters to vary continuously results in relativ ely well-behav ed Gaussian probabilty distributions.The MCMC trace is able to find the most likely parameter values and iterate around the most likely values around it.The results are listed in Table 3 .The uncertainties on most of the parameters are on the order of 50 per cent.The only parameter that is strongly constrained is the fraction of metals available for grain growth in the diffuse   Next, we look at the degeneracies between the different parameters.There are multiple positive and negative correlations between the parameters.The strongest positive correlations observed are (i) the correlation between SN red and f dif , and (ii) between k gg,dif and M destr .These positive correlation can be explained as the effects of these parameters cancel each other out.One of the terms increases the dust content (increasing f dif or k gg,dif ) and the other reduces the dust content (increasing SN red or M destr ).Fig. 6 demonstrates that k gg,dif and M destr have very similar, but opposite, effects.
The strongest ne gativ e correlation is between SN red and M destr , both of which decrease the dust content.Finally, there is also an anticorrelation between SN red and k frag (especially if the values with SN red > 5 and k frag > 0.1 are ignored).The latter correlation arises from the constraints imposed by the low-metallicity galaxies in our sample.In order to obtain a good fit to these unevolved galaxies with our chemical evolution models, SN red or k frag have to be increased (Fig. 6 ).Both reduced SN dust or photofragmentation can explain the low dust masses of low-metallicity galaxies.

Dust budget of best models
To impro v e on our interpretation, we look in more detail at the dust budget of our best 'family' of models.Here, we have run models with different initial masses and SFE and taken the best values from Table 3 as the rele v ant dust parameters.Fig. 9 (left) shows how the time-scales associated with the different dust processing mechanisms vary with metallicity for two different galaxy masses.We see that for our best models with M tot = 10 9 M , the cloud dust grain growth time-scale varies on a linear power-law from about 1 billion years for low metallicities till about 10 million years for high metallicities.The dif fuse grain gro wth is o v er an order of magnitude slower and initially also follows a power law with metallicity.However, towards high metallicities, the diffuse grain growth first levels off as the maximum dust-to-metal ratio in the diffuse ISM is reached (which reduces the diffuse grain growth time-scale as apparent in equation 11) and then increases again as a balance is formed between dust production and dust destruction.
The dust destruction by SN shocks is initially inefficient due to the low dust-to-gas ratios in the low-metallicity ISM.As the dust-togas ratio increases with increasing metallicity, the dust destruction time-scales for SN shocks become faster, but it never catches up with the increasing rate of the dust grain growth.The photofragmentation of large grains is the dominating process (i.e.short time-scales) at low metallicities and then becomes weaker (longer time-scales) as the galaxy evolves.This can be explained in the following way: in unevolved galaxies, there are many young stars and thus there is a very harsh radiation field.In addition, due to the low dust-to-gas ratio the dust does not self-shield and the photofragmentation of dust grains is thus initially very efficient.As galaxies evolve the radiation field becomes less harsh and the large dust grains can survive longer.
Next, we look at how much the different dust production mechanisms contribute to the total produced dust mass (i.e.ignoring all dust destruction) in Fig. 9 (right).We see that at low metallicities, SN contributes nearly all of the dust mass, and AGB stars and grain growth only produce a marginal fraction.Ho we ver, as the metallicity increases, the efficiency of both the cloud and diffuse grain growth increases.At a metallicity of about 12 + log (O / H) = 7 .75, cloud grain growth has produced half of the dust present in the galaxy, and from then on cloud grain growth is the dominant mechanism.This is a significantly lower transition value than for other models in the literature (Kuo & Hirashita 2012 ;De Vis et al. 2019 ;V ílchez et al. 2019 ), we note that these studies did not include photofragmentation.
In order to explain the low dust-to-gas and dust-to-metal ratios observed in multiple observational studies (R émy-Ruyer et al. 2014 ;De Vis et al. 2017a, 2019 ;Cigan et al, in preparation), previous modelling attempts (e.g Zhukovska 2014 ;Feldmann 2015 ;De Vis et al. 2017b ) have often invoked very low SN dust yields.In this work, we have found that including a dust destruction term for photofragmentation of grains can also explain the lower dust-togas and dust-to-metal ratios of low-metallicity galaxies.In this case, the SN dust contribution only needs to be lowered by a factor of 2 to about on average ∼0.5 M of dust per SN.This is more consistent with observations of SN remnants (Dunne et al. 2003 ;Morgan et al. 2003 ;Rho et al. 2008 ;Barlow et al. 2010 ;Matsuura et al. 2011Matsuura et al. , 2015 ; ;Temim et al. 2017 ;Rho et al. 2018 ;Chawner et al. 2019 ) than the more extreme values that are necessary without photofragmentation of large grains.

C AV E AT S
Finally, we want to point out a number of important caveats to this work.The first and most important being that we do not claim our model can be used to state that a given parameter (e.g.Salpeter versus Chabrier IMF) is definitively more realistic than another.Within our statistical framework, one parameter might give a better fit than another and this does count for something, but this is only valid within the assumptions that were made (from the assumption of the outflow prescription to the dust mass absorption coefficient).There are a number of additional factors that need to be considered.
The first consideration is that we have only modelled galaxies evolving in isolation, and have not taken into account their merging histories.In addition, the prescription we used for inflows and outflows might not be ideal for our type of galaxies.The Nelson et al. ( 2019 ) prescription was chosen as the most realistic available, but may not encompass the physical properties of the whole range of galaxies observed in our sample.Additionally, when calculating the redshifts for the outflows, it is assumed that each galaxy was formed shortly after the big bang.When comparing our models to the observ ations, we allo w the age of the galaxy to vary (i.e.we compare to the model tracks rather than just the model end point).This means we are ef fecti vely allo wing the formation time to be much later than the big bang.For models where this is the case, the model will have used outflows based on a redshift that is slightly too large.This introduces an error into our models, but we consider this error to be negligible compared to the o v erall uncertainty associated with the outflow prescription.
In Section 6.1, we found there is a de generac y between the IMF, f recy , and SN yield tables.Rather than varying f recy , it is also possible to change the outflow prescription.This strongly affects the de generac y with IMF and SN yield tables.In the preparation of this study, multiple outflow prescriptions were tried, with significantly different results for what parameters gave the best fit (e.g. the Chabrier IMF was found as best IMF in one iteration).Ho we ver, e ven though the MNRAS 505, 3228-3246 (2021) Downloaded from https://academic.oup.com/mnras/article/505/3/3228/6294485 by guest on 18 May 2023 results were different, a very similar de generac y between IMF, SN yield tables, and outflows was always reco v ered.
Another key caveat is that our 'one-zone' approach without resolution has inherent limitations.We do not resolve different regions within the galaxy (other than the separation of the diffuse and cloud ISM).A real galaxy does, of course, not evolve so uniformly and different regions evolve at different time-scales.This undoubtedly affects the chemical evolution of galaxies.Detailed comparisons between resolved hydrodynamical simulations (which can unfortunately not sample the dust processing parameter space in detail) and 'one-zone' models could shed light on where the shortcomings of the latter lie, and how they can be overcome.
In addition, in this work we have only considered large dust grains without separating between silicate and carbonaceous dust.Silicate and carbonaceous dust do not have the same origins nor the same evolutionary pathways.Furthermore, the carbonaceous dust component should ideally include information about the grain size distribution or be divided into two populations.There would then be one population with silicate dust (not sensitive to photoprocessing), one population with large a-C(:H) grains (prone to photofragmentation), and one with nano a-C grains (subject to photodissociative processing into their hydrocarbon or atomic constituents).Note that the two latter destructive routes are energetic e xtreme-UV photon-driv en but the microscopic processes are different and depend on the interstellar radiation field (photofragmentation, important at the grain sub-structure lev el, v ersus photodissociation at the C-H and C-C bond molecular level; Schirmer et al. 2020 ).Including these different populations could lead to a more detailed understanding of dust evolution (especially when combined with MIR-based observational constraints).
Another consideration is that the priors we have used might not be the most appropriate for this kind of work.Constraints from theoretical modelling could impro v e our prior knowledge of the possible parameter values, which could in turn change our MCMC results.Using a uniform prior (either in logarithmic or linear space depending on the kind of parameter), we have aimed to use the most uninformative prior available, yet impro v ements could be made in future work by constructing realistic priors consistent with theoretical predictions.
Finally, a further impro v ement could also be made to our statistical framework.Currently, we have added measurement and model uncertainties, where the model uncertainty is to account for systematic errors and the limited model sampling.Ho we ver, in our approach these uncertainties are treated as random uncertainties, and thus scale as σ ∝ √ N (where N is the number of observations).In reality, the systematic uncertainty component would result in the same error for each galaxy (e.g. if the dust mass absorption coefficient is biased, this bias would be the same for each galaxy) and would thus not scale with the number of observations.To account for this in our framew ork, we w ould need to use a full error co-variance matrix in the likelihood calculation in equation ( 19).Unfortunately, this would make this calculation too time consuming for the number of galaxies, models, and computational power we are using.
Instead, we have tried to address some of the effects of systematic biases by manually perturbing all the observations.We added a bias of 0.1 dex to each of the observed variables one-by-one, and repeated our analysis.For brevity, we do not show all these results, but only summarize that these perturbations had only a limited effect on the probability distributions and did not affect any of the conclusions in this work.We also experimented with changing the metallicity calibration, where we use the IZI metallicities from De Vis et al. ( 2019 ) instead of PG16 .The results using the IZI calibration are shown in Appendix A. These kind of perturbations can give us some idea of the systematic ef fects.Ho we ver, doing a more realistic calculation including a full covariance matrix, would make the uncertainties on the chemical and dust evolution parameters significantly more realistic, and would be the logical next step.

C O N C L U S I O N S
In this work, we have used a rigorous statistical framework to provide statistical constraints on the chemical and dust evolution parameters for nearby galaxies.The models are compared to 340 nearby LTG observations from the DustPedia, HIGH, and HAPLESS samples.The key effects of varying each of the parameters is illustrated using plots of how the metallicity increases with decreasing gas fraction, and how the dust-to-gas ratio increases with increasing metallicity.A Bayesian MCMC framew ork w as used to provide statistical constraints, where the relative probabilities were calculated from the χ 2 , and we marginalized o v er the different time-steps, 10 different galaxy masses and 6 SFHs.Metallicity-related parameters were studied first, and the best model was chosen for subsequently studying the dust parameters.From studying the statistical constraints, we conclude the following: (i) Our main conclusion from exploring the full parameter space is that there are multiple viable models that compare well to the observations.In other words, there are significant degeneracies between various chemical and dust evolution parameters.It is thus necessary to sample a large parameter space when trying to put statistical constraints on the model parameters.
(ii) There are multiple combinations of metallicity parameters that give a realistic fit to the build-up of metals as galaxies evolve.In particular, there is a de generac y between varying the IMF, f recy , and SN yield tables.The best-fitting combination is the one with a Salpeter IMF, f recy = 1, and Limongi & Chieffi ( 2018 ) SN yields with v rot = 150 km s −1 .We note this fit is particularly dependent on the outflow prescription.
(iii) The Karakas et al. ( 2018 ) AGB yields with high mass-loss rates give consistently the best-fitting results within our statistical framework.Older AGB yield tables such as the ones from van den Hoek & Groenewegen ( 1997 ) and Karakas ( 2010 ) provide a poor fit to the build-up of the nitrogen-to-oxygen ratio as galaxies evolve.
(iv) From varying the dust parameters continuously, it is possible to find the best-fitting values and uncertainties.The best-fitting model has SN red = 2.01, M destr = 15 .74 M , k frag = 0.030, k gg,cloud = 3820, k gg , dif = 5 .59Gyr −1 , and f dif = 0.204.The uncertainties on these parameters are typically quite large (on the order of 50 per cent or more), except for f dif .
(v) There are degeneracies between a number of dust parameters.In particular, there is a positive correlation between SN red and f dif and between k gg,dif and M destr .There is also a ne gativ e correlation between SN red and M destr , and between SN red and k frag .In other words, both SN dust reduction or photofragmentation can explain the low dust content of unevolved galaxies.The best fit is for a combination of the two.
(vi) For the best-fitting models, the grain growth time-scales get shorter following a linear power law with increasing metallicity.The photofragmentation time-scales get longer with increasing metallicity.Nearly all the dust of low-metallicity galaxies was made by SN.The dust mass of high-metallicity galaxies is dominated by grain growth dust.Galaxies with PG16 metallicities around 12 + log (O/H) = 7.75 have similar amounts of SN and grain growth dust.This is a lower transition metallicity than most models in the literature.

Figure 1 .
Figure1.Star formation histories resulting from our empirical SFE prescription for a subset of chemical evolution models (one family of models -see Section 4).The colour of the model indicates its total galaxy mass, and the line type the reference SFE ( SFE 0 = 10 −8 . 5yr −1 : solid, SFE 0 = 10 −9 yr −1 : dashed and SFE 0 = 10 −9 . 5yr −1 : dotted).SFHs including b ursts ha ve been included using lighter shaded lines.The SFHs have the shape of a delayed SFH, which peaks at different times for different mass galaxies.

Figure 2 .
Figure 2. The build-up of metallicity with decreasing gas fraction for our observed samples (data points) and one family of chemical evolution models (the curves).The data are obtained from the HIGH survey (the blue squares), DustPedia (the green and purple hexagons for LTG and ETG, respectively) and HAPLESS (the cyan circles).Top: Changing M tot (sum of initial cloud mass + pristine inflows mass) shows the spread of models in a family [shown by the different colour curves (see legend) with the reference model of 10 10 M in yellow].Bottom: Changes in the build-up of metals due to changing the SFH (with the reference 'average' SFE model shown in red).

Figure 3 .
Figure 3.The variation in the chemical evolution models due to changing the IMF ( top left ), the SN metal yields ( top right ), the recycling scaling factor f recy ( bottom left ), and the AGB metal yields ( bottom right ).Symbols for the observed data points are the same as in Fig. 2 .The red curve shows the reference model in each subplot, see legend for more details on the parameter values used.

Figure 4 .
Figure 4.The variation of the Nitrogen-to-Oxygen ratio with metallicity for the different input AGB metal yields (Table1).Changing the input AGB metal yields changes the N/O ratio in spite of not changing the metallicity (see also Fig.3).Symbols for the observed data points are the same as in Fig.2.

Figure 5 .
Figure 5.An illustration of the various kind of SFHs used in this work.Note that the bursty and non-bursty models o v erlap during the interbursts periods.Symbols for the observed data points are the same as in Fig. 2 .Different coloured curves illustrate different star formation efficiencies.

Figure 6 .
Figure6.The build-up of the dust-to-gas mass ratio with metallicity for our observed samples (data points, see Fig.2) and models (the coloured curves) illustrating the variation in the chemical evolution models due to changing the SN dust yields ( top left ), the diffuse grain growth rate ( top right ), the cloud grain growth rate ( middle left ), the fraction of metals available for grain growth ( middle right ), the mass of dust destroyed by SN shocks for an MW-type galaxy ( bottom left ), and the dust photofragmentation rate ( bottom right ).The reference model for each parameter is shown in red, see legend for more details on the parameter values used.

Figure 7 .
Figure 7. Corner plot showing the probability distributions (see colour bar) on the grid used for studying the metallicity-related parameters.The probability distributions reveal multiple likely combinations (i.e.degeneracies).

Figure 8 .
Figure 8. Corner plot showing the probability distributions (probability indicated by the colour bar) for the continuous MCMC run studying the dust-related parameters.

Figure 9 .
Figure 9. Left: Evolution in the time-scales associated with the different dust processing mechanisms with increasing metallicity.These include grain growth in the diffuse ISM (cyan) and dense ISM (blue), dust destruction in SNe (orange), and photofragmentation (red).The solid and dotted lines indicate different M tot values.Right: Similar to the left image, but now following the fraction of dust mass contributed by different dust sources including SNe (blue), AGB stars (cyan), grain growth in the dense (red), and diffuse ISM (green).Galaxies transition from being SN-dust dominated to cloud grain growth dust dominated around a metallicity of 12 + log (O/H) = 7.75.

Table 1 .
Downloaded from https://academic.oup.com/mnras/article/505/3/3228/6294485 by guest on 18 May 2023 Free parameters in the chemical evolution models together with the grid values used to sample the parameter space.

Table 2 .
Priors on the MCMC parameters constrained in this work.,AGB , and f recy ) in a Bayesian manner, sampling the posterior probability distribution functions using the EMCEE (F oreman-Macke y et al. 2013 ) MCMC package for PYTHON .Ef fecti vely, for the metal results we are sampling the indices for the model grid described in Section 3.5.Since the indices of a grid are discrete, we round each number in each proposed MCMC sample before selecting the model.

Table 3 .
Best values and uncertainties on the dust parameters from the percentiles of their distribution in the continuous MCMC sample.