Main

Young, vigorously star-forming galaxies have been identified in the very early Universe1,2,3. These galaxies should be excellent sources of Lyman-α emission (Lyα; wavelength (λ) = 1,215.67 Å)—the intrinsically brightest emission line4—which stems from the recombination of hydrogen that has been ionized by their young stellar populations. However, deep in the epoch of reionization, 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 redshift6,7, and this neutral gas is expected to resonantly scatter Lyα emission. Hence, as a result of ‘local’ attenuation by a gas-rich interstellar medium (ISM) and scattering by a neutral IGM, Lyα emission should be detectable only towards the end of the reionization era, about one billion years after the Big Bang4,8,9,10.

Although the decreasing observability of Lyα emission with increasing redshift has been repeatedly observed11,12,13, this picture has been challenged by the occasional, surprising detection of Lyα emission in several galaxies deep in the reionization era14,15. It has been suggested that Lyα can escape through the neutral IGM if the galaxies reside in sufficiently large ionized bubbles embedded in 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 objects19,20,21,22,23,24,25. However, a solution to the escape of Lyα emission through the ISM and circumgalactic medium of what are expected to be very gas-rich galaxies remains elusive. Before the advent of the James Webb Space Telescope (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 redshift (z) > 7 with high-resolution spectrographs and that 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 in the GOODS-North, GOODS-South, EGS and COSMOS fields and have been observed as part of five different programmes: Public Release IMaging for Extragalactic Research (PRIMER, primary investigator (PI) Dunlop), the First Reionization Epoch Spectroscopic Complete Survey (FRESCO, PI Oesch)27, the Cosmic Evolution Early Release Science survey (CEERS, PI Finkelstein)28, the JWST Advanced Deep Extragalactic Survey (JADES, PI Eisenstein)29 and Director’s Discretionary Time (DDT) programme 4426 (providing NIRSpec integral field unit (IFU) observations of GN-z11, PI Maiolino)30. Six of these nine Lyα-emitting galaxies are known to lie within overdensities21,23,24,25,31, three of which also probably host accreting black holes16,17,18. Moreover, it has been shown that the remaining galaxies are incapable of alone blowing a large enough ionized bubble to facilitate the escape of Lyα through the neutral IGM, indicating that an underlying population of faint galaxies must surround these LAEs25.

Table 1 The properties of Lyα-emitting galaxies

However, the presence of spectroscopically confirmed galaxies within known ionized bubbles that do not show Lyα emission13,23,25 indicates that there must be further local processes at play driving Lyα emission deep into the epoch of reionization. With this in mind, we use the unparalleled sensitivity and resolution of JWST/NIRCam to accurately determine the properties of these nine LAEs with equivalent width greater than 10 Å, at z > 7.

Remarkably, NIRCam images of these galaxies show the presence of several components of all these LAEs, as shown in Fig. 1. Given the tight redshift constraints provided by strong [O iii]4959,5007 emission falling into medium-band filters, we estimate the probability of these being unrelated, high-redshift objects and find this is always less than 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 Data Figs. 1 and 2. We then fit each SED by assuming several parametric star formation histories (SFHs: burst, constant, delayed, exponential), allowing a large range for all parameters32. In all cases, we deduce that the photometric redshifts of the companions are consistent with that of their spectroscopically confirmed main component.

Fig. 1: An abundance of galaxy mergers seen with NIRCam.
figure 1

Cutouts of JWST images (as different surveys use different filters, the filter in 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.7671 and z7-13433) showing the 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 subpanels in the top right of each panel to emphasize that the companion is always visible before smoothing. We use different color-bar ranges for each system to make the companion galaxy well defined, so some unsmoothed images have what seems like a nosier background than others. Orange arrow, position of the main target (*A in the text); blue arrow, position of the first companion (*B); green arrow, position of the third companion, if any. The name of each candidate is indicated in the bottom-right corner. Scale bars, 0.5″.

To spectroscopically confirm their identity as companions, we exploit existing X-shooter/Very Large Telescope33, MOSFIRE/Keck34,35, NIRSpec on-board JWST36 and NIRCam wide-field slitless spectroscopy (WFSS) observations37 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/NIRSpec observations. For the remaining six systems in our sample, we are able to spectroscopically confirm the companion galaxy as it falls within the field of view (FoV) of the observations of the main LAE. For two companions in our sample, we see tentative (3σ) evidence of Lyα emission in Extended Data Fig. 3 (one of which has recently been confirmed by NIRSpec IFU observations38). For a further two companions, we observe the [O iii]5007 emission line at >8σ (as seen in Extended Data Fig. 4, with fluxes reported in Table 2), and NIRSpec IFU observations spectroscopically confirm the companions of the final two systems (refs. 30,39 and Carniani et al., manuscript in preparation). Therefore, for all systems for which we have resolved spectroscopic observations of the companion (67%), we spectroscopically confirm its nature as a physical companion (see Methods for further details of the spectroscopic analysis).

Table 2 Properties inferred from NIRCam Grism spectroscopy

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 detected13. 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., manuscript in preparation). 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 continuing interactions drive an increase in the detectability of Lyα emission during the epoch of reionization, we explore comparable galaxy mergers using the Azahar simulation suite (Martin-Alvarez et al., manuscript in preparation). This simulation suite was performed before we obtained the NIRCam data, but remarkably, we find many simulated galaxies that match the photometry and spatial geometry of our objects. Azahar is a cosmological, high-resolution, zoom-in simulation that uses a magnetohydrodynamical solver together with state-of-the-art cosmic ray feedback (see the pathfinder Pandora project40 and Methods). Most importantly for this work, Azahar also features on-the-fly radiative transfer41,42,43 capable of self-consistently reproducing reionization while fully modelling the ISM of galaxies (with a maximum spatial resolution of 10 pc) and, crucially, resolving the propagation and escape of ionizing radiation on the ISM scales44,45. We post-process our simulations with the publicly available RAdiative SCattering in Astrophysical Simulations (RASCAS) code46 to account for both the production and resonant scattering of Lyα photons as well as scattering or absorption of Lyα photons by dust47.

