V Patricio
email: patricio@dark-cosmology.dk

J Richard

D Carton

T Contini

B Epinat

J Brinchmann

K B Schmidt

D Krajnovic

N Bouché

P M M Weilbacher

V Patrício

D Krajnović

R Pelló

J Caruana

M Maseda

H Finley

F E Bauer

J Martinez

G Mahler

D Lagattuta

B Clément

G Soucail

L Wisotzki

Kinematics, Turbulence and Star Formation of z ∼1 Strongly Lensed Galaxies seen with MUSE

Keywords: galaxies, kinematics and dynamics -gravitational lensing, strong 1 INTRODUCTION

HAL is

which makes them significantly more massive than giant clouds in local galaxies, with masses of 10 5 -10 6 M . While it was initially speculated that these features had their origin in merging episodes, the advent of Integral Field Spectrographs showed that most clumpy galaxies have smooth velocity fields and are rotationally supported (e.g. Genzel et al. 2006;Bouché et al. 2007;Cresci et al. 2009;Förster Schreiber et al. 2009;Wisnioski et al. 2015;Contini et al. 2016), suggesting that these clumps may be part of the secular evolution of disc galaxies.

The physical picture that explains why rotating discs are more clumpy and turbulent at z ≥1 is still being investigated, both observationally and theoretically. One possible scenario is that the highly turbulent interstellar medium is stirred by radiation pressure and winds of strong star formation taking place in these young galaxies (e.g. Lehnert et al. 2009;Green et al. 2010;Lehnert et al. 2013). Another possibility is that the inflow of gas into the galaxy from the cosmic web provides enough energy to sustain these high velocity dispersions, although some works seem to point that this energy is insufficient to maintain the high-velocity dispersions over long timescales (e.g. Elmegreen & Burkert 2010). Finally, a third scenario is gravitational instability: the high gas fractions of these galaxies make them gravitationally unstable, causing gas to spiral down to the centre of the galaxies, converting gravitational energy into turbulent motions and driving galaxies to stability (e.g. Bournaud et al. 2007;Ceverino et al. 2010;Krumholz & Burkert 2010). Possibly, a combination of factors is at play. Recent simulations by Krumholz et al. 2017, that include both stellar-feedback and gravitational instability (transport driven) turbulence, show that while both processes contribute to the gas turbulence, transport-driven mechanisms dominate, especially for z > 0.5 and more massive galaxies.

From the observational side, numerous surveys have provided some insights into these hypothesis. Recent work by Johnson et al. 2018 studied the velocity dispersion of ∼450 star forming galaxies at z ∼ 0.9 from the KROSS survey, finding that the data are equally well explained by a scenario where turbulence is driven by stellar feedback or increased gas fractions. Previous studies have also not provided a definitive answer. In a analysis of z ∼ 3 galaxy analogues (the DYNAMO sample), Green et al. 2014 found a correlation between the star formation rates and velocity dispersions that supports the idea that turbulence is driven by stellar feedback. On the other hand, at high redshift (z ∼ 2), the analysis of KMOS 3D survey data by Wisnioski et al. 2015 found only a weak correlation between gas velocity dispersion and SFR or gas fractions at each redshift.

Another possible path to understanding the properties of these young discs is to study their resolved properties at sub-kiloparsec scales. This is particularly challenging when trying to derive intrinsic velocity dispersions. Owing to the limited resolution of observations, the velocity gradient present in rotating discs artificially increases the observed velocity dispersion of galaxies. The effects of beam smearing are particularly problematic in the centres of galaxies, where the velocity field rapidly changes. A few works used adaptive optics to measured the resolved intrinsic velocity dispersion maps of high-z disc galaxies, overcoming some of the issues caused by beam-smearing (e.g Genzel et al. 2011;Swinbank et al. 2012a;Newman et al. 2013). Genzel et al. 2011 find that the velocity dispersion maps are broadly compatible with a constant distribution, as seen in local disc galaxies, despite its higher value (σ ∼60 km s -1 compared with 10-20 km s -1 in local galaxies). However, other studies find a relation between local star formation rate surface densities (Σ SFR ) and velocity dispersions, which seems to point to some degree of structure (e.g. Swinbank et al. 2012a;Lehnert et al. 2009).

Studying turbulence in the high-z star-forming clumps could also provide some insights into how they differ from their local counterparts. High spatial resolution is necessary to study both the morphology of the velocity dispersion and the clumps' turbulence. However, many studies of highredshift galaxies are still severely hampered by the relatively low spatial resolution achieved with the current facilities. A possible strategy to improve this is to target lensed galaxies, where even distant z ∼ 1 galaxies can be resolved down to a few hundred parsecs (e.g. Jones et al. 2010;Livermore et al. 2015;Leethochawalit et al. 2016;Yuan et al. 2017).

Here, we target a sample of typical z ∼1 clumpy discs that, due to strong gravitational lensing, appear as extremely extended objects in the sky. Their high magnification factors, as well as the high quality MUSE 1 data acquired with excellent seeing, allows us to alleviate resolution issues, providing a sub-galactic view of the properties of these typical galaxies.

This paper is organized as follows. In Section 2, we present the sample and the MUSE observations used in this work as well as other ancillary data and the lensing models used to recover the intrinsic properties of these galaxies. In Section 3 we derive the integrated physical properties of this sample from the MUSE spectra and HST photometry, and in Section 4 we study the resolved kinematic properties of the galaxies. We conclude with a discussion of the results and comparison with other samples and works in Section 5. Throughout this paper, we adopt a Λ-CDM cosmology with Ω=0.7, Ω m =0.3 and H 0 = 70 km s -1 Mpc -1 . Magnitudes are provided in the AB photometric system (Oke 1974). We adopt a solar metallicity of 12 + log(O/H) = 8.69 (Allende Prieto et al. 2001) and the Chabrier 2003 initial mass function.

OBSERVATIONS AND DATA REDUCTION

In this work, we analyse a sample of eight strongly lensed galaxies that lie behind eight different galaxy clusters (see Table 1). All of these clusters were observed with MUSE: Abell 370 (A370), Abell 2390 (A2390), MACSJ0416.12403 (M0416), MACSJ1206.2-0847 (M1206), Abell 2667 (A2667) and Abell 521 (A521) within the MUSE Guaranteed Time Observations (GTO) Lensing Clusters Programme (PI: Richard); Abell S1063/RXJ2248-4431 (AS1063) during MUSE science verification (PI: Caputi & Clément); and MACSJ1149.5+2223 (M1149) during Director's Discretionary Time (PI: Grillo), targeting the newly discovered supernova in a lensed galaxy [START_REF] Kelly | The Astronomer's Telegram[END_REF].

We base our selection criteria on the apparent (i.e. due to gravitational lensing magnification) size of the targets with the goal of resolving these galaxies at sub-kiloparsec a) see also Karman et al. 2015 (b) see also Grillo et al. 2016 (c) measured in F606W.

scales. Within the sample of lensing clusters observed with MUSE, we select the highly magnified (µ > 3) and exceptionally extended (> 4 arcsec 2 ) objects, generally dubbed gravitational arcs, that allow to probe properties at physical scales down to a few hundred parsec at z = 1. For comparision, without lensing, at this redshift, a good seeing of 0.6" results in a resolution of ∼5 kpc. In order to investigate the physical properties of these gravitational arcs, such as kinematics and star formation rates, we require that strong, non-resonant emission lines are detected. Therefore, we target objects with visible [O ii ]λ3726, 29 emission in the range 4650 and 9300 Å. This limits our sample to eight galaxies between 0.6 < z < 1.5.

For clarity, we refer to these galaxies by the cluster name plus their respective number in the lensing models used in this work (see Table 1), for example "A370-sys1" or "A2667-sys1". Exceptions to this are the objects in cluster AS1063 and A2390, which do not have multiple images and are referred to as "AS1063-arc" and "A2390-arc".

In the following subsections, we discuss the observations and data reduction and briefly summarise the ancillary data available for these targets, four of which are part of the Frontier Fields HST programme (Lotz et al. 2017). We also present the lensing models and magnification factors used throughout the paper to derive the intrinsic properties of the gravitational arcs.

MUSE observations and Data Reduction

The targets were observed with MUSE between 2014 and 2017. MUSE is an Integral Field Unit instrument with a field of view of 1'×1', sampled at 0.2" per pixel, and covers the optical wavelength range between 4650 and 9300 Å with a spectral sampling of 1.25 Å (Bacon et al. 2014). The GTO targets have a variety of integration depths (between 2 and 6h). Each observation is comprised of individual 1800 s exposures. The position angle is rotated by 90 degrees between each exposure, minimising the stripe pattern that arises from the image slicer. In order to further minimise this pattern, a small dither (< 1 ) was added to each pointing. The observing strategy for the non-GTO targets (AS1063, MACS1149 and partially A370) is described in Karman et al. 2015 andGrillo et al. 2016 respectively. They follow very similar strategies to the one used within the GTO sample, including a small dithering and a 90 degrees rotation between exposures. A521 was observed in Oct. 2017 with the newly-commissioned Adaptive Optics Facility.

The data from all eight clusters was reduced using the ESO MUSE reduction pipeline version 1.2 [START_REF] Weilbacher | MUSE-DRP: MUSE Data Reduction Pipeline[END_REF]. The calibration files used to perform bias subtraction and flat fielding, including illumination and twilight exposures, were chosen to be the ones closest to the date of the exposure. To perform flux calibration and telluric correction, all 15 available standard stars were reduced and their flux and telluric response derived using the esorex recipe muse standard. The response curves were visually inspected and extreme responses (e.g. very high or low, or with unusual features compared with others) were removed before producing a median flux calibration and telluric correction with the 6 remaining response curves. These median calibration curves were applied to all individual cubes.

The individual cubes were aligned with HST F814W images using a combination of sextractor (Bertin & Arnouts 1996), to identify bright sources in individual the cubes and HST image, and scamp (Bertin 2006), to calculate the offset between the two.

Sky subtraction was performed using the ESO pipeline and the remaining sky subtraction residuals were removed from each individual exposure using the Zurich Atmosphere Purge tool (ZAP version 1; Soto et al. 2016), a principal component analysis that isolates and removes sky line residuals. Finally, the individual exposures were combined, rejecting voxels more than 3 σ from the median, in order to eliminate cosmic rays. To normalize the exposures, we estimate the sky transparency using the median of the fluxes estimated with sextractor. For each cluster, the exposure with the highest median flux was taken as the reference and the mean flux of the other exposures was rescaled to match this reference, before combining the data. A second sky residual subtraction using ZAP was then performed on the combined cubes, and a median filter along the wavelength axis (box of 100 Å) was applied in order to further smooth the background. During both operations, ZAP and median subtraction, the brightest sources were masked when estimating the background. Finally, Bacon et al. 2014 used the same data reduction and found that the variance propagated by the pipeline is underestimated by a factor of 1.6. We correct the variances in out cubes by this factor.

To determine the PSF size, the final cubes were compared to the available HST data. By weighting the MUSE cubes with the F814W filter transmission curve, we produced MUSE F814W pseudo broad band images. We model the PSF following the Bacon et al. 2017 method, which consists of assuming a Moffat profile, with a fixed power index of 2.8, and fitting the Full-Width Half Maximum (FWHM) by minimising the difference between the MUSE pseudo broad band and the HST F814W image convolved with the Moffat kernel (F606W was used for A521). With the exception of AS1063, all targets were observed in excellent seeing conditions, always better than a FWHM of 0.70 (see Table 1).

Ancillary Data and Lensing Models

The gravitational arcs used in this work were selected to be strongly magnified and are either multiply-imaged or very close to the multiple-image region near the cluster core. In order to recover their intrinsic (i.e. corrected for the lensing magnification) properties and morphology, accurate lensing mass models are required.

