David Carton
email: david.carton@univ-lyon1.fr

Jarle Brinchmann

Thierry Contini

Benoit Epinat

Hayley Finley

Johan Richard

Vera Patricio

Joop Schaye

Themiya Nanayakkara

Peter M Weilbacher

Benoît Epinat

Vera Patrício

Lutz Wisotzki

First gas

Keywords: galaxies, evolution -galaxies, abundances -galaxies, ISM 1 INTRODUCTION

come

redistributed within galaxies. Equally, we can also study how galaxies accrete and lose their gas.

In the classical inside-out picture of galaxy evolution, the inner regions of galaxies formed first from low angular momentum gas. And with the increase of angular momentum over time, the radial scale-length of star-formation has progressed outwards in galaxies (Larson 1976). Inside-out growth can explain why the centre of the Milky Way is more chemically evolved (metal-rich) than its outskirts (Portinari & Chiosi 1999). Moreover it can also explain why exponentiallydeclining radial metallicity profiles are ubiquitous in isolated and relatively massive ( ∼ > 10 8 M ) low-redshift galaxies (e.g. Vila-Costas & Edmunds 1992;Zaritsky et al. 1994, and references therein).

Interestingly, not only do all low-redshift star-forming galaxies present negative metallicity gradients, they also present similar slopes (when the metallicity gradient is normalized to the size of the galaxy). For example, Sánchez et al. (2014) and Ho et al. (2015) find the 1σ scatter in the metallicity gradients between galaxies is approximately the same magnitude as the mean of the metallicity gradient.

That said, while the radial metallicity profile of the Milky Way appears to be well-described simply by an exponentiallydeclining function (Esteban et al. 2017), there is some debate as to whether this is really true for all galaxies. Indeed, as noted by Sánchez et al. (2014) and many others (e.g. Rosales-Ortega et al. 2011;Bresolin et al. 2012;Marino et al. 2016), a significant fraction of galaxies exhibit shallow/flat metallicity profiles in their innermost and/or outermost regions. More specifically, in a recent study of 102 massive ( ∼ > 10 10 M ) galaxies, Sánchez-Menguiano et al. (2018) report that only 54% galaxies can be described by a single metallicity gradient, with the remainder exhibiting flattening either within ∼ < 0.5 r e (half an effective radius) and/or beyond ∼ > 1.5 r e . However, while this appears to be the prevailing picture, it should be noted that Carton et al. (2015) report that the metallicity profile of some galaxies may actually steepen in the outer (gas dominated) regions.

Nevertheless, although it appears that in some situations there are deviations from an exponential metallicity profile, it appears also that (apart from the innermost regions) most galaxies show similar exponential metallicity profiles. At higher redshift, the assumption of exponential metallicity profiles is hard to test (due to the limited spatial resolution), but it does seem to give a reasonable description of the data.

However, it is intriguing that the common (similar slope) metallicity gradient, which is seen in low redshift galaxies, is not observed at higher redshift (z ∼ > 0.6). Studies have found that not only was average metallicity gradient previously flatter than today, but there was also a large amount of scatter in the observed metallicity gradients (Stott et al. 2014;Wuyts et al. 2016).

Indeed, most striking is the fact that many high-redshift galaxies have positive (inverted) metallicity gradients (e.g. Queyrel et al. 2012). Galaxies with centres more metal poor than their outskirts are rarely, if ever, observed in the present-day Universe. The prevailing explanation for this phenomenon is that metal-poor gas is flowing (or has flowed) into the inner regions of these galaxies. The inflowing gas dilutes the metals, suppressing the metallicity. The acquisition of extra gas is subsequently expected to trigger intense star formation in the galaxy. In this regard Stott et al. (2014) identified a weak trend for elevated star-formation rates in the galaxies with flatter and inverted metallicity gradients.

There are two mechanisms that have been proposed to cause the inflow of metal-poor gas: galaxy-galaxy interactions and cold flows. Firstly, galaxy-galaxy interactions might trigger radial flows within a galaxy's disc, transporting metal-poor gas from the outskirts to the inner regions. At low-redshift there is observational support for this idea. Indeed, while it is true that there is a common metallicity gradient in isolated galaxies, it has been found that non-isolated (interacting) galaxies possess significantly flatter metallicity gradients (Rich et al. 2012). Furthermore, this mechanism (where mergers flatten metallicity gradients) has been demonstrated in numerical simulations (Rupke et al. 2010;Torrey et al. 2012). It appears, however, that galaxy-galaxy interactions are merely capable of flattening the metallicity gradient of galaxies, but not inverting it.

On the other hand, cold flows (the other mechanism proposed for producing inverted gradients), may be more successful (Cresci et al. 2010). These flows are cold streams of gas which can penetrate through a galaxy's hot halo to reach the galaxy itself (Kereš et al. 2005;Dekel & Birnboim 2006). However, it is not clear that this metal poor gas would flow directly to the centre of the galaxy, as would be needed to invert the metallicity gradient. Indeed it has been suggested that these streams may build an extended gas disc (Stewart et al. 2011;Danovich et al. 2015), which itself could, in turn, contract to form a compact star-forming clump at the centre of the galaxy (Dekel & Burkert 2014;Zolotov et al. 2015). There is some observational support for this, with the tentative identification of such cold-flow discs in a couple of high redshift galaxies (Bouché et al. 2013(Bouché et al. , 2016)).

Simulations predict that these hypothesized cold flows would dominate the gas supply of galaxies in the early Universe (z ∼ > 1.8) but we should, however, expect them to become increasingly rarer at late times (e.g. van de Voort et al. 2011;Woods et al. 2014). So, while cold flows may explain why galaxies observed at z ≈ 3.4 present inverted metallicity gradients (Troncoso et al. 2014), it is harder to invoke cold flows as an explanation for inverted gradients observed at z ≈ 1.

To summarize briefly, there is a disparity between metallicity gradients in the high-redshift and low-redshift Universe. While not necessarily contradictory, the high-redshift results point to stochastic processes dominating galaxy evolution, whilst the low-redshift results suggest a secular evolution of galaxies. There are few, if any, observations of metallicity gradients in galaxies between 0.1 ∼ < z ∼ < 0.6. Clearly, bridging this gap is a necessary step towards understanding the disparity between the high-and low-redshift results. Here, with the Multi Unit Spectroscopic Explorer (MUSE; [START_REF] Bacon | Ground-based and Airborne Instrumentation for Astronomy III[END_REF]Bacon et al. , 2017)), we will provide for the first time a large sample of metallicity gradients in intermediate-redshift galaxies (0.08 < z < 0.84).

Analysing these observations presents several challenges. The first challenge is to correct for the effects of seeing on our data. As demonstrated by both Yuan et al. (2013) and Mast et al. (2014), failing to correct for seeing effects will produce systemically flatter metallicity gradients. But the challenge posed by seeing is not unique to our work, and accordingly other recent studies have used simulated observations to apply an a posteriori correction to infer the true metallicity gradient. In Carton et al. (2017, hereafter C17), however, we presented an alternative forward-modelling approach. This method is better able to quantify the degeneracies that arise from seeing-corrupted data, and therefore yields formal estimates for uncertainty in the recovered metallicity gradient.

A second challenge we face is that we derive metallicities from a combination of nebular emission-lines. Depending on a galaxy's redshift, different emission-lines fall within the wavelength range of a spectrograph. It is a well-documented issue that different metallicity calibrations (especially when using different emission-lines) produce different results (e.g. Kewley & Ellison 2008). With a forward-modelling approach we can overcome these limitations and thereby self-consistently infer metallicity gradients independently of redshift.

With our observations of intermediate-redshift galaxies we will attempt to reconcile the high-and low-redshift pictures of galaxy evolution, particularly with respect to the gas supply in these systems. We structure the paper as follows. In Section 2 we describe our observations and outline our galaxy sample selection. We detail our methodology in Section 3, where we also include a sensitivity analysis for our model. Section 4 is dedicated to presenting the results on the derived metallicity gradients. In Section 5 we provide a discussion of these results. Finally we conclude our findings in Section 6.

Throughout the paper we assume a ΛCDM cosmology with H 0 = 70 km s -1 Mpc -1 , Ω m = 0.3 and Ω Λ = 0.7.

DATA

We wish to spatially resolve metallicity gradients in distant (0.1 ∼ < z ∼ < 0.8) galaxies. Using integral-field spectroscopy (IFS) we can map the nebular emission from star forming regions in these galaxies, and therefore measure radial metallicity variations.

Here we will use observations taken with the MUSE instrument situated at UT4 of the Very Large Telescope (VLT). We will construct our galaxy sample by combining data from both Guaranteed Time Observations (GTO) programmes and commissioning activities. However, because of the differing observing strategies employed in these observing campaigns, our data sample is rather inhomogeneous. Galaxies were observed with a variety of integration times (between 1 -31 h) and in a variety of seeing conditions. We will describe these datasets fully in Section 2.1.2.

MUSE Observations

Instrument Description

MUSE is an integral-field spectrograph that employs an image slicing technique at optical wavelengths. In normal (non-extended) wide-field mode, MUSE provides spectra over a continuous 1 × 1 Field of View (FoV) with continuous spectral coverage (4750 Å -9300 Å). These spectra have a wavelength resolution of ≈ 2.5 Å full-width half-maximum (FWHM), although this is not entirely constant with wavelength (see Bacon et al. 2017). The spatial sampling of the data is 0.2 × 0.2 , but in actuality the spatial resolution of our data limited by the seeing.

Field Description

Given MUSE's large contiguous field, we do not target individual galaxies, rather we target collections of galaxies with limited pre-selection. While each field was chosen to optimize the scientific objectives of the different observing programmes, in general the galaxy selection is essentially blind. There is one exception where one field (CGR30) targets a galaxy group at z ≈ 0.7. But even then there will be foreground and background galaxies in this field that are blindly selected. In total there are 35 MUSE pointings, covering ∼ 35 arcmin 2 on the sky.

While the parent galaxy sample selection is essentially blind, there are no straightforward criteria for selecting the galaxies for which we can measure metallicity gradients. That said, we would expect that we can measure metallicity gradients in the largest and brightest galaxies at a given redshift. We therefore do not impose a priori selection criteria, and analyse all MUSE detected galaxies that have known redshifts z < 0.9, rejecting those with insufficient signal-to-noise (S/N). We describe this S/N cut in Section 2.3.1. We present a post hoc description of the final sample in Section 2.3.

We will now outline the data used in our analysis as follows (a summary is displayed in Table 1):

