A massive interacting galaxy 510 million years after the Big Bang

Springer Nature 2021 LATEX template 2 A massive interacting galaxy 510 million years after the Big Bang A massive interacting galaxy 510 million years after the Big Bang Kristan Boyett1,2*, Michele Trenti1,2*, Nicha Leethochawalit3 , Antonello Calabr´o4 , Benjamin Metha1,2,5, Guido Roberts-Borsani5 , Nicol´o Dalmasso1,2, Lilan Yang6 , Paola Santini4 , Tommaso Treu5 , Tucker Jones7 , Alaina Henry8,9, Charlotte A. Mason10,11, Takahiro Morishita12, Themiya Nanayakkara13, Namrata Roy9 , Xin Wang14,15,16, Adriano Fontana4 , Emiliano Merlin4 , Marco Castellano4 , Diego Paris4 , Maruˇsa Bradaˇc17,7, Matt Malkan5 , Danilo Marchesini18, Sara Mascia4 , Karl Glazebrook13, Laura Pentericci4 , Eros Vanzella19 and Benedetta Vulcani20 1*School of Physics, University of Melbourne, Parkville 3010, VIC, Australia. 2ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia. 3National Astronomical Research Institute of Thailand (NARIT), Mae Rim, Chiang Mai, 50180, Thailand. 4 INAF Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Rome, Italy. 5Department of Physics and Astronomy, University of California, Los Angeles, 430 Portola Plaza, Los Angeles, CA 90095, USA. 6Kavli Institute for the Physics and Mathematics of the Universe, The University of Tokyo, Kashiwa, Japan 277-8583. 7Department of Physics and Astronomy, University of California Davis, 1 Shields Avenue, Davis, CA 95616, USA. 8Space Telescope Science Institute, 3700 San Martin Drive, Baltimore MD, 21218. 9Center for Astrophysical Sciences, Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD, 21218. 10Cosmic Dawn Center (DAWN), Denmark. Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 3 11Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark. 12IPAC, California Institute of Technology, MC 314-6, 1200 E. California Boulevard, Pasadena, CA 91125, USA. 13Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia. 14School of Astronomy and Space Science, University of Chinese Academy of Sciences (UCAS), Beijing 100049, China. 15National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China. 16Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal University, Beijing 102206, China. 17University of Ljubljana, Department of Mathematics and Physics, Jadranska ulica 19, SI-1000 Ljubljana, Slovenia. 18Physics and Astronomy Department, Tufts University, 574 Boston Avenue, Medford, MA 02155, USA. 19INAF – OAS, Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti 93/3, I-40129 Bologna, Italy. 20INAF Osservatorio Astronomico di Padova, vicolo dell’Osservatorio 5, 35122 Padova, Italy. *Corresponding author(s). E-mail(s): [email protected] ; [email protected] ; Abstract JWST observations spectroscopically confirmed the existence of galaxies as early as 300 million years after the Big Bang and with a higher number density than what was expected based on galaxy formation models and Hubble Space Telescope observations. Yet, the majority of sources confirmed spectroscopically so far in the first 500 million years have rest-frame UV-luminosities below the characteristic luminosity (M∗ UV ), limiting the signal to noise ratio for investigating substructure. Here, we present a high-resolution spectroscopic and spatially resolved study of a bright (MUV = −21.66 ± 0.03, ∼ 2M∗ UV ) galaxy at a redshift z = 9.3127 ± 0.0002 (510 million years after the Big Bang) with an estimated stellar mass of (1.6 +0.5 −0.4 ) × 109 M⊙, forming 19+5 −6 Solar masses per year and with a metallicity of about one tenth of Solar. The system has a morphology typically associated to two interacting galaxies, with a two-component main clump of very young stars (age less than 10 million years) surrounded by an extended stellar population (120 ± 20 million years old, identified from modeling of the Springer Nature 2021 LATEX template 4 A massive interacting galaxy 510 million years after the Big Bang NIRSpec spectrum) and an elongated clumpy tidal tail. The observations acquired at high spectral resolution identify oxygen, neon and hydrogen emission lines, as well as the Lyman break, where there is evidence of substantial absorption of Lyα. The [O ii] doublet is resolved spectrally, enabling an estimate of the electron number density and ionization parameter of the interstellar medium and showing higher densities and ionization than in analogs at lower redshifts. We identify evidence of absorption lines (silicon, carbon and iron), with low confidence individual detections but signal-to-noise ratio larger than 6 when stacked. These absorption features suggest that Lyα is damped by the interstellar and circumgalactic medium. Our observations provide evidence of rapid and efficient build up of mass and metals in the immediate aftermath of the Big Bang through mergers, demonstrating that massive galaxies with several billion stars are in place at early times. Keywords: galaxies: high-redshift, galaxies: formation, galaxies: interactions, galaxies: ISM 1 Main The first generations of stars and galaxies in the Universe formed in physical conditions different to those of the modern Universe. In fact, gas is expected to have nearly primordial composition, with low levels of chemical enrichment and dust content (e.g., [1]). Gas cooling is further limited by the higher Cosmic Microwave Background radiation, possibly altering the characteristic fragmentation mass of protostellar clouds [2]. In addition, early forming galaxies are expected to experience an elevated merger rate [3], affecting their morphology and stellar populations. Hubble Space Telescope observations have been able to identify galaxy candidates at redshift z ∼ 8−11 (∼ 600−400 million years after the Big Bang; e.g. see [4]), and follow-up with the Spitzer Space Telescope provided evidence of relatively old stellar populations, suggesting that star formation started at z > 15 [5, 6]. Yet, comparison to theoretical and numerical modeling has been restricted to number counts and luminosity functions by the limited angular resolution and dearth of spectroscopic data. With JWST commencing science operations in July 2022, progress has been rapid and transformational. Already in the first Cycle of science programs, JWST has built a convincing sample of z > 10 candidates based on NIRCam photometry [7–13], and a growing number of sources at z ≳ 8 are being confirmed spectroscopically with NIRSpec [14–27]. Surprisingly, the number density of (spectroscopically confirmed) high-redshift galaxies has been higher than expected by most models, in particular at the bright end of the luminosity function, possibly suggesting that we are missing key physical processes connected to the formation of first galaxies [28, 29], though observations do not suggest inconsistency with ΛCDM [27, 30, 31]. The superior spectral Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 5 and spatial resolution of JWST in the near-infrared offers a path to investigate these initial findings through detailed stellar population studies in early galaxies. An initial spectroscopic census of JWST galaxies above z > 7 paints a picture of elevated ionization parameters, low metallicities, low dust content, high specific star formation rates (sSFR), and potentially higher ionizing radiation escape fraction (fesc) than seen locally (e.g., [18, 22, 32–36]). Many of these galaxies show resemblance to the extreme interstellar medium (ISM) conditions observed in metal-poor actively star forming galaxies, both local (blueberries, greenpeas [37–40]) and at moderate redshifts (z ∼ 2, [41–43]). The spatial resolution of JWST also means the resolved stellar populations of galaxies in the epoch of reionization can be studied as well, revealing that these systems, once considered compact with HST, exhibit spatial variations in their physical properties [44–46]. In this work we extend the frontier of detailed investigations of individual galaxy properties at very high redshift by reporting on imaging and spectroscopic observations of one of the brightest among the galaxy candidates at z ≳ 9 observed with JWST, with flux ∼ 0.35 µJy in F444W (corresponding to mAB ∼ 25.0). These observations include 6-band NIRCam imaging (F115W, F150W, F200W, F277W, F356W, F444W, presented in Extended Data Table 1,) and NIRSpec high resolution (R ∼ 2700) Multi-object Spectroscopy (see 2.2). The galaxy - which we call here Gz9p3 - was initially identified as a potential F105W or F125W dropout based on HST observations, and then confirmed as a NIRCam F115W drop-out with a probable z ∼ 9.45 redshift [47]. The NIRSpec data we use are part of the GLASS-JWST program ERS-1324 [48] centered on the foreground galaxy cluster Abell 2744, which is gravitationally lensing the high-redshift background sources. Gz9p3 is located in the outskirts of the cluster and relatively far away from high-magnification regions, with a best estimate of the lensing magnification µ = 1.66±0.02 based on [49] (see 2.3). The JWST photometry for Gz9p3 is shown in the top panel of Figure