As a result of a substantial overdensity of galaxies in our simulated volume, the main progenitor undergoes repeated mergers with other galaxies brought in by the cosmic filaments, often involving several companions or mergers in rapid succession as has been observed in ref. 48. To highlight one such merger occurring at z ≈ 7.3, Fig. 2b–e 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 Data Fig. 5 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 × 108M, M*,2 = 3 × 108M and M*,3 = 8 × 107M at this redshift, which are comparable to the stellar masses of the EGSY8p68 galaxies (see Table 1). Furthermore, the simulated star formation rates (SFRs) vary between 1 and 10 M yr−1 (shown in Fig. 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 although it varies, its average value is 1043 erg s−1 (shown in Fig. 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 EGSY8p6821. We thus conclude that the simulated system matches all key observed photometric and spectroscopic properties of EGSY8p68 (see also Methods).

Fig. 2: Comparison of Azahar to an observed galaxy merger.
figure 2

a, NIRCam F150W (top) and HST F160W (bottom) imaging of the LAE EGSY8p68. The NIRCam imaging shows three components of the system that were previously unresolved by HST. be, analogue galaxy merger from the Azahar simulation: large-scale view of the filaments encompassing the galaxy merger (b; blue, H i density; grey-black, H ii density; yellow-white, F150W intensity), simulated NIRCam F150W observation of the Azahar merger (c); fully resolved simulated observation in the NIRCam filters (d) and hydrogen gas and ionizing radiation properties (e; blue, H i density; grey-black, H ii density; yellow-white, F150W intensity; purple, LyC radiation energy density; orange, He ii ionizing radiation energy density). The appearance of the observed system is very well reproduced by a continuing merger of three galaxies with an extended escaping LyC radiation field.

Fig. 3: The evolution of galaxy properties with a continuing merger.
figure 3

a,b, The time evolution of our simulated merger system at z = 8.1 (close-up and large-scale views) (a), 7.3 and 6.1 (b). a, Maps of H i density (blue), H ii density (grey-black), NIRCam F150W intensity (yellow-white), LyC radiation energy density (purple) and He ii ionizing radiation energy density (orange). b, RASCAS intrinsic Lyα flux (blue) and scattered Lyα flux (orange). c, SFR for each of the three galaxies and the final system (black line). Dot symbols show mergers with companions (filled markers) and other secondary systems (open markers). d, Intrinsic Lyα luminosity of the three galaxies and the entire system. e, 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 LOSs (lighter band)). Distance between the main progenitor and its companions is shown as well. ce, Grey vertical lines correspond to the redshifts shown in a and b. Right, key physical relations for the observations (black) and the simulations (pink).

The observations shown in Fig. 2a clearly indicate the superior resolution and sensitivity of JWST (top) compared to the Hubble Space Telescope (HST, bottom). The tricomponent nature of EGSY8p68 is entirely unresolved in existing HST observations but clearly identifiable in the F150W NIRCam imaging. Fig. 2b shows the large-scale view (150 kpc across) of the three filaments encompassing the merging system, where large amounts of H i gas (blue) along the filaments are feeding the star formation in these systems, resulting in substantial NIRCam F150W emission (yellow-white). This active star formation is driving ionized hydrogen bubbles (grey) away from the filaments and into the low-density regions. Fig. 2c shows the simulated JWST-like 150W observation (with dust extinction modelled as an absorption screen along the line of sight (LOS)), where we select the LOS 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 probable physical separation of the observed galaxies at z ≈ 7. Interestingly, the mock F150W emission shows compact galactic cores analogous to these JWST observations, with diffuse emission (blending with the background) emerging from the extended, low-density stellar disks. Fig. 2d shows a synthetic colour-composite observation using the JWST filters at the full resolution of our simulation. Although the two companion galaxies are more diffuse, with stars actively forming in their disks, the main galaxy has a much more compact core due to the preceding rapid growth through repeated mergers. Fig. 2e is a colour composite of various simulated properties relevant for Lyα detectability. The emission from LyC ionizing photons (purple) and He ii ionizing photons (orange) is distributed on a scale of a few (physical) kiloparsecs. Importantly, there is a clear separation of the stellar emission and the H i gas during the merger event, with the ionizing radiation escaping from star-forming regions.

To provide a quantitative confirmation and in-depth physical understanding of this effect, we explore in Fig. 3 the time evolution of the three merging galaxies in detail. In Fig. 3a, the colours encode the same quantities as in Fig. 2e, and Fig. 3b shows 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 view 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 continuing merger of our galaxies 1, 2 and 3 at z ≈ 7.3 (with a different orientation with respect to Fig. 2), where bursty star formation drives ionized channels through the ISM, with a substantial amount of gas also tidally stripped from the galaxies. The fourth column shows a post-merger view of the newly formed disk-dominated galaxy at z ≈ 6.1. Here, the H i is confined to the merger remnant galaxy, and the radiation can only escape through the channels opened by collimated, cosmic-ray-driven biconical outflows. Focusing now on the Lyα radiation, all three galaxies have substantial intrinsic Lyα emission due to their continuing star formation. However, at this particular viewing angle, because of 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, Fig. 3c shows 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 Fig. 3d, 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−18erg 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 Fig. 3e, 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 LOSs. As the fraction of escaping Lyα emission has several high peaks, we conclude that the Lyα emission should be detectable for several of the redshifts studied here. Interestingly, during close galactic encounters (as shown by the distance trajectories of three main galaxies shown in Fig. 3e by coloured lines) and mergers, substantial gas tidal stripping and bursty star formation feedback lead to the opening of more low-column density sightlines, 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 the Lyα escape fraction. Our findings about higher intrinsic Lyα emission and enhanced Lyα escape fraction in mergers are robust to the assumed dust model (Methods and Extended Data Fig. 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 AGN49.

The rightmost panels of Fig. 3c–e show various physical relations for the observations (black points) and simulations (pink ‘violin’ symbols, denoting the distribution of Lyα escape fraction over 12 LOSs) using the same flux filter as above, demonstrating a very good correspondence between our simulations and our observed LAEs. Interestingly, for both observations and simulations we do not find clear correlations between the Lyα escape fraction and any other galaxy properties (distance, SFR, stellar mass), further corroborating our finding that ‘favourable’ LOSs with evacuated neutral hydrogen channels lead to directional Lyα photon leakage.

In this Article, 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 reionization: galactic mergers driving high intrinsic Lyα emission in the host galaxy, a ‘favourable’ LOS cleared of local neutral hydrogen in the host galaxy by tidal interactions with companions and by star formation feedback and a sufficiently large ionized bubble facilitating the escape of Lyα emission through the IGM.

Methods

Data

The emergence of JWST with its substantial technical advancements over both ground- and space-based telescopes has the potential to provide new hints as to the driving forces behind the observed Lyα emission discussed in this Article. We therefore decided to analyse all spectroscopically confirmed LAEs at z > 7 with JWST/NIRCam imaging. This search originally showed 11 such LAEs; however, further constraints on the SED of these galaxies from NIRCam imaging shows that 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 have companions that are identifiable in publicly available NIRCam imaging, and GN-z11 has spectroscopically confirmed companions that do not show an ultraviolet (UV) continuum30,39. These galaxies are found in the GOODS-North, GOODS-South, EGS and COSMOS fields, which have been observed with NIRCam by FRESCO (PI Oesch, ID 1895), JADES (PI Eisenstein, ID 1180), CEERS (PI Finkelstein, ID 1345) and PRIMER (PI Dunlop, ID 1837), respectively.

The PRIMER survey covers 150 arcmin2 of the COSMOS field with the F090W, F115W, F150W, F200W, F277W, F356W, F410M and F444W filters all at a 5σ limiting depth of AB magnitude (mAB) > 28. The CEERS survey took imaging of 100 arcmin2 of the EGS field in the F115W, F150W, F200W, F277W, F356W, F410M and F444W filters, reaching a 5σ limiting depth of mAB > 29 in the short-wavelength filters and mAB > 28.4 in the long-wavelength (LW) filters50. The June 2023 JADES data release surveys a 36 arcmin2 area of the GOODS-S field using the F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M and F444W filters, reaching 5σ limiting depths of mAB ≈ 30 (ref. 51).

The FRESCO survey’s primary aim is to obtain NIRCam 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 the F182M and F210M filters, and direct imaging obtained in the F444W filter also provides us with a total of three filters for our galaxies in GOODS-S/GOODS-N. FRESCO’s LW WFSS brings substantial advantages as the F444W filter covers the wavelength range of 3.9–5 μm, therefore observing both [O iii]4959,5007 and Hβ at z ≈ 7–9. These emission lines ordinarily contaminate imaging in the F444W filter, so constraints on the fluxes of these emission lines allow us to reduce the uncertainty in galaxy properties such as stellar mass, SFR and age. We use these constraints on the emission lines in F444W to create a mock F430M filter. Furthermore, when we fit the SED (discussed further below), we check the fluxes of the [O iii]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.

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. 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 v.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 and (Stage 3) producing the final mosaic. We take more 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 SEXTRACTOR52 on point-spread function (PSF)-matched images, extracting all objects with more than 4 px (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 use the photometry of the LAE and its companion reported in ref. 53, and for EGSY8p68, we model the shape of each component using GALFIT54 assuming a Sersic profile and NIRCam PSFs, simulated using WEBBPSF. We first run GALFIT on one of the NIRCam images in which 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 Data Fig. 5 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 Multi-Shutter Array (MSA) 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 use all existing ground-based observations of our sample of LAEs, four of which 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 calibration, sky subtraction and flat fielding as well as further standard data-reduction steps 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 ≈ 3,380, 0.1798per px and R ≈ 8,900, 0.158″per 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 LW Grism as part of the FRESCO observing programme. FRESCO’s LW WFSS covers the wavelength range of 3.9–5 μm using the F444W filter, therefore observing both the [O iii]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 data27,55. Moreover, the substantially superior spatial resolution of NIRCam LW WFSS (60 Mas allows us to resolve emission lines where MOSFIRE and X-shooter would fail to do so.

We reduce the Grism data following the steps described in ref. 56, including flat fielding, background subtraction, 1/f noise subtraction and World Coordinate System 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 one-dimensional spectrum at the position of the source with an aperture of 3 px is extracted. We also perform a continuum subtraction from the one-dimensional spectrum to remove any contaminant continuum, and we focus on measuring [O iii]4959,5007 and Hβ emission line fluxes.

Redshift determination

We initially take care to ensure that for galaxies that have only 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). To confirm that this emission line is Lyα, we use the SED-fitting code BAGPIPES32 on the available NIRCam data to produce an improved photometric redshift of each LAE. One of the key diagnostics in this fitting becomes the bright [O iii]4959,5007 emission lines expected from high-redshift objects, which fall into either the medium-band filter F410M (for z ≈ 6.8–7.6) or 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 to be Lyα in all but two of these galaxies (z8-19326 and z7-20237). The photometric redshift, constrained by a 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. Boosts in both F356W and F444W strongly indicate that z7-20237 is in fact at z ≈ 6. The original HST photometry reported in ref. 31 also slightly favoured a z ≈ 6 solution, making the observed emission line probably 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 an important advancement in this area, given NIRSpec’s sensitivity, allowing the detection of several lines, but also given the frequency of programmes that exploit the NIRCam medium-band filters around 4 μm. All the galaxies in our final sample show a clear [O iii]4959,5007 excess in NIRCam imaging. We can therefore be confident that the emission line detected is indeed Lyα.

Following the previous step, to confirm that all components of each system are at a similar redshift, we again use 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 use all available ground- and space-based spectroscopic observations that cover the companions 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.76716

These galaxies constitute the subsample 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.76716 were made with the NIRSpec MSA, and the positioning of the shutter misses the companion; as such, spectroscopic constraints cannot be placed 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, but as a result of seeing 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 show EGSY8p68 is probably an AGN, given the presence of high-ionization emission lines and a broadening of the Hβ emission line17. Although 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 NIRSpec MSA flux is seven times fainter than that observed in the MOSFIRE spectrum57. Given that the MSA shutter size is notably smaller than the MOSFIRE 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; as such, this component may also host Lyα emission at the same redshift as the main component. Moreover, ref. 17 claims that the system is consistent with being involved in a major merger, and recent NIRSpec IFU observations of EGSY8p68 spectroscopically confirm that they are physical companions (Carniani et al., manuscript in preparation). However, because these results are not yet published and the data are not yet public, we continue our analysis of the redshift constraints assuming that we do not have this spectroscopic confirmation.

Despite the lack of spectroscopic constraints on this subsample, 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 [O iii]4959,5007 and Hβ emission, in the same filters, as can be seen in Extended Data Figs. 1 and 2. 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. Although the redshifts of these components are very confidently constrained to within this Δz = 1.4, further constraining their redshift is challenging without spectroscopy. However, we use the UV luminosity function from ref. 58 to estimate the expected number of galaxies within the volume to which the [O iii]4959,5007 emission redshift constraint constrains the companion galaxies. We note that the photometric redshift from BAGPIPES (reported in Table 1) is much more tightly constrained than that provided by the [O iii]4959,5007 boost alone, but 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× 3″ region, we use this as the area for this estimate and the [O iii]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 FoV that are consistent with Hβ and [O iii]4959,5007 producing the excess F444W flux and with observed UV magnitude is \(0.0{9}_{-0.04}^{+0.16}\) for JADES-GS+53.15682-27.76716, \(0.00{7}_{-0.002}^{+0.009}\) for CEERS-1027 and \(0.00{3}_{-0.001}^{+0.006}\) 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

JADES-GS-z7-LA is an extremely high equivalent width LAE (400 Å) discovered in JADES NIRSpec observations of galaxies in the GOODS-S field53. The presence of a photometric candidate companion galaxy to this LAE has already been reported53, 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 [O iii]4959,5007 emission is in the F410M filter, limiting the redshift to between z = 6.8 and z = 7.6, and find the expected number of galaxies to be \(0.1{3}_{-0.05}^{+0.14}\). 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 galaxies53, supporting the conclusion that these two galaxies are at the same redshift.

COSY and z7-13433

COSY, z7-13433 and their companions show excess fluxes in the F410M filter driven by [O iii]4959,5007 emission. SED-fitting of their companions strongly prefers a z > 7 solution, and the presence of [O iii]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.0{3}_{-0.01}^{+0.03}\) for COSY and \(0.0{4}_{-0.01}^{+0.03}\) for z7-13433, and thus we conclude that 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 Data Fig. 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 px, which is consistent with the 0.5″ offset between the two components in the slit. The Lyα emission coincident with the position of the companion is detected at 2.8σ and corresponds to z = 7.479, and 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−1. In the case of COSY, we expect a similar spatial offset between the two components, but the 0.7″ seen 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−1 between COSY-A and COSY-B. Although these Lyα detections are clearly tentative, the fact that they are coincident with the expected position in the 2D spectra of strong-photometric candidate galaxies, together with their proximity to confirmed LAEs, boosts their likelihood of being a true detection.

Finally, we note that recent NIRSpec IFU observations of COSY spectroscopically confirm its companion galaxy38. However, as with EGSY8p68, because the data is not yet public, we continue our analysis without this spectroscopic information.

z7-GSD-3811 and GSDY

Both galaxies in this sample are in the FoV of the FRESCO NIRCam WFSS survey. Their final 2D WFSS spectra (Extended Data Fig. 4) show [O iii]4959,5007 emission lines for z7-GSD-3811, GSDY and their companions. For z7-GSD-3811, both components of the [O iii]4959,5007 emission are present, and there is an 3σ faint detection of Hβ. Moreover, [O iii]4959,5007 emission is also detected in the companion galaxy. This emission is offset from the central galaxy by 170 ± 70 km s−1. The 2D spectrum of GSDY shows a clear detection of [O iii]5007 in both the central and companion galaxies, and there is a very faint indication of [O iii]4959 and Hβ in the companion galaxy, but at a less than 3σ confidence for both. The emissions detected in the central and companion galaxies are offset by just 26 km s−1, but given the low spectral resolution of the NIRCam WFSS, this is smaller than the 70 km s1 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 Te = 10,000 K and Ne = 104 cm−3, and following ref. 59, we estimate the instantaneous SFR reported in Table 2. Given the lack of constraints on the reddening of these galaxies, we do not dust-correct the Hβ flux, and as such, these measured instantaneous SFRs are lower bounds.

The observation of several spatially resolved components of 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 LOS separation or, more likely, a difference in the LOS velocities of the two objects as they interact with each other. We therefore consider this further evidence that the systems are indeed interacting, as opposed to being two components of the same stable system.

GN-z11

Recent observations of GN-z11 with the NIRSpec IFU have shown the presence of three regions emitting either He ii (GN-z11-C)30 or Lyα emissions (GN-z11-B and a more tentative candidate)39. Although 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-C30 indicate that 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 several LAEs in the local vicinity of GN-z11 indicates that peculiar velocity is not playing a vital role in the escape of Lyα through the neutral IGM (see ref. 39 for a detailed analysis of the Lyα emission present in this system). The density of objects in this field also makes GN-z11 probably a protocluster core39.

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 Data Figs. 1 and 2. The best-fit model is then obtained by taking the SFH that leads to the smallest Bayesian information criterion, as described in ref. 60. We report the galaxy properties returned from this method in Table 1.

Close-companion fraction of non-LAEs

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}}}}}_{\odot })] < 9.3\) from refs. 23,61,62,63. 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 it is incredibly rare to have spectroscopic information on the local environment of such galaxies. Thus, although a galaxy may be confirmed as non-LAE, there is no information about 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 high-resolution G140M grating, that show no Lyα emission13. This does not provide any Lyα information on any would-be companion galaxies; however, the companion fraction for this sample, 43%, is again consistent with both our larger sample and that of Puskás et al. (manuscript in preparation). 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 system64 with an absence of detected Lyα emission. The reported specific SFR of this galaxy (approximately−8) is consistent with the high specific SFR 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. Although we see 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 although most viewing angles facilitate the escape of Lyα, several viewing angles exist for which the only Lyα emission escaping the galaxy is diffuse and therefore unobservable. Moreover, without the presence of a large ionized 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 reionization. A combination of a large-scale ionized bubble and a preferable LOS, stripped of neutral hydrogen by the interaction with the companion, are also required.