Hubble Deep Field South (HDFS): As one of the commissioning activities of MUSE, a single deep (26.5 h = 53 × 1800 s) field in the HDFS was acquired. The average seeing conditions were good (FWHM = 0.66 at 7000 Å). [START_REF] Bacon | [END_REF] present a full description of the data. Here we use a slightly improved data reduction to the one presented therein. This new reduction (version 1.24) includes improvements to the sky subtraction and slice normalization (quasi flat-fielding), see Bacon et al. (2017). Hubble Ultra Deep Field (UDF): The MUSE-Deep GTO survey has observed a 9 field mosaic that covers the UDF. This 3 × 3 field has been observed to a depth of ≈ 10 h (in exposures of 1500 s each). In addition, there is also an extra-deep 1 × 1 portion of the mosaic that reaches ≈ 31 h. During the observations the average seeing conditions were good, resulting in a final combined PSF with FWHM = 0.61 -0.67 at 7000 Å. This data will appear in Bacon et al. (2017) and the full catalogue, including redshifts and line fluxes etc., in Inami et al. (2017). Here we are using a slightly older data reduction (version 0.31) than that presented in the aforementioned papers. The subsequent improvements in the reduction primary focus on removing low-level systematics, however, since we concern ourselves here with only the brightest galaxies these improvements are by and large irrelevant for our science.

Chandra Deep Field South (CDFS):

The MUSE-Wide GTO program is surveying a portion of the Chandra Deep Field South (amongst other fields). In the end this will produce a 60 tile mosaic of the CDFS at 1 h depth (using exposures of 4 × 900 s). Here we will use only the first 24 fields that have been observed. These observations were performed in moderate and poor seeing conditions, resulting in a FWHM = 0.7 -1.1 at 7000 Å) This dataset (version 1.0) is described by Herenz et al. (2017) and in a forthcoming data release (Urrutia et al., in prep.). COSMOS Group 30 (CGR30): A third GTO program is surveying galaxies in group environments. In our analysis here we will use observations of one of these galaxy groups Table 1. Summary of galaxy observations. The final sample of galaxies was obtained from various targeted fields (with differing exposure depths and seeing conditions). We list the number of galaxies obtained from each field.

Field

Depth PSF FWHM # of galaxies (namely Group 30 as identified in the zCOSMOS 20k Group Catalogue Knobel et al. (2012)). The deepest portion of the field covers slightly less than the full 1 × 1 FoV and reaches a 9.75 h depth (39 × 900 s). However, one galaxy in our final sample (COSMOS2015 ID: 506103) lies in a shallow "snapshot" region and therefore was only observed to 1 h depth. The average seeing conditions were good (FWHM = 0.60 at 7000 Å). This field (version 1.0) was presented in Epinat et al. (2018).

Data Reduction

Above we described fields from four different observing programmes and as is to be expected, there are differences in the specifics for each of the data reductions. However, all reductions use the standard data reduction pipeline (Weilbacher, in prep.) 2 to produce calibrated datacubes. In all fields sky subtraction is performed using the Zurich Atmospheric Purge (ZAP; Soto et al. 2016) software, which employs a principal component analysis technique developed specifically for MUSE data.

The largest difference between the reductions are the implementations (or lack thereof) of slice normalization. Slice3 normalization improves the uniformity (flatness) of the field and is primarily required because the flat-field calibrations are not taken at the exact same time as the science exposures. Small changes in the instrument alignment due to thermal variations can alter the throughput to the slits. Slice normalization is essentially a secondary flat-fielding, that selfcalibrates using the individual science exposures. Because multiple exposures are combined, these semi-random slice systematics contribute to the effective noise in the final datacube. As a result, the application of slice normalization is very important for faint galaxies, but will have little impact on the bright galaxies that we study here. Thus the fact that the various data reductions implement the normalization differently will not affect our results.

The final datacubes are constructed with equal sized voxels4 (0.2 × 0.2 × 1.25 Å), mirroring the native pixel size at the charge-coupled device (CCD) level. With this voxel size, the typical seeing-limited point-spread function (PSF) is well sampled, whilst the 2.5 Å line-spread function (LSF) is only just critically sampled.

PSF Determination

A critical part of our analysis is to forward model the seeing effects on our data. It is therefore necessary to measure the final PSF directly from our datacubes. MUSE commissioning activities established that the PSF is relatively spatially invariant across the FoV; a PSF model fit to a bright star somewhere within the FoV can be applied across the whole field. Unfortunately, not all fields contain such bright stars, and therefore we use a variety of PSF determination techniques in the different fields:

HDFS: This field contains a bright star to which [START_REF] Bacon | [END_REF] fit a Moffat function. The FWHM of the Moffat profile is allowed to vary as a function of wavelength, but the Moffat β parameter is not. We describe the FWHM as a piecewise linear function with three knots {4750 Å, 7000 Å, 9300 Å}. UDF: Most of the MUSE UDF pointings do not contain any bright stars, in which case the PSF must be inferred from non-point source objects (i.e. galaxies). Hubble Space Telescope (HST) images are convolved with a Moffat function and fit5 to MUSE pseudo-broadband images, see Bacon et al. (2017). We obtain a best-fit Moffat profile as a linear function of wavelength. The Moffat β parameter is assumed to be constant. The accuracy of this method has been verified by comparing the results in those fields that do contain bright stars. CDFS: As with the UDF, many of the CDFS fields are also devoid of bright stars. For these fields, Herenz et al. (2017) determined the PSF using bright, compact galaxies. Choosing only galaxies with minimal substructure, they parametrized the HST images of these galaxies as 2D elliptical Gaussian functions. By convolving the galaxy models with a circular Gaussian PSF, they then derived the PSF that best matches a series of MUSE pseudo-broadband images (constructed at various wavelengths). From this they were able to construct a PSF model for each field, where the Gaussian FWHM varied linearly as a function of wavelength. There were, however, some fields that also contained a bright star, enabling the PSF to be additionally derived by a direct fit to this point source. In these fields, and where the FWHM was a closer match to the telescope autoguider measurement, they adopted the stellar PSF fit in place of that fit to the galaxies. CGR30: This field contains four relatively faint stars. We perform a simultaneous fit to all stars using a Moffat PSF. The FWHM is assumed to have a 3 rd order polynomial dependence as a function of wavelength, whilst for the Moffat β parameter we allow only a linear dependence with wavelength.

Of the four targeted areas, determining a reliable PSF model for CDFS fields was the most challenging. Thus, to check the validity of our procedure, we applied the same method as used for the UDF to the CDFS data; we found that both methods yielded similar results.

To summarize, we model the PSF as an axisymmetric function (either a Moffat or Gaussian function). The FWHM is free to vary as function of wavelength, and is typically found to be ∼ 25% larger in the blue than in the red. These wavelength-dependent PSF models are directly used in the forward-modelling of our observations.

Derived global properties

As part of our analysis, we study metallicity gradients as a function of global galaxy properties, e.g. stellar mass, star formation rate, and disc size. We now outline how these quantities are derived.

Stellar Mass

There exists extensive broadband photometry for all of the fields that we study here. Fall (2000) dust model, and the absorbed energy is self-consistently reradiated in the infrared. The photometry used in each of the four fields is derived from various sources, all of which are approximations of the total magnitude: HDFS: For this field we use the four-band HST photometry {F300W, F450W, F606W, F814W} from Casertano et al. (2000). UDF: Extensive deep multi-band HST photometry is provied by Rafelski et al. (2015) in {F225W, F336W, F435W, F606W, F775W, F850LP, F105W, F125W, F140W, F160W}. Where possible, we use all filters. CDFS: Here we use the photometric catalogue of Guo et al. (2013) using exclusively the HST photometry {F606W, F775W, F814W, F850LP, F105W, F125W, F160W}, using all where available. CGR30: In this field we adopt the photometric catalogue of Capak et al. (2007) using {Subaru B j V j g + r + i + z + NB816, SDSS u g r i z, CFHT u * i * , HST F814W, CTIO/KPNO K s }.

To test the sensitivity of our results to the choice of stellar mass estimates, we also derive them using fast (Kriek et al. 2009), assuming both an exponentially declining starformation history and adopting a Calzetti et al. (2000) dust law. And although we identified some small differences between the FAST and MAGPHYS stellar mass estimates ( ∼ < 0.2 dex), we find that neither our results nor conclusions are impacted by our choice of stellar mass estimate.

Star Formation Rate (SFR)

We derive global star formation rates directly from the MUSE data, taking the spectrum integrated across the whole galaxy.

On this we perform a full spectral-fitting using platefit (Tremonti et al. 2004;Brinchmann et al. 2004). We will describe the spectral fitting in Section 3.1.2. Here it suffices to say that we obtain the Hα, Hβ and Hγ emission-line fluxes, accounting for the underlying stellar absorption.

For low-redshift galaxies (z ∼ < 0.4) we use Hα and Hβ to compute the SFR. At higher redshifts Hα is redshifted beyond the MUSE wavelength range, so we compute the SFR in these galaxies using Hβ and Hγ instead.

To correct for dust, we adopt the Charlot & Fall (2000) birth-cloud absorption curve which attenuates the luminosity, L(λ ) as follows

L ext (λ ) = L(λ )e -τ(λ ) , (1)
with

τ(λ ) = τ V λ 5500 Å -1.3 , ( 2
)
where τ V is the V-band optical depth. Depending on the redshift we use either the observed Hα/Hβ or Hγ/Hβ ratios to calculate τ V . For this we assume intrinsic Case B Balmer recombination ratios of j Hα / j Hβ = 2.86 and j Hγ / j Hβ = 0.468. These values are appropriate for H ii regions of temperatures, T e = 10, 000 K, and electron densities, n e = 100 cm -3 [START_REF] Dopita | Astrophysics of the diffuse universe Dopita[END_REF]. Finally, we convert the inferred dust-corrected Hα luminosities to SFRs using a scaling relation between Hα and SFR log 10 SFR M yr -1 = log 10 L(Hα) erg s -1 -41.27, (3) (Murphy et al. 2011;Hao et al. 2011;Kennicutt & Evans 2012). This assumes a Kroupa IMF (Kroupa 2001), and is therefore effectively consistent with the Chabrier IMF we adopted when obtaining stellar masses. For z ∼ > 0.4 galaxies (where Hα is not available), we assume L(Hα) = 2.86 L(Hβ ) (i.e. the Case B ratio).

Galaxy Morphology

In C17 we presented a method for modelling the metallicity gradients in our galaxies. As inputs, this method requires four basic morphological parameters: the galaxy centre (Right Ascension, RA, and Declination, Dec.), the inclination of the galaxy (inc.) and the position angle of the major axis on the sky (PA). For our discussion we also need galaxy size, which we shall express as the exponential disc scale-length, r d .

All the fields we present here are well studied and have existing morphological catalogues.

HDFS: For this field Contini et al. (2016) provide a detailed morphological analysis of the resolved galaxies. They perform a bulge-disc decomposition on the HST F814W imaging, which yields all the necessary morphology information. UDF + CDFS: For both these fields van der Wel et al.

