The VANDELS survey: dust attenuation in star-forming galaxies at z = 3-4 Free

Abstract

We present the results of a new study of dust attenuation at redshifts 3 < z < 4 based on a sample of 236 star-forming galaxies from the VANDELS spectroscopic survey. Motivated by results from the First Billion Years (FiBY) simulation project, we argue that the intrinsic spectral energy distributions (SEDs) of star-forming galaxies at these redshifts have a self-similar shape across the mass range 8.2 ≤ log (M_⋆/M_⊙) ≤ 10.6 probed by our sample. Using FiBY data, we construct a set of intrinsic SED templates which incorporate both detailed star formation and chemical abundance histories, and a variety of stellar population synthesis (SPS) model assumptions. With this set of intrinsic SEDs, we present a novel approach for directly recovering the shape and normalization of the dust attenuation curve. We find, across all of the intrinsic templates considered, that the average attenuation curve for star-forming galaxies at z ≃ 3.5 is similar in shape to the commonly adopted Calzetti starburst law, with an average total-to-selective attenuation ratio of R_V = 4.18 ± 0.29. In contrast, we find that an average attenuation curve as steep as the SMC extinction law is strongly disfavoured. We show that the optical attenuation (A_V) versus stellar mass (M_⋆) relation predicted using our method is consistent with recent ALMA observations of galaxies at 2 < z < 3 in the Hubble Ultra Deep Field (HUDF), as well as empirical A_V − M_⋆ relations predicted by a Calzetti-like law. In fact, our results, combined with other literature data, suggest that the A_V–M_⋆ relation does not evolve over the redshift range 0 < z < 5, at least for galaxies with log(M_⋆/M_⊙) ≳ 9.5. Finally, we present tentative evidence which suggests that the attenuation curve may become steeper at lower masses log(M_⋆/M_⊙) ≲ 9.0.

1 INTRODUCTION

Interstellar dust absorbs and scatters the ultraviolet (UV) to near-infrared (near-IR) radiation emitted by stars in galaxies, and re-radiates the absorbed energy in the infrared (IR) and far-infrared (FIR). The observed spectral energy distribution (SED) of a galaxy is therefore strongly affected by the amount, and spatial distribution, of the dust within it. Understanding exactly how dust influences the SEDs of galaxies is crucial to our interpretation of observations.

Characterizing the line-of-sight extinction towards individual stars has been one method used to determine the wavelength-dependent effect of dust in the Milky Way (MW) and other nearby galaxies, referred to as the extinction curve (e.g. Prevot et al. 1984; Cardelli, Clayton & Mathis 1989; Gordon et al. 2003; Hagen et al. 2017). In the case of an extinction curve, we are sensitive to both absorption and scattering out of the line of sight by dust. Uniformly across all curves, the extinction peaks in the UV and falls off with increasing wavelength, a shape which is related to the size distribution and composition of the dust grains (e.g. Draine 2003).

For distant galaxies it is, in general, not possible to measure the extinction of stars along individual sight-lines. Exceptions do exist, for example in the case of gamma-ray burst afterglow spectra (e.g. Zafar et al. 2011; Fynbo et al. 2014; Heintz et al. 2017); however, in the vast majority of cases, we are observing the average effect of dust across an extended region. In this case, the effect of dust is a combination of extinction plus scattering back into the line of sight, which implies a strong dependence on the dust geometry, and is referred to as an attenuation curve (see Calzetti 2001, for a review). Furthermore, since knowing the dust-free SED shape is crucial for constraining the normalization and wavelength dependence of any extinction or attenuation curve, an additional complication in the case of galaxy attenuation is that, unlike for individual stars, the intrinsic shape of a galaxy SED is less well constrained, being a complex function of a number of unknown parameters (e.g. star formation history, stellar metallicity, initial mass function etc.).

Nevertheless, attenuation curves have been derived for both local starburst galaxies (e.g. Calzetti, Kinney & Storchi-Bergmann 1994; Calzetti et al. 2000; Wild et al. 2011; Battisti, Calzetti & Chary 2017a) and, more recently, for star-forming galaxies at higher redshifts (e.g. Reddy et al. 2015, z ≃ 1.5–3). These studies have found, for the most part, that the shape of the attenuation curve is similar to the extinction curve, although, due to the different geometry, the increase towards the UV is not as steep, and attenuation curves are often referred to as being greyer than extinction curves in this respect. A number of studies have attempted to characterize the attenuation curve in terms of physical and structural parameters of galaxies (e.g. inclination, specific star formation rate), and investigated the presence/absence of the |$2175\,\rm {\AA }$| UV ‘bump’ feature common to many extinction curves, but frequently absent in attenuation curves (e.g. Noll et al. 2009; Conroy, Schiminovich & Blanton 2010; Buat et al. 2011; Wild et al. 2011; Buat et al. 2012; Kriek & Conroy 2013; Salmon et al. 2016). Nevertheless, by far the most commonly adopted attenuation law in the galaxy formation research community is the Calzetti starburst law (or Calzetti attenuation curve; Calzetti et al. 2000). For example, our current understanding of the global star formation rate density and stellar-mass density evolution out to z ≃ 8 is generally based on the assumption that the Calzetti law is applicable for all star-forming galaxies at all epochs (e.g. Madau & Dickinson 2014).

A key question in galaxy evolution is whether, or to what extent, the shape of the galaxy attenuation curve changes, either as a function of galaxy properties (e.g. stellar mass), or as a function of redshift. Interestingly, a number of recent results in the literature have suggested that the form of the attenuation curve evolves at z ≳ 2. For example, several recent sub-mm/mm measurements of the infrared excess (IRX) versus UV spectral slope (i.e. f_λ ∝ λ^β) appear to support a scenario in which the attenuation curve becomes much steeper at these redshifts, and more similar in shape to the SMC extinction curve (Capak et al. 2015; Bouwens et al. 2016; Reddy et al. 2017). Indeed, some evolution in attenuation properties may be intuitively expected given the growing consensus that star-forming galaxies at z ≳ 2 are characterized by systematically different physical conditions with respect to their low-redshift counterparts. In particular, in combination with an increase of the specific star formation rates of galaxies with redshift (e.g. Mármol-Queraltó et al. 2016), an evolution in the physical properties of the interstellar medium (ISM) is also observed, including an increase in ionization parameters and decrease in gas-phase metallicities (e.g. Cullen et al. 2014; Kewley et al. 2015; Cullen et al. 2016; Steidel et al. 2016; Strom et al. 2017). Moreover, at early epochs, possible changes in the relative importance of different dust production channels (e.g. supernova, asymptotic giant branch stars, ISM growth) may become important (e.g. Michałowski, Watson & Hjorth 2010; Rowlands et al. 2014; Watson et al. 2015; Popping, Somerville & Galametz 2017). These observations suggest plausible mechanisms (e.g. more extreme radiation fields, different dust production channels) by which the size distribution and composition of dust grains may be different at these epochs.