The selected clusters have all been well studied in the past using numerous and deep HST observations. All clusters except A370, A2390 and A521 are part of the Cluster Lensing And Supernova survey with Hubble (CLASH, Postman et al. 2012), for which photometry is available in 12 filters (see Table B1 for full list). Furthermore, four of them -AS1063, A370, MACS0416 and MACS1149 -are part of the Frontier Fields initiative (Lotz et al. 2017), that provided deep imaging in the F435W, F606W, F814W, F105W, F125W, F140W and F160W bands.

Only A2390 and A521 were not part of any of these programmes. For A2390, abundant HST data is also available (see Richard et al. 2008;Olmstead et al. 2014, as well as Table B1). For A521 only WFPC2/F606W band is available plus a NIRC2/Ks band ground-based image retrieved from the Keck archive. This wealth of data has allowed a large number of multiple systems to be identified (between 5 and 60, creating up to 200 multiple images, e.g. Caminha et al. 2017) resulting in very well-constrained lensing models. The MUSE-GTO programme confirmed, as well as discovered, a significant number of these systems, improving the mass models in the process (Richard et al. 2015;Bina et al. 2016;Lagattuta et al. 2017;Mahler et al. 2018).

The mass models used in this work were all constructed using the Lenstool2 software (Jullo et al. 2007), following the same methodology described in detail in earlier works (e.g. Richard et al. 2009). To summarise, we model the 2Dprojected mass distribution of the cluster as a parametric combination of cluster-scale and galaxy-scale dPIE (double Pseudo Isothermal Elliptical, Elíasdóttir et al. 2007) potentials. To limit the number of parameters in the model, the centres and shapes of the galaxy-scale components are tied to the centroid, ellipticity and position angle of cluster members as measured on the HST images. The selection of cluster members is performed using the red sequence in a colour-magnitude diagram (see e.g. Richard et al. 2014) and confirmed with MUSE spectroscopy. Constraints on the lens model are robustly identified multiple images selected based on their morphology, colours and/or spectroscopic redshifts. Although previous Lenstool models produced by our team were published for some of the clusters (e.g. Richard et al. 2010bfor A2390 and A2667, Richard et al. 2010a for A370), we update them using the most recent spectroscopic information, in particular coming from MUSE (Jauzac et al. 2016;Lagattuta et al. 2017). We summarize the references of the lensing models used in Table 1.

SAMPLE CHARACTERISATION

We start by describing the lensed morphology of these galaxies and then proceed to derive their integrated properties, such as mass, gas-phase metallicity3 and star formation rate.

Morphology

Most galaxies studied in this work have multiple images, i.e., the same object can be seen in two (or more) different positions in the sky with different distortions due to lensing effects (see Fig. 1). Often, the most extended objects do not contain the entire image of the galaxy, that is, only part of the galaxy is lensed into multiple images, which typically appear mirrored with respect to the critical line. One advantage of these multiple images is that the integrated signal of these gravitational arcs is very high, since we observe the same galaxy twice or more, depending on the multiplicity of the image. However, a potential drawback is that properties measured in gravitational arcs, such as mass or morphology, may not be a correct representation of the full galaxy, since the galaxy images may be only partially lensed. To overcome this issue, less magnified but complete images, usually called counter-images, can be used to recover the full morphology of the galaxy and to rescale the physical properties measured in the high signal-to-noise ratio data from gravitational arcs.

A good example of all these effects is A370-sys1, dubbed "the dragon" (see Fig. 1, upper right column). The "head" of the dragon is the complete (or counter) image of this galaxy, where the entire galaxy is lensed with an average magnification factor of 9 and the distortion is small. The "body" of the dragon (the elongated image) is composed of two to four multiple images (Richard et al. 2010a), each of them containing only a part of the galaxy. The magnification in these multiple images is much higher than in the complete image, reaching factors of almost 30. It is in these higher magnification images that the smallest physical scales can be resolved, down to 100 pc for A370-sys1. Similar structures can be seen in A2667-sys1, M1206-sys1, M0416-sys28 and A521-sys1, while AS1063-arc and A2390-arc are single images.

M0416-sys28

M1149-sys1

Figure 1 -continued MNRAS 000, 1-28 (2018) We apply the lensing model to each set of multiple images to obtain reconstructed images, or source plane images images. From these source plane images, we then estimate the effective radius of these galaxies. We do this by fitting a 2D exponential profile to the reconstructed F160W HST image of the complete images. We leave as free parameters the centre of the galaxy, the effective radius, the inclination and position angle. The fit was done using the AstroPy 2D modelling package (Astropy Collaboration et al. 2013). We caution that these galaxies display complex morphologies, with visible spiral arms and bulges, that were not modelled in the fit. However, we masked the central bright region in A370-sys1, AS1063-arc, M1149-sys1 and A2390-arc, where a bulge is possibly present, to ensure it did not bias the fit. We list the fits parameters in Table 2 and plot ellipses with the best fit position angle, inclination and centre at one effective radius with magenta markings in the lower right hand panels of Fig. 1.

Spectra Extraction and Photometry

The rest of the analysis is performed in the image plane, working directly in data cubes and HST images. For the spectral characterisation of these galaxies, we are mainly interested in extracting high signal-to-noise ratio (S/N) spectra, where emission lines and continuum can be best analysed. Since the sources have an irregular shape, we do extract spectra from circular or elliptical apertures in the MUSE cube, but instead define extraction areas by imposing a flux threshold on MUSE pseudo narrow bands. To achieve this, we produce MUSE pseudo narrow band images centred on the [O ii ] doublet, where we choose a spectral width that maximizes the S/N in the narrow band. We measure the background flux level and variance in these images, and select pixels that are above a 3σ threshold, summing the spectra of the respective spaxels (see Fig. 1). For M0416 and M1149 we exclude multiple images heavily contaminated by cluster members. For M1206, we did not include the counterimage, since it does not significantly increase the signal-tonoise ratio.

When coadding the spectra, magnification effects are corrected on a spaxel by spaxel basis, scaling the flux in each spaxel by 1/µ. This is done in order to avoid differential magnification issues. Nonetheless, we note that applying a global (i.e. averaged) magnification factor leads to negligible differences on the overall sample. The average µ value for each spectrum is listed in Table 1. The variance propagated during data reduction is extracted in a similar way.

Photometry was derived using the publicly available HST data. Fluxes are measured over the multiple image which gives a complete coverage of the source in each system. Therefore, for A370-sys1, A2667-sys1, M1206-sys1 and M0416-sys28, the photometry was measured in the counter image and not in the arc. To minimise contamination from cluster members that are too closet to the gravitational arcs, we performed a 2D fit to the light of close cluster galaxies. We do this separately for each photometric band, by modelling each cluster member with a Sérsic profile and subtracting it from the image. The residuals of this subtraction were measured and included in the error budget of the photometry. Next, the F814W image, that matches the MUSE wavelength range, is used to define an extraction mask, by measuring the mean background and variance and placing a 3σ threshold. For A521, the F606W image was used, since there no F814W observations available. We also mask strong residuals arising from the cluster member subtraction. The flux inside the aperture is summed and the background of the image subtracted. The luminosity-weighted magnification factors over the corresponding HST apertures are calculated and the observed fluxes corrected. Uncertainties on the magnification factors are obtained from the Markov-Chain Monte-Carlo samples produced by Lenstool as part of the optimisation of the mass model. The photometry and magnification factors of the HST apertures are listed in Table B1 and the spectra are shown in Fig. A1, in Appendix.

The magnification-corrected photometry was used to renormalise the MUSE spectra, correcting both for aperture size and image multiplicity. We do this by measuring the flux in the extracted spectra (applying the F814W transmission curve), and normalising it to the value obtained from the HST F814W image (discussed below). This process yields high signal to noise spectra, corrected for lensing and aperture effects, that we use to determine the integrated properties of the sample.

Emission Line Measurements

Besides strong [O ii ]λ3726,29 present in all galaxies due to the selection criteria, six out of the eight galaxies display [Ne iii ]λ3869 emission. Other strong lines, such as [O iii ]λ4959,5007 and Balmer emission and absorption lines, as Hβ and Hγ, can also be found in the spectra depending on the object redshift. The narrow line profile of the emission lines, particularly at its centre (see in Section 4), as well as the absence of [Mg II] emission, make the presence of broad-line AGN unlikely, although we cannot completely rule out this possibility.

In order to better constrain the properties of the emission lines (peak, flux and line width) we first fit and subtract the spectral continuum. To do this, we use the pPXF code version 6.0.2 (Cappellari 2017) and a sample of stellar spectra from the Indo-US library (Valdes et al. 2004), which includes all stars with no wavelength gaps (448 stars). This library covers the entire MUSE observed wavelength range with a spectral resolution of 1.35 Å (FWHM) at rest frame. This corresponds to a FWHM of 2.7 Å in the ob-Table 3. Kinematic properties from integrated spectra. Redshift (z) was measured from the emission lines. Uncertainties were calculated using the Monte-Carlo technique described in subsection 3.3 and the error on the last digit is indicated. The uncertainties do not include systemic errors on the wavelength calibration (∼0.030 Å). The velocity dispersion (σ EL ) was measured from the FWHM of the emission lines, corrected for the instrument Line Spread Function (LSF), but not for beam smearing (which is done in Section 4). The velocity dispersion of the stellar population (σ ), measured with pPXF was also corrected for instrumental broadening, with a constant value of 2.5 Å. Values for M1149-sys1 are not reported since there are no strong stellar features in the spectrum (see Appendix C). Errors of σ are formal errors of the χ 2 minimisation. Both velocity dispersions are presented in rest-frame.

Object z σ EL σ (km s -1 ) (km s -1 ) AS1063-arc 0.611532[7] 69±5 108±6 A370-sys1 0.72505[1] 61±2 83±6 A2390-arc 0.91280[4] 88±10 78±1 M0416-sys28 0.939667[8] 48±2 374±39 A2667-sys1 1.034096[2] 51±1 27±10 M1206-sys1 1.036623[3] 98±3 69±2 A521-sys1 1.04356[2] 64±2 58±2 M1149-sys1 1.488971[1] 55±2 -
served frame at z = 1. This is comparable to the MUSE LSF, about ∼2.5 Å, with a wavelength dependency, as determined in Bacon et al. 2017. 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. We manually tune the degree of the additive polynomial by evaluating the quality of the fit as well as the distortion it causes in the shape of the original best fit (i.e., the combination of the stellar templates), since the polynomial can change the original continuum shape. Results of the fit are shown in Fig. A1. The continuum-subtracted spectra were used to measure emission lines fluxes. We use a Gaussian model implemented in the mpdaf package [START_REF] Bacon | MPDAF: MUSE Python Data Analysis Framework[END_REF]. To estimate the error of the fit, we fit 500 realisations of the spectra, drawn randomly from a Gaussian distribution the with mean and variance corresponding to the observed spectrum flux and variance. We then take the error of the measurements (flux, FWHM and z) as the half-distance between the 16th and 84th percentile of the 500 fits. The flux values are presented in Table B2, and the redshift and line width of the spectra are listed in Table 3.

Metallicity and Star Formation Rates

From the emission line fluxes, we simultaneously derive the global metallicity and attenuation factor of the galaxies. We do this by fitting multiple metallicity and attenuation sensitive line ratios. This method is a modification to that presented by Maiolino et al. 2008 (herein M08), whereby we extend their line ratio list and place it on a more formal Bayesian footing. We use the following line ratios: (Storey & Zeippen 2000) 7 Hγ/Hβ, 0.466 (Osterbrock & Ferland 2006, case B) 8 Hδ/Hβ, 0.256 (Osterbrock & Ferland 2006, case B) 9 Hδ/Hγ, 0.549 [START_REF] Osterbrock | Astrophysics of gaseous nebulae and active galactic nuclei[END_REF], case B) 10 H7/Hγ, 0.339 (Osterbrock & Ferland 2006, case B) The line ratios [1][2][3][4][5] are the same as those adopted by M08. These line ratios provide constraints on metallicity. But, because line ratios [1,3,5] compare lines at disparate wavelengths, they are also sensitive to the attenuation. So, to help alleviate the metallicity-dust degeneracy, we include ratios that are predominantly independent of metallicity [6][7][8][9][10].

