Resolved scaling relations and metallicity gradients on sub-kiloparsec scales at z=1

TL;DR

This study uses gravitational lensing and advanced modeling to analyze resolved star formation and metallicity relations in galaxies at z=1, revealing that these relations are already established at sub-kiloparsec scales and are consistent with local galaxy patterns.

Contribution

First detection of the spatially resolved Mass-Metallicity Relation at z=1 on sub-kpc scales using lensing and 2D gradient modeling.

Findings

Confirmed the presence of the rSFMS at z=1 on sub-kpc scales.

Established the existence of the rMZR at z=1.

Measured metallicity gradients consistent with linear profiles.

Abstract

The existence of a spatially resolved Star-Forming Main Sequence (rSFMS) and a spatially resolved Mass-Metallicity Relation (rMZR) is now well established for local galaxies. Moreover, gradients with metallicity decreasing with radius seem to be common in local disc galaxies. These observations suggest that galaxy formation is a self-regulating process, and provide constraints for galaxy evolution models. Studying the evolution of these relations at higher redshifts is still however very challenging. In this paper, we analyse three gravitationally lensed galaxies at z = 0.6, 0.7 and 1, observed with MUSE and SINFONI. These galaxies are highly magnified by galaxy clusters, which allow us to observe resolved scaling relations and metallicity gradients on physical scales of a couple of hundred parsecs, comparable to studies of local galaxies. We confirm that the rSFMS is already in place at these redshifts on sub-kpc scales, and establish, for the first time, the existence of the rMZR at higher redshifts. We develop a forward-modelling approach to fit 2D metallicity gradients of multiply imaged lensed galaxies in the image plane, and derive gradients of -0.027+/-0.003, -0.019+/-0.003 and -0.039+/-0.060 dex/kpc. Despite the fact that these are clumpy galaxies, typical of high redshift discs, the metallicity variations in the galaxies are well described by global linear gradients, and we do not see any difference in metallicity associated with the star-forming clumps.

Tables

Table 1. Table 1: Sample properties derived by Patrício et al. ( 2018 ) using MUSE and HST data. From left to right: instrument, observation program identification, point spread function FWHM measured using a Moffat profile, redshift, stellar mass, magnification-corrected star-formation rate from dust-corrected Balmer lines and effective radius, calculated from the disc length (R d ) measured in Patrício et al. ( 2018 ) (table 2) using the F160W HST source plane images as R e =1.67835 R d .