Nevertheless, there remain a number of uncertainties related to these recent IR measurements which have yet to be fully characterized. First, IR luminosities carry significant systematic uncertainties related to the unknown dust temperatures in high-redshift galaxies (e.g. Schaerer et al. 2015; Faisst et al. 2017), and biases related to the treatment of β may not have been properly accounted for in all cases (McLure et al. 2018). Secondly, other independent sub-mm/mm observations of the IRX-β relation have shown consistency with grey Calzetti-like attenuation curves (e.g. Coppin et al. 2015; Dunlop et al. 2017; Fudamoto et al. 2017; Laporte et al. 2017; McLure et al. 2018). To further compound this confusing picture, similarly contradictory results from simulations exist, with some studies supporting a steep UV attenuation curve (e.g. Mancini et al. 2016) and others finding that Calzetti-like laws are preferable (e.g. Cullen et al. 2017). However, since the majority of the inferences regarding the shape of the attenuation curve at z ≳ 3 remain indirect, one potential avenue of progress is to attempt to directly measure the attenuation curve shape at high redshift.

Predominantly, the method used to derive the attenuation curve is to compare the observed SEDs of a sample of galaxies for which an independent indicator can be used to identify those that are more, or less, affected by dust. In practice, the ratio of hydrogen emission lines (most commonly H α/H β) is used to rank galaxies in increasing order of their Balmer decrement (Calzetti et al. 1994; Reddy et al. 2015). ¹ The primary advantage of this approach is that, in not having to assume the shape of a dust-free SED, potentially large systematic errors are avoided. However, one disadvantage is that since the dust-free intrinsic SED shape is not known, this method is sensitive to the total-to-selective attenuation law, rather than the attenuation law directly. An alternative approach is to explicitly assume the underlying intrinsic SED of all sources. In this case, it is possible to directly probe the wavelength dependence of the attenuation curve by taking the ratio of intrinsic to observed spectra (e.g. Scoville et al. 2015). This method suffers from an increase in systematic errors related to the choice of the intrinsic SED shape. However, these problems will be mitigated somewhat at high redshifts where star-formation histories are predicted to become self similar across the galaxy population (e.g. Finlator, Oppenheimer & Davé 2011).

In this paper, we implement this type of direct analysis to derive the shape of the attenuation curve for a preliminary sample of 236 galaxies at 3 < z < 4, selected from the VANDELS survey (McLure et al. 2017). Crucially, building on the study of Scoville et al. (2015), we construct a set of intrinsic, dust-free, SED templates that contains both simple constant star formation rate models as well as physically motivated SEDs derived from the First Billion Years (FiBY) simulation (Paardekooper, Khochfar & Dalla Vecchia 2015). Furthermore, we utilize a comparable sample size over a much narrower redshift range, mitigating any redshift-dependent systematics. Our aim is to investigate whether we find any evidence to suggest that the average shape of the attenuation curve is evolving with redshift, or whether the commonly adopted attenuation curve of Calzetti et al. (2000) is still, on average, applicable for star-forming galaxies at z ≃ 3.5.

The structure of the paper is as follows. In Section 2, we describe the VANDELS data and the ancillary photometry used in our analysis, and detail the assumptions and methods used to generate the intrinsic galaxy SEDs. In Section 3, we outline our method for deriving the attenuation curve before presenting our results in Section 4, including a investigation of the mass-dependence of both the total attenuation and the shape of the attenuation curve. Finally, we summarize our conclusions in Section 5. We adopt the following cosmological parameters throughout this paper: Ω_m = 0.3, |$\Omega _\Lambda =0.7$|⁠, H₀ = 70 km s⁻¹ Mpc⁻¹. All stellar masses are calculated using a Chabrier (2003) initial mass function (IMF).

2 DATA AND MODELS

In this section, we provide a brief overview of the spectroscopic and imaging data sets used in the current study and define the final sample of star-forming galaxies. We also provide a full description of how the observed SEDs were derived from the rest-frame optical photometry, and outline how the intrinsic template SEDs were generated.

2.1 Vandels survey

VANDELS is a public ESO spectroscopic survey targeting |$\rm {N} \approx 2000$| galaxies at 1.0 < z < 7.0 with the VIMOS spectrograph on the VLT in the CDFS and UDS extragalactic survey fields. Observations have been performed using the medium resolution grism, which provides a spectral resolution of R ≈ 600 over the wavelength range |$4800\hbox{--}10\,000\,\rm {\AA }$| at a dispersion of |$2.5\,\rm {\AA }$| per pixel. Details of the survey and data reduction will be provided in two upcoming papers (McLure et al., in preparation; Pentericci et al., in preparation). The survey is due to be completed in 2018 January, therefore in this work we use the current sample of 1502 fully or partially observed spectra, from which we select 242 galaxies with secure redshifts in the range 3.0 ≤ z ≤ 4.0 (see Section 2.3 for more details). Galaxies at the redshift of our sample are selected to have I_AB ≤ 27.5 and H_AB ≤ 27, with integration times optimized to provide 1D spectra with a signal-to-noise ratio (S/N) of 15–20 per resolution element. Since not all the spectra used this paper have been fully observed, this S/N is not achieved for all galaxies in our sample.

2.2 Photometry

All VANDELS targets have been selected from the CDFS and UDS survey fields. Within each field, |${\simeq } 50\, {\rm per \, cent}$| of galaxies fall within the CANDELS region (Grogin et al. 2011; Koekemoer et al. 2011) and therefore benefit from deep optical to near-IR Hubble Space Telescope (HST) imaging, while the other |$50 \, {\rm per \, cent}$| fall within the wider survey regions and relies on predominantly ground-based imaging. For the remainder of this paper, we reference these four separate subsets by their field (CDFS or UDS) and whether they are covered by primarily space-based (HST) or ground-based (GROUND) imaging (i.e. CDFS-HST, CDFS-GROUND, UDS-HST and UDS-GROUND).