1 [O ii ]λ 3727/Hβ, M08 2 [O iii ]λ 5007/Hβ, M08 3 [O iii ]λ 5007/[O ii ] λ 3727 , M08 4 [Ne iii ]/[O ii ]λ 3727, M08 5 R23 4 , M08 6 [O iii ]λ 5007/[O iii ]λ 4959, 2.98
We use all line ratios for which the lines are observed. However, since H β is not available for galaxies at z >0.9, we replace the [O ii ]/Hβ in ratio [1] with [O ii ]/Hγ (assuming case B Hγ/Hβ ratio of 0.466, for T e = 10000 K and low electron density). In Table 4, we list the specific line ratios used for each galaxy.

Another modification to the M08 method is that we adopt the Charlot & Fall 2000 dust model to correct the observed line fluxes (F obs ) for dust attenuation, assuming an index of -1.3 for the exponential attenuation law, proposed by the authors and found to be appropriate for star-forming regions (e.g. Brinchmann et al. 2013). The intrinsic fluxes are corrected in the following manner:

F int = F obs e τ v λ 5500 Å -1.3 (1)
We derive metallicity (Z) and attenuation (τ v ) from emission line ratios, fitting several ratios simultaneously in a Bayesian framework, using the emcee (Foreman-Mackey et al. 2013) Markov chain Monte Carlo Sampler to maximise the following a Gaussian (log-)likelihood function:

ln p = - 1 2 r (M r (Z) -O r (τ)) 2 σ 2 r + ln(2πσ 2 r ) (2)
where O r (τ) are the observed line ratios corrected for attenuation factor τ and M r (Z) are the predicted ratios using the above calibrations for the metallicity Z. σ 2 r is the quadratic sum of the observed error and an additional model uncertainty. We adopt a model uncertainty of 10% for the M08 calibrations, a 1% uncertainty for the Balmer line ratios and 2% for the [O iii ]λ 5007/[O iii ]λ 4959 ratio (Storey & Zeippen 2000). We assume a solar abundance of 8.69, as in M08. We use wide flat priors for all galaxies with metallicity 7 < 12 + log(Z) + 12 < 10 and attenuation 0 < τ v < 4. The marginalised metallicity and attenuation distributions for the AS1063-arc can be found in Fig. 2 and the remaining distributions in Appendix D. The median of these distributions are listed in Table 4, with the half-distance between the 16th and 84th percentiles adopted as the error.

M1149-sys1 was left out of this analysis, since only [O ii ] is present in the MUSE spectrum. All the galaxies analysed have well determined metallicities and attenuations, ranging from 8.71 to 9.05 and 0.09 to 2.86 respectively, with symmetrical marginalised distributions (except for attenuation in M0416-sys28, A2390-arc and A521-sys1). For A2390, the only metallicity diagnostic that could be used was [O ii ]λ 3727/Hγ, since the signal to noise ratio does not allow a robust measurement of the [Ne iii ] flux. As the [O ii ]λ 3727/Hγ diagnostic is degenerate, with both low and high metallicity solution, we investigate each branch by repeating the fit allowing only low or only high metallicity values ([7.0,8.6] and [8.6,10.0] priors, respectively). For solutions found in each branch, the [Ne iii ] flux would be ∼6.2 ×10 -18 erg cm -2 s -1 and ∼1.2 ×10 -18 erg cm -2 s -1 for the low and high metallicity solutions respectively. The low metallicity estimate approximately corresponds, for example, to the flux of Hδ seen in this galaxy (∼ 5.7 ×10 -18 erg cm -2 s -1 ), so the low metallicity branch seems less likely than the high metallicity solution. The high metallicity solution would yield lower [Ne III] fluxes that would likely be undetected given the ∼0.3 ×10 -18 erg cm -2 s -1 1σ spectral noise for this galaxy. Furthermore, the star formation rates derived via the Balmer lines and the [O ii ] line (corrected for metallicity) are closer to the high metallicity case than the low one. These indications do not rule out the low metallicity hypothesis, but we adopt the high metallicity values as the most probable solution for A2390-arc throughout the rest of this work.

The star formation rates (SFR) are calculated both from Balmer lines and the [O ii ] doublet. For the Balmer lines, the Kennicutt 1998a calibration is adopted, adapting this calibration to be used with Hβ or Hγ and assuming Table 4. Metallicity, attenuation and intrinsic SFR. Gas-phase metallicity and attenuation (τ v ) were calculated as described in subsection 3.4 and using the diagnostics listed in the 'Line Ratios' column. Both star formation rates derived from the Kennicutt 1998a (SFR Balmer ) and Kewley et al. 2004 4.

Stellar Mass

We fit both spectra and photometry using the FSPS (Flexible Stellar Population Synthesis, Conroy et al. 2009) models implemented within a Bayesian framework in the Prospector5 code (Johnson et al., in prep). The only exception is A521-sys1, for which only the spectrum was fit and used to normalise the spectra, since only one HST filter was available. Prospector allows us to fit models produced with FSPS 'on the fly', i.e., the fit are not based on a precomputed grid of models but are instead computed for each exact set of parameters tested. We assume a decaying star formation history and choose the Padova isochrones (Marigo & Girardi 2007;Marigo et al. 2008) with the Chabrier 2003 initial mass function and the MILES stellar library (Sánchez-Blázquez et al. 2007;Falcón-Barroso et al. 2011) to generate the models for all fits. We also adopt the Charlot & Fall 2000 dust attenuation here, with a fixed index of -1.3 for the star-forming regions, applied only to the star-light originating in stars younger than 10 Myrs, and an index of -0.7

for the global dust screen as suggested by Charlot & Fall 2000, applied to all starlight equally. In total, 6 physical parameters were fit: stellar mass, the e-folding time of the star formation history, stellar population age, stellar metallicity and optical depth of both birth clouds and interstellar medium (ISM) (see Table 5 for details and priors). Beside these physical parameters, a 3rd order polynomial is also fitted (similarly to what is done with pPXF). For the sake of simplicity on an already quite complex fit, we decided to mask the emission lines during these fits, and do not attempt to fit them during the continuum fit with PROSPECTOR and do not include nebular emission in the models.

The fits for each galaxy can be found in Appendix C. The best values and errors, estimated from the marginalised distributions as in the previous section, can be found in Table 5. The magnification errors were included at this point in the mass and SFR error estimates. Overall, the best fits provide a good match to both the spectral and photometric data. Most of the derived model parameters are well constrained, with the exception of the attenuation, which in some instances have large tails extending to high values. In principle, the value of τ BC obtained here should agree with the attenuation derived from the emission lines. In practice, since no FUV data is available, nor were the emission lines fitted with prospector, this value cannot be constrained and will be degenerate with τ I S M . The obtained stellar masses vary from ∼1 × 10 10 to ∼9 × 10 10 M , with a mean value of 3.65 × 10 10 M , and with an e-folding time that varies from 3 to 9 Gyrs and dominant ages from 3 to 7 Gyrs.

Comparison with other samples

We can now use the physical parameters previously derived both from emission lines and photometry to place this sample in context. For M1149-sys1, only the [O ii ] line is present in the MUSE wavelength range, so it is not possible to constrain metallicity and SFR from these data. Two metallicity estimates exist for this galaxy in the literature. Yuan et al. 2011 Firstly, we compare the sample with the mass-SFR relation (top panel of Fig. 3). Within errors, all the galaxies lie close to what is measured in main sequence star-forming galaxies. The only exception is A370-sys1 that has a slightly lower star formation rate than expected from its derived mass. We then compare this sample with the fundamental metallicity relation (Lara-López et al. 2010;Mannucci et al. 2010), that defines a tight relation between stellar mass, gas metallicity and SFR of normal star-forming galaxies (bottom panel of Fig. 3). Five out of the seven arcs -AS1063arc, A370-sys1, A2390-arc, M0416-sys28 and M1206-sys1differ by up to 0.1 dex from this fundamental relation, which is higher than the dispersion derived for the local SDSS sam- ple of about 0.05 dex, but still within the 0.5 dex dispersion observed for higher redshift samples (Mannucci et al. 2010).

We conclude that all the galaxies studied in this work are normal star-forming galaxies and are good representatives of the population of disc galaxies at z ∼ 1.

We can estimate the gas mass of these galaxies by using the Kennicutt-Schmidt law (Kennicutt 1998b) to convert SFR surface density to gas surface density. We take the Table 5. SED fitting best results and their error (16th and 84th percentiles) obtained with Prospector and including magnification corrections. M corresponds to the current mass of the galaxy; τ is the e-folding time of the star formation history and age the age of the composite stellar population. We estimate the average SFR (SFR SED ) as the ratio of the formed mass (higher than the current mass, which does not include gas that was recycled into the ISM) and the age of the galaxy (the age parameter). τ V , BC is the attenuation factor for the birth clouds, applied only to starlight coming from stars with less than 10 Myrs, and τ V , I S M is the global attenuation factor of the ISM and is applied to all starlight, respectively. The prior ranges of each parameter are shown in the first row. With the exception of τ, where a logarithmic prior was used, all other priors are flat. The maximum allowed age is set by the age of the Universe at the redshift of each galaxy.

Object M τ age SFR SED Z τ V , BC τ V , I S M (10 10 M ) (Gyr) (Gyr) (M yr -1 ) (log Z/Z ) Range [0.5 , 10] [0.01 , 10] [0.2 ,4-8] - [-1, 1] [0, 4] [0, 4]
AS1063-arc 8.74 global (dust and magnification corrected) SFR derived from the Balmer lines of the integrated spectrum and we calculate the SFR density by dividing it by the photometry aperture area (the area to which the spectra were normalised). Note that the area is calculated in the source plane. The SFR is then converted to gas surface density, and thereby also the total gas mass. These are listed in Table 6. The star formation rate densities, gas masses (M g, KS ) and gas fractions (M g, KS /(M g, KS + M )) are listed in Table 6. These galaxies have low gas fractions that range from ∼0.1 to 0.3, despite the star-forming clumps clearly visible in most of them (particularly AS1063-arc, A370-sys1 and M1206-sys1). These values are consistent with the values obtained from the DYNAMO survey, a sample of comparable local galaxies that covers a stellar masses of 10 9 to 10 11 M and star formation rates of 0.2 -100 M yr -1 , selected to be analogues to z ∼3 galaxies (Green et al. 2014). The gas fraction in our galaxies are also similar to what White et al. 2017 found from making CO measurements and inverting the KS law.

RESOLVED PROPERTIES

Observed Velocity Maps

The velocity maps (or kinematic maps) were derived measuring the [O ii ] emission lines present in all galaxies. The fits were performed with a slightly modified version of the camel6 code (Epinat et al. 2010), where we added the option of fitting binned data. camel fits emission lines with 1D Gaussian models, having as free parameters the redshift, the Gaussian FWHM and the flux of each emission line. In order to more robustly probe the outer parts of the galaxies, we bin the data using the Voronoi binning method7 presented in Cappellari & Copin 2003. We bin on the signal to noise of the [O ii ] pseudo narrow bands used in Sect. 3.2. We impose a signal to noise of 5 per bin in all galaxies, which we find to be a good compromise between high spatial sampling and robust spectral fits. The velocity dispersion was allowed to vary between 0 and 250 km s -1 . The 2D observed velocity maps of the 8 galaxies are presented in Fig. 4. These maps were inspected and the velocities values were corrected (adding or subtracting a constant value in the entire field) so that the velocity at the morphological centre (as given by the 2D exponential fit to the source plane image, see Table 2) corresponds to the kinematic centre and has velocity zero.

Fitting Velocity Maps

We fit the 2D velocity maps with the three following kinematic models: an arctangent model (Courteau 1997), an isothermal sphere model (Spano et al. 2008) and an exponential disc model (Freeman 1970). The first is an empirical model commonly used to fit local galaxies. The isothermal sphere model is a good approximation of the kinematics of a dark matter dominated galaxy, under the assumption that the dark matter halo has an isothermal profile. The exponential disc assumes a smooth exponentially decaying distribution of baryonic matter (and no dark matter) as seen in most disc galaxies (see Epinat et al. 2010 for a summary of these models and analytical expression implemented in this work).

