Deciphering Lyman-α Emission Deep into the Epoch of Reionisation

During the epoch of reionisation the first galaxies were enshrouded in pristine neutral gas, with one of the bright-est emission lines in star-forming galaxies, Lyman-α (Ly α ), expected to remain undetected until the Universe became ionised. Providing an explanation for the surprising detection of Ly α in these early galaxies is a major challenge for extra-galactic studies. Recent JWST observations have reignited the debate on whether residence in an overdensity of galaxies is a sufficient and necessary condition for Ly α to escape. Here, we take unique advantage of both high-resolution and high-sensitivity images from the JWST instrument NIRCam to reveal that all galaxies in a sample of z > 7 Ly α emitters have close companions . We exploit novel on-the-fly radiative transfer magnetohydrodynamical simulations with cosmic ray feedback to show that galaxies with frequent mergers have very bursty star formation which drives episodes of high intrinsic Ly α emission and facilitates the escape of Ly α photons along channels cleared of neutral gas. We conclude that the rapid build up of stellar mass through mergers presents a compelling solution to the long-standing puzzle of the detection of Ly α emission deep into the epoch of reionisation. Young, vigorously star-forming galaxies have been identified in the very early Universe 1,2,3 . These galaxies should be excellent sources of Lyman-α emission (Ly α , λ =1215.67˚A) – the intrinsically brightest emission line 4 – which stems from the recombination of hydrogen that has been ionised by their young stellar populations. However, deep in the epoch of reionisation galaxies are expected to be exceptionally gas-rich such that their stellar nurseries are enshrouded in copious amounts of neutral hydrogen, which leads to extreme damped absorption of Ly α 5 . Furthermore, the intergalactic medium (IGM) is increasingly neutral as we probe to higher redshift 6,7 and this

