Revealing the Nature of a Lyman-$\alpha$ Halo in a Strongly Lensed Interacting System at $z=2.92$

Spatially extended halos of H I Ly$\alpha$ emission are now ubiquitously found around high-redshift star-forming galaxies. But our understanding of the nature and powering mechanisms of these halos is still hampered by the complex radiative transfer effects of the Ly$\alpha$ line and limited angular resolution. In this paper, we present resolved Multi Unit Spectroscopic Explorer (MUSE) observations of SGAS J122651.3+215220, a strongly-lensed pair of $L^{*}$ galaxies at $z=2.92$ embedded in a Ly$\alpha$ halo of $L_{Ly\alpha}=(6.2\pm1.3)\times10^{42}$ erg s$^{-1}$. Globally, the system shows a line profile that is markedly asymmetric and redshifted, but its width and peak shift vary significantly across the halo. By fitting the spatially binned Ly$\alpha$ spectra with a collection of radiative transfer galactic wind models, we infer a mean outflow expansion velocity of $\approx 211$ km s$^{-1}$, with higher values preferentially found on both sides of the system's major axis. The velocity of the outflow is validated with the blueshift of low-ionization metal absorption lines in the spectra of the central galaxies. We also identify a faint ($M_{1500} \approx -16.7$) companion detected in both Ly$\alpha$ and the continuum, whose properties are in agreement with a predicted population of satellite galaxies that contribute to the extended Ly$\alpha$ emission. Finally, we briefly discuss the impact of the interaction between the central galaxies on the properties of the halo and the possibility of in situ fluorescent Ly$\alpha$ production.


INTRODUCTION
Galaxies are embedded in envelopes of a multiphase gas known as the circumgalactic medium (CGM). The CGM exists at an intermediate scale between the interstellar medium (ISM) and the intergalactic medium (IGM), and in the case of star forming galaxies (SFGs), it contains large reservoirs of cool (T ∼ 10 4 K) gas that feed and regulate the star formation activity in the host galaxy (see Tumlinson et al. 2017 for a review). Therefore, understanding the kinematics, physical conditions and spatial properties of the cool CGM is critical for answering the question of how galaxies evolve (e.g., Péroux & Howk 2020;and references therein). Historically, our knowledge of the CGM has been inferred from the statistical properties of intervening absorption arXiv:2206.02949v2 [astro-ph.GA] 16 Aug 2022 systems toward distant quasars (e.g., Churchill et al. 2000a,b;Péroux et al. 2003;Prochaska et al. 2003;Krogager et al. 2013;Nielsen et al. 2013;Sánchez-Ramírez et al. 2016), but this technique is not well suited for single-object studies, since all the information is obtained from one or very few lines of sight per galaxy. It is thus desirable to use spatially resolved spectroscopy of CGM emission lines to obtain a more complete picture of the gas properties in individual galaxies. Unfortunately, at z > 0.3 the H I 21 cm line, the canonical tracer of neutral gas, becomes too faint for individual detections. However, at z > 2.3 the H I Lyα transition at λ rest = 1215.67Å is shifted to visible wavelengths and has emerged as a powerful alternative for investigating the neutral phase of the CGM (e.g., Wisotzki et al. 2016). Besides its large intrinsic brightness, Lyα has gathered interest due to the discovery of diffuse emission extending beyond the stellar component of galaxies in deep narrowband (NB) imaging surveys (Steidel et al. 2000(Steidel et al. , 2011Cantalupo et al. 2012;Matsuda et al. 2012), which suggest that Lyα does indeed trace the CGM. These regions of extended Lyα emission around SFGs are also known in the literature as Lyα halos (LAHs; Hayashino et al. 2004;Steidel et al. 2011;Wisotzki et al. 2016;Leclercq et al. 2017;Chen et al. 2021b). With typical exponential scale lengths of 1 − 10 kpc and luminosities L Lyα 10 43 erg s −1 (Ouchi et al. 2020), LAHs should be distinguished from the larger and more luminous (but less abundant) Lyα "blobs" (e.g., Steidel et al. 2000;Matsuda et al. 2004;Ouchi et al. 2009;Borisova et al. 2016;Shibuya et al. 2018;Drake et al. 2020) that are typically associated with overdense regions containing several galaxies.
Due to its resonant nature, the Lyα line becomes optically thick at very low densities (N H > 10 13 cm −2 ; Ouchi et al. 2020), and thus in most environments it experiences complex radiative transfer effects that conceal the kinematics, column density, and ionization state of gas. As a consequence, researchers have attempted to infer the gas properties indirectly by combining spatially resolved spectroscopy (e.g., long-slit and first-generation integral-field unit (IFU) spectrographs) with models and simulations. However, this has been preferentially done on the brightest and most extended systems (i.e., Lyα blobs with L Lyα > 10 43 erg s −1 ), as fainter halos are more challenging to detect due to sensitivity limitations. In recent years, these limitations have been relaxed thanks to the availability of high-throughput integral-field spectrographs such as the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT) and the Keck Cosmic Web Imager (KCWI; Morrissey et al. 2018) at the Keck II telescope. With these instruments, galaxy-scale LAHs are now routinely detected and they are confirmed to be ubiquitous among high redshift SFGs (Wisotzki et al. 2016).
Theoretical models and simulations show that LAHs form by several mechanisms. On one hand, Lyα photons produced in central H II regions can propagate through the neutral gas out to the CGM in a series of resonant scattering events. On the other hand, Lyα photons in the CGM can also be produced in situ, by means of freebound collisional interactions, photoionization by internal or external sources, or star formation in satellite galaxies. Observationally, determining which mechanisms are at play is challenging, since the spectral shape of the line can be very similar for different underlying scenarios. Fortunately, simulations also predict LAHs to have a rich spatial substructure, featuring filaments, clumps, and satellites, as well as spatial variations in the spectral properties of the line (Mitchell et al. 2018;Behrens et al. 2019;Smith et al. 2019). For this reason, high angular resolution spectroscopy of LAHs is very valuable, since it can help us understand the physical nature of LAHs and potentially disentangle some of the degeneracies that affect the interpretation of the Lyα line.
In this context, a growing number of studies have exploited the power of strong gravitational lensing to resolve LAHs in a great level of detail (e.g., Swinbank et al. 2007;Karman et al. 2015;Caminha et al. 2016;Patrício et al. 2016;Smit et al. 2017;Erb et al. 2018;Claeyssens et al. 2019;Chen et al. 2021a;Claeyssens et al. 2022). This approach was pioneered by Swinbank et al. (2007), who obtained early IFU data of a giant lensed arc corresponding to a galaxy at z = 4.8. With MUSE, the efficiency of this technique was enhanced, providing the first tentative evidence of spatial variations in the line profile of a few lensed LAHs (Smit et al. 2017;Erb et al. 2018). Later on, Claeyssens et al. (2019) showed two examples of bright lensed LAHs with a robust measurement of the Lyα variation across the halo, finding in both cases a broader and redder line toward the outskirts of the halo. Although similar results have been obtained in a handful of nonlensed halos in the MUSE Hubble Ultra Deep Field (hereafter UDF; Leclercq et al. 2020), the spatial scales reached with lensing remain unrivaled.
In this paper we present MUSE observations of SGAS J122651.3+215220 (hereafter SGASJ1226), a lensed, multiply imaged pair of SFGs at z = 2.92. The main arc was discovered as a u-band dropout (Koester et al. 2010) and thanks to its high apparent brightness (r = 20.6 mag), it has been subject to several follow-up observa-tions (Wuyts et al. 2012;Saintonge et al. 2013;Malhotra et al. 2017;Gazagnes et al. 2018;Rigby et al. 2018;Chisholm et al. 2019;Solimano et al. 2021), becoming one of the best-studied Lyman-break galaxies (LBGs). Here, we report the discovery of an LAH associated with this galaxy and its close companion, which thanks to the lensing effect spans ∼ 20 on the sky. We take advantage of its extreme magnification (between µ ≈ 10 and µ ≈ 100 across the whole system) to spatially sample the Lyα line on subkiloparsec scales, thus offering a unique view of the CGM.
Throughout this paper, we adopt a flat ΛCDM cosmology with a matter density of Ω m, 0 = 0.3 and a Hubble parameter at z = 0 of H 0 = 70 km s −1 Mpc −1 . Unless otherwise specified, we will refer to physical (proper) distances rather than comoving distances. Also, all photometric magnitudes quoted in the paper are in the AB system. When relevant, we assume a universal Chabrier (2003) initial mass function.