(2012) provide a catalogue of single Sérsic fits. While this catalogue provides most of the relevant information, it does not provide estimates for the galaxy inclinations. The catalogue only provides the axis ratios of the galaxies. Since axis ratio is primarily a function of inclination, we can convert axis ratios into inclinations. However, axis ratios have a secondary dependence on the morphological type of the galaxy. Since the Sérsic index is a proxy for morphological type, we can reduce this secondary dependence by partitioning the catalogue into quintile bins of the Sérsic index. Within each of the these bins we order the galaxies by decreasing axis-ratio and use the rank order of the galaxy within the bin to estimate its inclination. For this we assume that we view galaxies from a random orientation on a 3D sphere, and therefore the probability of observing a galaxy with an given inclination, i, is P (i) ∝ sin (i). In other words we are more likely to observe a galaxy edge-on that we are to see it exactly face-on. To estimate the galaxy sizes we use the half-light radii, r e , reported in the catalogue. To convert these to disc scalelengths we simply assume that the galaxy profile is a bulgeless exponential disc (i.e. r d ≈ 0.596 r e ). Note that as a result of this assumption, we will likely underestimate the true disc size. In fact, because the most massive galaxies tend to have more prominent bulges, this might create an unwelcome trend such that the sizes of the most massive galaxies are more significantly underestimated. Fortunately, to a great extent the correlation between stellar mass and bulge fraction is driven by the large number of passive (and bulgy) galaxies at high masses (Bluck et al. 2014;Huertas-Company et al. 2015). And since we will be selecting galaxies with bright nebular emission, we preferentially select against such bulgy galaxies. As a result, we do not consider a potential mass trend to be of great concern, neither do we anticipate extremely bulge dominated systems. The catalogue provides the morphology derived from three HST bands (F105W, F125W and F160W). For any given galaxy we use the morphology of the band with the highest S/N. Finally, note that while there are other optical morphological catalogues available, we choose to use the van der Wel et al. (2012) catalogue because the near infrared data is less influenced by recent star formation activity and therefore should provide a better description for the morphology of the overall stellar population. CGR30: For this field we use the morphological assessment provided by the COSMOS 2005 Morphology Catalogue, which uses morpheus (Abraham et al. 2007) to measure the morphological parameters. As above, this catalogue also only reports a galaxy's axis ratio, not inclination. We again apply the rank ordering method to convert to axis ratio into an inclination. To avoid mixing morphological types, we use the concentration index6 as a proxy for morphological type. We divide the catalogue into decile bins of the concentration index, and perform the rank ordering within each. This catalogue provides the half-light radii of the galaxies, which we convert to exponential disc scale-lengths as above.

The different galaxy morphology catalogues are measured in different photometric bands. This will systematically affect the measured galaxy sizes; sizes of late-type galaxies measured at redder rest-frame wavelengths will systematically appear smaller. We correct for this using the parametrization of van der Wel et al. ( 2014) (their equations 1 & 2). This correction depends only on the galaxy's redshift and stellar mass. In this work all galaxy sizes are quoted as if they were measured at a rest-frame wavelength of 5000 Å.

Sample Description

As mentioned previously, we do not make an a priori selection for our sample. Nevertheless, there are many galaxies for which we cannot meaningfully constrain the metallicity gradient. Naturally we would expect that we can only measure metallicity gradients in the largest and brightest galaxies. However, it is non-trivial to map this into a clean cut on global properties (e.g. stellar mass, size and SFR). Therefore we build our final sample based upon data-driven criteria (i.e. S/N).

Selection Criteria

We extract spatially resolved emission-line fluxes from our parent sample of 590 MUSE detected galaxies (and with MUSE redshift measurements z < 0.85). The procedure for this extraction is described in Section 3.1. In many galaxies we fail to detect any emission-line component.

To meaningfully constrain metallicity (and distinguish its effects from dust) we need two strong forbidden lines and two Balmer lines to be detected at S/N ≥ 5. Exactly which emission lines are chosen depends on the galaxy. Using the globally integrated spectrum we choose the two forbidden lines with the highest S/N, and the two Balmer lines with the highest S/N. A typical choice for a low-redshift galaxy (z ∼ < 0.4) is {Hβ , [O iii]5007, Hα, [S ii]6717,6731}. And a typical choice at higher redshift is

{[O ii]3726,3729, Hγ, Hβ , [O iii]5007}.
Since we need to constrain the metallicity gradient of the galaxy, not just its metallicity, the line emission must be detected in multiple spatial bins. Explicitly we require that the four chosen emission-lines are all detected at S/N ≥ 5 in at least three spatial bins. (The spatial binning scheme is described in Section 3.1.1.)

Overall these criteria amount to a minimum S/N cut. How this selection maps into galaxy properties depends on the field (because our observations have different depths and seeing conditions).

In addition to this main selection cut, we apply three further criteria:

• In rare occasions the emission we detect may not be associated with the target of interest, i.e. the data is contaminated by a brighter neighbouring galaxy, at the same redshift. These cases can be identified through visual inspection. We manually excluded galaxies where significant contamination is apparent.

• When modelling the metallicity gradients in our galaxies we assume the galaxies to be infinitesimally thin discs. This approximation is acceptable for face-on galaxies, but, however, becomes increasingly questionable for more inclined systems. We therefore exclude galaxies with an estimated inclination inc. > 70 • .

• Galaxies with bright active galactic nuclei (AGN) will produce bright line-emission. Such emission would alter the observed emission-line ratios, and thus alter the inferred metallicities. In the following section we explain how we exclude AGN.

AGN Exclusion

It is difficult to determine the metallicity of galaxies when they are contaminated by emission from AGN or lowionization nuclear emission-line regions (LINERs). We could treat such galaxies in one of three ways. Firstly, one could attempt to model the flux contribution from a compact central source. Secondly, we could mask out the central regions of a galaxy, and derive the metallicity gradient from the outer regions of the galaxies. A third approach, and the one we adopt here, is to simply discard galaxies from our sample if they appear to have a significant AGN/LINER component. We classify our galaxies using standard emission-line ratios classifications, which we apply to a galaxy's globally integrated spectrum.

At low redshift (z ∼ < 0.4) we use both the [N ii]/Hα and [O iii]/Hβ ratios to classify galaxies. We follow the classification scheme of Brinchmann et al. (2004) to divide galaxies into three categories: pure star forming, those with significant AGN, and those that fall in between (i.e. intermediate). We exclude galaxies classified as AGN galaxies, but we do not automatically exclude the intermediate case galaxies. These intermediate galaxies are inspected manually and we exclude those that clearly possess broad emission-line velocity components, indicative of AGN galaxies.

At z ∼ > 0.4 the [N ii]/Hα is redshifted out of MUSE's wavelength range. We therefore adopt the Mass-Excitation (MEx) diagnostic [START_REF] Juneau | [END_REF] to classify the galaxies into the same three classifications (star forming, intermediate, AGN). As before, intermediate cases are manually inspected.

In Fig. 1 we show where our galaxies lie with respect to the two diagnostics. Out of 87 galaxies we exclude one galaxy below z ≈ 0.4 and two galaxies above.

Contamination from evolved stars

In some galaxies hot low-mass evolved stars (HOLMES) can contribute a significant number of ionizing photons, giving rise to LINER-like line emission (Stasińska et al. 2008). And, therefore, these stars present another possible source contamination which could affect our metallicity measurements.

In the following section (2.3.4) we shall demonstrate that our galaxy selection criteria ultimately result in a preferential selection of galaxies with high specific star formation rates. As a result, it is a priori unlikely that, on the whole, the evolved population would dominate over the young as the main source of line emission. However, if the star formation is not uniform, it is entirely plausible that local regions within the galaxies could be dominated by evolved stars.

Fortunately, Cid Fernandes et al. (2011) have proposed a simple cut based on the equivalent width (EW) of Hα that allows one to classify galaxies with emission dominated by the evolved stellar population; galaxies with EW (Hα) > -3 Å are to be considered dominated by evolved stars7 . Unfortunately, we are not able to observe Hα for the highest redshift galaxies in our sample (z ∼ > 0.4). However, all else being equal, and assuming a case B ratio of L(Hα) ∼ 3 L(Hβ ), we can define a new cut EW (Hβ ) > -1 Å.

In our analysis (Section 3.1) we measure the line emission in our galaxies divided into 2445 spatial bins. Of these, only 54 (∼ 2%) of the bins have EW (Hβ ) > -1 Å. Note that there is also a large uncertainly associated with these measurements (S/N (EW (Hβ )) < 2 for all 54 bins), implying that these bins are relatively faint.

Thus we are able to conclude that the potential impact of line emission from evolved stars is quite limited, and should not affect the metallicities we derive.

Sample Properties

In Fig. 2 we present the global properties of our final sample of 84 galaxies. Therein we show the distributions of stellar mass, SFR, rest-frame B-V colour and main-sequence offset (∆SFR)8 . In addition, for comparison, we also plot the parent sample of all galaxies with a MUSE detected redshift. This includes all 590 galaxies, even those that do not meet our selection criteria (Section 2.3.1).

It can be seen that our final sample preferentially selects the more massive and more strongly star forming galaxies. There is also a clear redshift dependence such that, in panel (a), the low-mass galaxies are almost exclusively lowredshift galaxies. The same is true for the SFR (panel b), where the effect appears even stronger.

In contrast, both B-V colour and ∆SFR show different trends. In panel (c) we see that our final sample is fairly representative of the parent sample. We note that the galaxy redshifts are relatively evenly distributed between each bin.

The same is true for the main-sequence offset parameter (panel d) where, above ∆SFR ∼ > 0 dex, the final sample traces the shape as the parent sample and the redshifts are evenly distributed. We can display this another way; in Fig. 3 we show the mass-SFR correlation for our galaxies. At high redshift our galaxies all lie above or on the main sequence. And at low redshift a large fraction of galaxies are found below the main sequence. However, note that the main sequence derived by Whitaker et al. (2012) uses ultraviolet and infrared photometry to calculate their SFRs, whereas we use emission lines to derive our SFRs. As a result the trend with redshift may not be purely a selection effect, since we could expect some systematic offsets and trends when comparing our SFRs relative to our adopted main sequence.

To summarize, at all epochs we are finding blue galaxies that lie on the upper half of the main sequence, however, at high redshift we are biased towards massive, strongly star-forming galaxies.

ANALYSIS

Many of the galaxies we observe are heavily corrupted by seeing. Furthermore, we must also aggregate (or "bin") spaxels 9 together to increase the S/N of our data, creating additional 8 David Carton et al. resolution loss. So, to recover the intrinsic metallicity gradient in our galaxies, we must model both the effect of seeing and binning on our data.

In C17 we demonstrated such a method for inferring both the central metallicity, log 10 Z 0 , and metallicity gradients, ∇ r (log 10 Z) in distant galaxies. While the central metallicity and metallicity gradients inferred this way are inevitably model dependent, we performed a series of tests (including mock observations) to validate our procedure. In this paper we apply this method to the sample of galaxies observed with MUSE. We have made a few minor modifications to the method presented in C17. For brevity here we only outline the method and highlight the changes.

In the next section we will first explain how we extract the emission-line fluxes (Section 3.1). We will then proceed to describe the fitting of our model to the data (Section 3.2). Finally, we shall investigate how sensitive our recovered model parameters are to particular model inputs (Section 3.3).

Emission line flux extraction

When extracting emission-line fluxes to measure metallicity gradients there is a trade off between the number of spatial bins and the S/N of the data within each bin. One must choose a S/N threshold that is sufficiently high to minimize systematic errors in the emission-line measurements, whilst avoiding losing too much spatial information.

Spatial Binning