The CDFS-HST photometry utilizes most of the same filters as the publicly available CANDELS catalogue described in Guo et al. (2013), with the difference that the Ks HAWKI photometry has been updated to include the final HUGS survey data (Fontana et al. 2014), and the CTIO U-band filter has been removed. The final catalogue contains photometry across the wavelength range |$0.3\hbox{--}8.0\,\mu \rm {m}$| (i.e. covers the rest-frame UV to near-IR SEDs of our sample out to |$\lambda _{{\rm rest}} \sim 2.0\,\mu \rm {m}$|⁠). Similarly, the UDS-HST photometry is essentially identical to the public CANDELS catalogue published by Galametz et al. (2013), and covers the same wavelength range. The CDFS-GROUND and UDS-GROUND photometry, whilst including some HST filters, are both composed primarily of ground-based observations. Crucially, these regions are not covered by deep Spitzer/IRAC observations, so the observed photometry only extends out the K band (⁠|${\sim } 2.2\,\mu \rm {m}$|⁠) and the rest-frame photometry therefore only covers the UV to optical SED of our sample out to |$\lambda _{{\rm rest}} \sim 0.6\,\mu \rm {m}$|⁠. In the attenuation curve analysis presented in this paper, we confine our analysis to the |$0.1 \lesssim \lambda \lesssim 0.6\,\mu \rm {m}$| region. A detailed description of the various photometric catalogues will be given in the VANDELS survey definition paper (McLure et al., in preparation).

2.3 VANDELS sample

The sample analysed in this paper is drawn from the current batch of 1502 VANDELS spectra. Of these, 480 have formal VANDELS redshifts in the range 3.0 ≤ z ≤ 4.0; however, for this initial study we restrict this sample to only those with the most secure redshift flags (FLAG = 3 or 4, i.e. a |$> 95 \, {\rm per \, cent}$| probability of being correct) to avoid any potential redshift contamination issues (see McLure et al., in preparation). This restricted sample contains 242 galaxies. We then applied a mass cut of log(M_⋆/M_⊙) ≤ 10.6, leaving a total of 236 galaxies. As described in Section 2.5, throughout this paper we compare our observed sample to a sample of simulated galaxies from the FiBY simulation, which only contains galaxies with stellar masses up to log(M_⋆/M_⊙) = 10.6, so we have applied the same upper mass limit to the observed sample. A description of how stellar masses were derived from the observations is given below. Therefore, the final selection criteria for the sample used in this paper is simply the original VANDELS selection criteria (I_AB ≤ 27.5 and H_AB ≤ 27 with a photometric redshift solution in the correct redshift range) plus the requirement of a secure spectroscopic redshift from current observations and an upper mass limit of log(M_⋆/M_⊙) = 10.6. Fig. 1 illustrates how the effect of these criteria is to select galaxies consistent with the star-forming main sequence at 3.0 < z < 4.0 (see also McLure et al., in preparation). The galaxies in our final sample span stellar masses in the range 1.6 × 10⁸ ≤ M_⋆/M_⊙ < 4.0 × 10¹⁰ with a median value of 4.5 × 10⁹ |$\rm {M}_{{\odot }}$|⁠, at a mean redshift of 〈z〉 = 3.49.

2.4 Observed SEDs

To determine the observed SED shape of each galaxy, we used two publicly  available SED-fitting codes: LePhare (e.g. Ilbert et al. 2008) and EAZY (Brammer, van Dokkum & Coppi 2008). EAZY is a redshift fitting code which uses linear combinations of a limited standard template set to best match the observed photometry. The template set is derived from fitting synthetic photometry of a semi-analytic model (SAM) with the PEGASE stellar population synthesis (SPS) library (Fioc & Rocca-Volmerange 1997). These standard templates can be thought of as the ‘principal component’ templates of all SED shapes in the SAM over the redshift range 0 < z ≲ 4 (Brammer et al. 2008). We ran EAZY on all galaxies in our sample using the default parameter set, including emission lines, fixing the redshift to the spectroscopic redshift measured from the VANDELS spectra. A major benefit of EAZY is the relatively small number of a-priori assumptions required; for example is it not necessary to assume a dust attenuation curve or star formation history. Examples of fits to three galaxies are shown in Fig. 2.

On the other hand, LePhare allows one to fit the observed photometry using a range of SEDs built using a chosen SPS model library along with a parametrized star formation history, dust attenuation law, colour excess, and metallicity. The main advantage of this approach is that the stellar mass of the galaxy can be reliably estimated, which we used in building the intrinsic template set, investigating the attenuation versus stellar-mass relation, and in splitting the sample by mass. To run LePahre, we built the SED library using constant star formation BC03 SPS models (Bruzual & Charlot 2003) at three metallicities (0.28, 0.56 and 1.4 Z_⊙)² with 0.0 ≤ E(B−V) ≤ 0.4 assuming the Calzetti et al. (2000) attenuation curve. Again, examples of the resulting fits are shown in Fig. 2. To check that our assumption of a constant star formation history was not biasing the resulting stellar masses, we also ran a set of exponentially rising star formation history models with e-folding times τ = 50–10 000 Myr using the Bayesian SED fitting code bagpipes (Carnall et al. 2017), adopting the same set of SPS models, metallicities and colour excesses. The masses and star formation rates returned from this test are fully consistent with the data in Fig. 1. We note that, given that the shapes of the SEDs are well constrained by the available photometry, the actual differences in the shapes of the resulting best-fitting templates are small and insensitive to model assumptions. Indeed, from Fig. 2 it can be seen that similar SED shapes would be obtained by a simple linear interpolation between the SED points. The median reduced χ² values for the LePhare and EAZY fits are 1.48 and 1.50, respectively; for the final observed SED of each galaxy, we averaged the best-fitting LePhare and EAZY SEDs.

2.5 Intrinsic template SED

As we describe in detail in Section 3, to derive the attenuation curve, the observed SED shape needs to be compared to an estimate of the intrinsic SED shape. In this paper, we make the assumption that all galaxies in the sample have the same intrinsic UV to optical SED shape.

2.5.1 Evidence for a homogeneous SED shape

