© A. Claeyssens et al. 2022

1. Introduction

The existence of bright Lyman-α radiation emitted by galaxies was originally predicted by Partridge & Peebles (1967) and progressively it 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 is known to produce extended halos surrounding galaxies up to 30 kpc (Matsuda et al. 2012; Momose et al. 2014), with the exception of quasars, however.

The origin of extended Lyman-α halos surrounding galaxies is still unknown, with two main hypotheses under consideration: (1) scattering of Lyman-α photons, produced mostly within star-forming regions, through the interstellar (hereafter ISM) and (2) circum-galactic medium or in situ photoionisation emission or collisional emission in the CGM (Mitchell et al. 2021). Lyman-α halos therefore represent a powerful probe of the hydrogen gas within the CGM, tracing both the spatial extent and velocity structure of the gas that surrounds galaxies and thus aiding in the investigation of 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-α halos is crucial for improving our understanding of the nature of the CGM at low and high redshifts. The LARS collaboration (for Lyman-α Reference Sample, Hayes et al. 2013, 2014; Östlin et al. 2014) characterized 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 typically lead to z > 2 galaxies appearing 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⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻². 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-α halos with typical exponential scale lengths of ∼5 − 10 kpc. However, a study of the possible diversity among individual Lyman-α halos 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. These authors presented a systematic morphological study of 145 of these Lyman-α halos 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; Chen et al. 2021b).

Scenarios explaining 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 and with the peak separation strongly dependent on the kinematic of the gas and neutral hydrogen column density (Verhamme et al. 2006, 2017; Henry et al. 2015; Gronke et al. 2016). In addition, emission from satellite galaxies and in situ emission probably contribute to the spatial extent and the clumpy morphology of the Lyman-α halos (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 these authors to distinguish among the aforementioned scenarios, by comparing these data with a zoom-in cosmological simulation, Mitchell et al. (2021) demonstrated that simulated Lyman-α halos 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 & Verhamme 2008; Verhamme et al. 2008; Schaerer et al. 2011; Gronke et al. 2015; Yang et al. 2016; Gronke 2017; Song et al. 2020).

One outstanding feature of Lyman-α halos is the spatial offset between UV continuum and Lyman-α emission peaks repeatedly reported in the literature (Shibuya et al. 2014a; Hoag et al. 2019; Lemaux et al. 2021; Ribeiro et al. 2020). Most studies so far have used long-slit spectroscopy to measure offsets between UV and Lyman-α emission, but these observations are restricted to one dimension (1D) only. Taking advantage of three-dimensional (3D) spectroscopy, IFUs are a more efficient tool in providing 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 offer a great opportunity to study the plausibility of different scenarios.

Multiple studies have focused on bright or extreme Lyman-α halos, such as Swinbank et al. (2015), who noted a variation of 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. (2020a,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 halos 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-α halos 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) in these studies has left the physical interpretation wide open, and some questions are still pending, particularly regarding: the origins of the Lyα photons as well as the physical mechanisms (e.g., outflows, inflows, satellites galaxies) responsible for the extent of Lyman-α halos 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, 2020; Smit et al. 2017; Claeyssens et al. 2019; Chen et al. 2021b). 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 and outflows). In addition, 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 any general conclusions.

With this in mind, we constructed 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), a total of 141 h were obtained mainly through MUSE Guaranteed Time Observation (GTO) on 18 different fields. We constructed a unique sample of 603 lensed LAEs (forming 959 images) at z = 2.9 − 6.6. Their strong lensing properties were characterized 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 or 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-α halos 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 compared 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 adopted a Λ cold dark matter cosmology with Ω_Λ = 0.7, Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. 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 here for the first time: MACS0451 and MACS0520 are completely new, while A2390, A2667 (de La Vieuville et al. 2019) and AS1063 (Mercurio et al. 2021) have been the subject of previous studies. This section briefly describes the method presented in R21 and the LAEs selections based on the global catalogs.

2.1. MUSE data

The MUSE data reduction procedure details are given in R21, largely following the prescription described in Weilbacher et al. (2020) along with some specific improvements for crowded fields. Of the 17 clusters we explore here, 5 (MACS0451, MACS0520, A2390, A2667 and AS1063) were analyzed 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 two 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 h 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²; 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).

2.2. HST data

To complement the MUSE data we use the available high-resolution ACS/WFC and WFC3-IR images in the optical and near-infrared (NIR) 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) surveys, and therefore have deep HST observations taken in 12 and 6 filters, respectively. For six additional clusters, HST observations were obtained as part of MACS survey (PI: Ebeling) as well as follow-up HST programs (PIs: Bradac, Egami, Carton).