Young, vigorously star-forming galaxies have been identified in the very early Universe 1, 2, 3 .These galaxies should be excellent sources of Lyman-α emission (Lyα, λ=1215.67Å) -the intrinsically brightest emission line 4 -which stems from the recombination of hydrogen that has been ionised by their young stellar populations.However, deep in the epoch of reionisation galaxies are expected to be exceptionally gas-rich such that their stellar nurseries are enshrouded in copious amounts of neutral hydrogen, which leads to extreme damped absorption of Lyα 5 .Furthermore, the intergalactic medium (IGM) is increasingly neutral as we probe to higher redshift 6,7 and this neutral gas is expected to resonantly scatter Lyα emission.Hence, due to 'local' attenuation by a gasrich interstellar medium (ISM) and scattering by a neutral IGM, Lyα emission should only be detectable towards the end of the reionisation era, about one billion years after the Big Bang 4,8,9,10 .
While the decreasing observability of Lyα emission with increasing redshift has been repeatedly claimed to be observed 11,12,13 , this picture has been challenged by the occasional, surprising detection of Lyα emission in several galaxies deep in the reionisation era 14,15 .It has been suggested that Lyα can escape through the neutral IGM if the galaxies reside in sufficiently large ionised bubbles embedded within the neutral IGM, driven either by active galactic nuclei (AGN) 16,17,18 or by an enhanced radiation field produced by an overdensity of associated objects 19,20,21,22,23,24,25 .However, a solution to the escape of Lyα emission through the ISM and CGM of what are expected to be very gas-rich galaxies remains elusive.Before the advent of JWST, the sensitivity and resolution of imaging instruments meant that studies of Lyα emitters (LAEs) at high redshift were spatially unresolved.Hence, it was not possible to probe the physical processes that could explain the escape of Lyα emission from the ISM.
To address this crucial issue we study a sample of nine galaxies that have been spectroscopically confirmed with the detection of Lyα emission at z > 7 with high-resolution spectrographs and which have been observed with publicly-available JWST Near Infrared Camera (NIRCam) 26 imaging (the properties of which are reported in Table 1).These galaxies fall within the GOODS-North, GOODS-South, EGS and COSMOS fields, and have been observed as part of five different programs: PRIMER (PI: Dunlop), FRESCO (PI: Oesch) 27 , CEERS (PI: Finkelstein) 28 , JADES (PI: Eisenstein) 29 and DDT program 4426 (providing NIRSpec IFU observations of GN-z11, PI: Maiolino) 30 .Six of these nine Lyα-emitting galaxies are known to lie within overdensities 23,31,21,24,25 , three of which also likely host accreting black holes 16,17,18 .Moreover, it has been shown that the remaining galaxies are incapable of alone blowing a large enough ionised bubble to facilitate the escape of Lyα through the neutral IGM, suggesting that an underlying population of faint galaxies must be surrounding these LAEs 25 .
However, the presence of spectroscopically confirmed galaxies within known ionised bubbles that do not show Lyα emission 23,13,25 indicates that there must be further local processes at play driving Lyα emission deep into the epoch of reionisation.With this in mind, we use the unparalleled sensitivity and resolution of JWST /NIRCam in order to accurately determine the properties of these nine LAEs, with equivalent width (EW) greater than 10 Å, at z > 7.
Remarkably, NIRCam images of these galaxies reveal the presence of multiple components to all of these LAEs, as shown in Figure 1.Given the tight redshift constraints provided by strong [OIII] 4959,5007 emission falling into medium band filters, we estimate the probability of these being unrelated, high-redshift objects and find this is always below 14% for each of our systems and is often just 1% (see Methods for further details).We estimate the redshift of each candidate companion by carefully extracting the Spectral Energy Distribution (SED) of both the main component and its companions in Extended Figure 2 and 3. We then fit each SED by assuming several parametric star formation histories (SFH -burst, constant, delayed, exponential), allowing a large range for all parameters 32 .In all cases, we deduce the photometric redshifts of the companions to be consistent with that of their spectroscopically-confirmed main component.Cutouts of JWST images (as different surveys make use of different filters, the filter within which the companion is most clearly present is shown: F182M for GSDY and z7-GSD-3811, F200W for EGSY8p68, F115W for CEERS-1027, JADES-GS-z7-LA and COSY and F150W for JADES-GS+53.15682-27.7671and z7-13433) showing the multiple components of each system.All images are smoothed with a Gaussian of FWHM ∼ 0.7 kpc to enhance the visibility of the companion, except for EGSY8p68 where the three components become unresolved if we smooth the images.We include the unsmoothed images as sub-panels in the top right of each panel in order to emphasise that the companion is always visible before smoothing.We use different colorbar ranges for each system in order to make the companion galaxy well defined, hence some  To spectroscopically confirm their identity as companions we exploit existing XSHOOTER/VLT 41 , MOSFIRE/Keck 42,43 , Near-Infrared Spectrograph (NIRSpec) on-board JWST 44 and NIRCam wide-field slitless spectroscopy (WFSS) observations 45 of these galaxies.For three systems in our sample we have no spectroscopic information on their companions given the limitations of seeing at ground-based telescopes and the small aperture size of JWST/NIR-Spec observations.For the remaining six systems in our sample we are able to spectroscopically confirm the companion galaxy as it falls within the FoV of the observations of the main LAE.For two companions in our sample we see tentative (∼ 3σ) evidence of Lyα emission in Extended Figure 4 (one of which has recently been confirmed by NIRSpec IFU observations 40 .For a further two companions we observe the [OIII] 5007 emission line at > 8σ (as seen in Extended Figure 5 with fluxes reported in Extended Table 1), and NIRSpec IFU observations spectroscopically confirm the companions of the final two systems 30,39 and [Carniani et al. in prep.].We therefore conclude that for all systems for which we have resolved spectroscopic observations of the companion (∼ 67%), we spectroscopically confirm its nature as a physical companion (see the Methods for further details of the spectrosocopic analysis).
To confirm that the presence of a close companion is the primary factor governing the visibility of Lyα, we determine the fraction of companions seen in a mass-matched sample of z > 7 galaxies with high-resolution spectroscopic data for which Lyα is not detected 13 .We find that 43% of these galaxies have photometric-candidate companions within 5 ′′ of the central galaxy which is consistent with the companion fraction determined by a more comprehensive study [Puskás et al., in prep].The lower fraction of close companions among samples that are not selected for Lyα emission is evidence that the 100% rate of companions for our sample of LAEs is atypical of z > 7 galaxies.
To reinforce the observational evidence supporting the idea that ongoing interactions drive an increase in the detectability of Lyα emission during the epoch of reionisation, we explore comparable galaxy mergers using the novel Azahar simulation suite (Martin-Alvarez et al., in prep.).This simulation suite was performed before we obtained the NIRCam data, but remarkably, we find many simulated galaxies that match the photometry and the spatial geometry of our objects.Azahar is a cosmological, high-resolution, zoom-in simulation which employs a magnetohydrodynamical solver together with state-of-theart cosmic ray feedback (see the pathfinder Pandora project 46 and Methods).Most importantly for this work, Azahar also features on-the-fly radiative transfer 47,48,49 , capable of self-consistently reproducing reionisation while fully modelling the ISM of galaxies (with a maximum spatial resolution of 10 pc), and crucially resolving the propagation and escape of ionising radiation on the ISM scales 50,51 .We post-process our simulations with the publicly available RASCAS code 52 to account for both the production and resonant scattering of Lyα photons as well as scattering or absorption of Lyα photons by dust 53 .d) hydrogen gas and ionising radiation properties (blue: HI density, grey-black: HII density, yellow-white: F150W intensity, purple: LyC radiation energy density, orange: HeII ionising radiation energy density).The appearance of the observed system is very well reproduced by an ongoing merger of three galaxies, with an extended escaping LyC radiation field.
Due to a substantial overdensity of galaxies in our simulated volume, the main progenitor undergoes repeated mergers with multiple other galaxies brought in by the cosmic filaments, often involving several companions or mergers in rapid succession as has been observed in 54 .To highlight one such merger occurring at z ∼ 7.3, the four rightmost panels of Figure 2 show different projections of Azahar.This particular interaction features three merging galaxies that will constitute the main progenitor of the spiral galaxy formed by z ∼ 1.These are very good analogues to the EGSY8p68 observations shown on the left (modelling of the observations with Galfit in Extended Figure 1 confirms the presence of three components), with very similar individual galaxy sizes, mutual distances and merger configuration.Specifically, the three simulated galaxies have stellar masses of M * ,1 = 2 • 10 8 M ⊙ , M * ,2 = 3 • 10 8 M ⊙ , and The SFR for each of the three galaxies and the final system (black line).Dot symbols show mergers with companions (filled markers) and with other secondary systems (open markers); the intrinsic Lyα luminosity of the three galaxies and the entire system; the fraction of escaped Lyα luminosity filtered for pixels with fluxes higher than a given limit (median (dark green curve), quartiles (darker band), and min-max of the distribution across 12 lines-of-sight (lighter band)).Distance between the main progenitor and its companions is shown as well.The grey vertical lines in Rows 3-5 correspond to the redshifts shown in Rows 1-2.The rightmost panels show key physical relations for the observations (black) and the simulations (pink).
M * ,3 = 8•10 7 M ⊙ at this redshift, which are comparable to the stellar masses of the EGSY8p68 galaxies (see Table 1).Furthermore, the simulated SFRs vary between 1 − 10 M ⊙ yr −1 (shown in Figure 3) and are consistent with the SFRs of the EGSY8p68 galaxies (7.62 M ⊙ yr −1 , 2.22 M ⊙ yr −1 , 4.26 M ⊙ yr −1 ).We also consider the total Lyα luminosity of the simulated system and find that, while it varies, its average value is ∼ 10 43 erg/s (shown in Figure 3), in agreement with the value for EGSY8p68, reported in Table 1.Finally we note that simulated volume hosts a substantial overdensity of galaxies, which is consistent with the local environments of the observed LAEs in our sample (as discussed earlier), including EGSY8p68 21 .We thus conclude that the simulated system matches all key observed photometric and spectrosopic properties of EGSY8p68 (see also Methods).
The observations shown on the left clearly indicate the superior resolution and sensitivity of JWST (top panel) over HST (bottom panel).The tri-component nature of EGSY8p68 is entirely unresolved in existing HST observations, but clearly identifiable in the F150W NIRCam imaging.Panel (a) shows the large-scale view (150 kpc across) of the three filaments encompassing the merging system, where large amounts of HI gas (blue) along the filaments are feeding the star formation in these systems, resulting in significant NIR-Cam F150W emission (yellow-white).This active star formation is driving ionised hydrogen bubbles (grey) away from the filaments and into the low density regions.Panel (b) shows the simulated JWST -like 150W observation (with dust extinction modelled as an absorption screen along the line-of-sight), where we select the line-of-sight to illustrate the resemblance with EGSY8p68.Note that the simulated merging system is displayed at a slightly earlier stage of the merger with respect to EGSY8p68, with our simulated galaxies approaching the likely physical separation of the observed galaxies at z ∼ 7. Interestingly, the mock F150W emission reveals compact galactic cores analogous to these JWST observations, with diffuse emission (blending with the background) emerging from the extended, low density stellar discs.Panel (c) shows a synthetic colour-composite observation using the JWST filters at the full resolution of our simulation.While the two companion galaxies are more diffuse, with stars actively forming in their discs, the main galaxy has a much more compact core due to the preceding rapid growth through repeated mergers.Panel (d ) is a colour composite of various simulated properties relevant for Lyα detectability.The emission from LyC ionising photons (purple) and HeII ionising photons (orange) is distributed on a scale of a few (physical) kpc.Importantly, there is a clear separation of the stellar emission and the HI gas during the merger event, with the ionising radiation escaping from star forming regions.
To provide a quantitative confirmation and in-depth physical understanding of this effect, we explore in Figure 3 the time evolution of the three merging galaxies in detail.In the top row the colours encode the same quantities as in panel (d ) of Figure 2, while in the second row we show the intrinsic Lyα emission (blue) and resonantly scattered Lyα emission (orange), produced by post-processing our simulations with RASCAS.The first column shows a larger scale view of the main progenitor (galaxy 2) and its approaching companions at z ∼ 8.1.The second column shows a close up views of the main progenitor at z ∼ 8.1 undergoing a minor merger, with a complex topology of neutral gas and channels through which Lyα photons are escaping.The third column shows the ongoing merger of our galaxies 1, 2 and 3 at z ∼ 7.3 (with a different orientation with respect to Figure 2), where bursty star formation drives ionised channels through the ISM, with a significant amount of gas also tidally stripped from the galaxies.The fourth column shows a post-merger view of the newly formed disc-dominated galaxy at z ∼ 6.1.Here, the HI is confined to the merger remnant galaxy, and the radiation can only escape through the channels opened by collimated, cosmic ray-driven bi-conical outflows.Focusing now on the Lyα radiation, all three galaxies have significant intrinsic Lyα emission due to their ongoing star formation.However, at this particular viewing angle, due to absorption by dust, Lyα scattered emission is only present around merging galaxy 2 at z ∼ 8.1 and around galaxies 1 and 3 at z ∼ 7.3.At z ∼ 7.3, the remaining emission, not absorbed by dust, is very diffuse on large scales and hence not observable.
To understand what drives the observable Lyα emission, in the third row we show the SFR evolution for our three merging galaxies, together with the SFR history of the merger remnant at z ∼ 6.1.SFRs of all three galaxies undergo repeated cycles of bursts, driven by numerous mergers ('major' mergers [filled dots] and 'minor' mergers [empty dots]) and a high SFR 'plateau' during the final merger episode from z ∼ 7 to 6.This is mirrored in the evolution of the intrinsic Lyα emission of each galaxy as well as the total emission integrated across our three simulated galaxies, as shown in the fourth row, where there is a clear correspondence between the SFR peaks and peaks of intrinsic Lyα.Our simulations hence demonstrate that high SFR triggered by gas-rich mergers results in brighter intrinsic Lyα emission, which enhances the probability of these systems being detected.
To constrain this more robustly, we select an observable Lyα flux threshold of 10 −18 erg s −1 cm −2 consistent with our observations, and compute the ratio of scattered to intrinsic Lyα emission over the spatial extent of our mock observations of the galaxies.This is shown in the fifth row, where the dark green curve denotes the median of the ratio of scattered to intrinsic Lyα flux and the shaded bands (quartiles and minimum-maximum values) show its distribution computed over 12 different lines-of-sight.As the fraction of escaping Lyα emission has a number of high peaks we conclude that the Lyα emission should be detectable for a number of the redshifts studied here.Interestingly, during close galactic encounters (as shown by the distance trajectories of three main galaxies shown in the same panel by coloured lines) and mergers, significant gas tidal stripping and bursty star formation feedback leads to the opening of more low-column density sight-lines, hence facilitating the escape of Lyα radiation.This effect can be clearly seen when comparing the most important merger instances to the peaks (high-tails) in the distribution of Lyα escape fraction.Our findings regarding higher intrinsic Lyα emission and enhanced Lyα escape fraction in mergers are robust to the assumed dust model (see Methods and Extended Figure 6).It is worth noting that inclusion of AGN feedback in our simulations would only reinforce our findings, as expected black hole fuelling during mergers would power AGN 55 .
The rightmost panels of Figure 3 show various physical relations for the observations (black points) and the simulations (pink 'violin' symbols, denoting the distribution of Lyα escape fraction over 12 lines-of-sight), using the same flux filter as above, demonstrating a very good correspondence between our simulations and our observed LAEs.Interestingly, both for observations and simulations we do not find any clear correlations between the Lyα escape fraction and any other galaxy properties (distance, SFR, stellar mass), further corroborating our finding that 'favourable' lines-of-sights with evacuated neutral hydrogen channels lead to directional Lyα photon leakage.
In this paper, we introduce a new interpretation to explain the unexpected detection of Lyα in z ≥ 7 galaxies in an epoch where the IGM is mostly neutral and gas-rich, early galaxies are heavily 'locally' Lyα damped.Combining the high-resolution and high-sensitivity of JWST data with state-of-the-art radiative transfer on-the-fly magnetohydrodynamics simulations, we demonstrate that three ingredients are key to making Lyα emission detectable from our sample of galaxies deep in the epoch of reionisation: galactic mergers driving high intrinsic Lyα emission in the host galaxy; a 'favourable' line-of-sight cleared of local neutral hydrogen in the host galaxy by tidal interactions with companions and by star-formation feedback; a sufficiently large ionised bubble facilitating the escape of Lyα emission through the IGM.