In C17 (appendix C) we designed a binning algorithm that attempts to maximize the number of spatial bins above a S/N threshold. In essence this algorithm performs successive passes over the data with an (initially small) fixed bin size. When it is longer possible to find spatial bins above the S/N threshold, the bin size is then increased and the process repeated. This continues up to a maximum allowed bin size, at which point we terminate the procedure and assign any remaining unbinned spaxels to a nearby bin.

When calculating the S/N of a spatial bin, we perform a full spectral fitting to the coadded spectrum. For a successful bin, we require set emission lines to be detected at S/N ≥ 5, where the set of tested emission lines is chosen on an objectby-object basis. This set typically consists of four lines, the two highest-S/N Balmer lines and two highest-S/N forbidden lines.

There are some peculiarities that arise from our binning strategy. To preserve as much radial information as possible we define our spatial bins in polar coordinates (in a plane inclined to the observer). However, a direct consequence of working in a non-Cartesian coordinate system is that the pixels within a bin are not necessarily all close to each other in Cartesian space. An additional oddity is caused by the successive passes with increasing bin size. This can sometimes result in bins that are (partially or entirely) enclosed within another. While neither of these effects are ideal, we remind the reader that we mirror the exact binning in the model (i.e. we give the binning segmentation map as a model input). As a result, we view the strangeness of our binning scheme to be a fair trade for optimizing S/N whilst still preserving radial resolution.

Finally we note that, because we do not impose a minimum bin size, our spatial bins can be much smaller than the PSF. Consequently emission-line fluxes of adjacent bins are not strictly statistically independent. Since we do not incorporate such covariances in the model likelihood function, this might lead us to underestimate the errors in the derived model parameters. However, in C17 (section 3.1) we showed that our model nevertheless provides reasonable error estimates. We therefore do not expect this issue to be of much concern.

Spectral Fitting

To extract the emission-line fluxes from the MUSE spectra we use the platefit spectral fitting code (Tremonti et al. 2004;Brinchmann et al. 2004). platefit applies a two-step process that first fits the stellar continuum (with emission lines masked) before fitting the nebular emission-line component (with the best-fit continuum subtracted). Note that the procedure we employ here is identical to that presented in C17 (section 4.2). We summarize it briefly here. The continuum fitting step of platefit does not fit either the redshift or velocity dispersion of the spectrum.

These two parameters must be provided in advance, and we do so as follows. The redshift of the spectrum is obtained using autoz (Baldry et al. 2014). We wish this to be robust, so if the value determined by autoz deviates by more than ±500 km s -1 from our initial redshift guess, then we default to that initial value. To estimate the velocity dispersion we use vdispfit10 . At low S/N, however, vdispfit can yield outliers beyond the realistic range [10 -300] km s -1 . If values outside this range are produced, we adopt a default value of 80 km s -1 .

With the values of the redshift and velocity dispersion predetermined, the stellar continuum is fit using a combination of Bruzual & Charlot (2003) SPS model templates. These templates form a grid of ten ages and four metallicites, spanning [5.2 Myr, 11 Gyr] and [0.004, 0.04], respectively. 10 David Carton et al. The SPS models are built with the MILES stellar libray (Sánchez-Blázquez et al. 2006), and therefore have a similar spectral resolution to that of our MUSE data (2.5 Å). This resolution is sufficient for the purpose of obtaining a good stellar continuum subtraction. If the stellar continuum is too faint / non-existent, the continuum fitting routine can fail.

In these cases we approximate the continuum by applying a running median filter to the spectrum (width 150 Å).

In the second platefit step, the best-fit continuum is subtracted from the observed spectrum. The emission lines are modelled as Gaussian functions. The velocity offset and velocity dispersions are the same for all emission lines. However, unlike the continuum fitting, these two velocity components are free parameters and need not be specified in advance.

The emission-line fluxes are determined from the spectral fitting. However, the formal emission-line flux errors are typically underestimated (see Brinchmann et al. 2013). We rescale these formal errors using a S/N-dependent correction to obtain better estimates of the true flux error. The correction factors were determined from duplicate Sloan Digital Sky Survey (SDSS; York et al. 2000) observations. Note that the SDSS spectra are of resolution comparable to MUSE, and the platefit correction factors are thus still applicable.

Inferring metallicities

We use a forward modelling method to correct for seeing effects allowing us to derive the true central metallicity and metallicity gradient of a galaxy. With this paper, we are also publicly releasing the code. 11

In C17 (section 2) we explained our method in detail 11 https://bitbucket.org/cartondj/metaldisc and described how it can be used to fit the observed 2D emission-line flux distribution of galaxies. Our method is almost identical to that presented therein, so we shall only briefly outline it here and highlight the changes. We approximate a galaxy as an infinitesimally thin disc, inclined to the observer. The disc is described by four fixed morphological parameters (RA, Dec, inc., PA). Our model contains five free parameters: total SFR of the galaxy, SFR tot , central metallicity, log 10 Z 0 , metallicity gradient ∇ r (log 10 Z), ionization parameter at solar metallicity, log 10 U , and the V-band optical depth, τ V .

The metallicity profile in the galaxy is assumed to be axisymmetric and is described by an exponential function

log 10 Z(r) = ∇ r (log 10 Z) r + log 10 Z 0 , ( 4
)
where r is the radius. This parametrization is commonly used in high redshift studies (e.g. Stott et al. 2014;Wuyts et al. 2016, amongst many others). However, as we discussed in the introduction, low redshift studies have suggested that the metallicity profile may be flat at both small ( ∼ < 0.5 r e ) and large ( ∼ > 1.5 r e ) radii (see Sánchez-Menguiano et al. 2018). That said, given the mass dependence they report for the inner flattening, and the fact that our galaxies are generally less massive than those presented therein, it is not obvious that we should expect inner flattening in our galaxies. As for the outer flattening, we simply do not expect to be sensitive to the low S/N outskirts of our galaxies. Therefore, in concordance with other high-redshift studies, we believe it is justified to adopt the simplistic parametrization of equation 4. So note that, although we do not attempt to prove this to be the exact metallicity profiles of our galaxies, it is still a useful description for the overall variation in metallicity across a galaxy.

As per C17, we also assume the ionization parameter to be anti-correlated with metallicity

log 10 U (Z) = -0.8 log 10 Z/Z + log 10 U , ( 5
)
where Z is the solar abundance and log 10 U is the ionization parameter at the solar abundance.

To predict the observed emission-line ratios we use the photoionization models of Dopita et al. (2013, hereafter D13), who tabulate emission-line fluxes over a grid of metallicities and ionization parameters. These models assume that most elemental abundances vary linearly with the oxygen abundance, except for carbon and nitrogen, whose dependence is empirically calibrated to observations of H ii regions.

At each spatial position we interpolate the D13 model grids to the appropriate values of metallicity and ionization parameter. The modelled emission-line ratios only depend on the radial coordinate; there is no azimuthal dependence.

We wish to include the Hδ and Hε emission lines in our model fit. However, these Balmer lines are not provided by the D13 photoionization models. To include these lines we need to extend the D13 photoionization models. We do this by tabulating the Case B recombination ratios j Hδ / j Hβ and j Hε / j Hβ as a function of j Hγ / j Hβ . By interpolating these at the D13 photoionization model values of L Hγ /L Hβ , we assign the appropriate L Hδ and L Hε for each photoionization model. The Case B recombination ratios were determined with PyNeb (Luridiana et al. 2015) using atomic data from Storey & Hummer (1995).

SFR Maps

To model the emission-line luminosities we need to model the SFR distribution in the galaxy. For this one could adopt a parametric model, c.f. Wuyts et al. (2016) who assume an exponential disc (a disc where SFR declines exponentially with increasing radius). However, since we have high-resolution HST imaging for all our galaxies, we prefer to relax this assumption and provide a 2D SFR map. Unlike the emissionline ratios, the modelled emission-line luminosities do have an azimuthal dependence.

We assume that the distribution of stellar light in the HST images provides a rough approximation for the relative SFR distribution. In the MUSE-Deep UDF fields we use the deep HST F775W imaging. In all other fields we use the F814W band. Although ideally we would like to use photometry that covers the rest-frame ultraviolet light in each galaxy, the necessary deep multi-wavelength HST imaging does not exist for all fields. To strike a good balance we use the F775W and F814W filters, since they providing the S/N required. This has the additional advantage that, because these filters have similar pivot wavelengths, we also ensure a degree of consistency between the fields.

We construct SFR maps of our galaxies by cropping the HST images to only include flux within an ellipse of radius 4× r d along the major axis. This ellipse has the same morphology (RA, Dec, inc., PA) as above. We inspect each image and alter the mask if necessary to ensure that we include all flux from the object and to remove other objects or defects. When necessary we interpolate over these rather than mask them. Negative flux values are set to zero. The final result is an SFR map that represents the relative spatial distribution of the SFR. The absolute SFR values are determined by normalizing the map to the total SFR, SFR tot , which is a free parameter in the model.

From experience with our C17 work, we found that it was quite important to include the complex spatial variations of the SFR in our model. Using a SFR map (as opposed to simply assuming an exponential disc) appeared to improve the accuracy of our model.

However, including a SFR map does not correct for all effects of clumpy star-formation. Indeed, a limitation of our model is that, by construction, we do not accommodate for galaxies with strong azimuthal metallicity variations12 ; we should anticipate that we may observe galaxies where the bright star-forming clumps have uncharacteristically low/high metallicities.

However, it should be noted that, for the low-redshift galaxies where azimuthal metallicity variations have been quantified, the variations are typically found to be quite small ( ∼ < 0.05 dex; e.g. Sánchez et al. 2015;Zinchenko et al. 2016;Sánchez-Menguiano et al. 2017), and there is only the odd example of a galaxy with variations as large as ≈ 0.2 dex (e.g. Vogt et al. 2017).

That said, even if there were large azimuthal variations in our galaxies, it does not mean that we cannot still make inferences in these cases because, in effect, we likely still mea- sure the radial metallicity profile of brightest star-forming regions (i.e. a quasi flux-weighted metallicity profile). Nevertheless, this would be in contrast with low-redshift observations, where we typically consider the measured metallicity profile to be describing the metallicity of the gas as a whole. This distinction is an important one to bear in mind when interpreting metallicity gradients reported by all studies of poorly resolved high-redshift galaxies.

Model Fitting

With this galaxy model we are able to mimic the resolution loss due to the seeing and spatial binning. Thus for every spatial bin we can generate a set of model fluxes that can be compared to those observed. We fit emission lines that are observed at S/N ≥ 513 . For clarity, we emphasize that this threshold is applied for all detected emission-lines, not just the four chosen in Section 2.3.1. In other words, some emission lines may only be detected in a few bins, but a critical subset (two Balmer, two forbidden) will be detected in all bins.

To fit our model to the data we use the MultiNest multi-modal nested sampling algorithm (Feroz et al. 2009;Feroz & Hobson 2008;Feroz et al. 2013) that we access through a Python wrapper (Buchner et al. 2014). With MultiNest we can calculate the posterior probability distributions (posteriors) of our five model parameters.