We fit these models directly to the 2D velocity maps produced by camel (i.e. in image plane). Note that we do not fit the velocity dispersion maps. To account for lensing distortions, we produce high resolution displacement maps using the optimised cluster mass models and the Lenstool software. These maps describe the geometrical transformation from source to image plane by predicting the spatial shifts (or displacement) that have to be applied to each source plane pixel to be projected to the image plane image.

The first step of the fit is to create a 2D kinematic map in source (unlensed) plane. We apply the lensing distortion to this kinematic model, using the displacement maps. We then convolve the lensed model velocity map with the seeing, using a Moffat kernel, flux-weighted by the [O ii ] flux image from camel. Finally, we apply the same binning to the models as the observations before comparing the modelled and the observed velocity fields. It is this final product that is compared with the data.

Each model has seven free parameters, that are explored using a MCMC sampler:

• X,Y position of the kinematic centre in image plane.

• V sys systemic velocity.

• Inclination of the galaxy (inc).

• Major axis kinematic angle in source plane (PA).

• Maximum Velocity (V max ).

• Transition radius (r t ).

The first three parameters (X,Y,V sys ) are mainly introduced to correct for small inaccuracies in the lensing correction (V sys is used to keep the v=0 km s -1 at the centre even if the centre changes slightly.). We assume that the kinematic centre is the same as the morphology centre (measured in the reconstructed image), and correct the observed velocity fields in order to have zero velocity at this position. To correct any possible errors when passing from the position in source plane to image plane, during the fit we allow the kinematic centre to vary within a 2×2 MUSE pixel box around the original position and let the velocity adjust up to ±100 km s -1 . The corrections for all galaxies were small with displacements of 1 pixel and velocity adjustments of up to ±30 km s -1 .

The maximum velocity and transition radius (V max , r t ) correspond to the maximum velocity (and the radius at which it occurs) for the exponential disc and isothermal sphere models. The arctangent model velocity does not reach a maximum at a finite radius, so the transition radius in the arctangent model corresponds to the radius at which 70% of the asymptotic velocity (V max ) value is reached, and it is therefore generally lower for the arctangent model than for the other two kinematic models.