Methods Data
The emergence of JWST with its significant technical advancements over both ground-and space-based telescopes, has the potential to reveal new hints as to the driving forces behind the observed Lyα emission discussed in this paper.We therefore decided to analyse all spectroscopically-confirmed LAEs at z > 7, with JWST /NIRCam imaging.This search originally revealed eleven such LAEs, however, further constraints on the SED of these galaxies from NIRCam imaging reveals the Lyα emission line detected in two of these LAEs is in fact not Lyα (discussed in more detail below).Therefore, the sample becomes nine LAEs (seen in Table 1), eight of which have companions that are identifiable in publicly available NIRCam imaging and GN-z11 that has spectrosocpically confirmed companions that do not show a UV continuum 30,39 .These galaxies are found in the GOODS-North and GOODS-South, EGS and COSMOS fields which have been observed with NIRCam by FRESCO (PI: Oesch, ID: 1895) and JADES (PI: Eisenstein, ID: 1180), CEERS (PI: Finkelstein, ID: 1345) and PRIMER (PI: Dunlop, ID: 1837), respectively.
The FRESCO survey's primary aim is to obtain NIRCam Long-Wavelength (LW) Grism spectra over the F444W filter for galaxies across GOODS-S and GOODS-N.However, the simultaneous imaging in the short-wavelength channel provides imaging in both F182M and F210M filters, and direct imaging also obtained in the F444W filter additionally provides us with a total of 3 filters for our galaxies in GOODS-S/GOODS-N.FRESCO's LW WFSS brings significant advantages as the F444W filter covers the wavelength range of 3.9 − 5µm, therefore observing both [OIII] 4959,5007 and Hβ at z ∼ 7 − 9.These emission lines ordinarily contaminate imaging in the F444W filter and hence constraints on the fluxes of these emission lines allow us to reduce the uncertainty in galaxy properties such as stellar mass, star-formation rate and age.We use these constraints on the emission lines present in F444W to create a mock F430M filter.Furthermore, when we fit the SED (discussed further below) we check the fluxes of the [OIII] 4959,5007 and Hβ emission lines in the best-fit SED model and confirm that these are consistent with the measured fluxes from the WFSS spectrum.

EGSY8p68 z7_20237
Extended Figure 1.Example of Galfit modelling for the EGSY8p68 system for which the separation between each component is too small for the photometry to be extracted by SEXtractor.