1. The source has a magnification-corrected apparent AB magnitude mAB = 25.56±0.06 in F444W, which for the cosmology adopted in this paper (see 2.1) corresponds to MAB = −21.77 ± 0.06 – approximately 50% brighter than the characteristic luminosity (M∗ ) of galaxies at this time [4]. The figure also includes a 1-D extraction of the NIRSpec spectrum in the middle panel, with spectral coverage in the [1.1 : 4.5] µm range (with some gaps due to the target location relative to the edge of the instrument field of view as discussed in 2.2.1). We present the full 2D spectrum in Extended Data Figure 1. The NIRSpec Multi-shutter Assembly (MSA) was configured based on HST imaging and as such the shutter covers the main body of the galaxy as shown in the bottom left panel of the figure. The spectrum shows a clear continuum detection and four emission lines are detected at [3.8 : 4.1] µm, which we identify as ([O ii], [Ne iii] λ3870, and [Ne iii] λ3969 blended with Hϵ) for a source Springer Nature 2021 LATEX template 6 A massive interacting galaxy 510 million years after the Big Bang [OII] [NeIII]�3870 [NeIII]�3969 + H� Fig. 1 JWST NIRCam and NIRSpec observations of Galaxy GLASS ID:Gz9p3. Top row: NIRCam direct imaging in broad-band filters. Second row: 1D standard extracted spectrum in units of fλ = 10−19erg s−1cm−2 ˚A−1 (with Npix=20 binning) from NIRSpec F100LP/G170H, F170LP/G235H & F290LP/G395H filter-dispersor configurations (from left-to-right in blue, orange and green). The observed frame location of the Lyman-break (λrest = 1215.67) and the [O ii], [Ne iii] λ3870 & [Ne iii] λ3969+Hϵ emission lines are overlaid as vertical dashed red lines for a zspec = 9.313. The Solid grey regions mask the location of contaminating emission lines from the dispersed spectrum of a galaxy falling in an open shutter of a separate quadrant of the NIRSpec MSA. Bottom left panel shows the color composite from NIRCam with NIRSpec slit positioning overlayed. The other three bottom panels are zoom in on a ±400˚A region centered on each emission line complexes ([O ii], [Ne iii] λ3870, [Ne iii] λ3969+Hϵ) with the continuum subtracted and the best-fit profile overlaid (again in units of 10−19 erg s−1 cm−2 ˚A−1 ). located at a redshift zspec = 9.3127 ± 0.0002 (see 2.5). In addition, the spectrum shows a Lyman break, and by modeling the stellar continuum as a step function around the break, we determine zbreak = 9.35+0.01 −0.05, consistent with the redshift measurement from emission lines (shown in Extended Data Figure 2, with further details in 2.5.1). The galaxy properties are derived from Spectral Energy Distribution (SED) modeling using both broad-band photometry and the 1-D spectrum as input (see Extended Data Figure 3, with further details in 2.8). The results are summarized in Table 1. Based on photometry of the whole galaxy (Kron fluxes, see [50] and 2.8), the modeling returns a magnification-corrected stellar mass of log10(M∗/M⊙) = 9.2 +0.1 −0.2 and MUV = −21.66 ± 0.03. (These are statistical uncertainties from photometric errors only). This makes Gz9p3 one of the most massive and intrinsically brightest galaxies confirmed in the epoch of reionization, and the brightest and most massive at z > 9 (see Figure 2 for a census for high-redshift spectroscopically confirmed galaxies). Even when compared Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 7 6 7 8 9 10 11 12 13 Spectroscopic redshift 22 21 20 19 18 17 M U V Wi22 JD1 GS-z10-0 GS-z11-0 GS-z12-0 GS-z13-0 Gz9p3 A23 CEERS_99715 CEERS_35590 MACS0647-JD GNZ11 CEERS2_588 Maisie’s galaxy CEERS2_7929 Prism Medium High Line emission This work 6 7 8 9 10 11 12 13 Spectroscopic redshift 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 S t ella r M a s s lo g 1 0(M

/ M ) Wi22 JD1 GS-z10-0 GS-z11-0 GS-z12-0 GS-z13-0 Gz9p3 CEERS_99715 GNZ11 CEERS2_588 Maisie’s galaxy CEERS_35590 MACS0647-JD CEERS2_7929 Fig. 2 Census of MUV and stellar mass in high-redshift galaxies. Left: Absolute UV magnitude of galaxies spectroscopically confirmed with JWST, corrected for lensing magnification where appropriate. Right: Stellar Mass (log10(M∗/M⊙)) distribution of confirmed galaxies, corrected for magnification where appropriate. Error bars derive from the measured mass or MUV of individual galaxies. In both panels, the point shape represents the resolution mode of the NIRSpec spectroscopy (low-resolution prism: circles; Medium resolution: squares; High resolution: diamonds and star). A black central dot marker indicates the detection of emission lines. Our target (red star with central dot) is one of the intrinsically brightest and most massive galaxies in the epoch of reionization among the current JWST samples from [15–27], and the highest and most massive at z > 9. It is also one of the highest redshift galaxies with emission line detections and the highest redshift one observed in the high resolution mode. Stellar masses were taken as quoted from each study. We note that we do not include MACS1149-JD1 [135] as a spectroscopically confirmed galaxy (at z=9.11) due to significant uncertainty on its lensing magnification, and hence also on its intrinsic MUV and Stellar Mass. Additionally, we note that only a subset of the [18] sample have their masses reported, and masses aren’t provided for [22]. We label all galaxies at z > 9. In these labels A23 refers to galaxy ID: 10058975 for [22] and Wi22 refers to [15]. against photometric galaxy candidates [51], Gz9p3 has one of largest masses known within the first 750Myr since the Big Bang. The SED modeling from the spectrum is restricted to the main region of the galaxy where the shutter was placed, and returns a stellar mass consistent with the photometric estimate (see 2.7). Both SED modeling approaches (photometry and spectrum fitting) also identify substantial ongoing star formation (9−19 M⊙yr−1 ), and limited evidence of dust due to the blue spectral slope β in the UV (−2.2 ≲ β ≲ −1.9), with robust results over a range of assumed star formation histories. Interestingly, the spectrum-based modeling shows evidence for older stellar populations in the central region of the galaxy (age 120±20 Myr), indicating that star formation started as early as z ≳ 15 to produce the average age observed (see 2.7). In contrast, modeling based on photometry infers younger ages, with integrated-light fits giving an age of 25+15 −12 Myr and spatially resolved analysis identifying regions with ages < 10 Myr (see 2.9.1 and Extended Data Table 2). Springer Nature 2021 LATEX template 8 A massive interacting galaxy 510 million years after the Big Bang The detection of rest-optical emission lines provides a window into the interstellar medium conditions in the galaxy by resolving the [O ii] doublet thanks to our high spectral resolution (Figure 1), measuring a line ratio of 0.94+0.14 −0.18 (see 2.7). The relative strength of these low-ionization lines is sensitive to the electron number density ne [52], leading to a measurement of ne = 590+570 −250cm−3 . This is marginally higher (at ≳ 1σ) than the median values of ne = 225cm−3 seen in galaxies at z = 2.3, and ne = 26cm−3 seen in local galaxies [52], qualitatively following the trend of ne increasing with redshift as reported in that study. From the spectrum we determine a Ne3O2 ([Ne iii] λ3870/[O ii]) ratio of 0.81 ± 0.09. As these two lines are close in wavelength, the ratio is insensitive to dust reddening. The measurement is higher than what is typically seen at lower redshift, indicating a high ionization parameter of log U = −2.13 ± 0.05 based on [53], and a low metallicity of 12 + log(O/H) = 7.6 ± 0.5 (depending on which one among the low-z Ne3O2 calibrations is used and including systematic uncertainties [54–57]; see 2.7.2 for further details). Together, these conditions indicate a sub-solar (Z ≲ 0.1Z⊙) metal-poor interstellar medium, with a high electron density and ionization parameter, exhibiting similar properties to other galaxies spectroscopically confirmed at zspec > 8 [19, 20, 22, 32–34, 36]. The ISM conditions are consistent with expectations from the young stellar ages derived from SED fitting, providing a self-consistent picture of the stellar populations and their surrounding gas. In Figure 3, the large stellar mass and low oxygen abundance place Gz9p3 below the mass-metallicity relations for z = 4 − 9 derived by [58], and marginally below the relation for z = 2 − 4 galaxies from [56, 59], even though systematic uncertainty may affect the robustness of these conclusions (see 2.7.2). The offset suggests Gz9p3 has a high gas fraction and potentially high accretion rates of pristine gas. This is qualitatively consistent with expectations from theoretical and numerical modeling of galaxies at these early times, given the short assembly times of their dark-matter halos [60]. Further insight into the stellar populations of Gz9p3 is offered by the detection - albeit at low confidence - of UV absorption lines, shown in Extended Data Figure 4. The spectrum shows a flux deficit associated with the SiIIλ1260, CIIλ1335 and FeIIλ2344 lines at a > 4σ significance individually, and at > 6σ when the lines are stacked together (see 2.7.7). The detection of metal absorption features reinforces the evidence from emission lines that the gas within the galaxy is not pristine and has been enriched by the older stellar population, again providing a consistent interpretation of the relatively evolved stellar ages inferred from SED modeling for the core of the galaxy where the slit is placed. The multi-component absorption profile hinted by the data suggests turbulence within the absorbing gas, as expected from models of star formation [61, 62], and/or from bulk gaseous flows that are associated to a merging or interacting system. However, higher signal to noise observations are required to quantify the robustness and confidence of these interpretations. Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 9 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 Stellar Mass log10(M * /M ) 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 1 2 +lo g 1 0(O / H) z 0 z 2.2 z 3.5 z 4 10 Maiolino+08 Andrews+13 Sanders+21 Nakajima+23 This work Fig. 3 Location of Gz9p3 on the Mass metallicity relation. Gz9p3 is shown in black, where the error bar presents the random and systematic uncertainty in the metallicity, based on Ne3O2 diagnostic calibrations from [54–57, 90]. The figure includes a comparison to mass metallicity relations covering 4 redshift epochs: z ∼ 0 (blue) from [56, 115], z ∼ 2.2 and z ∼ 3.5 (green and red) from [56, 59] and z = 4 − 10 (purple) from [58]. We additionally show the JWST z = 4 − 10 galaxies from [58] in purple, where the error bars derive from the mass and metallicity calculation for each individual galaxy. Interestingly, not only is there no evidence of Lyα emission, with a stringent limit on the equivalent width from Table 1, but also the stellar continuum at 1216˚A< λrest < 1240˚A (redward of the Lyman break) shows a deficit. Measurements over ∆λrest = 4˚A and 24˚A windows show a 80% and 40% deficit at a 5.7σ and 6.5σ significance, respectively (see 2.7.6). The softening of the spectral break supports the presence of absorption from Lyα damping wings (seen in the spectra of many z > 9 galaxies, [24, 63, 64]). One interpretation is that the damping is due to absorption by the intergalactic medium [65], which would indicate that the galaxy does not reside within a large ionized bubble. Such a scenario falls in line with the expected transmission due to damping wing in a neutral IGM from [66], with a predicted flux of ∼ 0.3× continuum and ∼ 0.8× continuum at 1000 km s−1 (∼ 4˚A) and 6000 km s−1 (∼ 24˚A) respectively. However, it is difficult to reconcile the lack of a large ionizing bubble with the high stellar mass and presence of relatively old (> 100 Springer Nature 2021 LATEX template 10 A massive interacting galaxy 510 million years after the Big Bang Myr) stars, especially because Lyα emission is detected in galaxies at higher redshift with lower star formation rates and stellar masses such as GN-z11 [23]. An alternative interpretation is that the interstellar and circumgalactic medium in Gz9p3 is primarily responsible for the lack of Lyα emission in the spectrum, irrespective of the IGM conditions. This scenario is supported by the detection of the SiIIλ1260 and CIIλ1335 absorption features with a restframe equivalent width of (−3.7 ± 0.8)˚A. In fact, it has been shown that their strength correlates with damping of Lyα and results in Lyα absorption when the equivalent width of low-ionization interstellar metal lines is ≲ −2˚A [67– 69]. Also, the presence of strong low-ionization ISM line absorption and the stringent upper limit on CIII] emission suggest that the galaxy is unlikely to be a Lyman continuum leaker [70]. In addition to providing detailed spectral insight on the stellar populations within the shutter aperture, the NIRSpec observations allow us to fix the redshift of Gz9p3 for photometric pixel-by-pixel modeling. This allows us to investigate spatial variations across the galaxy, following an approach similar to the one adopted by [46] at lower redshift, to create a 2D distribution of the galaxy’s physical properties (see 2.9). The analysis is carried out with BAGPIPES [71] to determine the following maps of resolved properties, shown in Figure 4: stellar mass surface density (SMD); 100Myr-time averaged star formation rate surface density (SFRD); mass-weighted stellar age; visual extinction (Av); and UV β slope (where fλ ∝ λ −β ). The native pixel resolution in F444W is 0.031”/pixel (with a FWHM F444W ∼ 4.5 native pixels, or ∼ 0.14”), but for our analysis we bin pixels 2 × 2 to improve the signal-tonoise ratio per pixel, generating the maps at a 0.27kpc/pixel resolution in the observed frame. This corresponds to 0.21kpc/pixel in the image plane after accounting for gravitational magnification. The galaxy shows a morphology comprised of an elongated tail and a central body, which in F444W appears as a single core. The peak of the star forming activity within the stellar population is situated in the central core and this traces the distribution of stellar mass. This central region of active star formation exhibits young stellar population ages (< 50Myr), as does the clump within the tail, whilst the surrounding regions on the outskirts of the central body show older populations, albeit with larger uncertainties (50 ≲ tage/Myr ≲ 125). The galaxy shows blue β slopes throughout and relatively low visual extinction, which implies low dust content (see 2.9). The spatially resolved modeling is suggestive of an interacting system undergoing (or having recently undergone) a major merger. To further investigate this scenario, we analyze photometry at rest-frame UV wavelengths, using a combined F150W+F200W image drizzled at 20 mas/pixel, presented in Figure 5. The data clearly show two distinct cores in the main region of the galaxy, and several components in the tail identified thorough a clump-finding algorithm (see 2.10). The morphology of Gz9p3 is described by morphological parameters that indicate the galaxy is a merger (Gini=0.61, M20=-1.29, A=0.35; see [72]). Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 11 SFRD 0.0 2.5 5.0 Age 0 20 40 60 SMD 5 6 7 SFRD 0.0 2.5 5.0 M yr 1kpc 2 Age 0 20 40 60 Myr SMD 0 1 2 log M kpc 2 F150W-F444W 0.5 0.0 0.5 1.0 Av 0.0 0.5 1.0 UV 1.5 2.0 2.5 F356W-F444W FWHM 0.27kpc/px 0.5 0.0 0.5 1.0 mag Av 0.0 0.5 1.0 Av UV 0.0 0.5 1.0 Fig. 4 2D color and physical parameter distribution of Gz9p3. Properties inferred from photometric spectral energy distribution fitting (NIRCam pixels matched to F444W and binned 2x2), with an observed frame resolution of 0.27kpc/pixel (0.21kpc/pixel in the image plane). From left to right: the top and second row present the star formation rate surface density, stellar age, and stellar mass surface density. The third and bottom row present the color, visual extinction and UV β slope. (Upper panels: Median value, Lower panels: Uncertainty based on 16th and 84th percentiles). The F150W 10σ contour is presented in orange in the F150W-F444W panel and the FWHM and pixel-scale are shown in the F356W-F444W panel. Springer Nature 2021 LATEX template 12 A massive interacting galaxy 510 million years after the Big Bang F150W+F200W ClumpMap Fig. 5 Morphology of Gz9p3. F150W+F200W direct image at a 20mas/px resolution in both panels. Left: direct imaging shows a double core within the central region and an elongated clumpy structure. Right: Overlaid Clump-map, showing 4 clumps detected within the tail of the system with a clumpiness parameter c = 0.56. Informed by the morphological analysis, we repeat the spatially resolved SED modeling by placing apertures over different stellar populations shown in Extended Data Figure 5 (innermost region, a surrounding annulus and the tail), clearly seeing a distinction between active regions of star formation and an underlying older stellar population, as expected from the merging scenario (see 2.9.1). The combined spectroscopic and imaging data paint a picture of a very bright and relatively massive interacting system just 0.5 Gyr after the Big Bang, raising the question of how likely such JWST observations should be. Figures 2-3 hint that the system could be an outlier. To quantify expectations we consider both analytical modeling of early galaxy formation and comparison to cosmological hydrodynamical simulations. We find that while the likelihood of capturing an interacting system is relatively high, i.e. ∼ 20% for a major merger, the stellar mass of Gz9p3 is higher than expected (see 2.11-2.12). This would indicate either the system is hosted in a very rare dark matter halo for that epoch, that we serendipitously observed, or more likely that the current recipes for star formation are missing some key ingredients at early times[73]. The latter interpretation would be consistent with the excess of sources identified by JWST at z > 10 through imaging programs [47] and with the high numbers of massive red galaxy candidates found at z ∼ 7.5 − 9 [51]. All these aspects make Gz9p3 an excellent target for further spectroscopic investigations, in particular through the Integral Field Unit mode on NIRSpec, that would shed further light on the kinematics of the system and on the complex interplay between assembly of dark matter halos, star formation and physical conditions in the interstellar, circumgalactic and intergalactic media at very early times. Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 13 Property Observed Value RA [Deg] 3.617193 DEC [Deg] -30.4255352 z a spec 9.3127 ± 0.0002 µ 1.66 ± 0.02 Line flux [10−19erg s−1cm−2 ] Lyα < 2.64† CIII]λ1908 < 1.45† MgII λ2804 < 1.07† [O ii] 6.7 +0.4 −0.5 [Ne iii] λ3870 5.4 +0.6 −0.5 ([Ne iii] λ3969 + Hϵ) 4.5 +0.6 −0.4 [Ne iii] λ3969 1.1 +1.7 −0.9 Hϵ 3.4 +1.1 −1.5 Line EWrest [˚A] Lyα < 7.6 † CIII]λ1908 < 1.0 † MgII λ2804 < 1.9 † [O ii] 25.8 +1.5 −1.9 [Ne iii] λ3870 21.4 +2.4 −2.0 ([Ne iii] λ3969 + Hϵ) 18.0 +2.4 −1.6 [Ne iii] λ3969 4.4 +6.7 −3.6 Hϵ 13.5 +4.4 −6.0 Full photometry SED fitb Stellar Mass [log10(M∗/M⊙)] 9.2 +0.1 −0.2 SFR† [M⊙yr−1 ] 19+5 −6 Stellar Age [Myr] 25+15 −12 β −1.94+0.05 −0.06 MUV,SED [AB mag] −21.66 ± 0.03 Spectrum+photometry SED fit of main componentc Stellar Mass [log10(M∗/M⊙)] 9.15 ± 0.04 SFR† [M⊙yr−1 ] 9.1 ± 0.6 Stellar Age [Myr] 120 ± 20 β −2.23 ± 0.04 MUV,SED [AB mag] −20.92 ± 0.02 Table 1 Physical properties for galaxy Gz9p3. †1σ upper limit. Line fluxes are not corrected for dust extinction or slit loss. b Photometry-only SED properties for the full system, corrected where appropriate for magnification. c Spectrum+Photometry SED properties for our main region aperture (see Extended Data Figure 5), corrected for magnification and adopting slit-losses based on aperture photometry. † The SFR is taken as the 100Myr average from the BAGPIPES SFH model. 2 Methods 2.1 Cosmology and conventions Where applicable, we use a standard ΛCDM cosmology with parameters H 0 = 70 km s−1 Mpc−1 , Ωm =0.3, and Ω∧ =0.7. All magnitudes are in the AB system [74]. Throughout this paper we refer to quantities as “observed” or “intrinsic” depending on whether the best available estimate of any gravitational lensing Springer Nature 2021 LATEX template 14 A massive interacting galaxy 510 million years after the Big Bang magnification correction has been applied. All wavelengths for emission lines are quoted as vacuum wavelengths. 2.2 Observations This work is based on JWST/NIRCam [75] direct imaging and JWST/NIRSpec [76, 77] high resolution multi-object spectroscopy as part of the GLASS-JWST survey (ERS 1324, PI Treu; [48]) and DDT program 2756 (PI Wenlei Chen). We refer the reader to [50, 78] for details of the observations and reduction strategy for NIRCam and to [14] for NIRSpec. As part of the GLASS JWST/NIRSpec observations, over 100 galaxies were assigned to slit apertures. Our target (GLASS ID: 10003, indicated as Gz9p3 hereafter) was included as a candidate, based on HST ACS/WFC3 photometry [79] which identified it as a potential F105W/F125W drop-out. This placed the object at a redshift zphot > 8. Without longer wavelength imaging at the spatial resolution needed to deblend near-proximity neighbors, the validity of the single broadband filter detection in HST/WFC3 F160W could not be confirmed and physical properties could not be determined. The photometric redshift solution has since been refined to zphot = 9.45 based on NIRCam observations [ID: DHZ1 in 47]. This object was one of 7 z ∼ 10 galaxies identified in the ABELL 2744 cluster region, suggesting an apparent over-density at this epoch in the field [see 47] and reaffirming the priority for spectroscopic follow-up. JWST NIRSpec observations now provide the spectral coverage needed to investigate this object further. 2.2.1 High-resolution spectroscopy We obtain high-resolution (R ∼ ∆λ λ = 2700) JWST/NIRSpec Multi-Object Spectroscopy of this target using the F100LP/G140H, F170LP/G235H and F290LP/G395H filter-disperser combinations. These spectroscopy configurations combined with the location of the object on the detector provide wavelength coverage from 1.00 − 1.48µm, 1.70 − 2.49µm and 2.88 − 4.19µm, respectively (970-4060˚A in the rest-frame, see zspec discussion in 2.5). The exposure time for each configuration is 4.9 hours. We note that while this work was undergoing peer-review and after an initial version of this manuscript had been posted on arXiv, observations of Gz9p3 using NIRSpec in the prism mode have been conducted by the UNCOVER team ([80, 81]) to extend the observed wavelength coverage out to ∼ 5.2µm, additionally detecting [O iii]+Hβ. In Figure 1, we overlay the location of the NIRSpec open shutters on the RGB direct image. The open shutters are placed over the main component but do not capture the total flux from the galaxy and therefore suffer from some slit losses. Due to the extended morphology and potential change in stellar populations across the galaxy, we do not correct the slit-losses to match the total photometry of the galaxy. Instead we adopt a slit-loss correction based only on the photometry of the main component (where the slit is placed), as part of our SED spectrum+photometry analysis of this region (see 2.7). Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 15 -1.5 -1.0 -0.5 0.0 0.5 1.0 yscale (arcsec) 1.1 1.2 1.3 1.4 -0.1 0 0.1 0.2 0.3 F Lyman-break 1.8 1.9 2.0 2.1 2.2 2.3 2.4 Wavelength - m 3.0 3.2 3.4 3.6 3.8 4.0 [OII] [NeIII] [N eIII] + H Extended Data Figure. 1 NIRSpec 2D spectrum of Gz9p3. Top: 2D observed-frame high resolution R ∼ 2700 spectrum in the f100lp/g140h, f170lp/g235h and f290lp/g395h filter-disperser configurations. Orange and white horizontal lines show the 1σ extraction trace for the optimal and narrow kernels. In the 2D extraction, the relative proximity of the dispersed light from other sources can be seen. These spectra are not associated with an additional source within our shutter but rather with targets in separate shutters in the NIRSpec MSA with a similar row number. The narrow kernel is introduced to ensure no contamination from the dispersed light of close proximity spectra is included for the fitting of the Lyman break. Bottom: the 1D extracted spectrum (as in Figure 1). We mark the location of the Lyman-break and detected emission lines in red in both the 1D and 2D spectra. The wavelength region contaminated by the additional source is marked in grey in the 1D. Therefore, for measurement of physical properties for the whole galaxy (e.g., Mass, SFR in 2.8), we rely on the photometric SED fitting (see [50]), and use the spectrum SED fitting as consistency check and for investigation of the properties of the core of the system. A detailed discussion of how the data were reduced is given in [14]. Briefly, we use the official STScI JWST pipeline (ver.1.8.2) for Level 1 data products, and the msaexp Python package for Level 2 and 3 data products. The orientation of the open-shutter configuration relative to Gz9p3 means the core of the galaxy is well-contained within the middle shutter, with the adjacent open shutters unoccupied, and so only sample the background. This orientation allows the adjacent open-shutters to be used for background subtraction and will not suffer from auto-subtraction. We extract the 1D spectrum following [14] and use the msaexp package to optimise our extraction using an inverse-variance weighted kernel. This kernel is derived by summing the 2D spectrum along the dispersion axis and fitting the resulting signal along the spatial axis with a Gaussian profile (σ ∼ 0.4”). This is then used to extract the 1D spectrum along the dispersion axis. Springer Nature 2021 LATEX template 16 A massive interacting galaxy 510 million years after the Big Bang We present the full 2D spectrum in Extended Data Figure 1 and the extracted 1D spectrum in Figure 1. Within the 2D image there is clear emission from an additional source below the trace associated to our target; this has been identified as a z ∼ 3 galaxy with prominent emission lines (λobs ∼ 2µm) and a faint continuum. Over the wavelength range 3-3.75µm the average continuum flux density of Gz9p3 is greater than that of the neighbouring source by a factor of 5 over the same window. The dispersed light from this second galaxy, which is another GLASS target, comes from an open shutter in a separate quadrant of the NIRSpec MSA. The dispersed light is spatially offset and the continuum is sufficiently faint as to not significantly contaminate our target, this can clearly be seen below the Lyman break where the continuum is consistent with zero (i.e. no contaminating emission). The weight of the neighbour’s continuum emission is lower in the optimal extraction, reducing its influence even more. However, the strong rest-optical line emission does leave a negative imprint on the continuum trace of our object due to our background subtraction process. These contaminated regions are masked in Figure 1 and when modeling the spectrum of our object. The close proximity of the dispersed light from other open shutters demonstrates the potential dangers of over filling the NIRSpec MSA. In the 2D and 1D spectrum of Gz9p3, several emission lines are visible at ∼ 4µm and we present a close up of this region in Figure 1 with individual panels focusing on each of the detected emission lines. In these panels it can be seen that the R ∼ 2700 resolving power can deblend the [O ii] doublet into the separate λrest = 3727 and 3730˚A lines. Finally, we perform a second extraction of the F100LP/G170H 2D spectrum with the intention of maximizing the signal to noise in the continuum around the observed spectral break (λobs ∈ 1.2 : 1.3µm), which we model in 2.5. For this purpose we restrict the width of the Gaussian profile (σ = 0.1”, ∼ 0.5px) to reduce the inclusion of pixels that contain only noise. 2.2.2 Direct Imaging Our spectroscopy is complemented by direct JWST/NIRCam imaging in the F115W, F150W, F200W, F277W, F356W and F444W broadband filters obtained by the DDT program 2756 (PI W. Chen), and reduced and discussed in [50]. We additionally use the HST ACS and WFC3 F606W, F814W, F105W, F125W and F160W broadband filters from [79]. All the direct imaging has been PSF-matched to the NIRCam F444W broadband filter, the longest wavelength filter with the coarsest resolution, using the ABELL 2744 region UNCOVER PSF models (epoch: 2022/11/07, PA: 41.2 deg, [50, 80]). We report the Kron flux measurements for each broadband filter in Extended Data Table 1, based on the analysis by [50]. In Figure 1, we present the direct imaging of the target in each of the available JWST filters. Here it can be seen that the galaxy exhibits a morphology comprised of a main body plus an extended tail. We examine the SED of the full broadband photometry in 2.8 and for specific regions within the galaxy using aperture photometry in 2.9.1. Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 17 Filter Fν (µJy) F606W 0.01 ± 0.02 F814W 0.02 ± 0.11 F105W −0.02 ± 0.04 F125W 0.07 ± 0.03 F160W 0.31 ± 0.05 F115W 0.03 ± 0.01 F150W 0.32 ± 0.01 F200W 0.34 ± 0.01 F277W 0.33 ± 0.01 F356W 0.32 ± 0.01 F444W 0.35 ± 0.02 Extended Data Table. 1 Observed broadband flux density (zero point: zp AB=23.9). HST photometry is presented above the horizontal line and JWST photometry is presented below. 2.3 Gravitational lensing Gz9p3 lies within the ABELL 2744 cluster region, albeit at a relatively large distance from the center of the cluster. Nonetheless this target will experience some degree of magnification. In this study we account for the cluster-induced magnification factor at the location and redshift (see Section 2.5) of Gz9p3 by adopting a magnification µ = 1.66±0.02 when appropriate, as in [47]. The magnification is derived from the latest strong lensing model of A2744 described in [82]. This model exploits 149 multiple images (121 with a spectroscopic confirmation at 1.03 ≤ z ≤ 9.76) over an area of 30 arcmin2 and takes advantage of newly identified multiple images based on NIRCam multiband photometry and Very Large Telescope(VLT)/MUSE spectroscopy. We refer to [82] for a detailed description of the modeling of the mass distribution, as well as for the catalogue of the multiple images. The magnification is sufficient to boost the observed magnitude, however is too low to affect the observed morphology of the galaxy, while the object size is magnified by ∼ 30% (p (µ). 2.4 Photometric properties This galaxy has a magnification corrected F444W magnitude of 25.56 ± 0.06. Under the adopted cosmology the distance modulus is 49.9, and using a kcorrection of −2.5log10(1 + z), the object has an absolute magnitude MAB = −21.77 ± 0.06. This brightness makes it comparable to the rare bright z ∼ 9 galaxies discovered with HST photometry in the large-area BoRG survey [e.g., 83]. The rest-frame 1500˚A is covered by the NIRCam/F150W broadband filter (with no detected emission lines that contribute to the flux density) and provides a magnification corrected estimate of MUV = −21.67 ± 0.04 using the same k-correction as above. These values place this object as one of the intrinsically brightest galaxies spectroscopically confirmed in the early Universe. We plot this object onto the growing sample of known galaxies at high redshift in MUV - redshift space in Springer Nature 2021 LATEX template 18 A massive interacting galaxy 510 million years after the Big Bang Figure 2. At the other luminosity extreme among z > 9 sources is a highly magnified z = 9.76 galaxy also from the GLASS-JWST survey [16] with an intrinsic MUV = −17.35 (µ ∼ 13), demonstrating the dynamical range capabilities of the GLASS program and of JWST in general. 2.5 Emission line detection and spectroscopic redshift determination The spectrum of this galaxy shows multiple emission lines and a spectral break in the UV, and both features enable the determination of a spectroscopic redshift. We detect emission from [O ii], [Ne iii] λ3870, [Ne iii] λ3969, Hϵ and we measure the corresponding line flux without correcting for slit-losses using the latest version of the specutils packages in Python. The stellar continuum in the high-resolution spectroscopy is well modeled by a simple polynomial after masking regions with detected emission line contribution or contamination from the dispersed light of neighboring targets) with a reduced Chi-square = 1.2. The continuum is subtracted before modeling the emission lines. Note that we neglect possible underlying absorption from stellar Balmer absorption lines, which could affect the Hϵ line because we use this line uniquely for redshift determination, which is not affected by any reduction in the equivalent width. We adopt a bootstrap method to represent the flux of each emission line, we model [Ne iii] λ3870 with a single Gaussian profile, and the [O ii] doublet and the [Ne iii] λ3969+Hϵ complex with a pair of Gaussians. The only constraint we place is that within the [O ii] doublet both emission lines share the same standard deviation (σ). For n = 10, 000 individual computations we vary each spectrum by applying random gaussian noise to each pixel following the associated uncertainty spectrum and fit gaussian profiles to each identified emission line. The n fits provide a resultant distribution of measured line fluxes from which we determine the median line flux as the 50th percentile and the associated uncertainty from the 16th and 84th percentiles. We report the measured line fluxes, which are not corrected for slit losses or reddening, in Table 1. We likewise use the bootstrap Monte Carlo to determine the best-fit redshift. We record λobs for each of the n Gaussian profile fits for the [Ne iii] λ3870 emission line and from the [16th, 50th, 84th] percentile we determine the best fit redshift and 1σ uncertainty to be zspec = 9.3127 ± 0.0002. The choice of using [Ne iii] λ3870 was made because it avoids the partial-blending in the [O ii] doublet and the [Ne iii] λ3969+Hϵ complex. 2.5.1 Spectroscopic redshift from the Lyman break We can also measure the redshift from the location of the Lyman break. In the f100lp/g140h spectrum, we attribute the observed spectral break feature (range λ ∈ [1.2, 1.3]µm) to be the Lyman break at λ = 1215.67˚A. To determine the best fit to the spectral break, we use the narrow 1D extracted spectrum Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 19 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 Wavelength ( m) 0.050 0.025 0.000 0.025 0.050 0.075 0.100 0.125 0.150 f ( J y) binned data (20px) raw data best fit fixed z 95% interval Extended Data Figure. 2 NIRSpec constraints on the Lyman break of Gz9p3. We present a best fit and 95% confidence interval of zbreak = 9.35+0.01 −0.05 using the narrow 1D extraction to minimize potential contamination (20px binning presented in blue at the midpoint of each bin, and raw data in gray). The orange region covers the 95% confidence interval with the solid black line showing the best fit model. The dashed lines shows the model fit using zspec = 9.313 derived from the emission lines. The independent best-fit of the Lyman break is consistent with the emission-line spectroscopic redshift solution. (to maximize the S/N) and randomize the data following the uncertainty spectrum by conducting Monte Carlo sampling a total of n = 10, 000 times. We model the Lyman break as a step function and determine a best-fit observed wavelength of 1.258+0.001 −0.006µm, with the uncertainty assigned by identifying the 95% confidence region when minimizing the χ 2 . This best-fit wavelength corresponds to a redshift of zbreak = 9.35+0.01 −0.05 (in full agreement with the redshift determination from the emission lines). The best fit and uncertainty regions are shown in Extended Data Figure 2. We note that while this break redshift is consistent with the line redshift within the errors, the redshift determined from the break may be systematically slightly higher because of the effect of the Lyα damping wing (see Section 2.7.6), since we do not consider when we model the spectrum with a step function break. 2.6 Constraints on other emission lines It is also of interest which emission lines fall within our wavelength coverage but are not detected. In the rest-optical we do not detect either HeI lines at λrest = 3890˚A or 4025˚A, placing 1σ upper limits of 0.50 and 0.39 ×10−19 erg s −1 cm−2˚A−1 based on the noise properties of the spectrum over a λrest ± 4˚A window (the extent of the pixel coverage of the detected emission lines). Neither do we detect any of the following commonly-studied rest-UV emission lines: Lyα, CIII](λ1907, 1909), MgII(λ2804), [NeV] (λ3427) which fall in our coverage (CIVλ1548, 1551, OIII]λ1661, 1666, HeIIλ1640, MgIIλ2799 are all out of our coverage). For each of these we place 1σ upper limits of [2.64, 1.45, 1.07, 0.54] ×10−19 erg s−1 cm−2˚A−1 . We note that many of these are intrinsically faint lines relatively to rest-optical [O ii]. Springer Nature 2021 LATEX template 20 A massive interacting galaxy 510 million years after the Big Bang 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Wavelength ( m) 0.0 0.1 0.2 0.3 0.4 F ( Jy) Extended Data Figure. 3 Spectral energy distribution fit to the spectrum and photometry of the central region of Gz9p3. Scaled-spectrum and aperture photometry of the main component of Gz9p3 with best fit BAGPIPES (Spec+Phot) model overlaid (red) on top of the spectrum (both binned at 20px) and the photometric flux densities (where the error bars indicate the filter width and uncertainty on the broadband imaging flux density). Wavelength regions affected by contamination are masked and shaded in gray. Each spectrum has been corrected for slit-losses to match the aperture photometry of the main component (see Extended Data Figure 5) in the F150W, F200W and F356W bands respectively. 2.7 Physical interpretation of the spectrum The detection of emission lines and the shape of the stellar continuum allows us to infer further galaxy properties of Gz9p3. Although, we note that the location of the MSA open shutter (as shown in Figure 1) means we capture the light from the core of the system only. 2.7.1 Spectrum SED fitting We utilize BAGPIPES [71] to fit a model SED provided both the spectrum and the photometry to determine the physical properties of the galaxy. For the combined analysis we apply a conversion factor to account for the slit-losses for each of the spectra such that the spectral flux density matches the broadband photometry in the F150W, F200W and F356W filters, respectively for each of the three non overlapping spectral regions covered by the data. Further discussion to validate this correction based on the resolved stellar population is discussed in 2.9. The SED fitting procedure is described in 2.8. Extended Data Figure 3 presents the best fit template overlaid onto the adjusted spectrum and photometry, and the best-fit properties are given in Table 1. We correct our measurements for the lensing magnification; (see 2.3) and we report that the core of the target covered by the shutter has an intrinsic log10(M∗/M⊙) = 9.15 ± 0.04 and a 100Myr time-averaged SFR of 9.1 ± 0.6M⊙/yr. The core also shows minimal visual extinction with Av = 0.14 ± 0.04, equivalent to E(B − V ) = 0.57 ± 0.20 for a Rv = 4.05 [84]. The mass-weighted stellar age of 120 ± 20Myr reflects the presence of an older stellar population and indicates that star formation started as early as z ≳ 15 to produce the average age observed. Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 21 2.7.2 Metallicity The detection of both [O ii] and [Ne iii] λ3870 allows us to study both the number density with the galaxy’s ISM, via the ratio of the [O ii] doublet line fluxes, and the ionization parameter and metallicity from the flux ratio of the [O ii] and [Ne iii] λ3870 lines. We measure a Ne3O2 (defined as [Ne iii] λ3870/[O ii]) line flux ratio of 0.81 ± 0.09. Considering how close these two lines are in wavelength, the measured ratio has only minimal sensitivity to either dust reddening or slitloss correction, which we therefore neglect. We determine a metallicity estimate for the galaxy from the measured Ne3O2 line ratio (where we take solar metallicity to be 12 + log10(O/H)=8.69, [85]). The Ne3O2 diagnostic exhibits an anti-correlation with Oxygen abundance [54, 55, 57, 86, 87] due to its relation with the ionization parameter (which is anti-correlated with metallicity, [88, 89]). Many studies have modeled the relation between the Ne3O2 diagnostic and oxygen abundance in samples of local galaxies. Using diagnostic calibrations from literature we determine metallicities of: 7.87 ± 0.02 [57] 7.86 ± 0.03 [54], 7.38 ± 0.07 [56], 7.17 ± 0.06 [55], with a mean and standard error from these diagnostics of 7.6 ± 0.2. Note that with the data available, we are unable to break the degeneracy in the non-monotonic behavior of the [90] relation and obtain two possible values of 7.01 ± 0.07 and 7.88 ± 0.07. Also, our Ne3O2 value appears to lie above the calibrated range of [59], placing an upper limit of < 7.4. All the values quoted here are subject to systematic uncertainties in the diagnostic calibration of the order 0.3dex, which dominate the random uncertainty. Our Ne3O2 line therefore supports an oxygen abundance of 10 % solar (12 + log10(O/H)= 7.6 ± 0.2 ± 0.3sys), in line with recent JWST studies of z ∼ 6−10 galaxies [19, 22, 36, 58]. However, an important caveat is that these metallicity results are based on diagnostics calibrated at low redshifts (SDSS z < 0.1). While criteria were used in these studies to select low-redshift EOR analogues, it has not yet been confirmed as to whether these relations still hold true for high-redshift galaxies at z ≳ 9. 2.7.3 Electron Density Following [91], we determine an electron density of ne = 590+570 −250cm−3 from a measured [O ii]λ3730/[O ii]λ3727 flux ratio of 0.94+0.14 −0.18. This has a large uncertainty as ne is highly sensitive to the flux ratio and the latter measurement has limited signal to noise. [91] report a median [O ii]λ3730/[O ii]λ3727 flux ratio of 1.18 in their z ∼ 2.3 sample (ne = 225cm−3 ). Gz9p3 shows stronger [O ii]λ3727 relative to [O ii]λ3730 than this lower redshift sample, indicating a higher electron density. 2.7.4 Ionization state The Ne3O2 ratio in this object is significantly higher than has been found in lower redshift samples. The stacked ratio from 66 MOSDEF galaxies of similar mass (log10(M∗/M⊙):8.23 − 9.51) at redshift z ∼ 2 is Ne3O2 = 0.19 Springer Nature 2021 LATEX template 22 A massive interacting galaxy 510 million years after the Big Bang [88], significantly lower than our value. This trend for a high Ne3O2 ratio has also been found in the recent JWST spectra of three z ∼ 8 galaxies (0.5 <Ne3O2< 2, [32, 33]) in the SMACS early-release observations (JWST program ID: 2736, [92]), and is indicative of a higher ionization parameter [89]. Using the [53] Ne3O2-U calibration, we determine an ionization parameter of log U = −2.13 ± 0.05. Due to a similar dependence on the ionization parameter, Ne3O2 can be used to trace the O32 ([O iii]/[O ii]) strong line diagnostic, with the advantage of being less sensitive to reddening [53, 57, 89, 93]. In a sample of z ∼ 0.8 galaxies, [57] demonstrate that the Ne3O2 diagnostic traces the dust-corrected O32, and can be reasonably approximated as [Ne iii] λ3870≈[O iii]/13, although we note scatter around this relation may hold a dependence on the hardness of the ionizing spectrum [94, 95]. Adopting this approximation would give Gz9p3 a dust-free O32 = 11.2 ± 1.3. This estimate is larger than the typical range of values found in local samples (O32 ∼ 0.3 − 4, for star forming galaxies selected from the Sloan Digital Sky Survey [96] Data-Release 7 [97]) and falls into the trend that star forming galaxy samples at higher redshift typically exhibit large O32 compared to the local Universe [e.g., 36, 94, 98]. Our estimate is comparable to recent JWST/NIRSpec observations of zspec ∼ 8 sources. We find consistency with the range of O32 values (O32∼ 8−20, [32, 33]) found in three SMACS galaxies, the range of O32 values (O32∼ 5 − 21, [19]) found in five GLASS galaxies, the lower O32 limits (O32> 5 − 8.2, [20]) placed on six galaxies from the Cosmic Evolution Early Release Science survey (CEERS, JWST Program ID: 1345, [99]), and the range of O32 values (O32∼ 6 − 30, [22]) found in five galaxies from the JWST Advanced Deep Extragalactic Survey (JADES, JWST program IDs: 1180, 1181, 1210, 1286 & 1287, [100]). The ISM conditions observed in Gz9p3 provide additional evidence for a higher ionization parameter at high redshift, likely due to the presence of a young stellar population and a harder ionizing spectrum, in line with the trend for an increase in the global specific star formation rate with redshift [101, 102]. 2.7.5 AGN constraints We report no evidence for Gz9p3 being an active galactic Nuclei (AGN). There is no broad-line component required in the model fitting of the detected emission lines, and we place an upper limit of 1.07×10−19 erg s−1 cm−2˚A−1 on the MgII line flux. This corresponds to an equivalent width limit EWrest < 1.9˚A. Based on a local sample of AGN, this limit would indicate mBH < 107 M⊙ and a bolometric luminosity of the source far exceeding the Eddington luminosity of the central BH [103]. The inference that an AGN is not the dominant source of light is also fully supported by the extended, spatially resolved nature of Gz9p3. Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 23 2.7.6 Lya absorption The shape of the Lyman break can be softened in the presence of a high neutral fraction in the Intergalactic medium (IGM) due to the strength of the Lyα damping wing, allowing the Gunn-Peterson absorption to extend redward of λrest = 1216˚A [65]. To assess whether there is any evidence for Lyα absorption we estimate the expected continuum level by fitting a fλ ∝ λ β model to the full stellar continuum in the narrow spectrum (to ensure we avoid any contamination from neighboring spectra) after masking the region of Lyα absorption, regions covering detected emission lines and regions of contamination from emission lines in close-proximity dispersed spectra. Based on the best-fit model we expect a continuum signal at λrest = 1216˚A of 1.01 ± 0.03 ×10−20 erg s −1 cm−2˚A−1 , which agrees with the mean signal over the wavelength range λrest = 1240˚A−1300˚A of < fλ >= 1.02 ± 0.05 ×10−20 erg s−1 cm−2˚A−1 . Over a λrest = 4˚A region redward of Lyα we measure a mean deficit against the expected continuum of −0.8 ×10−20 erg s−1 cm−2˚A−1 (missing 80% of the expected continuum flux) at a 5.7σ significance (based on the noise properties of the spectrum). Over a broader λrest = 24˚A region we measure a mean deficit of −0.4 × 10−20 erg s−1 cm−2˚A−1 at a 6.5σ significance (40% of the expect continuum flux). Our findings are at face value consistent with the expected transmission due to damping wing in a neutral IGM from [66], with a predicted flux of ∼ 0.3×continuum and ∼ 0.8×continuum at 1000 km s−1 (∼ 4˚A) and 6000 km s−1 (∼ 24˚A). From this, the flux deficit can be interpreted as Lyα absorption and thus as evidence that this galaxy does not sit within a large ionized bubble. Similar interpretations have been drawn for 4 galaxies at 10 < zspec < 14 which also show a softening of their Lyman break [17]. However, this physical interpretation is challenging to reconcile with the evolved stellar populations and large stellar mass estimates, which imply the reionization process of the IGM should be well underway. In this respect, we highlight that Lyα has been observed in emission in GN-z11, a younger, less massive galaxy at higher redshift [23]. An alternative interpretation of the Lyα damping is explored in the context of ISM and CGM absorption in 2.7.7. 2.7.7 UV absorption features We apply the same method used for Ly-α absorption to measure the significance of a deficit between the expected stellar continuum and the measured continuum for several common UV absorption lines which fall within our spectral coverage (SiIIλ1260, 1304, OIλ1302, CIIλ1335, FeIIλ2344, 2374, 2382). We examine the region of ±5000km s−1 centered on each line and measure the average continuum level outside the central ±500km s−1 . We measure the significance of any absorption within this central region against the continuum level and find greater than > 3σ detections in SiIIλ1260, CIIλ1335 and FeIIλ2344, with EWrest = (−3.6 ± 0.7)˚A, (−3.7 ± 0.8)˚A and (−6.2 ± 1.6)˚A, respectively. The remaining lines show < 2σ significance. To examine potential UV-absorption further we stack the ±5000km s−1 window for these 7 lines, Springer Nature 2021 LATEX template 24 A massive interacting galaxy 510 million years after the Big Bang 3000 2000 1000 0 1000 2000 3000 0 0.5 1 1.5 2 s t a c k e d F mean 1.02± 0.02 6.3 at -500:500 km s 1 3000 2000 1000 0 1000 2000 3000 -1 0 1 2 3 F 1 0 2 0 e r g s 1 c m 2 Å 1 SiII 1260 5.1 at -500:500 km s 1 mean=1.07±0.06 3000 2000 1000 0 1000 2000 3000 Velocity offset km s 1 -1 0 1 2 CII 1335 4.8 at -500:500 km s 1 mean=0.88±0.05 3000 2000 1000 0 1000 2000 3000 -0.2 0 0.2 0.4 0.6 FeII 2344 3.8 at -500:500 km s 1 mean=0.23±0.02 Extended Data Figure. 4 UV absorption features in the high-resolution spectroscopy of Gz9p3. Top panel: Stack of region ±3000km s−1 centered on common UVabsorption lines (SiIIλ1260, 1304, OIλ1302, CIIλ1335, FeIIλ2344, 2374, 2382). The orange filled region highlights the −500 : 500km s−1 window which shows a series of absorption features exhibiting a 40% reduction flux compared to the mean stellar continuum (red line) at a 6.3σ significance. Bottom panels: the individual SiIIλ1260, CIIλ1335 and FeIIλ2344 absorption features shown over the same velocity window. scaling the measured continuum level to unity and using a 1/σ2 weighting. In Extended Data Figure 4 we present the stacked spectrum plus the individual spectral regions for the detected absorption lines. The stacked spectrum shows a double absorption feature centered at ±100km s−1 . The deficit of the combine absorption feature over a [-500:500]km s−1 window is a 40% reduction flux compared to the stellar continuum at a 6.3σ significance (57% reduction at 7.5σ over a [−500 : 200]km s−1 window). The equivalent width of SiII and CII absorption lines is closely correlated with the Lyα equivalent width as these ions arise from HI gas, which is also responsible for the damping of Lyα [67–69, 104, 105]. From [69], we estimate that rest-frame equivalent width of ≲ −2 ˚A for low-ionization interstellar metal lines is associated with weakened Lyα emission and moderate absorption damping wings. Therefore, our measurement of SiII and CII absorption with EW ∼ (−3.6 ± 0.8)˚A can qualitatively explain the damped Lyα feature present in the spectrum redward of the Lyman break without necessarily requiring Gz9p3 to live within a partially neutral IGM region. While our data are inconclusive to determine the IGM ionization, we note that the combination of strong absorption for low-ionization ISM metal lines and weak CIII] emission (EWrest < 1.0˚A) suggest that the galaxy is unlikely to be a Lyman continuum leaker in its current state [70, 105]. In both the stack and the individual lines the overall shape of the deficit shows hints of a multi-component structure. The absorption is centered at a 100km s−1 offset to the expected location, either to one side or both. Such multi-component absorption may indicate the presence of turbulent gas within Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 25 the system, but the signal to noise of our data is not sufficient to draw robust conclusions on this potential interpretation. 2.8 Measurements of integrated physical parameters Given that the NIRSpec slit does not cover the full extend of the object due to its orientation (see Figure 1), we additionally use integrated-light photometry to derive an estimate for the stellar mass and the 100Myr time-averaged star formation rate (SFR) for the whole galaxy by fitting the spectral energy distribution using two alternative SED-fitting software tools; BAGPIPES [71] and ZPHOT [106]. The application of each tool for the GLASS photometric dataset is described in [107] for BAGPIPES and [108] for ZPHOT, and we refer the reader to these articles for a detailed discussion. Briefly, for BAGPIPES, we model the photometry using a lognormal SFH, a Kroupa [109] initial mass function (IMF), a Calzetti [84] dust attenuation law and a BPASS (v2.2.1 [110]) stellar population model that includes binary populations that is reprocessed with photoionization code CLOUDY (c17.03 [111]) generated in [34, private communication with Adam Carnall]. Consistent stellar mass and SFR estimates within 1σ are obtained using a delayed exponential SFH or when using the 2016 version of the BC03 [112] stellar population models. For ZPHOT, the important changes are the use of a delayed exponentially declining Star Formation Histories, the 2016 version of the BC03 [112] stellar population models and a Chabrier [113] IMF. In Section 2.7.2 we place constraints on the metallicity of Gz9p3 and as part of our SED fitting we apply a Z=[0.01-0.33]Z⊙ flat prior to match the measured range of metallicities determined from our Ne3O2 analysis and their associated systematic uncertainty. We note that we obtain consistent measurements (within 1σ) for our key parameters (Mass, SFR, Age) when we place no prior on the metallicity. We run both tools twice, first allowing the redshift to vary and then fixing the redshift at zspec = 9.313 (see 2.5). All variations of the photometric analysis return consistent estimates of the stellar mass and star formation rate, and we find the photometric redshift estimates for both tools are reasonably consistent with the spectroscopic value and place the redshift firmly at z > 9 (zphot: 10.2 and 9.7 ± 0.3 for ZPHOT and BAGPIPES respectively). By fixing the redshift at the solution derived from emission line analysis, BAGPIPES and ZPHOT estimate an observed log10(M∗/M⊙) of 9.43+0.11 −0.15 and 9.51+0.20 −0.14 and SFRs (M⊙/yr) of 31+8 −10 and 57+9 −21. The consistency between the tools is discussed in [108], who find good agreement in all but the lowest mass galaxies. We present the best-fit SED using the BAGPIPES fixed-redshift model in Extended Data Figure 5. We correct our measurements for the lensing magnification (µ = 1.658+0.016 −0.015; see 2.3) and we report that our target has an intrinsic log10(M∗/M⊙) = 9.2 +0.1 −0.2 (or 9.3 +0.2 −0.1 ) and SFR = 19+5 −6M⊙/yr (or Springer Nature 2021 LATEX template 26 A massive interacting galaxy 510 million years after the Big Bang 34+5 −13M⊙/yr) for BAGPIPES (or ZPHOT). These measurements are consistent with those obtained by [47] using ZPHOT, with [114] templates and a delayed-star forming history (SFR=40+70 −10M⊙/yr, log10(M∗/M⊙)=9.2 ± 0.3). 2.8.1 Mass metallicity relation This galaxy exhibits a large stellar mass and low oxygen abundance. In Figure 3 we present the location of Gz9p3 relative to measured mass-metallicity relations. The measured metallicity range for Gz9p3 [12 + log10(O/H) = 7.2−7.8] places it offset below local z = 0 [56, 115] and moderate redshift z = 2 − 4 mass-metallicity relations [56, 59] (which estimate 12+log10(O/H) in the range ∼ 8.5 − 8.8 and ∼ 7.7 − 8.3 at log10(M∗/M⊙)= 9.2, respectively). In a recent census of JWST galaxies at z ∼ 4−10, [58] report little evolution in the massmetallicity relation within their metallicity uncertainties (∼ 0.3dex) and we find Gz9p3 remains offset below the best-fit mass-metallicity relation at these higher redshifts. In fact, the relation from [58] provides 12+log10(O/H) = 8.01 at log10(M∗/M⊙)= 9.2. 2.9 Resolved stellar populations In Figure 1, Gz9p3 is clearly isolated from close neighbors and shows an elongated spatially-resolved structure with relatively high S/N in individual pixels. The direct images have been PSF-matched to the F444W image and while the PSF-matching means the information encoding within a pixel will be blended with its neighbors, our target is spatially large enough to provide regions far enough apart to be resolved (FWHM of F444W is ∼ 4.5 native pixels, or ∼ 0.14”). This offers the opportunity to carry out a resolved stellar population analysis. Here, we follow the method of [46] and perform a pixel-by-pixel fitting of the SED. By fitting the physical parameters for each pixel we can create a 2D map of the distribution of galaxy properties (Stellar Mass, SFR, Stellar Age). To ensure we have sufficient sensitivity to perform pixel-by-pixel SED model fitting, we first bin the direct image using a (2×2) kernel to improve the signal to noise (SNR) and then enforce a SNR > 2 threshold in the F150W and F200W NIRCam broadband filters (which have the highest noise levels). After this binning, there are ∼ 60 pixel that meet this threshold, which sample the full extent of the object with a (binned) pixel scale of 0.062”/px, which at z = 9.313 scales to 0.27kpc/px. For each pixel we fit the 11-band HST+JWST photometry using BAGPIPES (as discussed in 2.8) and we fix the redshift to zspec = 9.313. Our SNR > 2 threshold choice is sufficient to obtain robust fits, with the log of the Bayesian evidence for all pixels > 2. From BAGPIPES we obtain estimates for the star formation rate and stellar mass, which we convert to surface density properties using the 0.27 kpc/pixel scale, as well as an estimate for the mass-weighted stellar age. We Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 27 additionally measure the UV slope β (where fλ ∝ λ β ) using the F150WF200W colour, taken from the BAGPIPES best-fit SED magnitudes. We determine β = C × (m1 − m2) − 2 (e.g., [116]), where C −1 = 2.5 · log10( λ2 λ1 ). For F150W-F200W, the effective wavelength λef f = 15007˚A and 19886˚A trace the rest-frame UV for our galaxy’s redshift (λrest = 1455˚A−1930˚A), and set C = 3.27. Here we adopt the ‘pivot’ wavelength as the effective value to account for the transmission profile and detector sensitivity as a function of wavelength for each filter. In Figure 4 we present the 2D distribution of these physical properties, where for each pixel we display the 50th percentile of the posterior and in the supporting panels we show the uncertainty of the estimates based on the 16th and 84th percentiles. Within the 2D distributions, the presence of multiple components is apparent; the main body of the object, a tail and a compact region at the end of the tail (highlighted in the F150W-F444W panel of Figure 4 by the 10σ F150W contour). The SFR surface density peaks in the center of the main body, with minimal star formation found in the tail. The mass-weighted age for the main body and tail-clump is estimated to be < 50Myr, with the tail showing older ages albeit with a larger uncertainty. Across the whole 2D distribution of the galaxy the estimates for the visual extinction are below Av = 0.7, with a median of only Av = 0.2. It is clear that the SED modeling prefers solutions with little to no dust extinction. Like the SFRD, the stellar mass surface density peaks in the center of the main body with a SMD of ∼ 106M⊙/kpc−2 maintained throughout the tail. The β slope estimates range from −1.5 to −2.5 and trace the stellar age 2D distribution, with the steepest slopes (associated with the youngest stellar populations) correlating with the stellar age. 2.9.1 Properties of Integrated regions In addition to presenting the pixel-by-pixel 2D distributions of the physical parameters, we investigate the integrated properties of three distinct regions (with greater spatial separation than the FWHM). Extended Data Figure 5 presents the placement of three apertures which follow the 5 and 20σ contours of the F277W direct image. One aperture (r = 3px) is placed over clump in the tail, which traces the region of low stellar ages and blue β slope seen in Figure 4. Two apertures are placed over the main component of the source, a compact aperture (r = 3px) following the 20σ contour and a larger ellipse (6px/7px semi-minor/major axis) following the 5σ contour, the placement of these two creates an annulus. The inner aperture traces the compact peak of the SFRD and SMD distributions while the annulus covers the surrounding lower SFRD and SMD region. To compare different regions of the target, we use the photutils package in Python to measure the aperture photometry. We perform BAGPIPES SED modeling of the aperture photometry from each region following the method in 2.8, obtaining tighter constraints on the physical parameters than in the pixel-by-pixel fitting. The photometry and best-fitting SED models are shown in the right panel of Extended Data Figure 5. Springer Nature 2021 LATEX template 28 A massive interacting galaxy 510 million years after the Big Bang Property Main (phot-only) Inner Tail Annulus Stellar Mass [log10(M∗/M⊙)] 8.3 +0.3 −0.2 7.9 ± 0.1 7.3 +0.3 −0.2 8.2 +0.3 −0.2 Star formation rate [M⊙ yr−1 ] 2.2 +2.5 −0.7 0.9 +0.4 −0.3 0.2 +0.3 −0.1 1.5 +1.6 −0.6 Stellar Age [Myr] 7+15 −3 8 +7 −4 8 +20 −4 7 +20 −4 Av 0.1 ± 0.1 0.03+0.04 −0.02 0.3 +0.3 −0.2 0.3 ± 0.2 β Slope −2.3 ± 0.1 −2.48+0.05 −0.03 −2.2 +0.4 −0.2 −2.1 ± 0.2 Muv [AB mag] −21.1 ± 0.1 −20.34+0.04 −0.03 −18.2 +0.3 −0.2 −20.4 +0.2 −0.1 Extended Data Table. 2 Galaxy properties for regions within Gz9p3, fit using aperture photometry shown in Extended Data Figure 5, corrected for magnification when appropriate. The spectrum+photometry SED fitting results for the main component is shown in Table 1. The Inner region and the surrounding annulus have comparable F150W and F200W flux densities, however diverge at longer wavelengths, with the inner region showing a steeper UV-slope (β = −2.48+0.05 −0.03 and −2.1 ± 0.2, respectively) and a fainter rest-optical stellar continuum. The clump in the tail shows similar properties to the annulus region, with a β = −2.2 +0.4 −0.2 UV slope. All three regions exhibit consistent sSFR (∼ 10Gyr−1 ) and Stellar ages (∼ 8Myr, with 16th/84th percentiles of ∼ 3 and ∼ 20Myr), with the distinguishing property being the visual extinction which is significantly lower in the Inner region (Av = 0.03+0.04 −0.02) than in the annulus and tail (Av = 0.3 ± 0.2 and 0.3 +0.3 −0.2 , respectively). In each region the extinction is minimal and for the Inner region it appears negligible, possibly suggesting inflow of pristine gas to the core due to a merger event (see Section 2.10). The inferred properties of these regions from the photometry alone are presented in Extended Data Table 2 and all show young stellar ages (< 10Myr), whereas the fit to the spectrum and photometry of the main component (the combined inner+annulus region) in 2.7 estimated a higher age (130 ± 20Myr). These differences highlight how the light from the young stellar populations dominates the photometry in the rest-UV and masks the older population, which can be seen in the spectra based on the emerging Balmer break. Likewise, the stellar mass estimate for the main component using the spectrum+photometry (log10(M∗/M⊙)=9.15 ± 0.04) is higher due to the inclusion of an older stellar population. From the photometry-only SED fitting, the stellar mass of the main component accounts for about one fifth of the fullphotometry estimate, but for half of the MUV flux density. Correspondingly, the outskirts of the main component (outside the aperture) and the extended tail exhibit a high stellar mass-to-UV-light ratio. This follows the pattern of SMD, stellar age and SFRD properties in the 2D distributions of Figure 4. We also note that the very young stellar population ages (< 10Myr) preferred in the fits of the individual regions have lower stellar mass-to-UV-light ratios, with the fit outcome dominated by the rest-frame UV light. The discrepancy in the star formation rates and stellar masses between SED fitting from broadband photometry and SED fitting from the spectrum highlights the importance of spectroscopy to study galaxies during the epoch of reionization. Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 29 Extended Data Figure. 5 Spectral energy distribution for Gz9p3 from integrated light in key regions. Left Panel: Overlaid onto the F277W direct imaging for the galaxy, white dashed lines show the 5σ and 20σ F277W contours. Three apertures are placed to approximately trace these contours. These include an aperture (green) over the tail tracing the 5σ contour of a clump and an inner and an outer aperture over the main component tracing the 20σ and 5σ contours (red, blue). Right Panel: Photometry and BAGPIPES SED fit (16th − 84th percentile shown by shaded region) for the full galaxy (orange, see Extended Data Table 1) and key regions from the left panel (blue = Inner, red = the annulus region created between the inner and outer apertures on the main source, and green = Tail, which we scale by ×4 to aid the viewer). Error bars derive from uncertainty on the broadband imaging flux density and the SED fitting. 2.10 Morphological properties and evidence for a merging system The non-negligible mass in the tail of Gz9p3 suggests the extended structure is a result of a merger rather than being due to gas stripping, which would show very low masses in a tail [117]. This is further supported by the presence of a double nucleus (reminiscent of a merger remnant) in the F150W-F444W panel (and in Figure 1) and by star formation concentrated to the core, which may indicate a nuclear starburst triggered by gas infall during a merger [118– 120]. On the other hand, the pixel-by-pixel SED modeling returns a relatively smooth spatial distribution of properties, which could be interpreted as evidence that the system is observed a few dynamical times post-merger and has therefore largely completed its relaxation process. Spectroscopic coverage over the tail, which is unavailable due to the placement of the slit, would likely allow stellar and gas kinematic measurements to further distinguish between a merger history and any gas stripping from either ram pressure stripping within the z ∼ 10 over-density reported by [47] or due to cosmic web stripping [121]. Lacking such spectroscopic data, to examine further whether this system is a merger, we utilize a tailored data reduction of the NIRCam direct imaging at short wavelengths, PSF matching to F200W (the longest wavelength filter with the coarsest resolution we consider in this Section, instead of F444W). This allows greater spatial resolution in the short wavelength F150W Springer Nature 2021 LATEX template 30 A massive interacting galaxy 510 million years after the Big Bang and F200W filters, which are drizzled to a 20mas/px resolution. We weightcombine the F150W and F200W direct images to improve the signal to noise ratio before investigating the morphology. The combined image is presented in Figure 5 and visually shows two cores within the center of the galaxy. We fit both core components (Left and Right as viewed) with Sersic profiles and measure consistent profiles with Remajor = 0.032”, 0.033”(corresponding to 0.14kpc), nsersic = 0.44, 0.55 and axis ratios of 0.81, 0.76. The central region (excluding the tail) shows a greater axis ratio of 0.44, with a Remajor = 0.15” (0.65kpc) for a nsersic = 0.50 profile. The best-fit F150W+F200W profile, which traces the rest-UV and hence active star forming regions, shows a smaller size and shallower gradient than the F444W best-fit profile to the central region (Remajor = 0.26”, nsersic = 3.54, axis ratio= 0.26) which traces the rest-optical and the older stellar population. 2.10.1 Morphological indexes Several quantifiable indicators exist for describing the morphology of a galaxy, with thresholds to indicate whether a system is likely to have gone, or be going through, a merger: Gini - M20 and asymmetry (A). These morphological parameters are defined in [72, 122] and [123], respectively. [72] define mergers as galaxies with Gini > 0.33 −0.14× M20, while [124] adopt the condition A ≥ 0.35 as a merger criterion. These parameters were already adopted to study the morphology of galaxies at the epoch reionization with JWST observations [125]. We thus follow the same measurement procedure on the combined F150W+F200W galaxy cutout shown in Figure 5. We determine each morphological parameter following the literature definitions with our inhouse JWSTmorph package. We obtain Gini = 0.61, M20 = −1.29, and A = 0.35, which robustly classify our system as a merger. 2.10.2 Clumpiness The system appears also very clumpy in Figure 5, with two bright and distinct clumps clearly detected in the central region. These might represent the two merger components. To be more quantitative, we calculate the clumpiness parameter c, which is defined as the fraction of light of the galaxy residing in clumps. This traces the relative importance of small scale structures inside a galaxy, and can shed light on star-formation episodes outside of the nucleus or material (gas+stars) that has been stripped away from the galaxy by tidal phenomena. At intermediate redshift, the presence of bright clumps is tightly associated to merger events [126]. To compute the clumpiness, we follow a similar approach to [126]. In brief, we smooth the original image (F150W + F200W combined) using a gaussian filter with a size of 0.2 ′′, corresponding to ∼ 1 kpc at the target redshift. We then subtract the smoothed image from the original image, and select the pixels with a flux at least 2σ above the background, in order to reduce the noise contamination. We also impose 0 for all the pixels with a negative value in the residual [124]. Finally, the clumpiness is calculated as the ratio between Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 31 the flux in those pixels and the total flux of the galaxy. We obtain a final value of c = 0.56, and the clump-map is shown in the right panel of Figure 5. In the clump-map we can see that, while most of the flux is coming from the two central nuclei, there are 4 clumps detected in the left tail of the system. If we remove the nuclei, we still find that 6.3% of the total flux (which roughly traces the SFR as we probe the rest-frame UV) resides in these external structures, which is completely different from the compact, nucleated morphology of typical star-forming galaxies seen at these redshifts [125]. We can also notice the elongated shape of the main clump in the center, which contributes to increase the asymmetry of the system. All this evidence suggests that a merger is ongoing and that the elongated, clumpy structure might be the result of the tidal forces generated during the interaction, as observed in lower redshift mergers. 2.11 Comparison to theoretical modeling The properties of Gz9p3 can be compared to basic theoretical modeling of galaxy formation in ΛCDM as a first step to assess the ability of current model to reproduce the observations. To estimate the number density and relation between stellar mass and MUV , we consider a set of models where DM halos are populated with stellar populations based on simple analytical recipes [28, 60, 127]. In summary, the modeling relies on three elements (1) the stellar-to-halo mass ratio is mass dependent but redshift independent; (2) the age of the stellar population is proportional to the DM halo assembly time and its scatter derives from the probability distribution of the halo assembly time based on extended Press Schechter theory; (3) the star formation efficiency is calibrated via abundance matching at a single reference redshift where a robust determination of the UV luminosity function is available. Using the results from [28] we can see from their Fig. 3 that Gz9p3 lies very close to the modeling predictions for the stellar mass versus MUV relationship at the median age of ∼ 80 Myr, indicating it should be a typical object for its mass and luminosity (and dust content). Based on that model, the inferred DM halo mass of the system is MDM ∼ 1011.4 M⊙ and the inferred number density is ∼ 10−6 M pc−3 (comoving). It is challenging to estimate precisely the parent survey volume to associate to this object (given it was included as a NIRSpec target based on pre-existing HST photometry), however bounding within a factor of two is possible. The lower limit is to consider the NIRSpec FoV of ∼ 3 × 3 arcmin2 and ∆z ∼ 1, while the upper limit is the area of the HST photometry that includes WFC3/IR data (needed given the redshift of the source), which is approximately 13 arcmin2 [125]. This gives a selection a volume of ∼ (2 − 3) × 104 M pc−3 (comoving) (see e.g. the cosmic variance calculator from [128]), and it implies only a few percent probability of having such a massive and bright galaxy in our observations. Given that the relationships between stellar mass, UV luminosity, and dust are in agreement with theory, the natural degree of freedom to change is an increased star formation efficiency relative to the halo mass as the redshift increases. This has been Springer Nature 2021 LATEX template 32 A massive interacting galaxy 510 million years after the Big Bang proposed in similar modeling, e.g. by [129], and an increased star formation efficiency could also explain the higher than expected density of photometric candidates at z ≳ 10 from NIRCam observations (e.g. [47]). Another possibility is that the source is hosted by a dark matter halo mass that is less massive than the average halo to stellar mass relationship implies, and thus we are observing an outlier with higher than average star formation efficiency and stellar mass. Based on previous work, it is possible that the efficiency in a single object is up to one dex above mean (i.e. a factor 10×) (e.g. [28, 127]). Theoretical modeling is also useful to estimate the probability of observing a disturbed morphology deriving from a galaxy merger. The merger rate of DM halos per unit redshift is nearly universal across primary halo mass and redshift and only depends on the mass ratio [3]. Estimating a mass ratio between 1/10 and 1/5 from the photometry of Gz9p3, [3] gives dNmerge/dz ∼ 1, which corresponds to one merger expected every ∼ 70 Myr at z ∼ 9.3. Assuming further that the disturbed morphology remains observable over a timescale of a few dynamical times for the system. The latter can be estimated to be tdyn ∼ 5 Myr considering a characteristic velocity dispersion of ∼ 200 km/s (from virial equilibrium) and a characteristic radius of 3 kpc (from imaging). Therefore, we conclude that catching galaxies in the act of merging at this epoch is relatively likely (at the level of ∼ 10 − 20% probability). 2.12 Comparison to cosmological hydrodynamical simulations To further assess whether this galaxy is in tension with predictions from the ΛCDM model, we searched for analogous galaxies with masses 109 ≤ M∗/M⊙ ≤ 1010 in the IllustrisTNG simulation suite [130–134]. We found 13 analogous galaxies in Snapshot 5 of TNG300-1 at a redshift of z = 9.39, with a box length of 302.6 cMpc, and an average mass resolution of 1.1 × 107M⊙ for baryonic particles. No analogous galaxies were found in TNG100, the second largest simulation of the suite with a box length of 110.7 cMpc, indicating that such massive galaxies are rare objects in the early Universe. While 13 simulated galaxies does not provide a statistical sample, we do find that these simulated sources do share similarities to Gz9p3. Each of the 13 galaxies shares a similar SFR, with the mean SFR = 30 ± 9M⊙yr−1 .