Object	$\alpha$	$\delta$	Inst.	Program ID	PSF	$z$	${\log }_{10}$ M⋆	SFRMUSE	Re
	J2000	J2000			["]		[M⊙]	[M⊙/yr]	[kpc]
AS1063-arc	22:48:42	-44:31:57	MUSE	060.A-9345(a)	1.03	0.6115	10.94$\pm$0.05	41.5$\pm$4.0	7.7$\pm$0.2
A370-sys1	02:39:53	-01:35:05	MUSE	094.A-0115, 096.A-0710	0.70	0.7251	10.40$\pm$0.02	3.1$\pm$0.3	12.0$\pm$0.7
M1206-sys1	12:06:11	-08:48:05	SINFONI	087.A-0700	0.78	1.0366	10.90$\pm$0.06	107.3$\pm$30.7	11.1$\pm$0.2

Table 2. Table 2: AS1063-arc and A370-sys1 metallicity gradient and morphology fit. GALFIT: results of the morphological fit to the reconstructed F160W HST band using galfit . Remaining rows: fit of the image plane metallicity gradient using the procedure described in 3.3 , fixing or letting the morphological parameters vary.

AS1063-arc
		$\nabla$Z	Z0	Centre RA	Centre Dec	q	$\theta$	${\chi }^2$/dof
		[dex/kpc]	[12+$\log$(O/H)]	J2000	J2000		[deg]
galfit	-	-	-	22h48m42.859s	-44d31m57.0464s	0.56	-32	37.25
Fixed morph.	prior	[-0.1:0.0]	[8.5:9.5]	22h48m42.859s	-44d31m57.0464s	0.56	-32
	fit	-0.034$\pm$0.002	8.985$\pm$0.007	-	-	-	-	6.04
Free par.	prior	[-0.1:0.0]	[8.5:9.5]	22h48m[41.634 : 41.871]s	-44d31m[55.169 : 57.708]s	[0.1:0.9]	[-90:90]
	fit	-0.042$\pm$0.002	9.038$\pm$0.008	22h48m41.750s	-44d31m56.016s	0.52$\pm$0.05	68$\pm$2	1.44
A370-sys1
		$\nabla$Z	Z0	Centre RA	Centre Dec	q	$\theta$	${\chi }^2$/dof
		[dex/kpc]	[12+$\log$(O/H)]	J2000	J2000		[deg]
galfit	-	-	-	02h39m53.716s	-01d35m03.55s	0.32	-52	42.04
Fixed morph.	prior	[-0.1:0.0]	[8.5:9.5]	02h39m53.716s	-01d35m03.55s	0.32	-52
	fit	-0.039$\pm$0.004	8.980$\pm$0.007	-	-	-	-	4.79
Free par.	prior	[-0.1:0.0]	[8.5:9.5]	02h39m[53.573 : 53.805]s	-01d35m[02.921 : 07.817]s	[0.1:0.9]	[-90:90]
	fit	-0.053$\pm$0.004	9.032$\pm$0.009	02h39m53.709s	-01d35m04.169s	0.39$\pm$0.04	-47$\pm$3	3.80

Table 3. Table 3: Comparison between metallicities derived in P18, using the full set of lines available in MUSE and the Charlot & Fall ( 2000 ) extinction law, and the metallicities derived using only the strongest lines and the Calzetti et al. ( 2000 ) extinction law. Z: the metallicity, in 12+ log \log (O/H) ; E(B-V): dust attenuation in magnitudes; τ v subscript 𝜏 𝑣 \tau_{v} : dust attenuation (adimensional).

Object	Line Ratios	Calzetti et al. (2000)		Charlot & Fall (2000)
		Z	E(B-V)	Z	${\tau }_v$
AS1063-arc	P18	-	-	8.82$\pm$0.02	1.09$\pm$0.12
AS1063-arc	O2, O3, O32, R23, H$\beta$/H$\gamma$	8.76$\pm$0.10	0.46$\pm$0.09	8.75$\pm$0.10	1.11$\pm$0.20
A370-sys1	P18	-	-	8.88$\pm$0.02	0.44$\pm$0.11
A370-sys1	O2, O3, O32, R23, H$\beta$/H$\gamma$	8.81$\pm$0.17	0.38$\pm$0.19	8.80$\pm$0.16	0.88$\pm$0.48
M1206-sys1	P18	-	-	8.91$\pm$0.06	0.74$\pm$0.33
M1206-sys1	O2, Ne3O2, N2, H$\gamma$/H7, H$\gamma$/H$\delta$	8.89$\pm$0.05	0.92$\pm$0.11	8.91$\pm$0.05	1.86$\pm$0.21
M1206-sys1	O2, Ne3O2, N2, H$\gamma$/H7, H$\gamma$/H$\delta$ , H$\alpha$/H$\delta$ , H$\alpha$/H$\gamma$	8.87$\pm$0.07	1.00$\pm$+0.00	8.88$\pm$0.08	2.00$\pm$0.00
M1206-sys1	N2	8.94$\pm$0.07	-	8.94$\pm$0.07	-

Equations

$$lnp=-\frac{1}{2}\sum _r\left[{\left(\frac{{\mathrm{M}}_{\mathrm{r}}(Z)-{\mathrm{O}}_{\mathrm{r}}(E(B-V))}{{\sigma }_r^2}\right)}^2+\ln (2\pi {\sigma }_r^2)\right]$$

Full text

††pubyear: 2019††pagerange: Resolved scaling relations and metallicity gradients on sub-kiloparsec scales at $z\approx 1$–C

1 Introduction

It has now been well established that the masses, star-formation rates, and gas metallicities of star-forming galaxies are tightly correlated by two relations: the Star-Forming Main Sequence (SFMS), that relates stellar mass and star-formation rates, and the Mass-Metallicity Relation (MZR), relating mass and metallicity. These scaling relations have been observed from $z=0$ up to $z=6$ (e.g. Brinchmann et al., 2004; Tremonti et al., 2004; Erb et al., 2006; Whitaker et al., 2012; Speagle et al., 2014). It has even been argued that these three properties are connected by a single plane, the fundamental mass-metallicity relation (Lara-López et al., 2010; Mannucci et al., 2010), that does not evolve with redshift, although its existence is still controversial (e.g. Sánchez et al., 2013; Erroz-Ferrer et al., 2019).

These scaling relations are well explained by "reservoir" models. In these analytical models, after an initial phase of gas accretion, galaxies self-regulate their star-formation rates, evolving in a quasi-steady state (e.g. Schaye et al., 2010; Dutton et al., 2010; Bouché et al., 2010; Davé et al., 2012; Lilly et al., 2013). These models can successfully predict the SFMS and MZR using only a couple of fairly simple "regulators" of star-formation, such as gas infall rates, outflow rates, and gas recycling rates, without involving any details about star-forming processes. Since it is possible to reproduce these scaling relations without specifying any details on star-forming or stellar feedback processes, additional observables are needed to further our understanding of galaxy evolution.

In recent years, with the increasing number of Integral Field Unit (IFU) spectrograph surveys of local disc galaxies (e.g. CALIFA Sánchez et al. 2014, MaNGA Bundy et al. 2015, SAMI Croom et al. 2012 and MAD Erroz-Ferrer et al. 2019), it has also been established that both the Star-Forming Main Sequence and the Mass-Metallicity Relation exist on sub-galactic scales (e.g. Sánchez et al., 2013; Erroz-Ferrer et al., 2019). Moreover, it has been argued that these resolved relations are in fact more fundamental than the integrated ones (Rosales-Ortega et al., 2012; Barrera-Ballesteros et al., 2016), i.e., that the galaxy wide scaling-relations are a consequence of the local relations between stellar mass surface density, star-formation surface density and metallicity.

It is unclear if the reservoir models can be extended to explain these resolved relations, since they are based on isolated galactic systems, rather then contiguous and possibly interacting kpc-scale regions (but see, for example, Ho et al. 2015; Carton et al. 2015; D’Eugenio et al. 2018). It is also not clear what the reservoir would correspond to in this case and how the equilibrium phase would be reached. New and additional observables are needed to advance these simple but powerful models of galaxy evolution, as well as to test complex simulations that include sub-grid recipes for smaller-scale physical processes (e.g. Trayford & Schaye, 2018).

Another area of rapid development thanks to recent IFU surveys is the study of metallicity gradients. In the local Universe, disc galaxies are commonly observed to have higher metallicities in the centre than in the outskirts (a negative metallicity gradient) (e.g. Pilyugin et al., 2015; Ho et al., 2015; Carton et al., 2015; Sánchez-Menguiano et al., 2016; Belfiore et al., 2017), possibly with a universal slope when normalised to the galaxy size (e.g. Sánchez et al., 2014; Ho et al., 2015). The negative metallicity gradients can be explained by the ’inside-out’ disc growth scenario, where the inner parts of galaxies are formed at earlier times and are, consequently, more metal enriched than the outskirts (Larson, 1976). However, models that predict metallicity gradients compatible with the ones observed locally, make different predictions for gradients at earlier epochs, predicting either a steepening of the gradient at earlier epochs (e.g. Pilkington et al., 2012), or a flattening (e.g. Mott et al., 2013).

Deriving metallicity gradients at high-$z$ remains challenging. While in the local Universe metallicity gradients are generally negative, at high-$z$ a wide range of gradients, from negative to positive, has been measured. Wuyts et al. (2016) measured the metallicity gradients of star-forming galaxies at $z=0.6-2.7$, finding that they are, on average, flat. At slightly lower redshifts, $z=0.1-0.8$, Carton et al. (2018) find a negative median gradient, but with a large scatter (8% of their sample have significant positive gradients and 31% are consistent with flat gradients).

The evolution of the resolved scaling relations with cosmic time also remains difficult to probe, since it requires both a high signal-to-noise ratio and a high spatial resolution. Wuyts et al. (2013), and more recently Abdurro’uf & Akiyama (2018), have measured the resolved Star-Forming Main Sequence on kilo-parsec scales in massive galaxies (M⋆ > 1010 M⊙) at $z=0.7-1.8$ using multi-band high-resolution HST imaging, finding that the rSFMS was already in place at those redshifts. On the other hand, the resolved Mass-Metallicity Relation has still not been studied outside the local Universe until now.

In this work, we combine IFU optical and IR data from MUSE (Bacon et al., 2010) and SINFONI (Eisenhauer et al., 2003) observations of strongly gravitationally lensed arcs at z$\approx$1 to derive metallicity using multiple line-ratio diagnostics, and dust-corrected SFR from emission lines at physical scales of only a couple of hundreds parsecs. Using these data, we probe the metallicity gradients, resolved Star-Forming Main Sequence and the resolved Mass-Metallicity Relation of typical z$\approx$1 star-forming disc galaxies.

We analyse a sample of 3 strongly lensed galaxies in the background of the Abell S1063/RXJ2248-4431 (AS1063), Abell 370 (A370) and MACSJ1206.2-0847 (M1206) lensing clusters. These gravitational arcs were selected for their large size in the image plane (i.e. as seen in the sky). Despite their high magnification, these galaxies are quite typical of $z=1$ rotating discs. We have presented their basic properties derived from MUSE and HST data in a previous paper, Patrício et al. (2018). Here, we combine MUSE and SINFONI data to derive the resolved metallicity maps for three of those objects.

This paper is organised as follows. In Section 2 we present the MUSE and SINFONI data used in this work. In Section 3 we describe our method to derive metallicity and extinction from line fluxes. In Section 4 we analyse the local scaling relations and in Section 5 we derive the resolved metallicity maps and describe how we account for lensing. We discuss and summarise our results in Section 6.

Throughout this paper, we adopt a $\Lambda$-CDM cosmology with $\Omega$=0.7, $\Omega_{m}$=0.3 and H0 = 70 km s*-1* Mpc*-1*. We adopt a solar metallicity of $12+\log$(O/H) = 8.69 (Allende Prieto et al., 2001) and the Chabrier (2003) stellar initial mass function.

2 Sample and Data Reduction

The physical properties of the three galaxies analysed in this work were derived in a homogeneous way from HST and MUSE data in Patrício et al. (2018) (see Table 1 for a summary). They have redshifts between 0.6 and 1.0, stellar masses around 1010 M⊙, and are compatible with the Fundamental Mass-Metallicity relation (Lara-López et al., 2010; Mannucci et al., 2010) up to 0.1 dex. The stellar masses were derived fitting multiple HST bands and the MUSE integrated spectra using prospector111prospector (https://doi.org/10.5281/zenodo.1116491), a SED fitting code, Conroy et al. (2009) stellar models and the Chabrier (2003) initial mass function. Emission lines were masked during this fit. Dust corrected star-formation rates were calculated from emission lines from the MUSE data, making use of H$\gamma$ or $H\beta$. From the kinematic analysis of the [O ii] $\lambda$3727,29  emission, we concluded that these are typical rotating discs, representative of the population of star-forming galaxies at $z\approx 1$.

For the two lowest redshift galaxies analysed here, AS1063-arc and A370-sys1 (lensed galaxies in the clusters AS1063 and A370), we use MUSE data to derive metallicity maps from optical emission lines ([O ii] $\lambda$3727,29, [O ii] $\lambda$3727, H$\gamma$ , H$\delta$  and [O iii] $\lambda$5007). M1206-sys1 was also observed with MUSE, and its global metallicity can be derived from the integrated spectrum using [Ne iii] $\lambda$3869  and [O ii] $\lambda$3727,29  emission lines. However, [Ne iii] $\lambda$3869  is too faint to derive the resolved metallicity of this galaxy using MUSE data, and we use instead H$\alpha$  and [N ii] $\lambda$6585  emission from SINFONI data.

We cannot rule out the presence of an AGN in any of these three galaxies, since none has all the required emission lines to compare it with widely used criteria such as the BPT diagram. However, as we argued in Patrício et al. (2018), none of these galaxies has [Mg II] emission, and the emission lines are generally narrow, particularly at the centre, which makes the presence of broad-line AGNs unlikely, although not impossible.

2.1 Optical IFU data

The MUSE data and their reduction, were already presented in Patrício et al. (2018) and we provide here only a short summary. AS1063 and A370 were observed for 3.25 and 6 hours, respectively. We used the ESO MUSE reduction pipeline version 1.2 (Weilbacher et al., 2016) with the usual calibrations (bias, flat, illumination and twilight). The pipeline sky subtraction was improved by using the Zurich Atmosphere Purge tool (ZAP version 1; Soto et al. 2016), a principal component analysis that isolates and removes sky line residuals, on the individual data cubes.

To determine the Point Spread Function (PSF), the final cubes were compared with HST data covering the MUSE wavelengths. We assume a Moffat profile, with a fixed power index of 2.8, and fit the Full-Width Half Maximum (FWHM) by minimising the difference between a MUSE pseudo F814W image and the HST F814W image convolved with the Moffat kernel (see Bacon et al. 2017 for details).

2.2 Infrared IFU data

MACS1206-arc was targeted with SINFONI in 2011 with a total exposure time of 6 hours. The SINFONI data were reduced with the pipeline developed by MPE (SPRED, Abuter et al. 2006; Förster Schreiber et al. 2009) together with custom codes for the correction of detector bad columns, cosmic ray removal, OH line suppression and sky subtraction (Davies, 2007) and flux calibration.

The main steps of the procedure are as follows. Master bias and flat images were constructed using calibration cubes taken closest in time to the science frames and used to correct each data cube. The science frames were pair-subtracted with an ON-OFF pattern to eliminate variation in the infra-red sky background. The wavelength calibration is based on the Ar lamp. For each set of observations, a flux standard star was observed at approximately the same time and airmass and was reduced in the same way as the science data. These flux standard stars were then used for flux calibration by fitting a black-body spectrum to the O/B stars or a power law to the cold stars and normalising them to the 2MASS magnitudes. These spectra were also used to remove atmospheric absorption features from the science cubes. The different observations were then combined spatially by using HST images with a larger field of view and good astrometry taken in a similar band as the SINFONI cube, and aligning the SINFONI cube relative to that image. Given that the lenses have such distinctive morphologies, this technique provides reliable coordinates. After these steps, voxels (3D pixels) with flux levels more than 7 standard deviations from the median of the neighbouring voxels were rejected, using a sigma clip algorithm.

At this point, we inspected the quality of the data. In Fig. 1 we present the H$\alpha$  pseudo-narrow band image obtained from the M1206-sys1 data cube by integrating the flux in a spectral window of 12 pixels centred on H$\alpha$  (which corresponds to 3$\sigma$, assuming a Gaussian shape for the H$\alpha$  line profile). The continuum was estimated from two close spectral windows of 6 Å  width each and subtracted from this pseudo-narrow band. Beside M1206-sys1, other two highly magnified $z\approx 1$ galaxies from the sample of Patrício et al. (2018) have been observed with SINFONI: A2390-arc and A521-sys1. However, only M1206-sys1 is bright enough to derive metallicity maps.

Finally, we adjusted the flux calibration and determined the PSF of the SINFONI M1206-sys1 data by comparing a SINFONI F125W pseudo-broad band image with the HST F125W band. The SINFONI field of view is too small to apply the same procedure of image convolution as done with MUSE data, so we fit the two cluster members visible in the SINFONI data. We assume a 2D Moffat profile and, using the astropy package (Astropy Collaboration et al., 2013), fit the cluster members both in HST and in the SINFONI F125W pseudo-broad band image. We then measure the photometry in both images in the same aperture, subtracting the background noise. We find that our nominally reduced SINFONI data overestimate the flux by $\approx 11\%$ when compared to HST and we correct the SINFONI data for this offset. We obtained a PSF FWHM of 0.75" and 0.80" for each cluster member, and we take the mean as the seeing of the SINFONI data throughout this work.

3 Deriving Metallicity

3.1 Data binning and spectral extraction

We start by producing a white light image (the sum along the wavelength axis of the data cubes) for AS1063-arc and A370-sys1 and bin these images using the Cappellari & Copin (2003) method of Voronoi tessellation. We opt to use the white light image as opposed to the [O ii] $\lambda$3727,29  pseudo-narrow band because, in the case of AS1063-arc, using this pseudo-band resulted in a tessellation highly biased towards the strong H ii  south region. For M1206-sys1, due to the higher noise in the SINFONI cube and the fact that we do not detect significant continuum, we use the H$\alpha$ pseudo-narrow band.

We choose a low, arbitrary target signal-to-noise ratio for the tessellation, extract the spectrum from each resulting Voronoi bin, and measure the emission line fluxes in each spectrum (details in the following sub-section). We then check the signal-to-noise ratio of the emission lines in each Voronoi bin. We repeat the process, increasing the target signal-to-noise ratio of the tessellation, until we obtain a signal-to-noise ratio of at least 3 in all bins and for all emission lines. Once this condition is met, we check the quality of the fits of the emission lines for each bin and reject problematic fits. We use the fluxes measured in each Voronoi bin to derive metallicity, dust attenuation and dust-corrected SFRs.

We check the size of the final bins by summing the number of pixels of each bin and converting this area to physical pc2, using the local value of the magnification to correct for lensing magnification. We then take the square-root of this area as an approximation of the size of the bin and plot this in Fig. 2. Most of the resulting bins have sizes smaller than 1 kpc.

3.2 Emission line measurements

The first step in measuring the emission line fluxes is to subtract the continuum, which is especially important for the Balmer lines, since the absorption features are quite significant in these galaxies. We make use of the pPXF routine (Cappellari 2017, version 6.0.2) and a sample of stellar spectra from the Indo-US library (Valdes et al., 2004). The continuum fit is performed masking emission lines. To improve the fit, we add a low-order polynomial to the templates and multiply by a first order polynomial.

After this, the continuum is subtracted from the spectrum and the emission lines are measured using the Automated Line Fitting Algorithm (alfa) from Wesson (2016). Comparing results obtained using alfa and the method of Patrício et al. (2018), we obtained flux differences of less then 8% in the integrated spectra for fainter Balmer lines (H$\delta$  and H7) and less than 1% for strong emission as [O ii] $\lambda$3727,29.

We present the resulting emission line maps, as well as the maps of the ratios used to derive metallicity in the following subsection, in Appendix A.

3.3 Determining metallicity, SFR surface density, and extinction

We use the following diagnostics to derive the metallicity in our sample:

O2

[O ii] $\lambda$3727,29/H$\beta$

O3

[O iii] $\lambda$5007/H$\beta$

O32

[O ii] $\lambda$3727,29/[O iii] $\lambda$5007

R23

([O ii] $\lambda$3727,29+[O iii] $\lambda$4959+[O iii] $\lambda$5007)/H$\beta$

N2

[N ii] $\lambda$6585/H$\alpha$

with the O2, O3, O32, R23 diagnostics and H$\beta$/H$\gamma$  being used for AS1063-arc and A370-sys1, and N2 for M1206-sys1. We use the Maiolino et al. (2008) strong line calibration to derive metallicities from these line ratios. Since H$\beta$  is not available for M1206-sys1, the O2 diagnostic was derived by extrapolating the H$\beta$  flux using the intrinsic H$\gamma$  (i.e. dust corrected, see below) from the MUSE data, assuming the H$\beta$/H$\gamma$  ratio of 2.135, for T${}_{e}=10000$ K and low electron density and case B recombination (Storey & Hummer, 1995).

We make use of the H$\beta$/H$\gamma$  ratio to derive the redenning correction in the case of AS1063-arc and A370-sys1. For M1206-sys1, no correction is applied to the N2 ratio, due to the large uncertainties when combining MUSE (H$\beta$, H$\gamma$) with SINFONI data (H$\alpha$). Moreover, the proximity of H$\alpha$  and [N ii] $\lambda$6585  makes the differential dust attenuation between these two lines small enough that it is reliable to derive metallicities without including dust-correction.

We do not correct for Galactic extinction. This correction would be very small in the case of AS1063-arc and A370-sys1 ($E(B-V)$ = 0.012 mag and 0.032 mag respectively), and with a variation of less 0.001 mag within the full length of the gravitational arcs (Schlafly & Finkbeiner, 2011). For M1206-sys1, the galactic extinction is higher ($E(B-V)$ = 0.063 mag), but for the reason mentioned above, we do not apply this correction either.

Finally, we derive metallicity ($Z$) and attenuation (E(B-V)) from several emission line ratios ($r$) in a Bayesian framework, fitting multiple strong line metallicity diagnostics and extinction simultaneously. We use the emcee Markov chain Monte Carlo Sampler (Foreman-Mackey et al., 2013) to maximise the following Gaussian ($\log$-)likelihood function:

[TABLE]

where ${\rm O_{r}}(E(B-V))$ are the observed line ratios corrected for attenuation and ${\rm M_{r}}(Z)$ are the respective expected ratios, obtained from the Maiolino et al. (2008) calibrations. $\sigma_{r}^{2}$ is the quadratic sum of the observed error and an additional model uncertainties. We adopt an uncertainty of 0.1 dex for the metallicity calibrations and a 1% uncertainty for the case B Balmer line ratios. We use a wide flat prior on metallicity, between 7.0<\mbox{12+\log(O/H)}<9.2, the range of the data analysed in Maiolino et al. (2008) (see their figure 5), and a wide flat prior on attenuation of $0<$E(B-V)$<1$ mag.

The star-formation rates densities are calculated by taking the H$\beta$  intrinsic fluxes and calculating the expected intrinsic H$\alpha$  fluxes, assuming case B, a temperature of T=10 000K and low electron density, and applying the Kennicutt (1998) calibration. Since we calculate SFR densities, no magnification corrections are needed because gravitational lensing conserves surface density brightness (the increased flux due to lensing covers a larger area). The dust attenuation also does not depend on lensing correction, since it is derived from line ratios of each pixel.

We adopt this Bayesian approach as a systematic way to combine different indicators, which has the advantage of having a self-consistent dust and metallicity treatment. However, we do not claim that this will necessarily yield statistically meaningful uncertainty estimates, since the line ratios used in the likelihood function are not independent from each other.

In order to estimate uncertainties in an alternative way, we compute the metallicity using each diagnostic independently and calculate the dispersion of values obtained for each bin (see Appendix A). We did not include any dust correction in these calculations. For AS1063-arc, we obtain a mean standard deviation between metallicity values of 0.09 dex and a maximum dispersion of 0.24 dex, compared with a mean and maximum of 0.03 dex and 0.04 dex obtained using our Bayesian approach. For A370-sys1, we obtain a mean and maximum dispersion of 0.05 dex and 0.12 dex from the individual diganostics, compared with 0.03 dex and 0.07 dex from the Bayesian likelihood maximisation. We notice that amongst the four diagnostics included – R23, O3, O2 and O32 – the latter is the one that most deviates from the mean for both galaxies.

Since the O32 ratio is sensitive to the ionisation parameter (e.g. Kewley & Dopita, 2002), it is possible that differences in local ionisation parameter are driving the dispersion in metallicity. For AS1063-arc, this diagnostic deviates most from the metallicities calculated with the other 3 diagnostics in the H ii  south region, where the highest SFR densities are also found (see Fig. 3) and the highest ionisation parameters is expected due to recent star-formation, which seems to further confirm this hypothesis. It is worth noticing however that the relation between SFR and ionisation parameter is not fully established. For example, Paalvast et al. (2018) do not find a relation between sSFR and the O32 ratio. Furthermore, Shirazi et al. (2014) suggest that high-z galaxies with elevated O32 ratios have high electron densities, not necessarily higher ionisation parameters.

3.4 Metallicity Maps

For both AS1063-arc and A370-sys1, [O ii] $\lambda$3727,29, [O iii] $\lambda$4959, [O iii] $\lambda$5007, H$\gamma$, and H$\beta$  can be measured with a signal-to-noise of at least 3 in each bin. We derive the metallicity and extinction maps using a total of 6 line ratios (O2, O3, O32, R23, [O iii] $\lambda$5007/4959, H$\beta$/H$\gamma$). Weaker lines, such as [Ne iii] $\lambda$3869  and H7 are also well detected in the integrated spectra of these galaxies, but cannot be used to derive resolved properties due to their faintness. We present the comparison between the metallicities derived using different line sets, with and without these fainter lines, in Appendix B. The resolved maps of metallicity, SFR and extinction for these two galaxies are shown in Fig. 3 and 4.

For M1206-sys1, we derive the metallicity using the [N ii] $\lambda$6585/H$\alpha$  ratio and not do include dust-correction. We show the metallicity map of M1206-sys1 in Fig. 5.

We did not account for Diffuse Ionised Gas (DIG) in this analysis, which might impact the values of metallicity. In their sample of local disc galaxies with resolutions of $\approx 100$ pc, Erroz-Ferrer et al. (2019) found that the DIG regions have metallicity on average 0.1 dex lower than the H ii  regions, so we might assume our metallicity values to be upper limits. On the other hand, Erroz-Ferrer et al. (2019) found that the radial gradient of both metallicities (H ii  regions and DIG) were similar, so this caveat in our analysis might not impact the derived gradients, if this result is also valid at $z\approx 1$.

4 Resolved Scaling Relations

We start our analysis by checking whether the resolved Star-Forming Main Sequence (rSFMS) and the resolved Mass-Metallicity Relation (rMZR) are in place for these galaxies. Since these relations only involve surface density quantities (or metallicity) which are conserved by gravitational lensing, we can investigate these correlations regardless of lensing correction.

4.1 Resolved Mass-SFR relation

We derive stellar mass surface densities ($\Sigma_{\star}$) by measuring the photometry in multiple HST bands (F105W, F110W, F125W, F140W, F160W, F435W, F606W, F625W, F775W, F814W and F850W) for each bin defined in the MUSE data. We then use FAST222We included both spectra and photometry to derive the total mass using prospector in our previous work. Since in this case we only use photometry (the spectral continuum signal-to-noise ratio is too low to further constrain the fit), we opted to use FAST, since the computational time to derive masses is substantially smaller. (Kriek et al., 2009), with the Bruzual & Charlot (2003) stellar synthesis models, the Chabrier (2003) IMF and an exponentially decaying star-forming history, and a Calzetti et al. (2000) dust attenuation law. We convert the output masses into mass surface densities, which, as stated before, is independent of lensing.

Using these mass densities and the star-formation rate densities derived from the H$\beta$ lines ($\Sigma_{SFR,H\beta}$) for AS1063-arc and A370-sys1, we plot the rSFMS in the first row of Fig. 6. For MACS1206-sys1, we use the flux of H$\alpha$ as a proxy for SFR, although this is merely indicative.

We fit the rSFMS using a hierarchical Bayesian model, linmix333https://linmix.readthedocs.io (Kelly, 2007), that fits a linear model taking into account uncertainties in both variables involved in the relation. We fit a linear model in the form $\log_{10}\Sigma_{SFR}=a+b\,\log_{10}\left(\frac{\Sigma_{\star}}{2.0}\right)$ for AS1063-arc and A370-sys1 and \log_{10}\mbox{H\alpha}=a+b\,\log_{10}\left(\frac{\Sigma_{\star}}{2.5}\right) for M1206-sys1, placing the pivot point of the linear relation in the middle of the data. We plot the resulting fits in Fig. 6. We obtain slopes of b = 1.03${}^{+0.32}_{-0.20}$ and 1.08${}^{+0.56}_{-0.18}$ for AS1063-arc and A370-sys1, confirming that the SFMS is locally present in these two galaxies. These uncertainties were calculated by taking values of the slope and intercept from several steps of the limix MCMC chain, and placing them in histograms. Since some of the resulting distributions are asymmetric, we take the upper and lower errors as the minimum and maximum values that contain 68% of the values centred in the maximum of the histogram. For a Gaussian distribution, this corresponds to the 1$\sigma$ error.

These slopes agree, within uncertainties, with what was obtained by Wuyts et al. (2013) using 473 massive star-forming galaxies at $0.7<z<1.5$ at kilo-parsec resolutions (slope of 0.95, in yellow dotted line in Fig. 6). In a recent work, Abdurro’uf & Akiyama (2018) also analysed the rSFMS at 1 kpc resolution for massive disc galaxies at $0.8<z<1.8$ (slope of 0.88, in green dotted line), calculating SFRs from broad band SED fitting, finding similar results to Wuyts et al. (2013) and the ones derived here.

For M1206-sys1, there seems to be no correlation between the mass density and the H$\alpha$ flux (the slope is compatible with zero), which might be an indication that the dust attenuation is not the same in the entire galaxy.

4.2 Resolved Mass-Metallicity relation

We plot the metallicity derived for each bin and the corresponding stellar surface density masses in the middle row of Fig. 6 in order to study the rMZR. For AS1063-arc and A370-sys1, we find that metallicity and stellar mass density are correlated, with higher density bins having higher metallicities. Although at lower redshifts (and with substantial more data) this relation is fit with a more complex function, given the small range explored by our data (2 orders of magnitude in $\Sigma_{\star}$), we fit the relation with a linear model, as done for the rSFMS.

We obtain different slopes of 0.28${}^{+0.04}_{-0.02}$, 0.38${}^{+0.15}_{-0.08}$ for AS1063-arc and A370-sys1, which are compatible within uncertainties. For M1206-sys1, we obtain a slope of 0.07${}^{+0.45}_{-0.21}$ between metallicity and stellar mass density, which indicates a very weak relation between these two quantities.

We also plot the relations obtained in the local Universe using the PINGS (Rosales-Ortega et al., 2012) in red, CALIFA (Sánchez et al., 2013) in orange (we use their equation (1) with the parameters a = 8.74, b = 0.018, c = 3.05 (Sánchez, private com., also used in Barrera-Ballesteros et al. (2016)), and MaNGA (Barrera-Ballesteros et al., 2016) in green, using more complex functional forms to fit this relation. Our data points generally fall above all these local relations, i.e. they all have higher metallicities for the same mass surface density than what is found in the local Universe.

However, since determining absolute calibrations is still challenging, it is difficult to directly compare results obtained in different works. Both Rosales-Ortega et al. (2012) and Sánchez et al. (2013) use the O3N2 ratio (([O iii] $\lambda$4959/H$\beta$) / ([N ii] $\lambda$6585/H$\alpha$)) and the calibrations of Pettini & Pagel (2004) (PP04), while Barrera-Ballesteros et al. (2016) use the same ratio but with the Marino et al. (2013) (O3N2-M13) calibrations. Sánchez et al. (2017) investigate these differences calculating the MRZ using different metallicity calibrators for a sample of 613 galaxies observed in the CALIFA survey, obtaining for the same mass, differences of up to 0.4 dex between calibrations. In this analysis, the O3N2 calibrations of Pettini & Pagel (2004) and Marino et al. (2013) are included as well as the R23 from Maiolino et al. (2008) (M08), that we will take as a good approximation to the results derived here combining R23, O3, O2 and O32.

The O3N2-M13 calibration gives results up to 0.2 dex lower than the ones with M18 (figure 3 of Sánchez et al. (2017)), which might explain why the rMZR of Rosales-Ortega et al. (2012) (in red in Fig. 6) predict lower metallicities for the mass-densities analysed here (see Fig. 6). However, the PP04 calibrations give similar results as the M08 calibrations ($\approx 0.02$ dex), while the results from Barrera-Ballesteros et al. (2016) (in orange in Fig. 6) are the ones that most deviate from our results.

It is then difficult to say with certainty if there is an evolution with redshift of the rMRZ or if the discrepancies seen here arise due to the differences in metallicity calibrations.

Trayford & Schaye (2018) used the EAGLE simulation to study the evolution of the rMZR with redshift. They find a strong evolution in the shape of this relation when AGN feedback is included, while it remains fairly similar from $z=0.1$ to 2 when no AGN are present. However, even in this last case, the normalisation (i.e. intercept) of the rMZR shows a strong evolution of about 0.4 dex for stellar mass densities of 102 M⋆/pc2, with higher-$z$ having lower metallicity values.

For the same range of stellar-mass densities studied here, we find metallicity values that are $\approx 0.4-0.5$ dex higher than predicted by Trayford & Schaye (2018) for $z=0.5$ and 1. As for the observational studies, it is not clear if this difference is driven by the choice of metallicity calibration.

4.3 Resolved Fundamental Mass-Metallicity relation

Finally, we investigate the correlation between the residuals of the rSFMS and rMRZ. We plot this in the lower panel of Fig. 6 and calculate the Spearman correlation test for these two residuals.

For AS1063-arc we measure a correlation of $\rho=0.19$ (with corresponding $p$ value of 0.027), corresponding to a weak correlation. For A370-sys1 we obtain a strong correlation of $\rho=0.67$ ($p<$0.0001). For M1206-sys we compare the residuals of the rMZR with the ones from the stellar mass density vs H$\alpha$  flux, that we denoted as r($\Sigma_{\star}$-H$\alpha$), and found no clear correlation between these residuals ($p=0.537$).

Excluding M1206-sys1 from the analysis, given our lack of $\Sigma_{SFR}$ for this galaxy, we measure a positive correlation between $\Delta$rSFMS and $\Delta$rMZR for AS1063-arc and A370-sys1, although weak in the case of AS1063-arc. This might indicate that a relation between resolved $\Sigma_{\star}$, $\Sigma_{SFR}$ and Z is present at higher$-z$. However, given the different values we obtained for this correlation in these two galaxies of comparable mass and metallicities, it might indicate that this relation is not fundamental, in the sense that it is not the same for all galaxies at all redshifts.

We notice also that we find a positive correlation between the two residuals, with higher residuals in rSFMS corresponding to higher residuals in rMRZ, instead of the negative correlation between SFRs and metallicity, for fixed stellar mass, found in other works (e.g. Lara-López et al., 2010; Mannucci et al., 2010), with higher residuals in rSFMS corresponding to higher residuals in rMRZ. However, we base these conclusions in only two objects, and a larger sample with wider redshift range is needed in order to confirm these results.

5 Metallicity Gradient

We now turn our attention to the metallicity distribution within each galaxy, deriving its gradient and inspecting the residuals. We start by describing how gravitational lensing affects the galaxy properties, and proceed to describe how we model the data with a simple 2D radial gradient, taking into account lensing and seeing effects with a forward-modelling approach.

5.1 Lensing distortion

AS1063-arc is the least magnified galaxy, with a mean magnification of $\mu=4$, and also only a small distortion. Using Lenstool (Jullo et al., 2007) and the respective lensing model, we can reconstruct the morphology in source plane, i.e. corrected for lensing magnification (see Fig 14). This process does not account for seeing effects, and the PSF in the source plane is not circular, with a smaller FWHM in the direction where the galaxy is more magnified, where effectively we can probe smaller spatial scales (see the second panel in Fig 14, in appendix). This means that spatial resolution is not homogeneous in this galaxy, which we will explore in the next section.

A370-sys1 and M1206-sys1 have higher magnification factors, reaching $\mu=30$ in some regions, and more complex lensed morphologies, with multiple images of the same regions, which makes the reconstruction process more challenging. The lensed image of A370-sys1 contains one complete image of the galaxy, plus 3 other partial images, i.e., only a portion of the galaxy was imaged into those multiple images. This is also the case for M1206-sys1, where 4 multiple images can be seen in the SINFONI data. However, unlike A370-sys1, the SINFONI data do not contain the full image, and only about half of the disc is available.

Each of these multiple images can be traced back to the source plane using the lensing models. However, this leads to different PSFs in the source plane, since their lensing distortions are different. For AS1063-arc, the FWHM of the PSF measures 2.3 kpc in the direction of highest magnification and 5.69 kpc in the lowest and for A370-sys1, between 0.73 and 3.10 kpc. This means that combining several multiple images in the source plane, without including seeing deconvolution, can produce misleading results. Strategies to deal with this issue have been developed (Sharma et al., 2018), but here we choose a simpler approach, and perform most of our analysis in the image plane, keeping the multiple images separated.

5.2 Forward-modelling metallicity gradients

In order to fully use the spatial information provided by the IFU observations, we fit the metallicity maps assuming a simple 2D axisymmetric gradient, where the metallicity depends on the deprojected galactocentric distance to the centre of the galaxy (corrected for inclination and lensing), the assumed gradient ($\nabla Z$) and central metallicity value ($Z_{0}$).

We build our gradient model in the source plane, calculate the lensing distortions using the lensing models, convolve the lensed gradient with the instrument seeing, and finally compare it with the data, minimising the difference between the two. This approach is similar to the one presented in Carton et al. (2017) for field galaxies, but includes the lensing correction.

We start by producing a deprojected galactocentric distance 2D map in the source plane, using the centre of the galaxy (cx, cy), the ratio between the minor and the major axis (q), and the position angle ($\theta$). Using the lensing model, we forward-lense this deprojected galactocentric distance map to the image plane and align it with the data, rescaling the pixel sizes to match the IFU observations. We then multiply this map by the gradient and add the central metallicity value (Z(x,y) = Z0 + $\nabla$ Z r) to produce a metallicity gradient in the source plane. Finally, we convolve it with the seeing, and apply the same binning as used to derive the metallicity maps. We compare this model gradient with the measured metallicity maps using a Gaussian log-likelihood function and the emcee sampler to maximise the likelihood and obtain the best-fit parameters (as done Section 3.3). We have made this method publicly available444The code, FRApy, for Fitting Resolved Arcs with Python, is available at https://frapy.readthedocs.io..

5.2.1 AS1063-arc

We start by producing a source plane image of the F160W HST band and fit it with galfit (Peng et al., 2010), in order to assess what values the morphological parameters of the metallicity gradient model – q, $\theta$ and centre – could have. We used a global Sérsic profile plus two more compact components for the bulge and the large H ii  southern region. We report the relevant results of the fit in Table 2.

We fit the data first keeping q, $\theta$ and the centre fixed to the values obtained with galfit, and then letting them vary within large intervals. The morphological parameters obtained in the second case are very different from what was obtained with galfit. The centre is offset by about 0.4 arcsec and the position angle $\theta$ differs by $\approx 90$ degrees. This difference arises from the fact that the two spiral arms (and the major axis of the galaxy derived with galfit) are aligned with the direction of the highest stretch caused by gravitational lensing, that together with the poor seeing at which this galaxy was observed ($\approx 1$"), makes it challenging to derive the correct morphology.

Following these two approaches, we obtain gradients of -0.034$\pm$0.002 and -0.042$\pm$0.002 dex/kpc, respectively and central metallicities (8.99$\pm$0.01 and 9.04$\pm$0.01 in 12+log(O/H)). We plot the 1D profiles for both these fits in Fig. 7.

5.2.2 A370-sys1

We fit the A370-sys1 metallicity map with the same technique, starting by fitting the morphology using the F160W HST band. Due to the difficulties in combining different multiple images (see subsection  5.1), we use only the complete multiple image to perform the galfit fit. Since this galaxy also has a complex morphology, we use several components in the fit (disc, bulge plus strong star-forming regions), and report the values for the disc in Table 2.

We then proceed to fit A370-sys1 fixing the morphology to the values found with galfit and also letting q, $\theta$ and the central position free. The results are listed in Table 2. In this case, we obtain axis ratios and $\theta$ closer to what was obtained with galfit, but still inconsistent with this method.

The central metallicities obtained in both fits are also close (8.98$\pm$0.01 and 9.03$\pm$0.01 in 12+log(O/H)), and although not formally compatible, they are well within the typical uncertainty of metallicity calibrations. We also obtain different gradients, -0.039$\pm$0.004 and -0.053$\pm$0.004 dex/kpc, respectively. The 1D profiles obtained with both fits are shown in Fig. 8.

5.2.3 M1206-sys1

Because of the complexity of the lens model and the low(er) number of metallicity measurements, which do not allow us to reliably constrain the parameters of the metallicity gradient model, we performed only a simple 1D analysis for M1206-sys1.

We produce a source plane deprojected distance map, using the ellipticity, position angle and centre from an elliptical fit to the F160W HST image of the complete multiple image of the galaxy. We forward-lens this map using lenstool, and define 1 kpc annular apertures starting at $r=0$, measuring the average metallicity in these annuli. This approach does not include any correction for seeing, which it is known to flatten gradients (Yuan et al., 2013). We fit the data with the linmix555https://linmix.readthedocs.io/en/latest/ package. We obtain a slope of -0.039$\pm$0.060 dex/kpc, a central metallicity 9.06$\pm$0.25 in 12+log(O/H). The data and fit are shown in Fig. 9.

5.3 Comparison with the literature

At high-redshift, a wide range of metallicity gradients have been derived from lensing studies, which range from quite steep negative gradients (e.g. Jones et al., 2013; Wang et al., 2017) to positive gradients (e.g. Leethochawalit et al., 2016) that are usually not observed in the local Universe. However, these previous lensing studies focused on galaxies at considerably higher redshifts (1.2 $\leq z\leq$ 2.3) than the three objects analysed here ($z=0.6$, 0.7 and 1.0).

A better match in redshift to our sample are the Wuyts et al. (2016) and Carton et al. (2018) surveys of field galaxies. Wuyts et al. (2016) analyse a sample of 180 star-forming galaxies from the KMOS3D survey, from $z=0.6$ to 2.7, with stellar masses between 109.5 and 1011.5 M⊙ and SFR between 0.1 and 1000 M$\odot$/yr , measuring the metallicity in annuli using the N2 indicator. Most of their sample have flat gradients, with only $\approx 7\%$ of the sample exhibiting positive gradients. Carton et al. (2018) analyse a sample of 84 galaxies from several MUSE GTO programmes, with stellar masses between $10^{7}$ and $10^{10.5}$ M⊙ and SFR between 0.01 and 10 M⊙/yr at $z=0.2-0.8$, combining several metallicity diagnostics in a 2D forward-modelling approach. They obtain a mean negative gradient of -0.039${}^{+0.007}_{−0.009}$ dex/kpc, but with a larger spread in gradients than found by Wuyts et al. (2016). AS1063-arc and A370-sys1, with redshifts of 0.6 and 0.7, are at the intersection of these two studies, and are compatible with the mean values of both. We compare M1206-sys1, at $z=1$, only with Wuyts et al. (2016). We obtain a gradient more negative than most galaxies between $z=0.9-1.1$ (-0.006 dex/kpc) , but still compatible with Wuyts et al. (2016) within uncertainty.

There are strong indications for the existence of a characteristic metallicity slope in low-$z$ galaxies, when the physical slope (dex/kpc) is normalised to the size of the galaxies. Both Sánchez et al. (2014) and Sánchez-Menguiano et al. (2018) find a characteristic (scaled) slope of -0.1 dex/Re, when the gradient is normalised to the effective radius Re (see also Ho et al. 2015 for a R25 normalisation). At higher redshift, Carton et al. (2018), find a steeper slope of -0.34 dex/Re (for galaxies with Rd>3 kpc, as the ones presented here, and converting Rd in Re), although with a higher spread than found at lower redshift ($\sigma_{\mathrm{int}}$ = 0.1 dex).

We normalise the gradients with the values of Re obtained from morphological fits (see Table 1), obtaining $\nabla\,Z$ of -0.323$\pm$0.007, -0.636$\pm$0.011 and -0.407$\pm$0.658 dex/Re for AS1063-arc, A370-sys1 and M1206-sys1, respectively. These are all significantly steeper scaled gradients than what is found for low redshift galaxies (-0.1 dex/Re), or for galaxies between 0.1$\leq z\leq$0.8 as in Carton et al. (2018). Part of the discrepancy might be explained by errors in the Re, derived using galfit.

5.4 Deviations from radial gradients

Here we analyse the residuals of the metallicity maps after subtracting the fitted gradients, which we refer to as metallicity residuals.

For AS1063-arc, when using the morphological parameters obtained with galfit, the radial residuals are as high as 0.1 dex, when radially binned in 0.5 dex metallicity bins, but without a clear radial trend (see the bottom panel of Fig. 7). For the fit where all variables are allowed to vary, the residuals are very low ($\leq$0.02 dex) up until 6 kpc ($\sim 0.8R_{e}$). After this, there seems to be a trend of increasing residuals with radius. This could be caused by a flattening of the metallicity gradient at outer radii (between 0.5 to 3 Re), as observed in some cases in the local Universe (Sánchez-Menguiano et al., 2018), but it would be necessary to probe the metallicity gradient further out in order to confirm this.

As for A370-sys1, both models, with fixed or free morphological parameters, result in residuals of about $\leq$0.05 dex, when the data is radially binned in bins of 0.5 dex.

In the 2D analysis of the metallicity residuals, we consider only the gradient modelled with free parameters, for simplicity. In Fig. 10 we plot the 2D residuals, as well as the residuals versus the stellar mass surface density and star-formation density. We do not see any trend with morphological features of the galaxies. We note that Erroz-Ferrer et al. (2019) in their analysis of local discs, found a metallicity increase of about $\approx 0.2-0.25$ dex in H ii  regions when compared with the surrounding metallicity. This does not appear to be the case for these z$\sim$1 galaxies, despite the fact that they do contain giant H ii  regions, typical of high-$z$ disc galaxies.

We investigate this further by plotting the residual metallicity versus the star-formation density, also in Fig. 10, and computing the Spearman rank correlation coefficient between these two quantities. We obtain values of $\rho=-0.1$ and -0.07, with $p$ values of 0.24 and 0.48, showing no clear correlation between the residual metallicity and the star-forming rates densities. One possible explanation for not observing the same increase in metallicity as noted in Erroz-Ferrer et al. (2019), is the difference in spatial scales probed. Although the work presented here probes sub-kiloparsec regions, which are at $z\approx 1$ only possible to study in lensed galaxies, Erroz-Ferrer et al. (2019) observe galaxies at $<$100 pc scales, an order of magnitude smaller.

6 Summary and Conclusions

In this work, we made use of HST, MUSE and SINFONI data to analyse the spatially resolved properties of 3 lensed galaxies at redshifts 0.6, 0.7 and 1, at exceptionally high spatial resolution (see Fig. 2). We derive the stellar-mass surface density using multiple HST bands. For the two lower-redshift targets, AS1063-arc and A370-sys1, we derive the gas metallicity using the line ratios (O2, O3, O32, R23, [O iii] $\lambda$5007/4959, H$\beta$/H$\gamma$) and the Maiolino et al. (2008) metallicity calibration. For M1206-sys1 only N2 was available. Using these results, we examine the resolved Star-Forming Main Sequence (rSFMS) at $z\approx 1$ at sub-kiloparsec resolution, at a physical scale unattainable with un-lensed galaxies. We also explore, for the first time at $z\approx 1$, the resolved Mass-Metallicity Relation (rMZR) and the resolved Fundamental Mass-Metallicity Relation (rFMZ).

In order to fit the 2D metallicity gradients, we develop a forward-modelling method that fits data in the image plane, correcting for seeing and lensing distortions, avoiding issues arising from combining different multiple images.

Our main results from this analysis are:

• •

We find that both the rSFMS and rMZR are in place for galaxies AS1063-arc ($z=0.6$) and A370-sys1 ($z=0.7$), although with different slopes as the ones observed in the local Universe (Fig. 6).

• •

For these two galaxies, we also find a correlation between the residuals of the rSFMS and the rMZR ($\rho$=0.19 and 0.67, Fig. 6), which might indicate the presence of a rFMZ. We notice however, that we find the opposite correlation (with higher rSFMS residuals corresponding to higher rMZR residuals) to what is found in other works. Moreover, the correlations are different for the two galaxies tested, which suggests that the relation evolves with redshift. A larger sample is needed in order to confirm these results.

• •

We measure metallicity gradients of $-0.027\pm 0.003$, $-0.019\pm 0.003$ and $-0.039\pm 0.060$ dex/kpc for our three targets (Table 2). This is in agreement with what was derived for surveys at similar redshifts.

• •

We find no significant deviations from an exponentially decreasing metallicity gradient (Fig. 7 and 8). In particular, we find no increase or decrease of the metallicity with star-formation rate density (Fig.10). We find a mean dispersion of the metallicity residuals of 0.01 dex for AS1063-arc and of $\approx 0.05$ dex for A370-sys1.

We conclude that, although the galaxies analysed are typical high-$z$ disc galaxies, with several large H ii  regions (clumps) and highly turbulent ionised gas, the relation between stellar mass surface density, star-formation rate surface density and metallicity at sub-kiloparsec scales observed at in local discs is already in place at $z\approx 1$. Moreover, a negative metallicity gradient is already established, although with steeper scaled gradients than seen in local disc galaxies, and there are no significant metallicity deviations from a linear gradient, either due to morphological structures such as spiral arms or star-forming regions.

The data and analysis done for this work can be found in https://github.com/VeraPatricio/Resolved_Metallicity.

Acknowledgements

We thank Nicole Nesvadba for reducing the SINFONI data. We also thank Tiantian Yuan, Lisa Kewley and Ayan Acharyya for useful and insightful comments on how to improve this work. Finally, we thank the referee for helpful suggestions that made this work clearer.

VP is supported by the grant DFF - 4090-00079. DC acknowledges support from the ERC starting grant 336736-CALENDS. CP thanks the Alexander von Humboldt Foundation for the granting of a Bessel Research Award held at MPA. CP is also grateful to the ESO and the DFG Cluster of Excellence "Origin and Structure of the Universe" for support. JB acknowledges support by FCT/MCTES through national funds by this grant UID/FIS/04434/2019 and through the Investigador FCT Contract No. IF/01654/2014/CP1215/CT0003.

This research made use several open source python libraries: numpy (van der Walt et al., 2011), scipy (Jones et al., 01), matplotlib (Hunter, 2007), and astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013). This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in A&AS 143, 23. This work has made use of dust extinction maps from the NASA/IPAC Infrared Science Archive.

Appendix A Emission line and line ratios maps

We present the signal-to-noise ratio maps of the emission lines used to derive metallicity ([O ii] $\lambda$3727,29, H$\beta$, [O iii] $\lambda$4959, [O iii] $\lambda$5007, [N ii] $\lambda$6585  and H$\alpha$) in the top panels of Fig. 11 and 12. For [O ii] $\lambda$3727,29, we plot the sum of the doublet. The signal to noise ratio was calculated using the flux and uncertainties measured with alfa, as detailed in Section 3.

Using these maps, without including any dust correction, we calculate the individual line ratios used in this work (middle rows of Fig. 11 and 12). Using these and the Maiolino et al. (2008) calibrations, we calculate the metallicity maps for each individual diagnostic. We notice that we obtain the largest discrepancies with O32, an ionisation sensitive diagnostic. We also measure the dispersion in metallicity for each bin, calculating the standard deviation in each bin between of all metallicity maps.

For M1206-sys1, since we have only one line ratio available, we present only the signal-to-noise ratio maps of the two lines used (H$\alpha$  and [N ii] $\lambda$6585) and the ratio of the two in Fig. 13.

Appendix B Comparison of metallicity derived from different line sets

We compare the metallicity derived in this work using only the strongest lines, with the one obtained in Patrício et al. (2018) (hereafter P18) from integrated spectra, where faint lines were also included ([Ne iii] $\lambda$3869, H$\gamma$, H$\delta$,  and H7). For M1206-sys1, we also test the consistency of the results derived using MUSE and SINFONI data or just SINFONI data.

Besides all the metallicity dependent rations presented in Section 3.3, we also included here the Ne3O2 ([Ne iii] $\lambda$3869/[O ii] $\lambda$3727,29) ratio and the following metallicity independent ratios:

H$\alpha$****/H$\gamma$

6.113

H$\alpha$****/H$\delta$

11.057

H$\alpha$****/H7

18.004

H$\beta$****/H$\gamma$

2.135

H$\gamma$****/H7

6.288

H$\gamma$****/Hd

1.809

H$\delta$****/H7

1.628

[O III]** $\lambda$5007/4959**

2.98

In P18, 10 line ratios were used to derive the integrated metallicity of AS1063-arc and A370-sys1 (see Table 3). In this work, only 5 ratios (O2, O3, O32, R23, and H$\beta$/H$\gamma$) are available to study the resolved metallicity and we re-derived the integrated metallicity using only those 5 ratios and compare it with the previous values. The new metallicity and extinction are presented in Table 3.

We obtain slightly lower metallicities for AS1073-arc – from 8.82$\pm$0.02 in P18 using 10 line ratios, to 8.75$\pm$0.10 in 12+$\log$(O/H)  – and A370-sys1 – 12+$\log$(O/H) = 8.88$\pm$0.02 in P18 and 8.83$\pm$0.15 in this work – but that are compatible within uncertainty. Indeed, the uncertainty of the metallicities derived in this work are considerably higher (and more realistic) than in P18, reflecting both the use of less constraints and the addition of the continuum subtraction uncertainty to the line flux errors. A similar trend is seen with the values of $\tau_{v}$, the extinction factor obtained with the Charlot & Fall (2000) law, that are higher than in P18. As previously described, the chosen extinction law has a very small impact in the metallicity derived, about 0.01 dex, much smaller than the associated uncertainties. It seems then possible to obtain metallicities comparable as the ones derived using a larger set of line ratios, using only the 6 line ratios involving the strongest lines, although with a higher associated uncertainty.

For M1206-sys1, the MUSE data only covers [O ii] $\lambda$3727,29  and [Ne iii] $\lambda$3869, as well as several weak Balmer lines, H$\gamma$, H$\delta$  and H7. However, using SINFONI, both H$\alpha$  and [N ii] $\lambda$6585  can be utilised. We first start to test whether the N2 diagnostic gives compatible results with the ones presented in P18, using Ne3O2 and O2 (see Table 3). We obtain a metallicity of 8.89$\pm$0.8 in 12+$\log$(O/H), compatible with what was previously derived not including N2. However, the $\tau_{v}$ obtained is quite higher, indicating some possible remaining issues with the flux calibration between MUSE and SINFONI data (we remind the reader that the method used here to determine metallicity uses all lines to determine extinction). Indeed, if we add more line ratios involving H$\alpha$  and other Balmer lines in the MUSE data, the dust attenuation values obtained are clustered around our highest allowed extinction, much higher than what is obtained with only the MUSE data, and surprisingly high (E(V-B)$>1$ mag). We conclude that our flux calibration between the two data sets is not accurate enough to allow to robustly determine the extinction combining H$\alpha$  with other Balmer lines. However, relying only on SINFONI data and the N2 metallicity diagnostic, without any Balmer ratios, we obtain a slightly higher global metallicity (12+$\log$(O/H)= 8.94$\pm$0.07) but that it is still compatible with what is derived using only Ne3O2 and O2. The proximity of H$\alpha$  and [N ii] $\lambda$6585  makes the differential dust attenuation between these two lines small enough that it is still reliable to derive metallicities not including dust correction.

Appendix C 2D Maps in Source Plane

We use lenstool to correct the image plane maps of metallicity, extinction and SFR densities for lensing distortions and plot the results in Figures 14, 15 and  16. For A370-sys1 and M1206-sys1, we reconstruct the different multiple images separately. We can see that in the case of A370-sys1 (Fig. 15) the results from the different multiple images are sightly different, as it is expected since they come from different voxels in the data cube, but show a global agreement, with higher metallicities, E(B-V) and SFRs in the centre of the galaxy. AS1063-arc also displays higher metallicity and E(B-V) values at the centre of the galaxy. However, E(B-V) is also high in the region of higher star-formation rates, at the edge of a spiral arm.