2.3. 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 as follows. We begin with the 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 source catalogue is then 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. Next, a spectrum for 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 optimize the signal-to-noise ratios (S/N) 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. Then, all MUSELET and PRIOR sources down to a limiting S/N in the MUSE continuum (averaged over the full MUSE wavelength range) are inspected individually and their redshift is 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). Lastly, the resulting catalogs are 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 a redshift confidence of > 1 (more details on the Lyman-α line identification can be found in R21 Sect. 3.5)¹. 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).

2.4. 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). It allows us to generate parametric models of each cluster’s total mass distribution, using numerous multiple images identified in the catalogs as constraints. The final model parameters and constraints are presented in Appendix B of R21. Each cluster mass model is optimized with between 7 and 100 multiple 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 used 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². The relative contribution of each cluster to the full survey covolume is provided in the last columns of Table 1.

2.5. 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 of 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 Sect. 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 the observation information on each field.

2.6. Global properties

Figure 1 shows the global properties of the LLAMAS: redshift, lensing magnification, UV magnitude, and UV star formation rate (SFR) distributions; the last two are corrected for lensing magnification. The grey area of each panel represents sources selected for the spatial fitting (see Sect. 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 Fig. 1 are the values computed in R21, based on the central position of each image. The UV absolute magnitude at rest-frame 1500 Å (M₁₅₀₀) 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₁₅₀₀ 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⁻¹, and the median value is SFR = 0.55 M_⊙yr⁻¹.

Fig. 1. Global properties of the galaxies from the LLAMAS. From left to right: distribution of redshift, lensing magnification, UV magnitude, and SFR (both corrected from lensing magnification). The grey area in each panel represents the sources selected for the spatial fitting (see Sect. 3.4). When more than one image is detected for the same source, the values of the most magnified images are presented in panels 3 and 4. The magnifications presented in this figure are the total magnification for each source, i.e., the sum of the global magnifications of each multiple image of the same galaxy.

3. Analysis

3.1. Narrow-band image construction and global spectral properties

To characterize the Lyman-α emission of each LAE, we first constructed a 5″  ×  5″ Lyman-α narrow-band (hereafter, NB) image of each object from the MUSE datacube. When a source produced multiple images, we studied 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 (Shibuya et al. 2014b) given by:

$$\begin{array}{c}\hfill f(\lambda )=A\exp (\frac{-{(\lambda -{\lambda }_{\mathrm{peak}})}^2}{2{(a(\lambda -{\lambda }_{\mathrm{peak}})+d)}^2}),\end{array}$$(1)

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 optimized: 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 of [λ_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 10 000 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 optimized, 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 optimized (from λ_right − λ_left = 2.5 to 20 Å) in order to maximize the S/N as measured in a 0.7″ radius circular aperture (typical MUSE PSF size) centered on each image (on the catalogue source position). The start and end positions of the spectral band for continuum level estimation were chosen to be close to the Lyman-α line while not containing any Lyman-α emission ([λ_left − 34 Å :λ_left − 6.25 Å] on the left side and [λ_right + 6.25 Å : λ_right + 34 Å]) on the right side. Following this procedure, we obtained NB images with spectral bandwidths ranging from 3.75 to 17.5 Å depending on the galaxy (which represent 3 to 14 MUSE spectral pixels, respectively). The mean spectral width is 6.9 Å which represents 345 km s⁻¹ at z = 4.

Third, when 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-α optimized 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⁻¹⁹ erg s⁻¹ cm⁻² arcsec⁻² in any of the tested NB images, the object was rejected from the sample (in total, 12 images for 5 rejected objects).

3.2. Forward modeling

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-modeling approach based on parametric source models, including both lensing and instrumental effects. It consists of generating parametric source models (based on user assumptions on the emission profile) in the source plane, then lensing it by the best cluster model, convolving with the point spread funtion (PSF), and re-gridding it to the spatial sampling of the observation to compare it with the observation. In this study, we only used Sersic profiles, and specifically exponential (Sersic index n = 1) ones, as described in Wisotzki et al. (2016) and Leclercq et al. (2017). The free parameters are: the 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 the additional effects of lensing model uncertainties on the source reconstruction. We applied this fit both to 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 manually constructed a contour around each image, extended enough to cover all significant flux pixels (see Fig. 2) and a large (at least 10 MUSE pixels around the image) empty area around it. Only the pixels inside this region were considered for the χ² calculation. A weight image was associated with 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/σ². 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 χ² estimate in LENSTOOL is then:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{p,q}{(I_{p,q}-(M_{p,q}+\mathit{sky}))}^2{W}_{p,q},}\end{array}$$(2)

Fig. 2. Examples of LLAMAS galaxies and the best spatial 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 arcsec. We indicate 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.

with I_p, q, M_p, q, and W_p, q, respectively, the value in pixels [p, q] of the observed model and weight images.

3.3. PSF estimation

Because our aim is to model the morphological properties of the LAEs, we have to start by obtaining a very good knowledge of the PSF values in both the 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 Eq. (1) in R21. We constructed MUSE PSF following the same procedure described in Bacon et al. (2017). The approach to determine PSF FWHM variations with wavelength for each cluster is described in R21. For the HST images, we modeled the PSF in each filter used for spatial fitting (F555W, F606W, F814W, F110W and F125W). We used at least five 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.

3.4. 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 (four 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 modeling 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). One particular region of the plot stands out visually, where 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 Fig. 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 characterized both in UV and Lyman-α emission (representing 206 images).

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 Sect. 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).

