© The Authors 2025

1. Introduction

The James Webb Space Telescope (JWST, Gardner et al. 2023) is transforming our understanding of galaxy formation, particularly during the so-called Epoch of Reionization (EoR). The JWST has photometrically detected potential galaxy candidates up to redshift 25 (Pérez-González et al. 2025) and spectroscopically confirmed galaxies up to redshift 14.44 (Naidu et al. 2025). At these redshifts, the rest-frame optical emission lines and continuum are redshifted into the mid-infrared (5–25 μm), which is the spectral range covered by the Mid-Infrared Instrument (MIRI, Rieke et al. 2015; Wright et al. 2015, 2023). At redshifts of 7 and above, this instrument allows us to explore the strong optical emission lines, such as Hα.

Mergers are expected to have a greater impact on galaxy evolution at higher redshifts than at lower redshifts, as the merger rate increases with redshift (Duan et al. 2025). The gravitational forces from mergers can cause significant disruptions in galaxy structures and movements, leading to features such as tidal tails and asymmetric stellar envelopes (e.g., Toomre & Toomre 1972; Mihos & Hernquist 1996; Bushouse et al. 2002; Conselice et al. 2003; Spilker et al. 2022). These disturbances can cause gas to flow from the outer regions into the center of galaxies, triggering the formation of new stars in central areas (Pérez-González et al. 2024a). Simulations such as FirstLight (Ceverino et al. 2017) show that about 10% of galaxies at redshifts of 6–10 exhibit clumpy structures, which are primarily formed through mergers. The central clumps are mainly composed of older stellar populations exceeding 50 Myr, while off-center clumps originate from gas debris along tidal tails during major mergers. These mergers have a short life, typically merging into the main galaxy within a few tens of million years (Nakazato et al. 2024). Observational studies also support the presence of mergers in the EoR, with recent works identifying clear merger signatures in high-redshift galaxies (e.g., Álvarez-Márquez et al. 2023; Arribas et al. 2024; Harshan et al. 2024; Hu et al. 2025).

High-z star-forming galaxies (SFGs) are characterized by optical lines with high (rest-frame) values of the equivalent widths (EW up to 3000 Å) in extreme emission line galaxies (Boyett et al. 2024). There is also a trend of EW increasing with redshift, which indicates an evolution in the galaxy-specific star formation rate (Rinaldi et al. 2023, 2025) and young stellar populations (< 10 Myr). Studies at low redshifts have shown that galaxies in close pairs exhibit an SFR increase compared to isolated galaxies, highlighting the impact of gravitational interactions on gas dynamics and star formation (Barton et al. 2000). This enhancement is also observed at higher redshifts (6.5–8.5), with an excess in SFR of 0.26 ± 0.11 dex found in very close pairs, with no significant increase at larger separations (Duan et al. 2024; Puskás et al. 2025).

B14-65666 (also known as the Big Three Dragons galaxy) was discovered based on wide-field imaging data of the UltraVISTA survey of the Cosmological Evolution Survey (COSMOS) field in Y, J, H, and K_s near-infrared bands (Bowler et al. 2012, 2014). It was then spectroscopically detected in Lyα (Furusawa et al. 2016). An image captured with the Hubble Space Telescope (HST) showed that the galaxy consists of two clearly separated components in the rest-frame ultraviolet (UV), thus suggesting that the object could be undergoing a merger of two galaxies (Bowler et al. 2017). Observations performed with the Atacama Large Millimeter/submillimeter Array (ALMA) with Band 6 reported the detection of an ≈160 μm dust continuum at the position of the components (Bowler et al. 2018). It also detected [O III] 88 μm, [C II] 158 μm, and underlying dust continuum emission at ≈90 μm in Bands 6 and 8 (Hashimoto et al. 2019). Furthermore, they showed that [O III] 88 μm and [C II] 158 μm can be decomposed into two components that are spatially resolved and kinematically separated. These components are associated with the two UV emission peaks identified in the HST images. Additional observations with ALMA Band 7 did not detect the presence of the far-infrared [N II] 122 μm fine-structure emission line, nor did ALMA Band 3 observations detect other lines such as CO (6-5), CO (7-6), and [C I] (2-1) (Sugahara et al. 2021; Hashimoto et al. 2023a).

Recent JWST Near Infrared Camera (NIRCam) imaging has revealed the complex morphology of the two galaxy components: galaxy E consists of a compact core, surrounded by diffuse, extended, rest-frame optical emission, and galaxy W has an elongated morphology interpreted as tidal tails (Sugahara et al. 2025). Follow-up observations with NIRSpec IFU from 2.9–5.3 μm directly map the strength of the rest-frame optical emission lines across the source (Jones et al. 2024). The profile of these lines revealed the presence of both narrow and broad line emission, with the latter indicating tidal interaction or outflows. This analysis showed that both galaxies have subsolar metallicities (0.24 Z_⊙ and 0.3 Z_⊙, E and W, respectively) and a high rate of star formation. The observations do not show evidence for the presence of an active galactic nucleus (AGN) in B14-65666.

B14-65666 is one of the few galaxies identified as undergoing a major merger during the EoR. This system allows the investigation of the physical properties of an early major merger, including its stellar structure, as well as the different phases of the interstellar medium (ISM). Additionally, this system provides a reference for direct comparisons with simulations. In this paper, we present the first spatially resolved spectroscopy of the Hα emission in the system B14-65666 at z = 7.15 based on new MIRI/JWST integral field spectroscopy. The paper presents several of the physical properties of B14-65666 based on the combination of the new MIRI data with ancillary NIRCam imaging and ALMA spectroscopy. This paper is organized as follows. In Sect. 2, we present the observations and data calibrations. In Sect. 3, we outline our analysis methods of spectroscopy and photometry. In Sect. 4, we present the results of the analysis. In the subsections, we study the spectral energy distribution of the galaxy (Sect. 4.1), the attenuation (Sect. 4.2), the star formation rate and burstiness (Sect. 4.3), the Hα EW and ionizing photon production efficiency (Sect. 4.4), the kinematics (Sect. 4.5), the stellar and dynamical mass (Sect. 4.6), the mass-size relation (Sect. 4.7), the starburst nature (Sect. 4.8), and finally the origin and evolution of the system according to simulations (Sect. 4.9). We discuss and analyze the spectral and spatial distribution of the Hα line and the derived properties of the galaxies, and we study the morphology and SFR of the galaxy and compare it with local galaxies. Our conclusions are presented in Sect. 5.

In this paper, we assume a flat ΛCDM cosmology with Ω_m = 0.310, and H₀ = 67.7 km s⁻¹ Mpc⁻¹ (Planck Collaboration VI 2020). With these parameters, 1″ is equal to 5.275 kpc at z = 7.15 and the luminosity distance is D_L= 72276.9 Mpc. For these cosmological parameters and redshift, the age of the Universe corresponds to 740 Myr. We use vacuum emission-line wavelengths throughout the paper.

2. Observations and data calibrations

2.1. MIRI

B14-65666 was observed with MIRI on 30 April, 2023 as a part of the MIRI European Consortium guaranteed time observations (program ID 1284, PI: L. Colina). The observations are composed of MIRI imaging (MIRIM, Bouchet et al. 2015; Dicken et al. 2024) and integral field spectroscopy using the Medium Resolution Spectrograph (MRS; Wells et al. 2015; Argyriou et al. 2023).

2.1.1. MIRI spectroscopy

The MRS observations were performed with the SHORT band for channels 1–4. We only used Channel 1 in this study. This channel covers the wavelengths 4.9–5.74 μm, thus targeting the Hα emission line at redshift 7.15 and other fainter emission lines such as [N II] 6550,6585 Å and [S II] 6718,6733 Å. The total on-source exposure time of the MRS observations is 22 743 seconds. The positive and negative four-point dither patterns for point sources were used. For each dither position, a total of six integrations were obtained using the SLOWR1 readout mode with 19 groups each. When the MRS targeted B14-65666, a parallel field was imaged using the MIRI F770W filter. This imaging data is utilized exclusively for the astrometric calibration of the MRS (see Sect. 2.4).

We used version 1.15.1 of the JWST calibration pipeline and context 1276 of the Calibration Reference Data System (CRDS) to process and calibrate the MRS observations of B14-65666. The standard MRS pipeline procedure (Bushouse et al. 2024) was used with additional customized steps to improve the quality of the final MRS-calibrated products and background (see Álvarez-Márquez et al. 2023, for details). Additionally, we implemented the algorithm developed by Spilker et al. (2023) to correct the striping pattern produced by the cosmic-ray showers in the spectral cube. The final Channel 1 cube has a spatial and spectral sampling of 0.13″ × 0.13″ × 0.8 nm (Law et al. 2023), and the resolving power (R  ≡  λ/Δλ) is about 3500 (Labiano et al. 2021; Jones et al. 2023).

2.1.2. MIRI imaging

The MIRIM observations of B14-65666 were performed using the F560W filter covering the rest-frame optical continuum at ≈0.7 μm and the Hα emission line. The MIRIM F560W image has a total on-source integration time of 2675 seconds. The observational setup consists of a four-point dither pattern, with the dither composed of two integrations of 120 groups using the FASTR1 readout mode.