The merger tree history of the 13 systems indicates that two had undergone a major merger (Mass ratio > 0.2, 15%) and four had undergone a minor merger (Mass ratio < 0.2, 31%) within the cosmic time interval between z = 10 and z = 9.3 (Snapshots 4 and 5 of the simulation). Such a rate of major mergers falls in line with our prediction from theoretical modeling. Unfortunately, the spatial and mass resolution of the simulation is not sufficient to obtain a detailed morphology for a quantitative comparison. However, qualitatively, two simulated systems show a double-core structure (two distinct mass peaks) which reflects the structure of Gz9p3. Therefore, while current state of the art modeling and simulations suggest that the probability of finding an object such Springer Nature 2021 LATEX template A massive interacting galaxy 510 million years after the Big Bang 33 as Gz9p3 in the cosmological volume probed by this survey is low (approximately 6% considering 13 objects in the TNG300 and assuming a survey area of GLASS-ERS+DDT+UNCOVER of ∼ 60 arcmin2 ), setting aside number density, galaxies with properties similar to Gz9p3 are expected to exist in the early Universe. This may suggest that star formation efficiency during the epoch of reionization may need to be revised upwards to improve the data-model comparison.

Declarations • Data availability: All data used in this paper are publicly available through the Mikulski Archive for Space Telescopes (MAST) server with the relevant program IDs (ERS-1324 for the NIRSpec spectroscopy and DDT-2756 for the NIRCam imaging). The reduced NIRCam imaging utilised in this work from the GLASS collaboration [50] is available at https://doi.org/10.17909/ kw3c-n857. All other data generated throughout the analysis are available from the corresponding author on request. • Code availability: Our analysis makes use of several publicly available codes. The NIRSpec data were reduced using the msaexp code, which can be found here: https://github.com/gbrammer/msaexp . The data reduction of the NIRCam images were performed with the official STScI JWST pipeline, which can be found here: https://github.com/spacetelescope/jwst . The SED fitting analyses were performed with BAGPIPES, the latest version of which (including the templates used here) is available at https: //bagpipes.readthedocs.io/en/latest/. We modelled the observed spectral emission lines using the specutils packages within Python, which can be found at https://specutils.readthedocs.io/en/stable/ . We performed aperture photometry on the direct imaging using the photutils packages within Python, which can be found at https://photutils.readthedocs.io/en/stable . Galaxy morphological parameter were measured using the GLASS inhouse JWSTmorph package, which is publicly available at https://github. com/Anthony96/JWSTmorph.git. All other code generated throughout the analysis are available from the corresponding author on request. • Acknowledgements: This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program JWST-ERS-1324 and JWST-DDT-2756. We acknowledge financial support from NASA through grant JWST-ERS-1324. KB, MT, BM, ND acknowledge support from the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. KG and TN acknowledge support from Australian Research Council Laureate Springer Nature 2021 LATEX template 34 A massive interacting galaxy 510 million years after the Big Bang Fellowship FL180100060. BM acknowledges support from Australian Government Research Training Program (RTP) Scholarships and the Jean E Laby Foundation. We acknowledge financial support through grants PRINMIUR 2017WSCC32 and 2020SKSTHZ. MB acknowledges support from the ERC Grant FIRSTLIGHT and from Slovenian national research agency ARRS through grants N1-0238 and P1-0188. CM acknowledges support by the VILLUM FONDEN under grant 37459 and the Carlsberg Foundation under grant CF22-1322. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. We acknowledge support from the INAF Large Grant 2022 “Extragalactic Surveys with JWST” (PI Pentericci). EV acknowledges support from the INAF GO Grant 2022 ”The revolution is around the corner: JWST will probe globular cluster precursors and Population III stellar clusters at cosmic dawn”. MC acknowledges support from INAF Minigrant “Reionization and fundamental cosmology with high-redshift galaxies”. PS acknowledges INAF Mini Grant 2022 “The evolution of passive galaxies through cosmic time”. DM acknowledges financial support from program HST-GO-17231, provided through a grant from the STScI under NASA contract NAS5-26555. • Authors’ contributions: KB identified the emission lines from the NIRSpec data, led the overall data analysis activities, produced all figures, and was primarily responsible for writing the Methods section. MT provided advice on the data analysis and on its physical interpretation, carried out the comparison to theoretical modeling and contributed associated text in Methods, and was primarily responsible for writing the Abstract and Main sections. NL led the SED fitting and contributed associated text in Methods. AC led the clumping analysis and contributed associated text in Methods. BM led the comparison to hydrodynamical simulations and contributed associated text in Methods. GRB led the NIRSpec data reduction. ND led the Lyman Break modeling. LY led the light profile fitting from imaging data. TT led the GLASS/ERS survey conception, design and execution as the Principal Investigator of the program, and contributed advice on the paper preparation. TJ and AH contributed to physical interpretation of the absorption lines. AH, CM, TM, TN, and XW contributed to the NIRSpec data reduction and to the development of the NIRSpec pipeline. AF, EM, CM, DP contributed to the NIRSpec data reduction and to the development of the NIRCam pipeline. All authors contributed comments during the research activities and the manuscript preparation.