Table 7. 2D kinematic fit priors and results for the three kinematic models tested: arctangent ("Atan"), exponential disc ("Exp") and isothermal sphere ("Iso"). Fitted parameters are the inclination ("Inc"), the position angle (PA), the maximum velocity (V ma x ) and the transition radius (r t ). The centre of the models was aligned with the morphological centre measured in Sect 2.2. For AS1063, the inclination was fixed, since the seeing of the observations (∼1") did not allow to constrain this quantity. The position angle is given in source plane. The logarithmic evidence for each model is listed in the last column and the lowest evidence model is highlighted in bold in the first column. We have rejected and masked (in grey) values with unreliable beam-smearing corrections (σ smear >σ obs ).

Inc

MNRAS 000, 1-28 ( 2018)

Figure 4 -continued MNRAS 000, 1-28 (2018)

MUSE Gravitational Arcs 15

We use the emcee (Foreman-Mackey et al. 2013) MCMC (Monte Carlo Markov Chain) implementation, with 100 walkers and let the chain evolve for 1000 steps, after an initial burn-in phase of 100 steps. The chains are inspected to confirm convergence. Also here, we adopt a Gaussian likelihood function:

ln p = - 1 2 r (Model(inc, P A, V max , r t ) -V el) 2 σ 2 + ln(2πσ 2 ) (3)
where Model(inc, P A, V max , r t ) is one of the kinematic models (arctangent, isothermal sphere or exponential disc) with inclination inc, position angle P A, maximum velocity V max and transition radius r t , after being lensed, convolved with the seeing and binned. The model is compared with the observed velocity field (V el) and weighed by the errors, σ. For these, we only consider the errors of the measured velocity, as given by camel, corresponding to the formal errors of the χ 2 minimisation of each bin. We assume flat priors to all the parameters being fit.

We choose a wide prior for the maximum velocity and transition radius that is the same for all galaxies. The inclination was given a wide prior of ±20 degrees around the inclination derived from the morphological fit. The only exception is AS1063-arc, observed with the worst seeing in this sample (∼1"), where the inclination was unconstrained by the kinematic data and was set to the value measured in the HST image. The position angle was also given a wide prior (±20 degrees) around the position angle measured in the exponential fit to the HST images. Finally, we present the Bayesian evidence of each model, a relative measure of how well a model performs, as calculated by emcee.

Results

The results for each model are presented in Table 7, where the best value of each parameter corresponds to the median of the marginalised distributions. Also here, the half distance between the 16th and 84th percentiles was taken as the error for each parameter. From the Bayesian evidences values, the preferred kinematic model in this sample is the arctangent model, with galaxies AS1063-arc, A2667-sys1, M0416-sys28 and M1206-sys1 being better fit by this model. (However, notice that the fit of M0416-sys28 is not constrained, with the transition radius hitting the upper limit allowed). The remaining galaxies, A370-sys1, A521-sys1 and M1149-sys1 are better fit by an exponential disc model and only A2390arc by an isothermal sphere. And although there is a clearly preferred model for each galaxy, none of the models is globally preferred. The position angles (PA), that are expected to be the same among the three models, agree quite well for each galaxy. The inclination shows a larger spread, but still displays a reasonable agreement. The 2D maps for the best model with the parameters of Table 7 are plotted in the middle panel of Fig. 4.

As previously stated, intrinsic velocity dispersions are difficult to estimate due to the beam-smearing effect. Here we benefit from high magnification of the lensing and the good seeing of the MUSE observations that reduce this problem, but nevertheless we still correct each pixel using the 2D Table 8. Kinematic properties of the sample. V max is the maximum velocity of the best model (repeated from Table 7); σ is the mean and standard deviation of intrinsic velocity dispersion (corrected for instrumental broadening and beam-smearing) measured in the 2D maps of Fig. 1. We also calculate the Gini coefficient of the intrinsic velocity dispersion maps (G σ ) and of the [O ii ] narrow band images (G [O ii ] ), as a comparison. The velocity dispersions are high but quite homogeneous. by: σ = σ 2 obs -σ 2 inst -σ 2 smear . We plot these 2D maps in the lower right panel in Fig. 4.

Object V max σ V max /σ G σ G O ii ( km s -1 ) ( km s -1 ) AS1063-arc 296±2 84±33
All galaxies are rotation-dominated, with regular 2D velocities and V/σ ratios between ∼2 and ∼10, estimated from the maximum velocity V max and the mean of the σ maps (see Table 8). These values are in agreement with V/σ ratios found for other high-z samples such as the KMOS 3D kinematics survey (Wisnioski et al. 2015), with V/σ between ∼3-11 for galaxies at z ∼ 1, the MUSE sample of low mass galaxies (Contini et al. 2016), with the range ∼1-5 and the KLASS survey of lensed galaxies at 0.7 < z < 2.3 (Mason et al. 2017), with V/σ ∼8 at z ∼ 1.

DISCUSSION

In this section, we discuss the implications of the kinematic fit results and the morphology of the velocity dispersion maps. We also explore the connection between star formation and turbulence in these typical z ∼1 discs. A summary of the individual properties of the galaxies, and comparison with previous studies of these objects, can be found in Appendix E.

Best Kinematic Model

Recent work by Genzel et al. 2017 andLang et al. 2017 found that both the individual velocity curves of massive galaxies and the stacked velocity curve of 101 star-forming galaxies in the range 0.6 < z < 2.6, display a significant drop in the velocity curve beyond the turnover radius. This would (Mason et al. 2017) are plotted as black diamonds (we have excluded possible mergers and galaxies with only upper limits). The KMOS Redshift One Spectroscopic Survey Sample at z∼1 from Johnson et al. 2018 is also plotted in grey dots. The gravitational arcs sample velocity dispersions are in good agreement with all these samples. imply a higher baryonic mass fraction in high-z galaxies than previously though. We briefly investigate this issue.

As stated before, four out of eight galaxies are better fit by the arctangent model, an empirical model. The exponential disc model, which describes rotation velocity of a gas disc, has lower Bayesian evidences on three of the galaxies and the isothermal sphere model, an approximation to what is expected for an isothermal dark matter halo, on only one. Most importantly, all models are remarkably similar in the inner parts of the galaxies. For most galaxies, the three models (with the best parameter estimations) are visually indistinguishable up to 1 R e (e.g. A370-sys1 and A2667-sys1), or even 2 R e (e.g. AS1063-arc).

In none of the galaxies analysed here do we see a drop in the velocity curve at large radii, which is confirmed by the fact that the arctangent model (the only one that does not decrease velocity with radius) is the most preferred model. However, with the exception of AS1063-arc and A2667-sys1, our data does not extend out as far as the stacked sample of Lang et al. 2017, so we cannot rule out the possibility of a drop in the velocity beyond 2 effective radii. We also cation the reader that our analysis is too simplistic to explore this subject, since a combination of the two componentsbaryons and dark matter -is expected, and a joint fit should be done in order to gain some more insight on this matter, as done in Lang et al. 2017.

Velocity Dispersion Morphology

We now turn our attention to the velocity dispersion. The mean intrinsic velocity dispersion of these disc galaxies ranges from 15 to 84 km s -1 , with a mean value of 48 km s -1 . We plot these values in Fig. 5 and compare them with other samples. The mean velocity dispersions are comparable to the global values derived for the KMOS 3D , the MUSE samples for 0.2 < z < 1.4 galaxies, and the lensed sample from the KLASS survey (Wisnioski et al. 2015, Contini et al. 2016and Mason et al. 2017, respectively). Our velocity dispersions are also in agreement with the values derived for the DYNAMO sample of local analogues Green et al. 2014;Oliva-Altamirano et al. 2017. When compared to the GHASP survey (Epinat et al. 2010), a sample of nearby latetype galaxies, these lensed galaxies have velocity dispersions higher by a factor of 1.4 to 4.2.

Several studies have investigated the kinematics of lensed galaxies (e.g Jones et al. 2010;Livermore et al. 2015;Leethochawalit et al. 2016;Mason et al. 2017), mostly focusing on the velocity fields of these galaxies and not so much in the resolved velocity dispersions. The sample presented in the current work was selected to be extremely extended and with high physical resolution and the high signal to noise of our data results in relatively small formal errors per bin in the 2D velocity dispersion maps. These are derived with camel and, for AS1063-arc, A370-sys1, A2390-arc, M0416-sys28, A2667-sys1 and M1206-sys1, have a median value that varies between 5 and 8 km s -1 (but can be as high as 50 km s -1 ). The remaining galaxies, A521-sys1 and M1149-sys1, have higher median errors, 24 and 16 km s -1 , respectively. Another source of uncertainty is the LSF correction, that we perform using expression 7 of Bacon et al. 2017. For our lowest redshift galaxy, the LSF dispersion is about 1.1 Å (see figure 15

M1149-sys1

Figure 6. Radial velocity and velocity dispersion in source plane for individual bins. Velocity was extracted along the major kinematic axis with an aperture of 0.5" both for the data (back circles) and for the models (thick coloured lines) and is inclination corrected. Velocity dispersion corresponds to the intrinsic velocity dispersion (beam-smearing corrected) and was extracted in the full galaxy. A vertical bar marks the effective radius (R e ).

a dispersion of ∼50 km s -1 . For the higher redshift galaxies this value is lower. Neither of these sources of uncertainty directly depends on the position of the bin in the velocity map. On the other hand, the beam-smearing correction does depend on the position of the bin in the 2D map. Calculate the errors associated with this correction is computationally expensive (since a reasonable number of velocity dispersion maps would have to be produced), so we estimate the uncertainties associated with the beam-smearing correction in a simplified way. We produce 16 velocity dispersion models with combinations of perturbed V m ax, r t , P A and the inclination. Since this is a simplified estimation of the error, we perturbed these parameters by ±20 km s -1 , ±1 kpc and ±5 degrees, values that are considerably higher than our typical errors. Comparing these perturbed models with the best model for each galaxy, we obtain median differences between 5 and 10 km s -1 . Overall, our major source of error is the LSF estimation. Nevertheless, Bacon et al. 2017 measure the LSF to be largely constant in the MUSE FOV, with an average standard deviation of ∼ 0.05 Å, therefore the effects on the study of the velocity dispersion structure are expected to be small.

The observed, σ smear and intrinsic velocity dispersion 2D maps of the 8 galaxies analysed in this work can be seen in the lower panels of Fig. 4. Both the observed and the σ smear (the dispersion expected due to observational effects) maps have a higher velocity dispersion at the centre of the galaxies (exceptions are A2667-sys1 and M0416-sys28). The intrinsic 2D velocity dispersion maps reveal no significant 2D structure at our precision limit, with the dispersion of the σ maps varying between 7 and 33 km s -1 (see Fig. 4). Our intrinsic velocity dispersion maps are broadly consis-tent with a constant velocity dispersion value throughout the entire galaxy.

In order to better quantify how homogeneous the velocity dispersion maps are we also measure their Gini coefficient [START_REF] Gini | Variabilità e mutabilità: contributo allo studio delle distribuzioni e delle relazioni statistiche[END_REF]. This is a measure of inequality in a data set and can vary from 1, absolute inequality, to 0, absolute equality (where, in this case, all the pixels have the same value). We obtain values between 0.14 and 0.30 with a mean value of 0.23 (see Table 8). This suggests a fairly homogeneous distribution of the values of the velocity dispersion and it is significantly lower than the values obtained doing the same exercise with the [O ii ] pseudo-narrow bands, that map the starforming clumps, where we measure Gini coefficients between 0.82 and 0.90. The relatively homogeneous velocity dispersion maps, compared with clumpy [O ii ] and UV morphologies, suggest that the star-forming clumps do not strongly disturb the local velocity dispersion, i.e., the clumps do not seem to locally increase turbulence. We caution however that we did this exercise based on individual pixel values, not deriving velocity dispersions per clump. Recent work analysing the DYNAMO sample by Oliva-Altamirano et al. 2017 found statistically significant structure in 6 of these local analogue galaxies with increased velocity dispersion possibly associated with clumps. These structures generally show an increase of the order of ∼20 km s -1 . This is very close to the precision we can reach with our data, and it is unlikely that we would detect this, even if this structure were present.

Turbulence and star formation

As stated in the introduction, the reasons why high redshift discs are more clumpy and turbulent, and how these two properties relate to each other, are still under discus- sion. While some works point to the higher gas fractions of these younger discs as the reason for their clumpy morphologies and turbulence (e.g. Genzel et al. 2011;Livermore et al. 2015), others found a correlation between SFR and higher turbulence (e.g. Green et al. 2010;Swinbank et al. 2009). We investigate this in the two lowest redshift galaxies of this sample, where the resolved extinction and extinction corrected SFR maps can be robustly derived. We will present and discuss the metallicity maps of these galaxies in a separate paper (Patrício et al., in prep.) and focus here on the SFR density.

To produce the star formation density (Σ SFR ) dust corrected maps, we first subtract the continuum from the observed cube, since this could bias the measurement of the Balmer lines. This is done by fitting the local continuum with pPXF. The data were spatially binned using Hγ pseudonarrow bands, which guarantees that both Hγ and Hβ can be reliably measured in each bin. The line fluxes are fit in the continuum subtracted cube using camel, as done for the observed velocity maps. We calculate the extinction and SFR from the emission lines in each bin as done in the integrated spectrum (see Section 3.4), with the exception that only the

[O iii ]λ 5007/Hβ, [O iii ]λ 5007/[O ii ]λ 3727, [O iii ]λ 5007/[O iii ]
λ 4959 and Hγ/Hβ ratios were used and that the parameter space explored this time was smaller (ex-tinction between 0 and 2.5) based on the results obtained in the integrated spectra. The resulting maps can be seen in Fig. 7. Since velocity dispersion and SFR density were derived with a different binning, in the following analysis we only consider the binning in Σ SF R and obtain the corresponding velocity dispersion by binning the intrinsic velocity dispersion maps to the same grid.

There are different results in the literature on the relation between star-formation rate density and velocity dispersion. While Genzel et al. 2011 found little correlation between the two in their sample of massive giant clumps at z ∼2, Swinbank et al. 2012b, using a sample of galaxies observed with Adaptive Optics, found it was well described by the relation σ ∝ Σ 1/n SFR + k. Since clump selection depends heavily on the criteria and data used (e.g. Dessauges-Zavadsky et al. 2017) we use the full 2D information of these two galaxies and plot in Fig. 8, for each Voronoi bin, the star formation density and the respective velocity dispersion. We also gather data on the velocity dispersion and star formation density of the high-z clumps available in the literature: 5 giant star forming clumps at z ∼ 2 from Genzel et al. 2011, from the z ∼ 1 lensed galaxies clumps from Livermore et al. 2015 and from the clumps of nine 0.8 < z < 2.23 galaxies from the HiZELS survey (Swinbank et al. 2012b). for the Table 9. Relation between σ 0 and Σ SFR for different subsets of data in Fig. 8: Only Clumps for individual clumps from Genzel et al. 2011, Swinbank et al. 2012band Livermore et al. 2015; AS1063-arc and A370-sys1 for the bin-per-bin values derived here and All Data for all the above. We list the Spearman rank correlation (ρ) and the and γ values derived from the fit and the reduced χ 2 of the fit. We also report the fit values from Genzel et al. 2011. Overall, σ last two, the Σ SFR were computed from the SFR and the radii of the clumps. We first investigate if these two quantities, Σ SFR and σ are correlated by computing the Spearman rank correlation (ρ) between the two (see Table 9). Since the data on the star-forming clumps gathered in the literature can be quite different from the individual bins studied here, we do this analysis separately.

The full data set on properties of clumps from Genzel et al. 2011, Swinbank et al. 2012band Livermore et al. 2015 shows very little correlation between Σ SFR and σ , with ρ=-0.12. Taking only the data of AS1063-arc we find a Spearman correlation of -0.21 and for A370-sys1 a correlation of 0.20, once more both very weak.

We proceed to fit both these quantities with a physically motivated relation, following Lehnert et al. 2009 andSwinbank et al. 2012b analysis. Assuming a marginally stable disc, with a Toomre parameter8 equal to unity, and the Kennicutt-Schmidt relation, we would expect that the velocity dispersion and the star formation density are related by:

σ 0 = Σ γ SFR (4)
Also here we analyse the data separately, fitting data only from individual star-forming clumps from the literature, only from AS1063-arc or A370-sys1, and finally all of it. The fit results and their reduced χ 2 are listed in Table 9. We generally find very low γ values, all below 0.1 except A370-sys1 (0.307±0.041), far from the ∼0.74 found by Swinbank et al. 2012b and implied by the local Kennicutt-Schmidt law. Lehnert et al. 2009 also present pixel-per-pixel data on 11 z ∼2 star-forming galaxies observed with AO and find that γ=0.33 is a good match to their data, which is in agreement with what we see for A370-sys1, but not to AS1063-arc or the individual clumps.

Because we made use of fully resolved (as opposed to clump) data of highly magnified objects, we reach lower star formation rate surface densities. Indeed, there is only a small overlap in Σ SFR between AS1063-arc and A370-sys1 when compared with the clump sample gathered from the literature. Other resolved studies, such as Lehnert et al. 2009 andSwinbank et al. 2009, generally reach star formation rate surface densities as low as 0.1 M yr -1 kpc -2 , while much of the our data sits below this limit (although see Cresci et al. 2009 for other Σ SFR < 0.1 M yr -1 kpc -2 measurements of high-z galaxies). However, for Σ SFR < 0.1 M yr -1 kpc -2 , velocity dispersions are expected to almost not vary with star formation surface density, even in a stellar-feedback driven turbulence scenario (Krumholz et al. 2017), so it is not surprising not to find a correlation at low star-formation rates. We can nevertheless look at the expected velocity dispersions for different scenarios for the measured stellar formation rate densities. For Σ SFR between 0.02 and 1 M yr -1 kpc -2 , where our data is more reliable, we find values of σ of 20 to 100 km s -1 . Comparing with theoretical predictions from Krumholz et al. 2017 (their figure 4) we see that these high velocity dispersions for low star formation density are only reached for models that include transport (i.e. gravitational instability) driven turbulence. Our data best compares with the transport+feedback and no-feedback predictions for local dwarf galaxies. These predictions are however done for integrated properties of the galaxies, and not resolved as presented in this work, which makes a direct comparison difficult.

We try to extend this study beyond these two galaxies, where we can correct each bin for extinction, to the other galaxies of the sample, but in a more simplistic way. Since we cannot derive resolved extinction maps for the remaining 6 objects, we test only if there is a correlation between the observed [O ii ] surface brightness (that probes star formation and is unaffected by magnification since gravitational lensing conserves surface brightness) and the velocity dispersion. A similar analysis, using Hα surface brightnesses, was presented in Lehnert et al. 2009 that found a positive correlation between the two. We plot these two quantities and the Spearman coefficients in Fig. 9. To make sure that the Spearman coefficients were not biased by unreliable measurements we have removed data points for which i) formal errors were 0 km s -1 ; ii) in which the velocity dispersion reached the maximum allowed value of 250 km s -1 ; or iii) that have S/N lower than 2. These unreliable values are plotted in black squares in Fig. 9.

This crude approach has limitations. Comparing the ρ coefficients obtained for AS1063-arc and A370 using both methods, we see they they are extremely different. However, in neither analysis do we find a strong correlation between star-formation rates (or [O ii ] surface brightness) and turbulence. All the 7 considered galaxies have ρ ≤ 0.26, denoting a very faint dependence between these two quantities (M0416-sys28 does not have enough data to perform this exercise). Bearing in mind the limitations of this test, we conclude that there is no strong evidence of local velocity dispersions being sustained by (low) star formation, which is corroborated by the more careful analysis of AS1063-arc and A370-sys1, and is in agreement with previous results from Genzel et al. 2011 andLivermore et al. 2015, but in contrast with what was found by Lehnert et al. 2009, Swinbank et al. 2012band Lehnert et al. 2013.

Energy Budget

In the previous section we argued that in this sample we do not see strong evidence that local SFR and velocity dispersion are correlated. There might be natural explanations for this: the sample probes low star formation rates, where it is not possible to distinguish if the turbulence is mostly driven by star formation. Another possible test to investigate the origin of the higher turbulence of high-z discs is to study the energy necessary to sustain those high velocity dispersions and where it could come from. Here we use (Krumholz et al. 2017) equations 1 and 5 to derive the energy associated with turbulence in a disc and star formation, respectively. We reproduce both expressions here:

dE dA turb = 3.1 × 10 9 Σ g [10 M pc -2 ] σ [10 km s -1 ] 2 (5) dE dA SF = 3.1 × 10 9 Σ g [10 M pc -2 ] σ [10 km s -1 ] r [10 kpc] V [200 km s -1 ] (6)
where energy associated with star formation is calculated over a galactic dynamical time t dyn = r V -1 ). We calculate both these energies within the photometric aperture (to which the total spectra were normalised), making use of the gas densities from Table 6. The velocity in expression 6 was calculated for the photometric aperture radius using the best model and parameters for each galaxy. We warn that Σ g and SFR are not independent, since Σ g was derived using the inverse Kennicutt-Schmidt law. We report these values in Table 10.