Data Reduction
Throughout our analysis we use the final reduced data products created by the PRIMER team for the COSMOS field, the publicly available reduced data products from the CEERS team for the EGS field, the publicly available reduced data products from the JADES team for the deep HST region of GOODS-S and our own reduction of the GOODS-S FRESCO imaging for regions outside of the publicly available JADES imaging.In order to produce the GOODS-S images, we take all available NIRCam imaging data of GOODS-S, available from the FRESCO survey, and reduce them using version 1.8.5 of the JWST pipeline.We first download the uncalibrated files from the MAST archive and follow the steps of the JWST pipeline: (Stage 1) detector-level corrections and ramp fitting, (Stage 2) instrument-level and observing-mode corrections producing fully calibrated exposures, (Stage 3) producing the final mosaic.We take additional care after Stage 1 to ensure the removal of horizontal and vertical striping that can be present in the rate files.
We extract the photometry of our objects using SEXtractor 58 on PSF-matched images extracting all objects with more than 4 pixels (DETECT MINAREA 4) detected at more than 1.5σ (DETECT THRESH and ANALYSIS THRESH 1.5).We used large deblending parameters to allow efficient separation of each component (DEBLEND NTHRESH 64 and DEBLEND MINCOUNT 0.000001), however, the small separation between the components of the EGSY8p68 and JADES-GS-z7-LA systems makes the extraction of each component's photometry impossible with SEXtractor.For JADES-GS-z7-LA we employ the photometry of the LAE and its companion reported in 35 while for EGSY8p68 we model the shape of each component using Galfit 59 assuming a Sersic profile and NIRCam PSFs, simulated using WebbPSF.We first run Galfit on one of the NIRCam images where all the components are clearly resolved (F115W), and we then fix the resultant position, Sersic index, axis ratio and angle for each component when extracting the photometry for the remaining images.Extended Figure 1 shows the modelling of the components for this system.
As discussed earlier, it is desirable to spectroscopically confirm as many companion galaxies as possible.The use of the Mutli-Shutter Array with NIRSpec to observe three of the systems (EGSY8p68, CEERS-1027 and JADES-GS+53.15682-27.76716)means that, given the small shutter size, the companion is not coincident with the shutter position.Instead, we make use of all existing ground-based observations of our sample of LAEs, of which four of the sample have been observed by MOSFIRE (z7-13433 -ID: N199, PI: Jung; EGSY8p68 -ID: C228M, PI: Zitrin; GSDY -ID: U069, PI: Treu; z7-GSD-3811 -ID: C182M, PI: Scoville) and one with X-shooter (COSY -ID: 097.A-0043, PI: Ellis).
We reduce the spectroscopic data using the standard MOSFIRE data reduction pipeline and EsoReflex.These pipelines perform wavelength callibration, sky subtraction and flat fielding as well as further standard data reduction steps in order to produce two-dimensional (2D) spectra.The resultant 2D spectra span the length of the original slit used for the observations and have spectral and spatial resolutions of R ∼ 3380, 0.1798 ′′ /px and R ∼ 8900, 0.158 ′′ /px for MOSFIRE and X-shooter respectively, but are limited by the seeing of each observation which can be as poor as ∼ 0.9 ′′ .Two of our sample are located in the GOODS-S field (z7-GSD-3811 and GSDY).This field was the target of WFSS observations by the NIRCam Long-wavelength (LW) Grism as part of the FRESCO observing program.FRESCO's LW WFSS covers the wavelength range of 3.9 − 5µm using the F444W filter, therefore observing both the [OIII] 4959,5007 and Hβ emission lines for z ∼ 7 − 9 galaxies, which are expected to be bright in high-redshift galaxies.These lines have previously been observed in z > 7 galaxies using the publicly available FRESCO data 60,27 .Moreover, the significantly superior spatial resolution of NIRCam LW WFSS (60 milli-arcesconds) allows us to resolve emission lines where MOSFIRE and X-shooter would fail to do so.
We reduced the Grism data following the steps described in 61 including flat fielding, background subtraction, 1/f noise subtraction and WCS assignment.We also perform a careful astrometric analysis to avoid any offsets between the LW direct imaging in F444W and the LW channel Grism spectra in F444W.This allows us to extract the 2D Grism spectra associated with the astrometric position of any source, these spectra are then stacked and the 1D spectrum at the position of the source with an aperture of 3 pixels is extracted.We also perform a continuum subtraction from the 1D spectrum to remove any contaminant continuum and we focus on measuring [OIII] 4959,5007 and Hβ emission line fluxes.