We were able to assess the validity of this assumption using data from the FiBY cosmological hydrodynamical simulation. A sample of 628 galaxies, with |$8.2 \le \rm {log(M/\mathrm{M}_{{\odot }})} \le 10.6$| (matching the mass range of the observed sample) at z = 4 (the minimum redshift of the simulation), was extracted from two simulation boxes with a combined comoving volume of |${\approx } 3.7 \times 10^4\,\rm {Mpc}^3$| (see Johnson, Dalla Vecchia & Khochfar 2013; Cullen et al. 2017, for simulation details). Fig. 3 shows the normalized mean star formation history and ±1σ scatter for the full sample of simulated galaxies, as well as the mean star formation histories split into three M_⋆ bins. To construct Fig. 3, we first normalized the star formation history of each galaxy using its mean star formation rate, which is the total stellar mass formed at z = 4 divided by the age of the oldest star particle in the galaxy (i.e. an approximation of the formation age of the galaxy). Fig. 3 illustrates how the typical star formation histories are self-similar across the whole population with a relatively small scatter. For example, over a lookback time of ∼300 Myr (from z = 4), where star formation rates are typically 3 × the mean, the scatter is ∼0.15 dex (a factor ∼ 1.4). Furthermore, in the inset panel of Fig. 3, we show the UV-weighted metallicity of the galaxies (an approximation of the mean metallicity of the galaxies formed over the last 100 Myr) as a function of stellar mass. Again, the distribution is relatively tight, with a mean UV-weighted metallicity of 〈Z/Z_⊙〉_Fe = 0.21 ± 0.07. Based on the similarity of the star formation histories, and the narrow range in metallicity, it is reasonable to assume that the underlying SED shape will be similar across all masses. In Section 4.2, we provide further empirical evidence for the validity of this assumption.

It is important to note here that we are defining metallicity as the iron abundance relative to solar (⁠|$[\rm {Fe}/\rm {H}]$|⁠) rather than the total abundance of all metals. This is motivated by the fact that stellar opacity in the UV, and hence the UV spectral shape, is dominated by Fe and is relatively insensitive to |$[\rm {O}/\rm {H}]$| (e.g. Rix et al. 2004). This is important because, for galaxies with rising star formation histories, |$[\rm {Fe}/\rm {H}] \ne [\rm {O}/\rm {H}]$| (e.g. Steidel et al. 2016; Cullen et al. 2017). This can be seen in the FiBY sample, where the mean UV-weighted total metallicity of all 628 simulated galaxies is 〈Z/Z_⊙〉_Tot = 0.76 ± 0.30 (i.e. a factor ≈3.5 larger than the Fe metallicity). As a consequence, since most standard SPS models assume solar abundance ratios, in practice one must choose which of the two abundance definitions to use as the ‘true’ metallicity when constructing SEDs from SPS models. For this paper, we adopt the UV-weighted Fe abundances for constructing the SEDs.

2.5.2 Building the intrinsic SEDs

The simulated galaxy SEDs from FiBY were built as described in Cullen et al. (2017). Briefly, each star particle associated with a given galaxy was assigned an instantaneous starburst SPS model which best matched its age and metallicity, and the final galaxy SED was then constructed by summing up the SEDs of all individual star particles. For the SPS models, we used BPASSv2 (e.g. Eldridge & Stanway 2016; Stanway, Eldridge & Becker 2016), and considered their four fiducial models, which we refer to according to the upper mass cutoff of the IMF and whether or not binary evolution is included. The BPASSv2-100bin models include binary star evolution with an IMF cutoff of 100 M_⊙, and BPASSv2-100 are the equivalent single-star evolution models; similarly, BPASSv2-300bin models include binary evolution with an IMF cutoff of 300 M_⊙ and BPASSv2-300 are the equivalent single-star evolution models. All models have an IMF index of −1.3 between 0.1–0.5 M_⊙ and −2.35 above 0.5 M_⊙. Finally, the nebular continuum contribution is included using cloudy (Ferland et al. 2017) as described in Cullen et al. (2017), assuming maximal nebular contribution (i.e. assuming escape fraction |$f_{{\rm esc}}=0\, {\rm per \, cent}$|⁠).

To generate the set of intrinsic template SEDs used in our analysis we constructed, for each of the four BPASS models, a stack of FiBY SEDs over a similar mass range to the observed sample. To generate each stack, we paired each VANDELS spectrum with the FiBY galaxy, which matched it closest in mass, avoiding duplications. Thus, the mass distribution of the stacked FiBY SED was the closest possible match to the mass distribution of the VANDELS sample. Fig. 4 shows the stacked SEDs for all FiBY models. For a given stack (e.g. FiBY-BPASSv2-100bin), the standard deviation in each pixel, across the wavelengths of interest (⁠|$0.125\,\mu \rm {m} \le \lambda \le 0.63\,\mu \rm {m})$|⁠, ranges from |${\approx } 1\hbox{--}20 \, {\rm per \, cent}$| with an average of |$9\, {\rm per \, cent}$|⁠. We note that, given the similarity in star formation history and metallicity across all masses, this method of mass-matched stacking yielded very similar results to simply stacking all 628 FiBY SEDs.

We also considered two constant star formation rate models from the latest version of Starburst99 (Leitherer et al. 2014) assuming the weaker-wind Geneva tracks without stellar rotation, a Kroupa (2001) IMF and an age of 100 Myr. The two models had metallicities Z_* = 0.002 and Z_* = 0.008, respectively. These models assume single-star evolution and an upper-mass IMF limit of |$100\,\rm {M}_{{\odot }}$|⁠, include nebular-continuum emission, and are referred to as S99-v00-z002 and S99-v00-z008 (see Table 1). The Z_* = 0.002 (⁠|$\rm {\mathrm{Z/Z}_{{\odot }}} = 0.14$|⁠) model was chosen as the closet match to the mean UV-weighted Fe/H based metallicity of the simulated sample (⁠|$\langle \rm {[Fe/H]} \rangle = 0.20$|⁠), while the Z_* = 0.008 (⁠|$\rm {\mathrm{Z/Z}_{{\odot }}} = 0.56$|⁠) was chosen as the closest match to the mean UV-weighted total metallicity (⁠|$\langle \rm {\mathrm{Z/Z}_{{\odot }}} \rangle = 0.70$|⁠). Comparing the two Starburst99 models allowed us to test how sensitive our results are to the assumed metallicity. Finally, an age of 100 Myr was chosen to match the typical UV and optical-weighted ages of the simulated galaxies. The mean ages weighted by |$0.15{}\,\mu \rm {m}$| and |$0.55{}\,\mu \rm {m}$| luminosity are 44 ± 5 Myr and 150 ± 6 Myr for the FiBY sample (averaged across the four BPASSv2 models). An age of 100 Myr is therefore a reasonable compromise between the estimated age of the stellar populations dominating the UV and the optical luminosity at z ≃ 3.5.

Fig. 4 shows the six intrinsic SED templates described above normalized at |$0.55\,\mu \rm {m}$|⁠. A summary of all of the observed and intrinsic template sets is given in Table 1.

3 DERIVING THE ATTENUATION LAW