The prior probability distributions (priors) that we place on our model parameters are outlined in Table 2. Except for two differences, these priors are mostly similar to the priors that were adopted in C17.

One difference is that we now adopt a logarithmic prior on SFR tot , where previously (for technical reasons) we had adopted a linear prior of [0, 100] M yr -1 . For a normalization parameter we believe that a logarithmic prior is more appropriate. Changing this prior should have little effect on the derived metallicity profiles, provided that the dust attention, τ V , is well constrained by the data.

A second and more significant change is that we now adopt a narrower prior on the metallicity gradient, ∇ r (log 10 Z) (previously we had adopted [-0.5, 0.5] dex/kpc). In C17 (appendix B) we found that the inferred metallicity gradients could be sensitive to the choice of prior. We attributed this sensitivity to the finite range of metallicities spanned by the photoionization model grid. Metallicity profiles of galaxies with low central metallicities and steep negative metallicity gradients will bottom-out at the lowest metallicity allowed by the photoionization model grid. As a direct result, negative metallicity gradients become somewhat degenerate in models with low central metallicities. Likewise, positive metallicity gradients become degenerate in models with high central metallicities.

In C17 we found that allowing an unnecessarily broad range of metallicity gradients could skew our inferences towards extreme values of the metallicity gradient. Hence we have now chosen to limit our metallicity gradient prior to the range [-0.3, 0.3] dex/kpc. This spans the range of currently measured metallicity gradients (e.g. Marino et al. 2016;Wuyts et al. 2016;Leethochawalit et al. 2016;Belfiore et al. 2017;Sánchez-Menguiano et al. 2018). In the following section we will assess the sensitivity of our results to various parameters, including the width of this prior, and present a way to flag galaxies where the aforementioned degeneracy may be of concern.

Sensitivity analysis

With our model we perform a fit with five free parameters. However, in our modelling we make particular choices, and fix certain parameters. We should attempt to assess the sensitivity of our results to these, considering both the magnitude of systematic errors we could expect, and the significance of unaccounted random uncertainties.

Here we will investigate four possible sources of additional errors: the width of the metallicity gradient prior, a misestimation of the PSF, a miscalculation of the galaxy's inclination, and a misidentification of the galaxy centre. To estimate the size of these effects we re-analyse our galaxies after altering one of the parameters (e.g. increasing the FWHM of the PSF by 10%). By comparing the difference between this alternate model and the fiducial model, we can assess the impact of unaccounted systematics.

Sensitivity to the width of the metallicity gradient prior

First we wish to identify the galaxies that are sensitive to the choice in metallicity gradient prior. In Fig. 4(a) we compare the differences between two flat priors with different widths: [-0.3, 0.3] (our fiducial prior) and [-0.5, 0.5] dex/kpc (the alternate prior). We see that most galaxies lie on the 1:1 line, indicating that they are robust against the choice of prior. However, we note that some of the galaxies with steep metallicity gradients deviate significantly, showing a high degree of sensitivity to the prior. This sensitivity to the prior indicates an unreliable estimate of the metallicity gradient. These galaxies should be treated with care in any analysis. Indeed, such behaviour might indicate that we have insufficient resolution (due to seeing and/or spatial binning effects) to constrain the metallicity profile of these galaxies.

To assess this sensitivity we need a metric to quantify the difference between the fiducial and alternate models. For this we use the Jensen-Shannon distance applied to the metallicity gradient posteriors derived from the two models

JS dist = 1 2 ∑ i P i ln P i M i + 1 2 ∑ i Q i ln Q i M i , (6)
with

M i = P i + Q i 2 , (7)
where P i and Q i are the discretized posteriors of the fiducial and alternate models, respectively. We place an arbitrary cut on galaxies with large Jensen-Shannon distances (JS dist > 0.3). From Fig. 4(a) this can be seen to be a good identifier of galaxies deviating from the 1:1 line. In the remainder of this work we consider the derived metallicity gradients for these 19 galaxies to be suspect, and we will flag them as potentially unreliable. For transparency we do not discard them and, unless otherwise stated, we include them in our statistics.

Sensitivity to PSF errors

The assumed PSF is perhaps one of the greatest sources of unaccounted error in our analysis. In our model we assume to know the PSF perfectly, however, the true PSF is somewhat uncertain, particularly in the fields that do not contain bright stars. Uncertainties in the PSF lead to additional random errors in our result. More importantly, if we incorrectly parametrize the PSF, we could systematically bias our inferred metallicity gradients. While it is challenging to truly characterize the effects of PSF errors, we can nevertheless attempt to estimate the magnitude of the problem. We re-analyse our galaxies adopting two different PSFs, one where the FWHM is systematically 10% smaller than the fiducial model, and another where it is 10% larger. We show the latter case in Fig. 4(b) and, as is to be expected, if we overestimate the PSF width then we will systematically overcorrect for seeing effects and infer systematically steeper metallicity gradients.

There is also noticeable scatter; a moderate uncertainty in the PSF will displace some galaxies significantly from the 1:1 line. In general, (although not shown here) we find that the smallest galaxies are the most sensitive. We can use the vertical displacement from the 1:1 line as an indication of the expected additional uncertainty. We calculate the standard deviation (weighted by the inverse of the mean 1 σ error), and find the additional uncertainty to be ≈ 0.026 dex/kpc, after excluding the unreliable galaxies14 . Naturally this value is only indicative and should not be treated as exact. Not least because it will vary from galaxy to galaxy with the smaller galaxies being more affected.

Sensitivity to galaxy inclination errors

Another possible source of uncertainty is the inclination of the galaxy. Whilst it is not so likely that we systematically miscalculate the inclination, there is certainly some random error. In general the inclination is most uncertain for faceon galaxies, however, this is counteracted by the fact that the metallicity gradients derived from face-on galaxies are probably the most robust against errors in the inclination.

As a first guess we consider inclination uncertainties . the errors of the two models are not independent and are likely to be highly correlated, consequently the errors appear much larger than the scatter. Each data point is coloured by the Jensen-Shannon distance between the 1D marginal posteriors of the fiducial and alternate models, highlighting the datapoints that deviate most from the 1:1 relation (shown as a red dotted line). Squares indicate a subset of galaxies which have large Jensen-Shannon distances (JS dist > 0.3) when comparing the fiducial [-0.3, 0.3] and the [-0.5, 0.5] metallicity gradient priors. We deem these galaxies as unreliable (see Section 3.3.1). In each panel a solid black line shows the best fit to the data, excluding unreliable galaxies and weighted by the inverse of the formal data errors. A red shaded region indicates the (weighted) vertical scatter about this 1:1 relation, again excluding the unreliable galaxies. Equations for the best fit and scatter are given in the bottom-right corner of each panel.

Table 3. Estimated additional uncertainty in the metallicity gradient. These are calculated for both with the unreliable galaxies excluded and also with them included. Again we stress the approximate nature of these results and there may be other unaccounted sources of uncertainty. of ±10 • and construct two alternate models. We show one of these models in Fig. 4(c). This plot indicates that if we overestimate the inclination, we will derive systematically smaller metallicity gradients. However, we note that most galaxies have small Jensen-Shannon distances and clearly lie on the 1:1 line (with minimal scatter). This suggests that the best-fit is in fact dominated by a few outliers, and therefore conclude that, in general, our model is robust against inclination uncertainties. As before, from the vertical scatter we estimate the additional uncertainty to be ≈ 0.016 dex/kpc.

Source of

Sensitivity to galaxy centre misidentification

Finally, we address the impact of misidentifying the centre of the galaxy (the point about which the radial metallicity profile is defined). We re-analyse the data, but shift the galaxy centre by 0.06 along the direction of the galaxy's major axis. This distance is approximately one tenth the size of the MUSE PSF. However, we remind the reader that the galaxy centre is actually defined from higher-resolution HST imaging. While an absolute 0.06 shift will impact the smallest galaxies in our sample the most, it is arguably easier to define the centre of small compact galaxies than it is to define the centre of a large irregular galaxies.

If galaxies have a radial metallicity profile, we might naïvely expect that shifting away from the true centre would result in a flatter metallicity gradient. However, in Fig. 4(d) we see no systematic trend towards flatter metallicity gradients. Misidentifying the galaxy centre appears to add no systematic bias, but it could add moderate scatter to the inferred metallicity gradients. We estimate this additional uncertainty to be ≈ 0.017 dex/kpc.

Total additional uncertainty

In summary we have considered three sources of additional uncertainty that are not incorporated by the model: misestimation of the PSF, miscalculation of the inclination and misidentification of the galaxy centre. While, these are not necessarily the only contributions to the total additional uncertainty, they are certainly likely candidates.

In Table 3 we summarize our findings and add sources of uncertainty in quadrature to roughly estimate the total additional uncertainty. When calculating the additional uncertainties we excluded the unreliable galaxies, so for full transparency we also present the estimated additional uncertainties if we also include the unreliable galaxies. When considering additional uncertainties induced by misestimation of the PSF and misidentification of the galaxy centre, we note that by excluding the unreliable galaxies we decrease the additional uncertainty by about 23% and 48%, respectively. Indeed it can be seen that some of the more discrepant galaxies in Fig. 4(b &d) are also those flagged as unreliable (as defined via the Jensen-Shannon distance metric, see equation 6). Thus it appears that a simple cut on galaxies that are sensitive to the prior, also partly helps to remove galaxies which are sensitive to other aspects of the model input.

RESULTS

In this section we study the metallicity gradients of our galaxies and search for trends with galaxy size, mass and SFR. Our main findings are as follows:

• The average galaxy in our sample has a negative metallicity gradient. However, there is considerable scatter and a few galaxies have positive metallicity gradients.

• If we select only the largest galaxies (r d > 3 kpc), the scatter reduces. There are, in fact, no large galaxies with inverted metallicity gradients.

• We do not find significant trends between a galaxy's metallicity gradient and either its mass or its SFR. The lack of correlation with the SFR is in contrast with other studies (as we shall discuss here and in §5.2).

From our MUSE observations we have constructed a sample of 84 galaxies between 0.08 < z < 0.84. At a 2σ significance level we identify 32 galaxies with negative metallicity gradients and 7 with positive gradients. We classify 26 galaxies as having metallicity gradients that are not inconsistent with zero. Hereafter we shall term these galaxies as having flat metallicity gradients, a name which should not be overinterpreted; we cannot (and do not) claim galaxies with flat metallicity gradients to have no metallicity gradient whatsoever, but rather our data/analysis lacks the statistical power to say otherwise. A further 19 galaxies are flagged as unreliable for being sensitive to the metallicity gradient prior (Section 3.3.1). In Table A1 we provide the inferred central metallicities and metallicity gradients of our galaxies. Additionally, we provide an image atlas of the galaxies in Fig. A1. (As electronic-only supplementary material we provide plots showing both the model fit quality and parameter constraints.

Examples of each are given in Figs. A2 & A3, respectively)