Redshift Determination
We initially take care to ensure that, for galaxies that only have Lyα emission detected in their spectra, the emission line detected is indeed Lyα.There are initially five LAEs at z > 7 in the literature, imaged by NIRCam, for which Lyα is the only observed emission line (z7-GSD-3811, z8-19326, z7-13433, z7-20237 and GSDY).In order to confirm that this emission line is Lyα, we use the SED-fitting code Bagpipes 32 on the available NIRCam data Extended Figure 2. Best SED-fits for all systems with strong Lyα.The yellow dots are the observed photometry, the grey line shows the best fit.The physical parameters of the best fit are also indicated (see text for details).
Extended Figure 3. Same as Extended Figure 2.
to produce an improved photometric redshift of each LAE.One of the key diagnostics in this fitting becomes the bright [OIII] 4959,5007 emission lines, expected from high redshift objects, that fall into either the medium-band filter F410M (for z ∼ 6.8 − 7.6)/F430M (for z ∼ 7.3 − 7.8) or a lower flux consistent with continuum emission in this filter but a boost in the F444W filter (for z ∼ 6.8 − 9).This photometric-redshift fitting results in a redshift that is consistent with the redshift of the emission lines observed being Lyα in all but two of these galaxies (z8-19326 and z7-20237).The photometric redshift, constrained by a significant boost in F410M (z ∼ 7.1) for z8-19326, becomes inconsistent with the emission line detected being Lyα and therefore this galaxy is not included in our sample.While boosts in both F356W and F444W strongly indicate that z7-20237 is in fact at z ∼ 6.The original HST photometry reported in 31 also slightly favoured a z ∼ 6 solution making the observed emission line likely to be C IV.The removal of these two galaxies as contaminants underlines the challenges in determining galaxy redshifts with just one emission line and wide-band HST filters pre-JWST.The JWST offers significant advancement in this area given NIRSpec's sensitivity allowing the detection of multiple lines, but also given the frequency of programs that exploit the NIRCam medium-band filters around ∼ 4µm.All of the galaxies in our final sample show a clear [OIII] 4959,5007 excess in NIRCam imaging.We can therefore be confident that the emission line detected is indeed Lyα.
Following this previous step, in order to confirm that all components of each system are at a similar redshift, we again employ Bagpipes to estimate their photometric redshifts.We then compare these with the spectroscopic redshift of the main component.This fitting confirms that all the companion galaxies preferentially have high photometric redshifts, all of which are consistent with the spectroscopic redshift of the massive galaxy.In combination with the prior of the companion being closely located to a confirmed high-redshift galaxy, these become strong-photometric candidate companion galaxies.We also exploit all available ground and space based spectroscopic observations that cover the companions and use these to spectroscopically confirm the companion galaxy where possible.Below we discuss the available redshift constraints on our sample.
EGSY8p68, CEERS-1027 and JADES-GS+53.15682-27.76716constitute the sub-sample of LAEs that have photometric candidate companions with no spectroscopic confirmation.This lack of confirmation is due to a combination of factors.The observations of CEERS-1027 and JADES-GS+53.15682-27.76716were made with the NIRSpec MSA and the positioning of the shutter misses the companion, and as such there is no ability to place any spectroscopic constraints on these companions.EGSY8p68 has been observed by both the NIRSpec MSA and the ground-based telescope MOSFIRE/Keck.The MOSFIRE observations of the system cover the companion galaxies, however due to the seeing of the observations (∼ 0.7 ′′ ) and the small separation between the components of the system (∼ 0.1 ′′ ) resolving the components was not possible.Follow-up observations with the NIRSpec MSA reveal EGSY8p68 is likely an AGN given the presence of high-ionisation emission lines and a significant broadening of the Hβ emission line 17 .While one of the companion galaxies lies in the shutter of this observation, the available spatial resolution and strength of the emission lines from EGSY8p68 makes it challenging to disentangle the emission from each component.However, the reported NIR-Spec MSA flux is seven times fainter than that observed in the MOSFIRE spectrum 38 .Given the MSA shutter size is significantly smaller than the MOS-FIRE slit, one solution to this discrepancy in flux is that some of the Lyα emission is extended outside of the MSA slit.We note that the component EGSY8p68-C sits outside of the MSA slit position and as such this component may additionally host Lyα emission at the same redshift as the main component.Moreover, 17 claim that the system is consistent with being involved in a major merger and recent NIRSpec IFU observations of EGSY8p68 spectroscopically confirm they are physical companions [Carniani et al. in prep.].However, given these results are not yet published and the data is not yet public we continue our analysis of the redshift constraints assuming we do not have this spectroscopic confirmation.
Despite the lack of spectroscopic constraints on this sub-sample, the similarities in the SEDs of the main and companion galaxies are distinct.The galaxies in each system show similar colours, Lyman-breaks that are consistent and an emission-line boost, driven by [OIII] 4959,5007 and Hβ emission, in the same filters, as can be seen in Extended Figure 2 and 3.This results in a photometric redshift for each of the companions that clearly favours a z > 7 solution, and given the emission line feature in the F444W filter, the photometric redshift is further constrained to z = 7.6−9.While the redshift of these components are very confidently constrained to within this ∆z = 1.4,further constraining their redshift is challenging without spectroscopy.However, we make use of the UV luminosity function from 62 in order to estimate the expected number of galaxies within the volume that the [OIII] 4959,5007 emission redshift constraint constrains the companion galaxies to.We note that the photometric redshift from BAGPIPES (reported in Table 1) is much more tightly constrained than that provided by the [OIII] 4959,5007 boost alone, however we use this wider redshift constraint (∆z = 1.4) to provide an upper limit on the expected number of galaxies.Given that we have searched for companions in a 3 ′′ x3 ′′ region, we use this as the area for this estimate and the [OIII] 4959,5007 redshift constraints to produce a volume.We then integrate the z = 7 or z = 8 (dependent on the spectroscopic redshift of the main component) UV luminosity function down to the 5σ magnitude limit of the NIRCam images.We estimate that the expected number of galaxies in this field-of-view that are consistent with Hβ and [OIII] 4959,5007 producing the excess F444W flux and with observed UV magnitude is 0.09 +0. 160.04for JADES-GS+53.15682-27.76716,0.007 +0.009 −0.002 for CEERS-1027 and 0.003 +0.006 −0.001 for EGSY8p68.It is therefore very unlikely that another galaxy that is not associated with the spectroscopically confirmed LAE would be present in this volume.
JADES-GS-z7-LA is an extremely high EW LAE (400 Å) discovered in JADES NIRSpec observations of galaxies in the GOODS-S field 35 .The presence of a photometric candidate companion galaxy to this LAE has already been reported 35 , but the companion lies outside of the NIRSpec MSA shutter, and hence is not spectroscopically confirmed.We perform the same analysis as above, noting that the [OIII] 4959,5007 emission is in the F410M filter, limiting the redshift to between z = 6.8 − 7.6, and find the expected number of galaxies to be 0.13 +0.14 −0.05 .Therefore, once again, it is unlikely that the observed companion is at a redshift that makes it unrelated to the spectroscopically confirmed LAE.Moreover, we note the observed excess in the F115W filter that may be driven by Lyα emission in both the main and companion galaxy 35 , supporting the conclusion that these two galaxies are at the same redshift.
COSY and z7-13433 have companions that show excess fluxes in the F410M filter driven by [OIII] 4959,5007 emission.SED-fitting of their companions strongly prefers a z > 7 solution, and the presence of [OIII] 4959,5007 emission confines the redshift to z = 6.8 − 7.6.We again perform the same analysis as above, estimating the expected number of galaxies in the FoV.The expected number of galaxies is 0.03 +0.03 −0.01 for COSY and 0.04 +0.03 −0.01 for z7-13433 and thus we conclude it is very likely that the central and companion galaxies are associated with each other.
Moreover, the companions of both COSY and z7-13433 are within the slits used for the observations with X-shooter and MOSFIRE, respectively.The two components of each system are shown in Extended Figure 3, where we can identify both the positive (white) and negative (black) traces of the objects, caused by the ABBA nodding pattern of the observations.We find that, for z7-13433, two emission lines are offset by ∼ 3 pixels, which is consistent with the ∼ 0.5 ′′ offset between the two components in the slit.The Lyα emission that is coincident with the position of the companion is detected at 2.8σ and corresponds to z = 7.479, while the Lyα flux associated with the main galaxy is at 4σ corresponding to z = 7.482.Therefore, these two components are offset by ∼ 125 km/s.In the case of COSY, we expect a similar spatial offset between the two components, but the ∼ 0.7 ′′ seeing during the observations makes distinguishing between the two components very challenging.However, there is a second, previously unidentified, component to the Lyα emission of COSY at ∼ 0.9898µm.We attribute this 2.7σ emission line to COSY-B and note a velocity offset of ∼ 280 km/s between COSY-A and COSY-B.While these Lyα detections are clearly tentative, given they are coincident with the expected position in the 2D spectra of strong-photometric candidate galaxies and given their proximity to confirmed LAEs this boosts their likelihood of being a true detection.
Finally, we note that recent NIRSpec IFU observations of COSY spectroscopically confirm its companion galaxy 40 .However, as with EGSY8p68, given these results are not yet published and the data is not yet public we continue our analysis without this spectroscopic information.
Extended Figure 4.The X-shooter (above) and MOSFIRE (below) Gaussian-smoothed 2D spectra of COSY and z7-13433, respectively.These cover the wavelength range of the previously spectroscopically confirmed Lyα emission of the system.White circular apertures indicate the proposed two components of the Lyα emission that are coincident with the expected positions of the main and companion galaxies in the slit.The positive trace of the emission is indicated by white pixels, the negative trace by black pixels, this effect of positive and negative traces is caused due to the ABBA nodding procedure employed by both instruments.z7-GSD-3811 and GSDY are both in the FoV of the FRESCO NIRCam WFSS survey.Their final 2D WFSS spectra, shown in Extended Figure 5 show [OIII] 4959,5007 emission lines for both z7-GSD-3811, GSDY and their companions.For z7-GSD-3811 both components of the [OIII] 4959,5007 emission are present and there is a ∼ 3σ, faint detection of Hβ.Moreover, [OIII] 4959,5007 emission is also detected in the companion galaxy.This emission is offset from the central galaxy by 170 ± 70 km/s.The 2D spectrum of GSDY shows a clear detection of [OIII] 5007 in both the central and companion galaxies, and there is a very faint indication of [OIII] 4959 and Hβ in the companion galaxy, but at a less than 3σ confidence for both.The emission detected in the central and companion galaxy are offset by just ∼ 26 km/s, but given the low spectral resolution of the NIRCam WFSS, this is significantly smaller than the ∼ 70 km/s uncertainty on this measurement.
We discuss in the main text the use of Hβ to estimate the intrinsic Lyα flux and hence the escape fraction of Lyα photons, however, it can also be used as a probe of the instantaneous SFR.We assume a Hα to Hβ ratio of 2.85 and Case B recombination at T e = 10, 000 K and N e = 10 4 cm −3 and following 63 we estimate the instantaneous SFR which is reported in Extended Table 1.Given the lack of constraints on the redenning of these galaxies, we do not dust-correct the Hβ flux and as such these measured instantaneous SFRs are lower bounds.
The observation of multiple, spatially resolved components to each emission line is taken as spectroscopic confirmation of the companion galaxies.Moreover, the velocity offset observed in all of these galaxies is indicative of either a line-of-sight separation, or more likely, a difference in the line-of-sight velocities of the two objects as they interact with each other.We therefore consider this as further evidence that the systems are indeed interacting, as opposed to being two components of the same stable system.Extended Table 1.Properties inferred from NIRCam Grism spectroscopy.The observed redshift and fluxes of emission lines (in units of 10 −18 erg/s/cm 2 ) detected in the NIRCam/WFSS spectra.The instantaneous starformation rate (SFR) is not dust-corrected given the lack of constraints on dust in these galaxies.
GN-z11 was recently observed with the NIRSpec IFU which has revealed the presence of three regions emitting either HeII (GN-z11-C) 30 or Lyα emission (GN-z11-B and a more tentative candidate) 39 .While none of these regions show evidence of a UV continuum in deep NIRCam imaging, a tentative indication of CIV]λλ1548,1551 emission in GN-z11-B and analysis of the emission in GN-z11-C 30 indicate they are indeed being driven by a UV-faint stellar population in these regions.As such, we conclude that GN-z11 does have spectroscopically-confirmed companion galaxies.Moreover, the presence of multiple LAEs in the local vicinity of GN-z11 suggests that peculiar velocity is not playing a vital role in the escape of Lyα through the neutral IGM (see 39 for a detailed analysis of the Lyα emission present in this system).Moreover, the density of objects in this field makes GN-z11 a likely protocluster core 39 .