Azahar simulations

The new Azahar simulations are a suite of several models spanning various combinations of canonical hydrodynamics, magnetic fields, radiative transfer and cosmic ray physics with the aim of understanding their effects on galaxy formation as well as their complex interplay. The simulations are generated using the magnetohydrodynamical code Ramses65, which uses a constrained transport method for divergenceless evolution of the magnetic field66. The main target of study for Azahar is a massive spiral galaxy with Mhalo (z = 1) ≈ 2 × 1012M in a relatively large zoom-in region, about 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 mDM, where mDM ≈ 4.5 × 105M is the mass of dark-matter particles in the zoom region. The stellar mass particle is m* ≈ 4 × 104M. Azahar follows the set of physics presented by the pathfinder Pandora40. We provide here a brief summary, and we refer the reader to that reference and Martin-Alvarez et al. (manuscript in preparation) for further details.

In this work, we investigate one of the most complete simulations in the Azahar suite, the RTnsCRiMHD model. In this Article, we refer to the RTnsCRiMHD simulation simply as Azahar. Radiative cooling is modelled both above and below 104 K (refs. 67,68). The model adopts a magneto-thermo-turbulent prescription for star formation69,70, mechanical supernova (SN) feedback44 and astrophysically seeded magnetic fields through magnetized SN feedback71. This feedback injects 1% of the SN energy as magnetic energy and intends to reproduce the approximate magnetization 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 solver72,73. In Azahar, cosmic rays are exclusively sourced by SN feedback, with each event injecting 10% of its total energy as cosmic rays. We assume a constant diffusion coefficient for cosmic rays of κCR = 3 × 1028 cm2 s−1 and do not account for cosmic ray streaming in this model. This model also includes on-the-fly radiative transfer41, with a configuration similar to that of the SPHINX simulations42, featuring three energy bins for radiation spanning the ionization energy intervals 13.6–24.59 eV (H i ionization), 24.59–54.42 eV (HeI ionization) and above 54.42–eV (He ii ionization). 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 are lower than in their absence74.