MUSE
We observed SGASJ1226 with the MUSE (Bacon et al. 2010) mounted at the Unit Telescope 4 (Yepun) of the VLT. The data were taken between 2018 April and 2019 February as part of program 101.A-0364 (PI: López), using the Wide Field Mode with adaptive optics and an extended wavelength range (WFM-AO-E). This setup provides a field of view (FoV) of 1 × 1 with a pixel scale of 0. 2, and a wavelength range of 4600Å-9350Å with a resolving power of R ∼ 1770 at 4800Å (the wavelength of Lyα emission at z ∼ 2.92). A subset of these data were already presented by Tejos et al. (2021), and in this paper we follow similar reduction steps. The main difference is that Tejos et al. only combined 20 exposures (six were discarded due to slightly suboptimal seeing) to create the final stacked data cube, while here we used the full set of 26 exposures to reach deeper into the Lyα surface brightness (SB). All exposures have an exposure time of ∼ 640 s, resulting in a total exposure time of 16,640 s (∼ 4.6 hr). We found that including these additional exposures provides a 14% reduction in the 5σ noise level of the Lyα pseudo narrowband image, while the PSF full width at half maximum (FWHM; (0.72 ± 0.03) , see below) stays within the uncertainties of the value measured in the original cube ((0.76 ± 0.04) ).
The data was reduced using the MUSE pipeline (Weilbacher et al. 2012) within the ESO Recipe Execution Tool (EsoRex) environment (ESO CPL Development Team 2015), using standard procedures and calibrations. The wavelength solution was calibrated to vac-uum. We measured and cleaned residual sky contamination from the standard pipeline using the Zurich Atmosphere Purge code (Soto et al. 2016). We then aligned the World Coordinate Ssytem (WCS) of the data cube to match that of the HST imaging (see below) using a bright star in the FoV as reference. The effective PSF FWHM is (0.72 ± 0.03) , as measured from a Moffat fit to the bright star in a spectrally stacked image around the wavelength of Lyα at z = 2.92.
The MUSE pipeline propagates uncertainties in the individual pixel tables during the data cube creation, yielding a "variance cube" as part of the output. However, variance cubes generated in this way are known to underestimate the true variances of the data, because the algorithm neglects the spatial correlation in the noise introduced by the dithering and the slicer patterns (e.g., Bacon et al. 2017). Since the analysis presented in this paper required a good knowledge of the uncertainties, we estimated the true variances following the method outlined by Urrutia et al. (2019) and Weilbacher et al. (2020), but with a small modification. While they assumed the noise is spatially uniform (and hence only wavelength-dependent), we kept the spatial structure of the original variance cube, but scaled it up according to a correcting factor. The rationale behind this choice is to account for the sources' contribution to variance, which becomes relevant in the brightest regions of the LAH. To obtain the correction factor, we created a new pixel table populated with noise following a Gaussian distribution (mean=0, variance=1), and then fed it to the muse_scipost routine in the EsoRex pipeline to produce a mock noise data cube. We then measured the variance in sky regions of the resulting data cube to be close to 0.4 (the actual value is slightly wavelengthdependent, but we assume it is constant). So the correction factor is 1/0.4 = 2.5, which is then multiplied by the original variance cube. After this correction, the signal-to-noise ratio (S/N) distribution of sky apertures is consistent with a variance of ∼ 1 while without the correction such values are ∼ 2.
After visual inspection of the data cube, it became apparent that the central galaxy of the lensing cluster (hereafter BCG) contributes significantly to the background light near the arcs (see Fig. 1). Due to the position, orientation, and angular extension of the arcs with respect to the BCG, modeling and subtracting the background contamination in apertures would have been impractical and prone to many systematic effects. Instead, we opted for the self-consistent approach of modeling the BCG with a parametric light profile as a function of wavelength. We achieved this using the Bulge-Disc Decomposition of IFU Data Cubes package (buddi; John- of an average channel map integrated from 4900 to 5700Å taken from the data cube. The lower panel shows the same view but the average was obtained from the BCG-subtracted residual cube created with buddi. The image does not contain significant features at the position of the BCG within the colormap cuts (±7σ). ston et al. 2017), which uses GalfitM (Häußler et al. 2013) to model the 2D light profile of a galaxy simultaneously across several wavelength slices. While buddi is usually used to cleanly extract the spectra of each component included in the model, it also provides a residual data cube where the light of the target has been subtracted and the foreground and background objects can be analyzed with minimal contamination. In this case, we modeled the BCG with a single Sérsic profile plus a central PSF component after masking the arcs and other unrelated sources. After extracting the best-fit model data cube, the residual data cube was obtained, in which the light of the arcs is free from contamination from the BCG (see bottom panel of Fig. 1). In what follows, we use this BCG-subtracted data cube to perform the Lyα analysis.

HST
Observations of SGASJ1226 were taken with the HST with the Advanced Camera for Surveys (ACS) in the broadband filters F606W and F814W, and with the infrared channel of the Wide Field Camera 3 (WFC3) in the F110W and F160W filters, as part of General Observer programs #12368 and #15378, respectively. We used DrizzlePac's AstroDrizzle routine to align and combine the calibrated exposures to a common grid with a pixel size of 0. 03 using a gaussian kernel with pixfrac=0.8. The images reach 5σ limiting magnitudes 1 at m 606 = 25.9, m 814 = 25.4, m 110 = 24.7, and m 160 = 24.6. To check the accuracy of the astrometric solution a posteriori, we crossmatched the point sources in the final ACS frames with the Gaia DR2 catalog (Gaia Collaboration et al. 2018) and found five sources with a mean shift of 0. 08 and no signs of rotation or major distortion.
A 30 ×20 cutout of the ACS F606W image is shown in the upper right panel of Fig. 2. The multiple images and arcs are labeled using the same nomenclature in Fig. 1. A.1 is the brightest (most highly magnified) component and corresponds to a twofold, almost symmetric pair of images of galaxy A on both sides of the lensing critical curve (dotted white line). A fainter counterimage (A.2) of the same galaxy is seen ∼ 15 to the east. The second galaxy of the system, B, has a single known image (arc B) stretching from east to west a few arcseconds below A.1. A zoomed-in view of arc B shown in the lower right panel of Fig. 2 reveals multiple knots of UV emission. Arc B is also the host of the unresolved 870 µm continuum detection (purple ellipse) reported in Solimano et al. (2021). Between arcs A.1 and B we indicate G1, the z = 0.77 foreground Mg II absorber studied in Mortensen et al. (2021) and Tejos et al. (2021).
At z = 2.92 the observed filters cover the UV continuum and the Balmer break at the source rest frame, sampling the spectral energy distribution (SED) of the young stellar populations. The high resolution imaging also provides exquisite morphological detail on the lensed arcs, revealing numerous clumps and substruc- tures. At the same time, it enables the identification of lensed image pairs and the positions of cluster members that are included in the development of the lens model (Dai et al. 2020;Tejos et al. 2021).
The complex morphology of the lensed galaxies in the SGASJ1226 system (see § 3.4.1) makes the standard photometric extraction techniques unsuitable. For this reason, we used manually defined ad hoc polygonal apertures to obtain the total flux from the arcs. To assess the systematic error associated with this method, four members of the team created several apertures for each arc based on the image with the broadest PSF (WFC3-IR/F160W). The flux standard deviation derived from 20 different contributed apertures is about 10%, three times larger than the statistical error inferred from linear propagation. The final image plane magnitudes are reported in Table 1 after correction from Milky Way reddening using the Schlafly & Finkbeiner (2011) extinction tables.

Lens model
A proper measurement of the intrinsic properties of strongly lensed galaxies requires accurate knowledge of the geometrical distortion and magnification produced by the lensing cluster. Good approximations of the amount and direction of the distortion at any given position are obtained by modeling the cluster as a collection of parametric dark matter profiles and then constraining the model to predict the positions of image pairs identified in the data. In this paper, we use the lens model presented in Tejos et al. (2021), and therefore we refer the reader to that work for further details. Briefly, the modeling was done using the Lenstool software (Jullo & Kneib 2009) following the procedure described in Sharon et al. (2020). The model was fitted using three images of galaxy A (two images in A.1 plus the counterimage A.2; see Fig. 2), with the position of individual clumps identified in the HST data serving as constraints. Both cluster-scale and galaxy-scale potentials were modeled as pseudoisothermal ellipsoidal mass distributions (PIEMDs), all located at z lens = 0.43. For the galaxy-scale potentials the center of each PIEMD was held fixed to match the optical centroids of the cluster members during the fitting process. The presence of a second strong-lensing cluster at the same redshift only 157 to the south of SGASJ1226 motivated the addition of a cluster-scale PIEMD at that position, contributing shear to the overall lensing potential. Finally, a perturber on top of the western part of arc A.1 and the z = 0.77 galaxy between A.1 and B (labeled G1 in Tejos et al. 2021;see Fig. 2) were included as individual components. The latter was treated as belonging to the same plane of the foreground cluster.
The resulting best-fit model provides a set of deflection matrices that prescribe the angular offsets (in both the R.A. and decl. axes) induced by the lens at any given position in the image plane. This model reproduces the positions of 26 constraints with an rms of 0. 08. Throughout this paper, we use the deflection matrices to reconstruct source plane positions and sizes.
In Fig. 3 we plot the lensing critical curve (i.e., the locus of maximal magnification) on top of the MUSE pseudo-NB image as a way to visualize the morphology of the lensing potential.