SED fitting
We then use Bagpipes, with the redshift fixed to that of the spectroscopic redshift of the system, using different SFHs (constant, delayed, burst and burst+constant) to fit the SED of all central and companion galaxies, as seen in Extended Figure 2 and 2. The best fit model is then obtained by taking the SFH that leads to the smallest Bayesian Information Criterion (BIC) as described in 64 .We then report the galaxy properties returned from this method in Table 1.

Close companion fraction of non-LAEs
In order to ascertain the companion fraction of a mass-matched sample of z > 7 galaxies, we take 30 galaxies with stellar masses from 7.4 < log(M * [M ⊙ ]) < 9.3 from 65,23,66,67 .We find that 47% of these galaxies have photometric-candidate companions within 5 ′′ of the central galaxy.However, many of these galaxies have not been spectroscopically followed-up and as such, we do not have constraints on their Lyα emission.
Unfortunately, examples of spectroscopically confirmed z > 7 non-Lyα emitting galaxies that have been observed by high-resolution spectrographs are rare in the literature.As such it is challenging to probe the companion fraction for known non-LAE galaxies.Moreover, the challenges of slit/aperture spectroscopy mean that it is incredibly rare to have spectroscopic information on the local environment of such galaxies and hence, while a galaxy may be confirmed as non-LAE, there is no information regarding the Lyα emission of any companion that may or may not be present.
The recent JADES data release of NIRSpec and NIRCam observations in GOODS-S includes seven z > 7 galaxies, observed using the NIRSpec highresolution G140M grating, that show no Lyα emission 13 .This does not provide any Lyα information on any would-be companion galaxies, however, the companion fraction for this sample, ∼ 43%, is once again consistent with both our larger sample, and that of Puskás et al. (in prep.).We therefore conclude that the 100% companion fraction of our LAE sample is clearly atypical compared to the companion fractions of samples not selected for Lyα emission.

Close companions are necessary but not sufficient for Lyα emission
We note the recently reported z = 9.3 merging system 68 with an absence of detected Lyα emission.The reported specific-star-formation rate of this galaxy (sSFR ∼ −8) is consistent with the high sSFR seen in our sample of merging LAEs, and therefore, one may naively posit that this is a contradiction to our results.We argue that non-detection of Lyα emission from a merging system is entirely reasonable and in fact our simulation modelling predicts this.While we see significant boosts in the SFR and hence intrinsic Lyα emission from merging systems, the escape of Lyα emission is only possible when the escape fraction of Lyα emission out of both the host galaxy and the surrounding IGM is non-zero.We explore the effect of viewing angle on the observed Lyα emission escaping our simulated merging systems and find that while the majority of viewing angles facilitate the escape of Lyα, several viewing angles exist for which the only Lyα emission escaping the galaxy is diffuse and is therefore unobservable.Moreover, without the presence of a large ionised bubble, the neutral IGM will not facilitate the escape and observation of Lyα emission.Ultimately we conclude that close companions are necessary but not sufficient for the observation of Lyα from galaxies deep in the epoch of reionisation.A combination of a large-scale ionised bubble and a preferable line-of-sight, stripped of neutral hydrogen by the interaction with the companion, are also required.

Azahar simulations
The new Azahar simulations are a suite of multiple models spanning various combinations of canonical hydrodynamics, magnetic fields, radiative transfer, and cosmic rays physics with the aim of understanding their effects on galaxy formation as well as their complex interplay.The simulations are generated using the magneto-hydrodynamical code ramses 69 , which employs a constrained transport method for divergence-less evolution of the magnetic field 70 .The main target of study for Azahar is a massive spiral galaxy with M halo (z = 1) ∼ 2 • 10 12 M ⊙ in a relatively large zoom-in region, approximately 8 cMpc across along its largest axis.Azahar has a maximum spatial resolution in cells with a full cell-width of ∆x ∼ 20 pc (or equivalently, an approximate cell half-size or radius of 10 pc), and refinement is triggered whenever its size is larger than a quarter of its Jeans length, or it contains a total mass larger than 8 m DM , where m DM ∼ 4.5 • 10 5 M ⊙ is the mass of dark matter particles in the zoom region.The stellar mass particle is m * ∼ 4 • 10 4 M ⊙ .Azahar follows the set of physics presented by the pathfinder Pandora 46 .We provide here a brief summary, and refer the reader to that reference and Martin-Alvarez et al. (in prep.) for further details.
In this work, we investigate one of the most complete simulations in the Azahar suite, the RTnsCRiMHD model.In this manuscript, we refer to the RTnsCRiMHD simulation simply as Azahar.Radiative cooling is modelled both above and below 10 4 K 71,72 .The model adopts a magneto-thermoturbulent prescription for star formation 73,74 , mechanical supernova (SN) feedback 50 and astrophysically-seeded magnetic fields through magnetised SN feedback 75 .This feedback injects 1% of the SN energy as magnetic energy, and intends to reproduce the approximate magnetisation of SN remnants.This particular simulation of Azahar also includes cosmic rays modelled as an energy density allowed to anisotropically diffuse, evolved with an implicit solver 76,77 .In Azahar, cosmic rays are exclusively sourced by SN feedback, with each event injecting 10% of their total energy as cosmic rays.We assume a constant diffusion coefficient for cosmic rays of κ CR = 3 • 10 28 cm 2 s −1 , and do not account for cosmic ray streaming in this model.This model also includes on-the-fly radiative transfer 47 , with a configuration similar to that of the SPHINX simulations 48 , featuring three energy bins for radiation spanning the ionisation energy intervals from 13.6 eV to 24.59 eV (HI ionisation), 24.59 eV to 54.42 eV (HeI ionisation), and above 54.42 eV (HeII ionisation).Finally, recent work has shown that when expanding simulation models to account for cosmic ray feedback physics at comparable resolutions, the escape fractions of LyC photons from galaxies is lower than in their absence 78 .
Therefore, Azahar features a realistic (yet computationally taxing) ISM model that considerably approaches the resolution required to converge in the propagation and escape of ionising photons from galaxies 50 .We note that, while close, our resolution still falls short from that regime and that even more sophisticated ISM models 79 may be important for a complete understanding of photon propagation in the ISM.Opportunely, a considerable part of our results relies on gas being ejected from galaxies during merger events, which implies that our estimate of LAE detectability during these events will be more resilient to these caveats.