In this section, we describe the method used to derive the shape and normalization of the UV–optical attenuation curve and assess the validity of the method using a simple simulation

3.1 Fitting method

The first step was to re-sample the observed and intrinsic SED templates on to a common rest-frame wavelength grid. We adopted a wavelength grid spanning |$0.12\hbox{--}0.63\,\mu \rm {m}$| with steps |${\rm d}\lambda =5\times 10^{-4}\,\mu \rm {m}$| (⁠|$5\,\rm {\AA }$|⁠). The upper and lower limits were chosen to mimic the wavelength window used in the derivation of Calzetti et al. (2000) attenuation curve at UV to optical wavelengths. Moreover, |$0.63\,\mu \rm {m}$| is a reasonable upper limit since it is close to the maximum rest-frame wavelength at which all galaxies in our sample have photometric coverage.

The second step was to decide which wavelength pixels to use when performing the fit. Not all pixels are suitable since some will be biased by strong stellar and nebular features. We adopted the following criteria to select good pixels: (i) those with an estimated nebular contribution (line or continuum) of |${<} 10 \, {\rm per \, cent}$|⁠; (ii) those free from stellar absorption/emission features. To implement (ii), we used the Calzetti et al. (1994) windows in the UV, and a set of custom windows at optical wavelengths masking the strong hydrogen absorption features. The resulting rest-wavelength regions that we used for fitting the attenuation curve are listed in Table 2 and shown in Fig. 5. An example fit to one galaxy in our sample is shown in Fig. 6.

3.1.1 The nebular continuum and emission spectrum

We decided to exclude regions of the spectrum with a significant nebular contribution due to the large uncertainties inherent in modelling the nebular emission. In particular, the nebular spectrum is sensitive to the assumed escape fraction (f_esc) of ionizing photons, which depends primarily on the covering fraction of neutral hydrogen and the dust content of a galaxy (e.g. Hayes et al. 2011).

Our nebular modelling is based on the latest determinations of physical parameters within H ii regions at z ≃ 2.5 (e.g. Steidel et al. 2016; Strom et al. 2017) and is described in detail in Cullen et al. (2017). By default we assume an escape fraction of f_esc = 0.0, but we illustrate in Fig. 5 how the resulting intrinsic SED is sensitive to different f_esc assumptions. It can be seen that at certain wavelengths (e.g. |$0.25 \lesssim \lambda \lesssim 0.35 \, \mu \rm {m}$|⁠) the effect of varying f_esc can be significant. Therefore, since the average value of f_esc at z ≃ 3.5 is unknown, with only a handful of individual measurements made to date (e.g. de Barros et al. 2016; Shapley et al. 2016; Vanzella et al. 2016) and, furthermore, variation from galaxy to galaxy might be expected (e.g. Paardekooper, Khochfar & Dalla Vecchia 2013; Paardekooper et al. 2015), the shape of the intrinsic spectrum at these wavelengths is highly uncertain.

Another issue is the potential difference in attenuation between nebular and stellar regions. For example, in the local Universe, and out to z ≃ 2, studies continue to find evidence for differential attenuation, such that the nebular regions suffer a larger amount of attenuation, on average, than the stellar continuum (e.g. Kashino et al. 2013; Price et al. 2014). In this case, equation (2) will not hold for wavelength regions with a large nebular contribution.

Indeed, properly accounting for the nebular emission when fitting the attenuation curve would require a model that incorporates variations in escape fraction, differential attenuation, and perhaps even a difference in shape of the attenuation curve between stellar and nebular regions. We decided, therefore, to restrict our analysis to wavelengths for which the nebular contribution is low. It is important to note, however, that the remaining fitting regions still cover both far-UV and optical wavelengths (Fig. 5), and we therefore remain sensitive to the UV to optical shape of the attenuation curve. Our method can still distinguish, for example, between shallow Calzetti-like laws, or steep SMC-like laws, regardless of the exact details of the nebular modelling/attenuation. We note, however, that by flagging these regions, our method is not sensitive to the presence or absence of the 2175 Å UV bump feature.

3.2 Simulating the fitting method

To assess how accurately the shape and normalization of an input attenuation law could be recovered using this fitting method we performed a simple simulation. Using the FiBY-BPASSv2-100bin SED as the intrinsic template, we constructed a sample of artificially reddened SEDs. First, we generated 1000 values of M_⋆ (equally spaced between |$8.2 \le \rm {log(M_{\star }/M}_{{\odot }}) \le 10.6$|⁠) and assigned a value of A_V using in the A_V–M_⋆ relation for a Calzetti law from McLure et al. (2018). This relation is shown as the dashed black curve in the bottom panel of Fig. 7. For each value of M_⋆ we constructed an attenuated SED in the following way. The purely stellar component of the FiBY-BPASSv2-100bin template was attenuated using a Calzetti law with its corresponding A_V value, while the purely nebular component was also attenuated with a Calzetti law, however this time assuming A_{V, neb} = 2 × A_V. The final simulated SED was taken as the sum of the two components, allowing us to mimic the effect of differential reddening in stellar continuum and nebular regions. Finally, we perturbed the flux values in the SED assuming a |$20\, {\rm per \, cent}$| error in each pixel to mimic the combined uncertainty in the photometry and choice of intrinsic SED shape.

After generating 1000 artificial SEDs, the next step was to recover the attenuation curve of the stellar continuum using the method described above. We first performed the fitting without masking any wavelength regions, the results of which are shown by the blue curves in Fig. 7. We found that the underlying Calzetti law attenuating the stellar continuum could not be recovered in this case because the fit is biased by pixels with a strong nebular contribution, and/or pixels with strong stellar absorption/emission features as described above. The recovered shape of the resulting attenuation curve (A_λ/A_V) is steeper than a Calzetti law, and the A_V–M_⋆ relation is heavily biased to low A_V values (by up to ≃ 0.6 dex at the largest values of M_⋆).

In contrast, by using the wavelength regions given in Table 2, it is possible to recover the input relations (red curves in Fig. 7). There are still small biases in the recovered A_V, which increase with M_⋆, however this bias does not strongly affect the recovered attenuation curve shape (A_λ/A_V). We note that these simulation results hold for all other intrinsic templates, and are also independent of the assumed attenuation law (i.e. similar results are found using an input SMC-like curve). Finally, we found that although it is common to parametrize the attenuation curve using a third-order polynomial (e.g. Calzetti et al. 2000; Reddy et al. 2015), the parametrization of A_λ given by equation (3) works better at recovering A_V and A_λ/A_V given the restricted wavelength range we used in our fitting. As can be seen, this parametrization nevertheless results in a very similar attenuation curve shape to the actual Calzetti law.