Host galaxy properties
In this section we give a characterization of the host galaxies in terms oftheir systemic velocity, mass, and luminosity. A crucial step to understanding the origins and kinematics of the Lyα-emitting gas is to secure the systemic redshifts for the galaxies within the halo. For galaxy A, numerous reports exist in the literature. First, in their discovery paper Koester et al. (2010) used UV absorption lines to get z ISM = 2.9233. Later, Wuyts et al. (2012) obtained Keck/NIRSPEC spectroscopy that allowed the measurement of rest-frame optical nebular lines to find z neb = 2.9257 ± 0.0004. The difference of ∼ 200 km s −1 between the emission and absorption solutions was then suggested as tentative evidence for outflows (Wuyts et al. 2012). For galaxy B, Gemini/GMOS spectroscopy has yielded z ISM = 2.9233 (Bayliss et al. 2011) but no nebular redshift is available. Here, we avoid the systematic effects arising from a mix of different instruments by measuring redshifts directly on the MUSE data cube. Unfortunately, the nebular lines accessible in the MUSE wavelength range are extremely faint, so we coadded all spaxels associated with each galaxy to improve the S/N, thereby losing Overview and details of the image plane morphology of the SGASJ1226 system. Upper right: 30 × 20 cutout of the HST ACS F606W image, with MUSE Lyα contours overplotted (pale red). The contours start at the 3σ = 2.27 × 10 −18 erg s −1 cm −2 arcsec −2 SB level of the smoothed NB image (see § 3.2) and increase in powers of √ 2. The thin, dotted white line segment approximately divides A.1 into the two symmetric halves of the arc. Lower right: Zoomed-in view of arc B that highlights the spatial offset between the UV and Lyα emission. For clarity, contours only show the top 12 √ 2, 24, 24 √ 2 and 48σ levels of Lyα SB. Top left: 6 . 1 cutout of the ACS image centered on a local maximum of Lyα SB, with contours starting at 6σ. Middle left: Galfit model of the two foreground elliptical galaxies near the local Lyα peak. Lower left: Residuals from the subtraction of the Galfit model. A compact excess of continuum appears near the center of the Lyα peak. The purple ellipse shows the location and deconvolved size of the dust continuum emission at 870 µm detected with the Atacama Compact Array (ACA; Solimano et al. 2021) any information on possible spatial velocity gradients. For each galaxy, we simultaneously fit the Si II* λ1533, [C III] λ1906 and C III] λ1908 lines on the continuumsubtracted spectrum with Gaussian profiles of a common width, varying only z and the line amplitudes. In this way we obtained z neb = 2.9257 ± 0.0001 and z neb = 2.9238 ± 0.0002 for A and B, respectively, corresponding to a radial velocity difference between A and B of (145 ± 17) km s −1 . The value for galaxy A is fully consistent with prior literature measurements (Wuyts et al. 2012;Rigby et al. 2018). We checked that our solutions are robust against the relaxation of the width constraint: letting Si II * λ1533 have a width parameter independent of the C III lines yields the same 16 th -50 th -84 th percentiles as the single-width run.
The general properties of the two brightest galaxies in the LAH were already reported in Solimano et al. (2021) by means of fitting SED models to the available broadband photometry. Briefly, the SED fits included photometry from the four HST bands mentioned above, Spitzer/IRAC 3.6 and 4.5 µm bands, and the ACA 870 µm continuum. The available Herschel bands were excluded due to severe blending and low S/N (see Saintonge et al. 2013). The magnified fluxes were then fit with the Magphys (da Cunha et al. 2008Cunha et al. , 2015 energy balance code to constrain the star formation rate (SFR), stellar mass, and infrared luminosity. The lack of midinfrared bands prevented the detection of continuum from dust heated by an active galactic nucleus (AGN).
Here, we reproduce the best-fit values and their uncertainties after demagnification in Table 2. The average magnification was computed as the ratio between the solid angle of the HST aperture in the image plane and the solid angle spanned by the smallest polygon that encloses all the delensed grid points in the source plane. The magnification uncertainty was fixed to 20% to account for typical systematic errors in the lens modeling (Raney et al. 2020;where statistical errors are less important). We tested this choice by computing the actual statistical uncertainties using a set of 100 realizations of the lens model sampled from a Markov Chain Monte Carlo (MCMC) chain. We found that for the apertures used here, the relative uncertainty in magnification is between 13% and 17%. We also estimated the rest-frame far-UV luminosity of the arcs directly from the ACS/F606W photometry tracing the stellar continuum at λ rest ≈ 1500Å at z = 2.92. As we did with all the available HST data, we divided the image plane fluxes by the average magnification factor of the most complete image. This means that despite A.1 having the largest S/N, our flux-dependent estimates for galaxy A come from the counterimage A.2 instead, since it has the largest footprint in the source plane. In this way, we used µ = 7.5 for A.2 and µ = 30.2 for the arc B to find log(L UV /L ) = log(λL λ ) = 10.80 ± 0.09 for galaxy A and log(L UV /L ) = 10.6 ± 0.1 for galaxy B, corresponding, respectively, to 1.0 ± 0.2 and 0.64 ± 0.15 times the typical UV luminosity at z = 3, L * = 6.2 × 10 10 L (e.g., Paltani et al. 2007;Reddy et al. 2008;Bian et al. 2013;Mehta et al. 2017).

Continuum Subtraction and Pseudo-NB Imaging
The first step in our analysis was to create NB Lyα images from the reduced MUSE cube. We first extracted a subcube between 4754 and 4810.3Å, a range that fully includes both the Lyα emission line and 10Å of the adjacent continuum on each side. In the UV-bright regions of the arc, the continuum shows a clear break at the wavelength of Lyα. At bluer wavelengths the continuum is strongly suppressed by a combination of damped self-absorption and decreased IGM transmission, while at redder wavelengths the emission is dominated by the UV power-law continuum but modulated by the H I damping wings. These effects produce an underlying continuum with a complex shape. In principle, one could model it as the superposition of the UV power law, a damped Voigt profile, and the break from the IGM transmission; however, the data lacks sufficient S/N to fit such a model on a spaxel-to-spaxel basis. Instead, we chose to model the continuum as a linear ramp between the blue and red continuum levels. This scheme has been successfully applied to similar datasets (e.g., Claeyssens et al. 2019). The functional form of this model is is the slope of the linear ramp. We estimated the mean blue (f blue ) and red (f red ) continuum levels by averaging the spectral channels to each side of the line from the subcube, at λ < λ 1 and λ > λ 2 respectively. We found that the values λ 1 = 4768.8Å and λ 2 = 4788.0Å produce robust continuum subtraction (average zero flux) on both sides of the line. Once the continuum was extracted from every spaxel, we integrated the resulting cube between 4769 and 4788Å (rest frame 1215.3 and 1219.6Å at z = 2.924) to obtain a pseudo-NB image of the Lyα emission. The integration limits were chosen to enclose the full spectral extent of the redshifted line. We excluded the blueshifted peak since it only appears in a limited image plane region and its contribution to the total flux is less than 1%. Fig. 3 shows the final image after smoothing with a Gaussian kernel of σ = 2 spaxels. The smoothing was applied only for the ease of visualizing the low-SB structure of the object, but in the subsequent analysis we used the unsmoothed version.

Image Plane Analysis
In order to obtain insights into the spatial Lyα properties of SGASJ1226 in a way that is independent of the lens model, we started our analysis in the image plane, rather than in the source plane. The NB image in Fig. 3 shows at least five strong peaks of Lyα SB that stand out from the diffuse emission. In aid of comparing these features to the UV continuum, we reproduce the Lyα SB contours of Fig. 3 in Fig. 2, placing them on top of the ACS F606W image, the band that traces the UV continuum at z = 2.92. The two northernmost local peaks are associated with images A.1 (west, bright) and A.2 (east, faint) of galaxy A. Toward the south, we observe the second brightest peak of Lyα SB, which is connected to a fainter arc-like structure extending about 5 to the east. The arc itself has two secondary peaks at the ≈ 34σ level. We associate the peak plus the arc with Arc B in the HST image, although there is an evident offset between the bright spots in Lyα and the location of the UV-bright clumps (see lower right panel of Fig. 2). In particular, the bright Lyα peak at B is offset by 1. 2 with respect to the UV centroid, which lies in the middle of the brightest knots of the arc. Also, the Lyα arc is offset by 0. 6 on average to the north of the UV arc B. Due to the achromatic nature of lensing, offsets in the image plane between two tracers (e.g., Lyα and UV here) imply intrinsic offsets in the source plane.
The other remarkable local peak is ∼ 5 southward from A.2 at 17σ above the sky background (see green box in Fig. 2). Interestingly, this peak is centered very close to two intervening foreground cluster members. Since no emission line is expected at λ obs = 4471Å at the redshift of the cluster (z lens = 0.43), we conclude that the signal comes from the same redshift of the LAH and is likely originates in a compact Lyα emitter (LAE) embedded in the halo (hereafter we refer to this Lyα source as SGASJ1226-LAE). Motivated by this hypothesis, we searched the HST data for a continuum counterpart, finding a candidate near the core of one of the cluster members in the F606W image. To confirm the presence of this counterpart, we employed Galfit (Peng et al. 2002(Peng et al. , 2010 to model the light of the two intervening galaxies with Sérsic profiles. After subtraction of the best-fit model, the residuals clearly show an excess emission within 0. 2 of the Lyα peak. The excess is also detected in the F814W image, but not in the near infrared filters. The limited spatial resolution in the F110W and F160W bands results in contamination the expected location of the galaxy by PSF subtraction residuals, preventing us from putting any constraint on its flux. We measured the flux in the residual images with SExtractor's (Bertin & Arnouts 1996) automatic apertures with the F606W residual as the detection frame, yielding (magnified) magnitudes of m 606 = 25.94 ± 0.17 and m 814 = 25.82 ± 0.30. From these two bands we inferred a UV slope of β = −1.6 ± 1.1. A Lyα flux of SGASJ1226-LAE of 7.5 × 10 −17 erg s −1 cm −2 µ −1 was measured by fitting an asymmetric Gaussian (AG) profile (see Appendix B) to the spectrum integrated in a 1 circular aperture. Extrapolating the continuum flux to λ rest = 1215.67Å we estimated a rest-frame Lyα equivalent width (EW) of (104 ± 19)Å, which falls in the classical definition of high-redshift LAEs (Ouchi et al. 2020). At the same redshift as the rest of SGASJ1226, this source is likely a satellite galaxy of the system. Details on the characterization of this source together with estimates of its contribution to the total Lyα luminosity are presented in § 4.2.3.

Integrated spectrum
In Fig. 4 we show the spatially integrated Lyα spectrum for three different apertures, one for each of the two arcs and another for the diffuse emission. To create the apertures, we first integrated the data cube at 1550Å < λ rest < 1650Å to obtain a map of the UV continuum at z = 2.92 at the same resolution, pixel scale, and astrometry of the Lyα NB image. Then, we applied a threshold of continuum S/N = 5 on this image to isolate spaxels containing bright UV emission. The resulting masks were manually inspected and tweaked to remove spurious spaxels that met the S/N threshold but were associated with unrelated sources. For galaxy A, the masks include all spaxels from both the arc A.1 and the counterimage A.2. The mask for the diffuse LAH was defined as all spaxels with Lyα SB above 3σ(smooth) = 2.2 × 10 −18 erg s −1 cm −2 arcsec −2 in the smoothed narrowband image. Then, the continuum masks were subtracted from the diffuse halo mask to exclude spaxels with continuum emission. The three apertures and their corresponding Lyα spectra are shown in Fig. 4.
The three profiles show a very strong red peak with no clear blue component. Also, the line profile is markedly asymmetric in the three cases, with a sharp drop from the peak to the blue and a broad red wing. The major difference between the three spectra is the width of the line and the location of the peak. We measured line properties such as the FWHM and the shift of the peak with respect to the systemic velocity of the system by fitting an AG profile to the spectrum after taking into account the line spread function (LSF; for details on this method see Appendix B). Galaxy A shows the narrowest line, with an observed-frame FWHM of 4.68 ± 0.60Å (equivalent to a rest-frame velocity of 294.2 ± 4.0 km s −1 ) and a 1.72 ± 0.03Å (108.2 ± 1.8 km s −1 ) offset between the peak and the systemic redshift z A sys = 2.9257, indicated by a dotted orange line. The Lyα profile emitted by galaxy B is broader (FWHM = 6.84 ± 0.10Å, 430.1 ± 6.5 km s −1 ), and, despite its having the same peak wavelength as galaxy A, the systemic redshift is lower (z B sys = 2.9238; dotted green line) and thus the velocity offset is 4.63 ± 0.06Å (291.5 ± 3.5 km s −1 ). Finally, the diffuse halo emission has the broadest line width (8.10 ± 0.07Å, 509 ± 4 km s −1 ) and a peak that is slightly displaced redward with respect to the peaks of A and B (4.31 ± 0.05Å, or 270.6 ± 3.0 km s −1 if we assume (z A sys + z B sys )/2 = 2.92475). We warn the reader that the relatively low uncertainties on the fitted parameters reflect the high S/N of the spectra, and thus do not include the systematic uncertainties arising from the unknown intrinsic profile shape.
The broader and redder line profile we observe in the diffuse halo relative to the central component agrees with previous findings by Claeyssens et al. (2019) for two other lensed LAHs (SMACS2031 and MACS0940) and by Leclercq et al. (2020) in a sample of unlensed LAHs in the UDF. In a scenario where the extended Lyα emission is explained exclusively by neutral gas scattering, the width of the line is linked to the number of scatterings those photons had to experience before escaping the halo in the direction of the observer. The broader the line, the more reprocessed are the photons in the highvelocity tail. In this context, our results would indicate that Lyα photons coming "down the barrel" from the SFGs are less reprocessed on average than photons that escape from the outskirts of the halo. Such an effect can naturally arise if the close environments of the galaxies have a higher ionization fraction that effectively reduces the optical depth along the line of sight. However, the Lyα signal depends not only on the neutral hydrogen column density, but also on the gas kinematics.

Morphology
In this section we present the spatial properties of the Lyα SB in the source plane. As we discussed in § 3.3, the NB data in the image plane already reveals that the Lyα emission is both spatially offset and more extended than the UV continuum. Qualitatively, these properties should be preserved by the lensing, but we now verify it in a quantitative way using our lens model. In what follows, we use the deflection matrices resampled to MUSE resolution (0. 2) unless otherwise specified.
We employed a Bayesian forward-modeling approach similar to the one presented in Claeyssens et al. (2022). We modeled both the MUSE UV and Lyα data for each of the two galaxies (A and B) as Sérsic profiles and their parameter space was explored by a MCMC sampling scheme using the emcee library (Foreman-Mackey et al. 2013). For each proposed set of parameters in the chain, the model was evaluated and traced back to the image plane according to the lens model prescription, convolved with the PSF, and compared to the data under a Gaussian likelihood. For the sake of simplicity, we used a single Sérsic profile per component with all its six parameters free and set uniform priors.
We first fit the MUSE UV continuum image at λ rest ∼ 1600Å defined in § 3.3.1. We chose to fit the UV model in the MUSE data rather than the ACS F606W data to have a more comparable data quality between the UV and Lyα datasets. After MCMC convergence we obtained a best fit Sérsic index and circularized effective radius of n A UV = 0.57 ± 0.04 and r A 50, UV = (0.97 ± 0.02) kpc, respectively, for galaxy A, while for galaxy B, n B UV = 1.4 ± 0.2 and r 50, UV = (4.6 ± 0.5) kpc. For completeness and direct comparison with Claeyssens et al. (2022), we also quote the 90%light radius r A 90, UV = (1.85 ± 0.05) kpc and r B 90, U V = 12.4 +2.2 −1.9 kpc. The resulting median distance 2 between the centers of A and B is (14.3 ± 1.4) kpc.
We In both the UV and Lyα cases the models cannot fully reproduce all of the image plane features, as revealed by the high-significance residuals of the fit (see Fig. D4 for the UV and Fig. 5 for Lyα). This can be a result of the clumpy nature of the galaxies and their halos, making the Sérsic profile unsuitable. A full morphological analysis of the individual clumps is outside the scope of this paper and will be presented elsewhere. Nevertheless, the fitted Sérsic parameters can be informative of the sizes and overall light distribution of the galaxies. For example, the high Sérsic indices indicate that the sources have a very compact core and extended tails, similar to the double-exponential profiles often invoked for describing LAHs (Wisotzki et al. 2016;Leclercq et al. 2017;e.g.,). The Lyα half-light radii, on the other hand, put the two sources in the top 10% of the Leclercq et al. (2017) sample and above any measurement in the Claeyssens et al. (2022) sample. In this context, the size of the SGASJ1226 system approaches the lower end of the size range of Lyα blobs (Ouchi et al. 2020).
Using the centroids of the models, we computed the intrinsic spatial offsets between Lyα and the UV to be ∆(Lyα − UV) A = (0.55 ± 0.09) kpc for galaxy A and ∆(Lyα − UV) B = (3.7 ± 0.2) kpc for galaxy B. Compared to the Claeyssens et al. (2022) sample, the offset of galaxy B ranks the largest at face value. However, when normalized to the size of the UV model (using the definition of "elliptical distance," ∆ ell , in Equation (3) of Claeyssens et al. (2022)), the offsets of A and B both qualify as "internal spatial offsets" for ∆ A ell ≈ 0.06 and ∆ B ell ≈ 0.02 are well below unity.

Spatial Distribution of Line Parameters
In this section we describe the procedure to extract and characterize the spectral line profile in resolved regions of the system. Following Claeyssens et al. (2019), we started from our best source plane model of Lyα emission (see § 3.4.1) and evaluated it in a grid covering the delensed coordinates, at a pixel size of 0. 03. The resulting image was then fed to the vorbin package (Cappellari & Copin 2003), which produces a tessellation of the source plane into regions of roughly equal flux.
With this method, we constructed two different tessellations: first, a high-resolution tessellation of 50 bins with a median S/N of ∼ 20 that we used for fitting the AG profiles and second, a lower-resolution tesselation of 24 spatial bins with a median S/N of ∼ 30, used to fit the galactic wind models (see § 3.4.3), which require higher S/N for producing reliable results. Then, for a given tessellation, we extracted the spectrum from each region by combining all the image plane spaxels that trace back to it and coadded the corresponding continuum-subtracted spectra. We made sure that the traced image plane regions have at least four contiguous spaxels and they are larger than the PSF at least in one direction.
We characterized the spectral properties of the redshifted Lyα line in each bin by fitting an AG profile (see Appendix B). An analysis of physically motivated models for the observed Lyα profiles is postponed to § 3.4.3.
Before fitting, each bin was assigned a systemic redshift based on its spatial overlap with UV emission of the galaxies. Bins having more than 50% of their UV flux inside the continuum masks defined in § 3.3.1 were assigned the systemic redshift of the corresponding galaxy (i.e., z A = 2.9257 for galaxy A and z B = 2.9238 for galaxy B). For all the other bins we set z mean = (z A + z B )/2 = 2.92475.
The fit is well behaved and the marginalized posterior distributions are not multimodal, except for a single bin of the 50-bin tessellation located at the outskirts of the halo whose spectrum does not have enough S/N for MCMC convergence. We used the integrated flux yielded by the fitted model to estimate both the intrinsic SB of each bin (dividing by the bin solid angle) and the total luminosity of the system (dividing by bin magnification 3 and taking the sum). In this way, we obtained L Lyα = (6.2 ± 1.3) × 10 42 erg s −1 for the whole SGASJ1226 system 4 based on the 50-bin tessellation. Fig. 6 shows the source plane tessellation with the bins color-coded according to the fitted SB, peak shift velocity, and line width. The leftmost panel shows the source plane distribution of SB, which peaks near the position of galaxies A and B (cyan contours) and rapidly declines toward the outskirts. However, instead of being two distinct halos, the two resolved Lyα SB peaks are connected by a low-SB "bridge" or filament. The peak velocity approximately ranges from 50 to 450 km s −1 across the halo, as seen in the middle panel of Fig. 6. The lowest values are the ones associated with galaxy A, with an uncertainty-weighted mean and standard deviation of 75 and 20 km s −1 respectively. The highest values are instead associated to either galaxy B or to bins at the outskirts of the halo. The right panel also shows a large spread in line FWHM, spanning from 265 to 690 km s −1 . We observe that the bins covering the UV continuum typically have samller widths than the bins of the diffuse halo, in agreement with the trends observed by Claeyssens et al. (2019) and Leclercq et al. (2020). Qualitatively, the presence of this pattern seems to confirm the results obtained above for the integrated apertures (see § 3.3.1). A remarkable exception, however, is the bin covering the faint companion SGASJ1226-LAE (indicated with a cyan star in Fig. 6), which exhibits the second largest line width of the system (≈ 670 km s −1 ). Similarities between the FWHM and peak shift maps should not come as a surprise. Recent studies have found that these two quantities are positively correlated, both among integrated measurements across different objects (Verhamme et al. 2018)   2020). From a theoretical perspective, the relation is expected to arise as a natural consequence of resonant scattering, as confirmed by radiative transfer calculations in simple geometries (e.g., Schaerer et al. 2011;Zheng & Wallace 2014;Song et al. 2020). In contrast, the relation is not found in more complex, high-resolution simulations that include full radiative hydrodynamics (Behrens et al. 2019;Mitchell et al. 2021). In Fig. 7 we show the FWHM versus peak shift values for the 49 spatial bins with good fits in the 50-bin tessellation, in comparison with the Verhamme et al. (2018) empirical relation. Our data also exhibits the FWHM-shift correlation with an uncertainty-weighted Pearson's coefficient of r = 0.76 +0.09 −0.21 (p < 10 −5 ), where the errors have been estimated with bootstrapping. Motivated by this result, we performed orthogonal linear regressions using SciPy's odr module on the data separated in the three systemic redshift groups. The best-fit lines are plotted in the left panel of Fig. 7 along with their corresponding 68% confidence intervals. Each line is described by a slope a and a zero-point b. We found that the slopes for the diffuse halo and galaxy B bins agree at a ≈ 1.25 within a 1σ error, while their zero-points differ significantly by about 120 km s −1 . These offsets put the resolved relations for B and the halo slightly below Verhamme et al.'s empirical relation (a = 0.9±0.14; b = −34 ± 60), but the slopes are still consistent within 1σ. Galaxy A's data, on the other hand, favors a much shallower relation (a = 0.51 ± 0.09), which cannot be brought into agreement with Verhamme et al.'s, but is within the broad range of slopes observed by Leclercq et al. (2020) in resolved individual halos. Additionally, we investigated the effect of SB in the line profile. Just as in the two objects studied by Claeyssens et al. (2019), we observe that above a certain SB threshold, brighter bins have smaller peak velocity shifts, as shown in the right panel of Fig. 7.

FLaREON Models
Here, we analyze the binned spectra with physically motivated models. We used the publicly available The simulations only include two radiative transfer effects, namely resonant scattering and dust absorption. Hence, for the interpretation of these models one has to assume that all the Lyα photons are produced in the central galaxies and then scattered away from resonance in the CGM. Other mechanisms for the production of extended Lyα emission, such as fluorescence and gravitational cooling, are discussed in § 4. Three geometries are currently available within FLaREON: The first two are a spherical, expanding thin shell of isothermal gas around a point source of Lyα photons (hereafter TS) and a galactic wind geometry (GW), which is identical to the TS geometry except for the density distribution of the gas. While the TS assumes all the gas is concentrated in a thin shell at a fixed distance of the source, in the GW geometry the source lies in an empty spherical cavity surrounded by a spherical distribution of isothermal gas with radially declining gas density. Finally, FLaREON also offers a biconical outflow geometry, which combines a static uniform medium with a bicone of lower-density, expanding gas. All geometries are governed by the same three parameters: the expansion velocity V exp , the neutral hydrogen column density log N H I and the pure absorption (dust) optical depth log τ .
While these models were originally developed to describe the profiles of spatially unresolved spectra of LAEs, here we use them to model individual bins of our source plane tessellation. The main caveat of this approach is that the emission source will not be, in most cases, at the center of the binned region. This means that a typical Lyα photon produced by the central galaxies will need to travel a larger distance (thus with a higher probability of absorption or scattering) to escape from any given bin than what the models imply. However, Chen et al. (2021a) proposed that if the Lyα photons do not originate inside the cloud, the problem can still be described with just half of the expanding sphere. Following their argument, the isotropy of the TS and GW geometries implies that a signal arising from a single hemisphere would have the same line profile as the full spectrum but with a reduction in amplitude. This property suggests that one can approximate the full halo as a collection of half-expanding clouds, where now the expansion velocity is measured relative to a reference point inside each cloud along the line of sight. Since the biconical wind geometry is not isotropic, we did not try that geometry in our modeling.
The code requires the input of a systemic redshift in order to transform observed wavelengths into rest-frame Marker shapes and colors separate the three components of the SGASJ1226 system, namely the cores of galaxies A (orange squares) and B (green crosses) and the halo diffuse emission (blue circles). For the latter, the systemic velocity was set to zmean = (zA + zB)/2 = 2.92475. The best fitting straight line is also displayed for each of these components, along with their 68% confidence intervals shown as a shaded area around the lines. The dashed gray line traces the empirical relation found by Verhamme et al. (2018). Right panel: Lyα peak shift velocity versus SB. Horizontal lines mark the uncertainty-weighted average peak shift velocity ∆v for each of the three components.
wavelengths. However, the choice of systemic redshift is known to have a significant influence on other parameters (Gurung-López et al. 2019). Due to this, one could in principle infer the systemic redshift from the Lyα profile shape alone. To explore this idea, we started by fitting FLaREON models with z sys as a free parameter to the integrated spectra of galaxies A and B (see § 3.3.1), because their systemic redshifts are known and their S/N is the highest. We ran an MCMC fitting scheme similar to the one used for fitting the AG for the two isotropic geometries, TS and GW. The total number of free parameters for each model is five, since we fit an amplitude scale factor in addition to the three main parameters and the systemic redshift. Fig. 8 shows the best-fit GW models and residuals for the integrated spectra of A and B in the left and right panels, respectively. Remarkably, the systemic redshift order (z A > z B ) is correctly inferred by the models, despite the uniform priors and the fact that the line centroid of B is actually redder than A's. Moreover, the fitted redshifts z A = 2.9262 ± 0.0001 and z B = 2.9222 ± 0.0001 are only (50 ± 10) km s −1 and (−163 ± 10) km s −1 from the nebular redshift solution of A and B, respectively, consistent with the scatter of the different solutions ( § 3.1). The inferred expansion velocities and neutral column densities are V A exp = 224.3 km s −1 and V B exp = 187.2 km s −1 ; and log(N A H I /cm −2 ) = 18.9 and log(N B H I /cm −2 ) = 20.3. These two crucial parameters are constrained by the data, whereas the optical depth is not.
The best-fit velocity for galaxy A is higher than the best value for galaxy B. This already gives an interpretation to the line profile being narrower in A: since the medium is moving at higher velocities, the Lyα photons need to experience on average fewer scattering events to shift their frequency out of resonance, because the required shift is smaller. This effect is coupled with the lower column density of neutral hydrogen in A with respect to B. A lower density of neutral hydrogen atoms reduces the number of scattering events that a photon will experience before escaping. At a fixed Doppler shift and gas velocity, a Lyα photon has higher chances of escaping if there are fewer atoms along its path. We also fitted the diffuse-only spectrum in this fashion, that is, by having the systemic redshift as a free parameter. Interestingly, the best-fit redshift is z sys = 2.9222, extremely close to the best value given to B's model. The column density inferred by the model is also similar to B's (N diffuse H I ≈ 10 20.4 cm −2 ), but the expansion velocity of 198.3 km s −1 is in between V A exp and V B exp . These results are, however, more difficult to interpret since we have explicitly excluded the MUSE spaxels with UV continuum emission (a proxy for a high SFR), so there is no source of Lyα photons along the line of sight under the assumption of a pure scattering scenario. In this sense, the high column density should be regarded as an average integrated along random optical paths of the escaping photons, rather than along the line of sight between the observer and the source. This is the meaning that we will give to our subsequent results on the binned regions' spectra.
The models using the TS geometry yield similar results, but the χ 2 is larger in all the three fits and even though the z A > z B property is recovered, the predicted redshifts are less accurate. For this reason, plus the fact that the line shape does not vary dramatically across the halo at the MUSE resolution, we adopted the GW as our fiducial geometry to model the spectra from the individual bins.
Finally, given that our FLaREON GW model recovered the systemic redshifts reasonably well and given our lack of systemic redshifts for most of the binned regions, we also let that parameter vary freely in the fitting procedure for the binned regions. Here we used the 24-bin tessellation, since we required a higher S/N per spectrum. The source plane maps of the best-fit GW parameters are shown in Fig. 9. Again, the dust absorption optical depth is not constrained by the data, and hence we are unable to make inferences on this parameter.
The expansion velocities inferred by FLaREON under the GW geometry range from 80 to 360 km s −1 , with an uncertainty-weighted average of V exp = (211 ± 3) km s −1 and a standard deviation of (62.5 ± 3.0) km s −1 . Interestingly, Fig. 9 suggests some degree of spatial correlation between the expansion velocity values across the halo: Low-velocity bins are clustered along the north-south axis (vertical in the figure, and hereafter referred to as the "major axis" of the halo), including bins associated with the continuum of the central galaxies. Similarly, high-velocity bins are associated with diffuse-only regions at the outskirts of the halo, and they are found on both sides of the major axis and are also clustered in the same direction. These results suggest that on halo scales the medium is not expanding isotropically, but rather with a preferred direction. In other words, they might indicate that the gas is traveling faster in a direction perpendicular to the line connecting the two galaxies (e.g., due to decreased resistance from the environment).
In terms of column density, the models are distributed around log N H I /cm −2 = 20, with a 3σ-clipped stan- dard deviation of 0.4 dex. After discarding three outlier bins (∆ log N H I 1 dex), we found that the pathintegrated column density is remarkably uniform across the halo. Under the interpretation proposed above, this would mean that the actual column density declines as a function of the radial distance to the Lyα sources. Finally, in the right panel of Fig. 8 we plot the systemic redshifts predicted by FLaREON for each bin. If we interpret them as the zero velocity of the reference point inside of each GW model, its complex spatial structure and a range of values exceeding differences of 700 km s −1 would indicate that the outflow is modulated by complex underlying kinematics. Such complexity would naturally arise in a scenario where the two main galaxies are subject to dynamical interactions, for example, those observed in a merger. Further discussion of the potential impact of the interaction between A and B is made in § 4.1.

Low-Ionization Absorption Lines
Independent insight into the kinematics of the system can be obtained from the study of absorption features in the spectrum, although they only probe intervening gas along the line of sight to the central galaxies. As mentioned in § 3.1, previous studies of SGASJ1226 have obtained redshift solutions based on interstellar absorption lines that are ∼ 200 km s −1 bluer than the nebular emission based solutions (Koester et al. 2010;Wuyts et al. 2012). This result was later secured using the higher-resolution Magellan Echellette (MagE) spectrum of arc A.1 obtained as part of the MegaSaura Survey (Rigby et al. 2018). Moreover, Gazagnes et al. (2018) fitted Voigt profiles to the Lyman-series absorption lines (from Lyβ 1025.7Å to Ly6 930.8Å), obtaining a central velocity of −264 ± 21 km s −1 (−218 km s −1 if we convert to our nebular redshift solution).
Here, we complement these results by measuring the absorption velocity in our MUSE data, including galaxy B. But since the Lyman-series lines lie below the MUSE wavelength cutoff, we opted for the low-ionization Al II λ1670 line as an alternative tracer. We extracted spatially integrated spectra from the global apertures for A and B defined in § 3.3.1. Both spectra show a very similar profile that includes an asymmetric blue tail, and thus were modeled with an absorption AG profile (see Fig. 10 and Appendix B). We found v A = (−197 ± 4) km s −1 and v B = (−142 ± 13) km s −1 . Also, despite the line being narrower in A than in B (FWHM A = 300 ± 7 km s −1 versus FWHM B = 373 ± 15 km s −1 ), their EWs are consistent with each other, for EW A 0 = (5.8 ± 0.5)Å and EW B 0 = (5.6 ± 0.2)Å. The combination of blueshifted absorption lines and redshifted Lyα emission is very common among z ∼ 3 LBGs (e.g., Steidel et al. 2003), and it is often interpreted as a telltale signature of (spherical) galactic outflows (Verhamme et al. 2006).

Kinematics of SGASJ1226
A major goal of this paper is to understand the physical configuration of the SGASJ1226 system, by building a picture that includes both its spatial and kinematic properties. Firstly, the lens model reveals that the UV centroids of the two main galaxies of the system are separated only by 14.3 ± 1.4 kpc in projection (see § 3.4.1).  This fact, together with the small systemic velocity offset between the two galaxies (∆v ≈ 145 km s −1 ), suggests that they are gravitationally interacting. The interaction hypothesis is further motivated by the reconstructed source plane continuum morphology (contours in Fig. 6), which exhibits a very distorted galaxy B resembling the tidal tails or collisional rings seen in galaxy interactions at low redshift (e.g., Darg et al. 2010). The source plane morphology, however, should be interpreted with care, because the distorted appearance of galaxy B can also be due to (1) strong dust attenuation in some parts of the ISM that conceal the full shape of the galaxy, or to (2), geometrical artifacts that arise from the lens modeling. With the available data, we cannot distinguish whether the galaxies are experiencing a single visit flyby or the early stage of a major merger. On one hand, the lack of an ongoing starburst in both galaxies (their specific SFRs are consistent with the main sequence, Solimano et al. 2021) indicates that the interaction (if any) is not currently boosting the SFR, as it would be expected in some stages of a major merger (e.g., Hopkins et al. 2008). On the other hand, the extended Lyα emission implies a significant amount of gas between and around the pair of galaxies, thus favoring a merger scenario in which the interaction was able to strip gas away from the ISM and push it into the CGM (see, for example, Yajima et al. 2013). Also, we observed in § 3.4.3 that the fitted systemic redshift map exhibits a complex morphology and a large range of values, which provide hints of complex underlying kinematics similar to the expectation for galaxy mergers (e.g., Sparre et al. 2022). But in order to confirm or reject the merger activity as a driver of the gas motions, it would be necessary to obtain resolved spectroscopy of a nonresonant line, such as Hα or [C II] 158 µm, to trace the internal kinematics of the ISM where the effects of a merger would be more evident.
The global kinematics of the halo are clearly dominated by outflow motions, as revealed by the characteristic redshifted and red asymmetric Lyα profile. Under simple isotropic outflow geometries such as the TS and GW considered here, RT calculations always predict an enhanced red peak with a broadened red tail starting from relatively low expansion velocities. Extra evidence for outflows is seen in the absorption signature of lowionization metal lines (see § 3.5), which are thought to trace the same gas phase as Lyα. The signature is seen toward the UV continuum of both galaxies, A and B, but it remains unclear whether the outflows are launched independently by each galaxy or the outflow originates preferentially in one galaxy and the absorption appears in the second galaxy by a projection effect. The first case seems to be more likely, since the two galaxies have similar stellar masses and SFRs. Alternatively, the spatially coherent outflow signature could be linked to a large stream of receding gas stripped from the galaxies by the interactions described above.
Finally, we did not find any evidence of CGM-scale rotation, although such a signal would be severely smeared by the RT effects. In fact, rotation-like gradients in the Lyα peak shift velocity of LAHs are very rare, with only one candidate out of six in the sample presented by Leclercq et al. (2020), the only one among the few resolved LAHs in the literature (Claeyssens et al. 2019;Chen et al. 2021a).

Powering mechanism
We have thus far assumed that the extended Lyα emission of SGASJ1226 is exclusively explained by resonant scattering of photons produced in the central galaxies. In this section we explore alternative mechanisms that could also drive the observed properties of SGASJ1226.

Gravitational cooling
Galaxies need to accrete gas from the environment in order to sustain their growth over gigayear timescales (e.g., Bouché et al. 2010;Davé et al. 2012;Scoville et al. 2017;Tacconi et al. 2018;Walter et al. 2020). Simulations show that gas accretion can occur, for example, through cold streams of pristine gas from the IGM (e.g., Dekel et al. 2009;L'Huillier et al. 2012) or through inflows of previously ejected, metal-enriched gas in the CGM (Springel & Hernquist 2003;Oppenheimer et al. 2010;Brook et al. 2014;Anglés-Alcázar et al. 2017). In either case, collisional interactions in the gas are expected to transform its gravitational binding energy into Lyα radiation. This mechanism is known as "gravitational cooling," and it may contribute significantly to the development of extended LAHs (e.g., Haiman et al. 2000;Furlanetto et al. 2005). In this scenario, Lyα photons are created in situ, and thus their peak shift should trace the line-of-sight velocity of the gas. Then, if there is inflowing gas between the observer and the center of the halo, the Lyα spectrum emitted from it will have a strong blueshifted peak and a blue tail (e.g., Haiman et al. 2000;Dijkstra et al. 2006). Such profiles can also be produced by scattering-only scenarios (on a collapsing sphere of gas: e.g., Verhamme et al. 2006) provided that the intervening gas has a velocity in the opposite direction of the propagation of the Lyα photons. In other words, the presence of an enhanced blue peak is a signature of inflowing motion, but not necessarily of gravitational cooling.
In SGASJ1226, neither the integrated ( § 3.3.1) nor the resolved spectra ( § 3.4.2) exhibit significant emission blueward of the systemic velocity. Only a few tentative blue peaks can be seen in the Lyα profiles of some halo regions (see regions 0, 2, and 5 in Fig. C3). The lack of blueshifted emission cannot be entirely explained by intervening absorption in the IGM, because at z = 2.9 the average IGM transmission at −200 km s −1 from the Lyα rest-frame wavelength is approximately 85% (Laursen et al. 2011), high enough for a blue peak to be detectable. The caveat is that this particular line of sight can have a higher-than-average neutral gas column density. We therefore conclude that the presence of significant inflowing motion along the line of sight is unlikely but not completely ruled out.

Fluorescence
Cool hydrogen gas in the CGM can be momentarily ionized after being exposed to Lyman-continuum radiation escaping the inner parts of the galaxy (from either an AGN or a starburst region) or coming from the cosmic UV background (UVB). If the density is high enough, the atoms rapidly recombine and fall to the ground state, where a Lyα photon is emitted in a process called fluorescence. For this mechanism to produce extended Lyα emission, the gas needs to be distributed in high-density clumps with a low covering fraction, so the radiative transfer occurs preferentially in the surfaces of the clumps, resulting in a reduced number of scattering and absorption events with respect to the case of a homogeneous medium. In the approximation where scattering and dust absorption are negligible, this mechanism generates an intrinsic Lyα luminosity proportional to the production rate of ionizing photons. Here, we follow Valentino et al. (2016) to estimate the ionizing photon rate Q required to power the emission under the assumption that fluorescence is the only mechanism at play. From an observed, delensed f Lyα esc L Lyα = (6.2 ± 0.2) × 10 42 erg s −1 we obtained where η Lyα = 0.68 is the fraction of ionizing photons converted into Lyα (Spitzer 1978) and hν Lyα = 10.19 eV is the energy of a single Lyα photon. Since we cannot directly measure the production rate of ionizing photons (due to extreme ISM opacity at λ rest < 912Å) we need to rely on a longer-wavelength proxy such as the far UV luminosity (λ rest = 1500Å). But for a given L 1500Å the actual production rate of ionizing photons Q depends on the properties of the stellar population, particularly on the luminosity-weighted age and metallicity (Steidel et al. 2001;Smith et al. 2002). Fortunately, the stellar population synthesis analysis presented by Chisholm et al. (2019) provides estimates of the ionizing photon production efficiency, ξ ion ≡ Q/L 1500Å , for all galaxies in the MegaSaura sample including SGASJ1226's arc A.1. Therefore, assuming that ξ ion is uniform across galaxy A and is the same for galaxy B, we can solve for Q multiplying ξ ion by the total reddening-corrected UV luminosity In Table 5 of Chisholm et al. (2019), the best-fit Starburst99 model for A.1's MagE spectrum implies a photon production efficiency of log ξ ion = 12.74 ± 0.16 while the best-fit BPASS model favors a slightly higher value of log ξ ion = 13.04 ± 0.16. Then, the total rate of production of ionizing photons by the galaxies in SGASJ1226 is Q = (1.9 ± 0.7) × 10 54 s −1 for the Star-burst99 model and Q = (3.7 ± 1.5) × 10 54 s −1 for the BPASS model. Taken at face value, these results imply that photoionization from young stellar populations would only account for 20% − 60% of the total photon rate required to power the Lyα luminosity if we assumed f Lyα esc = 0.082, the value predicted by the Sobral & Matthee (2019) empirical relation based on the observed EW. Conversely, a larger escape fraction (between 0.15 and 0.30) would be needed to match the photon rates. Now, these calculations only considered young stars as the source of ionizing photons, but we cannot rule out the presence of an AGN that is obscured along the line of sight with our current data.
Also, we did not expect a significant contribution of metagalactic ionizing photons from the cosmic UVB, since the latest observational constraints imply that the UVB produces Lyα profiles at z ≈ 3 with peaks at the 2 × 10 −20 erg s −1 cm −2Å −1 arcsec −2 SB level (Gallego et al. 2021), at least 3 orders of magnitude fainter than the observed SB peak of SGASJ1226 on halo scalesalthough its contribution to the profile can become relevant at large radii, in the interface with the IGM.
In any case, we warn the reader that these ionizing photon budget considerations cannot constrain the contribution of fluorescence in the absence of an independent measure of the escape fraction and evidence of CGM clumpiness. The most direct test for fluorescence as a major contributor to the extended Lyα emission is the concomitant presence of extended Hα emission (Mas-Ribas et al. 2017). This is because under case B recombination (Osterbrock & Ferland 2006), the Lyα emissivity of the gas is 8.7 times the Hα emissivity. In other words, the same regions of photoionized gas should glow in Hα by a proportional amount, simultaneously solving the Lyα escape fraction and the question of in situ Lyα production. Promisingly, SGASJ1226 has been selected as a NIRSpec IFU target for the JWST Early Release Science "TEMPLATES" program (Rigby et al. 2017) and thus resolved Hα observations of the arc A.1 will become available. While the FoV of the instrument will not cover the whole LAH, it will certainly tell us if the Hα extends beyond the ISM.

SGASJ1226-LAE
In § 3.3 we reported the discovery of a continuum counterpart to a local maximum of Lyα SB in the MUSE NB image labeled SGASJ1226-LAE. In this section we argue that SGASJ1226-LAE is a satellite of the main system composed of galaxies A and B, rather than another UV clump in the ISM of galaxy B. According to our best-fit lens model, this source has an average magnification of µ = 14 and lies ∼ 3 kpc away from the UV centroid of galaxy B (see Fig. 6). This is approximately twice the exponential UV scale length of the con-tinuum of galaxy B, further away than any of the UV clumps identifiable in that galaxy. Also, the candidate counterpart appears spatially resolved in the F606W image, with an exponential scale length of 250(14/µ) 1 2 pc, unlike the clumps of galaxy B, which are all pointlike. SGASJ1226-LAE has an intrinsic Lyα luminosity of (1.0 ± 0.2) × 10 41 erg s −1 , representing the 2% of the parent LAH luminosity. The F606W photometry implies an absolute UV magnitude of M 1500 ≈ −16.7, which lies at the faint end of the luminosity function of LAEs at this redshift (Ouchi et al. 2008;Kusakabe et al. 2020). In fact, the UV continuum at this luminosity is extremely difficult to detect, and only a handful of objects have been robustly detected in lensed fields (Claeyssens et al. 2022) or in deep stacks from the UDF (Maseda et al. 2018). Although with very low significance, the UV slope of this source is steeper than the mean slope of galaxies A and B, which is an expected property of strong LAEs due to their young ages (e.g., Nakajima et al. 2012;Hagen et al. 2014) and low dust attenuation (e.g., Ono et al. 2010;Stark et al. 2010;Kojima et al. 2017).
Theoretical models and simulations predict the presence of several satellites populating the outer parts of LAHs, and some authors propose that they contribute a significant fraction of the Lyα SB at large radii (r 0.25R vir ) (Mas-Ribas et al. 2017;Mitchell et al. 2021). It is thus plausible that SGASJ1226-LAE is indeed a satellite of the SGASJ1226 system, made detectable by the chance alignment of the lensing caustic boosting its flux above the background.
Unfortunately, we did not detect additional lines in the MUSE spectrum of SGASJ1226-LAE due to strong contamination from the foreground galaxy light at λ obs 4800Å. Without a systemic redshift for this galaxy, interpreting the line profile becomes even more difficult. Nevertheless, the line profile shows a broad red tail, suggesting that Lyα photons coming out of SGASJ1226-LAE are also scattered in the CGM.

CONCLUSION
We have analyzed the spatial and spectral properties of the diffuse Lyα emission associated with a pair of lensed LBGs at z ≈ 3. The remarkable brightness and extension of this system, together with high-quality MUSE observations, allowed us to probe in detail its physical nature.
The system is composed of two main-sequence galaxies (labeled A and B) of similar stellar mass (≈ 10 10 M ) and SFRs (≈ 10 M yr −1 ) that are separated by less than 15 kpc when projected in the source plane, and by 145 km s −1 in the velocity space, suggesting an interact-ing pair. The galaxies are associated with a single LAH of L Lyα = (6.2 ± 1.3) × 10 42 erg s −1 which we decomposed into two Sérsic profiles in the source plane with the largest component having a circularized half-light radius of 19.4 kpc. Despite its apparent Lyα brightness, the whole system has a rest-frame Lyα EW of only (17.0 ± 2.7)Å. Globally, the Lyα line exhibits a redshifted peak with an asymmetric red tail, typical of CGM-scale outflows.
We found significant ±200 km s −1 spatial variations of the line FWHM and peak shift velocity across the halo. The lowest values of the FWHM and peak velocity shift were preferentially found on top of the central galaxies. We also recovered a correlation between these two spectral properties, in line with recent results of resolved LAHs.
We divided the source plane emission into 24 spatial bins and fitted them with radiative transfer models in isotropic galactic wind geometries. At an average expansion velocity of (211 ± 3) km s −1 and standard deviation of (62.5 ± 3.0) km s −1 we found tentative evidence of structured gradients along the minor axis, which suggests the outflow has a preferred direction. Also, the best-fit models imply that a typical Lyα photon in the halo encounters a neutral column density of ∼ 10 20 cm −2 integrated along its path. The existence of the outflow is further confirmed by the presence of blueshifted, asymmetric absorption lines of low-ionization metal species in the UV spectrum of the central galaxies. In particular, we measured the Al II λ1670 absorption central velocity at (−197 ± 4) km s −1 and (−142 ± 2) km s −1 in A and B, respectively, in broad agreement with the velocities inferred from the wind models. Finally, the recovered systemic redshifts for the different source plane regions show a complex structure that could be also explained by interaction processes between the galaxies.
We explored different mechanisms that could be producing the extended emission besides the resonant scattering of photons produced in the central galaxies. Lyα production in situ by gravitational cooling is disfavored since we found no indication of infalling gas mo-tion (e.g., a dominant blue peak) assuming an average IGM transmission. However, a major contribution of fluorescent radiation is allowed by energy budget arguments but it otherwise remains unconstrained due to uncertainties in the Lyα escape fraction and the clumpiness of the CGM. Upcoming Hα observations with JWST will be key to establishing the contribution of extended fluorescent radiation in SGASJ1226.
The boost in spatial resolution provided by the lensing effect allowed us to detect the continuum counterpart of a faint (M 1500 ≈ −16.7) satellite that is a strong LAE (EW 0 = (104 ± 19)Å) and contributes 2% of the total Lyα luminosity. Moreover, this source is resolved with an exponential scale length of ≈ 250 pc. This is one of the few observational hints that such satellites do indeed exist, and contribute to the Lyα luminosity.
We thank the anonymous referee for helpful suggestions that which have contributed to improving the quality of this manuscript. We also warmly thank Lucia Guaita for her insightful comments and discussion. This  Another important observable is the rest-frame Lyα EW, W Lyα , since it has been shown to correlate very strongly with the Lyα escape fraction f Lyα esc of a galaxy (e.g., Harikane et al. 2018;Sobral & Matthee 2019), defined as the ratio between the observed and intrinsic Lyα luminosies. However, measuring the EW is subject to some complexities when the line profile is composed of both emission and absorption components (e.g., Kornei et al. 2010;Erb et al. 2019), as is the case of SGASJ1226. Here, we calculate the net or total EW dividing the Lyα flux by the expected continuum level at λ rest = 1215.67Å, which was extrapolated from a powerlaw fit to the continuum at λ rest 1270Å. Using the image plane apertures for A and B, we extracted the full MUSE spectra and fitted a power law described by f λ ∝ λ β to all line-free channels at λ rest 1270Å. We used uniform priors on β and the normalization factor. This resulted in slopes β of −1.18 ± 0.15 for galaxy A and −0.79 ± 0.21 for galaxy B. We used these values to extrapolate the demagnified F606W magnitudes (λ pivot = 5921Å) to the redshift Lyα wavelength (4771Å) continuum, finding a total (A+B) intrinsic flux level of f 4771 = (1.03 ± 0.16) × 10 −18 erg s −1 cm −2Å −1 . Then, we computed the EW by dividing the total delensed Lyα flux of the LAH, F Lyα = (6.9 ± 0.2) × 10 −17 erg s −1 cm −2 (see § 3.4.2), by the underlying continuum level obtained above. After propagating all uncertainties, this ratio yields W Lyα = (66 ± 10)(1 + z) −1Å = (17.0 ± 2.7)Å. Such a value is typical of LBGs (e.g., Steidel et al. 2003) and, according to the empirically calibrated relation of Sobral & Matthee (2019), it implies a global f Lyα esc = 0.082 ± 0.018. However, with the available data it is also possible to investigate spatial variations of the EW across the arcs. The distribution of the EW can inform us about the homogeneity of the intervening gas. For example if the EW is uniform across the source, it could mean that the Lyα signal is processed by an approximately homogeneous slab of gas. If the EW distribution is clumpy or has gradients, a more complex geometry can be in place. We constructed a map of the EW as the ratio between the Lyα NB and the extrapolated UV flux density at λ Lyα . To avoid discrepancies in the spatial resolution at different wavelengths, we estimated the UV slope in each spaxel by taking the ratio between PSF-homogenized images at λ obs = 5300Å and λ obs = 7060Å (rest-frame 1350Å and 1800Å, respectively) before extrapolating to λ Lyα . The images were obtained by averaging over the spectral axis on a ≈ 200Å window centered at these two wavelengths and masking the channels with line absorption or emission. We created models of the MUSE-AO PSF by fitting a circular Moffat 2D profile to a single star near the center of the field and used them to convolve all images to a common PSF with FWHM=0.72 . We computed the UV slope array from the ratio between the blue and red continuum images and then used it to extrapolate the continuum flux density to the Lyα wavelength. The EW map was finally constructed from the resulting continuum image at λ Lyα and the NB image (see Fig. A1). This operation was restricted to the same continuum apertures defined in § 3.3.1.
We observe that W Lyα is mostly uniform toward galaxy A, with an average value of 4Å, whereas galaxy B has values that range from 3Å to 42Å. The increased W Lyα toward the northern and eastern edges of arc B can be explained by the large spatial offset between the continuum and the brightest knots of Lyα emission. As we mentioned above ( § 3.3), the UV arc is offset by 0. 6 from the Lyα arc, and there is also a bright knot of Lyα emission to the east of the arc that has little overlap with the continuum. In Fig. A1 we also provide the correspondence between EW and f Lyα esc as calibrated by Sobral & Matthee (2019), but the interpretation of the EW map as an f Lyα esc map would require the relation to hold also for resolved regions.

B. AG FITTING
The AG profile as parameterized by Shibuya et al. (2014) serves as a tool to measure basic properties of typical Lyα spectra, namely the amplitude, peak shift velocity, FWHM, and an asymmetry parameter that quantifies the skewness of the curve. The use of the AG profile is becoming increasingly common in the resolved LAH literature (Claeyssens et al. 2019;Leclercq et al. 2020;Claeyssens et al. 2022). In this paper we fit Lyα spectra either from the integrated apertures ( § 3.3.1), from individual resolved regions ( § 3.4.2) or from absorption lines ( § 3.5) using a standard MCMC posterior sampling scheme. For each fit, we imposed a Gaussian likelihood for the residuals, which were computed from the difference between the data and the AG profile convolved with the LSF. We sampled the posterior distribution using 32 walkers and 2000 steps within the emcee library (Foreman-Mackey et al. 2013) to obtain robust estimates of the parameter uncertainties and covariances. We set uniform priors on each of the five parameters (the four named above plus an amplitude parameter) and let them vary freely over their domain ranges. During optimization, we also kept track of the integrated flux of the model. The value for the peak shift velocity depends on the systemic redshift input by the user. An example of an AG fit is displayed in Fig. B2.  UV surface brightness (sigma) Figure D4. Upper panel: MUSE broadband image at λrest ∼ 1600Å. Lower panel: Same as above but with the best-fit model for galaxies A and B subtracted (after convolution with the MUSE PSF). In both panels the SB units are normalized to the background RMS and the dashed curve indicates the "diffuse halo" aperture (see § 3.3.1).