We find that in our sample, star formation can provide the necessary energy to sustain turbulence, requiring only 6% to 46% of the available energy. This is somehow at odds with what is found by Krumholz et al. 2017 in their disc evolution model where at z > 0.5 disc turbulence is dominated by transport-driven mechanisms (e.g. from gas accretion) and not stellar-feedback, especially for more massive galaxies, although both processes co-exist. Such a difference can come from the approximations made to calculate these energies, particularly in estimating the gas surface density.

SUMMARY AND CONCLUSIONS

We discussed the physical properties of seven clumpy discs at 0.6 < z < 1.5, that are highly magnified by strong lensing and can therefore be resolved to sub-kiloparsec scales. From high-quality MUSE data (with <0.7" seeing) and abundant HST archival data we derive the following:

• We derive stellar masses, global metallicities and extinctions from multi-band HST photometry and multiple line diagnostics for the full sample and for the first time for 3 of the 8 objects (AS1063-arc, A370-sys1 and A521-sys1). Comparing these objects with the main-sequence of star-forming galaxies, we confirm that these are typical starforming galaxies for the considered redshifts (Fig. 3).

• The observed velocity fields, derived from the MUSE data cubes, all reveal smooth velocity gradients with a velocity dispersion peak at the centre of rotation, typical of rotating discs. These are rotation-dominated galaxies with V/σ between ∼2 and ∼10, also typical of what was found by IFU surveys (Fig. 5). We do not see any evidence of a drop in the velocity curve within the sample, even for galaxies where 2 effective radii are reached (Fig. 6).

• The observed velocity maps were fit with three different kinematic models (arctangent, isothermal sphere and exponential disc model), with a novel method that allows to fit several multiple images at the same time, accounting for lensing distortions. The arctangent model was a better fit to 4 galaxies of the sample, the exponential disc model to 3 and the isothermal sphere model to the remaining galaxy. However, all models deliver very similar results, particularly in the central part of the galaxies, which seems to indicate that the velocity is not probed far enough to reliably infer which model is the best (and therefore to verify if the velocity curves decay at larger radii).

• The intrinsic velocity dispersion maps (that we correct for beam-smearing effect) do not display any strong 2D structure, revealing fairly isotropic 2D maps with a mean Gini coefficient of 0.23 (compared with 0.87 for the [O ii ] maps, Table 8).

• There is no strong evidence of increased (or decreased) velocity dispersion at the positions of the star-forming clumps. We also do not find a strong correlation between the local (resolved) SFR density and the velocity dispersion for the two galaxies where the dust attenuation corrected SFR maps could be derived in a resolved way (Fig. 8). Comparing the resolved results with theoretical predictions, we conclude that the high velocity dispersion found for relatively low star formation rate surface densities are more in agreement with scenarios where the turbulence is driven by gravitational instability (Krumholz et al. 2017).

• However, from an energetic point of view, the observed star-formation is sufficient to explain the observed turbulence.

• We find that the star-forming clumps share the dynamical properties of the ionised gas surrounding them (rotation velocity and turbulence). In a future work, we will study the chemical composition of these galaxies and the clumps within them, to further test the in-situ clump formation hypothesis.

Finally, here we only focus on the ionised gas, traced by emission lines. The analysis of the properties of molecular gas, made from high-resolution ALMA data obtained for a number of these objects, will allow to probe the turbulence of the cold gas and to have a more complete view of the process of formation of star-forming clumps.

M0416-sys28

This galaxy has the lowest mass on this sample 9.5±0.25×10 9 M and it is also the most compact, making it difficult to fit of the gas kinematics. For example, the transition radius (r t ) is not correctly fit for any model, reaching limit values of 12-14 kpc (we impose a limit of 15). We estimate a maximum velocity of 157±7 km s -1 and a V/σ ratio of 6, best fit with an arctangent model. This galaxy is included in the Mason et al. 2017 sample of lensed galaxies (ID 394) and from Hα kinematics the authors derive a maximum velocity of 133±22 km s -1 , compatible within errors with our estimate. The authors also derive an intrinsic velocity dispersion of 12±15 km s -1 from the integrated spectra, which is in good agreement with our intrinsic dispersion velocity value derived from the 2D map (15±7 km s -1 ).

A2667-sys1

The redshift of this galaxy was firstly derived by Covone et al. 2006 using VIMOS IFU observations. Using HST and IRAC data, Cao et al. 2015 derived its SED and source plane morphology revealing a disc with visible spiral arms, and a negative radial colour gradient, most likely due to a stellar metallicity gradient. Yuan et al. 2012 used OSIRIS data to investigate the resolved metallicity of this galaxy through the [N ii ]/Hα ratio. The central part of the galaxy, that seems to be unaffected by shocks, yields a metallicity of 12 + log(O/H) = 8.57±0.03, with a global value of 12 + log(O/H) = 8.606±0.051 for the entire galaxy, using the Pettini & Pagel 2004 calibrations. The authors also measure a radially increasing [N ii ]/Hα ratio, but observe that there is evidence of galactic outflows in the outskirts of the galaxy and estimate that shock ionisation (originated by the outflows) contributes significantly to the excitation of [N ii ]. Yuan et al. 2012 conclude that the observed ratio gradient is probably due to the presence of shocks in the outer parts of the galaxy and not by a positive gas-phase metallicity gradient. From the MUSE data, we derive a high-metallicity of 12 + log(O/H) = 9.08 , in strong disagreement with the 12 + log(O/H) = 8.6 value previously derived by Yuan et al. 2012. Some of this difference can be explained by the different calibrations used. Averaging the [N ii ]/Hα ratio from Yuan et al. 2012 over the several spatial regions (Box 1 to 4 in their figure 4) we obtain a value of 0.37, that converted to metallicity using the Maiolino et al. 2008 calibrations yields a value of 12+log(O/H) = 8.8 , which is still lower than the value derived in this work, but considerably higher (∼0.2 dex) than the one obtained with Pettini & Pagel 2004 calibrations. The remaining difference may be explained by the different dust corrections applied in both works. Yuan et al. 2012 derive a mass of 1.90±2.0×10 10 M , that agrees within errors with the value of 1.37 ± 0.21 × 10 10 M derived in this work.

On the kinematic side, Yuan et al. 2012 derive a V max sin(i)/σ ratio of ∼0.73, which, using the inclination derived in this work (∼21), gives an intrinsic ratio of about 1.05 -placing it at the limit of rotation-dominated galaxies. In contrast, we obtain a V max /σ of 3. As for A2390, this difference can come from the fact that MUSE covers a larger area than OSIRIS (see figure 1 of Yuan et al. 2012). Comparing the best 2D kinematic model with the observed one, it seems that the small twist on the zero-velocity curve is not due to lensing effects and might indicate the presence of a small bar. On the velocity dispersion maps corrected for beam smearing ( Fig. 4) some residuals are seen (up to ∼ 60 km s -1 ) but show no significant 2D structure.

M1206-sys1

MACS1206, nicknamed the cosmic snake (Cava et al. 2018), is another example of the "head and tail" structure: it has a long and highly magnified (although not complete) multiple image and a less magnified counter-image, that contains the total spatial information of this galaxy (see Eichner et al. 2013 for details). From the MUSE and HST data, we derive a stellar mass (7.87±0.49 × 10 10 M ) and, for the first time, the gas-phase metallicity (12 + log(O/H) = 8.90 ± 0.06 ) and the star formation rate of (38-43 M yr -1 ). This is also a rotation-dominated galaxy, displaying a smooth velocity field. Cava et al. 2018 performed SED fitting to the HST data and IRAC photometry, deriving a stellar mass of 4.0 × 10 10 M and a SFR SED of 30 M yr -1 .

A521-sys1

We present here the first study of the resolved properties in this system. This is a remarkable object, that displays clear spiral arms and numerous star-forming clumps (Richard et al. 2010b). It is also a very good example of a typical disc galaxy at z ∼ 1, with mass and SFR that agree quite well with the main-sequence of star-forming galaxies at these redshifts.

M1149-sys1

Finally, MACS1149, the highest redshift galaxy in the sample, is a spectacular spiral seen in 4 different multiple images. From the MUSE data, only the [O ii ] doublet can be seen (besides the Mg ii and Fe ii emission and absorptions, studied by Karman et al. 2016). Also here, the arctangent model provides a better fit to the observed kinematics of the 3 multiple images. Yuan et al. 2011, using OSIRIS data, derived the velocity field and metallicity gradient of this galaxy, concluding that it drops from 12 + log(O/H) = 8.54 ± 0.04 at the centre of the galaxy to 8.05 ± 0.05 in the edges, with a global value of 8.36 ± 0.04. More recently, Wang et al. 2017 also derived the metallicity of this galaxy using HST -grism data (ID 04054), arriving to a value of 12 + log(O/H) = 8.70 ± 0.11 , that we adopt in this paper.

This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .

 1 Figure 1. Gravitational arcs sample. For each galaxy: No frame: HST composite image with filters F435W, F814W and F160W in blue, green and red, respectively (except A521-sys1 where F606W is green and blue, and NIRC2-K band in red). The photometric aperture is plotted in dashed red. Yellow frame: corresponding MUSE [O ii ] pseudo-narrow band of the same sky region. The spectroscopic aperture is plotted in thick blue. Blue frame: F814W band source reconstruction. Both the photometric (HST ) and spectroscopic (MUSE) apertures are plotted in red and blue, respectively. Blue frame: F160W band source reconstruction. The magenta cross marks the morphological centre of the galaxy, and the ellipse is placed at 1 effective radius with the inclination and position angle derived from the morphological fit.

Figure 2 .

 2 Figure 2. Example of simultaneous determination of metallicity and extinction for galaxy AS1063-arc. The marginalised distributions for metallicity and extinction are displayed in the top left and bottom right panel, respectively. The bottom left panel displays the joint distribution. Similar plots for the remaining galaxies can be found in Appendix D.

  their ratios relative to Hα are given by Case B theory at electronic temperature T e =10000K and low electron density. For the [O ii ], the metallicity dependentKewley et al. 2004 calibration is used. The SFR from Balmer lines is obtained by drawing 500 random values of Hβ or Hγ fluxes (assuming a Gaussian distribution) and correcting them for dust attenuation using a sample of 500 τ v values drawn from the marginalised distribution obtained from the metallicity and attenuation fit. The dust corrected line fluxes are converted to SFR using the Kennicutt 1998a calibration and the SFR taken as the median of the sample. For the SFR estimated from [O ii ], an equivalent process is used. We produce a sample of [O ii ] fluxes, τ v and metallicity values and apply theKewley et al. 2004 calibrations. Finally, both SFR distributions are convolved with the magnification error, since they depend on the absolute flux of the emission lines. We list the calculated SFR in Table

SDSSFigure 3 .

 3 Figure3. Top panel: Gravitational arcs sample plotted against the main sequence. We plot in solid grey lines the parametrisation ofWhitaker et al. 2012 for this relation at z = 0.6 and z = 1, with the dispersion in grey. Bottom panel: Distance to the Fundamental Metallicity Relation (FMR) fromMannucci et al. 2010 (with the 1σ range of the SDSS sample form the same work shown in grey). Most galaxies lie close to this fundamental relation. The same colour scheme is used in both plots to identify objects. Both relations considered, these galaxies are typical star-forming galaxies for z ∼1.

Figure 4 .

 4 Figure 4. Observed and model velocity maps of the 8 gravitational arcs. For each galaxy, the 3 first upper panels display the velocity maps: the observed velocity derived from [O ii ] emission, the best model (with the model name) and the residuals. On the first panel, a black circle gives the seeing size for that observation. The 3 lower panels are the equivalent for the velocity dispersion. Both colour bars are in in units of km s -1 .

Figure 5 .

 5 Figure 5. Comparison of the velocity dispersion of the gravitational arcs with other samples. Velocity dispersion against stellar mass on the top panel and velocity dispersion against star formation surface density on the bottom (M1149-sys1 is omitted from this plot since we do not derive its SFR). The mean value for the local sample of rotating galaxies from Epinat et al. 2008 in plotted as a black line, with scatter masked by the grey area; the local (z = 0.1) DYNAMO sample of Green et al. 2014 in black squares (we plot only the rotating disc galaxies); the z ∼ 0.2 -1.4 sample from Contini et al. 2016 in black triangles; a sub-set of the z = 1 KMOS 3D galaxies, with velocity dispersions derived by Di Teodoro et al. 2016 as black crosses (only in the top panel); the lensed galaxies from the KLASS sample (Mason et al. 2017) are plotted as black diamonds (we have excluded possible mergers and galaxies with only upper limits). The KMOS Redshift One Spectroscopic Survey Sample at z∼1 from Johnson et al. 2018 is also plotted in grey dots. The gravitational arcs sample velocity dispersions are in good agreement with all these samples.

  of Bacon et al. 2017), which translates in