The green line on Fig. 3, represents the criterion used to select sources for which morphology best-fit parameters are reliably recovered. To determine whether 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, and so, 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_{\mathrm{eps}}=\frac{(e_1+e_2-2(e_1\cap e_2))}{e_1}$, with e₁ the area of the input ellipse, e₂ the area of the best-fit ellipse, and e₁ ∩ e₂ 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 to a sufficient extent when this R_eps < 0.3, that is, when the error on morphology of the source concerns less than 30% of the total source area. The S/N value used for this selection is integrated on the complete image. Images having both a good S/N (between 10 and 35) and a large number of pixels (between 200 and 1000) cannot be fit well as the S/N of each individual pixel is relatively low.

3.5. Modeling 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) have mainly characterized the morphological properties of the Lyman-α nebulae through one- or two-component circular exponential models. Leclercq et al. (2017) decomposed the spatial fitting of the Lyman-α emission in two steps: first, they fit the UV emission on HST images using a single circular exponential profile. Then they fit the Lyman-α emission on MUSE NB images using a combination of two circular exponential profiles, at the same location, one with the scale radius fixed to the UV one. This modeling reproduced the objects of these studies well, with good residuals. However, in our study, thanks to the lensing magnification, we are able to observe LAEs with an improved spatial resolution. After applying this two-component circular exponential 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, when the object is detected in HST with a enough S/N value to be selected for the spatial fitting (see Fig. 3), we performed two 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 fit nine models on the MUSE Lyman-α emission:

• M3. The classical two-component circular exponential profile based on UV modeling (five free parameters).

• M4. The same two-component exponential profile with elliptical profiles based on the elliptical UV modeling (seven free parameters).

• M5. The M1 profile with the possibility to adjust the UV-like component scale radius (six free parameters).

• M6. The M3 profile with the possibility to adjust the UV-like component scale radius, ellipticity, and position angle (10 free parameters).

• M7. Two-component circular exponential profiles allow for centroids to be obtained at two different locations (8 free parameters).

• M8. Two-component elliptical exponential profiles allow for centroids to be obtained at two different locations (12 free parameters).

• M9. A single-component circular exponential profile (4 free parameters).

• M10. A single-component elliptical exponential profile (6 free parameters).

• M11. A single-component Sersic elliptical profile (7 free parameters).

The complete description of free and fixed parameters used in each model is presented in Table 2.

To disentangle the type of modeling that is better adapted to each object, we took into account the number of constraints, the number of free parameters, and the final χ², using a Bayesian Information Criterion (BIC), defined as: BIC = −2ln(ℒ)−k ln(N) with ℒ the likelihood of the best-fit, k as the number of free parameters, and N the number of constraints. We the fit with the minimum BIC kept as the best for the purposes of modeling. Thanks to these different modeling types, 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 into three categories: two components fixed at the same location, two components free to vary, and the single-component models. Ultimately, 67% of the LLAMAS galaxies are well described with a two-component fixed model (but the two-component circular exponential model is chosen only for 15%). Then, 12% of the galaxies are best described by two free components, which correspond to LAEs presenting either two 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 the median S/N of the objects in each category of model and found that the two more complex fits (models 7 and 8, with two free components) have a median S/N that is twice as high as the median S/N of the other models (S/N = 25). Finally, we found that objects modeled 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 the S/N and magnification seem to impact the models 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 Sect. 3.1), we obtained a measurement 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 characterization of fainter LAEs than non-lensed studies. To estimate the more realistic intrinsic flux, we used 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 a magnification under 10, the two magnification values are very similar; for higher magnification values, the new measurement is, on average, two to ten times lower than the first estimate, which is 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.

Fig. 4. Lyman-α line flux distribution of the LLAMA galaxies and UDF Lyman-α halos presented in 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 Fig. 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-α halos presented in Leclercq et al. (2017).

4. Results

In this section, we present our results on the Lyman-α nebulae morphology in terms of the spatial extent and axis ratio of the Lyman-α halos, the spatial offsets between UV continuum, and the Lyman-α emissions distributions.

