The Lensed Lyman-Alpha MUSE Arcs Sample (LLAMAS) I. Characterisation of extended Lyman-alpha haloes and spatial offsets.

Aims. We present the Lensed Lyman-Alpha MUSE Arcs Sample (LLAMAS) selected from MUSE and HST observations of 17 lensing clusters. The sample consists of 603 continuum-faint (-23 < M UV < -14) lensed Lyman- α emitters (producing 959 images) with secure spectroscopic redshifts between 2.9 and 6.7. Combining the power of cluster magniﬁcation with 3D spectroscopic observations, we are able to reveal the resolved morphological properties of 268 Lyman- α emitters. Methods. We use a forward modelling approach to model both Lyman- α and rest-frame UV continuum emission proﬁles in the source plane and measure spatial extent, ellipticity and spatial o ﬀ sets between UV and Lyman- α emission. Results. We ﬁnd a signiﬁcant correlation between UV continuum and Lyman- α spatial extent. Our characterization of the Lyman- α haloes indicates that the halo size is linked to the physical properties of the host galaxy (SFR, Lyman- α equivalent width and Lyman- α line FWHM). We ﬁnd that 48% of Lyman- α haloes are best-ﬁtted by an elliptical emission distribution with a median axis ratio of q = 0 . 48. We observe that 60% of galaxies detected both in UV and Lyman- α emission show a signiﬁcant spatial o ﬀ set ( ∆ Ly α − UV ). We measure a median o ﬀ set of ∆ Ly α − UV = 0 . 58 ± 0 . 14 kpc for the entire sample. By comparing the spatial o ﬀ set values with the size of the UV component, we show that 40% of the o ﬀ sets could be due to star-forming sub-structures in the UV component, while the larger o ﬀ sets (60%) are more likely due to larger distance processes such as scattering e ﬀ ects inside the circumgalactic medium or emission from faint satellites or merging galaxies. Comparisons with a zoom-in radiative hydrodynamics simulation of a typical Lyman- α emitting galaxy show a very good agreement with LLAMAS galaxies and indicate that bright star-formation clumps and satellite galaxies could produce a similar spatial o ﬀ sets distribution.