In reality, there will of course be a scatter in the value of A_V at a given M_⋆, in the ratio of nebular to stellar attenuation, and even in the shape of the attenuation curve which we have not captured here; however, we argue that this simple simulation illustrates how the underlying shape and normalization of the attenuation curve, as well as the A_V–M_⋆ relation, can be robustly derived using the method outlined in Section 3.

4 Results: the UV–optical attenuation curve

In the section, we present the shape and normalization of the attenuation curve at z ≃ 3.5, derived using the method described above.

4.1 Average shape of the attenuation curve

We first investigated the average shape of the attenuation curve for the 236 galaxies in our VANDELS sample. In Fig. 8, we show the average attenuation curve shapes (A_λ/A_V) corresponding to each of the six intrinsic template SEDs listed in Table 1. For each template, we individually fitted every galaxy in the sample (an example of which is given in Fig. 6); and in Fig. 8, we plot the unweighted mean of all the individual curves along with the ±1σ scatter.

For each of the intrinsic templates, it was not possible to find an individual solution for every galaxy in the sample with that template. This occurred when the UV/optical ratio was larger in the observed SED than in the intrinsic SED (i.e. the observed SED was bluer), implying a negative attenuation was required to map the intrinsic to observed shape. In these cases, either the assumed intrinsic shape, or the observed shape derived from the photometry, is clearly incorrect. For all these cases, we assumed A_V = 0 and excluded these galaxies when averaging. The degree of success of each template at fitting individual galaxies varied, and the number of successful fits is indicated in Fig. 8. It can be seen that for half of the templates the number of unsuccessful fits was a negligible fraction of the sample (⁠|${\lesssim } 5 \, {\rm per \, cent}$|⁠).

It is clear from Fig. 8 that the majority of attenuation curves are shallow and Calzetti-like in shape. Steeper curves are suggested from the FiBY-BPASSv2 single-star models (100, 300), however these are still not, on average, as steep as the SMC extinction law or even the Reddy et al. (2015) curve. Furthermore, the single star models have the lowest success rates for individual fits (⁠|${\approx } 75{\rm -}80\, {\rm per \, cent}$|⁠) indicating that these intrinsic SED shapes, which have the lowest UV/optical ratios (Fig. 4), are perhaps not accurate representations of the intrinsic z ≃ 3.5 population. Both the Starburst99 templates and FiBY-BPASSv2 binary-star models have much higher individual success rates (⁠|${\approx } 90{\rm -}99 \, {\rm per \, cent}$|⁠) and clearly favour a shallower curve similar to the Calzetti et al. (2000) law, or the Charlot & Fall (2000) model. Overall, the simple constant star-formation rate Starburst99 and the FiBY-BPASSv2-300bin templates have the highest success rate for individual fits.

Unfortunately, this analysis alone cannot distinguish which of the intrinsic templates is most likely to be correct one. Nevertheless, it does demonstrate that, irrespective of intrinsic template choice, assumed star formation history, or metallicity, the average shape of the attenuation curve at z ≃ 3.5 is consistent with a grey Calzetti-like law within ±1σ. In other words, the ratio of UV (evaluated at |$\lambda =0.12\,\mu$|m) to optical (⁠|$\lambda =0.55\,\mu$|m) attenuation is consistently within the range A_UV/A_V ≃ 2.5–3.5. Extremely steep attenuation curves, similar in shape to the SMC extinction law (A_UV/A_V ≃ 6.5), are strongly disfavoured. Finally, it is interesting to note that our results are generally incompatible with the direct measurement of the attenuation curve at z ≃ 2 by Reddy et al. (2015), which has A_UV/A_V ≃ 4.0 (dotted line in Fig. 8). While the shape of the Reddy et al. (2015) curve in the UV (⁠|${\lesssim } 0.26\,\mu \rm {m}$|⁠) is consistent with Calzetti (see e.g. Reddy et al. 2015; Cullen et al. 2017), the full UV to optical shape is significantly steeper. The reason for this discrepancy is not immediately obvious, although it is plausibly related to the difference between the two methodologies. Reddy et al. (2015) do not assume an intrinsic SED shape but instead rank galaxies by increasing amount dust attenuation (using the Balmer decrement as a proxy), and use ratios of stacked SEDs, as a function of absolute attenuation, to derive the shape of the attenuation curve. As we discuss in detail below, growing evidence for a relation between stellar mass and attenuation opens up the possibility that a similar analysis to Reddy et al. (2015) can be carried out using stellar mass, rather than Balmer decrement, as a proxy for absolute attenuation. In future, utilizing the full VANDELS data set, it should be possible to directly compare the two different approaches.

4.2 The |$\mathbf {A_V{\rm -}M_{\star }}$| relation

One way to further discriminate between these different models is to investigate the predicted absolute attenuation as a function of stellar mass, which we discuss below. A major strength of our fitting method is that it returns an accurate estimate of A_V for each galaxy (Fig. 7), which can be used to investigate the relationship between absolute attenuation and stellar mass (M_⋆). In Fig. 9, we show the A_V–M_⋆ relation returned from our fitting method for the FiBY-BPASSv2-300bin templates. We compare this result to recent measurements of A_V–M_⋆ at 2 < z < 3 by McLure et al. (2018), based on a deep Atacama Larges Millimeter Array (ALMA) continuum mosaic of the Hubble Ultra Deep Field (HUDF; Dunlop et al. 2017). The McLure et al. (2018) data points shown in Fig. 9 are inferred from measuring the infrared excess (⁠|$\rm {IRX} \equiv \rm {{\it L}_{IR}}/\rm {{\it L}_{UV}}$|⁠) using a stacking analysis. Crucially, these measurements provide an independent comparison sample based on a direct detection of dust emission in the IR.

In Fig. 9, we also show two empirical A_V–M_⋆ relations, one corresponding to the assumption of a Calzetti et al. (2000) law, and the other to an SMC-like attenuation curve. These are derived by McLure et al. (2018) from an analysis of the M_⋆–β relation for a sample of 8407 galaxies at 2 < z < 3 across the stellar mass range 8.25 ≤ log(M_⋆/M_⊙) ≤ 11.50. We refer the reader to their paper for a full description. It is sufficient here to note that the functional form of these relations is dependent on the shape of attenuation curve and the assumed intrinsic UV continuum slope. McLure et al. (2018) use β_int = −2.3 at z ≃ 2.5, which is consistent with the UV continuum slopes of all of our intrinsic template SEDs (−2.41 ≤ β_int ≤ −2.31). As is evident from Fig. 9, McLure et al. (2018) found that their data follow an A_V–M_⋆ consistent with a grey Calzetti-like attenuation law. These data and empirical relations provide an interesting comparison for the template-dependent A_V–M_⋆ relations which result from our analysis.