Our first result is that on average our galaxies have a negative metallicity gradient with a median of -0.039 +0.007 -0.009 dex/kpc (using bootstrap and Monte Carlo resampling). We also apply a Wilcoxon signed-rank test and reject the null hypothesis that the average metallicity gradient is zero (p = 0.0014). Contrast this with the findings of Wuyts et al. (2016) who, in a study of 0.6 ∼ < z ∼ < 2.6 galaxies, find the average metallicity gradient to be flat.

In Fig. 5 we present the metallicity gradients of galaxies as a function of their size. Panels (a & b) show the metallicity gradient expressed in scaled units (dex/r d ) and physical units (dex/kpc), respectively. Throughout the remainder of this work we will present metallicity gradients in scaled units that normalize for the galaxy's size. At low redshift there is some consensus that when expressed this way, isolated massive . Metallicity gradients of galaxies as a function of their size. On the horizontal axes we plot the disc scale-length of the galaxies. On the vertical axes we show the metallicity gradient. In panel (a) we have normalized the metallicity gradient to the disc scale length, while in (b) we display the same data but with the metallicity gradient expressed in physical units. Symbols indicate our metallicity gradient classification: triangles pointing up/down are galaxies with 2σ positive/negative gradients, large circles represent galaxies that have metallicity gradients consistent we being flat, and small circles indicate galaxies flagged for their sensitivity to the prior (the squares in Fig. 4(a)). Data points are coloured according to the galaxy's redshift. We overplot three red crosses, which indicate the median trend of the metallicity gradient with size. Errors on the median are determined by bootstrapping the sample, and Monte Carlo sampling of the errors (this includes unreliable galaxies). In panel (a), a grey shading denotes the region disallowed by our prior on the metallicity gradient.

galaxies ( > ∼ 10 8 M ) share a common value for the metallicity gradient (Sánchez et al. 2014;Ho et al. 2015). While this may not be true at higher redshift, we will nonetheless use scaled units.

From Fig. 5 it is clear that the average galaxy has a negative metallicity gradient. However, it is also clear that there is considerable scatter. There is a noticeable trend for the median metallicity gradient to become more negative with increasing galaxy size. We see that amongst the large galaxies ( ∼ > 3 kpc) there are no galaxies with positive Table 4. Analysis of unexplained (intrinsic) scatter in the metallicity gradients. We model the metallicity gradients as if they were normally distributed with mean, µ, and standard deviation, σ int . The galaxies are divided into two groups: those smaller, and those larger than r d = 3 kpc. We show results for metallicity gradients expressed in both physical units (dex/kpc) and scaled units (dex/r d ).

The unreliable galaxies, which are sensitive to the prior, are not included in this analysis. Including them would not significantly alter the results. metallicity gradients, regardless of redshift. Indeed almost all large galaxies have negative metallicity gradients. Conversely, the small galaxies ( ∼ < 3 kpc) present a range of metallicity gradients, some negative and some positive.

Gradient

The scatter in the metallicity gradient appears to sharply increase below ∼ < 3 kpc. However, we see that the error bars of the small galaxies also increase. It is therefore important to ask whether there is a true increase in the intrinsic scatter in the small galaxies. In Table 4 we present an analysis of the amount of intrinsic scatter in both the small and large galaxies. For this we model the metallicity gradients as normally distributed with mean, µ, and standard deviation, σ int . This analysis accounts for the full posterior distributions on the metallicity gradient. We see that there is indeed an increase in the intrinsic scatter in the small galaxies ( ∼ < 3 kpc) and that the mean is overall negative.

However, we caution that although we call this intrinsic scatter, it is perhaps more honest to call it unexplained scatter. It could simply be that our model is underestimating the true error in the metallicity gradient. Indeed it is somewhat suspicious that we see an increase in the scatter towards the smaller galaxies. We naïvely would expect small galaxies to be more sensitive to model specification errors, i.e. errors in the PSF or inclination (see Section 3.3). Nevertheless, the sharpness of the increase in scatter for galaxies with sizes ∼ < 3 kpc is certainly worthy of note.

It is well known than there exists a correlation between the size and total stellar mass of a galaxy (see van der Wel et al. 2014, and references therein). So, given that the galaxies with positive metallicity gradients are exclusively small, we also might expect them to be also less massive. However, this does not appear to be the case in Fig. 6, where we plot the metallicity gradients against the total stellar mass of the galaxy. Galaxies with positive metallicity gradients span a similar mass range as those with negative gradients. That said, it does appear that on average the metallicity gradient decreases with increasing mass. This trend would be qualitatively similar to that found in studies of low-redshift galaxies, e.g. Belfiore et al. (2017) and (somewhat indirectly) the works of Pérez et al. (2013) and Ibarra-Medel et al. (2016). However, because of our galaxy selection, it is hard to interpret the trend we observe. With a larger sample one might be able to disentangle the coupled selection effects on mass, size and redshift.

It is easier, however, to discuss the more massive portion of our sample. For example, at masses above ≈ 10 9 M we notice there is a large range of metallicity gradients, with hints that the scatter is perhaps greater in the higher-redshift galaxies. This certainly is in direct contrast to the low-redshift (z 0.1) literature result of a common metallicity gradient. This discrepancy might be resolved if many of these galaxies are not isolated, but are instead interacting with other galaxies. Because, interacting galaxies typically have flatter metallicity gradients than would otherwise be expected (Rich et al. 2012).

Since interactions could trigger increased star-formation rates in interacting galaxies, Stott et al. (2014) suggested that the scatter in the observed metallicity gradients could be explained by the SFR. In Fig. 7 we compare the metallicity gradient to the SFR, however, we do not observe a trend towards more positive metallicity at high SFRs. Admittedly the association between interacting galaxies and higher SFR might be weak. However, in fact we see a slight trend in the opposite direction with the metallicity gradient becoming more negative at high SFRs. We caution that this trend could be a manifestation of our sample selection (i.e. because our sample is incomplete). Notably the SFR of a galaxy in our sample correlates with its redshift. This is in part because at higher redshifts we are biased towards the brighter, more strongly star-forming galaxies. However, it is also in part because we expect the galaxy main-sequence to be evolving (the average SFR is less today than it was at z ≈ 1). So to put the galaxies on an even footing, we normalize the SFR relative to the main sequence. We define the mainsequence offset, ∆SFR, as the difference between a galaxy's SFR and the SFR of a galaxy on the main sequence which has identical mass and redshift, adopting the parametrisation of the main sequence evolution from Whitaker et al. (2012). In Fig. 8 we show the metallicity gradients of our galaxies against ∆SFR. We see that the trend towards negative gradients with increasing SFR has now disappeared. In fact, this projection shows no trend whatsoever; the average gradient is more or less constant, and there is large scatter irrespective of a galaxy's position relative to the main-sequence. This is in contrast to the results of Stott et al. ( 2014) who find a positive correlation between specific star-formation rate (sSFR) and metallicity gradient. They find galaxies that are vigorously forming stars have flatter or even positive gradients. The authors suggest that an event causing infall of metal-poor gas would simultaneously trigger intense star formation and reduce the central metallicity (thus flattening or inverting the metallicity profile). Our findings challenge this explanation for the cause of inverted metallicity gradients. However, our results do not automatically preclude the mechanisms of galaxy-galaxy interactions and/or cold flows for triggering star formation.

In the next section we shall discuss our results in more detail, placing them in context with other observations and theoretical work.

DISCUSSION

Literature comparison

The work presented here is the first large systematic study of metallicity gradients in galaxies between 0.1 ∼ < z ∼ < 0.8. However, there have been numerous studies of metallicity gradients in galaxies at lower and higher redshifts. In Fig. 9(a) we plot our results alongside several of these studies. Many of the high-redshift observations have used IFS techniques similar to ours, although some use observations of gravitationally lensed galaxies. Ideally these lensed observations will have sufficient resolution such that seeing does not significantly affect the observed metallicity gradient. Therefore it is worth noting that the lensed and non-lensed studies are consistent with one another, with both having means close to zero. Furthermore, apart from some lensed galaxies with very steep negative metallicity gradients, the scatter in the two distributions is similar.

There is, however, a slight discontinuity between our observations and the other studies at 0.8 ∼ < z ∼ < 2.6. This transition is perhaps most visible in the increased scatter of our observations. But there is also a shift in the average gradient to slightly negative gradients in our sample. While these effects could be attributed to the real evolution of metallicity gradients with cosmic time, a discontinuity would argue for less astrophysical causes.

There are a few plausible explanations, of which the most concerning for us is that this could highlight an issue with our method for determining metallicity gradients. A systematic overcompensation for seeing effects certainly could produce steeper gradients. It is harder, however, to conceive of systematic effects that would produce a shift away from a non-zero average gradient. In support of our results, a slightly negative average gradient is in fact consistent with observations in the low-redshift Universe (e.g. Sánchez et al. 2014). In particular, Pérez-Montero et al. (2016) report positive metallicity gradients occurring in 10% of their z ∼ 0.02 sample. This is similar to the (8 ± 3)% we find here (7 out of 84 galaxies). Note, given the 2σ significance level that we choose when classifying galaxies as having positive metallicity gradients, we estimate a non-negligible number of false positives (0.77 galaxies) 15 .

There is another potential reason for the discontinuity between the observational studies. The vast majority the 0.8 ∼ < z ∼ < 2.6 galaxies have metallicity gradients determined using the N2 ratio ([N ii]6584/Hα) and the calibration proposed by Pettini & Pagel (2004, hereafter PP04). Whereas, we use a very different method for deriving metallicity. It is well known that there are large discrepancies between different metallicity determination methods (see Kewley & Ellison 2008). So it is worth noting, although not demonstrated here, that combining the D13 photoionization models with our metallicity and ionization parameter coupling (equation 5) tends to provide similar metallicity gradients to those produced by the method of (Maiolino et al. 2008, hereafter M08). The calibration adopted by M08 has a shallower dependence of the N2 ratio on metallicity than that used by PP04. Therefore, using M08 yields systematically steeper 15 (26 + 7) × 2.275% -The number of galaxies with flat or positive metallicity gradients, multiplied by the one-tailed 2σ probability. metallicity gradients than PP04. And, consequently, when comparing our work with the other studies that use the PP04 calibration, one might also expect our metallicity gradients to be systematically steeper.

An advantage of our method is that it produces a selfconsistent metallicity gradient analysis, independent of the available emission-lines. With our MUSE observations we lose emission lines redward of Hα at z ∼ > 0.4. However, we do not observe a systematic shift in our measurements at z ≈ 0.4, suggesting that our method is indeed self-consistent.

On a related note, it is questionable whether any of these metallicity methods are valid at high redshift, since most metallicity calibrations and photoionization models are designed for low-redshift interstellar medium (ISM) conditions. For example, it is generally accepted that the electron density of the ISM was higher at earlier times (Shirazi et al. 2014;Sanders et al. 2016), but the D13 photoionization models we use are computed only at low densities (≈ 10 cm -3 ). At earlier times the ionization parameter, the hardness of the ionizing spectrum, or the nitrogen-to-oxygen abundance ratio may also have been different, although there is little consensus on this (e.g. Onodera et al. 2016;Kashino et al. 2017;Steidel et al. 2016, and references therein). Our method is certainly not immune to these issues. However, since our observations are all below z ≈ 0.8 it may not be inappropriate to apply the same assumptions as we use at low redshift. It is worth noting that because we marginalize over galaxy wide variations in the ionization parameter, our method may partially mitigate against some of the variations in ISM conditions. This is because (to first order) the largest variations in nebular emission-line spectra are typically due to metallicity and ionization parameter [START_REF] Dopita | Astrophysics of the diffuse universe Dopita[END_REF].