4.1. 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 used half-light (r₅₀) and 90%-light radii (r₉₀) to characterize both the UV and Lyman-α emissions. For the total sample, 38% of the LAEs have an HST counterpart bright enough to be spatially modeled with LENSTOOL. We estimate r₅₀ and r₉₀ on a source-plane image produced from the parametric model (minimizing the BIC) of each source. Depending on the type of model, we measured the two radii from elliptical or circular rings. For objects with a circular best model (models M3, M5, and M9), we took ϵ = 0; for elliptical models M4, M6, and M11, we took the axis ratio value of the brightest component; and for the models M7, M8 (with two components at different locations) or M10, we used the mean axis ratio measured from the model M10 (one elliptical exponential component). We randomly produced 200 source-plane images selected from the LENSTOOL MCMC samples and the error bars were obtained by measuring r₅₀ and r₉₀ 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₅₀ and r₉₀ radius on the best-model images of the Lyman-α halos from Leclercq et al. (2017).

Based on the r_50, Lyα and r_90, Lyα, we measured the concentration parameter of the Lyman-α emission, which represents how compact the Lyman-α light profile is. We measure values ranging from c_Lyα = r₉₀/r₅₀ = 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 halos are also the smaller ones.

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 circularized 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₉₀ = 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 UV-like component. These values show that there is a large diversity of Lyman-α emission concentrations (very diffuse or peaked); however, 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 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₅₀ = r_50, Lyα/r_50, UV, with x_50,μ,LLAMAS = 12.0 and x_50,μ,UDF = 4.6.

Fig. 5. Comparison of Lyman-α and UV circularized 90%-light radii. The figure shows the Lyman-α 90%-light radius. LLAMAS and UDF (Leclercq et al. 2017) galaxies are represented in red and grey, respectively. The UDF and LLAMAS values are both correlated (Spearman weighted p-value of p0 = 0.02 and p0 = 0.05 for UDF and LLAMAS, respectively). Both datasets are roughly centered around rLyα = 10 rUV.

From the spatial modeling, we obtained 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 modeled in the UV emission, 38 % prefer a circular best model and 62% an elliptical one. For Lyman-α emission, 52% of the halos are better described by a circular model and 48% by 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 measured 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 the S/N weighted mean, as the magnification value is strongly correlated with the spatial integrated S/N). When we consider only images located inside the green contour in the second panel of Fig. 3, we find a median of ⟨q⟩ = 0.38. Applying the same procedure to measuring the distribution of Lyman-α emission axis ratios, we find that Lyman-α halos 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 Fig. 6. These results are consistent with the trend presented in Chen et al. (2021b) 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 used the axis ratio from the best source-plane model and applied 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. Halos with a small axis ratio value indicate that the CGM is structured along a preferred direction around the UV source.