The MIRIM F560W observations are calibrated with version 1.12.0 of the JWST pipeline (Bushouse et al. 2024) and context 1150 of the CRDS. We followed the general procedure to calibrate the data, but we also included more steps to correct for striping artifacts and background gradients (see Pérez-González et al. 2024b; Östlin et al. 2025, for details). The final MIRIM F560W mosaic has a pixel size of 0.06″, after resampling, and a field of view (FoV) of 74″ × 113″. The point spread function (PSF) full width at half maximum (FWHM) is 0.207″ (Libralato et al. 2024).

2.2. NIRCam

The NIRCam (Rieke et al. 2023) data used in this work were presented in a previous study (Sugahara et al. 2025). The observations were taken with filters F115W, F150W, F200W, F277W, F356W, and F444W (JWST GO1 program ID 1840; PIs: J. Álvarez-Márquez and T. Hashimoto). At redshift 7.15, F115W, F150W, and F200W sample the rest-frame UV spectrum, F277W and F356W cover the Balmer break region, and F444W targets the rest-frame continuum around 5000 Å and the [O III] 4960,5008 Å emission lines.

NIRCam and MIRIM images present large variations on the spatial resolution, with PSF FWHMs increasing with wavelength from 0.04″ to 0.21″. To obtain consistent multiwavelength aperture photometry in the different NIRCam and MIRIM images, we performed a PSF homogenization in the NIRCam images to match the spatial resolution of NIRCam images to the coarser MIRIM F560W image. However, we also used the NIRCam images with the native spatial resolution to investigate the structure of B14-65666 in the rest-frame UV (see Sect. 3.4).

The PSF homogenization procedure includes the following steps. We employed the observed effective PSFs (ePSFs) for MIRI (Libralato et al. 2024) and for NIRCam (Nardiello et al. 2022)¹. We compared these ePSFs with point-like sources from our observations and conclude that they are consistent with each other. The first step was to rotate the NIRCam ePSF to the position angle (PA) of the MIRI image (PA = 108°), ensuring orientation alignment between the two ePSFs. We rotate the MIRI ePSF to match the PA of the MIRI image. Next, we resampled the pixel scale of the MIRI ePSF (0.028″, Libralato et al. 2024) to match that of the MIRI imager and the NIRCam ePSF (0.008″ for F115W, F150W, and F200W; and 0.016″ for F277W, F356W, and F444W; Nardiello et al. 2022) to match the pixel scale of the MIRI imager. Using these aligned and resampled ePSFs, we then generated a matching kernel using the python package PHOTUTILS to bridge the differences between the NIRCam and MIRI ePSFs. Finally, to apply the kernel using the python function creatematchingkernel (Bradley et al. 2022), we rescaled it to the NIRCam image pixel size (0.03″ after resampling). Subsequently, we convolved the NIRCam images with the kernel, thus homogenizing the PSF.

2.3. ALMA ancillary data

The ALMA data used in this work were previously presented in other studies. ALMA Band-6 and Band-8 data were taken in Cycles 4 (ID 2016.1.00954.S, PI: A. K. Inoue) and 5 (ID 2017.1.00190.S, PI: A. K. Inoue) to target the [C II] 158 μm and [O III] 88 μm emission lines (Hashimoto et al. 2019). ALMA Band-7 data were taken during Cycle 7 (ID 2019.1.01491.S, PI: A. K. Inoue; Sugahara et al. 2021). In this article, we used the spectra of [C II] 158 μm and [O III] 88 μm and the emission-line and continuum (90, 120, and 160 μm) fluxes presented in Hashimoto et al. (2019) and Sugahara et al. (2021).

2.4. MIRI-NIRCam astrometry

We corrected the astrometry in the MIRI F560W image by using the centroid of two point-like sources within the FoV of the available GAIA DR3 (Gaia Collaboration 2023) coordinates. The uncertainty in the final positioning is 18 mas. With the corrected MIRIM F560W image, we aligned the simultaneous MIRIM F770W image with two point-like sources in the common FoV of both images. As the MRS data and the MIRIM F770W image were taken simultaneously, we used the astrometric offset found in the MIRIM F770W image to correct the MRS 1SHORT cube, obtaining an uncertainty of less than 60 mas (0.3 kpc). Finally, we aligned the NIRCam images using the astrometrically corrected MIRIM F560W image as we did for the simultaneous MIRIM F770W image. The NIRCam images have an uncertainty smaller than 35 mas; this was obtained with the available star and galaxy that have Gaia astrometry (Gaia Collaboration 2023).

3. Analysis

3.1. MIRI and NIRCam images and Hα line map

Previous HST images showed that B14-65666 is composed of two spatially separated galaxies in the rest-frame UV, a brighter galaxy in the northeast (galaxy E), and a fainter galaxy in the southwest (galaxy W) (Bowler et al. 2017). ALMA data further demonstrated that the [O III] 88 μm and [C II] 158 μm emission lines can be spatially and spectrally decomposed into two components associated with galaxies E and W (Hashimoto et al. 2018). NIRCam images revealed the complex morphology of the individual galaxy components (Sugahara et al. 2025): galaxy E has a compact core surrounded by diffuse extended optical emission, while galaxy W displays elongated clumpy morphology.

The MIRIM and MRS data cover the emission line of Hα together with the continuum at 0.7 μm. MIRI observations show a similar trend to NIRCam images. MRS shows an extended Hα emission associated with E and W galaxy components; additionally, each of these galaxy components is spectroscopically resolved as they were already presented in ALMA observations (Hashimoto et al. 2019). We present MIRI and NIRCam images in Fig. 1.

Fig. 1. Top left: MIRI Image F560W. The white circumference represents a 0.9″ radius aperture. The solid black elliptical line represents the E+W galaxies’ apertures for the Hα emission of the whole galaxy; the dashed black lines represent the apertures for galaxies E and W. The white filled circle area in the bottom left represents the spatial resolution (PSF FWHM) of the MIRI image F560W. Top right: Hα line map. The Hα line map is generated by integrating Hα line emission in the velocity range: –500 < v [km s−1] < 300. The white filled area represents the spatial resolution (PSF FWHM) of the MRS. Bottom: NIRCam cutouts in six filters, F115W to F444W (left to right and top to bottom). The images are homogenized to the PSF of the MIRI image, represented with the white filled circle in the NIRCam F115W image.

For this study, we defined a total elliptical aperture for the whole system (named galaxies E+W) with a semi-major axis of 0.456″, semi-minor axis of 0.343″, and position angle (PA) of 85.83°. In the case of the galaxy E we define an elliptical aperture with semi-major axis of 0.211″ and semi-minor axis of 0.201″, position angle (PA) of 112.83°. For the galaxy W, we extracted an ellipse with a 0.228″ and 0.177″ semi-major and semi-minor axis, respectively, and a PA of 112.83°. The coordinates of the center are shown in Table 1. The E and W elliptical apertures are compatible with the main emission of these galaxy components in our MIRIM F560W image. The elliptical aperture of galaxies E+W covers the extended Hα emission from the MRS Hα line map. The Hα line map is generated by integrating the Hα line emission, centered at the Hα peak of galaxies E+W (5.3510 μm) in the velocity range –500 < v [km s⁻¹] < 300 (see Fig. 1).

3.2. Hα spectra and line fluxes

We extracted the integrated Hα spectra of galaxies E+W, E, and W using the three defined apertures (Fig. 1). We corrected the two individual spectra for aperture losses, assuming that the two galaxies are unresolved sources for the MRS angular resolution taking into account the MRS PSF (Argyriou et al. 2023), which has a FWHM of 0.27″ × 0.35″ at 5.35 μm. The aperture correction is 1/0.57 for the galaxies E and W. In the case of the galaxies E+W, we also corrected for aperture losses, assuming that the flux is mainly dominated by galaxies E and W, which are unresolved. The aperture correction in this case is 1/0.81. We also corrected the spectrum for the contamination of flux from galaxy E in the W aperture by modeling galaxy E at the center of the aforementioned PSF and measuring how much of its flux would appear in the W aperture. We applied the same procedure for galaxy W to estimate its contribution to the E aperture. These corrections are 1/0.02 and 1/0.03, respectively.

We extracted twelve 1D spectra using the same apertures as the E+W, E, and W galaxies, respectively, at random positions of the MRS FoV clean of emission. We combined these spectra to generate the 1D median of the local background. The median is subtracted from the Hα spectra with the aim of removing any systematic residual feature left in the MRS calibration process. The three spectra are presented in Fig. 2.

Fig. 2. MIRI MRS 1 SHORT spectra for Hα observed at z = 7.1513. Top left: Hα spectra shown in velocity space for two spatially separated galaxies, E and W (red line and green line, respectively), identified in the Hα line map. The black line shows the integrated MRS spectrum of B14-65666 extracted using the galaxies E+W aperture. Top right: E+W galaxies’ integrated Hα spectrum. A vertical dashed line shows the wavelength of the weaker undetected lines [N II] 6550 Å and [N II] 6585 Å that are well separated in wavelength from Hα 6563 AA. Bottom left: Galaxy E’s integrated Hα spectrum. Bottom right: Galaxy W’s integrated Hα spectrum. The red dashed line shows the single Gaussian fits, with the residuals shown in the bottom panels. 1σ noise level is shown in yellow; it was calculated from the continuum adjacent to the emission line.