Simulated ionised bubble evolution in the neutral IGM
In the high resolution patch of the Universe probed by our simulation (≲ 8 Mpc on a side), centred on the main progenitor of our z ∼ 1 galaxy, the first ionised bubbles emerge at z ∼ 20.By z ∼ 15 several ionised bubbles exist that are a few (physical) kpc in radius, and rapidly begin merging as redshift approaches z ∼ 10.This drastically increases the mean free path of ionising photons, with the coalesced bubble reaching scales of 100 kpc.The resulting main bubble continues to grow, reaching a radius ∼ 0.5 Mpc at z ∼ 7.3.
For intrinsically bright LAEs, as studied here, a ≲ 0.5 Mpc local ionised bubble may be sufficient for these sources to be detectable 80 .We further note that due to ongoing mergers our simulated galaxies have considerable relative peculiar velocities of up to 230 km s −1 .Sufficiently large velocities relative to the intervening neutral IGM gas may reduce damping wing suppression for particular lines-of-sight, hence requiring smaller local ionised bubbles to facilitate the detection of Lyα.Our simulation also does not consider star formation outside of the high-resolution zoom-in region.Hence, our simulated ionised bubbles could be even larger in size, e.g.due to the presence of a neighbouring, large-scale overdensity or large-scale clustering 81 on spatial scales beyond our zoom-in region.
Nevertheless, knowing from our observations that the LAEs sit in sufficiently large ionised bubbles, here we primarily focus on the probability of Lyα photons escaping the host galaxy and its intermediate surroundings.We find that a large amount of gas remains neutral within the local ionised bubble, actively feeding the galaxies along the cosmic web.

Galaxy selection and measurements in the simulations
In order to select galaxies in our simulation, we employ the halomaker software 82 to detect and characterise dark matter halos.We identify the three progenitor galaxies of interest (as well as the galaxy merging with the main progenitor at z ∼ 8.1) and follow them through time by tracking their innermost stellar particles.Their centres are determined using a shrinking spheres algorithm applied to their stellar component 83 , and they are assigned their corresponding halo as obtained by halomaker.In order to select the system for study, we broadly reviewed the evolution of all galaxies with stellar masses M * > 10 7 M ⊙ in the redshift interval z ∈ [9.0, 6.0].The most promising candidate, the main progenitor of the main Azahar galaxy was selected for further investigation.
In order to assign measurements to individual galaxies, we measure values within their galactic region, defined by the radius r gal < 0.2 r rvir .During the merging stages, we employ r = min[r gal , 0.45D ij ], where D ij is the distance between the i and j progenitors, down to a minimum distance of 1.5 kpc.
We have briefly reviewed the observability of lower mass pairs within the simulated domain.While these present a similar behaviour during mergers to that of our main studied system, their observability and that of other galaxies with similar masses remains low.This is due to their low stellar masses generating low luminosities.Their low stellar masses also make them more vulnerable to disruption during powerful star formation bursts, which increases their escape fractions 49 .Despite this, most isolated systems present quiet star formation histories and in the absence of mergers, have low observabilities throughout their evolution.

RASCAS post-processing
We post-process the simulation with the publicly available, massively parallel code RASCAS 52 , for modelling the Lyα emission and resonant scattering in our simulations.RASCAS accounts for two sources for the Lyα emission, recombination and collisional excitation.For recombination, RASCAS adopts the case B recombination coefficient from 84 , and the fraction of recombinations producing Lyα photons from 85 .For collisional excitation, RASCAS uses the fitting function for the collisional excitation rate from level 1s to 2p from 86 .We cast N MC = 10 6 Monte Carlo photon packets to sample the real Lyα photon distributions, from recombination and collisional excitation, respectively.We sample N MC = 10 7 photon packets for the images along the LOS in Extended Figure 4.
After each scattering event, the Lyα photon changes its frequency and direction according to a phase function as implemented in RASCAS.RASCAS adopts the phase function in 87,88 for the scattering of Lyα photons around line centre and Rayleigh scattering for Lyα photons in the line wing.At high H i column density, Lyα photons will scatter many times locally until they shift in frequency large enough so that they can have a long mean free path.To reduce the associated computation cost, RASCAS adopts a core-skipping mechanism 89 to transit the photons to line wings while avoiding local scattering in space.RASCAS also implements the recoil effect and the transition due to deuterium with an abundance of D/H = 3 × 10 −5 .The dust distribution is modelled according to the dust model described in 53 .Dust can either scatter or absorb Lyα photons.The probability of scattering is given by the dust albedo a dust = 0.32 following 90 .The dust scattering with the Henyey-Greenstein phase function 91 and the asymmetry parameter is set to g = 0.73 following 90 .We generate synthetic images and spectra with a peeling algorithm described in 92,93,94 , along 12 lines-of-sight uniformly sampled using the healpix algorithm 95 .
Finally, we note that the dust distribution has a large effect on the Lyα radiative transfer as well as the [OIII] 4959,5007 and Hβ emission that is used to infer the intrinsic Lyα emission from observations.In Extended Figure 6, we show how reducing the amount of dust changes the observed Lyα images.We see the dust is able to completely mask the Lyα emission of some individual galaxies as resonant scattering of H i increases the probability of Lyα photons encountering dust grains.The dust can suppress the Hα emission by up to 20% around the center of the three merging galaxies, if we assume intrinsic Lyα can be converted directly to intrinsic Hα emission.We therefore conclude that inclusion of dust modelling, as done here, is fundamental for understanding observed Lyα emission.However, our results regarding the boosted intrinsic Lyα emission in star formation bursts and enhanced Lyα escape fraction in mergers, due to the channels cleared of local neutral hydrogen are robust to the assumed dust model.
We further consider the Lyα profiles escaping along the line-of-sight shown on the upper panel of Extended Figure 4 at z = 7.3 after applying a reasonable IGM damping wing from 96 .We note that the shape of the Lyα profile is very complex and sensitive to both the LOS 97 and IGM attenuation 34 .As such, further analysis of this is needed to constrain the expected Lyα profiles from simulations requiring a larger simulated galaxy sample with careful IGM attenuation modelling as well as a larger observational sample to compare to, which is beyond the scope of this work.However, as an immediate consistency test, we take the Lyα profile from our simulated system at z = 7.3 after IGM attenuation 96 and compare the velocity offset and FWHM to that reported in our observational comparison, EGSY8p68.We find that the observed Lyα profile 17 is well matched with the simulated profile's offset of ∼ 200 km/s and FWHM of ∼ 200 km/s, and shows a clear asymmetric line shape as observed.