It can be seen from Fig. 9 that the A_V–M_⋆ prediction for the FiBY-BPASSv2-300bin template is in excellent agreement with both the McLure et al. (2018) data, and the empirical relation for a Calzetti dust law at log(M_⋆/M_⊙) ≳ 9.5. This may seem unsurprising given the fact that the derived attenuation curve shape for this template is so similar to the Calzetti law (Fig. 8). However, it is important to note that this level of consistency is not necessarily guaranteed. First, the McLure et al. (2018) data points are derived from an independent sample of galaxies, at a slightly lower redshift, using a completely different method to the one presented here. Secondly, the empirical relation derived by McLure et al. (2018) depends only on the assumption that the shape of the UV continuum (⁠|$0.12{\rm -}0.25\,\mu \mathrm{m}$|⁠) is constant over the full range of stellar mass, while our method is based on the premise that the intrinsic UV–optical (⁠|$0.12{\rm -}0.63\,\mu \mathrm{m}$|⁠) shape is unchanging. Therefore, the excellent agreement between these separate analyses is both an argument in favour of (i) the Calzetti law being representative of the average attenuation curve at 2 < z < 4 (at least for log(M_⋆/M_⊙) ≳ 9.5), and (ii) our fundamental assumption being correct (i.e. that the intrinsic UV–optical SED shape for star-forming galaxies at these redshifts is relatively constant). Finally, despite only showing the FiBY-BPASSv2-300bin template in Fig. 9, we note that both the FiBY-BPASSv2-100bin and S99-v00-z008 templates give very similar results. In contrast, the S99-v00-z002 template and both the FiBY-BPASSv2 single-templates are incompatible with the McLure et al. (2018) data, as illustrated by the running median A_V–M_* relation for the FiBY-BPASSv2-100 template in Fig. 9.

4.2.1 Redshift dependence of A_V–M_⋆

In Fig. 9, we also show three A_V–M_* relations compiled from various literature sources and presented in McLure et al. (2018). The green curves are taken from Heinis et al. (2013), Pannella et al. (2015), and Álvarez-Márquez et al. (2016) and are based on IR detections of dust emission from star-forming galaxies at 1 < z < 3. McLure et al. (2018) also show how these relations are in remarkably good agreement with attenuation inferred for star-forming galaxies drawn from the SDSS presented by Garn & Best (2010). The consistency between our data and these studies is striking given the variety of different methods involved.

Finally, the square cyan data point shows the average A_V between 9.00 ≤ log(M_⋆/M_⊙) ≤ 9.5 taken from the A_V–M_⋆ relation inferred at z ≃ 5 by Cullen et al. (2017). In Cullen et al. (2017), we derived A_V–M_⋆ by comparing the observed z ≃ 5 luminosity function and colour–magnitude relation with predictions from the FiBY simulation. This stellar-mass range is chosen to cover the masses across which both this study and Cullen et al. (2017) benefit from robust statistics. Again, the consistency is remarkable, given the differences in the two methods. Interestingly, combining all of these literature results with our data suggests both that stellar mass is a good proxy for attenuation, and that the relationship between mass and attenuation, for normal star-forming galaxies with |$\rm {log(M_{\star } / \mathrm{M}_{{\odot }})} \gtrsim 9.5$|⁠, does not evolve between z = 0 to z ≃ 5.

In summary, the A_V–M_⋆ relations predicted by the FiBY-BPASSv2 binary star templates (100bin, 300bin), and S99-v00-z008 template, are fully consistent with observed data at 1 < z < 3. These results support the view that star-forming galaxies at these redshifts are reddened, on average, by a shallow attenuation law similar to the Calzetti starburst law. The consistency of our results with data from previous studies at different redshifts, using different techniques, is also remarkable. Current data appear to support a scenario in which the A_V–M_⋆ relation for star-forming galaxies does not evolve significantly from the local Universe out to z ≃ 5. Deeper IR data are needed in order to make robust statements at |$\rm {log(M_{\star } / \mathrm{M}_{{\odot }})} < 9.5$|⁠.

4.3 Parameterization of the attenuation curve

4.4 Attenuation curve at low masses

There is a suggestion from Fig. 9, both from our data and the measurements of McLure et al. (2018), that the average attenuation curve may be steeper than the Calzetti law at the lowest masses, particularly below log(M_⋆/M_⊙) ≃ 9.0. The idea that the attenuation curve could be mass dependent has been suggested by some studies in the literature, which argue that the attenuation curve steepens as the optical depth decreases (e.g. Seon & Draine 2016; Leja et al. 2017), which could be the result of a mass-dependent change in the dust geometry in galaxies (e.g. Paardekooper et al. 2015). On the other hand, the opposite trend has been suggested in other papers (e.g. Zeimann et al. 2015).

In Fig. 11, we show the average attenuation curve shapes for the 33 galaxies in our sample with log(M_⋆/M_⊙) < 9.0. To improve the accuracy of this fitting, we constructed stacked intrinsic spectra from FiBY with masses in the same mass range; however, the results are similar using the original stacks. It can be seen that we find some evidence for a slight steepening using our fiducial FiBY-BPASSv2-300bin template. The steepening is more pronounced for the S99-v00-z008 template. However, although we don't show the ±1σ scatter of the curves in Fig. 11 for clarity, both remain consistent with the Calzetti attenuation law within 1σ, and rule out and SMC-like curve at >1σ. Therefore, though we find some tentative evidence for a steepening of the attenuation curve at log(M_⋆/M_⊙) < 9.0, a larger sample of galaxies at these masses is needed to confirm this result. Moreover, we still find no evidence to support an average attenuation curve as steep as the SMC extinction law, even at these lower masses.

4.5 Evidence for a |$2175\,\rm {\AA }$| UV bump?