The Hα spectrum is modeled by a single Gaussian function, to fit the emission line, and a second-order polynomial, to fit any residual background gradient. Using a Monte Carlo simulation, we estimated the uncertainties for the derived emission line parameters (such as FWHM, flux, and central wavelength). To determine the noise level in the spectrum, the root mean square (rms) was calculated from the continuum adjacent to the emission line. Using this rms value as the 1σ standard deviation, we generated 1000 synthetic spectra by adding random Gaussian noise to the original spectrum. Each synthetic spectrum was then refit to obtain new measurements of the line parameters. The final uncertainty for each parameter is given by the standard deviation of these measurements across all simulations. The Gaussian fit is shown in Fig. 2, and the results are presented in Table 1.

From the Hα line, the E+W galaxies of B14-65666 are at a redshift of 7.1513 ± 0.0007. Galaxy E is redshifted to 7.1529 ± 0.0004; however, galaxy W is blueshifted to 7.1481 ± 0.0007. This implies that galaxy E and W present a velocity offset between the two components of 175 ± 28 km s⁻¹, compatible with that observed from [O III] 88 μm and [C II] 158 μm, 177 ± 16 km s⁻¹ (Hashimoto et al. 2019).

The Hα line is spectrally resolved with the MIRI-MRS. It presents an intrinsic FWHM of 407 ± 86 km s⁻¹ for galaxies E+W and 243 ± 41 km s⁻¹ and 312 ± 44 km s⁻¹, respectively, for E and W, after deconvolving with the line spread function of the MRS (Labiano et al. 2021). Due to the offset in velocities the component of the E+W galaxies presents a higher FWHM than that of galaxies E and W. The observed Hα fluxes are for E+W, E, and W: (23.1 ± 3.1) × 10⁻¹⁸, (11.9 ± 1.3) × 10⁻¹⁸, and (8.6 ± 1.2) × 10⁻¹⁸ erg s⁻¹ cm⁻², respectively. All results are presented in Table 1. The Hα flux for galaxies E+W is in agreement within the uncertainties of the sum of the individual galaxy fluxes. This suggests that the Hα emission is mainly dominated by the two galaxies excluding any relevant diffuse component.

3.3. MIRIM and NIRCam photometry

MIRIM F560W photometry for galaxies E+W and for galaxies E and W is derived with the same apertures as for the Hα spectra (see Fig. 1 and Sect. 3.1). The local background level for each aperture is estimated in an elliptical annulus centered at the same position as the E+W galaxies’ elliptical aperture and with an inner radius of 0.95″ and outer radius of 1.95″. Additionally, we used 300 apertures on blank regions of the MIRI image with the same size as the galaxy apertures to estimate the flux uncertainties. We also performed the NIRCam photometry using the same aperture as for MIRIM F560W. To do so, we utilized the NIRCam images (F115W, F150W, F200W, F277W, F356W, F444W) that are homogenized to the MIRI F560W PSF. The final observed MIRIM F560W and NIRCam fluxes are presented in Table 1. These NIRCam fluxes are lower by a factor of ≈0.13 than the ones obtained from previous results (Sugahara et al. 2025) because of the PSF homogenization from NIRCam F444W to MIRI F560W (see Sect. 2.2) and the different aperture sizes.

Galaxy E is an unresolved source in the NIRCam F150W image with an estimated physical size of 63 pc (see Sect. 3.4 for details). The selected aperture mainly contains emission from galaxy E; thus, we implemented aperture correction in the NIRCam and MIRIM photometry using the MIRIM F560W PSF. This aperture correction has a value of 1/0.67. On the contrary, galaxy W is fully resolved in the NIRCam F150W, and for that reason we did not implement aperture correction to its photometry.

3.4. Stellar structure and effective radius

We study the stellar structure of the two galaxies, E and W, individually. To achieve this, we utilized the F150W image, as it offers the best combination of signal-to-noise ratio (S/N) and spatial resolution. It presents the largest exposure time, 1675 s, arriving at a sensitivity of 7.7 nJy at 5 σ depth, which is a factor of 2–3 better than any other NIRCam band. Also, the PSF FWHM is 0.067″, allowing us to resolve structure up to 0.35 kpc at redshift 7.15 (Sugahara et al. 2025). The F150W image traces the UV emission, which is dominated by the young and attenuated stellar population. We worked with native-resolution images (before PSF-matching and resampling), and the analysis was done using the GALFIT code (Peng et al. 2002).

Galaxy E consists of a compact core, which appears unresolved, and a nondominant extended component. To obtain the effective radius of galaxy E, we modeled the F150W data in GALFIT with a Sérsic profile (with n = 0.5, which gives a Gaussian profile) to fit the diffuse emission and a point-like component (using our ePSF) to fit the more compact part of the galaxy (Fig. 3). The surface brightness of galaxy E is compatible with a nuclear, point-like source plus a diffuse emission component at radii larger than 0.05″. However, as the core remains unresolved in the F150W image, we extend our analysis to the F115W image, in which the S/N is lower, but the spatial resolution is higher than in the F150W image. The core of galaxy E also remains unresolved in the F115W image; thus, we obtained an upper limit for the effective radius of 63 pc (0.012″).

Fig. 3. NIRCam F150W images and encircled light profiles for B14-65666. Top panels: Galaxy E and residuals of the GALFIT model (modeled as one point source and one extended source). The white circular aperture represents 0.2 arcsec. Bottom panels: Galaxy W and residuals of the GALFIT model (modeled with three point sources and one extended one, marked with white crosses). The white circular aperture represents 0.28 arcsec. The white circles shown in the bottom left of the top and bottom images represent the NIRCam F150W PSF FWHM. Top and bottom right panels: Brown dashed curve shows the decline of the ePSF (from our observations) profile. The red and green solid lines represent the surface brightness profile of the individual galaxies E and W, respectively. The black line represents the surface brightness profile of the model of galaxies E and W, respectively, with GALFIT.

Galaxy W has a clumpy elongated morphology and is resolved in the F150W filter. It is well modeled by four different clumps. Three of them are unresolved point sources, which we modeled in GALFIT with point-like components. The other one is a more diffuse clump, associated with the brightest emission, modeled with point-like and Sérsic-profile (with an index n = 0.5) components (Fig. 3). When considering galaxy W as a single entity and circularizing the profile, the intrinsic effective radius is 348 pc (0.057″) after deconvolving with the NIRCam F150W PSF.

The merger nature of B14-65666 is supported by the morphological analysis using the statmorph code (Rodriguez-Gomez et al. 2019). We analyzed galaxies E and W separately and find that in both cases the shape asymmetry, A_S – a reliable tracer of merging galaxies (Rinaldi et al. 2024a) – is greater than 0.2, reinforcing the idea that each galaxy exhibits significant asymmetry. Specifically, we find A_S of 0.35 and 0.80 for galaxies E and W, respectively. We also tested the multimode-intensity-deviation (MID) statistics, which are useful for detecting multicomponent systems, such as the double-nucleus of a late-stage merger, highly disordered post-merger remnants, galaxies with bright star-forming clumps in rest-UV emission, or an apparently single galaxy with an extended emission component(s) (Freeman et al. 2013). M (multimode statistics) quantifies how a galaxy’s light is distributed into distinct regions, and D (deviation statistic) measures the offset between the brightest region and the overall centroid. M values approaching 1 are likely to indicate a double nucleus, while values close to zero are interpreted as indicating a single source; D values close to zero indicate symmetric and ordered morphologies such as spheroids and disks. In the case of galaxy E, M ≈ 0.0 and D ≈ 0.1. This corroborates the interpretation that galaxy E is a single (compact spheroidal) source. In the case of galaxy W, M ≈ 0.4 and D ≈ 0.3, indicating that galaxy W is composed of at least two clumps, with one of them dominating the emission.

4. Results and discussion

We detected and spatially resolved the rest-frame optical continuum and Hα emission in the system B14-65666 at z = 7.15 with MIRI imaging and spectroscopy. The new MIRI data are combined with ancillary NIRCam imaging and ALMA spectroscopy to obtain the physical properties of galaxies E and W, and of the entire system: star formation rate, instantaneous SFR, ionizing photon-production efficiency, Hα equivalent width, dynamical mass, stellar-mass surface density, and SFR surface density.

4.1. Spectral energy distribution and star formation history

We performed an SED fitting using CIGALE (Noll et al. 2009; Boquien et al. 2019) for galaxies E+W as well as for the two individual components. The SED covers the 1 μm to 1.3 mm observed wavelength range, combining JWST photometric data (NIRCam and MIRIM) with ALMA continuum imaging (Hashimoto et al. 2019; Sugahara et al. 2021), along with our Hα and [O III] 88 μm (Hashimoto et al. 2019) line fluxes. The redshift used in the analysis is the one derived spectroscopically from the Hα line.