Introduction
The existence of bright Lyman-α radiation emitted by galaxies was originally predicted by Partridge & Peebles (1967) and progressively became a prominent target in searches for high redshift galaxies. Then, the detection of extended Lyman-α emission around high redshift galaxies was predicted by Haiman et al. (2000). The Lyman-α emission is now also used as a neutral hydrogen gas tracer in the circumgalactic medium (CGM) given its resonant nature, which produce extended halos surrounding galaxies up to 30 kpc (Matsuda et al. 2012;Momose et al. 2014), with the exception of quasars.
The origin of extended Lyman-α haloes surrounding galaxies is still unknown with two main hypotheses being considered: scattering of Lyman-α photons, produced mostly within star-forming regions, through the interstellar (hereafter ISM) and circum-galactic medium or in-situ photoionisation emission or collisional emission in the CGM (Mitchell et al. 2021). Lyman-α haloes therefore represent a powerful probe of the hydrogen gas within the CGM, tracing both spatial extent and velocity structure of the gas surrounding galaxies and thus investigating the galaxy formation processes and the reionisation epoch at z = 6. These resonant scattering events increase the path length of Lyman-α photons and consequently the observed Lyman-α emission is potentially influenced along this path by a large number of physical parameters (column density, temperature, dust content, kinematics, covering fractions and clumpiness), which modify both the spectral profile of the line and its spatial distribution (Ouchi et al. 2020).
Observing and characterizing Lyman-α haloes is crucial to understand the nature of the CGM at low and high redshifts. The LARS collaboration (for Lyman-α Reference Sample, Hayes et al. 2013Hayes et al. , 2014Östlin et al. 2014) characterised the Lyman-α emission in low redshift star-forming galaxies (z < 0.5), demonstrating the presence of a complex structure in the outer parts of the disks. At high redshift (z > 2), conducting similar studies is harder because of limitations in sensitivity and spatial resolution. Narrow-band imaging observations noted that z > 2 galaxies appear more extended in Lyman-α than in the rest-frame UV continuum (Møller & Warren 1998;Fynbo et al. 2001). Hayashino et al. (2004) later detected the first extended halo by stacking 74 Lyman-α emitters (hereafter LAEs) at z = 3.1 using the Subaru Telescope, while Rauch et al. (2008) also observed 27 faint LAEs between z = 2.6 and 3.8 in very deep long slit exposures with ESO-VLT, finding that the majority of their LAEs had spatial profiles larger than the UV sources. More recently, Steidel et al. (2011) stacked the images of 92 bright Lyman Break Galaxies (LBGs) at z = 2.3 − 3, and demonstrated that Lyman-α emission is detected out to 10" from the UV continuum emission at a surface brightness level of ∼ 10 −19 erg s −1 cm −2 arcsec −2 . This stacking approach was adopted by other groups (Matsuda et al. 2012;Feldmeier et al. 2013;Momose et al. 2014;Xue et al. 2017;Wu et al. 2020) and confirmed the presence of extended Lyman-α haloes with typical exponential scale lengths of ∼ 5 − 10 kpc. However, a study of the possible diversity among individual Lyman-α haloes remained out of reach in these observations due to the lack of sensitivity.
The arrival on sky of the ESO-VLT instrument MUSE (the Multi Unit Spectroscopic Explorer, Bacon et al. 2010), an integral-field spectrograph with an unrivalled sensitivity, has substantially increased the number of observed LAEs at high redshift (z > 3) (Bacon et al. 2015). The sample was extended with the MUSE Hubble Ultra Deep Field (UDF, Bacon et al. 2015;Leclercq et al. 2017) and the MUSE Wide Survey (Wisotzki et al. 2016). Leclercq et al. (2017) reported the detection of individual extended Lyman-α emission around 80% of the LAEs detected, confirming the presence of a significant amount of hydrogen gas in the CGM. They presented a systematic morphological study of 145 of these Lyman-α haloes in the UDF, showing that the majority of the Lyman-α flux comes from the halo surrounding each galaxy, whose properties seem to be related to the UV properties of the galaxies. In parallel, the Integral Field Unit (IFU) instrument KCWI (Keck Cosmic Web Imager, Morrissey et al. 2018) at the Keck Observatory, recently started to confirm similar results in 2 < z < 3 LAEs (Erb et al. 2018;).
Scenarios about the origin of the extended Lyman alpha emission include cooling radiation from cold accretion, outflows, emission from satellites galaxies and resonant scattering of photons produced in the ISM and the CGM (Laursen & Sommer-Larsen 2007;Steidel et al. 2010;Zheng et al. 2011). Spatially integrated Lyman-α emission is almost always redshifted relative to the systemic velocity, indicating the presence of galactic outflows in both observational and theoretical works (Heckman 2001;Verhamme et al. 2006;Scannapieco 2017;Song et al. 2020). Cosmological simulations predict also the presence of inflowing filamentary streams of dense gas (Kereš et al. 2005;Dekel & Birnboim 2006;Mitchell et al. 2021;Byrohl et al. 2021) which could produce an overall blue-shifted Lyman-α line (Dijkstra et al. 2006;Verhamme et al. 2006;Mitchell et al. 2021;Garel et al. 2021). Galaxies with a relatively low Lyman-α optical depth typically exhibit double-peaked profiles, with a dominant red peak, with the peak separation strongly dependent on the kinematic of the gas and neutral hydrogen column density (Verhamme et al. 2006;Henry et al. 2015;Gronke et al. 2016;Verhamme et al. 2017). In addition, emission from satellite galaxies and in-situ emission probably contribute to the spatial extent and the clumpy morphology of the Lyman-α haloes (Mas-Ribas & Dijkstra 2016;Mas-Ribas et al. 2017;Mitchell et al. 2021;Byrohl et al. 2021 in simulations).
Although the spatial resolution of the MUSE deep fields in Leclercq et al. (2017) data did not allow them to distinguish between the aforementioned scenarios, by comparing these data with a zoom-in cosmological simulation, Mitchell et al. (2021) demonstrated that simulated Lyman-α haloes are likely powered by a combination of scattering of galactic Lyman-α emission, in-situ emission of the (mostly infalling) CGM gas and Lyman-α emission from small satellite galaxies. In their simulation, Mitchell et al. (2021) showed that each of the scenarios is dominant on a different scale. Another approach to study the physical processes influencing the Lyman-α emission is to use simple wind models, in which the central Lyman-α source is surrounded by an expanding neutral hydrogen medium associated with dust. Albeit simple in nature, these models can successfully reproduce the majority of observed Lyman-α line profiles (Ahn 2004;Schaerer et al. 2011;Gronke et al. 2015;Yang et al. 2016;Gronke 2017;Song et al. 2020).
One outstanding feature of Lyman-α haloes is the spatial offset between UV continuum and Lyman-α emission peaks repeatedly reported in the litterature Hoag et al. 2019;Lemaux et al. 2020;Ribeiro et al. 2020). Most studies so far used long-slit spectroscopy to measure offsets between UV and Lyman-α emission, but these observations are restricted to one dimension only. Taking advantage of the 3D spectroscopy, IFUs are a more efficient tool to provide a more thorough and detailed picture of the Lyman-α and UV emission in individual galaxies. The presence of spatial offsets could indicate that Lyman-α photons can preferentially be produced (or scattered) far away from the star forming regions that are responsible for the UV emission. Depending on the range of offsets observed they could preferentially support one of the two main scenarios mentioned previously. For example, an offset smaller than the UV source extent could indicate an off-center star-formation clump emitting a large amount of Lyman-α photons. On the contrary an offset larger than the UV source size will support the scenario of satellite galaxy emission or resonant scattering from escape channels. Measuring these spatial offsets very precisely and comparing them with the UV emission distribution represent a great opportunity to study the preponderance of different scenarios. hydrogen column density in the CGM at z = 4.1. Similar studies were performed by Erb et al. 2018 (a double peaked Lyman-α line at z = 2.3), Vernet et al. 2017 (at z = 3), Matthee et al. 2020,b (an LBG at z = 6.53 and one galaxy at z = 6.6) and Herenz et al. 2020 (one Lyman-α blob at z = 3.1). Leclercq et al. (2020) performed a similar study of six bright haloes from the total UDF sample, searching for resolved variations of the Lyman-α line profile across the halo and correlations with host galaxy properties. They showed that the Lyman-α line is in general broader and redder in the extended part of the halo, suggesting that Lyman-α haloes are powered either by scattering processes in an outflowing medium. However the lack of sufficient spatial resolution (typically ∼ 3 − 4.9 kpc at z = 4) of these studies left the physical interpretation widely open, and some questions are still pending, namely: What are the origins of the Lyα photons? Which physical mechanisms (e.g., outflows, inflows, satellites galaxies) are responsible for the extent of Lyman-α haloes and the spectral shapes of the Lyman-α line profiles?
In order to improve our understanding of the properties of the CGM and provide robust new constraints on theoretical models, we focus our observations on high-redshift lensed galaxies. Gravitational lensing boosts and magnifies the total observed flux of sources, causing them to appear physically larger and (in some cases) creating multiple images of a single object, making them ideal targets for spatially resolved studies. A small but growing number of highly magnified LAEs have already been individually studied in order to characterize the CGM gas at 2 < z < 7. One of the first is a lensed galaxy at z = 4.9 presented in Swinbank et al. (2007), showing extended Lyman-α emission. Following this study, many subsequent efforts have also targeted strongly lensed sources (Karman et al. 2015;Caminha et al. 2016;Patrício et al. 2016;Vanzella et al. 2016;Smit et al. 2017;Claeyssens et al. 2019;Vanzella et al. 2020;. These studies, focusing only on one or two objects, have provided the first evidence of variations in the Lyman-α line profile across the halo, revealing the complex structure of the neutral hydrogen distribution surrounding galaxies (in terms of covering fraction, column density, presence of inflows/outflows). Besides these studies have demonstrated that lensing observations represent a privileged field to study the CGM at high redshift. However, most studies of lensed galaxies only concern too few objects to draw general conclusions.
With this in mind, we construct a statistically large sample of lensed Lyman-α emitters named the Lensed Lyman-α MUSE Arc Sample (hereafter LLAMAS). Based on the recent MUSE lensing clusters data release presented in Richard et al. (2021) (hereafter R21), totalling 141 hours were obtained mainly through MUSE Guaranteed Time Observation (GTO) on 18 different fields; we construct a unique sample of 603 lensed LAEs (forming 959 images) at z = 2.9 − 6.6. Their strong lensing properties are characterised using well-constrained models of massive galaxy clusters, with magnification values ranging from 1.4 to 40. A partial sample of these LAEs was used previously to study the lens models (e.g. Mahler et al. 2018;Lagattuta et al. 2019) and the Lyman-α luminosity function in individual / several clusters (Bina et al. 2016;de La Vieuville et al. 2019). The LLAMA sample makes use of the unique combination of strong lensing, deep MUSE IFU and HST high resolution images to study the spatial and spectral properties of the Lyman-α haloes in the entire sample. In this first publication of the series, we present the general properties of this sample and focus on the morphological parameters, with emphasis on the spatial extent of the Lyman-α emission and spatial offsets between UV continuum and Lyman-α emission. In order to physically interpret these results and disentangle between the many Lyman-α halo production scenarios, we compare our results with a zoom-in radiation-hydrodynamical simulation. The paper is organized as follows: we describe our data and the sample selection in Sect. 2. Section 3 presents our procedure for image construction, spectral extraction and fitting, and modeling of the UV and Lyman-α spatial distribution in the source plane. We describe the results in Sect. 4. We discuss these results and compare with the zoom-in simulation in Sect. 5. Finally we present our summary and conclusions in Sect. 6. All distances are physical and we use AB magnitudes. We adopt a Λ cold dark matter cosmology with Ω Λ = 0.7, Ω m = 0.3 and H 0 = 70 km s −1 Mpc −1 .

Lyman-α sample
The data and redshift catalogs presented in this study have all been processed following the method described in R21. Among the 17 galaxy clusters catalogs used for this work, 12 were presented in R21 and 5 are presented for the first time (MACS0451, MACS0520 are completely new and A2390, A2667 (de La Vieuville et al. 2019) and AS1063 (Mercurio et al. 2021) had been the subject of previous studies). This section will shortly describe the method presented in R21 and the LAEs selections based on the global catalogs.

MUSE data
The MUSE data reduction procedure details are given in R21, largely following the prescription described in Weilbacher et al. (2020) with some specific improvements for crowded fields. Of the 17 clusters we explore here, five (MACS0451, MACS0520, A2390, A2667 and AS1063) were analysed after the publication of R21, but the final data products and galaxy catalogs were constructed following exactly the same procedure. The final output is given as a FITS datacube with 2 extensions containing the flux and the associated variance over a regular 3D grid at a spatial pixel scale of 0.2" and a wavelength step of 1.25 Å between 4750 and 9350 Å. The final seeing, defined as the FWHM of the point spread function at 7000 Å, varies from 0.52" to 0.79" among the fields. Every cluster in our sample has a redshift between 0.2 and 0.6 and are all known to be massive strong lenses. The integration times vary between 2 and 14 hours per field, using a combination of standard and adaptive-optics (for observations done after 2014) modes. The fields of view are centered on the core regions of each cluster, in order to maximise the number of strongly lensed LAEs. The MUSE field of view is 1 × 1 arcmin 2 ; for five clusters, multiple contiguous MUSE pointings were mosaicked to cover the complete multiple image area of the clusters (between 2 and 4 pointings, see Table 1).

HST data
To complement the MUSE data we use the available highresolution ACS/WFC and WFC3-IR images in the optical / nearinfrared covering the MUSE observations. Six clusters in the sample are included in either the CLASH (Postman et al. 2012) or Frontier Fields (Lotz et al. 2017)

Input redshift catalogs
Based on these observations, a complete spectroscopic catalogue was constructed for each cluster. The complete procedure is described in R21 and the main steps are: -Production of an input photometric catalogue of continuum sources that overlap with the MUSE field of view. The detection image (produced by combining the HST images into an inverse-variance-weighted detection image) is given as input to the SExtractor software (Bertin & Arnouts 1996). The photometry performed by SExtractor in each band is merged together into a final catalogue of HST sources (hereafter PRIOR sources). -Independent of the HST catalogue, a line-detected sources catalogue is produced directly from the MUSE datacubes, performed by running the muselet software which is part of the MPDAF python package (Piqueras et al. 2019) hereafter MUSELET sources. -A spectrum of each source (both PRIOR and MUSELET) is extracted from the MUSE datacube based on weighted images (created for each source by taking the flux distribution over the segmentation maps produced as part of the detection process) to optimise the signal-to-noise of the detections. A local background spectrum around each source was estimated and subtracted (to remove large-scale contamination from bright sources, such as stars and cluster members, as well as potential systematics in the background level remaining from the data reduction) from the weighted spectrum to compute an optimized spectrum for source identification and redshift measurement. -All MUSELET and PRIOR sources down to a limiting signalto-noise ratio (S/N) in the MUSE continuum (averaged over the full MUSE wavelength range) were inspected individually and their redshift were assessed. For each source we determine the redshift value and confidence (between 1 for insecure redshifts and 3 to most secure redshifts, the details of redshift classification can be found in R21 Sect 3.5), as well as any association between PRIOR and MUSELET sources (i.e. when the same source is detected in both HST and MUSE data). -The resulting catalogs are then tested with the corresponding lenstool (Jullo et al. 2007) mass model of the cluster to associate multiple images together and predict potential new multiple images. -The final catalogs are composed of 4020 secure redshifts with 0 < z < 6.7 and 634 unique LAEs with redshift confidence > 1 (more details on the Lyman-α line identification can be found in R21 Sect. 3.5) 1 . The redshifts of the LAEs are measured based on the Lyman-α emission line or the presence of a strong Lyman break or from nebular lines when detected.

Mass models
In order to study the intrinsic galaxy properties of the LAEs in the source plane, we used the parametric models of the cluster mass distribution presented in R21, generated using the public lenstool software (Jullo et al. 2007). lenstool allows to generate parametric models of each cluster's total mass distribution, using numerous multiple images identified in the catalogs as constraints. The final model's parameters and constraints are presented in Appendix B of R21. Each cluster mass model is 1 https://cral-perso.univ-lyon1.fr/labo/perso/johan.richard/MUSE_data_release/ optimized with between 7 and 100 multiples systems of images with secure spectroscopic redshifts. The precision of the lens models, which corresponds to the typical spatial uncertainty in reproducing the strongly-lensed images, is typically from 0.5" to 0.9". One crucial value for the study of lensed background source morphologies is the lensing amplification and shear. We use the value from R21 as a first estimate, based on the central location of each image in the catalogs, and refine it in the Section 3.5. As lenstool uses a Monte Carlo Markov Chain (MCMC) to sample the posterior probability distribution of the model, statistical errors were estimated for each parameter of the models. A second important measurement derived from the lens model is the equivalent source-plane area covered by the MUSE observations. The intrinsic survey volume for lensing studies differs from the image-plane area due to strong lensing effects. The effective survey volume is reduced by the same amount as the magnification factor. This value varies depending on the cluster (due to the different mass distribution and MUSE coverage) and redshifts of the sources. At z = 4, which is the median redshift of all the LAEs detected in this sample, the total source-plane area covered is about 1.8 arcmin 2 . The relative contribution of each cluster to the full survey covolume is provided in the last columns of Table  1.

LAE selection
The catalogs presented in R21 include all the spectroscopic redshifts measured in each field. In order to construct a sample of LAEs, we selected all the sources with a secure redshift (confidence 2 or 3 based on multiple emission lines, a clear asymetric Lyman-α emission line, a clear Lyman break or a lensing confirmation of the high-redshift nature of the image) between 2.9 and 6.7 (1031 Lyman-α images selected on 1269 detected). We flagged all galaxies for which no Lyman-α emission is detected (e.g. with an Lyman-α line integrated S/N>3). For these sources, we searched for extended Lyman-α emission in multiple NB images produced around the predicted location of the Lyman-α line (based on the galaxies' systemic redshift) with different velocity windows. We rejected 20 galaxies with no significant emission features around the galaxy location (S/N<3). After a visual inspection we rejected all images detected in close vicinity of a bright cluster galaxy (BCG) or bright cluster members. In such cases the Lyman-α line is too contaminated by the foreground galaxy continuum and the HST emission is not sufficiently isolated to be correctly spatially fitted. In fact, all the images rejected are either central poorly magnified images or small counter-images from a multiple system. The rejection of these images does not affect the final sample results since either the images would have a too low S/N to be spatially fitted (see Section. 3.4) or they are parts of a multiple system of which the other images are kept in the sample. The final LLAMAS catalogue is composed of 959 Lyman-α images from 603 unique objects. Among these sources, 341 have at least one image with an HST detection. The number of LAEs detected in each cluster is presented in Table 1 together with observation information on each field. Figure 1 shows the global properties of the LLAMAS: redshift, lensing magnification, UV magnitude and UV star formation rate (SFR) distributions, the last two corrected for lensing magnification. The grey area of each panel represents sources selected for the spatial fitting (see Section 3.4). The median redshift of the complete sample is z = 4 with 40 galaxies at z > 6. The median magnification value is µ = 5.4, with a range from 1 to 40. The magnification values presented in Figure1 are the values computed in R21, based on the central position of each image. The UV absolute magnitude at rest-frame 1500 Å (M 1500 ) and the UV spectral slope (β) are estimated by adjusting a single power law to the HST broad-band photometry available in each cluster (see R21) redward of Lyman-α in the UV. The UV SFR is derived from M 1500 using the Kennicutt (1998) relation. Among the sample, 39% of the objects are pure MUSELET detections (without any UV detection in HST). For these objects, no UV magnitude or SFR can be computed. Among the PRIOR sources, 65% of the galaxies have SFR < 1 M ⊙ yr −1 , and the median value is SFR = 0.55 M ⊙ yr −1 .

Narrow-band image construction and global spectral properties
To characterize the Lyman-α emission of each LAE we first construct a 5"×5" Lyman-α narrow-band (hereafter NB) image of each object from the MUSE datacube. When a source produces multiple images, we study each image independently. With the intention of maximising the S/N in the NB images and recovering most of the Lyman-α emission of each source, we applied an iterative process consisting in three steps: spectral fitting, NB image construction, and spectral extraction repeated three times.
First we fitted the Lyman-α line with the asymmetric Gaussian function  given by: with A the amplitude of the line, λ peak the peak of the line, a the asymetry parameter and d the typical width of the line. Four parameters were optimised: peak position (λ peak ) of the line, FWHM, flux F, and asymmetry a of the profile. The prior on the flux corresponds to the integrated flux of the last produced spectrum on the spectral range [λ peak − 6.25 : λ peak + 6.25] Å. The prior on λ peak is based on the catalogue redshift of each source. When two spectral peaks were detected in the Lyman-α line, the spatial fit only took the red peak into account. We applied a uniform prior of FWHM and asymmetry with mean values of 7 Å and 0.20, respectively. We applied this fit on a spectral window around the Lyman-α emission peak, in which each pixel has a minimum signal to noise ratio of 2.5. Each spectrum and its associated variance were fitted using the python package emcee (Foreman-Mackey et al. 2013). We performed the fit with 8 walkers and 10000 steps. We used the median values of the resulting posterior probability distributions for all the model parameters as best-fit parameters. Errors on the parameters were estimated using the 16th and 84th percentiles. In the first loop of the process, this spectral fit was applied on the continuum optimised, sky-subtracted spectrum produced in the data release of R21.
Secondly we constructed a Lyman-α narrow-band image for each Lyman-α image. The central wavelength is based on the results of an asymmetric Gaussian function fit of the Lyman-α line. The continuum level was measured on the left and right side of the Lyman-α line, over a 28 Å width band and subtracted from the collapsed Lyman-α line image. The NB bandwidth [λ left : λ right ], with λ left and λ right being the left and right wavelengths of the band where the spectral layers were summed, was optimised (from λ right − λ left = 2.5 to 20 Å) to maximise the S/N as measured in a 0.7" radius circular aperture (typical MUSE PSF size) centred on each image (on the catalogue source position). The start and end position of the spectral band for continuum level estimation were chosen to be close to the Lyman-α line while not containing any Lyman-α Article number, page 5 of 22 A&A proofs: manuscript no. output 5.1. Modélisation de l'émission dans le plan source 115 FIGURE 5.2 -De gauche à droite : distribution en redshift, amplification gravitationnelle, magnitude UV (corrigée de l'amplification) et taux de formation stellaire (corrigé de l'amplification) des galaxies sélectionnées pour être modélisées dans le plan source (en gris) par rapport à l'échantillon total (en couleur).
Third, the new Lyman-α image was obtained, a non-weighted spectrum was re-extracted from the MUSE datacube based on this new NB image. The purpose of this new extraction is to get a Lyman-α optimised spectrum containing the most of the Lyman-α emission from the galaxy and its halo. We spatially smoothed the best NB image with a Gaussian filter of FWHM smooth = 0.4". The total spectra is the sum of the spectra of each MUSE pixel with a flux in the NB image higher that the typical value of the dispersion measured in the image. To avoid external contamination, a mask was created manually for each object to isolate them from possible cluster galaxies or star residuals in the NB images. Finally the sky subtraction was performed using the same sky spectrum used to create the first spectrum.
The same process was performed two more times, each time based on the last NB image and extracted spectra. If an object presents a Lyman-α emission too low for the spectral extraction (i.e. no pixels with a smoothed SB level > 6.25 10 −19 erg.s −1 .cm −2 .arcsec −2 in any of the tested NB images, the object was rejected from the sample (in total 12 images for 5 objects rejected).

Forward modelling
To study intrinsic morphological properties of the LAEs, we modeled the Lyman-α and UV continuum emission in the source plane making use of the cleanlens function from the latest version of lenstool. The method used in this study is the forward modelling approach based on parametric source models, including both lensing and instrumental effects. It consists in generating parametric source models (based on user assumptions on the emission profile) in the source plane, lensing it by the best cluster model, convolving with the PSF and re-gridding it to the spatial sampling of the observation, to compare with the observation. In this study we used only Sersic profiles and specifically exponential (Sersic index n = 1) as used in Wisotzki et al. (2016) and Leclercq et al. (2017). The free parameters are spatial location (x and y) of the center, exponential scale radius (a), ellipticity (ϵ = (a − b)/(a + b) with b the minor axis of the distribution), position angle (θ), magnitude (m) and Sersic index (n, fixed to 1 in case of exponential profiles). The best-fit parameters were found by minimizing the residuals between input observed images and simulated image-plane observations. One or several Sersic components could be used to reproduce one object. The multiple images of the same object could be fitted together or separately. In this study we decided to fit each multiple image separately to avoid multiple images models uncertainties effect on the source reconstruction. We applied this fit both on Lyman-α emission from NB images and UV continuum emission from HST images for objects with an HST detection in the catalogue. To isolate each object from other features presented in the NB images, we constructed manually a contour around each image, extended enough to cover all significant flux pixels (see Figure  2) and a large (at least 10 MUSE pixels around the image) empty area around it. Only pixels inside this region are considered for the χ 2 calculation. A weight image is associated for each object. For MUSE observations the weight values were estimated as 1/Var [p, q] in each pixel (p,q) from the NB variance image associated. For the HST images, the standard deviation σ was measured in an empty region close to the object, and the values of all pixels of the weight image were fixed to 1/σ 2 . The model includes a contribution from a local background (sky) estimated from the median flux measured in a large empty region close to each source. The χ 2 estimate in lenstool is then: with I p,q , M p,q and W p,q respectively the value in pixel [p,q] of the observed, model and weight images.

PSF estimation
Because we aim at modelling morphological properties of the LAEs, we have to obtain a very good knowledge of the point spread function (PSF) in both HST and MUSE observations. Since the PSF varies with wavelength across the MUSE spectral range we determined a specific circular Moffat monochromatic PSF for each object with an FWHM estimated following equation Examples of LLAMAS galaxies and spatial best models. In each row, we show, from left to right, the MUSE NB image, the best model and residuals. The contours present smoothed surface brightness levels at 12.5 and 25.0 × 10 −19 erg s −1 cm −2 arcsec −2 in light yellow and at 62.5, 100 and 200 × 10 −19 erg s −1 cm −2 arcsec −2 in red. The scale in the left panels are all at 2 arcseconds. We indicated from left to right the magnification value (µ), the redshift and ID of the source. In the last column, the areas without pixels indicates the edges of the area located in the region where the fit is applied, when any area without pixels is visible, it means that the region used is larger than the image presented in the figure. In the last row, the yellow line represent the critical line.
filter used for spatial fitting (F555W, F606W, F814W, F110W and F125W). We used at least 5 non-saturated, bright and isolated stars detected in each cluster (in all filters). The HST images of all these stars were combined to create a 51 × 51 pixels average image of the PSF centred on the brightest pixel used as PSF by lenstool.

Validation and selection with lenstool
To estimate the robustness of the source-plane modeling as a function of the S/N and extent of images we fitted a range of simulated Sersic profiles with the same method applied on real LAEs (for both MUSE and HST images). We generated more than 4000 simulated sources Sersic profiles with randomly selected parameters (position, scale radius, ellipticity, positional angle and magnitude). Each source was lensed by a real cluster model (4 randomly selected clusters from R21) and convolved with the MUSE PSF. We added random realizations of the noise based on different noise measurements from the Lyman-α NB images. Once a simulated image was created, we detected the multiple images using the Python/ photutils package (Bradley et al. 2016). If multiple images of the source were detected, only one image was fitted (chosen randomly). After applying the forward modelling approach described before, the best-fit parameters were compared with initial source parameters. Figure 3 shows how the difference between input and best-fit magnitudes varies as a function of S/N and area of each simulated image. We consider that a difference < 0.3 in magnitude is enough to get a good representation of the flux distribution in the source plane (with a relative error lower than 5% and 10% on the scale length and ellipticity parameter respectively). A region of the plot stands out visually, we defined a contour at the level of ∆ mag = 0.3 that we used after as a selection function for the source-plane spatial study. The S/N and the number of pixels contributing to the spatial fit were measured on the optimized NB image of each object with exactly the same detection process applied on simulated sources. The total distribution of the complete sample of LAEs is represented with grey points in Figure 3. Finally we obtained 475 MUSE Lyman-α images and 271 objects selected for source-plane emission spatial characterisation. Among them, 142 objects have enough resolved HST data to be characterised both in UV and Lyman-α emission (which represent 206 images). The green line on Figure 3, represents the criterion used to select sources for which morphology best-fit parameters are reliably recovered. To determine if the morphology of a source is correctly recovered by the fit procedure we compared in the source plane the input and best-fit models. Both input and best-fit profiles are elliptical profiles, to measure the difference between the 2 ellipses, we determined the proportion of non-recovered morphology with respect to the area of the input ellipse by measuring: R eps = (e 1 +e 2 −2 (e 1 ∩e 2 )) e 1 with e 1 the area of the input ellipse, e 2 the area of the best-fit ellipse and e 1 ∩ e 2 the common area. The two ellipses are defined by their scale radius, ellipticity and positional angle and the position. We consider the morphological properties measured well enough when this R eps < 0.3, e.g. the error on morphology of the source concerns less than 30% of the total source area. The signal-to-noise ratio used for this selection is integrated on the complete image. Images with a good signal-to-noise ratio (between 10 and 35) but a very large number of pixels (between 200 and 1000) can not be well fitted as the signal-to-noise ratio of each individual pixel is relatively low.

Modelling the morphology in the source plane
The large majority of previous studies on the SB spatial distribution of LAEs (Steidel et al. 2011;Wisotzki et al. 2016;Leclercq et al. 2017) mainly characterised the morphological properties of the Lyman-α nebulae through 1 or 2 components circular exponential models. Leclercq et al. (2017) decomposed the spatial fitting of the Lyman-α emission in 2 steps: first they fitted the UV emission on HST images using a single circular exponential Article number, page 7 of 22 A&A proofs: manuscript no. output  Fig. 3. Distribution of residual magnitudes between simulated source and best-fit parameters for MUSE Lyman-α NB images in the left panel and HST images in the right panel (see section 3.4). The grey points represent the LLAMAS galaxies. All the objects in the blue dashed contour were selected for spatial fitting with lenstool based on the magnitude difference between input and best-fit parameters (∆ mag < 0.3). The green dashed contour represents where the spatial fit will recover both the flux distribution and morphological information on the sources as estimated by the fraction of non-overlap between the simulated and fitted sources (see Sect. 3.4).
profile. Then they fitted the Lyman-α emission on MUSE NB images using a combination of 2 circular exponential profiles, at the same location, one with the scale radius fixed to the UV one. This modelling reproduced well the objects of these studies, with good residuals. But in our study, thanks to the lensing magnification, we observe LAEs with an improved spatial resolution. After applying this two circular exponential components model on all the Lyman-α images selected for spatial fitting, it was obvious that this model was not suitable for a large fraction of the LLAMAS galaxies. Indeed a lot of lensed galaxies from our sample present either a two-components Lyman-α distribution or a strongly asymmetric fainter Lyman-α emission surrounding a bright core emission. With the intention of obtaining the best source-plane reproduction for the UV and Lyman-α emission distribution we chose to apply between 11 and 7 different models, respectively presented in Table 2.
First if the object is detected in HST with enough S/N to be selected for the spatial fitting (see Figure 3), we performed 2 fits on the UV image: the first is a circular exponential profile (model M1) and the second an elliptical exponential profile (model M2). We applied the fit on the F555W, F606W, F814W, F110W or F125W filter (depending of target redshift, quality of the detection in terms of S/N of the images and available HST data in each cluster) whichever was the closest to the 1500-2500 Å restframe. Then we fitted 9 models on the MUSE Lyman-α emission: The complete description of free and fixed parameters used in each model is presented in Table 2.
To disentangle which modelling is the more adapted for each object we took into account the number of constraints, the number of free parameters and the final χ 2 , using a Bayesian Information Criterion (BIC) defined as: BIC = −2ln(L) − k ln(N) with L the likelihood of the best-fit, k the number of free parameters and N the number of constraints. We kept as the best modelling the fit with the minimum BIC.
Thanks to these different modellings, we obtained for each LLAMAS galaxy the best source-plane spatial model of both UV (when it is detected) and Lyman-α emission. Table 3 shows how sources are distributed among the best-fit models for each cluster and for the complete sample. We divided the different Lyman-α emission fits in 3 categories: two components fixed at the same location, two components free to vary, and the single component models. In the end, 67% of the LLAMAS galaxies are well described with a two-components fixed model (but the twocomponents circular exponential model is chosen only for 15 %). 12 % of the galaxies are best described by two free components, which correspond to LAEs presenting either 2 Lyman-α emission spatial peaks or a very asymmetric Lyman-α distribution.
We looked for trends between the best-fit model and the source S/N or lensing magnification effects on the models distribution. We measured median S/N of the objects in each category of models and found that the two more complex fits (models 7 and 8 with 2 free components) have a median S/N of twice as high as the median S/N of the other models (S/N=25). Finally we found that objects modelled with one of the two free components fits have a median lensing magnification of 6.2 versus 3.7 for the two fixed components models and 4.3 for the one component models. Both S/N and magnification seem to impact the mod-  Leclercq et al. (2017). The light pink histogram shows the distribution of the observed flux of Lyman-α images in the LLAMA sample. The dark pink histogram represents the intrinsic flux distribution of the LLAMA galaxies, obtained by dividing the observed flux by the lensing magnification. This histogram shows only the images selected for spatial fitting (see Figure 3). The typical uncertainty at 1 σ for LLAMAS flux measurement is about 1.2 × 10 −18 erg.s −1 .cm −2 . The grey histogram shows the observed Lyman-α flux of the Lyman-α haloes presented in Leclercq et al. (2017). els distribution between the different fits. The fact that high S/N and high magnification objects tend to prefer a more complex fit to reproduce intrinsic emission distribution suggests that LAEs have a more complex structure than a circular Lyman-α halo surrounding a circular ISM Lyman-α emission component.
Based on the spectral fit performed on the Lyman-α line (see Section. 3.1) we obtain a measure of the total Lyman-α emission flux for each image. Figure 4 shows the distribution of observed and intrinsic Lyman-α flux of the LLAMAS galaxies selected for spatial fitting. The grey histogram shows the distribution of the total Lyman-α line flux of the 145 LAEs presented in Leclercq et al. (2017). Thanks to the lensing magnification, our sample allows the characterisation of fainter LAEs than non-lensed studies. To estimate the more realistic intrinsic flux, we use the best-fit of each object to update the lensing magnification value to account for the morphology of each source. The magnification was estimated by measuring the ratio between Lyman-α emission in the image plane (measured in the best-fit image plane results) and the source-plane emission (measured in the best-fit source-plane model). These new magnification values better represent the total amplification of the Lyman-α emission compared to the previous values measured at a specific UV location, especially for the most strongly magnified sources. We find that for magnification under 10, the two magnification values are very similar; for higher magnification values, the new measurement is on average 2 to 10 times lower than the first estimate, which was an expected outcome. Indeed, for highly magnified images, the magnification varies across the image and thus a value measured at one position in the image does not reflect the average magnification of this image.

Results
In this section we present results on the Lyman-α nebulae morphology: spatial extent and axis ratio of the Lyman-α haloes, spatial offsets between UV continuum and Lyman-α emissions distributions.

Extended emission properties
Once we obtain a good model fit for each source, we need a common measurement to compare the spatial extent between the different individual objects. We use half-light (r 50 ) and 90%-light radii (r 90 ) to characterise both the UV and Lyman-α emissions. For the total sample, 38% of the LAEs have an HST counterpart bright enough to be spatially modelled with lenstool. We estimate r 50 and r 90 on a source-plane image produced from the parametric model (minimizing the BIC) of each source. Depending on the type of models, we measure the two radii from elliptical or circular rings. For objects with a circular best model (models M3, M5 and M9) we take ϵ = 0, for elliptical models M4, M6 and M11 we take the axis ratio value of the brightest component, and for the models M7, M8 (with 2 components at different locations) or, M10, we use the mean axis ratio measured from the model M10 (one elliptical exponential component). We randomly produce 200 source-plane images selected from the lenstool MCMC samples and the error bars are obtained by measuring r 50 and r 90 on each image. Error bars are estimated from the 68 % confidence interval on each side of the best model value. We also measure in the same way the r 50 and r 90 radius on the best-model images of the Lyman-α haloes from Leclercq et al. (2017).
Based on the r 50,Lyα and r 90,Lyα , we measure the concentration parameter of the Lyman-α emission, which measures how compact the Lyman-α light profile is. We measure values ranging from c Lyα = r 90 /r 50 = 33.3 to c Lyα = 1.15 with a median value of 2.57. The median value of c Lyα = 2.57 corresponds to a Sersic index of n = 1.2 which is close to the exponential profile index value n = 1. The concentration of Lyman-α emission is only weakly correlated with the Lyman-α extent, the more compact haloes are also the smaller.
We used the 90-light radius to compare the spatial extent of the Lyman-α emission in LLAMAS and UDF galaxies (Leclercq et al. 2017), as the half-light radius of UDF galaxies will be overdominated by the bright continuum-like component and will not reflect the extended halo properties. Figure 5 shows the distribution of the circularised Lyman-α 90%-radius r 90,Lyα as a function of UV r 90,UV (for multiply-imaged systems we present the value of the most extended image for both UV and Lyman-α emission). Among the LLAMAS, 40% of sources have a Lyman-α halo substantially smaller than the vast majority (i.e. 97%) of the halos in Leclercq et al. (2017). We measured the inverse-variance weighted mean of the ratio x 90 = r 90,Lyα /r 90,UV in both UDF and LLAMAS and find that the LLAMAS galaxies present on average a higher value, with x 90,µ,LLAMAS = 18.0, compared to the UDF sample for which we find x 90,µ,UDF = 10.40. The value of x 90,µ,UDF is consistent to the values measured in Leclercq et al. 2017 between the halo component scale radius and the core UVlike component. These values show that there is a large diversity of Lyman-α emission concentrations (very diffuse or peaked), but in 97% of the cases we observe a Lyman-α halo more extended than the UV continuum emission, and 75% and 47% are significantly more extended at 1 and 3 σ respectively. This confirms that the Lyman-α emission is intrinsically more extended than Article number, page 9 of 22 A&A proofs: manuscript no. output Model description N components Fixed parameters Free parameters M1. One circular exponential 1 ϵ = 0 ; n=1 x; y; a; m M2. One elliptical exponential 1 n=1 x; y; a; m; ϵ; θ M3. Two circular exponential 2 ϵ 1 = ϵ 2 = 0; n 1 = n 2 = 1; x 1 = x 2 ; y 1 = y 2 ; fixed components a 1 = a UV a 2 ;m 1 ; m 2 M4. Two elliptical exponential 2 ϵ 1 = ϵ UV ; θ 1 = θ UV x 1 = x 2 ; y 1 = y 2 ; a 2 fixed components a 1 = a UV ; n 1 = n 2 = 1; m 1 ; m 2 ; ϵ 2 ; θ 2 M5. Two circular exponential 2 ϵ 1 = ϵ 2 = 0; n 1 = n 2 = 1 x 1 = x 2 ; y 1 = y 2 ; fixed components a 1 = a UV ± 0.1"; a 2 ;m 1 ; m 2 ; M6. Two elliptical exponential 2 n 1 = n 2 = 1 x 1 = x 2 ; y 1 = y 2 ; fixed components a 1 = a UV ± 0.1"; a 2 ;m 1 ; m 2 ; ϵ 1 ; ϵ 2 ; θ 1 ; θ 2 M7. Two free circular 2 n 1 = n 2 = 1; ϵ 1 = ϵ 2 = 0 x 1 ; x 2 ; y 1 ; y 2 ; a 1 ; a 2 ; m 1 ; m 2 components M8. Two free elliptical 2 n 1 = n 2 = 1 x 1 ; x 2 ; y 1 ; y 2 ; a 1 ; a 2 ; m 1 ; m 2 components ϵ 1 ; ϵ 2 ; θ 1 : θ 2 M9. One circular exponential 1 n 1 = 1; ϵ 1 = 0 x 1 ; y 1 ; a 1 ; m 1 component M10. One elliptical exponential 1 n 1 = 1 x 1 ; y 1 ; a 1 ; m 1 , ϵ 1 , θ 1 component M11. One Sersic circular exponential 1 x 1 ; y 1 ; a 1 ; m 1 ; n 1 component Table 2. Description of free and fixed parameters used for the different spatial models applied. x and y are the RA and DEC positions, m is the magnitude, a is the scale radius, n is the Sersic index, ϵ is the ellipticity and θ is the positional angle. The models M1 and M2 are applied only on the UV images and models from M3 to M11 on Lyman-α images.  Table 3. Summary of the spatial best-fits distribution in each cluster. The first column shows the total number of Lyman-α images selected for spatial fitting and the boldface numbers the number of unique objects. The second and third columns shows the number of images with both UV and Lyman-α detection and only a Lyman-α detection respectively. The three last columns show the distribution in the 3 main categories of Lyman-α emission models. The last row shows the repartition for the global sample.
the UV component as measured in previous studies (Steidel et al. 2010;Leclercq et al. 2017;Wisotzki et al. 2018). We find smaller values for x 50 = r 50,Lyα /r 50,UV with x 50,µ,LLAMAS = 12.0 and x 50,µ,UDF = 4.6. From the spatial modelling, we obtain a measure of the axis ratio of the UV and Lyman-α emission distribution (through the axis ratio q = b/a between the major (a) and minor (b) axes of the best-fit ellipse). This value is interesting to study as it should correlate with the spatial distribution of the neutral hydrogen in the CGM and possibly with the inclination of the galaxy. We use the best model of each source to determine the degree of ellipticity of the UV and Lyman-α distributions in the source plane. Among the galaxies selected to be spatially modelled in UV emission, 38 % prefer a circular best model and 62 % an elliptical one. For Lyman-α emission, 52 % of the haloes are better described by a circular model and 48% an elliptical one (we consider models M7 and M8 as elliptical models). The median UV axis ratio q, measured only on elliptical sources, is ⟨q⟩ ∼ 0.22. For multiple systems we measure a magnification weighted mean of the axis ratio values of each image of the system (e.g. the same results are found if we use signal-to-noise ratio weighted mean as the magnification value is strongly correlated with the spatial integrated signal-to-noise ratio). When we consider only images located inside the green contour in the second panel of Figure. 3, we find a median of ⟨q⟩ = 0.38. Applying the same procedure to measure the distribution of Lyman-α emission axis ratios, we find that Lyman-α haloes are on average less elongated (median value of 0.22 for UV and 0.48 for Lyman-α emission). The distributions of axis ratio values for the UV and Lyman-α emission distributions are presented in the Figure 6. These results are consistent with the trend presented in  and with the measurements for the high-z simulated LARS galaxies in Guaita et al. (2015) who found similar values of UV and Lyman-α axis ratio (with axis ratio of their Lyman-α emission between 0.4 and 0.9, with a mean value around 0.7). Wisotzki et al. (2016) measured also that 75% of their LAEs came out with a UV axis ratio smaller that 0.5. We use the axis ratio from the best source-plane model and apply the same procedure described previously for the multiple systems. The axis ratio is a useful indicator of the galactic disk inclination and hydrogen distribution morphology for the Lyman-α emission.
Haloes with a small axis ratio value indicate that the CGM is structured along a preferred direction around the UV source.
We find no significant variation of the axis ratio with redshift. However we found a significant correlation between the spatial extent of the emission in the image plane and the proportion of circular and elliptical best model. We measured the division between circular and elliptical in three equal-size bins based on number of pixels in the detection map. For Lyman-α emission, we measure 13%, 40% and 77% of elliptical best model for respectively 3.3 < area < 7, 7 < area < 9.6 and 9.6 < area < 37.5 with area in arcsec 2 . The more a Lyman-α halo is resolved (with a high detection map area which is strongly correlated to the lensing magnification and signal to noise ratio), the more it would prefer an elliptical model. The same effect is observed in S/N as detection map area and S/N are strongly correlated (see Figure 3). These results also indicate that the circular shape measured on the Lyman-α emission could be partly a limitation due to lower signal to noise ratios and incorrect PSF estimation.
Since the lensing models are used to measure all the properties presented here, lensing uncertainties in the mass model used in the lensing reconstruction could strongly impact the results. To estimate the impact of the lensing model on the spatial measurements we study the dispersion of the different measurements between the different images of the 80 multiple systems with at least two images selected for spatial fitting (including 22 systems with 3 or more images with both Lyman-α and UV detections). For each system we measured the magnification-weighted mean (⟨.⟩ µ ) and standard deviation (σ µ ) for each of the following parameters : r 50,Lyα , r 50,UV , UV (q UV ) and Lyman-α axis ratio (q Lyα ) and spatial offsets (presented in the following sections). All these measurements are presented in Figure 7.
For the spatial extent measurements (left panel of the Figure 7), we found that 19% of the multiple systems present a small dispersion between the different images (variation smaller than 20% of the mean value in ⟨r Lyα ⟩ µ and in ⟨r UV ⟩ µ ). For 30% of the systems, the variation between the different images is moderate (between 20 and 50% of the mean value) and the remaining 50% present a large variation between multiple images. Finally 4 multiples systems present a variation of the Lyman-α extent larger than the mean value. After a visual inspection, we found that these systems are the most magnified galaxies (with a total magnification between 20 and 50). Within these specific systems the variations in amplification and shear, and therefore in spatial resolution, between images are very large, which explains the variations observed in the measurements. In all the results presented in this study, we keep the values measured on the most extended image of each multiple system. We find the same kind of trends for the axis ratio (middle panel of the Figure 7).

Variations in the spatial extent of Lyman-α emission
In the Figure 5, we observe a correlation between Lyman-α and UV continuum r 90 for both UDF and LLAMAS datasets. The error-bars weighted Pearson coefficients (ρ UDF = 0.22 and ρ LLAMAS = 0.20) suggest that these correlations are weak but significant, with p-values of (p 0,UDF = 0.02 and p 0,LLAMAS = 0.05). The LLAMAS values seem more scattered (σ UDF = 7.6 kpc < σ LLAMAS = 22 kpc) and the correlation more marginal. This higher spread could be due to the larger uncertainties of the LLAMAS measurements only or also to the fact that lensing studies provide access to a new population of weaker and smaller LAEs than those characterised by the previous studies. We notice a strong effect of the detection isophotal area (number of pixels in the detection map of each Lyman-α image) on this correlation. We divided the sample in three equal-size bins of low, medium and high value of spatial extent and measured the weighted Pearson p-value and coefficient in each bin. We find that the less extended images (between 3.3 and 6.1 arcsec 2 ) present no correlations between Lyman-α and UV r 90 (p 0 > 0.2 and ρ < 0.1). On the contrary the two bins of medium (between 6.1 and 10.2 arcsec 2 ) and highly (between 10.2 and 44.2 arcsec 2 ) extended images show a significant correlation (p 0 = 0.033, ρ = 0.36 and p 0 = 0.05, ρ = 0.33 respectively). This shows that the correlation is stronger when we look only at the images with higher resolution and then suggests that the absence of correlation for the low spatial resolution is due to some bias. The same effect is measured with the magnification and the signal to noise ratio which are both strongly correlated with the image-plane resolution. Finally the different models used to reproduce the LLAMAS galaxies, of which 7 out of 9 are completely independent of UV spatial properties could probably play a role in measuring a lower correlation between Lyman-α and UV spatial extent.
We searched for correlations between the Lyman-α emission size (using the r 50,Lyα parameter) and the physical parameters governing the host galaxies (star formation rate, Lyman-α equivalent width, flux and luminosity, and properties of the Lyman-α line profile), and find three more or less significant correlations presented in Figure 8. First we find a very significant positive correlation between r 50,Lyα and the FWHM of the Lyman-α line (with weighted Pearson coefficient of ρ = 0.35 and p-value of p 0 = 0.0002). Secondly, we observe a negative trend between r 50,Lyα and the Lyman-α rest-frame equivalent width W 0 (with weighted Pearson coefficient of ρ = −0.25 and p-value of p 0 = 0.07; the more extended the Lyman-α halo, the smaller is W 0 . Finally we observe a weaker positive correlation between r 50,Lyα and the UV SFR (ρ = 0.18 and p 0 = 0.0002). We observe similar correlations between theses three parameters and r 90,Lyα .
Article number, page 11 of 22 A&A proofs: manuscript no. output

Spatial offsets between UV and Lyman-α emissions
Looking at the best model of each UV and Lyman-α source, we notice a significant spatial offset between the two peaks locations; some examples are shown in Fig 9. We measure this offset in each galaxy with a UV component observed in HST and selected for spatial fitting (see Section 3.4). We performed two different measurements on each galaxy. First we measured the two locations of peak emission (UV and Lyman-α) in the image plane. We measure the peak location after a spatial Gaussian filtering of FWHM=0.2 or 0.1" for Lyman-α NB and HST images, respectively. We ray-trace these values into the source plane to measure the intrinsic (physical) spatial offset between the two emission peaks, in kpc. The distribution of this "peak to peak" offsets measurements is presented in dark blue in Figure 10. For the total sample, we measure a median value of ∆(peaks) = 0.67 ± 0.13 kpc. When we keep only spatial offsets larger than 2 MUSE pixels (i.e. 0.4") in the image plane, as they are considered more significant, the median value became ∆(peaks) = 1.16 kpc. Second, we measure the spatial offset between UV and Lyman-α best model centroids (measured on the best-model source-plane image). For models with 2 Lyman-α emission peaks (M7 and M8 in Table 2), we measure the offset between the UV component centroid and the 2 Lyman-α peaks and kept only the smallest. The distribution of these offset measurements is shown in light blue in Figure 10, with the median value of the total sample is ∆(centroids) = 0.58 ± 0.14 kpc. Among the total sample, we identified 3 galaxies (two of them presented in Fig 9) which present a strong Lyman-α absorption feature (LBG) in their spectrum and a huge spatial offset between UV and Lyman-α component.
For these three sources we can attribute the large offset value to the absorption of Lyman-α at the location of the UV component. We compare our results with the three recent studies on spatial offsets performed by Hoag et al. (2019); Lemaux et al. (2020) and Ribeiro et al. (2020) on three samples with similar ranges of Lyman-α luminosities and UV magnitudes, and found a very good agreement. Hoag et al. (2019) measured the spatial offsets in 300 galaxies at 3 < z < 5.5 observed in slit data, and find a median value of 0.61 ± 0.05 kpc. Ribeiro et al. (2020) measured a similar value in a sample of ∼ 900 galaxies at 2 < z < 6 of 0.60 ± 0.05 kpc. When they selected 11% of galaxies with secure offsets (after a visual inspection), the median value become 1.9 ± 0.2 kpc. Finally Lemaux et al. (2020) measured for 64 objects (mix of lensed and non-lensed galaxies) with 5 < z < 7 a median offset of 0.61 ± 0.05 kpc. In LLAMAS galaxies, we do not measure any significant variation of the offset values with redshift. Finally, we measure the impact of the lensing model to the offset measurements (see Figure 7). We found that 75% of the multiple systems have a dispersion smaller than the mean value.

Significance of the offsets measurements
To measure the robustness of these measurements we estimate the probability of measuring such offsets if Lyman-α and UV peaks are supperposed. Using the MCMC optimisation result of the best model from lenstool, which provides a list of values of x and y positions in the source plane tested during the optimisation, for both UV and Lyman-α models. We apply the measured offset to the Lyman-α position sample to center the cloud of Lyman-α positions on the UV ones. We randomly draw 10.000 pairs of UV and Lyman-α positions and measured for each couple the offset value δ 95 , measured in this way, corresponding to a cumulative probability to randomly measure an offset smaller than the real one ∆ of 95%. In other words it corresponds to the offset value from which there is less than 5% of chance to randomly measure a similar or larger offset from the best models. We consider that an offset measurement is significant when δ 95 < ∆. We find that 63% of galaxies have a significant offset. Among the remaining 37%, we measure with the same method, if the offsets corresponding to the cumulative probability of 68% ( Best-fit oAEsets Peak-to-peak oAEsets 7. Dispersion measured in each multiple system of the three main parameters measured : UV and Lyman-α extent, axis ratio and spatial offset. We measured in each multiple systems the µ-weighted mean values and µ-weighted dispersion of values among the images. From left to right we represent the dispersion of half-light radius (r 50 ) with respect to µ-weighted mean values, the dispersion of axis ratio (q) with respect to µ-weighted mean values and the dispersion of spatial offset (∆ Lyα−UV ) with respect to µ-weighted mean values. In each panel the grey stars and squares represent 2 images systems and colored points the 3 and more images systems. In the left and middle panels, the stars represents UV r 50 and q respectively, and square Lyman-α values. In the right panel, the stars represent best model centroid offsets and green points represents the "peak to peak" offset values (see Sect. 4.3). In all panels we plot dashed lines to represent the 10%, 50% and 100% error levels.  at one sigma) which are smaller than ∆ (corresponding to offset significant at only one sigma). We consider that all the galaxies with ∆ < δ 68 are compatible with the non-offset scenario. For these galaxies, we use the δ 68 as an upper limit to the offset measurement. We measure this fraction of spatial offsets in 3 bins with 0 < ∆ < 1 kpc, 1 < ∆ < 2 kpc and ∆ > 2 kpc. We find respectively 48%, 92% and 92% of significant offsets. The fraction of measurements compatible with the non-offset hypothesis is of 37%, 5% and 2% in the 3 same bins, respectively. Note that non significant offsets could simply be due to small intrinsic values or indeed coincident UV and Lyman-α emission.
The spatial offsets measured in kpc in the source plane should be correlated to the UV size of the galaxy to be physically interpreted. With the aim to propose a more consistent measurement of spatial offset with respect to the UV continuum emission size, we measure the elliptical distance, ∆ ell (normalised by the UV emission 90% isocontour elliptical parameters) between the Lyman-α emission centroids and the UV emission centroid: with (X Lyα , Y Lyα ) the position of the Lyman-α emission centroid in the referential formed by the axis of the UV emission elliptical distribution and (R x,UV , R y,UV ) respectively the semi-major and semi-minor axis of the UV elliptical emission distribution. A value of ∆ ell < 1 means that the Lyman-α emission center is located in the UV continuum component and probably produced by substructures within the star-forming region. When ∆ ell > 1, the Lyman-α emission peak is produced outside of the stellar component, probably by satellite galaxies, or large-scale effects in the CGM. We show the distribution of ∆ ell values in grey in Fig 11. We measure that 40% of the galaxies present an Article number, page 13 of 22 external spatial offset as high as ∆ ell >2.
We measure no correlation (p-values > 0.05) between the values of elliptical distance and the physical parameters of the host galaxies (Lyα W 0 , flux and luminosity, ∆ in kpc). After dividing their sample of galaxies in 2 bins, based of W 0 , Hoag et al. (2019) found a significant trend: galaxies with lower W 0 values show higher offset values (mean of 1.92 ± 0.13 kpc) than galaxies with higher W 0 (mean of 1.51 ± 0.11 kpc). We do not confirm this trend, neither with offset value in kpc nor with the elliptical distance. Lemaux et al. (2020) observed that the spatial offset increases with the UV brightness of the galaxies. We find the opposite correlation when we compare UV SFR values and elliptical distances for all galaxies with secure spatial offset measurements (90 sources). Galaxies with an elliptical distance ∆ ell < 3.9 (45 sources) show on average a higher UV SFR value (⟨S FR⟩ = 1.74 ± 0.23 M ⊙ yr −1 ) than galaxies with ∆ ell > 3.9 (45 sources, ⟨S FR⟩ = 1.26 ± 0.28 M ⊙ yr −1 ). We notice also that the galaxies with a higher elliptical distance present a smaller UV size (⟨r 50,UV ⟩ = 0.08 ± 0.01 kpc) and larger Lyman-α versus UV emission extent (⟨ r 50,Lyα r 50,UV ⟩ = 25.5 ± 4.3), than galaxies with smaller elliptical distances (⟨r 50,UV ⟩ = 0.31 ± 0.4 kpc and ⟨ r 50,Lyα r 50,UV ⟩ = 7.0 ± 0.8). Thus, more than half of the Lyman-α peaks are located outside of the stellar body of the source which could be due to the presence of an extra Lyman-α emission source such as a satellite or merging galaxy.
Besides, it should be easier to distinguish a spatial offset in the Lyman-α emission produced by a satellite galaxy for smaller (low r UV ) galaxies, because these sources should emit less Lyman-α photons and thus the contribution of a satellite galaxy would be easily detectable in the global profile. This could explain why the galaxies for which we measured higher elliptical distances present on average a smaller UV component.

Production of mock observations
In order to physically interpret these results we compare the measured LAE properties with a cosmological radiation hydrodynamical (RHD) simulation of a high-redshift galaxy evolving from z = 6 to z = 3, described in Mauerhofer et al. (2021). This zoom-in simulation was produced using the ramses-rt code (Rosdahl et al. 2013;Rosdahl & Teyssier 2015). The simulation includes all the expected Lyman-α production mechanisms (photo-ionization and photo-heating of hydrogen by local sources, collisional excitation of hydrogen as well as contribution from the UV background) and so represents a powerful tool to study Lyman-α photons escape in both ISM and CGM. This simulated galaxy has been deliberately chosen to be representative of the faint-UV LAEs detected in the recent MUSE studies of LAEs at high redshift (Leclercq et al. 2017;Wisotzki et al. 2016) in terms of halo mass (M h = 6 × 10 10 M ⊙ at z = 3). We study mock observations of the Lyman-α line cubes and UV continuum emission at 1500 Å rest-frame. Each mock Lyman-α dataset consists in a 10"×10"×10.9 Å datacube centered on the Lyman-α line with 0.067"×0.067"×0.0625 Å pixels, which is three times better than the MUSE sampling both in spatial and spectral directions. The UV continuum 1500 Å rest-frame images are produced with 0.01"×0.01" pixels which are 5 times smaller than HST/ACS pixels. To represent the diversity of observed Lyman-α profiles, 12 mocks datacubes were produced for each of the 129 simulation timesteps, by projecting along 12 different line of sights defined by healpix Nside=1 (identical at all redshifts). Thus 12×129 mocks were produced at the 129 different redshifts ranging from z = 3.000 to z = 5.989 with a regular lookback time interval of 10 Myrs. Although this simulation focuses only on one galaxy, the fact that this galaxy is studied at 129 different redshifts in 12 different directions adds some diversity due to variations in SFR with time, galaxy growth, effects of radiative transfer into the CGM and line of sight projections. The global properties of the simulated galaxy sample and UDF and LLAMAS sample are presented in Figure 12. Altogether, this sample of mock observations has physical properties close to the LLAMAS properties in terms of redshift, Lyman-α luminosity and SFR. The median redshift value is z = 3.92 for the simulated data versus z = 4 for the LLAMAS sample. Note that the distribution of UV magnitudes in the simulation is narrower than that of LLAMAS galaxies, but it roughly covers a similar range of M UV with a median value slightly brighter than the median value of the LLAMAS sample (−17.8 vs. −17.0 respectively).
The SFR of the simulated galaxy varies with redshift, on average the SFR increases during the formation and evolution of the galaxy from 0.5 M ⊙ .yr −1 to 1.23 M ⊙ .yr −1 between z = 6 and z = 3 with few SFR peak episodes (up to 3.0 M ⊙ .yr −1 ) during its history. These SFR variations are strongly correlated with the variation of the total Lyman-α luminosity. The median value of the SFR on the simulated sample (SFR Sim = 0.80 M ⊙ .yr −1 ) is larger than the LLAMAS median value (SFR LLAMAS = 0.48 M ⊙ .yr −1 ), although lower and higher SFRs can also be found in the LLAMA sample. We notice a very good match in terms of Lyman-α luminosities between LLAMAS and the simulation, with Lyman-α luminosities ranging from log(Lyα/erg.s −1 ) ∼ 40 to 42.5. with a median value around 41.5.
We note that the UV size of the simulated galaxy increases from 0.15 to 0.45 kpc between z = 6 and z = 3 as a result of the continuous mass growth due to gas accretion and mergers over this period of ∼ 2 Gyr.
To compare with lensed and non-lensed MUSE observations of high redshift LAEs, we produced mock "UDF-like" and "LLAMAS-like" Lyman-α NB and HST images from the simulated galaxies. To produce Lyman-α NB images we collapse a cube containing all the Lyman-α emission without continuum. To reproduce non-lensed MUSE Lyman-α NB images of LAEs, we convolve the initial raw Lyman-α NB images by a typical MUSE UDF PSF, depending on the redshift of the galaxy, as described in Section 3.3. Finally we re-grid the images at the MUSE sampling of 0.2"×0.2" pixels and then add a random Gaussian noise based on the typical level observed in UDF Lyman-α NB images. We follow the same method for UV images; we convolve the raw UV images with UDF HST PSF, constructed following the method explained in Section 3.3, using HST F606W or F814W images (depending on the redshift). We finally re-grid the images at the HST sampling of 0.05"×0.05". We add a random Gaussian noise based on the typical level observed in the HST images of the UDF. To reproduce "LLAMAS-like" observations, we first chose randomly a cluster model and a specific source location from the R21 data release. We lens both UV and Lyman-α NB images by the lens model. We then follow the same procedure for PSF convolution, pixels sampling and addition of Gaussian noise as for UDF-like observations. We use typical MUSE PSF parameters and HST PSF measured in the cluster used as a lens. For both kinds of observations (UDF-like and LLAMAS-like) we use the python package photutils to detect the different images, following the same criteria used for the LLAMAS sample (see section 2.5). When the lensing produces multiple images, we keep only one of them for comparison with observations. We The light blue distribution represents spatial offsets measured between UV and Lyman-α emission centroids measured in the source plane, and the dark blue one the distribution of spatial offsets measured between UV and Lyman-α emission peaks projected in the source plane. The hatched distribution represents the galaxies with a spatial offset probability higher than 95%. The dashed and solid lines show the median values of centroids and peaks offsets distributions respectively. apply a S/N threshold of S/N=6 for Lyman-α NB "UDF-like" images as applied in Leclercq et al. (2017). For LLAMAS-like images, we apply the same selection for spatial fits as for real observations, highlighted by the blue contour in Figure 3. We finally obtain a sample of 1164 raw images (both UV and Lyman-α NB), 271 "UDF-like" images and 254 "LLAMAS-like" (the other produced UDF-like and LLAMAS-like images had too faint UV or Lyman-α images to be detected or spatially characterised).

Variations in Lyman-α extent
We apply the same spatial fit as presented in Section 3.5 to the three types of mock observations ('raw' simulation, UDF-like and LLAMAS-like observations). For UV measurements, we apply a single elliptical exponential component fit on the three different mock observations. For the raw data and UDF-like observations we apply on Lyman-α NB images a single elliptical exponential two components fit (M4 in Table 2). For LLAMAS-like mock observations, in order to fairly compare with the LLAMAS sample results, we applied 3 different fits on each selected image and compared BIC criteria to choose the best model of each galaxy. We choose to apply the models labelled M4, M6 and M8 in Table 2. We measured half-light and 90%-light radius on best-fit images for each source. The morphology of simulated LLAMAS haloes gives a similar range of concentrations (c = 1.16-32.5) as the real observations. Figure 13 shows the distribution of Lyman-α r 90 for UDF, LLAMAS, UDF-like and LLAMAS-like galaxies. We observe that the simulated UDF galaxies (empty grey points) are confined to a small region of the cloud of points from Leclercq et al. (2017). The All P o↵ > 95% Δ ell < 2 : 38 % Δ ell > 2 : 62 % Δ ell < 1 : 15 % Median value : 3.7 Fig. 11. Distribution of elliptical distance (cf Section 4.4) measured in the source plane between Lyman-α emission centroids and the ellipse formed by the UV emission distribution (using r 90,UV as radius). The grey distribution represents all the LLAMAS galaxies and the black hatched only the galaxies with a spatial offset probability higher than 95%. The vertical dotted line represents the separation between internal and external spatial offsets (∆ ell = 2). simulated LLAMAS-like galaxies are located in a different region, presenting both smaller and higher values of Lyman-α extent, closer to the real LLAMAS galaxy Lyman-α spatial extent values. In all LLAMAS-like and UDF-like galaxies, the Lyman-α emission is measured to be more extended than the UV central component, as for real LLAMAS galaxies. However the mean x = r 90,Lyα /r 90,UV ratio is smaller for the simulated galaxies whatever the simulated observation method (⟨x 90,UDF−sim ⟩ = 4.3 and ⟨x 90,LLAMAS−Sim ⟩ = 5.0). On average the simulated LAEs present a Lyman-α emission which is less extended, with respect to the UV component, than the observed ones. The fact that the same simulated galaxy leads to different physical parameters is indicative of an inconsistency between the two measurements which could be due as much to the measurement method as to the instrumental effects (PSF and noise). Figure 14 represents the relative errors (i.e. the ratio of LLAMAS-like or UDF-like measurements and original mock measurement) of UV and Lyman-α size measurements both for LLAMAS-like (in blue) and UDF-like sources (in green). We notice that UDF measurement tend to overestimate both UV and Lyman-α emission sizes. This effect could be due to the PSF smoothing, dominant in the UDF simulated images. The UV measurements on the LLAMAS-like simulated galaxies are in good agreement with the values estimated on the original high resolution images of the simulation. The Lyman-α extent is less constrained, however no systematic bias is observed. These different results confirm the gain from the lensed samples in the study of extended Lyman-α emission. LLAMAS galaxies represented here are only the objects selected for spatial fitting (see Sect. 2.5). All values are measured following the same procedure on UV and NB Lyman-α images.

Origin of the extended Lyman-α emission
We notice in the simulations a strong effect of the line of sight on the spatial extent measurements but we did not measure any significant correlation between the Lyman-α spatial extent and the physical parameters of the galaxies in the simulation both for the UDF-like and LLAMAS-like samples. Both in the UDF and LLAMA samples, the Lyman-α emission is almost always more extended than the UV component traced by the UV restframe emission thus dominated by the young stellar population emission. From the results of this study we propose two possible scenario to explain this result. First the Lyman-α haloes could be due only to the scattering of the Lyman-α photons from the source emission to the outskirts of the halo. This scenario is supported by strong correlation measured between r 50,Lyα and the width of the Lyman-α line and presented in the Figure 8. Indeed, assuming that more extended Lyman-α haloes trace optically thick media where the number of Lyman-α scatterings is increased, we expect from theoretical studies that these halos will exhibit broader line profiles as a result of resonant scattering (e.g. Verhamme et al. 2018) Secondly, the correlation measured between r 50,Lyα and the UV SFR (Figure 8) indicate that the spatial extent of the CGM may also depend on the UV stellar activity. Anisotropic outflows, as observed in star-forming galaxies at 2 < z < 6 (Steidel et al. 2010;Pelliccia et al. 2020;Ginolfi et al. 2020) can be produced by stellar feedbacks. These outflows could push the gas and thus let Lyman-α photons diffuse further away from the galaxy center by decreasing the covering fraction of the hydrogen gas (Lemaux et al. 2020); causing the halo expansion. This scenario could also explain the asymmetric and anisotropic Lyman-α emission distribution noticed in some LLAMAS galaxies. Rasekh et al. (2021) measured the same trend at low redshift (in 28 galaxies out of a sample of 45 galaxies at z < 0.24) where the Lyman-α emission extent is correlated with the stellar mass and the star formation regions sizes.

Offsets in mocks vs. observations
We measured a significant spatial offset between UV and Lymanα emission in 60% of the LLAMAS galaxies, ranging from 0.1 to 7 kpc. Following the same procedure used in observational data (Sect. 3.2 ), we measured spatial offsets between UV and Lyman-α emission in the three samples of simulated sources. The left panel of Figure 15 shows the values of spatial offsets measured in raw simulations (grey) and LLAMAS-like sources  Fig. 13. Lyman-α emission 90-light radius r 90,Lyα as a function of the UV emission 90-light radius r 90,UV for UDF galaxies (in grey), LLAMAS galaxies (in red), UDF simulated galaxies (empty grey circles) and LLAMAS simulated galaxies (empty red circles).
(yellow). The values shown in orange are the sources with a spatial offset larger than 0.4" in the image plane (which is often the limit given in the literature to claim significant offsets in MUSE observations). We notice that a high number of LLAMASlike sources (22%) present an offset larger than 0.4" in the image plane, thanks to lensing magnification. The distribution of spatial offset measured on LLAMAS-like simulated galaxies is close to the raw simulations measurements. (∆ UV−Lyα,simus = 0.40 kpc and ∆ UV−Lyα,LLAMAS−like = 0.64 kpc). In the UDF-like mocks, only 16 sources (6%) show an offset larger than 0.4", which can explain why no spatial offsets are usually reported with MUSE in nonlensed galaxies due to resolution limits and emphasizes the gain provided by lensing surveys. When we compare the distribution of spatial offsets measured in the LLAMAS sample with the distribution measured in the LLAMAS-like sample (right panel of Figure 15), we find a very good match, highlighting that these simulated galaxies incorporate physical mechanisms capable of producing similar spatial offsets to those observed. Figure 11 Article number, page 17 of 22 represents the values of ∆ ell for the simulated LLAMAS-like galaxies in gold. We notice that the simulated galaxies have on average lower elliptical distances but span anyway a large range of r ell values as measured in the LLAMAS galaxies. 72% of the simulated LLAMAS-like galaxies present an internal spatial offset with ∆ ell < 2. We do not measure very high values (>20) of ∆ ell in the simulated galaxies. In the LLAMAS, we identify the galaxies with ∆ ell > 20 as being LBGs with a strong absorption feature observed in their spectra (one example is presented in the Figure 10 at z = 4.69).

Origin of the spatial offsets in the simulation
Thanks to the high spatial resolution of the simulation, we can investigate the origin of the spatial offsets found in the mocks. Figure 16 shows that there are two regimes in the elliptical distance distribution of LLAMAS galaxies: sources with ∆ ell < 2 and ∆ ell > 2. Galaxies with ∆ ell < 2 represent 38% of the LLAMAS sample and 72% of the LLAMAS-like simulated galaxies. In this case, the spatial offset is likely due to an offsetted star formation clump emitting a high quantity of Lyman-α photons or due to the inhomogeneous neutral hydrogen distribution surrounding the galaxy. For low redshift galaxies, the LARS sample (Hayes et al. 2013Östlin et al. 2014) observed a high distinct clumpiness of the ISM emission in both UV and Hα emission as well as a complex structure of the Lyman-α emission. As we know that such sub-structures are present in high redshift galaxies (Elmegreen et al. 2013 andFörster Schreiber et al. 2018), they could explain the formation of small offsets between Lyman-α and UV emission at the scale of the continuum component. In the high resolution images of the simulated galaxy, we can observe a very clumpy UV emission (see two examples in Figure 17) and we are able to visually associate some small spatial offsets values with a clear UV emission clump in the outskirt of the galaxy (one example shown in the top row of Figure 17). Due to the resolution limits, it is hard to distinguish the different potential UV emission clumps in the real observed LLAMAS galaxies, except for some highly magnified (less than 10 in the LLAMAS) objects as for example the source at z = 4.03 in Figure 9.
Moreover we measured a significant (> 2 HST pixels) distance between the UV brighter pixel (i.e. the UV peak emission location) and the UV emission centroids in 27% of the LLAMAS sources (distribution shown in grey in Figure 16) which reveals the clumpy nature of some galaxies and could explain the formation of some small offsets. Among these galaxies, 54% present a spatial offset with ∆ ell < 1, the presence of a clumpy UV emission distribution seems to favour the measurement of an internal offset in the galaxy. We notice the same trend in the simulated LLAMAS-like galaxies: 22% of the objects present a significant UV-UV offset and among them 90% present a value of ∆ ell < 1.
On the other hand, 62 % of the LLAMAS galaxies have an elliptical distance too high to be explained by internal substructures of the UV emission. Many other scenarios are compatible with these large offsets such as inflows of gas, emission from faint satellite galaxies, outflows in the CGM or Lyman-α scattering effects. In the LLAMAS galaxies, without new observations of ISM lines (such as Hα or [O iii]) or deepest UV data, it is difficult to disentangle these scenarios in each individual galaxy. The possibility of bright Lyman-α emission from faint satellite or merger galaxies was suggested by Maiolino et al. (2015) and Mitchell et al. (2021) and could explain both the larger values of r 50,Lyα /r 50,UV and elliptical distances. Nevertheless, we did not clearly detect UV satellites coincident with this Lyman-α emission, even in the deepest HST fields. In the original simulation images, we can visually assign the spatial peak of the Lymanα to a faint UV emission component located outside the main UV component, as shown in Figure 17, in 24% of the LLAMAS simulated galaxies which represent the larger spatial offsets. The small measured offsets (< 10 kpc) suggest mainly cases of merger galaxies emission both in LLAMA and LLAMAS-like samples. We do not notice any significant trend between the spatial offset values and the physical properties of the LLAMAS-like simulated galaxies.

Other possible origins and future direction
By comparing lensed LAE observations and a zoom-in simulation, we identified two different ranges of ∆ values (cf. Figure 11) showing that different scenarios are at play in the formation of spatial offsets between UV and Lyman-α emission at z > 3: 1. Formation of small-scale offsets by bright UV clumps inside or in the outskirts of the main UV component. If most of the Lyman-α photons are produced in the ISM, the initial Lyman-α emission distribution should follow the ISM spatial morphology. In their simulated galaxy, Mitchell et al. (2021) showed that at r < 7 kpc the ISM contribution dominate the Lyman-α emission. This scenario is coherent with LARS galaxies (Östlin et al. 2014) at low redshift which are highly irregular and clumpy. The LARS study shows also that the Lyman-α W 0 is varying significantly as a function of the position in the galaxy. Recently, Rasekh et al. (2021) have measured the same type of spatial offsets between UV and Lyman-α emission centroids in low redshift galaxies (with offsets ranging from 0.8 to 2.25 kpc and a median value of . Left: Distribution of the spatial offsets measured in the simulated raw data (grey), simulated lensed observations (yellow). The black distribution show the spatial offsets which are observed in the image plane higher than 2 MUSE pixels. Right: Normalized spatial offset distribution for the LLAMAS galaxies (in blue) and the LLAMAS simulated galaxies (in yellow). The dashed lines show the mean value of the two samples. . Left: Distribution of elliptical distance measured in the source plane between the Lyman-α emission centroids and the ellipse formed by the UV emission distribution (using r 90,UV as radius) for the LLAMAS galaxies (in grey) and the simulated LLAMAS galaxies (in yellow). Right: Distribution of the UV-UV spatial offsets measured between the UV emission peak and the UV centroid position in LLAMAS galaxies (in grey) and simulated LLAMAS galaxies in yellow. We show here only the galaxies for which this UV-UV offset is higher than 2 HST pixels in the image plane.
1.13 kpc) and found a correlation between the offset and the stellar mass and the size of the star-formation regions. Guaita et al. (2015) showed that the UV component of these galaxies presents a morphology compatible with mergers or star-burst galaxies at high redshifts. Messa et al. (2019) measured that the LARS galaxies with a Lyman-α escape fraction higher than 10% have more than 50% of their UV luminosity which comes from UV stellar clumps. Finally, the turbulence in actively star-forming galaxies is strongly connected to ISM conditions that favour an escape of Lymanα radiation (Herenz et al. 2016;Kimm et al. 2019). In this scenario, the offset could be intrinsic or affected by dust effects (Behrens et al. 2019), which could locally obstruct the emission of Lyman-α photons and thus produce a small spatial offset. In their sample, Hoag et al. (2019) noticed that the less dusty galaxies present on average a larger spatial  Fig. 17. Two examples of simulated galaxies. From left to right: high resolution image (UV emission), simulated HST image, simulated HST image of the same galaxy lensed by a LLAMAS cluster and HST image of a real LLAMAS galaxy. On each image, the contours represent the Lyman-α emission. The orange stars show the location of the centroid of the UV emission and the pink circle the location of the Lyman-α centroid. The first row presents an example of spatial offset produced by an offsetted UV bright clump. The second row presents an example of a spatial offset produced by a faint UV component spatially coincident with a strong Lyman-α emission peak. offset. We do not report a clear trend between the UV slope β and the spatial offset in the LLAMAS galaxies but we observed that the galaxies with larger offsets present higher β. Finally, the Lyman-α photons are produced in the vicinity of young short-lived stars while the 1500−Å UV emission arises on longer timescales and could be dominated by more evolved massive stars which produce less Lyman-α photons. An external young SFR cluster in the outskirts of the main UV component will therefore produce a spatial offset with the UV total light centroid and a spatial offset between UV and Lyman-α emission. 2. Lyman-α emission from faint UV satellites producing larger offsets values. This scenario was already proposed to explain either the formation of very extended Lyman-α haloes (Mitchell et al. 2021;Byrohl et al. 2021), or spatial offsets at z ∼ 3 − 7 Maiolino et al. 2015). Lemaux et al. (2020) measured a correlation between the UV brightness and the spatial offsets in kiloparsec as we measured in the LLAMAS galaxies with the SFR. The UV brightest galaxies are also the most massive, they are therefore more likely to reside in a more massive dark matter halo and thus be surrounded by faint satellite or merger galaxies. 3. Scattering effects of the Lyman-α photons across an inhomogeneous medium in such a way that Lyman-α emission seems to be offseted from the UV counterpart. This is more likely for the small offsets observed (∆ ell < 2). However, as the brightness of scattered Lyman-α emission decreases as a function of 1/r 2 , scattering effects alone are unlikely to produce the largest spatial offset such as those measured in the LLAMAS galaxies. The presence of ionised or low column density channels in the ISM and CGM, produced for instance by stellar feedbacks (Rivera-Thorsen et al. 2017;Erb et al. 2019;Reddy et al. 2021), could also produce this type spatial offsets.
Each offset measured can also be produced by a combination of several of these phenomena, as proposed by Matthee et al. (2020b) to explain a spatial offset measured in a z = 6.6 galaxy. Another way to try to distinguish all these scenarios is to study the spatially resolved properties of the lines for the most extended objects. For example, Erb et al. (2019) measured in a spatial offset of 650 pc between UV and Lyman-α emission in a lensed galaxy at z = 1.84 extended by 12 arcseconds. They explained this offset by a significant variation of the neutral hydrogen column density across the object, which supports a model in which ionizing radiation escapes from galaxies through channels with low column density of neutral gas. In a similar way,  identified 2 Lyman-α nebulae spatially offset from the associated star-forming regions. The variation of the Lyman-α surface brightness suggests large spatial fluctuations in the gas properties, and their results on spatial variations of the Lyman-α line profile, support a scenario in which high column density gas is driven toward up to 10 kpc. They conclude that the Lyman-α photons originate from a combination of resonant scattering from the star-forming regions and recombination radiation due to escaping ionizing photons, but they were unable to determine the relative contribution of these two mechanisms. A detailed study of the spectral and spatial properties of most extended Lyman-α haloes, as performed in Leclercq et al. (2020); Smit et al. (2017); Claeyssens et al. (2019) will allow us to detect potential variations of the CGM gas properties, such as hydrogen column density and kinematics, across the halo. In Claeyssens et al. (2019), we studied the spatial variation of the Lyman-α line within a lensed halo at sub-kpc scale. We identified a region, in the outskirts of the halo, presenting a smaller spectral Lyman-α line shift (with respect to the systemic redshift) than the rest of the extended emission. The local emission of Lyman-α photons by a faint, non-detected, UV component could explain the presence of weakly scattered Lyman-α photons at this location.

Summary and conclusions
We presented the largest statistical sample of lensed Lyman-α emitters observed with MUSE and HST. We observed 603 sources (producing 959 images) lensed by 17 different massive Article number, page 20 of 22 A. Claeyssens et al.: The Lensed Lyman-Alpha MUSE Arcs Sample (LLAMAS) clusters. Thanks to the lensing magnification (ranging from 1.4 to 40 in the total sample), we characterized the spatial properties of this new population of small and faint LAEs. We observe that 97% of the LLAMAS galaxies present an extended Lyman-α halo. We measured that the spatial extent of Lyman-α emission seem to be correlated with the UV SFR and the FWHM of the line. We confirmed the correlation from Leclercq et al. (2017) between Lyman-α and UV emission spatial extent and extended it to fainter LAEs but with a higher dispersion. We measured also the axis ratio of the UV continuum and Lyman-α emission distribution and notice that the 48% of the halo present an elliptical morphology (this fraction increase if we consider only the most resolved haloes). The Lyman-α haloes are on average less elliptical than the UV emission. We measured secure spatial offsets between the UV and Lyman-α emissions for 63% of the sources, and found a distribution of values in very good agreement with Hoag et al. (2019); Lemaux et al. (2020); Ribeiro et al. (2020) with a median value of ∆ = 0.58 kpc. We found very small or no correlations between the offset measurements and the physical parameters of the host galaxies (UV star formation rate, Lyman-α equivalent width and Lyman-α line FWHM). We identified 2 regimes in the offset distribution. First, galaxies with small offsets values (38%) with respect to the UV emission distribution, are likely due to bright star formation clumps in the outskirts of the UV component, emitting a strong Lyman-α emission. For the 62% other sources showing larger offsets, many scenarios could explain the large offsets such as inflows of gas, scattering effects of the photons in the CGM, extinction, outflows or Lyman-α emission from faint satellite galaxies not detected in UV. This last scenario is supported by the fact that we found higher values of r Lyα /r UV for galaxies with higher elliptical distances. Finally we compared our results with a zoom-in RHD simulation, following one typical faint Lyman-α emitter from z = 6 to z = 3, by producing, both in UV and Lyman-α emission, high resolution, "UDF-like" and "LLAMAS-like" mocks. The simulated galaxy is representative of the LLAMAS sample in terms of UV magnitude and Lyman-α halo size. We measured a similar spatial offsets distribution for the 3 samples and the LLAMAS galaxies. The simulation favors the interpretation where substructures in star-forming galaxies account for the smaller offsets (with ∆ ell < 2) and satellites/merger galaxies explain the larger offsets (with ∆ ell > 2). The scattering of Lyman-α photons could also contribute to production of spatial offsets.
It is likely that future works on these galaxies, especially the study of the spatial variations of emission lines profiles, and future observations of lower redshift LAEs with BlueMUSE , will help us understand more about the nature and the origin of the spatial offsets and Lyman-α haloes.