Finally, we investigated whether there was any evidence for excess attenuation in our sample caused by the |$2175\,\rm {\AA }$| absorption feature (often referred to as the UV bump; Stecher 1965). The |$2175\,\rm {\AA }$| absorption feature is observed in the MW and Large Magellanic Cloud (LMC) extinction curves, and has recently been shown to be present in the extinction curves along the majority of sight-lines to the SMC (Hagen et al. 2017). Intriguingly, the feature is not present in the Calzetti starburst attenuation law, although a number of other studies, both in the local Universe and up to z ∼ 2, have found evidence for its presence (e.g. Noll et al. 2009; Conroy et al. 2010; Buat et al. 2011; Wild et al. 2011; Kriek & Conroy 2013; Battisti, Calzetti & Chary 2017b). The absence/presence of a UV bump, and its dependence on other physical parameters, is a potentially useful constraint on the grain properties and physical conditions within star-forming galaxies (e.g. Seon & Draine 2016).

As discussed in Section 3, our method for deriving the shape of the UV to optical attenuation curve explicitly excluded the UV bump wavelength region. However, the availability of deep rest-frame UV spectra for the galaxies in our sample provides an alternative route to determine whether the average attenuation curve at z > 3 contains a |$2175\,\rm {\AA }$| UV bump feature. The feature is not sampled by the VIMOS spectra for galaxies at z > 3.5, therefore we could only utilize the 122 galaxies at 3.0 < z < 3.5 for this analysis (with 〈z〉 = 3.288). A stacked composite UV spectrum was formed by shifting each spectrum to the rest-frame, re-sampling on to a common wavelength grid and taking the median at each dispersion point after rejecting 5σ outliers. To correct for redshift differences, the flux of each individual spectrum was scaled to the flux that would be observed at the mean redshift of the sample. The 1σ error spectrum was estimated by bootstrap re-sampling the flux values at each dispersion point.

The two intrinsic templates were attenuated using equation (11) for various values of E_b and compared to the composite spectrum after normalizing to the median flux in the wavelength range |$0.17 < \lambda < 0.18\mu \rm {m}$| (Fig. 12). The data robustly rule out the presence of a strong UV bump with E_b > 1.0. Interestingly, assuming a weak |$2175\,\rm {\AA }$| UV bump feature with E_b ≃ 1 results in a better match to the observed composite spectrum than a pure Calzetti law (E_b = 0). This result is consistent with the K-band selected sample of galaxies at 0.5 < z < 2.5 from Kriek & Conroy (2013), who find that the strength of the UV bump is correlated with the steepness of the attenuation curve, and that for Calzetti-like curves E_b = 0.85 ± 0.09. However, given the uncertainties on the composite spectrum, we consider this only marginal evidence for a UV bump. Nevertheless, the comparison does constrain the average strength of the UV bump to be, on average, relatively weak at z ≃ 3 (E_b ≲ 1). For comparison, the UV bump in the MW extinction curve has E_b ≃ 4. In a future work, with the full VANDELS sample, we will present a thorough analysis of the presence/absence of a UV bump at 2.5 < z < 3.5, and its dependence on various galaxy properties.

5 CONCLUSIONS

We have presented an investigation into the attenuation curve of galaxies at z = 3–4 using an initial sample of 236 star-forming galaxies from the VANDELS survey. By comparing the observed shape of the galaxy SEDs to a set of intrinsic SED templates we have been able to derive the shape and normalization of the attenuation curve based on the assumption that all galaxies across the mass range |$8.2 < \rm {log(M_{\star } / \mathrm{M}_{{\odot }})} < 10.6$| have the same underlying intrinsic SED shape. Using this method, we investigated the average attenuation curve at z ≃ 3.5 and its mass dependence, as well as the relation between absolute attenuation and stellar mass. Our main findings can be summarized as follows:

• From an analysis of the star formation histories and metallicities of a sample of 628 simulated galaxies extracted from the FiBY simulation at z = 4, we argue that the intrinsic shape of the UV–optical SED of star-forming galaxies at these redshifts is approximately constant across the full stellar-mass range of our observed sample. We construct a set of six intrinsic SED templates, both using physically motivated star formation histories and metallicities from FiBY, as well as simple Starburst99 constant star formation rate models.

• Based on the assumption of a constant intrinsic SED shape, we outline a method for fitting the shape and normalization of the attenuation curve for any individual observed galaxy SED. We demonstrate the feasibility of this method with a simple simulation. Applying this method to the full sample, we find that the average attenuation curve shapes derived from each of the intrinsic templates are consistent with a grey Calzetti-like attenuation law within ±1σ, and find no evidence for a steep SMC-like attenuation law.

• By comparing the A_V–M_⋆ relations predicted by each template to recent observations at 2 < z < 3 from McLure et al. (2018), we find that a subset of the intrinsic template sets yield results which are fully consistent with the observed data, and with the A_V–M_⋆ relation predicted for a Calzetti-like attenuation law. Again, we find no evidence for an A_V–M_⋆ relation consistent with the SMC extinction law. This result, combined with other literature data, suggest that the relationship between stellar mass and A_V does not evolve over the redshift range 0 < z < 5. Furthermore, across all these redshifts, we confirm previous results that suggest stellar mass is a good proxy for absolute attenuation.

• Using our fiducial intrinsic template (FiBY-BPASSv2-300bin), we transform our inferred attenuation curve (A_λ) into the more standard total-to-selective attenuation curve formulation (k_λ). We find a total-to-selective attenuation ratio of R_V = 4.18 ± 0.29, consistent with the original Calzetti et al. (2000) value of R_V = 4.05 ± 0.80.

• We present tentative evidence for steeper attenuation curve shapes at masses log(M_⋆/M_⊙) ≲ 9.0, in line with models which predict a dependence of the attenuation curve shape on the total optical depth. Better statistics are required to confirm this results, or alternatively deeper ALMA continuum imaging at these masses.

• Finally, using a stacked composite spectrum of 122 galaxies at 3.0 < z < 3.5, we find marginal evidence for a weak |$2175\,\rm {\AA }$| UV bump in the average attenuation curve at z ≃ 3, but robustly rule out the presence of a strong UV bump feature. A more detailed examination will be possible in future using the full VANDELS spectroscopic sample.

Acknowledgements

The authors would like to thank the anonymous referee whose comments and suggestions significantly improved the final version of this manuscript. FC, RJM, SK and JSD acknowledge the support of the UK Science and Technology Facilities Council. ACC acknowledges the grants ASI n.I/023/12/0 ‘Attività relative alla fase B2/C per la missione Euclid’ and MIUR PRIN 2015 ‘Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid’. The Cosmic Dawn centre is funded by the DNRF. This research made use of Astropy, a community-developed core python package for Astronomy (Astropy Collaboration et al. 2013), NumPy and SciPy (Oliphant 2007), Matplotlib (Hunter 2007), IPython (Pérez & Granger 2007), and NASA's Astrophysics Data System Bibliographic Services.

Footnotes

REFERENCES