We made different assumptions to represent the star formation history of the system. The best results are determined from CIGALE fitting by minimizing the χ² value and ensuring the best agreement with observed photometry and line fluxes. These are obtained when considering two components: an instantaneous burst with ages ranging from 30 to 500 Myr to model the mature stellar population; and a constant burst with ages lower than 10 Myr to model the young stellar population. We used the Bruzual & Charlot (2003) stellar population models and the Chabrier (2003) initial mass function (IMF). The nebular continuum and emission lines were modeled using electron densities from 100 to 1000 cm⁻³ and the ionization parameter (logU) ranging from –3.5 to –1.5. The metallicity was treated as a free parameter, ranging from 0.03 Z_⊙ to 0.7 Z_⊙, as previous studies showed that it varies in galaxies E and W (Sugahara et al. 2025; Jones et al. 2024). We used the dust attenuation law from Calzetti et al. (2000), and the dust emission was modeled using the Draine et al. (2014) models.

The best-fit models for galaxies E+W and for galaxies E and W separately are shown in Fig. 4. These models reveal that B14-65666 is dominated by a young stellar burst characterized with a total SFR of 127 ± 66 M_⊙ yr⁻¹ and a mass for the young stellar component of (8 ± 1) × 10⁸ M_⊙. Galaxies E and W have SFRs of 84 ± 49 and 50 ± 12 M_⊙ yr⁻¹, respectively, and show the same young stellar component with a stellar mass of (4 ± 1) × 10⁸ M_⊙. The SED fitting results indicate the presence of a mature stellar population in B14-65666 with a stellar mass of around 0.5 − 2 × 10⁹ M_⊙ and average ages about 200 Myr. SED fitting results are presented in Table 2. We conclude from the SED fitting that the B14-65666 system consists of two galaxies experiencing starbursts characterized by a similar age (see Sect. 4.9 for a more detailed discussion). Also, despite the compactness of galaxy E (see Sect. 4.7), the SED based on the NIRCam and MIRI imaging and covering the rest-frame 0.16–0.78 μm does not correspond to that of the average little red dots at redshifts ∼6–8 (Pérez-González et al. 2024b; Kokorev et al. 2024; Rinaldi et al. 2024a).

Fig. 4. Observed near-infrared-to-millimeter SEDs and best-fit models derived from CIGALE SED fitting analysis of galaxies E+W, E, and W. Gray line: Best-fit model for the galaxies E+W. Red line: Best-fit model for galaxy E. Green line: Best-fit model for galaxy W. The circles represent the observed photometric fluxes for galaxies E+W, E, and W, with the associated uncertainties given in Table 1. The crosses represent the modeled fluxes. Also indicated are the Hα observed and predicted fluxes (top right) and the corresponding [O III] 88 μm fluxes (bottom right).

4.2. Nebular and stellar attenuation

Dust attenuation plays a key role in shaping the observed properties of B14-65666, as revealed by MIRI, NIRCam, and ALMA observations. Combining the MRS Hα and the NIRSpec Hβ fluxes², the Hα/Hβ ratio gives observed values of 4.5 ± 0.6 and 2.9 ± 0.4 for galaxies E and W, respectively. We note that these values correspond to the narrow component of the line emission for both galaxies as the broad component present on the Hβ line (Jones et al. 2024) is not detected in Hα with the present data (see Sect. 4.5). Assuming a case B recombination, T_e ≈ 1.5 × 10⁴ K and n_e = 10³ cm⁻³, the theoretical ratio is 2.8. Using Cardelli et al. (1989), we obtain that the nebular visual extinction (A_V) is A_V = 1.5 ± 0.4 mag for galaxy E, and it is A_V = $0.1_{-0.1}^{+0.5}$ mag for galaxy W. Thus the nebular attenuation is higher for galaxy E, in agreement, within the uncertainties, with the values of A_V = 1.6 ± 0.3 mag for galaxy E and A_V = 0.7 ${}_{-0.2}^{+0.3}$ mag for W from Sugahara et al. (2025) (obtained with SED fitting and applying the conversion factor to nebular attenuation, f = 0.44; Calzetti 1997).

We also calculated the dust attenuation in the FUV, A_FUV, with the IR-to-UV luminosity ratio (IRX = L_IR/L_UV), as this is a robust tracer of A_FUV. Taking the value of log_IRX = 0.7 for the total galaxy (Sugahara et al. 2025), we adopted expression (5) from Meurer et al. (1999) assuming the bolometric correction for the total light emitted by dust is 1, and for the total light emitted by stars, UV, is 1.68 (Meurer et al. 1999). We obtain a value of A_FUV = 1.50, which is compatible with that obtained with our SED fitting for galaxies E+W, A_FUV = 1.41 ± 0.09 mag (Calzetti et al. 2000).

4.3. Star formation rate and burstiness

We derived the recent (i.e., less than 10 Myr) star formation rate directly from our Hα measurements, assuming solar (to facilitate comparison with other star formation tracers at other wavelengths) and subsolar metallicity (0.29 Z_⊙, Jones et al. 2024 and see Table 2). We applied the corresponding SFR-to-Hα conversion factor for each metallicity, which is approximately 1.7 times lower in the subsolar case (Reddy et al. 2018). The final SFR values are corrected for attenuation using the Hα/Hβ ratios (see previous section). We obtain SFRs of 76  ±  8 (46  ±  5) M_⊙ yr⁻¹ and 30  ±  4 (18  ±  3) M_⊙ yr⁻¹ for galaxies E and W and solar (subsolar) metallicities, respectively.

Additional estimates of the SFR can be obtained directly from far-infrared emission lines or using a combination of multiwavelength tracers. The [C II] 158 μm emission line is a tracer of SFR in low-redshift, star-forming, and composite SF+AGN galaxies (De Looze et al. 2014); it is confirmed for high-redshift galaxies (Schaerer et al. 2020; Carniani et al. 2020). Adopting the ancillary [C II] 158 μm luminosities derived from ALMA (Hashimoto et al. 2019), values similar to those derived from our Hα measurements are obtained, within the uncertainties (see Table 2). The large uncertainty in the value of the SFR derived from the [C II] 158 μm emission line is due to the nature of this emission, which could include extended or diffuse regions not directly linked to active star formation. Also, to trace less massive stars we used the SFR(UV) for solar metallicities and considering a constant SFH over 100 Myr (Calzetti 2013).

Other SFR estimates involve the use of hybrid indicators such as the UV+TIR (total infrared) luminosity (Kennicutt & Evans 2012), or Hα+TIR (Catalán-Torrecilla et al. 2015; Calzetti 2013), that capture both obscured and unobscured star formation over longer periods of time (up to 100 Myr). Assuming a simple energy balance (Kennicutt et al. 2009), and considering the far-UV luminosity derived from our SED fitting (L(UV) = (1.15 ± 0.11) × 10¹¹ L_⊙), and the TIR luminosity L(TIR) = (1.05 ± 0.2) × 10¹² L_⊙ (Hashimoto et al. 2019), we obtain an SFR consistent with our extinction-corrected Hα estimates. The SFR derived from the Hα and TIR luminosities are also in good agreement with those obtained from UV+TIR (see Table 2). The SFRs derived in this analysis (70–140 M_⊙ yr⁻¹) are lower than the SED-based SFR reported by Sugahara et al. (2025), SFR = 225 ${}_{-56}^{+71}$ M_⊙ yr⁻¹, and the (narrow+broad) Hβ-based SFR from Jones et al. (2024), SFR_Hβ = 308 ± 39 M_⊙ yr⁻¹. This difference is attributed to variations in aperture sizes where the contribution of the extended emission (Hβ) is significant, as well as other factors such as different estimates of the internal dust attenuation. Additionally, in the case of the Hβ-based SFR, the presence of both narrow and broad components in the Hβ line emission also contributes to the higher SFR value by a factor of 2 (Jones et al. 2024).

A useful diagnostic for assessing the recent star formation in high-z galaxies is the ratio of the SFR(Hα)-to-SFR(UV), which is referred to as the burstiness parameter and probes star formation over the last 10 Myr relatively to the last 100 Myr. If the ratios are above 1, it indicates that the galaxy is in a bursty phase (Atek et al. 2024). We obtain a value of the ratio SFR(Hα)-to-SFR(UV) for solar metallicities of 1.8 ± 0.1 and 0.8 ± 0.1 for galaxies E and W, respectively (we would expect similar values for subsolar metallicities). When compared to the range of burstiness ratios reported by Atek et al. (2024) (0.7 < z < 1.5), values above or close to unity are typical in low-mass galaxies (M_⋆ < 10⁹ M_⊙) and below unity in high-mass galaxies (M_⋆ > 10⁹ M_⊙). Thus, our values suggest that galaxy E shows signs of bursty activity, forming stars at nearly twice the average rate of the past 100 Myr; it is in a mild bursty phase, while galaxy W forms stars at a rate closer to that if this past 100 Myr; i.e., experiencing a more continuous burst. This aligns with the results from Navarro-Carrera et al. (2024), whose authors analyzed ∼4500 Hα emitters across z ∼ 2–6.5 using JADES data, observing a correlation between enhanced SFR(Hα)/SFR(UV) and EW₀(Hα), reinforcing the scenario where galaxies with above-unity SFR(Hα)/SFR(UV) experience a bursty phase in the last ∼50 Myr.

4.4. Equivalent width of Hα and ionizing photon-production efficiency