Therefore, Azahar features a realistic (yet computationally taxing) ISM model that considerably approaches the resolution required to converge in the propagation and escape of ionizing photons from galaxies44. We note that although close, our resolution still falls short of that regime and that even more sophisticated ISM models75 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 ionized 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 ionized bubbles emerge at z ≈ 20. By z ≈ 15, several ionized 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 ionizing 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 ionized bubble may be sufficient for these sources to be detectable76. We further note that because of continuing 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 LOSs, hence requiring smaller local ionized 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 ionized bubbles could be even larger: for example, as a result of the presence of a neighbouring large-scale overdensity or large-scale clustering77 on spatial scales beyond our zoom-in region.

Nevertheless, knowing from our observations that the LAEs sit in sufficiently large ionized 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 ionized bubble, actively feeding the galaxies along the cosmic web.

Galaxy selection and measurements in the simulations

To select galaxies in our simulation, we use the Halomaker software78 to detect and characterize 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 component79, and they are assigned their corresponding halo as obtained by Halomaker. To select the system for study, we broadly reviewed the evolution of all galaxies with stellar masses M* > 107M 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.

To assign measurements to individual galaxies, we measure values within their galactic region, defined by the radius rgal < 0.2 rrvir. During the merging stages, we use r = min(rgal, 0.45Dij), where Dij 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 in the simulated domain. Although these present behaviour during mergers similar to that of our main studied system, their observability and that of other galaxies with similar masses remains low. This is because their low stellar masses generate low luminosities. Their low stellar masses also make them more vulnerable to disruption during powerful star formation bursts, which increases their escape fractions43. 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 RASCAS46 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 ref. 80 and the fraction of recombinations producing Lyα photons from ref. 81. For collisional excitation, RASCAS uses the fitting function for the collisional excitation rate from level 1s to 2p from ref. 82. We cast NMC = 106 Monte Carlo photon packets to sample the real Lyα photon distributions from recombination and collisional excitation, respectively. We sample NMC = 107 photon packets for the images along the LOS in Extended Data Fig. 3.

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 refs. 83,84 for the scattering of Lyα photons around the 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 enough that they can have a long mean free path. To reduce the associated computation cost, RASCAS adopts a core-skipping mechanism85 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 ref. 47. Dust can either scatter or absorb Lyα photons. The probability of scattering is given by the dust albedo adust = 0.32, following ref. 86. The dust scattering with the Henyey–Greenstein phase function87 and the asymmetry parameter is set to g = 0.73, following ref. 86. We generate synthetic images and spectra with a peeling algorithm described in refs. 49,88,89 along 12 LOSs uniformly sampled using the healpix algorithm90.

Finally, we note that the dust distribution has a large effect on the Lyα radiative transfer as well as the [O iii]4959,5007 and Hβ emission used to infer the intrinsic Lyα emission from observations. In Extended Data Fig. 6, we show how reducing the amount of dust changes the observed Lyα images. We see that 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 centre 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, because the channels are cleared of local neutral hydrogen, our results about the boosted intrinsic Lyα emission in star formation bursts and enhanced Lyα escape fraction in mergers are robust to the assumed dust model.

We further consider the Lyα profiles escaping along the LOS shown in the upper panel of Extended Data Fig. 3 at z = 7.3 after applying a reasonable IGM damping wing from ref. 91. We note that the shape of the Lyα profile is very complex and sensitive to both the LOS92 and IGM attenuation93. 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 attenuation91 and compare the velocity offset and full-width at half-maximum (FWHM) to that reported in our observational comparison, EGSY8p68. We find that the observed Lyα profile17 is well matched with the simulated profile’s velocity offset of 200 km s−1 and FWHM of 200 km s−1 and shows a clear asymmetric line shape as observed.