The dust budget crisis in high-redshift submillimetre galaxies

We apply a chemical evolution model to investigate the sources and evolution of dust in a sample of 26 high-redshift ($z>1$) submillimetre galaxies (SMGs) from the literature, with complete photometry from ultraviolet to the submillimetre. We show that dust produced only by low-intermediate mass stars falls a factor 240 short of the observed dust masses of SMGs, the well-known `dust-budget crisis'. Adding an extra source of dust from supernovae can account for the dust mass in 19 per cent of the SMG sample. Even after accounting for dust produced by supernovae the remaining deficit in the dust mass budget provides support for higher supernova yields, substantial grain growth in the interstellar medium or a top-heavy IMF. Including efficient destruction of dust by supernova shocks increases the tension between our model and observed SMG dust masses. The models which best reproduce the physical properties of SMGs have a rapid build-up of dust from both stellar and interstellar sources and minimal dust destruction. Alternatively, invoking a top-heavy IMF or significant changes in the dust grain properties can solve the dust budget crisis only if dust is produced by both low mass stars and supernovae and is not efficiently destroyed by supernova shocks.


I N T RO D U C T I O N
The first blind submillimetre surveys discovered a population of highly star-forming (100-1000 M yr −1 ), dusty galaxies at high redshift (Smail, Ivison & Blain 1997;Barger et al. 1998;Hughes et al. 1998;Eales et al. 1999). These submillimetre galaxies (SMGs) are thought to be undergoing intense, obscured starbursts (Alexander et al. 2005;Greve et al. 2005;Tacconi et al. 2006;Pope et al. 2008), which may be driven by gas-rich major mergers (e.g. Tacconi et al. 2008;Engel et al. 2010;Riechers et al. 2011;Wang et al. 2011), or streams of cold gas (Dekel et al. 2009;Davé et al. 2010;van de Voort et al. 2011). Observational studies show that SMGs typically have stellar masses of ∼10 11 M (e.g. Hainline et al. 2011;Magnelli et al. 2012), large dust masses (∼10 8−9 M ; Santini et al. 2010;Magdis et al. 2012;Simpson et al. 2014), high gas fractions (30-50 per cent; Tacconi et al. 2008;Bothwell et al. 2013) and solar or sub-solar metallicities (Swinbank et al. 2004;Banerji et al. 2011;Nagao et al. 2012). E-mail: ker7@st-andrews.ac.uk The source of interstellar dust in SMGs is still a controversial issue, particularly whether it originates from supernovae (SNe) or from the cool, stellar winds of low-intermediate-mass stars (LIMS). Recent work has revealed a 'dust budget crisis' (Morgan & Edmunds 2003, hereafter ME03;Dwek, Galliano & Jones 2007;Michałowski et al. 2010a;Michałowski, Watson & Hjorth 2010b;Santini et al. 2010;Gall, Andersen & Hjorth 2011b;Valiante et al. 2011), whereby it is difficult to explain the high dust masses observed in high-redshift galaxies with dust from LIMS. 1 At z > 5 this is further compounded as there is little time for LIMS to produce significant amounts of dust (ME03; Di Criscienzo et al. 2013). The surprisingly constant dust-to-metals ratio measured in galaxies over a wide range of cosmic time also indicates that a rapid mechanism of dust formation is needed (Zafar & Watson 2013, and references therein), requiring dust formation time-scales to be the same order as the metal enrichment time-scale. Although Valiante et al. (2009) and Dwek et al. (2011) argue that asymptotic giant branch (AGB) stars may contribute significantly to the dust budget after only 150-500 Myr (and thus may be a significant source of dust at high redshift), the amount of dust produced is highly sensitive to the assumed initial mass function (IMF). Furthermore, in the former study a high star formation rate (SFR) in excess of 1000 M yr −1 sustained over ∼0. 3-0.4 Gyr is required to build up a significant mass of dust. Due to their short lifetimes, massive-star SNe have long been proposed as a potential source of dust at early times (Dunne et al. 2003bME03;Nozawa et al. 2003;Dwek et al. 2007;Gall, Hjorth & Andersen 2011a).
Observational evidence of dust formation in SN ejecta has come to light recently with SN1987A, Cas A and the Crab nebula remnants containing significant quantities of dust (0.1-1 M ; Dunne et al. 2009;Matsuura et al. 2011;Gomez et al. 2012b). There is now little doubt that dust is formed in SN ejecta (Dunne et al. 2003b;Sugerman et al. 2006;Rho et al. 2008; Barlow et al. 2010;Matsuura et al. 2011;Gomez et al. 2012b;Temim et al. 2012) though the amount of dust which will ultimately survive the SN shocks is still highly uncertain (Bianchi & Schneider 2007;Kozasa et al. 2009;Jones & Nuth 2011). Additionally, dust grain growth in the interstellar medium (ISM; Draine 2009) has been proposed as an extra source of dust in galaxies at both high redshift (Michałowski et al. 2010b;Hirashita & Kuo 2011;Valiante et al. 2011;Calura et al. 2014) and at low redshift (Dwek et al. 2007;Dunne et al. 2011;Boyer et al. 2012;Inoue 2011;Kuo & Hirashita 2012;Asano et al. 2013), which could make up the shortfall in the dust budget of galaxies. The difficulty with determining the origin of dust in galaxies and its lifecycle arises due to a lack of large samples of sources in which to test these issues. Previous authors including ME03, Valiante et al. (2009Valiante et al. ( , 2011, Dwek et al. (2011), Gall et al. (2011b) and Michałowski et al. (2010b) investigated the origin and evolution of dust in high-redshift galaxies, but these were limited to one or two (or, at most, a handful) extreme starbursting systems, selected in a non-uniform way and often missing crucial far-infrared (FIR) photometry spanning the peak of the dust emission.
In Rowlands et al. (2014), a sample of SMGs were carefully selected from the comprehensive data in Magnelli et al. (2012) and galaxy properties were derived for the population by fitting their spectral energy distributions (SEDs) from the UV to the submillimetre in a consistent way. Here, we investigate the origin of dust in these high-redshift SMGs using an updated version of the chemical evolution model of ME03 which incorporates realistic star formation histories (SFHs) for each galaxy, with greater complexity than previous chemical evolution studies have attempted. The sample properties and derivation of the observational parameters are described in full in Rowlands et al. (2014, see also Magnelli et al. 2012. We briefly comment on our sample selection and the SED fitting method in Section 2. In Section 3, we present the updated chemical evolution model which follows the build-up of dust over time, with comparison to the observed properties of SMGs in Section 4. Our conclusions are summarized in Section 5. We adopt a cosmology with m = 0.27, = 0.73 and H 0 = 71 km s −1 Mpc −1 .

D E R I V I N G P H Y S I C A L PA R A M E T E R S F O R S M G s
In order to investigate the physical properties of SMGs, we selected a high-redshift sample from M12 with 1.0 < z < 5.3. Full details of the sample selection and caveats/selection effects are pro-vided in M12 and Rowlands et al. (2014), but briefly, the SMGs are selected from blank field (sub)millimetre surveys (850-1200 µm) which have robust counterparts identified with deep radio, interferometric submillimetre and/or mid-infrared (MIR) imaging. The SMGs are located in fields which have a wealth of multiwavelength observations (GOODS-N, ECDFS, COSMOS and Lockman Hole), which is required in order to derive statistical constraints on galaxy physical properties using SED fitting. Rowlands et al. (2014) used a modified version of the physically motivated method of da Cunha, Charlot & Elbaz (2008, hereafter DCE08 2 ) adapted for SMGs to recover the physical properties of the galaxies in our sample. Briefly, the energy from UV-optical radiation emitted by stellar populations is absorbed by dust, and this is matched to that re-radiated in the FIR. Spectral libraries of 50 000 optical models with stochastic SFHs, and 50 000 infrared models, are produced at the redshift of each galaxy in our sample, containing model parameters and synthetic photometry from the UV to the millimetre. The model libraries are constructed from parameters which have prior distributions designed to reproduce the range of properties found in galaxies. The optical libraries are produced using the spectral evolution of stellar populations calculated from the latest version of the population synthesis code of Bruzual & Charlot (2003). The stellar population models include a revised prescription for thermally pulsing asymptotic giant branch (TP-AGB) stars from Marigo & Girardi (2007). A Chabrier (2003) Galactic-disc IMF is assumed. The libraries contain model spectra with a wide range of SFHs, metallicities and dust attenuations.

SED fitting
The infrared libraries contain SEDs comprised of four different dust components, from which the dust mass (M d ) is calculated. In stellar birth clouds, these components are polycyclic aromatic hydrocarbons, hot dust (stochastically heated small grains with a temperature 130-250 K) and warm dust in thermal equilibrium (30-60 K). In the diffuse ISM, the relative fractions of these three dust components are fixed, but an additional cold dust component with an adjustable temperature between 15 and 30 K is added. The dust mass absorption coefficient κ λ ∝ λ −β has a normalization of κ 850 = 0.077 m 2 kg −1 (Dunne et al. 2000;James et al. 2002). A dust emissivity index of β = 1.5 is assumed for warm dust, and β = 2.0 for cold dust. The prior distributions for the dust temperatures are flat, so that all temperatures within the bounds of the prior have equal probability in the model libraries.
The attenuated stellar emission and dust emission models in the two spectral libraries are combined using a simple energy balance argument, that the energy absorbed by dust in stellar birth clouds and the diffuse ISM is re-emitted in the FIR. Statistical constraints on the various parameters of the model are derived using the Bayesian approach described in DCE08. Each observed galaxy SED is compared to a library of stochastic models which encompasses all plausible parameter combinations. For each galaxy, the marginalized likelihood distribution of any physical parameter is built by evaluating how well each model in the library can account for the observed properties of the galaxy (by computing the χ 2 goodness of fit). This method ensures that possible degeneracies between model parameters are included in the final probability density function (PDF) Figure 1. Example best-fitting rest-frame SED of a high-redshift SMG, with observed photometry (red points) from the rest-frame UV to the submillimetre. The photometry is described in Rowlands et al. (2014). The black line is the best-fitting model SED and the blue line is the unattenuated optical model. The residuals between the best-fitting model photometry and the observed data are shown in the bottom panel. Table 1. Summary of physical properties for the z > 1 SMGs derived from stacking the PDFs derived from MAGPHYS. For each parameter, we use the first moment of the average PDF to estimate the mean of the population, with the variance on the population taken from the second moment of the average PDF minus the mean squared. The error on the mean is simply the square root of the population variance, normalized by the square root of the number of galaxies in the sample. We also list the full range of medianlikelihood parameters values. The parameters are stellar mass M * /M ; dust mass M d /M ; dust-to-stellar mass M d /M * and SFR averaged over the last 10 7 years. log 10 (M * ) log 10 (M d ) log 10 (M d /M * ) log 10 (SFR) Mean 10.80 ± 0.10 9.09 ± 0.09 −1.71 ± 0.10 2.59 ± 0.08 Range (9.87−11.74) (7.89−9.56) (−2.47 to −0.81) (1.05−3.34) of each parameter. The effects of individual wavebands on the derived parameters are explored in DCE08 and Smith et al. (2012a), but we emphasize the importance of using the FIR-submillimetre data from the Herschel Space Observatory 3 (Pilbratt et al. 2010) to sample the peak of the dust emission and the Rayleigh-Jeans slope in order to get reliable constraints on the dust mass and luminosity.
The SEDs of 26 high-redshift SMGs were fitted with MAGPHYS, producing model parameters for each source including stellar mass M * /M ; dust mass M d /M ; dust-to-stellar mass ratio M d /M * and the SFR averaged over the last 10 7 years ψ/M yr −1 . An example best-fitting SED is shown in Fig. 1, and the range of values derived for the SMG sample along with their average properties are listed in Table 1.
Evidence from X-ray studies suggests that many SMGs host an active galactic nucleus (AGN; Alexander et al. 2005), indeed six SMGs in the parent sample of Rowlands et al. (2014) show excess emission in the rest-frame NIR, which may be due to dust heated to high temperatures by an obscured AGN (Hainline et al. 2011). As the MAGPHYS SED models do not include a prescription for AGN emission, we follow Hainline et al. (2011) by subtracting a power-law component from the optical-NIR photometry of the subset of SMGs which exhibit excess emission in the NIR. Note that AGN contribute a negligible amount to the flux at wavelengths longwards of rest-frame 30 µm (Netzer et al. 2007;Hatziminaoglou et al. 2010;Pozzi et al. 2012). Subtracting the NIR power-law results in a reduction in the average stellar mass of the sample by 0.1 dex and a negligible change in the recent SFR. We have excluded two galaxies in the original M12 sample from our analysis as the uncertainties on the parameters due to the subtraction of the power law are too large to make them useful members of the sample (see Rowlands et al. 2014 for details).
MAGPHYS also allows us to recover an estimate of the SFH for each source. We note that whilst the exact form of the SFH cannot be measured using broad-band SED fitting, the best-fitting SFHs are consistent with the physical properties of each SMG. The sensitivity of our results to the SFH is explored in Section 3.4. The SFHs in the MAGPHYS models are parametrized by both exponentially increasing and decreasing models of the form exp(−γ t), where γ is the star formation time-scale parameter which is distributed uniformly between −1 and 1 Gyr −1 . The time since the start of star formation in the galaxy (t) is uniformly distributed between 0.01 Gyr and the age of the Universe at the galaxy redshift. Bursts of star formation are superimposed at random times on the underlying SFH, but with a probability such that 50 per cent of the model galaxies will have a burst in the last 2 Gyr. The strength of the burst is defined as the mass of stars formed in the burst relative to the mass of stars formed in continuous star formation over the lifetime of the galaxy; this parameter ranges from 0.1 to 100. The best-fitting SFHs for the SMGs are shown in Fig. 2. Many of the best-fitting SFHs are bursty: some galaxies have evidence of recent burst(s) of star formation producing a significant fraction of their stellar mass. Others appear to have a smoother (either exponentially declining or increasing) SFH, though as most of these occur over a short period of time (10-100 Myr), they are also 'bursty' in nature. As expected, SMGs are therefore likely to rapidly exhaust their gas supply within ∼100 Myr (Simpson et al. 2014, and references therein).

Figure 2.
Best-fitting SFHs of the 26 SMGs derived from MAGPHYS SED fitting (Section 2, see also Rowlands et al. 2014). The majority of SFHs can be described as 'bursts' of star formation, either because they have a short elevated SFR near the current age, or because their SFHs are so short and extreme they can be considered a burst. These SFHs are a key ingredient of our chemical evolution models. It is important to emphasize that unlike most chemical evolution models in the literature, we will not just use simple parametric SFHs, but we will use SFHs that are consistent with the observed galaxy SEDs.

Physical properties of SMGs
A detailed comparison of the properties of SMGs and low-redshift dusty galaxies is presented in Rowlands et al. (2014). In summary, this sample has an average stellar mass of 6.3 +1.6 −1.3 × 10 10 M (in agreement with Hainline et al. 2011 and M12) and an average SFR of 390 +80 −70 M yr −1 (∼120 times higher than a low-redshift galaxy sample matched in stellar mass, where the average SFR = 3.3 ± 0.2 M yr −1 ). This is consistent with the observed evolution in characteristic SFR of galaxies out to z ∼ 2. The SMGs harbour an order of magnitude more dust with mass 1.2 +0.3 −0.2 × 10 9 M compared to (1.6 ± 0.1) × 10 8 M for low-redshift dusty galaxies selected to have a similar stellar mass. The dust masses derived for the SMGs are consistent with those found in the literature (Santini et al. 2010;Magdis et al. 2012;Simpson et al. 2014). It is not surprising that a high-redshift submillimetre sample would have a higher average dust mass, since moderate dust masses would not be detectable at high redshifts. However, such a selection effect does not account for the much larger space density of galaxies with the highest dust masses at high redshift, since these would have been detected should they exist at lower redshift. This is consistent with the well-documented strong evolution in the dust content of massive, dusty galaxies with redshift, in agreement with Dunne & Eales (2001), Dunne, Eales & Edmunds (2003a), Eales et al. (2010), Dunne et al. (2011), Bourne et al. (2012 and Symeonidis et al. (2013).
In order to compare to chemical evolution models, we need to know the gas fraction and metallicity of the SMGs. As these values are not available for all SMGs in our sample, we compare to literature values derived for similar SMG samples as summarized in Table 2. The gas fractions (f g ) of SMGs are, on average, 30-50 per cent (Tacconi et al. 2008;Riechers et al. 2011;Bothwell et al. 2013) based on CO observations and assuming a conversion from CO luminosity to M H 2 of α CO = 0.8−1.0 M (K km s −1 pc 2 ) −1 . The typical metallicities (Z) of SMGs are found to be solar or slightly sub-solar, albeit with uncertainties due to possible AGN contamination of emission lines (Swinbank et al. 2004;Banerji et al. 2011;Nagao et al. 2012). The dust-tometals ratio η Z out to redshift 6 measured from absolute extinction and metal column densities for a sample of γ -ray burst afterglows and quasar foreground absorption systems is found to be approximately constant (Zafar & Watson 2013). Typical gas-to-dust ratios (η g ) for SMGs are estimated at = 46 ± 25 (Swinbank et al. 2014) assuming an average gas mass of M H 2 = (3.6 ± 1.0) × 10 10 M Table 2. Summary of typical physical properties for SMGs derived from the literature (see main text for details). The parameters are the final gas fraction f gas , metallicity in units of solar metallicity (Z; the ratio of metal mass to gas mass, with Z = 0.019), the dust-to-metal mass ratio (η Z = M dust /M Z ) and the gas-to-dust ratio (η g = M gas /M dust ). (Bothwell et al. 2013). Similar values of η g = 28 +14 −22 were found by Kovács et al. (2006). 4

T H E C H E M I C A L E VO L U T I O N M O D E L
In order to investigate the origin of dust in our sample of highredshift galaxies, we compare the observed dust masses of SMGs to predictions using an updated version of the chemical evolution model of ME03. The model is based on chemical evolution models in the literature (Tinsley 1980;Pagel 1997;Dwek 1998;Calura, Pipino & Matteucci 2008). By relaxing the instantaneous recycling approximation to account for the lifetimes of stars of different masses, the model tracks the build-up of heavy elements over time produced by stars (LIMS and SNe) where some fraction of the heavy elements will condense into dust. Given an input SFH, gas is converted into stars over time, assuming an IMF. The total mass of the system is given by where M g is the gas mass and M * is the stellar mass. The gas mass changes with time as described in equation (2), as gas is depleted by the SFR, ψ(t), and returned to the ISM as stars die, e(t): The first two terms in equation (2) on their own describe a closed box system, the third term describes gas inflow with rate I and the fourth term describes outflow of gas with rate O. Inflows and outflows are discussed further in Sections 4.5 and 4.6 and are (in the simplest form) parametrized as a fraction of the instantaneous SFR. Assuming that mass-loss occurs suddenly at the end of stellar evolution, the ejected mass, e(t) from stars is and the remnant mass is Prantzos, Casse & Vangioni-Flam 2003). τ m is the lifetime of a star of mass m from Schaller et al. (1992), m U is 100 M and m τm is the mass of a star whose age is that of a system where a star formed at t − τ m has died at time τ m . For consistency with the SED fitting method in Rowlands et al. (2014), we adopt a Chabrier (2003) IMF, unless stated otherwise. This takes the form: where m c = 0.079 and σ = 0.69. The IMF is normalized to 1 in the mass range 0.1−100 M . The choice of a Chabrier IMF results in higher stellar dust production than the Scalo or Salpeter IMFs (i.e. compared to the results in ME03 and Dwek 1998), since fewer stars with m < 1 M are produced in a given population, hence less metals are locked up on time-scales of the order of the Hubble time.
The evolution of the mass of metals in the ISM (M Z ) is described by where Z is defined as the fraction of heavy elements by mass in the gas phase, i.e.
The first term of equation (4) describes the metals locked up in stars, and the second term describes the metals returned to the ISM via stellar mass loss (as described in equation 6). Together these two terms describe the evolution of metals in a closed box system. The third term of equation (4) describes an inflow of gas with metallicity Z I and the fourth term of equation (4) describes an outflow of gas with metallicity Z O . The final term M Z, i allows for pre-enrichment from Pop III stars. We set this to zero, but adding pre-enrichment at the expected level of Z i ∼ 10 −4 Z (see Bromm & Yoshida 2011 for a review) does not change any of our results. The mass of heavy elements ejected by stars at the end of their lives is described by where mp z is the yield of heavy elements from a star of initial mass m and metallicity Z, interpolated from Maeder (1992) for massive stars, and van den Hoek & Groenewegen (1997) for LIMS (for progenitor masses up to 8 M ). The integrated yield (p z ) is defined as the mass fraction of stars formed in the mass range m 1 −m 2 which are expelled in the form of heavy element z in equation (7), The evolution of the dust mass will depend on (i) the IMF, (ii) the SFH (iii) the amount of heavy elements produced in stars (the yield), (iv) the dust destruction efficiency and (v) whether or not dust can be formed in the ISM in addition to stellar winds or explosions. Heavy elements are produced by both LIMS and SNe, therefore in the model, the fraction of metals turned into dust is parametrized by a dust condensation 'efficiency' for both SN and LIMS yields. In general, the evolution of dust mass (M d ) with time is described by The first term within the parentheses describes metals locked up in stars, a fraction of which are then recycled into dust through stellar winds. The second term accounts for dust produced from freshly synthesized heavy elements in stars (LIMS and SNe) with initial stellar mass 1 ≤ m i ≤ 40 M . The third and fourth terms account for dust lost in forming stars (astration) and dust lost via destruction processes (parametrized by δ dest , see Section 4.3), respectively. The fifth and sixth terms represent grain growth in the ISM (parametrized by δ grow , see Section 4.4) and dust produced by Pop III stars (set to zero, see equation 4). The final two terms describe dust mass gained or lost via inflow and outflow of gas i.e. (M d /M g ) I and (M d /M g ) O .

Dust produced by stars -δ dust
The dust condensation efficiency in equation (8) (δ dust ) describes the fraction of heavy elements which are incorporated into dust for newly synthesized elements. This can be split into the dust efficiency from SNe (δ sn ) and from the stellar winds of LIMS (δ lims ).

Dust from LIMS -δ lims
In this work, we take the dust condensation efficiencies for LIMS with mass 1 ≤ m i ≤ 8 M from ME03 and apply them to the metal yields of van den Hoek & Groenewegen (1997) Since the amount of dust produced by a population of stars depends not only on the dust condensation efficiency but also the chosen metal yields, and given the wide range of theoretical yields in the literature, here we take a moment to discuss the effect the chosen yields and condensation efficiencies will have on our results. In Fig. 3(a), we compare the dust masses from LIMS assumed here (black solid line) to other literature studies. We compare these with the dust masses from Dwek (1998, also used in Calura et al. 2008 where the dust yield is simply assumed to be 1.0 × mp Z for Mg, Si, S, Ca, Fe and 16 O where C/O < 1 and 1.0 × mp C where C/O > 1; with p Z taken from the metal yields of Renzini & Voli (1981). This rather unrealistic assumption therefore assumes that 100 per cent of the available carbon or oxygen formed in LIMS will condense into dust. It is clear that for all progenitor masses in the range 1 < m i < 5 M , the dust masses from LIMS in Dwek (1998) are an order of magnitude higher than the dust masses predicted by this work and by Ventura et al. (2012). The difference stems from the high dust condensation efficiency assumed in Dwek (1998) and the metal yields used. The van den Hoek & Groenewegen (1997) theoretical yields used in this work are more physical than those of Renzini & Voli (1981), and importantly include the TP-AGB phase consistent with our SED fitting technique (see the comprehensive description in Romano et al. 2010 for a detailed comparison of the available yields in the literature).
The oft-used theoretical AGB dust formation model of Ferrarotti & Gail (2006, see Zhukovska, Gail & Trieloff 2008) also predicts higher dust yields from LIMS than this work, particularly in the range 2 < m i < 6 M ; indeed, using the Ferraroti & Gail dust yields ( Fig. 3a) in our model would lead to 3.6 times more dust from LIMS within 0.1 Gyr. This disagreement is a combination of a choice of (different) input yields and a higher equivalent 'δ lims '. We note that the Ferraroti & Gail dust yields (derived from synthetic stellar evolution models) also disagree with recent results from the self-consistent full stellar evolution model in Ventura et al. (2012). In some cases, the Ferraroti & Gail dust yields for a given AGB star exceeds the total yield of heavy elements produced by the van den Hoek & Groenewegen yields at solar metallicity.
Finally, we have compared our choice of yields (van den Hoek & Groenewegen 1997) with Marigo (2000) and updated LIMS models  Dwek (1998) and Calura et al. (2008, blue dotted line). The purple dashed line shows the yields from the theoretical AGB dust formation models from Ferrarotti & Gail (2006) and Zhukovska et al. (2008, also Valiante et al. 2009), with dust masses from Ventura et al. (2012) in red. The minimum average dust yield per AGB star required to explain observations of high-redshift SMGs (Michałowski et al. 2010b) is shown in the shaded light blue region. (b) We compare the theoretical dust yields from core-collapse SNe (δ SN ) with observations. The yields from TF01 used in this work (black solid line), yields from Dwek (1998) and Calura et al. (2008;blue dotted), and Sarangi & Cherchneff (2013, red solid); the range of expected yields from the theoretical SN dust formation model of Bianchi & Schneider (2007) which includes dust destruction (purple dashed lines). The two lines compare the dust mass which has 'survived' the SN shock expanding into two different densities. The average dust yield per SN required to explain high-redshift SMGs (Michałowski et al. 2010b; ME03) is shown in the shaded light blue region. Observed dust masses from Galactic and nearby young SNRs are indicated in the shaded purple regions (Rho et al. 2008;Dunne et al. 2009;Barlow et al. 2010;Otsuka et al. 2010;Matsuura et al. 2011;Gomez et al. 2012b). The boxes indicate the range of dust mass values derived from IR-submillimetre data as well as uncertainties in the mass of the progenitor stars. from Karakas (2010, which as discussed in Romano et al. 2010 provides the best fit to a range of observations) and find that the van den Hoek & Groenewegen (1997) yields are often lower in the 3 < m i < 6 progenitor mass range. However, we argue in favour of keeping these yields since they are the only uniformly calculated set that includes the super-AGB phase between 5 and 8 M (important for producing much of the chemical enrichment in galaxies, Romano et al. 2010), they still explain many of the observational tests which are well fitted by the yields from Karakas (Romano et al. 2010) and they agree with the updated AGB dust model of Ventura et al. (2012). Ultimately, the lower dust yields from LIMS used here compared to Dwek (1998), Calura et al. (2008) and Ferrarotti & Gail (2006) are compensated for by our choice of Chabrier IMF compared to the Scalo/Salpeter IMFs used in these studies since the Chabrier IMF produces approximately four times more interstellar metals (and therefore dust) after 0.5 Gyr of galactic evolution compared to Scalo/Salpeter. Note that if we were to use the Ferraroti & Gail models combined with the Chabrier IMF, we would still not be able to solve the dust budget crisis. Thus, the major results of this work are robust to changing the stellar yields from LIMS, i.e. our conclusions would not be significantly different if we used the same yields as previous literature studies, though Fig. 3(a) suggests that previous works may have overestimated the contribution from LIMS.

Dust from SNe -δ sn
The dust masses formed per core-collapse SN (m d = mp z δ sn ) are taken from the theoretical model by Todini & Ferrara (2001, hereafter TF01) who predict ∼0.1−1.0 M of dust per SN, depending on the progenitor mass and metallicity. Table 3 (see also Fig. 3b) lists these dust masses for each massive-star SN at solar metallicity for the progenitor range 11 < m i < 25 M . Table 3 also compares the TF01 dust masses with the metal yields per SN expected from stellar evolutionary models in the literature (e.g. Woosley & Weaver Table 3. A list of the predicted dust mass m d from TF01 produced in each SN event for different stellar progenitor mass m i at Z . Also shown are the total metal yields (mp z ) available in the ejecta (WW95: Woosley & Weaver 1995; N06: Nomoto et al. 2006;M92: Maeder 1992) and an estimate of the SN dust efficiency parameter in the model δ sn = m d /(mp z ) (the range indicates the difference in yields WW95, N06, M92) a . Observational values for m d from the literature are also shown; note that only SNRs with Herschel observations have been included here since SN dust masses estimated from MIR photometry of remnants are likely to be lower limits. a The metal mass will be overestimated since this includes all metals produced except for light elements, and it is likely that the metals available to form dust are less (depending on the carbon/oxygen ratio and compounds formed). b The Large Magellanic Cloud remnants (N49, SN 1987A) are in a lower metallicity environment (0.5 Z ) compared to the Galactic SNRs, but the predicted metal and dust yields for these SNe do not change significantly between 0.1 and 1.0 Z .
Theoretical SN dust and metal yields from core-collapse SNe  1995; Nomoto et al. 2006) with indicative values for δ sn for each stellar mass, to allow for comparison with the model. Note the wide range of predicted stellar yields from SNe and the range of expected dust masses, this makes δ sn difficult to pin down observationally.
It has only recently become possible to compare total theoretical dust masses from SNe with observations. FIR and submillimetre observations with Herschel have detected cool (T d ∼ 30-40 K) dust in supernova remnants (SNRs) with masses of ∼0.1 M (Rho et al. 2009;Barlow et al. 2010;Gomez et al. 2012b). There is also evidence that the Cas A and SN1987A SNRs have a more massive colder population of dust (T d ∼ 20 K, m d ∼ 0.4-1.0 M ; Dunne et al. 2003bDunne et al. , 2009Matsuura et al. 2011), and these dust masses are close to the higher end of the range predicted by TF01 (see Table 3 and Fig. 3). However, little is known about how much dust will survive the passage through the shockfront (e.g. Bianchi & Schneider 2007). Assuming that dust destruction in the reverse shock is not efficient, the TF01 model dust yields appear to explain the highest observed dust masses in nearby SNRs at solar metallicity (Table 3, Fig. 3b). Therefore, we set δ sn to reproduce the TF01 models in the mass range 9 ≤ m i ≤ 40 M (elsewhere it is set to zero). In Fig. 3(b), we compare the SN dust masses predicted per progenitor star used here (based on TF01 and Herschel observations) with those used in other chemical evolution studies, in particular the works of Dwek (1998) and Calura et al. (2008) who assume dust condensation efficiencies of 0.8 × mp Z for Mg, Si, S, Ca, Fe and 16 O and 0.5 × mp C (where p C and p Z is taken from the published yields by Woosley & Weaver 1995). We also compare the range of expected SN dust yields from the dust formation model in Bianchi & Schneider (2007), an updated version of TF01, which also includes destruction by SN shock waves expanding into the surrounding ISM. The range of dust masses obtained from FIR/submillimetre observations of Galactic and nearby SNRs are indicated via the shaded purple boxes.
The average dust yield per SNe required to explain dust masses in a handful of high-redshift SMGs (from Michałowski et al. (2010b) and ME03) is also highlighted by the light blue shaded region on this plot. The highest observed dust masses from Herschel and the TF01 models agree well with the their estimates. Finally, we ignore dust formation in Type Ia SNe as recent Herschel observations suggest these events are not contributing a significant mass of dust to the ISM (Gomez et al. 2012a).

Dust destruction -δ dest
Since dust is thought to be removed from the ISM by the sputtering and shattering of dust grains by SN shocks (McKee 1989), it needs to be accounted for in this model. The efficiency of dust destruction is highly uncertain but is assumed to depend on the density and composition of the ISM, and the SN shock velocity (McKee 1989;Jones 2004;Dwek et al. 2007). Depending on the adopted destruction efficiency, the predicted dust mass in a galaxy from chemical evolution models can vary by a factor of 10 for a Salpeter IMF (Dwek et al. 2007). We follow Dwek et al. (2007) by parametrizing the dust destruction as a function proportional to the SN rate. In this case, the dust destruction time-scale τ dest is described by equation (9): where M g is the gas mass and R SN is the SN rate: The parameter m ISM is the effective mass of ISM cleared by each SN event, usually assumed to be 1000 M for typical Galactic interstellar densities of 0.1−1 cm −3 (e.g. Dwek et al. 2011;Gall et al. 2011b, see Section 4.3 for more details). The dust destruction time-scale varies over time depending on the gas fraction and SN rate, with a value of <0.09 Gyr for a gas fraction of 0.5 and SFR of >60 M yr −1 . The destruction parameter in our model is given by δ dest = τ −1 dest .

Dust growth -δ grow
We also include a prescription which accounts for accretion of atoms on to dust grain cores in the cold, dense regions of the ISM (Dwek & Scalo 1980;Tielens 1998;Draine 2009). The shortfall of dust from stellar sources in galaxies supports evidence that grain growth is a significant contributor to the dust budget (e.g. Zhukovska et al. 2008;Draine 2009;Michałowski et al. 2010a;Pipino et al. 2010;Dunne et al. 2011;Gall et al. 2011a;Jones & Nuth 2011;Valiante et al. 2011;Boyer et al. 2012;Kuo & Hirashita 2012;Asano et al. 2013). We follow the prescription of , where the rate of grain growth is linked to the metallicity (as grains accrete metals in order to grow), and the SFR, assumed to be proportional to the amount of molecular gas in a galaxy since dense regions (where molecular gas is present) is also the environment where grain growth is likely to take place. The time-scale for grain growth τ grow in the ISM is then given by where η d is the dust-to-gas ratio and Z is the metallicity. τ 0 is defined in  as where ψ is the SFR and is an efficiency parameter (which is unconstrained). The growth parameter in our model is given by δ grow = τ −1 grow (see Section 4.4 for more details). We define τ −1 o as proportional to the SFR because it is mathematically convenient, but we note that in a more correct physical model of τ −1 o may be proportional to the molecular gas density. Note that the parametrization of the grain growth does not make much difference in the present context.

A more realistic treatment of SFH
The detailed treatment of the lifetimes of stars of different stellar masses is important for SMGs which may have bursts of star formation which occur on short time-scales. Previous studies of chemical evolution in SMGs (e.g. ME03) often assumed an SFR proportional to the gas mass which decreased smoothly with time. One of the main differences between this work and ME03 (among others) is the incorporation of a more realistic SFH with bursts of variable strength and duration, and with an underlying SFH which can be either exponentially rising or declining (indeed 15 of our SMGs have exponentially rising SFHs, see Fig. 2). Incorporating the SEDderived SFHs allows us to carry out chemical evolution modelling in a manner that is consistent with the SED fitting method.
Degeneracies between the MAGPHYS-derived SFH parameters, such as the timing, strength and duration of bursts, means that the best-fitting SFH may not be a unique solution. Since the SFHs are constrained predominantly by the UV-optical light which is emitted mostly by stars <100 Myr old, there is a large uncertainty on the form of the SFH at times prior to this. Our approach, however, is correct in the statistical sense for the SMG sample (if not for each individual galaxy), such that each best-fitting SFH produces the observed UV-submillimetre SEDs of these galaxies. Whilst differences between the SFHs can cause some variation in the dust mass, the overall mass of dust produced will ultimately depend on the mass of metals formed, which is governed by the total mass of stars formed (Edmunds & Eales 1998). The adopted SFH gives a physically plausible and self-consistent representation of the SFH which we can use as an input to our chemical evolution models. Although a small number of studies have used physically realistic SFHs consistent with galaxy physical properties (Valiante et al. 2009(Valiante et al. , 2011, our method represents an improvement over previous works that use arbitrary SFHs (or indeed just one SFH to describe all galaxies) which may not be appropriate.
In addition to the derived SFHs for each individual SMG, we also use four fiducial models to explore the changes in dust mass when key parameters such as the SFH and IMF are varied. These fiducial models represent the whole range of continuous SFHs as derived by the MAGPHYS SED-fitting and thus provide a test of the range of entire dust masses expected to be formed. The fiducial models are parametrized by (i) an exponentially declining SFH with initial SFR of 150 M yr −1 (ii) an exponentially increasing SFH with final SFR of 150 M yr −1 , (iii) a constant SFR of 150 M yr −1 and (iv) a burst SFH with all the stellar mass produced in a burst of 1000 M yr −1 . Each model has an initial gas mass set to the median initial gas mass of the SMG sample (see Section 4). The fiducial SFHs reach the mean stellar mass of the SMG sample (Section 4) at 0.9, 4.4, 0.6 and 0.1 Gyr after the onset of star formation, respectively. The short-burst fiducial SFH produces a much lower mass of dust at a given gas fraction compared to the other three SFHs. In total, the variation between the dust masses built up by the model for the different fiducial SFHs is a factor of 29 (if dust is formed only by LIMS) and a factor 1.5 if dust is contributed by both SNe and LIMS.

T H E O R I G I N O F D U S T I N S M G s
We now explore how well different chemical evolution models can reproduce the dust masses of our SMG sample, and other observational properties of SMGs as a population (Section 2.2, Table 2). To model these sources, in the first instance, we consider a closed box model, assuming no inflow or outflow of gas or metals. The initial gas mass is set at 2× the best-fitting stellar mass derived from the SED fitting, such that at the end of the SFH history, ∼50 per cent of the total galaxy mass ends up in stars. 5 The initial gas masses for this sample (with median M g (0) = 1.25 × 10 11 M ) are therefore tuned to reproduce the observed gas fractions of SMGs. In the closed box model, the average final gas mass of the SMGs (5.8 × 10 10 M ) is in agreement with observations. Furthermore, by design, the final stellar masses are in close agreement with the best-fitting stellar masses (mean of the PDF is 6.3 +1.6 −1.3 × 10 10 M ) derived from the SED fitting (Table 1; Rowlands et al. 2014).
We then use our model to go beyond the simple closed box by including the effects of inflows and outflows on the gas, metals and dust. We also explore the effect of dust destruction and grain growth on the dust mass. A summary of all of the observational results derived from the different chemical evolution models considered in this work are given in Table A1, with the 'final' dust masses summarized in Table 4. Table 5 shows a set of 'good' models which reproduce many of the observed properties of the SMG sample. The

Dust production in stars -LIMS only
In the first instance, we consider dust production in a closed box from LIMS only for each SMG in the sample, assuming no dust destruction as an optimistic case. The dust produced by LIMS only is indicated by the solid black line in Fig. 4. The delay between the onset of star formation (as traced by the build-up of stellar mass shown as the solid grey line) and significant dust production by LIMS is evident in these plots, requiring more than a few hundred Myr to reach dust masses greater than 10 7 M . In Fig. 5(a), we show the difference between the median-likelihood dust mass from the SED fitting and the final dust mass derived from the chemical evolution modelling. The dust masses calculated from the chemical evolution model with dust from LIMS fall far short of the observed dust masses for the majority of SMGs. On average, the theoretical dust masses are 5.0 × 10 6 M , which is a factor 240 lower than the average observed dust mass in the SMG sample (1.2 +0.3 −0.2 × 10 9 M ). This provides definitive evidence that (without changing the IMF) the majority of dust in SMGs must come from a source other than LIMS. Although noted by previous authors, this was based previously on smaller samples and/or simple parametrized SFHs (e.g. ME03; Matsuura et al. 2009;Michałowski et al. 2010b;Dunne et al. 2011;Gall et al. 2011a). Whilst different SFHs can produce Table 5. Summary of the properties derived from different chemical evolution models which best describe most of the observed average properties of the 26 z > 1 SMGs (see Table A1 for a list of all the model results). The properties are the final gas fraction f gas , metallicity in units of solar metallicity (Z; the ratio of metal mass to gas mass, with Z = 0.019), the dust-to-stellar mass ratio (M d /M * ), the dust-to-metal mass ratio (η Z = M d /M Z ) and the gas-to-dust ratio (η g = M g /M d ). For reference, the average f g of SMGs is 30-50 per cent (Tacconi et al. 2008;Riechers et al. 2011); the typical Z is solar or slightly subsolar (Swinbank et al. 2004;Banerji et al. 2011;Nagao et al. 2012); the mean log 10 (M d /M * ) is −1.71; the typical η Z is ∼0.5 (Zafar & Watson 2013); and average η g values are ∼30−50 (Kovács et al. 2006;Swinbank et al. 2014). For each chemical evolution model, we list the median value of the sample in bold and the 16th and 84th percentiles to indicate the range of values in the sample. differences in the dust mass of up to a factor of 29 (Section 3.4), the dust deficit seen here is far larger than this, and occurs for all SMGs in our sample regardless of the SFH. Thus, this conclusion is robust for our ensemble of SMGs.
For the closed box model, the median metallicity for the SMGs reaches 0.9 Z ; this will be discussed further in Section 4.5. The median fraction of metals in the ISM in the form of dust (η Z = M d /M Z ) in the LIMS-only model is 0.4 per cent; this is well below the dust-to-metal ratios observed in local galaxies and out to redshifts of 6 (Whittet 1992;Pei, Fall & Hauser 1999;James et al. 2002;Watson 2011;Zafar & Watson 2013) which are typically 50 per cent. The closed box model reproduces the observed average gas fraction in SMGs (30-50 per cent), but in reality in a sample of galaxies there will be a range of gas fractions. We therefore briefly explore the effect of increasing the initial gas mass in each galaxy by a factor of 2 to M g (0) = 2.5 × 10 11 M . This results in a decrease in the median final dust mass by a factor of 1.3 to M d 3.9 × 10 6 M . The median final metallicity of the sample in this case is 0.5 Z , lower than observed metallicities of some SMGs (see Section 4.5). The metallicity and the dust mass are decreased compared to the original model because the same mass of stars enriches a larger mass of gas. The stellar mass is unchanged as this is set by the SFH, resulting in a final average gas fraction of 0.75. Although this is larger than the average observed values in SMGs, some high-redshift systems have been found with f g ∼ 0.8 (Tacconi et al. 2010;Riechers et al. 2011) which is consistent with this model. In summary, even with realistic SFHs including bursts and a larger SMG sample, we can definitively rule out a model with LIMS stardust alone, since this cannot explain the observed dust mass or the dust-to-metals ratio of the SMG population.

Dust production in stars -adding SNe
If we include dust production from both SNe (using the yields from TF01) and LIMS, dust builds up more rapidly in SMGs with a delay of only tens of Myr between the highest mass stars forming and evolving to the SN phase. This is evident in Fig. 4, as the dust produced by both SNe and LIMS closely tracks the stellar mass build-up over time, with bursts of star formation resulting in an almost instantaneous increase in the dust mass. Adding dust from SNe accounts for more than an order of magnitude increase in the dust mass of SMGs (with a median mass of 1.9 × 10 8 M for the LIMS+SNe model) compared to the dust mass from LIMS only (5.0 × 10 6 M ). The model dust masses using LIMS and SNe match the observed values (accounting for the ±0.2 dex uncertainty in the MAGPHYS-derived dust masses) in ∼19 per cent of cases (Table 4). The median metallicity of the SMGs in this model is the same as with LIMS only (0.9 Z ), but with the inclusion of SN dust the median fraction of metals in the ISM in the form of dust is higher (η Z = 16 per cent). In Fig. 5(b), it can be seen that the predicted dust masses for the majority of the SMGs falls short of the observed dust masses, which indicates additional sources of dust, or even higher SN dust yields than TF01 (where δ sn ∼ 0.3-0.8, Table 3) are required.
In Fig 5(c), we consider the extreme case of maximal dust production from SN e.g. δ sn = 1, such that all metals ejected in this phase are incorporated into dust. Sufficiently high dust masses are achieved in this scenario (a median of 7.6 × 10 8 M ) to account for the observed dust in 15/26 SMGs (see the black dashed lines in Fig. 4). The gas-to-dust ratio for the maximal SN dust model is in agreement with observed values for SMGs of 28 +14 −11 and 42 ± 25 6 ( Kovács et al. 2006;Swinbank et al. 2014). The resulting median fraction of metals in the ISM in the form of dust is η Z ∼ 68 per cent, somewhat higher than the values observed in nearby galaxies (James et al. 2002;Watson 2011, and references therein) and in galaxies out to z ∼ 6 (Zafar & Watson 2013). Herschel observations of SNRs do indicate high condensation efficiencies with δ sn ∼ 0.1−1.0 (see e.g. Table 3), suggesting that at some point, the majority of heavy elements produced in SN ejecta can condense into dust. However, Figure 4. The stellar and dust mass evolution over time derived from chemical evolution modelling for the sample of 26 SMGs. Note that the y-axis labels of each panel are different on the right and left sides. The stellar mass growth from the input MAGPHYS SFH is represented by the grey line and corresponds to the left axis. All of the other lines represent different dust models and correspond to the right axis. The black solid line is the dust mass produced by LIMS only, the black dotted line is LIMS and SN dust and the black dashed line is LIMS and maximal SN dust production. The red line represents the dust mass in a model where dust is produced by LIMS and grain growth, and the blue line shows the dust mass if dust produced by LIMS is efficiently destroyed by SN shocks. At early times, dust destruction and grain growth models have a dust mass track similar to that with dust from LIMS only. Horizontal dot-dashed grey and blue lines represent the observed best-fitting stellar masses and median-likelihood dust masses, respectively, with the blue shaded region indicating the 84th-16th percentile range from the SED fitting. evidence for efficient dust production is observed in only a handful of SNRs, and it is unclear how much dust will survive passage into the ISM. Furthermore, uncertainties in the progenitor mass makes δ sn difficult to estimate observationally. It is therefore still possible that SN dust is not able to account for all the dust in galaxies, i.e. either a significant mass of dust must form rapidly in the ISM (see Section 4.4) or theoretical metal yields from SNe obtained from stellar evolution models are systematically underestimated. Whilst differences between SFHs can cause variations in dust mass, as explored in Section 3.4, this is only a factor of 1.5 for models with LIMS and SN dust. Since this is much smaller than the dust shortfall our conclusions are robust to differences between SFHs.

The effects of dust destruction
There are large uncertainties about the effectiveness of dust destruction in the ISM, yet it is often used in chemical evolution analyses. If SN shocks are efficient in destroying dust, then assuming 1000 M of ISM is cleared by each SN event, the dust mass we obtain is reduced by a factor of 6-10 on average if dust is produced by LIMS and LIMS+SNe, respectively. This compounds the dust budget crisis further. The effect of dust destruction on the build-up of dust is shown as the dark blue line in Fig. 4. With LIMS as the only source of dust, the dust destruction efficiency in our model is similar to the maximum dust destruction case in Gall et al. (2011b) with m ISM = 800 M . If dust is produced by both LIMS and SNe, then any increase in the dust mass by including SNe is effectively cancelled out by the dust destroyed, resulting in a median dust mass of 1.8 × 10 7 M , approximately 3.6 times the mass obtained with the LIMS only model (5 × 10 6 M ).
If the dust is shielded in cold, dense regions of the ISM, then it is possible that the dust destruction efficiency of SN shocks will be reduced. Lower dust destruction rates have been suggested by Dwek et al. (2007Dwek et al. ( , 2011 and Gall et al. (2011b), who struggle to produce the dust masses of high-redshift galaxies with efficient dust destruction. These authors suggest that m ISM = 100 M may be more appropriate given the increased density of the ISM gas (n > 10 3 cm −3 , e.g. Alaghband-Zadeh et al. 2013, see fig. 2 in Dwek et al. 2007) compared to the Milky Way. Including dust destruction in the model with m ISM = 100 M now lowers the dust mass by a factor of 1.2-1.6 (compared to the LIMS only and LIMS+SNe models). A cautionary note here is that theoretical dust destruction models are not well understood or appear to be too efficient at destroying dust grains, the very fact that we observe so much dust in galaxies, including those with very little recent star formation (e.g. Rowlands et al. 2012), implies the destruction rate must be balanced by the injection rate from stars and another source of dust e.g. grain growth in the ISM (see also Michałowski et al. 2010b;Dunne et al. 2011;Boyer et al. 2012;.

Adding grain growth
In order to obtain a minimum τ grow ∼ (20-200) Myr in our fiducial models, in line with expected grain growth time-scales and to relieve the 'budget crisis' (Zhukovska et al. 2008, Mattsson, Andersen & Munkhammar 2012, we set = 500 (Section 3.3, equation 11). If the value of is lower, then the grain growth time-scale is longer which reduces the dust mass produced by grain growth at a given age. It is possible that is larger than the value adopted here, however, all of the metals in each SMG are rapidly incorporated into dust grains.
The effect of adding grain growth to the LIMS only model (with no dust destruction) on the dust mass is shown in Fig. 4 by the solid red line. However, this model fails to adequately reproduce the observed dust masses for the majority of SMGs. In Fig 5(d), we therefore consider dust produced by LIMS, SNe (using TF01 yields) and grain growth in the ISM. We find that including grain growth on average increases the dust mass by a factor of 200 to 9.8 × 10 8 M , compared to a model in which dust is contributed by LIMS only. Grain growth can therefore easily make up the shortfall in the predicted dust masses for 60-70 per cent of SMGs in the sample (Table 4), but only if grain growth is the dominant form of dust production in SMGs (see also Mattsson 2011). Indeed, this model can easily account for the observed dust masses of the majority of SMGs in our sample; however, in this scenario, a large fraction (77 per cent) of the metals is locked up in the form of dust which is higher than typically observed. At first glance the large fraction of metals locked up in dust is in conflict with Zafar & Watson (2013) who found a relatively constant dust-to-metals ratio of 0.5 out to z ∼ 6. However, the uncertainties on the dust and metal masses are likely to be a factor of a few; therefore, the average SMG dust-to-metals ratio is not strictly in disagreement with the results of Zafar & Watson (2013).
In order to explain all of the dust in every SMG in our sample, we run the most optimistic models where dust is produced by LIMS, SNe (with all metals incorporated into dust) and grain growth, with no destruction. Even with such high dust condensation efficiencies, we cannot account for the dust in all SMGs in our sample, only reproducing the dust masses for ∼69 per cent of SMGs. Whilst the average dust mass and gas-to-dust ratio in this scenario agree well with observed values, around 96 per cent of metals are in the form of dust, which is much higher than observed. Given that the majority of metals are already in dust grains, increasing the grain growth efficiency does not substantially improve the dust yield and therefore grain growth cannot explain all of this shortfall. The close agreement in the average dust mass for the SMG sample (1.2 +0.3 −0.2 × 10 9 M ) and the average mass of metals (assuming a metallicity of ∼Z and a gas mass of 5 × 10 10 M ) further suggests that metal yields of stars may be systematically underestimated (see also Mattsson 2011).
Ultimately, if efficient dust destruction is included along with LIMS, SN dust and grain growth, then the dust produced in these galaxies is not enough to account for the observed dust masses, with the median dust mass reaching 6.0 × 10 7 M . For the SMGs whose predicted dust masses fall short of the observed value, it is possible that dust destruction is less efficient than that assumed in the literature. Whilst there are considerable uncertainties in the sources of dust production and destruction in galaxies, we can definitively state that LIMS cannot be the only source of dust in SMGs, and show that this result is robust to larger samples and bursty SFHs. In order to explain the observed dust masses of 19 per cent of the SMGs in this work, we require dust from LIMS and ∼20 per cent of metals produced in SN to condense into dust grains and survive (δ dest = 0). If (as is more likely) δ dest > 0, another dust source is required to account for the observed masses. In the next sections, we look at the effects of inflows and outflows on the dust properties of SMGs.

Inflows
The closed box chemical evolution model for the SMGS produces a (median) final metallicity of 0.9 Z (Tables 5 and A1) in line with solar/sub-solar metallicities observed in Swinbank et al. (2004) and Banerji et al. (2011). However, although the closed box model is the simplest approach to chemical evolution, in reality, galaxies are unlikely to be closed systems (e.g. Erb 2008). Examples of this include the well-known G-Dwarf problem, which requires infall of material in the Milky Way (e.g. van den Bergh 1962;Searle & Sargent 1972;Pagel & Patchett 1975;Tinsley 1980), and the observed wide range in stellar metallicity for galaxies at a fixed gas-phase oxygen abundance (Gallazzi et al. 2005) which also requires inflows and outflows of gas. As yet, there is limited direct observational evidence for gas inflows, but recent studies at high redshift provide further indirect evidence that inflows are required, including the need to sustain high SFRs (Giavalisco et al. 2011;Reddy et al. 2012;Tacconi et al. 2013) as well as being an essential ingredient in galaxy formation simulations at these epochs (e.g. Dekel et al. 2009). Due to the importance of gas accretion in galaxy evolution models, we therefore investigate the effect of inflows on our results from the chemical evolution model. In general, inflows will act to decrease the gas-phase metallicity (as the enriched interstellar gas is diluted by the metal-poor inflow, Edmunds 1990Edmunds , 2001Edmunds & Eales 1998), though the dust mass contributed by stellar sources is largely unchanged. Conversely, inflows will have a more significant effect on the dust produced by grain growth by decreasing the grain growth time-scale (equation 11). To determine the 'inflow prescription' to include in this work, we use an inflow rate of the order of the SFR to be consistent with the semi-analytic model of Dutton, van den Bosch & Dekel (2010) and Erb (2008), who find that the rate-of-change of the gas mass (inflow-outflow) is in a steady state with the SFR.
We initially assume that an inflow delivers metal-free gas (Z I = 0, equation (4) and (M d /M g ) I = 0, equation 8) to the galaxy at a rate proportional to the SFR throughout the lifetime of each SMG. We tune the initial gas mass such that the SMGs have the same final gas fraction as the closed box model ∼0.5 to provide a consistent comparison. To demonstrate the effect of inflow on the chemical evolution of SMGs, we run the model including dust from LIMS only, with no destruction or grain growth. We find that an inflow rate equal to the SFR reduces the median metallicity of the SMG sample to 0.7 Z , whilst the median dust mass is increased by a factor of 1.3 (Tables 5 and A1). This is because the inflow model starts with a smaller mass of gas which is enriched to a higher metallicity than the closed box model, which results in a higher dust mass. As the inflow of pristine gas continues, this dilutes the metallicity of the gas but does not affect the mass of dust in the galaxy.
The inflow rates we require to match the observed properties of SMGs are <1 × SFR (for the SMG sample the 16th-84th percentile range of the median SFR over time is 6-600 M yr −1 ), which is consistent with indirect observational support from some studies of high-redshift galaxies, but higher than the gas accretion rates required in simulations (typically 40-60 M yr −1 ; Kereš et al. 2005;Davé et al. 2010;van de Voort et al. 2011). Unfortunately, the large range in the observed metallicities of SMGs and uncertainties due to possible AGN contamination of the emission lines do not allow us to discriminate between models which have different gas inflow rates.

Outflows
Outflows of gas are thought to be common in actively star-forming galaxies at all epochs (Heckman et al. 2000;Weiner et al. 2009;Rubin et al. 2010;Diamond-Stanic et al. 2012;Bradshaw et al. 2013), and may be either driven by stars (stellar winds and SN), or by AGN. Significant outflows of material are implied by the results of Ménard et al. (2010), who found evidence for dust in galaxy haloes with a mass comparable to that of dust in the disc. Furthermore, Erb et al. (2006) suggest that the mass-metallicity relation at z ∼ 2 is modulated by metal-rich outflows from galaxies, with rates of up to four times the SFR. The upper end of their range agrees with results from Dunne et al. (2011), who used a simple analytic chemical evolution model to demonstrate that outflow rates of four times the SFR best describes the evolution of the dust mass function of Herschel-Astrophysical TeraHertz Large Area Survey (H-ATLAS) galaxies at 0 < z < 0.5. Outflows could therefore be responsible for the significant metal enrichment of the IGM and are an important component of the chemical evolution of galaxies.
In this work, we assume that the gas and dust in the ISM are well mixed so that Z O = Z (the so-called unenriched outflow, Dalcanton 2007, see also equation 4), and that outflows remove material (including gas, metals and dust) from the galaxy at a rate proportional to the SFR. To demonstrate the effect of outflows on the chemical evolution of SMGs, we again run the model including only dust from LIMS, with no destruction or grain growth. In general, outflows decrease the overall dust mass (though not the contribution from stellar sources) and decrease the ISM metallicity. In this work, as with the inflow model, the initial gas mass is tuned such that the SMGs have a final gas fraction of 0.5. Outflows of 1× the SFR reduce the dust mass in SMGs on average by a factor of 1.4, and the ISM metallicity is reduced to 0.7 Z . This is well within the large range of the observed dust masses and metallicities of SMGs; therefore, outflows of gas equal to the SFR can be accommodated in our chemical evolution models.
It is also possible that both inflows and outflows occur simultaneously (Sakamoto et al. 2013), or in short succession (Dalcanton 2007). By allowing simultaneous inflow and outflow in our model, with rates equal to the SFR, the metallicity is decreased to 0.6 Z . The dust mass from LIMS is reduced by a factor of 1.3 compared to the closed box model, which is well within the observational range of SMG dust masses. The amount of dust removed from the galaxy in the outflow model is much lower than that suggested by the Ménard et al. (2010) results, which imply a higher outflow rate is needed to remove half of the dust mass. However, outflow rates significantly larger than the rate of gas inflow would serve to decrease the dust mass, which increases the tension between model and observed dust masses and compounds the dust budget crisis further.

Variations in the IMF
The amount of dust formed in galaxies is strongly linked to the amount of metals produced by stars. Increasing the yields and decreasing the amount of low-mass 'dead' stars produced in a stellar population by varying the IMF is therefore one way to possibly solve the dust budget crisis. Previous works (e.g. Dwek et al. 2011;Gall et al. 2011a;Valiante et al. 2011) have shown that by invoking a top-heavy IMF, one can easily reproduce the observed dust masses in some high-redshift galaxies, but it is unclear whether the IMF in other galaxies is Milky Way like and invariant with time and location (see the review in Bastian, Covey & Meyer 2010). Many studies have suggested that a top-heavy IMF is a natural consequence of the extreme environment in high-redshift galaxies, for example due to bursty SFHs and the denser ISM in comparison to local galaxies (Dabringhausen, Kroupa & Baumgardt 2009;Papadopoulos et al. 2011;Kroupa 2012). Indeed, Gunawardhana et al. (2011) found evidence for a strong relationship between SFR and IMF slope for z < 0.3 galaxies, such that galaxies with higher SFRs form more massive stars in a given stellar population. Here, we investigate the sensitivity of the derived dust mass to the IMF in the models and whether this allows us to predict the slope of the IMF required to resolve the dust budget crisis in SMGs.
We increase the power-law slope of the Chabrier IMF (φ(m) ∝ m −α where α is the slope) from −1.3 to −0.67, but leave the low-mass end (<1 M ) unchanged. The value of the high-mass slope is found by extrapolating the relationship of Gunawardhana et al. (2011) between IMF slope and SFR for low-redshift starforming galaxies to the average SMG SFR (390 M yr −1 ). Using the fiducial SFHs introduced in Section 3, i.e. an exponentially declining SFH with ψ(0) = 150 M yr −1 , an exponentially increasing SFH with ψ(f) = 150 M yr −1 , a constant SFR of 150 M yr −1 and an instantaneous burst, a top-heavy IMF does not increase the dust mass enough to account for the observed dust masses of SMGs with an LIMS-source of dust only, even with the increase in the number of 'super-AGBs' formed. Considering dust production from both LIMS and SNe and using the fiducial SFHs, we find that an IMF slope of −0.67 reproduces the average observed SMG dust masses (1.2 +0.3 −0.2 × 10 9 M ) within a factor of 2. This slight shortfall in dust mass can be alleviated by including a small amount of grain growth which allows the average SMG dust mass to be easily reached.
A top-heavy IMF also leads to a higher destruction rate because of the increased SN rate (equations 9 and 10; Gall et al. 2011b;Mattsson 2011). Therefore, the increase in the dust mass from LIMS and SNe with efficient destruction (m ISM = 1000 M ) achieved with a top-heavy IMF is on average a factor of 1.3 lower compared to a Chabrier IMF with the same dust sources and destruction rates (at a time of 0.5 Gyr after the onset of star formation). In summary, invoking a top-heavy IMF with no dust destruction and no grain growth can solve the dust budget crisis in SMGs, but given the uncertainties involved, the high dust masses do not provide unequivocal evidence for a top-heavy IMF.

The dust emissivity caveat
On a final note, it is also possible that dust produced by SNe in high-redshift galaxies has different properties to that assumed in this work, where the dust emissivity κ used to determine the dust mass is calibrated from observations of the Milky Way and other nearby galaxies (Section 2, see also Valiante et al. 2011). If, for example, the dust emissivity in high-redshift SMGs was systematically higher due to changes in the dust grain size distribution, shape or grain composition, this would serve to decrease the observed dust masses thereby alleviating the tension between observed and predicted properties of SMGs. One possibility to increase κ is if the dust grains are in dense gas environments and have amorphous structures (Ossenkopf & Henning 1994). Although there is no observational evidence for different grain properties in SMGs, recent FIR observations of the Magellanic Clouds found that dust properties may be different from those in the Milky Way (Meixner et al. 2010) due to a recent increase in the Type II SN rate. It is therefore plausible that the dust grain properties may be different in environments where dust is produced and/or reprocessed by SNe. Indeed, Bianchi & Schneider (2007) derive κ values for freshly formed SN dust and SN-processed dust, with the latter case predicting a κ 850 value 7 which is 2.46 times higher than the value assumed in this work. If dust is produced by LIMS only, κ would need to differ by a factor of 240 compared to the MW in order to resolve the dust budget crisis in this sample. Adding an SN dust source with no grain growth or dust destruction, the crisis can be solved with κ increasing on average by a factor of 7 at 1 < z < 5. Note that to explain the dust shortfall in individual SMGs by dust emissivity variation alone, κ 850 would need to vary from 0.074−5.9 m 2 kg −1 . If we include efficient dust destruction in our models, then κ 850 would need to be a factor of 70 higher (i.e. 5.4 m 2 kg −1 ) at 1 < z < 5 to solve the dust budget crisis. In comparison, the range of κ 850 in the literature (Valiante et al. 2011, see references therein) is 0.02−2.0 m 2 kg −1 where the higher end of this range refers to the results from theoretical simulations of dust processed by SN shocks (Bianchi & Schneider 2007). Although variations in the dust emissivity by up to a factor of 10 are theoretically possible (Ossenkopf & Henning 1994), there is currently no observational evidence to suggest the large variations in the dust emissivity which would be required to solve the crisis in this sample. Furthermore, Rowlands et al. (2014) suggest that due to the consistency between observations of gas and dust in individual SMGs and the gas-to-dust ratios implied by the ratio of FIR to CO luminosity, κ does not evolve strongly between low-redshift dusty galaxies and SMGs.

Reproducing SMG properties
In summary, the models which best reproduce the dust-to-stellar mass and gas-to-dust properties of SMGs are LIMS+SNe+grain growth (62 per cent of the sample) and LIMS+maximal SN dust production (explains 58 per cent of dust masses in the SMGs sample, although this is unphysical as it requires all of the SN metal yields to be in the form of dust). The most plausible models are summarized in Table 5. In reality, a mixture of these models will most likely best describe the properties of high-redshift SMGs.
Based on these findings, more complex models were run with dust created by LIMS, SNe and grain growth with moderate dust destruction (m ISM = 100 M ) to reduce the dust-to-metals ratio slightly to better match observations of SMGs. For a closed box model, the average dust-to-stellar mass and gas-to-dust ratios agree well with observed values for SMGs (see Table 5). This shows that a modest amount of dust destruction can be accommodated if dust is produced by both stellar and interstellar sources. In reality, galaxies are unlikely to be closed boxes, we therefore run the same models but with an inflow and outflow with a rate equal to the SFR. Compared to the closed box these models predict slightly lower average final M d /M * (by a factor of 2) and slightly higher gas-todust ratios (by a factor of 2), but these differences are still well within the observational range of these parameters for SMGs.

C O N C L U S I O N S
In this paper, we have used an updated chemical evolution model to reproduce the properties of a submillimetre selected sample of 26 massive, dusty galaxies in the redshift range 1.0 < z < 5.3. Our chemical modelling for the first time utilizes complex SFHs derived from SED fitting to the UV-submillimetre photometry and a detailed treatment of the dust sources and sinks in galaxies.
We can rule out a number of models (Table A1) which result in dust-to-stellar masses and/or gas-to-dust ratios which are inconsistent with observations of SMGs. These models include those with dust produced by LIMS only, and those which have efficient dust destruction (mass of ISM m ISM = 1000 M cleared of dust). The models which best match the observed gas-to-dust ratios include rapid dust build-up from grain growth and SN dust sources. Our main results are as follows.
(i) We find that dust produced only by LIMS falls a factor 240 short of the observed dust masses of SMGs. Adding an extra source of dust from SNe can account for the dust mass in SMGs in only 19 per cent of cases. Even after accounting for dust produced by SNe, the remaining deficit in the dust mass budget suggests that higher SN metal yields, and/or substantial grain growth are required in order for the dust mass predicted by the chemical evolution models to match observations of SMGs.
(ii) Efficient destruction of dust grains by SN shocks (m ISM = 1000 M ) on average decreases the dust mass from LIMS+SNe by a factor of 6-10. Additional sources of dust are required in order to account for the additional shortfall of dust in SMGs caused by dust destruction. Alternatively, dust destruction may be less efficient if dust grains are shielded from SN shocks in dense regions of the ISM. A small amount of dust destruction (m ISM = 100 M ) can be accommodated in our models only if dust is produced efficiently by both stellar and interstellar sources.
(iii) The average metallicity in the closed box model reaches 0.9 Z , which is consistent with the metallicity measured in SMGs. If inflows of pristine gas occur with a rate equal to the SFR the metallicity is reduced to 0.7 Z ; a similar metallicity is reached with enriched gas outflows. Inflows and outflows result in a modest decrease of a factor <1.5 in the dust mass of SMGs. Given the current large range in observed gas-phase metallicities in SMGs, and uncertainties due to possible AGN contamination of the emission lines, we cannot currently distinguish between different inflow and outflow rates. Measurements of gas-phase metallicities which are not affected by the presence of an AGN are required for larger samples of SMGs.
(iv) A top-heavy IMF cannot account for the observed dust masses if dust is produced by LIMS only. With no dust destruction we found that a top-heavy IMF with dust produced by both LIMS and SNe can produce the average dust mass observed in SMGs (within a factor of 2) therefore resolving the dust budget crisis. Yet, given the uncertainties involved (e.g. in the dust destruction rate and metallicity in SMGs), this does not provide unequivocal evidence for a top-heavy IMF in dusty high-redshift galaxies.
(v) Increasing the dust emissivity on average by a factor of 7 can solve the dust budget crisis if dust is produced by LIMS and SNe and is not destroyed by SN shocks. Variations in the dust emissivity are theoretically predicted to be a factor of <3, and, currently there is no observational evidence to suggest a large variation in emissivity occurs in high-redshift SMGs. Finally, an alternative explanation for the dust budget crisis is that the metal yields of stars may be systematically underestimated. This paper has been typeset from a T E X/L A T E X file prepared by the author.