To derive the rest-frame equivalent width of Hα (EW₀(Hα)) for galaxies E+W and galaxies E and W, we make use of the F560W and Hα fluxes measured in our MIRIM and the MRS observations, respectively. First, we obtain rest-frame optical continuum emission by subtracting the contribution of the Hα line (Sect. 3.2) to the F560W flux (Sect. 3.3). The EW₀(Hα) is then derived as the ratio of the Hα to the continuum flux under the Hα line. We obtain rest-frame equivalent widths of 832 ± 100 Å and 536 ± 78 Å for galaxies E and W, respectively. The difference in the EW₀(Hα) of the two galaxies is attributed to a different star formation history in the two galaxies with ionizing young stars and more mature, non-ionizing stellar populations contributing in different proportions to the continuum emission. According to our SED analysis, galaxy E, with a larger EW₀(Hα) value, has a slightly younger stellar population (6 Myr) relative to galaxy W (8 Myr).

The ionizing photon production efficiency, log(ζ_ion), is the ratio of the ionizing to the non-ionizing UV flux. Assuming recombination, log(ζ_ion) is empirically obtained as the ratio of the number of ionizing photons (N_LyC) derived from the Hα line (after extinction correction) and the luminosity at 1500 Å. Assuming zero escape fraction, a dust attenuation as derived in Sect. 4.2, and a temperature of 1.5 × 10⁴ K (Álvarez-Márquez et al. 2024) the integrated Hα luminosity provides a total number of ionizing photons (N_LyC) equal to 1.1 × 10⁵⁵ ph s⁻¹ and 4.1 × 10⁵⁴ ph s⁻¹ for galaxies E and W, respectively. Different approaches to dust correction can yield significantly different values of ζ_ion introducing uncertainty in its determination and affecting the interpretation of its physical implications (Matthee et al. 2017). Thus, correcting the UV and Hα luminosities, we obtain ionizing photon-production efficiencies log(ζ_ion) equal to 25.1 ± 0.1 and 25.5 ± 0.1 Hz erg⁻¹ for galaxies E and W, respectively. This difference in the ζ_ion between the two galaxies of the system can be interpreted as partly due to the age difference in the young stellar populations in the two galaxies. Models indicate that in addition to the stellar continuum, there is a strong nebular contribution due to the two-photon continuum at 1500 Å. This contribution can be close to 50% of the total for very young populations (1 Myr), decreasing with age and becoming almost negligible for 10-Myr-old populations (Katz et al. 2025). Thus, the impact of this nebular continuum would be to increase the 1500 Å luminosity and consequently decrease the corresponding log(ζ_ion) by factors of 0.15–0.04 for ages of 1–8 Myr. In addition, residual differential extinction in the two galaxies could also play a role due to uncertainties in the derivation of the extinction in the UV.

We present the previous results in the log(ζ_ion)–EW₀(Hα) plane (see Fig. 5). Additional samples of galaxies at redshifts 4–9 are also presented for comparison. These include Lyman-alpha emitters (LAE) at redshifts 3.8–6 (Maseda et al. 2020; Ning et al. 2023; Lam et al. 2019; Harikane et al. 2018), Hα and [O III] emitters at redshifts 6.1–8 (Sun et al. 2023; Rinaldi et al. 2023), and the mean of Lyman-break galaxies at redshift 8 (Stefanon et al. 2022). Individual galaxies with resolved integral-field spectroscopy are also shown (RXCJ0600-2007 at z = 6.07, Giménez-Arteaga et al. 2024 and MACS1149-JD1 at z = 9.1, Álvarez-Márquez et al. 2024). While galaxies E and W forms stars at a rate much higher than the average galaxy at redshift 7, they have log(ζ_ion) and EW₀(Hα) in the range of high-z galaxies. These values are consistent with the predicted values for young stellar populations. The Binary Population and Spectral Synthesis (BPASS; Eldridge & Stanway 2020) model demonstrates that the log(ζ_ion) and the EW₀(Hα) of a single-aged, instantaneous-starburst binary population increases as the age of the population decreases (see Fig. 5), with values compatible with those observed for young populations with ages from 1 to 20 Myr. The production of ionizing photons also depends on the presence of massive stars and therefore the IMF (Stanway et al. 2016). Values for the logarithm of the photon production efficiency above 25.5 Hz erg⁻¹ can only be obtained with massive bursts, 10⁷ M_⊙ due to the stochasticity of the formation of massive stars (Stanway & Eldridge 2023). The galaxies in the B14-65666 system and the high-z galaxies identified with the JWST should have higher young stellar masses.

Fig. 5. log(ζion) as function of Hα equivalent width, both of them corrected from dust attenuation. The red star represents galaxy E of B14-65666 and the green star galaxy W of B14-65666. The blue line represents the relation for galaxies at redshifts 3–7 (Prieto-Lyon et al. 2023). B14-65666 is compared with the galaxy MACS1149-JD1 (Álvarez-Márquez et al. 2024), RXCJ0600-z6-3 (Giménez-Arteaga et al. 2024), mean of 102 Lyman-break galaxies at redshift 8 (Stefanon et al. 2022), MIDIS Hα-emitters at redshifts 7–8 (Rinaldi et al. 2023), faint Lyα emitters (LAEs) at redshifts 4–5 (Lam et al. 2019), LAEs at 3.8–5 (Maseda et al. 2020), a LAE subsample at redshift 4.9 (Harikane et al. 2018), line emitters at redshifts 6.11–6.35 (Sun et al. 2023), and LAEs at redshift 6 (Ning et al. 2023). We also present the BPASS models of a binary stellar population for an instantaneous burst with ages from 1 to 20 Myr (crosses; from red to blue, the ages are 20, 16, 13, 10, 8, 6, 4, 3, 2, and 1 Myr and are connected with the solid black line) (Stanway & Eldridge 2023).

4.5. Kinematics of the ionized ISM

The kinematics of the ionized gas can be measured with the Hα emission-line profile in our MRS data. Based on the integrated Hα line, we measure a redshift of z = 7.1513 ± 0.0007 for the B14-65666 system, compatible with the values obtained with the FIR emission lines [O III] 88 μm and [C II] 158 μm (z = 7.1521 ± 0.0004, Hashimoto et al. 2019) and optical [O III] line (z = 7.1517 ± 0.0001, Jones et al. 2024). The two galaxies of the B14-65666 system have a relative velocity of 175 ± 28 km s⁻¹ along our line of sight as derived from the centroid of the Hα line emission, showing it is a galaxy merger.

The kinematics of the ionized gas traced by Hα appears to show differences between the two galaxies. Specifically, we find that the FWHM for galaxy W is higher than for galaxy E, with values of 312 ± 44 km s⁻¹ and 243 ± 41 km s⁻¹, respectively. These values agree within the uncertainties (from those measured) with the far-infrared emission lines, [C II] 158 μm and [O III] 88 μm. The [C II] 158 μm emission has similar kinematics in both galaxies (FWHM of 288 ± 41 km s⁻¹ for galaxy W and 276 ± 34 km s⁻¹³ for galaxy E). For [O III] 88 μm, galaxy E shows a broader line (322 ± 41 km s⁻¹) than galaxy W (267 ± 42 km s⁻¹). In Fig. 6, we present the comparison of these three lines for galaxies E+W and E and W. The observed differences in the two galaxies indicate the complex structure of the velocity field in the different phases of the ISM. While the Hα line traces the overall kinematics of the ionized gas, the [C II] 158 μm line is closely associated with the neutral gas and the [O III] 88 μm with the high excitation gas. Simulations also show that different tracers are sensitive to different kinematic properties (Kohandel et al. 2024). In addition, the continuum imaging (Sugahara et al. 2025) of both galaxies shows a different morphology; galaxy E as a very compact source, while galaxy W shows a more extended structure, in particular in the rest-UV light, reminiscent of tidal tails. Moreover, low-spectral-resolution (R ≈ 1000) integral-field spectroscopy with NIRSpec suggests the presence of a narrow and a broad kinematic component in both galaxies and in all optical emission lines tracing the ionized gas (Jones et al. 2024). These two components are potentially associated with the motions in the host galaxy (narrow) and the presence of high-velocity gas outflows (broad). The width of the narrow component is only 144 ± 11 km s⁻¹ for galaxy E and 227 ± 6 km s⁻¹ for galaxy W. These are significantly narrower than the values derived from our Hα line and far-infrared emission lines, obtained with a better spectral resolution but lower S/N. We tested whether the observed Hα line profile could be compatible with the combination of the narrow and broad (653 ± 10 km s⁻¹ for galaxy E and 755 ± 23 km s⁻¹ for galaxy W) line components identified with NIRSpec. Our tests (see Fig. 7) indicate that the narrow+broad (N+B) component could be slightly compatible with our observed Hα line profiles. However, the simulated Hα N+B line profile appears to have a flux excess in the spectral range of the line (i.e., residuals mostly below zero) for galaxy E, while the faint broad component in galaxy W is close in velocity to the peak emission in galaxy E. Since the two galaxies are separated by only 0.44 arcsec, this secondary emission in W could be due to a contamination from the brighter emission from E (Jones et al. 2024). Further higher spectral resolution data with NIRSpec, and higher S/N MRS spectra, would be required to confirm the presence of such a broad line emission in the two galaxies.