Fig. 6. Distribution of axis-ratio (q = b/a) of UV (blue) and Lyman-α (pink) emission in the source plane. For multiple systems, we measured the magnification weighted mean of the axis ratio values of only the multiple images located inside the green contour of Fig. 3. The mean values of the two subsets are indicated on the plot.

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². The more a Lyman-α halo is resolved (with a high-detection map area that is strongly correlated to the lensing magnification and S/N), the more an elliptical model would be preferred in this case. The same effect is observed in S/N, since the detection map area and S/N are strongly correlated (see Fig. 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 Fig. 7.

Fig. 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 (r50) 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 two-image systems and the colored points show the three (or more) image systems. In the left and middle panels, the stars represents UV r50 and q, respectively, along with the 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 the dashed lines to represent the 10%, 50%, and 100% error levels.

For the spatial extent measurements (left panel of the Fig. 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, four multiple 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 the images are very broad, 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 Fig. 7).

4.2. Variations in the spatial extent of Lyman-α emission

In Fig. 5, we observe a correlation between Lyman-α and UV continuum r₉₀ for both UDF and LLAMAS datasets. The error-bar 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 characterized 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 values 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²) present no correlations between Lyman-α and UV r₉₀ (p₀ > 0.2 and ρ < 0.1). On the contrary, the two bins of medium (between 6.1 and 10.2 arcsec²) and highly (between 10.2 and 44.2 arcsec²) extended images show a significant correlation (p₀ = 0.033, ρ = 0.36, and p₀ = 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 S/N, 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 the properties of the Lyman-α line profile); thus we found three more or less significant correlations presented in Fig. 8. First, we found a very significant positive correlation between r_50, Lyα and the FWHM of the Lyman-α line (with a weighted Pearson coefficient of ρ = 0.35 and p-value of p₀ = 0.0002). Secondly, we observe a negative trend between r_50, Lyα and the Lyman-α rest-frame equivalent width W₀ (with weighted Pearson coefficient of ρ = −0.25 and p-value of p₀ = 0.07; the more extended the Lyman-α halo, the smaller is W₀. Finally we observe a weaker positive correlation between r_50, Lyα and the UV SFR (ρ = 0.18 and p₀ = 0.0002). We observe similar correlations between these three parameters and r_90, Lyα.

Fig. 8. r50, Lyα with respect to (from left to right) the UV SFR, Lyman-α rest-frame equivalent width (W0), and the velocity dispersion of the global Lyman-α line (σ) on the x-axis. We indicate the weighted Pearson coefficient, ρμ, and the associated p-value, p0, at the top of each panel.

4.3. 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 Sect. 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 measured the peak location after a spatial Gaussian filtering of FWHM = 0.2 or 0.1″ for Lyman-α NB and HST images, respectively. We ray-traced these values into the source plane to measure the intrinsic (physical) spatial offset between the two emission peaks (in kpc). The distribution of these “peak to peak” offset measurements is presented in dark blue in Fig. 10. For the total sample, we measured a median value of Δ(peaks) = 0.67 ± 0.13 kpc. When we kept only spatial offsets larger than two 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 two Lyman-α emission peaks (M7 and M8 in Table 2), we measured the offset between the UV component centroid and the two Lyman-α peaks and kept only the smallest. The distribution of these offset measurements is shown in light blue in Fig. 10, with the median value of the total sample is Δ(centroids) = 0.58 ± 0.14 kpc. Among the total sample, we identified three 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 were able top 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. (2021) 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 found 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, for 64 objects (mix of lensed and non-lensed galaxies) with 5 < z < 7, Lemaux et al. (2021) measured median offset of 0.61 ± 0.05 kpc. In LLAMAS galaxies, we did not measure any significant variation of the offset values with redshift. Finally, we measured the impact of the lensing model to the offset measurements (see Fig. 7). We found that 75% of the multiple systems have a dispersion smaller than the mean value.

Fig. 9. Presentation of 12 objects from the LLAMA sample presenting different values of spatial offsets. For the four galaxies on the top row, we show also the source-plane reconstruction. The observed images are labelled with a green “I” and the source-plane modelisations with a blue “S”. In each panel, we show HST images (F814W or F606W depending of the redshift) or source-plane reconstruction. The MUSE Lyman-α emission is represented by the red, orange, yellow, and white contours which correspond at a smooth surface brightness level of 250, 100, 50, and 12.5 10−19 erg s−1 cm−2 arcsec−2, respectively. The ID, redshift, and magnification values are given for each object. The Δ value is the value of the spatial offset measured in the source plane between the UV and Lyman-α best model centroids. The blue cross indicates the UV position in the image plane and the pink circle the Lyman-α one. In all panels we add a scale (white line) of 1″ in the image-plane and of 1 kpc in the 4 source-planes images.

Fig. 10. Distributions of measured offsets in kpc in the source plane. 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 centroid and peaks offset distributions, respectively.

4.4. Significance of the offsets measurements

To measure the robustness of these measurements, we estimated the probability of measuring such offsets if Lyman-α and UV peaks are superposed. We used 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 the UV and Lyman-α models. We applied 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 pair, the offset value δ₉₅; measured in this way, it corresponds 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 δ₉₅ < Δ. We find that 63% of galaxies have a significant offset. Among the remaining 37%, we measure (using the same method) whether the offsets corresponding to the cumulative probability of 68% (δ₆₈ at one sigma) are smaller than Δ (corresponding to offset significant at only one sigma). We consider that all the galaxies with Δ < δ₆₈ are compatible with the non-offset scenario. For these galaxies, we used the δ₆₈ as an upper limit to the offset measurement. We measured this fraction of spatial offsets in three bins with 0 < Δ < 1 kpc, 1 < Δ < 2 kpc, and Δ > 2 kpc and we found 48%, 92%, and 92% in significant offsets, respectively. The fraction of measurements compatible with the non-offset hypothesis is of 37%, 5%, and 2% in the three same bins, respectively. We 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 measured the elliptical distance, Δ_ell (normalized by the UV emission 90% isocontour elliptical parameters), between the Lyman-α emission centroids and the UV emission centroid:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{ell}}=(\frac{X_{\mathrm{Ly}\alpha }}{R_{x,\mathrm{UV}}}{)}^2+(\frac{Y_{\mathrm{Ly}\alpha }}{R_{y,\mathrm{UV}}}{)}^2,\end{array}$$(3)

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, as 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 obtained a result confirming that 40% of the galaxies present an external spatial offset as high as Δ_ell > 2.

Fig. 11. Distribution of elliptical distance (cf Sect. 4.4) measured in the source plane between Lyman-α emission centroids and the ellipse formed by the UV emission distribution (using r90,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).