Figure 7 .

 7 Figure 7. Resolved extinction and SFR density maps for AS1063-arc (top) and A370 (bottom).In each row, from left to right: attenuation (τ) map, star formation rate surface density, from Hβ (units of [M yr -1 kpc -2 ]) and 1D radial profiles of both. These were obtained by averaging the values of extinction and SFR within distorted elliptical annuli of increasing radius (see Fig.1). These distorted elliptical annuli correspond to circular apertures in source plane. from the centre of the galaxy. The position of the strongest star forming clumps are marked with a black arrow.

Figure 8 .

 8 Figure 8. Velocity dispersion vs. SFR density for several subgalactic samples: in purple and blue circles, the results from AS1063 and A370, per Voronoi bin; in brown squares, the Genzel et al. 2011 giant star-forming clumps at z ∼2; in green triangles the sub-kiloparsec clumps at 1 < z < 4 of Livermore et al. 2015 and in dark triangles the clumps from Swinbank et al. 2012b. The lines show the results of the fit. Purple and blue full lines are fits to AS1063-arc and A370-sys1, respectively, dotted line the fit only to the Livermore et al. 2015, Swinbank et al. 2012b and Genzel et al. 2011, and dashed dotted line fit to all data.

Figure 9 .

 9 Figure 9. [O ii ] surface brightness and velocity dispersion for the full sample except M0416-sys28 where not enough data was available.The Spearman coefficients of these two quantities are plot in the upper left corner of each panel. Unreliable points, (in black circles in the plots, see text for details) were not included in the calculation of the Spearman coefficients. There is no correlation between [O ii ] surface brightness and velocity dispersion in the sample.

Figure C1 .

 C1 Figure C1. SED Fits. Best model in black (both circles and line) and data in colour. Flux in unit is maggies, linear flux units used in SDSS-III Collaboration (1 maggie = 3.631×10 3 Jy). Wavelengths are in Å and rest frame.

Table 1 .

 1 List of gravitational arcs. Coordinates α and δ correspond to the complete image position (see Section 3.1). The Point Spread Function (PSF) FWHM was obtained by fitting a MUSE pseudo F814W image with a seeing convolved HST F814W image. The magnification factor µ MUSE is the mean amplification factor within the MUSE aperture, predicted by the lensing model in Ref µ. This is merely indicative, since spectra were corrected pixel by pixel.

	Object	MUSE	α	δ	Exp. PSF	z	Ref z	Size µ MUSE	Ref µ
		Program	(J2000.0)		(hrs)	( )			( 2 )
	AS1063-arc 060.A-9345 (a) 22:48:42 -44:31:57 3.25 1.03 0.611 Gómez et al. 2012	33	4±1 Clément et al. in prep.
	A370-sys1	094.A-0115,	02:39:53 -01:35:05 6.0	0.70 0.725 Soucail et al. 1988	30 17±1 Lagattuta et al. 2017
		096.A-0710
	A2390-arc	094.A-0115	21:53:34 +17:41:59 2.0	0.56 0.913	Pello et al. 1991	31 10±1	Pello et al. in prep.
	M0416-sys28 094.A-0115	04:16:10 -24:04:16 2.0	0.45 0.940 Caminha et al. 2017 15 29±2 Richard et al. in prep.
	M1206-sys1 095.A-0181,	12:06:11 -08:48:05 3.0	0.51 1.033 Ebeling et al. 2009 50 18±1	Cava et al. 2018
		097.A-0269
	A2667-sys1 094.A-0115	23:51:39 -26:04:50 2.0	0.62 1.033 Covone et al. 2006	89 30±2	Pello et al. in prep.
	A521-sys1	100.A-0249	04:54:06 -10:13:23 1.67 0.57 (c) 1.043 Richard et al. 2010b 33 40±3	Richard et al. 2010b
	M1149-sys1 294.A-5032 (b) 11:49:35 +22:23:45 4.8	0.57 1.491 Smith et al. 2009	30	9±1	Jauzac et al. 2016

Table 2 .

 2 Best fit morphology parameters, from fitting a 2D exponential profile to the lensed reconstructed F160W HST images.

	Centre position (α, δ), effective radius (R e ), inclination (inc.) and
	position angle (PA) for the reconstructed source. Position angle
	0 corresponds to North and +90 to East.
	Object	α	δ	R e	inc.	PA
		(deg)	(deg)	(kpc) (deg) (deg)
	AS1063-arc 342.17851383 -44.53256249 4.6±0.1 52±15 -33
	A370-sys1	39.97151460 -1.57945147 7.2±0.4 41±12 -48
	A2390-arc	328.39577581 17.69923321 10.4±0.9 51±27 -77
	M0416-sys28 64.03737998 -24.0710587 4.1±0.2 57±15	72
	A2667-sys1 357.91536922 -26.08347133 2.4±0.1 44±21	54
	M1206-sys1 181.55001582 -8.80075740 6.6±0.1 64±9	10
	A521-sys1	73.52840603 -10.22305071 10.2±0.8 51±24	47
	M1149-sys1 177.40341077 22.40244261 15.4±0.9 51±43 -28

  (SFR [O ii ] ) relations are listed. M1149-sys1 was not included in this analysis, since only [O ii ] is present in the MUSE spectra.

	Object	Line Gas Metallicity Ratios (12 + log(O/H))	τ v	SFR Balmer SFR [O ii ] (M yr -1 ) (M yr -1 )
	AS1063-arc 1...10	8.82±0.02	1.09±0.12 41.5±4.0 50.3±10.1
	A370-sys1	1...10	8.88±0.02	0.44±0.11 3.1±0.3	3.1±0.6
	A2390-arc	1,9,10	9.00±0.11	0.60±0.40 7.3±2.5	7.9±6.1
	M0416-sys28 1,4,9,10	8.72±0.6	0.14±0.18 2.0±0.7	1.8±1.0
	A2667-sys1 1,4,9,10	9.04±0.04	0.53±0.23 15.7±3.7 9.8±4.2
	M1206-sys1 1,4,9,10	8.91±0.06	0.74±0.33 107.3±30.7 85.1±55.5
	A521-sys1 1,4,9,10	9.05±0.08	2.86±0.50 17.4±8.3 30.2±31.8

  derive a metallicity of 12 + log(O/H) = 8.36 ± 0.04 using the [N ii ] / Hα ratio, while Wang et al. 2017 measure 12 + log(O/H) = 8.70 ± 0.10 using HST grism spectra. Wang et al. 2017 use M08 calibrations, including line ratios with [O ii ], [O iii ] and Hβ emission lines, and calculate metallicity using a Bayesian inference framework, comparable to what is done here. We make use of both their metallicity and SFR (16.99 ± 4.3 M /yr) values derived from emission lines throughout this work.

Table 6 .

 6 Star formation rate surface density (Σ

			+0.39 -0.44	8.86 +0.27 -0.78	4.74 +0.46 -0.37	29.59 +3.12 -4.11	-0.48 +0.04 -0.02	0.00 +0.01 -0.00	0.91 +0.05 -0.06
		A370-sys1	2.49 +0.04 -0.04	7.10 +0.40 -1.47	7.27 +0.02 -0.24	5.68 +0.20 -0.18	-0.49 +0.01 -0.01	0.00 +0.00 -0.00	0.65 +0.02 -0.07
		A2390-arc	2.39 +0.50 -0.30	3.06 +1.65 -0.61	3.28 +0.44 -0.32	11.58 +2.86 -1.93	-0.41 +0.05 0.13	3.45 +0.41 -0.47	0.40 +0.16 -0.10
		M0416-sys28	0.95 +0.25 -0.25	6.16 +0.67 -0.47	6.46 +0.03 -0.07	2.47 +0.65 -0.64	-0.98 +0.03 -0.01	0.21 +0.03 -0.03	0.02 +0.03 -0.01
		A2667-sys1	1.37 +0.21 -0.19	9.57 +0.32 -0.82	3.72 +0.07 -0.09	6.36 +1.14 -0.93	-0.98 +0.01 -0.01	0.86 +0.03 -0.03	0.00 +0.01 -0.00
		M1206-sys1	7.87 +0.49 -0.47	2.07 +0.96 -1.31	4.38 +2.62 -2.29	39.91 +2.80 -2.69	0.07 +0.19 -0.19	0.24 +0.06 -0.16	0.26 +0.10 -0.08
		A521-sys1	3.21 +0.77 -0.64	9.63 +0.30 -4.73	6.00 +0.00 -0.01	8.81 +2.05 -1.82	-0.50 +0.01 -0.01	0.00 +0.01 -0.00	1.27 +0.07 -0.19
		M1149-sys1	1.76 +0.11 -0.11	7.10 +0.40 -1.47	7.27 +0.02 -0.24	18.83 +2.41 -2.27	-1.00 +0.01 -0.01	0.00 +0.01 -0.00	0.65 +0.02 -0.07
	Object	Σ SFR	Σ g	M g, KS	f g, KS
		(M yr -1 kpc -2 )	(M pc -2 ) (10 10 M )
	AS1063-arc	0.22±0.01	124.66±5.25 2.53±0.11 0.22±0.03
	A370-sys1	0.03±0.01	33.87±1.58 0.31±0.01 0.11±0.01
	A2390-arc	0.04±0.01	107.72±26.98 0.42±0.10 0.15±0.13
	M0416-sys28	0.02±0.01	21.49±5.63 0.35±0.09 0.27±0.15
	A2667-sys1	0.23±0.09	139.24±38.30 0.58±0.16 0.30±0.09
	M1206-sys1	0.70±0.39	246.00±93.41 2.16±0.83 0.22±0.07
	A521-sys1	0.03±0.01	29.78±10.03 1.74±0.59 0.35±0.12

SFR ), gas surface density (Σ g ) and gas mass (M g, KS ) estimated through the Kennicutt-Schmidt law and the respective gas fractions ( f g, KS ) calculated as M g, K S /(M g, KS + M )

  0 and Σ SFR are weakly correlated.

	Data	ρ		γ	χ 2 /do f
	Genzel et al. 2011	-	66	0.1±0.04	-
	Only Clumps	-0.12 68±2 -0.044 ± 0.023	2.83
	AS1063-arc	-0.21 61±1	0.004±0.012	6.64
	A370-sys1	0.18 106±15 0.314±0.042	0.20
	All Data	0.24	62±1	0.006±0.008	4.56

Table 10 .

 10 Energy associated with turbulence d E

				d A turb	and star
	formation d E d A SF 4. The gas densities used are listed in Table 6. The radius cor-following Krumholz et al. 2017 equations 1 and
	responds to the radius of the circularised areas used to calculate
	the gas densities. The velocities (V ) at that radius were calculated
	for the best kinematic model of each galaxy. M1149-sys1 was not
	analysed, since the SFR cannot be reliably calculated from MUSE
	data.
	Object	radius (kpc) ( km s -1 ) (erg cm -2 ) (erg cm -2 ) d E d E V d A SF d A turb	E turb E SF
	AS1063-arc	17	280	6.83e+12 2.48e+12 0.36
	A370-sys1	11	162	6.37e+11 1.69e+11 0.26
	A2390-arc	7	217	2.00e+12 8.40e+11 0.42
	M0416-sys28 15	105	1.76e+11 6.67e+09 0.04
	A2667-sys1	8	159	3.25e+12 6.97e+11 0.21
	M1206-sys1	11	205	7.94e+12 1.22e+12 0.15
	A521-sys1	30	55	6.00e+12 3.32e+11 0.06