Fig. 6. Comparison of integrated Hα emission-line profile, with the ALMA spectrum of the [O III] 88 μm line and [C II] 158 μm line (Hashimoto et al. 2019) for galaxies E+W, for galaxy E, and for galaxy W, respectively, from top to bottom. The systemic velocity corresponds to redshift 7.1513.

Fig. 7. Direct comparison of MRS-observed Hα emission-line profile for galaxies E (top) and W (bottom), with the rescaled Hβ narrow- and broad-line components identified in lower spectral resolution NIRSpec IFS (Jones et al. 2024). The presence of the narrow plus broad components could be compatible with the observed Hα line profile (see Sect. 4.5 for a detailed discussion). The systemic velocity corresponds to redshift 7.1513. The residuals are shown in the bottom plots, together with the 1σ noise level (yellow).

The consistency of the trends between our results and the ancillary data provides robust evidence of the dynamic differences between galaxies E and W. The broader profiles observed in galaxy W for both Hα and [C II] 158 μm may indicate enhanced turbulence or kinematic effects, potentially driven by interactions or feedback processes, as well as a possibility of a connection with outflows. On the other hand, the higher FWHM in [O III] 88 μm for galaxy E suggests that highly ionized gas dynamics differ from those of lower ionization states, pointing to localized physical conditions or energetic phenomena. Also, dust attenuation could impact the interpretation, in particular in galaxy E where a large visual extinction is measured. The optical lines are more susceptible to dust extinction, which can affect measurements of line widths and kinematic features, whereas far-infrared lines are less impacted by such effects.

4.6. Stellar and dynamical mass

Although the line profiles of the different emission lines suggest that the gas is not fully (dynamically) relaxed, here we obtain an upper limit of the dynamical mass of each of the galaxies assuming the system is close to relaxation and characterized by the random motions traced by the ionized gas velocity dispersion; i.e., not dominated by rotational motions. Under this working hypothesis, we followed previous works on similar systems such as the low-z luminous infrared galaxies (e.g., Eq. (1) in Bellocchi et al. 2013). The dynamical masses we derive are (1.5 ± 0.3) × 10⁹ M_⊙ and (8.5 ± 2.4) × 10⁹ M_⊙ for galaxies E and W, respectively, assuming velocity dispersions of 103 ± 17 km s⁻¹ and 132 ± 19 km s⁻¹ (obtained with the FWHM from our Hα line), and the effective radii derived from our NIRCam images (see Sect. 3.4 for details). The total stellar masses are derived from the multiwavelength SED fitting and correspond to (15 ± 9) × 10⁸ M_⊙ and (8 ± 4) × 10⁸ M_⊙ for galaxies E and W. Thus, while there is a good agreement between the stellar and dynamical masses for galaxy E, an order of magnitude difference appears for galaxy W. The difference in galaxy W could be due to the presence of a dominant molecular gas component in this galaxy amounting to close to 90% of the total mass. However, the stellar structure shows the presence of extended emission consisting of star-forming clumps along tidal tails (Sugahara et al. 2025). This suggests that this galaxy is far from relaxed, and therefore the estimated dynamical mass is highly uncertain.

4.7. Mass-size relation

The two galaxies of the B14-65666 system occupy a very different location in the stellar mass-size plane (see Fig. 8). Galaxy W, with an effective radius of 348 pc and a stellar mass of (8 ± 4) × 10⁸ M_⊙, is consistent with the expected values according to the relations derived for high-redshift galaxies from the JADES program (redshifts 4–10; Langeroodi & Hjorth 2023) and from the CEERS, PRIMER-UDS, and PRIMER-COSMOS surveys (redshifts 3–9; Allen et al. 2025). Therefore, galaxy W appears as an average galaxy at approximately redshift 7. On the other hand, the second component of the system, galaxy E, appears to be deviating from the standard relations, showing a very small and compact size (upper limit of 63 pc to the effective radius) for a stellar mass of 15 (± 9) × 10⁸ M_⊙. This translates into a high-stellar-mass surface density (6 × 10⁴ M_⊙ pc⁻² for E, versus 1 × 10³ M_⊙ pc⁻² for W), close to the limit of the highest values (10^5.5 M_⊙ pc⁻²) measured in the stellar clusters and nuclei of low-z galaxies (Grudić et al. 2019). Galaxies as compact as galaxy E and with similar stellar surface densities are not uncommon at high-z (Álvarez-Márquez et al. 2025). Star-forming clumps and galaxies with an effective radius in the 30–200 pc range and a wide mass range (4 × 10⁷ to 2 × 10¹⁰ M_⊙) are already known at redshifts above 6 (see Fig. 8; Claeyssens et al. 2023; Baggen et al. 2023). Thus, B14-65666 is identified as a system that consists of a pair of interacting galaxies, with similar stellar masses (within a factor of 2) but with one of the galaxies (E) being extremely compact and already showing stellar-mass surface densities close to those measured in the nuclei of galaxies at low redshifts. The galaxies of the B14-65666 system are much less massive (≈ × 20–30) and far more compact (≈ × 10–60) than the main-sequence and most luminous star-forming galaxies at low redshifts, with effective radii ≈ 3.5 kpc (Leroy et al. 2021) and ≈2.2 kpc (Bellocchi et al. 2022), respectively (see Fig. 8). Thus, galaxies E and W appear to be building up their stellar mass as they go through the advanced phases of the interaction/merging process.

Fig. 8. Mass-size relation for galaxies E (red star) and W (green star) in B14-65666, together with clumps of SMACS0723 at redshifts 6.4–8.5 (yellow squares, Claeyssens et al. 2023) and a subsample of massive and compact galaxies from the CEERS program (purple circles, Baggen et al. 2023). The blue line represents the mass-size relation for galaxies (blue dots) at redshifts 4–10 (Langeroodi & Hjorth 2023), while the solid black line indicates the mass-size relation for redshifts 6–9 derived for the NIRCam filter F150W (Allen et al. 2025). As B14-65666 is an interacting system, we also represent the mean values of the mass and size for local LIRGs and MS SFGs (Bellocchi et al. 2022) and (Leroy et al. 2021). Dotted gray lines represent constant stellar mass surface density for 10, 102, 103, 104, and 105 M⊙ pc−2. The red line indicates the observed maximum value of stellar-mass surface density in clusters and nuclei of galaxies, 105.5 M⊙ pc−2.

4.8. Starburst nature

The analysis of the SFR derived with the Hα-versus-stellar-mass plane reveals that both galaxies E and W are well above the main sequence (MS) of galaxies at redshift 7 (Heintz et al. 2023). Galaxies E and W of our system are situated between the region defined by Caputi et al. (2017) for Hα emitters, the starburst (SB) region, and the MS (see Fig. 9). Such extremely compact starburst galaxies are expected to have short depletion times (< 100 Myr, see Franco et al. 2020), and therefore galaxy E may quickly quench out of gas exhaustion – at least temporarily – which is in line with a bursty/stochastic SFH.

Fig. 9. Stellar-mass recent SFR diagram. The large red star represents galaxy E, and the large green star represents galaxy W. We use the SFR traced by Hα for solar metallicities. The orange line delimits the starburst area (Caputi et al. 2017). The blue line represents the MS for galaxies derived for redshift 7, and the green line shows the one for galaxies at redshift 0 (Heintz et al. 2023; Elbaz et al. 2011). The dashed gray lines represent the lines of constant sSFR. We also present galaxies of redshifts 6–8 (Endsley et al. 2023) (blue dots), the main-sequence galaxies of ALPINE at redshifts 4.4–5.9 (Faisst et al. 2020; Le Fèvre et al. 2020) (olive dots), strong Hα-emitters (HAEs) at z ≈ 7–8 (Rinaldi et al. 2024b) (yellow dots) and the mean for the LIRG sample (Bellocchi et al. 2022) (pink dot).

We compared the position in the SFR stellar-mass plot with more galaxies at similar redshifts (Fig. 9). The sSFRs of E and W (51 and 38 Gyr⁻¹, respectively) are similar to those of low-mass galaxies, but also to some very massive ones, such as REBELS-25 (sSFR = 25 Gyr⁻¹), a relaxed, dynamically cold disc galaxy at a similar redshift, z = 7.31 (Rowland et al. 2024). A recent large compilation of NIRSpec prism data and stacked spectra showed that the average galaxy becomes progressively more bursty with increasing redshift, as well as with increasing sSFR (Langeroodi & Hjorth 2024). Based on this, B14-65666 seems slightly less bursty (but consistent within error bars) than the average galaxy at similar redshift and sSFR. UV-faint galaxies at redshifts from 6 to 8 from the Cosmic Evolution Early Release Science (CEERS) JWST/NIRCam survey also show a bursty star formation history (Endsley et al. 2023). Galaxies from the ALMA-ALPINE [C II] survey at redshifts 4–6 (Le Fèvre et al. 2020; Faisst et al. 2020) show normal star formation following the MS of galaxies at these redshifts (Speagle et al. 2014). LIRGs at low redshifts, with typical stellar masses in the 10¹⁰ − 10¹¹ M_⊙ and SFR 30 M_⊙ yr⁻¹ range (Bellocchi et al. 2022), are placed above the low-z MS plane (Elbaz et al. 2011) and indicate that at high redshifts LIRGs are less massive but more compact.