We measured no correlation (p-values > 0.05) between the values of elliptical distance and the physical parameters of the host galaxies (LyαW₀, flux and luminosity, Δ in kpc). After dividing their sample of galaxies into two bins, based on W₀, Hoag et al. (2019) found a significant trend: galaxies with lower W₀ values show higher offset values (mean of 1.92 ± 0.13 kpc) than galaxies with higher W₀ (mean of 1.51 ± 0.11 kpc). We did not confirm this trend, either with an offset value in kpc nor with the elliptical distance. Lemaux et al. (2021) observed that the spatial offset increases with the UV brightness of the galaxies. We found 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 of Δ_ell < 3.9 (45 sources) show, on average, a higher UV SFR value (⟨SFR⟩ = 1.74 ± 0.23 M_⊙ yr⁻¹) than galaxies with Δ_ell > 3.9 (45 sources, ⟨SFR⟩ = 1.26 ± 0.28 M_⊙ yr⁻¹). 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 ($⟨\frac{r_{50,\mathrm{Ly}\alpha }}{r_{50,\mathrm{UV}}}⟩=25.5\pm 4.3$) than galaxies with smaller elliptical distances (⟨r_50, UV⟩ = 0.31 ± 0.4 kpc and $⟨\frac{r_{50,\mathrm{Ly}\alpha }}{r_{50,\mathrm{UV}}}⟩=7.0\pm 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.

In addition, 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 are expected to emit fewer 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.

5. Discussion

5.1. Production of mock observations

In order to physically interpret these results, we compared the measured LAE properties with a cosmological radiation hydrodynamical (RHD) simulation of a high-redshift galaxy evolving from z = 6 to z = 3, as 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, along with the collisional excitation of hydrogen as well as contribution from the UV background) and so it represents a powerful tool for studying Lyman-α photon 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¹⁰ 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 five 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 lines of sight 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 Myr. 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 Fig. 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. We 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).

Fig. 12. Global physical properties of the simulated, UDF and LLAMAS LAEs. From left to right: redshift, Lyman-α luminositie, UV magnitude, and SFR distributions. The Lyman-α luminosity and UV magnitude of the LLAMAS galaxies are intrinsic values measured (i.e., unlensed values). 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.

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⁻¹ to 1.23 M_⊙ yr⁻¹ between z = 6 and z = 3, with few SFR peak episodes (up to 3.0 M_⊙ yr⁻¹) 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⁻¹) is larger than the LLAMAS median value (SFR_LLAMAS = 0.48 M_⊙ yr⁻¹); however, 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⁻¹) ∼ 40 to 42.5 and 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 make comparisons 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 collapsed a cube containing all the Lyman-α emission without continuum. To reproduce non-lensed MUSE Lyman-α NB images of LAEs, we convolved the initial raw Lyman-α NB images by a typical MUSE UDF PSF, depending on the redshift of the galaxy, as described in Sect. 3.3. Finally we re-gridded 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 followed the same method for UV images: we convolved the raw UV images with UDF HST PSF, constructed following the method explained in Sect. 3.3, using HST F606W or F814W images (depending on the redshift). We finally re-gridded the images at the HST sampling of 0.05″  ×  0.05″. We added 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 followed 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 used the python package PHOTUTILS to detect the different images, following the same criteria used for the LLAMAS sample (see Sect. 2.5). When the lensing produces multiple images, we kept only one of them for comparison with the observations. We applied a S/N threshold of S/N = 6 for Lyman-α NB “UDF-like” images, as shown in Leclercq et al. (2017). For LLAMAS-like images, we applied the same selection for spatial fits as for real observations, highlighted by the blue contour in Fig. 3. We finally obtained a sample of 1164 raw images (both UV and Lyman-α NB), 271 “UDF-like” images and 254 “LLAMAS-like” (the othesr produced UDF-like and LLAMAS-like images that had too faint UV or Lyman-α images to be detected or spatially characterized).

5.2. Extended Lyman-α halos

5.2.1. Variations in the Lyman-α extent

We applied the same spatial fit as presented in Sect. 3.5 to the three types of mock observations (“raw” simulation, UDF-like and LLAMAS-like observations). For UV measurements, we applied a single elliptical exponential component fit on the three different mock observations. For the raw data and UDF-like observations we applied to the Lyman-α NB images a single elliptical exponential two-component fit (M4 in Table 2). For LLAMAS-like mock observations, in order to fairly compare them with the LLAMAS sample results, we applied three different fits to each selected image and compared BIC criteria to choose the best model of each galaxy. We chose to apply the models labeled 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 halos gives a similar range of concentrations (c = 1.16 − 32.5) as the real observations. Fig. 13 shows the distribution of Lyman-αr₉₀ 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 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).

Fig. 13. Lyman-α emission 90-light radius r90, Lyα as a function of the UV emission 90-light radius r90, UV for UDF galaxies (in grey), LLAMAS galaxies (in red), UDF simulated galaxies (empty grey circles), and LLAMAS simulated galaxies (empty red circles).

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.

Fig. 14. Distribution of the ratio of the measurements of r50, Lyα and r50, UV for UDF simulated (in grey) and LLAMAS simulated (in red) galaxies to the same measurements done on raw simulated data. The middle panel shows the distribution of the points for the 2 samples and the top and right panel the distribution of each parameter. The grey lines indicate δ50, Lyα = 1, δ50, UV = 1.