Interpretation

In the work presented here we find that, on average, galaxies between 0.1 ∼ < z ∼ < 0.8 have a negative metallicity gradient. However, there is considerable scatter, with some galaxies exhibiting a positive metallicity gradient. As reported by Pilkington et al. (2012), numerical simulations predict that there is intrinsic scatter in the metallicity gradient. This can be seen in Fig. 9(b) when comparing the scatter within the MUGS galaxies (Stinson et al. 2010) and within the RaDES galaxies (Few et al. 2012). One should note that the difference between the MUGS and RaDES simulations are, however, attributable to differences in the numerical recipes for star formation and feedback in the various simulations. Flatter metallicity gradients arise from the model prescriptions with more intense feedback. For example, the MaGICC simulations (Gibson et al. 2013) re-simulate two of the MUGS galaxies using an enhanced feedback recipe. The increased feedback produces galaxies with metallicity gradients that are consistently flatter at all redshifts. Gibson et al. (2013) attribute this to central gas that is lost in outflows and later recycled back into the galaxy, but at larger radii.

The intrinsic diversity in metallicity gradients that we observe (and which is also present in the simulations) would indicate that at earlier times there is no common metallicity gradient. It is perhaps surprising that we observe a common abundance gradient in the Universe today. That said, as we shall later suggest, this contradiction is not necessarily as Troncoso et al. (2014) and Wuyts et al. (2016). We also overplot a red trend line which indicates the median metallicity gradient in equally spaced redshift bins (excluding both the unreliable and z ≈ 0 galaxies). The shaded region around this line represents the 1σ scatter of the median, derived from a bootstrap and Monte Carlo resampling of the errors. (Right) Observational results (ours and others) compared with various numerical simulations. The grey shaded region indicates the ±1σ spread of metallicity gradients from 19 RaDES simulated galaxies (Few et al. 2012). Thin coloured lines indicate four individual MUGS simulated galaxies (Stinson et al. 2010). Two of these galaxies (the orange and green lines) were re-simulated within the MaGGIC simulations (Gibson et al. 2013) and are shown as thick orange/green lines.

acute as it would first appear. A common metallicity gradient may only exist in large galaxies.

Even though these simulations can produce a large range of metallicity gradients, they typically do not reproduce the same inverted metallicity gradients that we and others observe. The simulations tend to produce galaxies with negative metallicity gradients (both steep and shallow). Some of these simulations have next to no redshift dependence, while others suggest that the metallicity gradient was steeper at earlier times. Despite this variety, the simulations never produce galaxies with positive metallicity gradients. Indeed, with the classical understanding of inside-out growth one expects negative metallicity gradients (Portinari & Chiosi 1999). Radial mixing of gas could flatten the metallicity gradient, but it is hard to conceive of secular processes that could produce positive metallicity gradients. Consistent with other higher redshift studies, we identify some of our galaxies that have significantly positive metallicity gradients. Contrast this with galaxies at much higher redshift (z ≈ 3.4), which all have centres that are systematically more metal poor than their outskirts (Troncoso et al. 2014). In related work, Cresci et al. (2010) attribute this metallicity gradient inversion to cold flows (e.g. Dekel et al. 2009). Cold flows are streams of cold gas that penetrate the hot galaxy halo and fuel star formation in galaxies. It is argued that if this metal-poor cold gas can reach the innermost regions of a galaxy (the deepest part of the gravitational potential) then this would explain the metallicity gradient inversion.

It should be noted that cold flows are not the only way to transport metal-poor gas to the inner portions of a galaxy. As suggested by Rich et al. (2012), galaxy-galaxy interactions could cause radial inflow of gas though the galaxy's disc. Metal-poor gas from the outskirts is therefore deposited in the galaxy centre. This is supported by simulations, which find that interacting galaxies have flatter (although never positive) metallicity gradients (Rupke et al. 2010;Torrey et al. 2012;Sillero et al. 2017). However, work by Divoy (2014, Ch. 8) suggested that merging systems can perhaps exhibit positive metallicity gradients once simulations have been downgraded to mimic the effects of spatial resolution loss and/or noise. Nevertheless, it is not clear whether galaxy interactions can truly invert the metallicity gradient.

In contrast to the inflow mechanisms above, outflows could also provide an alternative explanation for positive metallicity gradients. Intense centrally-concentrated star formation could produce significant outflows that entrain metalrich gas. Using simple analytical chemical evolution arguments, Cresci et al. (2010) disfavoured a wind scenario as it would require mass outflow rates far in excess of the observed SFR. That said, Troncoso et al. (2014) point out that if this metal-rich gas falls back preferentially onto the outer regions of the galaxy, it could raise the outer metallicity. This "fountaining" could then enhance the ability of winds to produce positive metallicity gradients and hence reduce the required outflow rates down to more realistic levels.

The inflow and outflow scenarios need not be mutually exclusive. Gas accretion could trigger intense star formation that drives the winds. In either case (inflows and/or outflows) we should expect to observe elevated star-formation rates. In support of this Stott et al. (2014) and Wuyts et al. (2016) find weak correlations between the specific star-formation rate (sSFR) and metallicity gradient (galaxies with more intense star-formation have flatter or even positive metallicity gradients). In contrast, we do not find systematically different metallicity gradients in galaxies with elevated SFR (Fig. 8). A difference in the sample selection might plausibly account for the difference between our results and previous studies. Alternatively, one should also consider that if the intense star-formation is confined to a galaxy's centre, the global SFR of a galaxy may not be a very sensitive indicator of inflow or outflow events. Furthermore, the timescale for observing an inverted metallicity gradient might be different from that to observe an elevated SFR. As such, some of the galaxies observed with elevated SFRs are likely moving vertically down the main sequence (i.e. have a declining SFR at the instant we observe them). Therefore timescale differences might further weaken any correlation between the specific star-formation rate and metallicity gradient.

There is perhaps another way to reconcile the lack of correlation between metallicity gradient and star-formation intensity. Contrary to our earlier assertion, it has been recently proposed by Schönrich & McMillan (2017) that inverted metallicity gradients can actually arise mostly within the framework of secular inside-out growth. This is achieved in conjunction with outflows that recycle enriched gas from the galaxy centre, transporting it to larger radii. Schönrich & McMillan (2017), however, show that these positive metallicity gradients do not persist throughout a galaxy's life, and in fact are only expected to exist at early times in a galaxy's evolution. Unfortunately, given the large number of model uncertainties, they are unable to make rigorous predictions for how long this phase may last.

We find the results of Schönrich & McMillan (2017) intriguing as they may shed light on another curious result of our study. In the small galaxies (r d < 3 kpc) we observe a large scatter of the metallicity gradient, whilst larger galaxies present negative metallicity gradients with minimal scatter (Fig. 5). If the smaller galaxies are comparatively less evolved systems, then this secular evolution (with metallicity gradient inversion at early times) may account for the diversity of metallicity gradients observed. Consequently the larger galaxies, with their negative metallicity gradients, would be emblematic of the more classical understanding of inside-out growth which produces a common metallicity gradient (e.g. Prantzos & Boissier 2000). This would reconcile our results with the low-redshift Universe. And as such we could expect to find inverted metallicity gradients in small galaxies today.

The idea of that negative metallicity gradients only exist in galaxies with well-established discs is broadly similar to the findings of Ma et al. (2017). With their numerically simulated galaxies they find that galaxies with the steepest negative metallicity gradients all have well-ordered, rotating gas discs. (N.B. by itself, this is a necessary, but insufficient condition.) And conversely, the galaxies with perturbed discs all have flat metallicity gradients.

Nonetheless we strongly caution that their results are only conceptually similar to ours, not least because their work does not support a dependence between a galaxy's metallicity gradient and its size, which our work would suggest.

CONCLUSIONS

Using MUSE data we present metallicity gradients for a sample of 84 intermediate-redshift galaxies (0.1 ∼ < z ∼ < 0.8). We infer the true metallicity gradient using a forward-modelling technique that corrects for the seeing effects. We search for trends of the observed metallicity gradient with global properties such as galaxy mass, size and SFR. From this we conclude the following:

• The average galaxy in our sample has a negative metallicity gradient. Nevertheless, there is significant scatter and we classify 7 of the 84 galaxies as having positive metallicity gradients.

• We do not identify any significant correlation of the metallicity gradient with either total SFR or stellar mass (Figs. 6 &7). And we find no correlation whatsoever once the SFR is normalized relative to the main sequence (Fig. 8).

• This lack of correlation in the latter runs contrary to previous studies. These studies, predominantly at higher redshifts (0.6 ∼ < z ∼ < 2.6), cite such a correlation as being suggestive of sudden gas inflow or merger events that might trigger star formation while simultaneously flattening/inverting the metallicity gradient. Our results, however, do not support this interpretation at the intermediate redshifts we study.

• The largest galaxies in our sample (r d > 3 kpc) are found to have almost exclusively negative metallicity gradients. On the contrary, the smaller galaxies (r d < 3 kpc) present a range of metallicity gradients (Fig. 5). The small galaxies exhibit a larger intrinsic scatter in the metallicity gradient.

• We liken the large galaxies as being similar to galaxies observed in the Universe today, where galaxies present a common metallicity gradient. In addition, we speculate that the size dichotomy may be related to a secular understanding of inside-out growth, where a common metallicity gradient is only established in large (well-evolved) systems.

• However, we advise caution on this last point. Intrinsically smaller galaxies are also more affected by seeing-induced resolution loss. Therefore, for these small galaxies our inferred metallicity gradients (and their errors) are presumably more dependent on the accuracy of our modelling (see Section 3.3).

The MUSE GTO surveys are currently ongoing and in the future will provide a larger sample of metallicity gradients. This data will allow us to separate potential selection effects and biases.

In the future we will also study metallicity gradients in conjunction with the gas kinematics of the galaxies. By also quantifying the local environments of our galaxies, we can study the impact of galaxy-galaxy interactions on the metallicity gradient, without relying on SFR as an indirect tracer. . Plot showing derived model parameters for HDFS 0003. This is an illustrative example and similar plots for all galaxies can be found as online supplementary material. We show both 1D and 2D marginalized histograms for all 5 parameters: the total star-formation rate, log 10 (SFR), central metallicity, log 10 Z 0 , metallicity gradient, ∇ r (log 10 Z), ionization parameter at solar metallicity, log 10 U , and V-band optical depth, τ V . In each 1D histogram the vertical lines indicate the median (solid), ±1σ quantiles (dashed) and ±2σ quantiles (dash-dotted). In each 2D histogram we plot 1σ and 2σ contours. All axes span a [-4σ , 4σ ] interval in their respective parameters.