4.9. B14-65666: A merger of starburst galaxies in the early Universe

The picture that emerges from our analysis is that B14-65666 is a system of two galaxies in the process of interacting and merging. We identified one of the components (E) as an extremely compact source, characterized by a stellar-mass surface density close to that of the nuclei of low-z galaxies. On the other hand, the second component of the system (galaxy W) shows a more elongated structure with four fainter clumps (Sugahara et al. 2025), reminiscent of star-forming regions along the tidal tails produced in gas-rich interacting galaxies. A significant velocity offset of 175 km s⁻¹ between the two galaxies supports the merger scenario. Furthermore, the estimated stellar masses for each component are on the order of ∼10⁹ M_⊙, which is three-to-four orders of magnitude higher than typical stellar masses of star-forming clumps in galactic disks (10⁵–10⁶ M_⊙, Adamo et al. 2024; Mowla et al. 2024), reinforcing the interpretation that we are observing two galaxy-scale systems. The stellar-mass ratio between the components is ∼2:1, which is potentially consistent with a major merger. Additionally, we find differences in metallicity between the two components, with values derived using the SED-fitting of 0.24 Z_⊙ for E and 0.31 Z_⊙ for W; this is in agreement with line diagnostics from NIRSpec IFS (Jones et al. 2024). Collectively, the presence of tidal features, velocity offset, distinct stellar populations and metallicities, and the high stellar masses all provide compelling evidence that B14-65666 is an example of a merging system of two galaxies. In addition, galaxies E and W are experiencing a bursty phase with an SFR well above that of the main-sequence star-forming galaxies of the same mass at those redshifts.

Recent simulations of the formation of galaxies during the EoR (Nakazato et al. 2024, and references therein) show that systems identified as clumpy – i.e., with two or more high-surface-density star-forming clumps as traced by emission lines – represent 10% of all simulated sources. Moreover, most of the clumpy systems are produced by mergers, with the young stellar populations forming along the gas-rich tidal tails, and with specific SFR of the order of 50 Gyr⁻¹. However, these clumpy systems appear to be temporal and short-lived, merging into a single galaxy in about 50 Myr (Hashimoto et al. 2023b). Under this scenario, the compact, and highly dense, galaxy E could already be a protobulge produced by the merger of previous clumps, while galaxy W would be in the process of merging of smaller and less massive clumps before the two main components, E and W, merge together in a single, more massive galaxy at a later time; i.e., at lower redshifts (see Figure 5 in Nakazato et al. 2024).

B14-65666 is not a unique system at those redshifts. Recent JWST observations show merging systems of star-forming galaxies such as the interacting system SPT0311-58, which resides in the core of an extremely massive protocluster at z = 6.9 (Arribas et al. 2024); the source A2744-z7p9OD, the most distant protocluster at z = 7.88 (Hashimoto et al. 2023b); or the JWST’s Quintet, a merger of at least five galaxies containing over 17 clumps, at z = 6.7 (Hu et al. 2025). The detection with the JWST and further studies of new systems such as these will allow us to understand the formation process of massive galaxies in the early Universe.

5. Conclusions and summary

We present mid-infrared imaging and spectroscopy of the z = 7.15 Lyman-break galaxy B14-65666 taken with MIRI/JWST. The data provide the first detection of the Hα emission line of this source, a system of two merging galaxies at the EoR. The study also combines the new MIRI data with NIRCam data (Sugahara et al. 2025) and archival high-angular-resolution ALMA imaging (Hashimoto et al. 2018) to derive additional physical properties. The main results are summarized below.

• The Hα-emitting gas shows a structure dominated by two spatially resolved sources (E and W) separated by a projected distance of 0.4 arcsec (i.e., 2.2 kpc) and relative velocity of 175 km s⁻¹. The rest-frame UV light shows galaxy E as a very compact, unresolved source (upper limit for the effective radius of 63 pc), while galaxy W is more extended (348 pc) with the presence of fainter sources along an elongated structure.

• The SFR derived from the Hα luminosity is 76 and 30 M_⊙ yr⁻¹ (for solar metallicity as reference), for galaxies E and W, respectively. Other measurements of the SFR using hybrid tracers; Hα and the luminosity TIR; and the UV-plus-TIR luminosities give values of 50 and 54 M_⊙ yr⁻¹ and 66 and 88 M_⊙ yr⁻¹ for galaxies E and W, respectively. These values are compatible within the uncertainties with the SFR derived from Hα. The ratio of Hα-to-UV-based SFR (burstiness) indicates that while galaxy E is in a mild bursty phase (ratio of 1.8), galaxy W is experiencing a more continuous burst (ratio of 0.8).

• The EW of the Hα line (832 ± 100 and 536 ± 78 Å for galaxies E and W, respectively) and the ionizing photon production efficiency (log(ζ_ion) = 25.1 ± 0.1 Hz erg⁻¹ and log(ζ_ion) = 25.5 ± 0.1 Hz erg⁻¹) in the two galaxies are within the range of values measured in other samples of galaxies at these redshifts. Their different locations in the production efficiency – equivalent width plane, with galaxy W above and galaxy E below the relation for high-z galaxies is interpreted as being due to slightly different ages in their young stellar populations.

• The specific star formation rate (40–50 Gyr⁻¹) indicates that the two galaxies involved in the B14-65666 system are experiencing an intense starburst phase, well above the typical values for galaxies of similar masses at those redshifts and similar to less massive (× 10) galaxies.

• The stellar-mass surface density of galaxy E is extreme (> 6 × 10⁴ M_⊙ pc⁻²) and close to the value measured in the nuclei of low-z galaxies, while galaxy W has a surface density (10³ M_⊙ pc⁻²) similar to that of a high-z galaxy.

• The Hα emission-line profile shows a different gas kinematics, with galaxy W showing larger velocities (i.e., broader line) than galaxy E. The ionized gas kinematics traced by the Hα line shows differences with respect to the neutral and highly ionized gas traced by the far-infrared [C II] 158 μm and [O III] 88 μm lines, respectively. This indicates a complex structure in the ISM, likely due to the internal substructure of the young and mature stellar populations and dust distribution in these galaxies on (sub)kiloparsec scales.

In summary, B14-65666 is a system consisting of two galaxies that are experiencing an extreme starburst phase due to the interaction/merging process. One of the galaxies (E) is already very compact and massive, with a stellar density close to that of the mature nucleus/bulge of low-z galaxies, while the other one (W) is consistent with a more extended star formation, distributed in a few fainter clumps along a tidal tail. These results agree with recent cosmological simulations of clumpy merging systems in the early Universe (Nakazato et al. 2024). According to these simulations, this clumpy phase identified in B14-65666 is short-lived, with the bursty clumps merging on a short timescale (50 Myr) and forming a massive, single-nucleus galaxy at later times. Detailed studies of systems similar to B14-65666 will help our understanding of the physical conditions and mechanisms involved in the formation process of massive galaxies and their evolution in the early Universe, just a few hundred Myr after the Big Bang and during the EoR.

Acknowledgments