5.2.2. Origin of the extended Lyman-α emission

We notice in the simulations a strong effect on the part 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 rest-frame emission (i.e., 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-α halos 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 Fig. 8. Indeed, assuming that more extended Lyman-α halos 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).

Second, the correlation measured between r_50,Lyα and the UV SFR (Fig. 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. 2021); causing the halo expansion. This scenario could also explain the asymmetric and anisotropic Lyman-α emission distribution noticed in some LLAMAS galaxies. Rasekh et al. (2022) 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.

5.3. Spatial offsets between UV and Lyman-α emissions

5.3.1. 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 Fig. 15 shows the values of spatial offsets measured in raw simulations (grey) and LLAMAS-like sources (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 LLAMAS-like 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 non-lensed 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 Fig. 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 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 did not measure very high values (>20) of Δ_ell in the simulated galaxies. In the LLAMAS, we identified the galaxies with Δ_ell > 20 as being LBGs with a strong absorption feature observed in their spectra (one example is presented in the Fig. 10 at z = 4.69).

Fig. 15. Spatial offset distributions. 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.

5.3.2. 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, 2014; Ö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; Fö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 Fig. 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 Fig. 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 Fig. 9.

Fig. 16. Spatial offsets measured in the simulated galaxies. 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 r90,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.

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 Fig. 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 favor 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 Fig. 17, in 24% of the LLAMAS simulated galaxies which represent the larger spatial offsets. The small measured offsets (< 10 kpc) are suggestive mainly of cases of merger galaxies emission both in LLAMA and LLAMAS-like samples. We did not notice any significant trend between the spatial offset values and the physical properties of the LLAMAS-like simulated galaxies.

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. In 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.

5.3.3. Other possible origins and future direction

By comparing lensed LAE observations and a zoom-in simulation, we identified two different ranges of Δ values (cf. Fig. 11) showing that different scenarios are at play in the formation of spatial offsets between UV and Lyman-α emission at z > 3:

First, we considered the 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₀ is varying significantly as a function of the position in the galaxy. Recently, Rasekh et al. (2022) 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 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 coming from UV stellar clumps. Furthermore, the turbulence in actively star-forming galaxies is strongly connected to ISM conditions that favor the 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 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; whereas the 1500–Å UV emission arises on longer timescales and could be dominated by more evolved, massive stars that 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.

Second, we considered Lyman-α emission coming from faint UV satellites and producing larger offset values. This scenario was already proposed to explain either the formation of very extended Lyman-α halos (Mitchell et al. 2021; Byrohl et al. 2021), or spatial offsets at z  ∼  3 − 7 (Shibuya et al. 2014a; Maiolino et al. 2015). Lemaux et al. (2021) 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, to be surrounded by faint satellite or merger galaxies.

Next, we considered the 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², scattering effects alone are unlikely to produce the largest spatial offset such as those measured in the LLAMAS galaxies. The presence of ionized or low column density channels in the ISM and CGM that is produced, for instance, by stellar feedback (Rivera-Thorsen et al. 2017; Erb et al. 2019; Reddy et al. 2022) could also result in this type of spatial offset.

Each offset measured may 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 arcsec. These authors explained this offset by means of a significant variation in 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, Chen et al. (2021a) identified two Lyman-α nebulae spatially offset from the associated star-forming regions. The variation of the Lyman-α surface brightness suggests a large spatial fluctuations in the gas properties, and their results on spatial variations of the Lyman-α line profile support a scenario in which a high column density gas is driven toward up to 10 kpc. These authors 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-α halos, as performed in Leclercq et al. (2020), Smit et al. (2017), Claeyssens et al. (2019), would 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 on the outskirts of the halo that presents 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.

6. 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 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 observed that 97% of the LLAMAS galaxies present an extended Lyman-α halo. We found 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 also measured the axis ratio of the UV continuum and Lyman-α emission distribution and noticed that the 48% of the halo presents an elliptical morphology (this fraction increases if we consider only the most resolved halos). The Lyman-α halos 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. (2021), 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 two 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 three 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 satellite or 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 studies focused on these galaxies, especially the study of the spatial variations of emission lines profiles, along with future observations of lower redshift LAEs with BlueMUSE (Richard et al. 2019) will help us better understand the nature and the origin of the spatial offsets and Lyman-α halos.

Acknowledgments

A.C. and J.R. acknowledge support from the ERC starting grant 336736-CALENDS. T.G. and A.V. are supported by the ERC Starting grant 757258 “TRIPLE”. H.K., F.L. and A.V. acknowledge support from SNF Professorship PP00P2_176808. F.E.B. acknowledges support from ANID-Chile Basal AFB-170002, FONDECYT Regular 1200495 and 1190818, and Millennium Science Initiative Program – ICN12_009. H.K. acknowledges support from Japan Society for the Promotion of Science (JSPS) Overseas Research Fellowship. Tran Thi Thai was funded by Vingroup JSC and supported by the Master, PhD Scholarship Programme of Vingroup Innovation Foundation (VINIF), Institute of Big Data, code VINIF.2021.TS.041.

References

All Tables

All Figures

Fig. 1. Global properties of the galaxies from the LLAMAS. From left to right: distribution of redshift, lensing magnification, UV magnitude, and SFR (both corrected from lensing magnification). The grey area in each panel represents the sources selected for the spatial fitting (see Sect. 3.4). When more than one image is detected for the same source, the values of the most magnified images are presented in panels 3 and 4. The magnifications presented in this figure are the total magnification for each source, i.e., the sum of the global magnifications of each multiple image of the same galaxy.

In the text

Fig. 2. Examples of LLAMAS galaxies and the best spatial 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 arcsec. We indicate 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.

In the text

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 Sect. 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).