Figure 1 .

 1 Figure 1. Diagnostic plots for AGN classification. Our final sample of star forming galaxies and included intermediate galaxies are plotted as circles and squares, respectively. Excluded intermediate and AGN/LINER galaxies are indicated by triangles and stars, respectively. Data points are colour coded by redshift. (Left) The [O iii]/Hβ versus [N ii]/Hα diagnostic (a BPT (Baldwin et al. 1981) diagram) for the low-redshift (z ∼ < 0.4) portion of our sample. As a solid line we indicate the Kewley et al. (2001) theoretical maximum limit for star formation alone. With a dashed line we show the Kauffmann et al. (2003) curve, the empirical division between emission from star formation and AGN/LINER emission. (Right) The Mass-Excitation (MEx) diagram for whole our sample, including also the Juneau et al. (2014) demarcations as solid black lines. Note that the same galaxies that appear in the left panel (the BPT diagram) also appear in the right panel (the MEx diagram).

Figure 2 .

 2 Figure2. The distributions of various global galaxy properties of our final sample are shown as coloured histograms. The grey histograms show the distribution of galaxies from the parent sample. Galaxies in the final sample are coloured by their redshift. If the histogram bins were independent of redshift then each bin would be similarly coloured. We show four properties: (a) stellar mass, (b) star-formation rate, (c) rest-frame Bessel B-V colour and (d) offset from the star-forming main sequence, accounting for redshift evolution. Rest-frame colours are determined from the best-fit magphys spectrum to the photometry (see Section 2.2.1). While most galaxies in the parent sample have reliable photometry, allowing us to derive masses and colours, fewer galaxies in the parent sample have detected line-emission. Consequently, there are fewer galaxies in the parent sample in panels (b & d).

Figure 3 .

 3 Figure3. Mass versus SFR for our final sample plotted as coloured circles. For comparison, we display the main sequence at four different redshifts as solid lines, adopting the parametrization ofWhitaker et al. (2012). Note that the SFR of the low redshift galaxies are derived from Hα and Hβ lines, whilst those for the high redshift galaxies are derived from Hβ and faint Hγ lines. Consequently the SFR errors are much smaller for the low redshift galaxies (z ∼ < 0.4).

Figure 4 .

 4 Figure 4. Sensitivity analysis of model inputs. On the horizontal axes we plot the metallicity gradients derived from our fiducial model. On the vertical axes we compare the metallicity gradients derived from the same data, but with the alternate model. The alternative models are identical to the fiducial model, except for one fixed parameter that has been adjusted slightly (see text for details). A different alternative model is plotted in each of the four panels, showing the effects of: (a) the width of the metallicity gradient prior, (b) a misestimation of the PSF, (c) a miscalculation of the inclination and (d) misidentification of the galaxy centre. The error bars show the ±1σ quantiles of the metallicity gradient posterior probability distribution. N.B. the errors of the two models are not independent and are likely to be highly correlated, consequently the errors appear much larger than the scatter. Each data point is coloured by the Jensen-Shannon distance between the 1D marginal posteriors of the fiducial and alternate models, highlighting the datapoints that deviate most from the 1:1 relation (shown as a red dotted line). Squares indicate a subset of galaxies which have large Jensen-Shannon distances (JS dist > 0.3) when comparing the fiducial [-0.3, 0.3] and the [-0.5, 0.5] metallicity gradient priors. We deem these galaxies as unreliable (see Section 3.3.1). In each panel a solid black line shows the best fit to the data, excluding unreliable galaxies and weighted by the inverse of the formal data errors. A red shaded region indicates the (weighted) vertical scatter about this 1:1 relation, again excluding the unreliable galaxies. Equations for the best fit and scatter are given in the bottom-right corner of each panel.

  Figure5. Metallicity gradients of galaxies as a function of their size. On the horizontal axes we plot the disc scale-length of the galaxies. On the vertical axes we show the metallicity gradient. In panel (a) we have normalized the metallicity gradient to the disc scale length, while in (b) we display the same data but with the metallicity gradient expressed in physical units. Symbols indicate our metallicity gradient classification: triangles pointing up/down are galaxies with 2σ positive/negative gradients, large circles represent galaxies that have metallicity gradients consistent we being flat, and small circles indicate galaxies flagged for their sensitivity to the prior (the squares in Fig.4(a)). Data points are coloured according to the galaxy's redshift. We overplot three red crosses, which indicate the median trend of the metallicity gradient with size. Errors on the median are determined by bootstrapping the sample, and Monte Carlo sampling of the errors (this includes unreliable galaxies). In panel (a), a grey shading denotes the region disallowed by our prior on the metallicity gradient.

Figure 6 .Figure 7 .

 67 Figure 6. Metallicity gradients of galaxies as a function of their stellar mass. See Fig. 5 for a description of the symbols and colours.

Figure 8 .

 8 Figure 8. Metallicity gradients of galaxies as a function of ∆SFR, the SFR of a galaxy relative to what is expected for a galaxy on the main-sequence with identical mass and redshift. We use the main-sequence parametrization of Whitaker et al. (2012). See Fig. 5 for a description of the symbols and colours.

Figure 9 .

 9 Figure9. Comparison of our results with the literature, including both observations and numerical simulations. (Left) We compare our observational results against other studies. These metallicity gradients come from studies of galaxies that have been gravitationally lensed (star symbols), and non-lensed galaxies (squares). Where possible we colour data points to indicate the SFR relative to the main sequence (∆SFR). The z ≈ 0 data point represents the mean and scatter in metallicity gradients reported inSánchez et al. (2014). We compile lensing results from Yuan et al. (2011), Jones et al. (2013), Jones et al. (2015), Leethochawalit et al. (2016), Wang et al. (2017) and Patrício et al. (in prep.). The non-lensed results are gathered from Queyrel et al. (2012), Swinbank et al. (2012), Stott et al. (2014), Troncoso et al. (2014) andWuyts et al. (2016). We also overplot a red trend line which indicates the median metallicity gradient in equally spaced redshift bins (excluding both the unreliable and z ≈ 0 galaxies). The shaded region around this line represents the 1σ scatter of the median, derived from a bootstrap and Monte Carlo resampling of the errors. (Right) Observational results (ours and others) compared with various numerical simulations. The grey shaded region indicates the ±1σ spread of metallicity gradients from 19 RaDES simulated galaxies(Few et al. 2012). Thin coloured lines indicate four individual MUGS simulated galaxies(Stinson et al. 2010). Two of these galaxies (the orange and green lines) were re-simulated within the MaGGIC simulations(Gibson et al. 2013) and are shown as thick orange/green lines.

Figure A1 .Figure A2 .

 A1A2 Figure A1. Image atlas showing HST cutouts of galaxies in our sample, shown in order of increasing metallicity gradient. For the UDF galaxies with display the F775W filter images, for all other galaxies we use the F814W filter. Each image is centred on the galaxy of interest, where North is up and East is left. The image size is 16 kpc (physical) edge-to-edge. Two red ellipses indicate the derived galaxy morphology, drawn at the 50% (1.68 r d ) and 90% (3.89 r d ) light radii. (Note that for the largest galaxies the 90% radius lies outside of the cutout region shown.) In the top-left corner of each panel we show the galaxy's redshift. In the bottom-left corner we display the inferred metallicity gradient (units are dex/kpc). The border colour displays the metallicity gradient type: negative gradients (blue), positive gradients (orange), flat gradients (green). Galaxies flagged as unreliable have dashed borders.

Table 2 .

 2 Priors on model parameters. For each parameter we detail the type of prior and the range of values covered.

	Parameter	Prior type	Range
	SFR tot	Logarithmic	[0.01, 100] M yr -1
	log 10 Z 0	Linear	≈ [-1.30, 0.70] dex
	∇ r (log 10 Z)	Linear	[-0.3, 0.3] dex/kpc
	log 10 U	Linear	≈ [-5.02, -1.42] dex
	τ V	Linear	[0, 4]

MNRAS 000, 1-27 (2018)

A short description of the pipeline can be found inWeilbacher et al. (2012).

Slice refers to the optical image slicers within MUSE, not the wavelength layers (channels).

Volumetric pixels.

The fit is performed in the Fourier space.MNRAS 000, 1-27 (2018)

The ratio between the radii that contain 50% and 90% of a galaxies light.

We adopt the sign convention that expresses emission as negative equivalent widths.

We define ∆SFR to be the difference in the observed SFR relative to what would be expected for a galaxy on the main sequence, with the same stellar mass and at the same redshift. Here and throughout this paper we adopt the main sequence parametrization ofWhitaker et al. (2012).

Spatial pixels.MNRAS 000, 1-27(2018)

http://spectro.princeton.edu/idlspec2d_install.htmlMNRAS 000, 1-27 (2018)

This limitation is by no means exclusive to our model. To properly treat poorly resolved galaxies, one would need a forward model with some non-parametric description of the 2D metallicity profile.

We fit [O ii]3726, [O ii]3729, [O ii]3726,3729, [Ne iii]3869, Hε, Hδ , Hγ, Hβ , [O iii]4959, [O iii]5007, [N ii]6548, Hα, [N ii]6584, [S ii]6717, [S ii]6731 and [S ii]6717,6731. We exclude redundant lines, i.e. if [O ii]3726,3729 coadd is detected at S/N ≥ 5, then we do not also fit [O ii]3726 and [O ii]3729, even if they individually have a S/N ≥ 5 MNRAS 000, 1-27 (2018)

This estimate is the mean uncertainty obtained after averaging the two alternative models (one where the PSF is 10% smaller and the other where the PSF is 10% larger). Therefore this number differs slightly from the one presented in Fig.4(b)MNRAS 000, 1-27(2018)

ACKNOWLEDGEMENTS

Thank you to Josephine Kerutt and Rikke Saust for their helpful discussions and providing PSF measurements for the CDFS fields.

This work was performed using the python and idl programming languages. In particular we heavily used the following packages: NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001-), matplotlib (Hunter 2007) and astropy (Astropy Collaboration et al. 2013). 20 David Carton et al. DC, JR, and VP acknowledge support from the ERC starting grant 336736-CALENDS. JB acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through national funds (UID/FIS/04434/2013) and by FEDER through COM-PETE2020 (POCI-01-0145-FEDER-007672). During part of this work, JB was supported by FCT through Investigador FCT contract IF/01654/2014/CP1215/CT0003. TC acknowledges support of the ANR FOGHAR (ANR-13-BS05-0010-02), the OCEVU Labex (ANR-11-LABX-0060), and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the "Investissements d'avenir" French government program. JS acknowledges support from the Netherlands Organisation for Scientific Research (NWO), through VICI grant 639.043.409, and the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement 278594-GasAroundGalaxies.

APPENDIX A: ADDITIONAL FIGURES AND TABLES

In Fig. A1 we provide an HST image atlas of our final galaxy sample. With each cutout image we indicate the inferred metallicity gradient. In Figs. A2 & A3 we provide example plots for inspecting the quality of the model fit. Similar plots for each galaxy in our final sample are provided in the online supplementary material. In Table A1 provide the derived metallicity gradients and other data pertaining to our galaxy sample.