We thank Michele Perna and Gareth Jones for discussions on B14-65666. CPJ, JAM, LC acknowledge support by grant PIB2021-127718NB-100 from the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. ACG acknowledges support by JWST contract B0215/JWST-GO-02926. AB, JM and GÖ acknowledge support from the Swedish National Space Administration (SNSA). KIC acknowledges funding from the Dutch Research Council (NWO) through the award of the Vici Grant VI.C.212.036 and funding from the Netherlands Research School for Astronomy (NOVA). AAH acknowledges support from from grant PID2021-124665NB-I00 funded by the Spanish Ministry of Science and Innovation and the State Agency of Research MCIN/AEI/10.13039/501100011033 and ERDF A way of making Europe. This work was supported by research grants (VIL16599,VIL54489) from VILLUM FONDEN. JPP and TVT acknowledge financial support from the UK Science and Technology Facilities Council, and the UK Space Agency. The project that gave rise to these results received the support of a fellowship from the “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/PR24/12050015. PGPG and LC acknowledge support from grant PID2022-139567NB-I00 funded by Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033, FEDER Una manera de hacer Europa. EI acknowledges funding from the Netherlands Research School for Astronomy (NOVA). IJ acknowledge support from the Carlsberg Foundation (grant no CF20-0534) and the Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. SG acknowledges financial support from the Villum Young Investigator grant 37440 and 13160 and the Cosmic Dawn Center (DAWN), funded by the Danish National Research Foundation under grant DNRF140. AE and FP acknowledge support through the German Space Agency DLR 50OS1501 and DLR 50OS2001 from 2015 to 2023. The work presented is the effort of the entire MIRI team and the enthusiasm within the MIRI partnership is a significant factor in its success. MIRI draws on the scientific and technical expertise of the following organisations: Ames Research Center, USA; Airbus Defence and Space, UK; CEA-Irfu, Saclay, France; Centre Spatial de Liége, Belgium; Consejo Superior de Investigaciones Científicas, Spain; Carl Zeiss Optronics, Germany; Chalmers University of Technology, Sweden; Danish Space Research Institute, Denmark; Dublin Institute for Advanced Studies, Ireland; European Space Agency, Netherlands; ETCA, Belgium; ETH Zurich, Switzerland; Goddard Space Flight Center, USA; Institute d’Astrophysique Spatiale, France; Instituto Nacional de Técnica Aeroespacial, Spain; Institute for Astronomy, Edinburgh, UK; Jet Propulsion Laboratory, USA; Laboratoire d’Astrophysique de Marseille (LAM), France; Leiden University, Netherlands; Lockheed Advanced Technology Center (USA); NOVA Opt-IR group at Dwingeloo, Netherlands; Northrop Grumman, USA; Max-Planck Institut für Astronomie (MPIA), Heidelberg, Germany; Laboratoire d’Etudes Spatiales et d’Instrumentation en Astrophysique (LESIA), France; Paul Scherrer Institut, Switzerland; Raytheon Vision Systems, USA; RUAG Aerospace, Switzerland; Rutherford Appleton Laboratory (RAL Space), UK; Space Telescope Science Institute, USA; Toegepast- Natuurwetenschappelijk Onderzoek (TNO-TPD), Netherlands; UK Astronomy Technology Centre, UK; University College London, UK; University of Amsterdam, Netherlands; University of Arizona, USA; University of Cardiff, UK; University of Cologne, Germany; University of Ghent; University of Groningen, Netherlands; University of Leicester, UK; University of Leuven, Belgium; University of Stockholm, Sweden; Utah State University, USA. A portion of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We would like to thank the following National and International Funding Agencies for their support of the MIRI development: NASA; ESA; Belgian Science Policy Office; Centre Nationale D’Etudes Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luft-und Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economía y Competitividad; Netherlands Research School for Astronomy (NOVA); Netherlands Organisation for Scientific Research (NWO); Science and Technology Facilities Council; Swiss Space Office; Swedish National Space Board; UK Space Agency. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to the Author Accepted Manuscript version arising from this submission. 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 (AURA), Inc., under NASA contract NAS 5-03127 for JWST; and from the European JWST archive (eJWST) operated by the ESDC. These observations are associated with program ID 1284. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via https://doi.org/10.17909/gftt-2t47. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support to MAST for these data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2022).

References

All Tables

All Figures

Fig. 1. Top left: MIRI Image F560W. The white circumference represents a 0.9″ radius aperture. The solid black elliptical line represents the E+W galaxies’ apertures for the Hα emission of the whole galaxy; the dashed black lines represent the apertures for galaxies E and W. The white filled circle area in the bottom left represents the spatial resolution (PSF FWHM) of the MIRI image F560W. Top right: Hα line map. The Hα line map is generated by integrating Hα line emission in the velocity range: –500 < v [km s−1] < 300. The white filled area represents the spatial resolution (PSF FWHM) of the MRS. Bottom: NIRCam cutouts in six filters, F115W to F444W (left to right and top to bottom). The images are homogenized to the PSF of the MIRI image, represented with the white filled circle in the NIRCam F115W image.

In the text

Fig. 2. MIRI MRS 1 SHORT spectra for Hα observed at z = 7.1513. Top left: Hα spectra shown in velocity space for two spatially separated galaxies, E and W (red line and green line, respectively), identified in the Hα line map. The black line shows the integrated MRS spectrum of B14-65666 extracted using the galaxies E+W aperture. Top right: E+W galaxies’ integrated Hα spectrum. A vertical dashed line shows the wavelength of the weaker undetected lines [N II] 6550 Å and [N II] 6585 Å that are well separated in wavelength from Hα 6563 AA. Bottom left: Galaxy E’s integrated Hα spectrum. Bottom right: Galaxy W’s integrated Hα spectrum. The red dashed line shows the single Gaussian fits, with the residuals shown in the bottom panels. 1σ noise level is shown in yellow; it was calculated from the continuum adjacent to the emission line.

In the text

Fig. 3. NIRCam F150W images and encircled light profiles for B14-65666. Top panels: Galaxy E and residuals of the GALFIT model (modeled as one point source and one extended source). The white circular aperture represents 0.2 arcsec. Bottom panels: Galaxy W and residuals of the GALFIT model (modeled with three point sources and one extended one, marked with white crosses). The white circular aperture represents 0.28 arcsec. The white circles shown in the bottom left of the top and bottom images represent the NIRCam F150W PSF FWHM. Top and bottom right panels: Brown dashed curve shows the decline of the ePSF (from our observations) profile. The red and green solid lines represent the surface brightness profile of the individual galaxies E and W, respectively. The black line represents the surface brightness profile of the model of galaxies E and W, respectively, with GALFIT.

In the text

Fig. 4. Observed near-infrared-to-millimeter SEDs and best-fit models derived from CIGALE SED fitting analysis of galaxies E+W, E, and W. Gray line: Best-fit model for the galaxies E+W. Red line: Best-fit model for galaxy E. Green line: Best-fit model for galaxy W. The circles represent the observed photometric fluxes for galaxies E+W, E, and W, with the associated uncertainties given in Table 1. The crosses represent the modeled fluxes. Also indicated are the Hα observed and predicted fluxes (top right) and the corresponding [O III] 88 μm fluxes (bottom right).

In the text

Fig. 5. log(ζion) as function of Hα equivalent width, both of them corrected from dust attenuation. The red star represents galaxy E of B14-65666 and the green star galaxy W of B14-65666. The blue line represents the relation for galaxies at redshifts 3–7 (Prieto-Lyon et al. 2023). B14-65666 is compared with the galaxy MACS1149-JD1 (Álvarez-Márquez et al. 2024), RXCJ0600-z6-3 (Giménez-Arteaga et al. 2024), mean of 102 Lyman-break galaxies at redshift 8 (Stefanon et al. 2022), MIDIS Hα-emitters at redshifts 7–8 (Rinaldi et al. 2023), faint Lyα emitters (LAEs) at redshifts 4–5 (Lam et al. 2019), LAEs at 3.8–5 (Maseda et al. 2020), a LAE subsample at redshift 4.9 (Harikane et al. 2018), line emitters at redshifts 6.11–6.35 (Sun et al. 2023), and LAEs at redshift 6 (Ning et al. 2023). We also present the BPASS models of a binary stellar population for an instantaneous burst with ages from 1 to 20 Myr (crosses; from red to blue, the ages are 20, 16, 13, 10, 8, 6, 4, 3, 2, and 1 Myr and are connected with the solid black line) (Stanway & Eldridge 2023).

In the text

	Fig. 6. Comparison of integrated Hα emission-line profile, with the ALMA spectrum of the [O III] 88 μm line and [C II] 158 μm line (Hashimoto et al. 2019) for galaxies E+W, for galaxy E, and for galaxy W, respectively, from top to bottom. The systemic velocity corresponds to redshift 7.1513.
In the text

Fig. 7. Direct comparison of MRS-observed Hα emission-line profile for galaxies E (top) and W (bottom), with the rescaled Hβ narrow- and broad-line components identified in lower spectral resolution NIRSpec IFS (Jones et al. 2024). The presence of the narrow plus broad components could be compatible with the observed Hα line profile (see Sect. 4.5 for a detailed discussion). The systemic velocity corresponds to redshift 7.1513. The residuals are shown in the bottom plots, together with the 1σ noise level (yellow).

In the text

Fig. 8. Mass-size relation for galaxies E (red star) and W (green star) in B14-65666, together with clumps of SMACS0723 at redshifts 6.4–8.5 (yellow squares, Claeyssens et al. 2023) and a subsample of massive and compact galaxies from the CEERS program (purple circles, Baggen et al. 2023). The blue line represents the mass-size relation for galaxies (blue dots) at redshifts 4–10 (Langeroodi & Hjorth 2023), while the solid black line indicates the mass-size relation for redshifts 6–9 derived for the NIRCam filter F150W (Allen et al. 2025). As B14-65666 is an interacting system, we also represent the mean values of the mass and size for local LIRGs and MS SFGs (Bellocchi et al. 2022) and (Leroy et al. 2021). Dotted gray lines represent constant stellar mass surface density for 10, 102, 103, 104, and 105 M⊙ pc−2. The red line indicates the observed maximum value of stellar-mass surface density in clusters and nuclei of galaxies, 105.5 M⊙ pc−2.

In the text

Fig. 9. Stellar-mass recent SFR diagram. The large red star represents galaxy E, and the large green star represents galaxy W. We use the SFR traced by Hα for solar metallicities. The orange line delimits the starburst area (Caputi et al. 2017). The blue line represents the MS for galaxies derived for redshift 7, and the green line shows the one for galaxies at redshift 0 (Heintz et al. 2023; Elbaz et al. 2011). The dashed gray lines represent the lines of constant sSFR. We also present galaxies of redshifts 6–8 (Endsley et al. 2023) (blue dots), the main-sequence galaxies of ALPINE at redshifts 4.4–5.9 (Faisst et al. 2020; Le Fèvre et al. 2020) (olive dots), strong Hα-emitters (HAEs) at z ≈ 7–8 (Rinaldi et al. 2024b) (yellow dots) and the mean for the LIRG sample (Bellocchi et al. 2022) (pink dot).

In the text