Figure 1 .
Figure 1.An abundance of galaxy mergers seen with NIRCam.Cutouts of JWST images (as different surveys make use of different filters, the filter within which the companion is most clearly present is shown: F182M for GSDY and z7-GSD-3811, F200W for EGSY8p68, F115W for CEERS-1027, JADES-GS-z7-LA and COSY and F150W for JADES-GS+53.15682-27.7671and z7-13433) showing the multiple components of each system.All images are smoothed with a Gaussian of FWHM ∼ 0.7 kpc to enhance the visibility of the companion, except for EGSY8p68 where the three components become unresolved if we smooth the images.We include the unsmoothed images as sub-panels in the top right of each panel in order to emphasise that the companion is always visible before smoothing.We use different colorbar ranges for each system in order to make the companion galaxy well defined, hence some unsmoothed images can have what appears like a nosier background than others.The position of the main target (*A in the text) is displayed by an orange arrow, the position of the first companion (*B) is indicated by a blue arrow and the third companion, if any, by a green arrow.The name of each candidate is indicated in the bottom right corner, and in the bottom left of each panel is a black scale showing 0.5 ′′ .
Figure 1.An abundance of galaxy mergers seen with NIRCam.Cutouts of JWST images (as different surveys make use of different filters, the filter within which the companion is most clearly present is shown: F182M for GSDY and z7-GSD-3811, F200W for EGSY8p68, F115W for CEERS-1027, JADES-GS-z7-LA and COSY and F150W for JADES-GS+53.15682-27.7671and z7-13433) showing the multiple components of each system.All images are smoothed with a Gaussian of FWHM ∼ 0.7 kpc to enhance the visibility of the companion, except for EGSY8p68 where the three components become unresolved if we smooth the images.We include the unsmoothed images as sub-panels in the top right of each panel in order to emphasise that the companion is always visible before smoothing.We use different colorbar ranges for each system in order to make the companion galaxy well defined, hence some unsmoothed images can have what appears like a nosier background than others.The position of the main target (*A in the text) is displayed by an orange arrow, the position of the first companion (*B) is indicated by a blue arrow and the third companion, if any, by a green arrow.The name of each candidate is indicated in the bottom right corner, and in the bottom left of each panel is a black scale showing 0.5 ′′ .

Figure 2 .
Figure 2. Comparison of Azahar to an observed galaxy merger.(Left panels) NIRCam F150W (above) and HST F160W (below) imaging of the LAE EGSY8p68.The NIRCam imaging reveals three components to the system that were previously unresolved by HST.(Right panels) Analogue galaxy merger from the Azahar simulation: a) large-scale view of the filaments encompassing the galaxy merger (blue: HI density, grey-black: HII density, yellow-white: F150W intensity); b) simulated NIRCam F150W observation of the Azahar merger; c) fully resolved simulated observation in the NIRCam filters; and d) hydrogen gas and ionising radiation properties (blue: HI density, grey-black: HII density, yellow-white: F150W intensity, purple: LyC radiation energy density, orange: HeII ionising radiation energy density).The appearance of the observed system is very well reproduced by an ongoing merger of three galaxies, with an extended escaping LyC radiation field.

Figure 3 .
Figure 3.The evolution of galaxy properties with an ongoing merger.(Rows 1-2) The time evolution of our simulated merger system at z = 8.1 (close up and large scale view), 7.3 and 6.1.First row: maps of HI density (blue), HII density (grey-black), NIRCam F150W intensity (yellow-white), LyC radiation energy density (purple) and HeII ionising radiation energy density (orange).Second row: the RASCAS intrinsic Lyα flux (blue) and scattered Lyα flux (orange).(Rows 3-5)The SFR for each of the three galaxies and the final system (black line).Dot symbols show mergers with companions (filled markers) and with other secondary systems (open markers); the intrinsic Lyα luminosity of the three galaxies and the entire system; the fraction of escaped Lyα luminosity filtered for pixels with fluxes higher than a given limit (median (dark green curve), quartiles (darker band), and min-max of the distribution across 12 lines-of-sight (lighter band)).Distance between the main progenitor and its companions is shown as well.The grey vertical lines in Rows 3-5 correspond to the redshifts shown in Rows 1-2.The rightmost panels show key physical relations for the observations (black) and the simulations (pink).

Figure 5 .
NIRCam WFSS Grism spectra of z7-GSD-3811 (above) and GSDY (below).Within each panel, the continuum subtracted 2D spectrum (above) and 1D spectrum (below) of the central galaxy are shown.The upper left sub-panel shows the F444W image of the system where two components can be identified.The horizontal dashed lines indicate the position of the central galaxy in both the image and 2D spectrum.The three dashed lines in the 1D spectrum denote the expected positions of the [OIII] 4959,5007 doublet (central and right dashed lines) and Hβ (left dashed line) emission lines given the object's spectroscopic redshift, while the grey shaded region denotes the one-sigma noise level.

Extended Figure 6 .
Dust opacities superimposed on synthetic Lyα images at z = 7.3.The top row shows the intrinsic Lyα emission, the middle row shows the scattered Lyα emission generated with a default dust model, and the bottom row shows the scattered Lyα emission with practically no dust (dust opacity reduced by a factor of 10 6 ).The three different columns show the galaxy from three different directions.Red and orange iso-contours correspond to a dust optical depth of 0.2 at the Lyα wavelength and Hα wavelength, correspondingly.Note that the dust does not act as a screen for the Lyα emission because of the resonant nature of Lyα scattering.
39,30ompanions of GN-z11 are spectroscopically confirmed39,30however given that they are very UVfaint we cannot estimate their properties.† The companions of both COSY and EGSY8p68 have recently been spectroscopically confirmed by NIRSpec IFU observations 40 ;[Carniani et al. in prep.].
30Table1: The properties of Lyα-emitting galaxies.Physical properties of all the systems studied in this paper.The columns are : (1) spectroscopic (italic) or photometric redshift of the candidate, (2) stellar mass, (3) Star formation rate averaged over at most the last 100 Myr, (4) separation between the main component (A) and each companion, (5) luminosity of observed Lyα emission (6) escape fraction of Lyα photons and (7) reference for the original spectroscopic confirmation and the Lyα escape fraction measurement.* ‡ JADES-GS+53.15682-27.76716has been abbreviated to JADES-GS.