dE d A tur b is the energy per unit area in [erg cm -2 ] associated with turbulence in a disc and dE d A SF the energy associated with star formation; Σ g is the gas surface density; σ the velocity dispersion; Σ SFR the star formation rate surface density; r the radius and V the rotation velocity. (The MUSE Gravitational Arcs 21

  Integrated MUSE spectra of the sample. Intrinsic flux, corrected for magnification, is shown in thick back and pPXF fits to the continuum in red.APPENDIX A: SAMPLE INTEGRATED SPECTRASample integrated spectra and pPXF fit. Details about the pPXF fit are given inside the inset box: reduced χ 2 of the fit, number of templates used of a total of 448 and degree of multiplicative (mdegree) and additive (degree) polynomials added (see Section 3.3 for details). M1149-sys1 is not shown since in the MUSE spectra only [O ii] is seen (spectra and absorption line analysis inKarman et al. 2016).

	5000 2 /dof : 0.1768 templates: 12 / 448 mdegree: 1 degree: 6 AS1063-arc 5000 2 /dof : 1.901 templates: 13 / 448 6000 [OII] H10 [NeIII] Ca-K H7 H8 H9 6000 mdegree: 1 degree: 6 A2390-arc 5000 6000 FeII FeII Fe*II FeII MgII MgI 2 /dof : 1.924 H templates: 20 / 448 7000 H 7000 [OII] 7000 mdegree: 1 degree: 8 A2667-sys1 5000 6000 7000 obs [Å] 2 /dof : 1.927 templates: 14 / 448 mdegree: 1 degree: 6 A521-sys1 FeII* FeII 0 2000 4000 6000 0 100 200 300 400 500 600 0 200 400 600 800 0 20 Flux (10 20 erg cm 2 s 1 Å 1 ) Flux (10 20 erg cm 2 s 1 Å 1 ) Flux (10 20 erg cm 2 s 1 Å 1 ) 40 60 80 Flux (10 20 erg cm 2 s 1 Å 1 ) FeII MgII 4000 5000 6000 0.0 0.2 0.4 0.6 0.8 Maggies 1e 8 AS1063-arc (z=0.611) 8000 H [OIII] [OIII] 8000 [NeIII] H 8000 [OII] [NeIII] H8 H7 8000 [OII] [NeIII] H7 7000 8000 3800 4000 3000 4000 5000 6000 7000 H H H 0.00 0.25 0.50 0.75 1.00 1.25 Maggies 1e 9 A2390-arc (z=0.912) 3800 4000 2000 3000 4000 5000 6000 0 1 2 3 4 5 6 Maggies 1e 10 A2667-sys1 (z=1.034) 3800 4000 2000 3000 4000 5000 6000 7000 8000 Rest Wavelength [ ] 0.0 Figure A1. 3000 0.5 1.0 1.5 2.0 2.5 Maggies 1e 10 A521-sys1 (z=1.043) 3800 4000	9000 9000 9000 H 9000 H 9000 4200 4400 8000 4200 4400 7000 4200 4400 9000 10000 H 4200 4400	5000 2 /dof : 1.505 templates: 19 / 448 mdegree: 1 degree: 9 A370-sys1 MgII 5000 2 /dof : 1.305 templates: 8 / 448 mdegree: 1 degree: 4 MACS0416-sys28 5000 FeII FeII* FeII FeII* MgII 6000 6000 F435W 0 200 400 600 800 0 200 400 600 0 500 1000 1500 F606W F625W F775W F814W F850LP F105W F110W F125W 0.25 0.50 0.75 1.00 1.25 1e 9 F140W F160W 3000 0.00 F555W F814W F110W F125W F160W 2000 3000 0 1 2 3 4 1e 10 F450W F606W F814W F850lp F110W F160W 2000 3000 6000 templates: 15 / 448 7000 [OII] [NeIII] H7 H H8 7000 7000 [OII] Observed Wavelength (Å) H [NeIII] H8 [OII] H7 2 /dof : 1.003 mdegree: 1 degree: 6 MACS1206-sys1 A370-sys1 (z=0.725) 8000 H 8000 H H [OIII] [OIII] 8000 H9 [NeIII] H7 H H8 3800 4000 4000 5000 6000 7000 4000 5000 6000 7000 H 9000 9000 9000 4200 8000 M0416-sys28 (z=0.940) 3800 4000 4200 4000 5000 6000 7000 0.0 0.5 1.0 1.5 2.0 2.5 1e 9 M1206-sys1 (z=1.036) 3800 4000 4200 2000 3000 4000 5000 Rest Wavelength [ ] 0 1 2 3 1e 10 M1149-sys1 (z=1.488) 3800 4000 4200	4400 9000 8000 4400 4400 6000 4400	F435W F475W F606W F625W F814W F105W F125W F140W F160W F435W F475W F606W F625W F775W F814W F850LP F105W F110W F125W F140W F160W F435W F475W F606W F625W F775W F814W F850LP F105W F110W F125W F140W F160W F435W F606W F814W F105W F125W F140W F160W

MNRAS 000, 1-28 (2018)

https://projets.lam.fr/projects/lenstool/wiki

Unless stated otherwise, we refer to gas-phase oxygen abundance simply as metallicity.MNRAS 000, 1-28 (2018)

([O iii ]λ

+ [O iii ]λ 4959 + [O ii ]λ 3727)/HβMNRAS 000, 1-28 (2018)

Prospector publicly available at: https://github.com/bd-j/ prospector MNRAS 000, 1-28 (2018)

camel is available at https://bitbucket.org/bepinat/camel. git

Voronoi binning code: http://www-astro.physics.ox.ac.uk/ ~mxc/software MNRAS 000, 1-28 (2018)

Q = σκ/πGΣ gas , where σ is the velocity dispersion; κ the epicyclic frequency; G the gravitational constant and Σ gas the gas surface density)

MNRAS 000, 1-28 (2018)20 Patrício et al.

ACKNOWLEDGEMENTS

We warmly thank Benjamin Johnson for helping using prospector and updating the code to fit spectra and photometry simultaneously. We thank Jorge Sánchez Almeida and Andrea Ferrara for insightful suggestions about clumps, turbulence and energy. We thank Mark Swinbank, Philippe Amram, Miroslava Dessauges and Mathieu Puech for useful comments on this work. We thank David Bina for reducing A2390 and A2667 MUSE data. We also thank the referee for a constructive report with useful suggestions.

VP, DC, JR, JM, DL and BC acknowledge support from the ERC starting grant 336736-CALENDS. TC, NB and

BE acknowledge support from the ANR FOGHAR (ANR-13-BS05-0010-02) and the OCEVU Labex (ANR-11-LABX-0060). TC acknowledge support from and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the "Investissements d'avenir" French government program. KBS was supported by the BMBF grant 05A14BAC. JB acknowledges support by Fundação para a Ciência e Tecnologia (FCT) through national funds (UID/FIS/04434/2013) by FEDER through COMPETE2020 (POCI-01-0145-FEDER-007672), and by FCT

MUSE Gravitational Arcs 25

Table B1. AB magnitude of the available HST bands for each object. All are amplification corrected magnitudes. This factor µ, calculated in the same aperture, is also listed in the table. Errors on the magnitudes include cluster member subtraction residuals and a systematic error of 0.03 mag. Magnification errors were not added to these, since they are achromatic, and are listed separately in the err µ column. Table B2. Emission line fluxes in units of 10 -18 erg cm -2 s -1 . These were measured in the continuum subtracted spectra and errors were derived from the MC technique described in section 3.2. We note that the A2390-arc spectrum is more noisy and does not allow to accurately determine [Ne iii ] flux and fainter Balmer lines.

APPENDIX B: PHOTOMETRY AND LINE EMISSION FLUXES

Here we list the photometry derived from the publicly available HST data (programmes listed in the Table ) with respective errors and magnification errors. This data was used, together with the integrated spectra, to derive the mass and stellar population properties of this sample. We also list the observed emission lines fluxes, measured in the continuum subtracted spectra and corrected for magnification but not for attenuation, used to derive the metallicity and attenuation factors in Section 3.5.

APPENDIX C: SED FITTING RESULTS

The photometry and integrated spectra of the galaxies was simultaneously fit using the SED fitting code PROSPECTOR (see Section 3.5 for details). Emission lines were masked. Here we presented the comparison of the photometric and spectral data with the model predictions for both. The inset shows a zoom of the Balmer series.

APPENDIX D: METALLICITY AND ATTENUATION FITS

We present here the marginalised distributions and joint distribution of metallicity and attenuation for each galaxies A370-sys1, A2390-arc, M0416-sys28, A2667-sys1, M1206-sys1 and A521-sys1. The results for AS1063-arc are presented in the main body of this work. M1149-sys1 has only [O ii ] available within the MUSE spectral range, reason for which we do not perform this fit in this galaxy.

APPENDIX E: INDIVIDUAL PROPERTIES

Here we summarize the properties of the individual galaxies derived in this work and compare them with the available studies in the literature for each specific galaxy.

AS1063-arc

This galaxy is a beautiful example of a grand design spiral galaxy whose most striking feature is a gigantic and bright H ii region at the edge of one spiral arm. The [O ii ] kinematics of the galaxy was briefly studied by Karman et al. 2015, based on the same MUSE observations used in this paper. They found a smooth velocity field with a symmetrical rotation with maximum velocities of up to 200 km s -1 . We confirm this result and also derive the velocity dispersion map that, after being corrected for beam-smearing effects, still displays a higher central velocity dispersion. Despite the presence of several star forming clumps visible in the HST imaging, none of these coincide with any increases or decreases in the velocity dispersion (Fig. 4). This is true with or without the beam smearing correction.

A370-sys1

A370-sys1, sometimes called "the dragon", was the first gravitational arc to be spectroscopically confirmed (Lynds & Petrosian 1986;Soucail et al. 1988). It extends over 23" in the sky and is composed of several multiple images: the "head" of the dragon is the complete image of the galaxy, and has a somewhat low (µ ∼10) magnification factor, while the "tail" (most elongated structure) is composed by two highly magnified (µ ∼35) multiple images, mirroring each other, and that are only a fraction of the complete galaxy (details in Richard et al. 2010b). Blue spiral arms are clearly seen in the complete image, with a dozen of strongly star forming regions populating them. As for AS1063-arc, the map shows no particular features at the location of these clumps. The analysis of kinematics also reveals a regular gas rotation, with a twist of the zero-velocity usually associated to the presence of a bar, unknown before this analysis.

A2390-arc

A2390-arc has been extensively studied in the past both spectroscopically [START_REF] Osterbrock | Astrophysics of gaseous nebulae and active galactic nuclei[END_REF]) and photometrically (Cowie et al. 2002;Rigby et al. 2008;Olmstead et al. 2014).

Its kinematics was studied by Swinbank et al. 2006, using VIMOS IFU data. The authors derived a maximum velocity (corrected for inclination) of ∼187 km s -1 , with an inclination of 69 degrees. However, the authors divide this arc in two parts (arcA and arcB, in their figure 3), with arcB corresponding to the z ∼ 1 object seen overlapping with the arc (see Fig. 1). The MUSE data also allows to confirm that there is a foreground object at z ∼ 1, but we can still reliably measure the velocity field of A2390-arc at this point.

The full arc has a smooth velocity field that is well described by a single velocity model, so we do not divide this arc in two objects as in Swinbank et al. 2006. The derived kinematic parameters are at odds with the previous analysis by Swinbank et al. 2006: we derive a lower inclination, ∼57 degrees, and consequently higher maximum velocity of ∼220 km s -1 . The better quality of the MUSE data, that allows larger radii to be explored as well as the fact that we do not divide this object in arcA and arcB, may explain this difference. We derive a mass of 2.39 ± 0.5 × 10 10 M , with the same order of magnitude as the one derived by Olmstead