In the text

Fig. 4. Lyman-α line flux distribution of the LLAMA galaxies and UDF Lyman-α halos presented in 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 Fig. 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-α halos presented in Leclercq et al. (2017).

In the text

	Fig. 5. Comparison of Lyman-α and UV circularized 90%-light radii. The figure shows the Lyman-α 90%-light radius. LLAMAS and UDF (Leclercq et al. 2017) galaxies are represented in red and grey, respectively. The UDF and LLAMAS values are both correlated (Spearman weighted p-value of p0 = 0.02 and p0 = 0.05 for UDF and LLAMAS, respectively). Both datasets are roughly centered around rLyα = 10 rUV.
In the text

	Fig. 6. Distribution of axis-ratio (q = b/a) of UV (blue) and Lyman-α (pink) emission in the source plane. For multiple systems, we measured the magnification weighted mean of the axis ratio values of only the multiple images located inside the green contour of Fig. 3. The mean values of the two subsets are indicated on the plot.
In the text

Fig. 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 (r50) 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 two-image systems and the colored points show the three (or more) image systems. In the left and middle panels, the stars represents UV r50 and q, respectively, along with the 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 the dashed lines to represent the 10%, 50%, and 100% error levels.

In the text

	Fig. 8. r50, Lyα with respect to (from left to right) the UV SFR, Lyman-α rest-frame equivalent width (W0), and the velocity dispersion of the global Lyman-α line (σ) on the x-axis. We indicate the weighted Pearson coefficient, ρμ, and the associated p-value, p0, at the top of each panel.
In the text

Fig. 9. Presentation of 12 objects from the LLAMA sample presenting different values of spatial offsets. For the four galaxies on the top row, we show also the source-plane reconstruction. The observed images are labelled with a green “I” and the source-plane modelisations with a blue “S”. In each panel, we show HST images (F814W or F606W depending of the redshift) or source-plane reconstruction. The MUSE Lyman-α emission is represented by the red, orange, yellow, and white contours which correspond at a smooth surface brightness level of 250, 100, 50, and 12.5 10−19 erg s−1 cm−2 arcsec−2, respectively. The ID, redshift, and magnification values are given for each object. The Δ value is the value of the spatial offset measured in the source plane between the UV and Lyman-α best model centroids. The blue cross indicates the UV position in the image plane and the pink circle the Lyman-α one. In all panels we add a scale (white line) of 1″ in the image-plane and of 1 kpc in the 4 source-planes images.

In the text

Fig. 10. Distributions of measured offsets in kpc in the source plane. 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 centroid and peaks offset distributions, respectively.

In the text

Fig. 11. Distribution of elliptical distance (cf Sect. 4.4) measured in the source plane between Lyman-α emission centroids and the ellipse formed by the UV emission distribution (using r90,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).

In the text

Fig. 12. Global physical properties of the simulated, UDF and LLAMAS LAEs. From left to right: redshift, Lyman-α luminositie, UV magnitude, and SFR distributions. The Lyman-α luminosity and UV magnitude of the LLAMAS galaxies are intrinsic values measured (i.e., unlensed values). 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.

In the text

	Fig. 13. Lyman-α emission 90-light radius r90, Lyα as a function of the UV emission 90-light radius r90, UV for UDF galaxies (in grey), LLAMAS galaxies (in red), UDF simulated galaxies (empty grey circles), and LLAMAS simulated galaxies (empty red circles).
In the text

	Fig. 14. Distribution of the ratio of the measurements of r50, Lyα and r50, UV for UDF simulated (in grey) and LLAMAS simulated (in red) galaxies to the same measurements done on raw simulated data. The middle panel shows the distribution of the points for the 2 samples and the top and right panel the distribution of each parameter. The grey lines indicate δ50, Lyα = 1, δ50, UV = 1.
In the text

Fig. 15. Spatial offset distributions. 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.

In the text

Fig. 16. Spatial offsets measured in the simulated galaxies. 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 r90,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.

In the text

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. In 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.

In the text