Formation of a Low-Mass Galaxy from Star Clusters in a 600-Million-Year-Old Universe (PDF)

Article Formation of a low-mass galaxy from star clusters in a 600-million-year-old Universe https://doi.org/10.1038/s41586-024-08293-0 Lamiya Mowla1,2,15 ✉, Kartheik Iyer3,15 ✉, Yoshihisa Asada4,5, Guillaume Desprez4,...

Article Formation of a low-mass galaxy from star clusters in a 600-million-year-old Universe https://doi.org/10.1038/s41586-024-08293-0 Lamiya Mowla1,2,15 ✉, Kartheik Iyer3,15 ✉, Yoshihisa Asada4,5, Guillaume Desprez4, Vivian Yun Yan Tan6, Nicholas Martis7, Ghassan Sarrouh6, Victoria Strait8,9, Received: 11 February 2024 Roberto Abraham10, Maruša Bradač7, Gabriel Brammer8,9, Adam Muzzin6, Camilla Pacifici11, Accepted: 28 October 2024 Swara Ravindranath12, Marcin Sawicki4, Chris Willott13, Vince Estrada-Carpenter4, Nusrath Jahan14, Gaël Noirot4, Jasleen Matharu8,9, Gregor Rihtaršič7 & Johannes Zabl4 Published online: 11 December 2024 Open access The most distant galaxies detected were seen when the Universe was a scant 5% of its Check for updates current age. At these times, progenitors of galaxies such as the Milky Way were about 10,000 times less massive. Using the James Webb Space Telescope ( JWST) combined with magnification from gravitational lensing, these low-mass galaxies can not only be detected but also be studied in detail. Here we present JWST observations of a strongly lensed galaxy at zspec = 8.296 ± 0.001, showing massive star clusters (the Firefly Sparkle) cocooned in a diffuse arc in the Canadian Unbiased Cluster Survey (CANUCS)1. The Firefly Sparkle exhibits traits of a young, gas-rich galaxy in its early formation stage. The mass of the galaxy is concentrated in 10 star clusters (49–57% of total mass), with individual masses ranging from 105M⊙ to 106M⊙. These unresolved clusters have high surface densities (>103M⊙ pc−2), exceeding those of Milky Way globular clusters and young star clusters in nearby galaxies. The central cluster shows a nebular-dominated spectrum, low metallicity, high gas density and high electron temperature, hinting at a top-heavy initial mass function. These observations provide our first spectrophotometric view of a typical galaxy in its early stages, in a 600-million-year-old Universe. The Firefly Sparkle is a gravitationally lensed arc identified with the redshifts from the CANUCS dataset. The model shows magnification Hubble Space Telescope (HST) in the CLASH survey of the galaxy cluster factors between 16 and 26. A NIRSpec Prism slitlet, placed on the high- MACS J1423.8 + 2404 (hereafter MACS 1423) and reported as a z > 7 can- est magnification region at the centre of the arc, shows strong [Oiii] didate2. Follow-up spectroscopy using MOSFIRE on the Keck telescope emission, dominating the F444W flux and making the object appear suggested a redshift of z = 7.6 based on a possible Lyman-α (Lyα) detec- red in the composite image (Fig. 1). The projected half-light size of the tion3. The Canadian Unbiased Cluster Survey (CANUCS)1 revisited the arc in the source plane is only 0.3 ± 0.1 kpc, with most bright clusters field with JWST4, using NIRISS5, NIRCam6 and NIRSpec7. Imaging in 11 near the centre. Eight of the ten unresolved clusters (FF-3–FF-10) are bands (0.8–5 μm) showed a long magnified arc with distinct star clusters near the centre, whereas two others (FF-1 and FF-2) lie along an elon- embedded in a low surface brightness component extending up to 4′. gated arm, with a distance of 1.4 kpc between FF-1 and the central clus- NIRSpec Prism spectroscopy, covering the central brightest region, ter FF-5 (all distances quoted in the paper are projected distances in confirmed the high redshift (zspec = 8.296 ± 0.001) through multiple the source plane). Neighbour FF-BF is even more strongly magnified emission lines, with no Lyα emission detected. The Firefly Sparkle has ( μ = 28.014.4 −4.7), located within 2 kpc of FF-1 and also exhibits strong [Oiii] two neighbours: Firefly-Best Friend (FF-BF) at zspec = 8.2996 ± 0.0008 emission, whereas FF-NBF is at a distance of 13 kpc with very faint [Oiii] and Firefly-New Best Friend (FF-NBF) at zspec = 8.2967 ± 0.0016. All three emission. galaxies are shown in Fig. 1; this article focuses on the Firefly Sparkle We use NIRCam and NIRISS imaging to study the resolved Firefly and its star clusters. Sparkle and explore the stellar mass distribution in the clusters ver- The Firefly Sparkle resides in a highly magnified region lensed by the sus the diffuse arc. Photometry, derived by joint modelling of the 10 MACS 1423 cluster, enabling us to resolve the galaxy down to its indi- clusters and the diffuse arc using GALFIT10, shows that nine clusters vidual star clusters. We created a magnification model using Lenstool8,9, are unresolved, even in the highest resolution F115W images. Only the constrained by three multiple image systems3 with spectroscopic central cluster (FF-4) exhibits an elongated component and is fit with 1 Whitin Observatory, Department of Physics and Astronomy, Wellesley College, Wellesley, MA, USA. 2Center for Astronomy, Space Science, and Astrophysics, Independent University Bangladesh, Dhaka, Bangladesh. 3Columbia Astrophysics Laboratory, Columbia University, New York, NY, USA. 4Department of Astronomy and Physics, Saint Mary’s University, Halifax, Nova Scotia, Canada. 5Department of Astronomy, Kyoto University, Kyoto, Japan. 6Department of Physics and Astronomy, York University, Toronto, Ontario, Canada. 7University of Ljubljana, Department of Mathematics and Physics, Ljubljana, Slovenia. 8Cosmic Dawn Center (DAWN), Copenhagen, Denmark. 9Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark. 10 David A. Dunlap Department of Astronomy and Astrophysics, University of Toronto, Toronto, Ontario, Canada. 11Space Telescope Science Institute, Baltimore, MD, USA. 12Astrophysics Science Division, NASA Goddard Space Flight Center, Greenbelt, MD, USA. 13NRC Herzberg, Victoria, British Columbia, Canada. 14Shahjalal University of Science and Technology, Sylhet, Bangladesh. 15 These authors contributed equally: Lamiya Mowla, Kartheik Iyer. ✉e-mail: [email protected]; [email protected] 332 | Nature | Vol 636 | 12 December 2024 a F115W + F150W F200W + F277W F356W + F444W b F115W F277W F444W 20 1 25 2 5 4 10 kl3 35 ark Firefly Sparkle 6 3 5 25 8 7 9 10 25 35 20 25 20 10 10 20 20 25 35 25 25 25 c 20 1 10 FF B F BF FF-BF 20 FF-NBF F 2 4 6 3 5 7 8 9 10 F115W + F150W + F200W Fig. 1 | The Firefly Sparkle is a redshift zspec = 8.296 ± 0.001 gravitationally lines of lensing magnifications (μ = 15, 20, 30 and 40). b, RGB (F444W, F277W magnified arc lensed by the MACS J1423.8 + 2404 cluster. a, Full field with and F115W) image of the Firefly Sparkle showing different colours of the star the three objects of interest—Firefly Sparkle (centre), FF-BF (bottom left) and clusters. c, Combined short wavelength (F115W + F150W + F200W) image of FF-NBF (bottom right)—shown in boxes and circles. The contours show the the Firefly Sparkle, in which the distinct clusters can be seen. Scale bars, 1″. a Gaussian ellipse (reff = 0.01). The diffuse arc is also fit with a Gaussian suggests a hotter source than typical massive type O stars. This implies ellipse. All 11 components are simultaneously fit in each filter to derive a higher upper mass limit for the IMF or a top-heavy IMF12. the total flux (see section ‘Photometry of Firefly Sparkle’). We derived We derive the physical properties of the arc and the star clusters an upper limit on the half-light radii (Reff < 0.02), which is 0.5 times by performing SED fitting using various models, including simple the FWHM of the PSF in the F115W image, for all 10 clusters, including stellar population (SSP) models and non-parametric star formation FF-4, whose deconvolved size is smaller than the PSF of the image. histories using the Dense Basis method (see Methods for description). As the tangential magnifications of the clusters range from μtan = 12 SSP models are typically used in studies of star clusters in the local to μtan = 21, this results in an upper limit on the half-light sizes of less universe, as both observational and numerical works find that they than 4–7 pc. can be approximated to single bursts, whereas other models allow for We derive properties of the ionization sources in Firefly Sparkle using extended star formation histories. Where spectrophotometry is avail- NIRSpec Prism spectra from two adjacent shutters (slit 1 and slit 2) able, we include NIRSpec Prism spectra along with NIRISS and NIRCam covering the central brightest region of the Firefly Sparkle. Spectra photometry in the two slits (Extended Data Fig. 4), followed by photo- and properties from slit 1 are shown in Fig. 2, whereas spectra from metric modelling of the properties of the individual clusters. Extended all slits (including BF and NBF) are shown in Extended Data Fig. 2. The Data Table 2 shows the properties of the 10 clusters from the four spectrum of slit 1 (including light from FF-6 and contributions from different fits. FF-5 and the diffuse arc) shows a Balmer jump at λobs ≈ 3.5 μm (ref. 11) The demagnified stellar masses of the 10 clusters are about 105–106M⊙ and a smooth turnover at λobs ≲ 1.4 μm, possibly because of two-photon when fit with SSPs, similar to those of globular clusters. The surface continuum12. The absence of this feature in the spectrum at slit 2 (see density of the star clusters ranges between 103 and 104M⊙ pc−2, similar section ‘Photometry of Firefly Sparkle’) makes damping wing of neutral to Milky Way globular clusters13 (see Fig. 3b adapted from ref. 14). Dense hydrogen absorption unlikely. We infer significant nebular continuum Basis fits, accounting for extended star formation histories, indicate contribution to the overall SED for FF-6 and possibly other clusters with higher masses and specific star formation rates (sSFR) of about 10−7 yr−1, similar high-EW line emission. The ionizing source effective tempera- showing a sharp rise in the past 10–100 Myr. The precise nature of the ture of Teff = 105.1 K, obtained by modelling the nebular continuum with clusters depends on the interpretation of their star formation histories. CLOUDY and from line ratios, after confirming the Balmer decrement If seen as star clusters, their masses lead to crossing times of 1–4 Myr. is dust-free (see section ‘Spectroscopy extraction and spectral fitting’), Combined with their age estimates, this puts them at tage/tcross ~ 1–2, Nature | Vol 636 | 12 December 2024 | 333 Article a b F115W + F277W + F444W Slit 2 Slit 1 0.3 Slit 1 [NeIII] [NeIII] [OIII] [OIII] [OIII] [OIII] [OII] Lyα Hβ Hδ Hγ 0.2 Slit 2 fν (μJy) 0.1 Slit 1 0 0.7 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 Observed wavelength (μm) 2 c Likelihood (∝e–͜ ) d e 10 −6 10−4 10−2 Firefly Sparkle (slit 1) Top heavy (SSP) 4.8 3 Firefly Sparkle (slit 2) 0.8 (0.82 < D< 1.88) CEERS at z = 2–9 Topping + 24 Top heavy (DB) log([OIII]5007/[OII]3727) SDSS 0.6 (0.94 < D< 2.25) 2 Likelihood log[Telectron (K)] Canonical 4.4 Kroupa 0.4 (D= 2.3) 1 0.2 0 Top light (D > 3) 4.0 0 4.6 5.0 5.4 0 1 2 3 1 2 3 4 log[Teff (K)] log([OIII]4959+5007/[OIII]4363) IMF high-mass slope D Fig. 2 | Physical properties of the Firefly Sparkle. The physical properties continuum (1σ, 3σ and 5σ shown by contours). The slit 1 region exhibits electron measured from the NIRSpec and Prism spectra, which include light from FF-6 temperature of Telectron ~ 20,000 K, and ionizing source effective temperature of as well as fractional contributions from FF-3, FF-4 and FF-5 (Fig. 1). a, Positions Teffective ~ 105 K. d, Emission line diagnostics estimated from the fitted line ratios of the slitlets on the arc. b, The 2D NIRSpec spectrum for slit 1 (middle line, for RO3 and O32 hints at a metal-poor stellar population. Error bars show 1σ containing light from FF-6 with minor contributions from FF-5 and the diffuse uncertainties on the line ratios as derived in section ‘Spectral fitting in Firefly arc) and slit 2 (top line, containing light from FF-4 with minor contributions Sparkle slit 1’. e, Marginal likelihood for the high-mass IMF slope from joint from FF-3 and the diffuse arc) and the 1D spectrum of slit 1 only (see Extended spectrophotometric fitting with DENSE BASIS and SSP fits suggest a top-heavy Data Fig. 2 for slit 2 spectrum). c, Likelihood of the effective black-body IMF (α < 2), full posteriors shown in Extended Data Fig. 4. temperature and electron temperature from CLOUDY modelling of the nebular which indicates that they are marginally bound. By contrast, if they are photoionization data, affecting estimates of stellar masses, IMF and star nuclear star clusters or the remnants of dwarf galaxies that have previ- formation histories (SFH). Future JWST observations will provide bet- ously merged with the system, their ages are consistent with having ter constraints. In the meantime, we have mitigated the impact on our survived several crossing times, and they are likely to remain bound interpretation by using four independent SED models and focusing on until ejected from the system or integrated into the nucleus. The aspects that are common to the models, and by providing independent smooth component of the arc has more demagnified mass than any measurements where possible (for example, electron temperature). individual cluster, at log(M⋆ /M⊙) = 6.7 +0.9 −0.8 and an sSFR similar to the Improvements to population synthesis models and refinements in the star clusters (Fig. 4). The total demagnified mass of Firefly Sparkle is lens magnification model would help, although the estimated ages, log(M⋆ /M⊙) = 7.0+1.0 −0.3 , one of the lowest stellar mass objects observed sSFR and surface densities would be mostly unaffected. at this epoch, similar in stellar mass to the progenitor of a Milky Way Irrespective of these limitations, the Firefly Sparkle provides insights mass galaxy at z ~ 8 (Fig. 3). into the early galaxy formation. With massive star clusters exhibiting The slit 1 region shows extremely low metallicity (log(Z /Zgas) = high surface density, low metallicity, high electron temperature and +0.13 +0.22 −0.56−0.27 , 12 + log(O /H ) = 7.05−0.37 ) , among the lowest observed at hints of a top-heavy IMF, the Firefly Sparkle exhibits the hallmarks z > 6 (refs. 9,15,16). Our analysis of slit 1 using varying power-law slopes of star formation in extreme environments, consistent with sce- of the Kroupa IMF in FSPS indicates an excess of high-mass stars. Both narios such as pressure-regulated feedback dominated star forma- SSP and Dense Basis fits show a preference for top-heavy IMFs in the tion17–19, although further observations of the gas mass are needed MILES + MIST fits (Fig. 2). The high-mass star excess results in a domi- to ascertain this. The stellar mass of the galaxy is consistent with nant nebular continuum and high equivalent width emission lines12. progenitors of Milky-Way-like galaxies, derived using the abundance The fits rule out top-light IMFs (α = 2.3) and prefer top-heavy slopes matching method20,21 and TNG50 simulation22. The Firefly Sparkle +0.9 (αslit1 = 1.7 −0.7 ), consistent with a high ionizing source effective tem- suggests that early galaxy assembly can occur by dense star clusters perature of more than 40,000 K. as well23–25. The current analysis has several limitations. The spectrophotometric The Firefly Sparkle is the farthest spectroscopically confirmed galaxy models are influenced by star formation history, stellar population and with well-resolved star clusters, made visible by gravitational lensing 334 | Nature | Vol 636 | 12 December 2024  a b 1011 2.5 5.0 7.5 106 Age (Myr) 1010 105 Surface density (Mത pc–2) 109 104 M (Mത) 103 108 102 Firefly Sparkle 10 7 Ref. 40 Firefly Sparkle (zspec = 8.3) Ref. 39 101 Cosmic Gems (zphot = 10.1) Ref. 37 Sunrise (zphot = 6.2) 106 Ref. 49 Ref. 38 Sunburst (zspec = 2.4) 100 TNG50 (ref. 22) Milky Way globular clusters Milky Way (ref. 48) Young star clusters in nearby galaxies 105 0 2 4 6 8 100 101 z Radius (pc) Fig. 3 | The Firefly Sparkle in context. a, Progenitors of Milky Way galaxy star clusters in the Firefly Sparkle and their sizes are shown, colour-coded by analogues in the TNG50 simulation22 and from observations through abundance their ages from the SSP-MIST + MILES model. Nine of the ten star clusters are matching 20 applied to observed stellar mass functions19,35–37. The orange star unresolved; hence an upper limit on the demagnified sizes (HWHM/μtan) and a and error bars show the median and 1σ uncertainties on the stellar mass derived lower limit on the surface density are shown. Other high-redshift star clusters from SED fitting. Given its stellar mass, the Firefly Sparkle is within 1σ of the in magnified galaxies (Sunburst28, Sunrise29,41 and Cosmic Gems14) are shown median mass of a progenitor at z ~ 8.3. For comparison, stellar masses and for comparison. Milky Way globular clusters42 and young star clusters in redshifts of galaxies from previous JWST observations are shown (filled star-forming spiral galaxies in the Local Volume (distance less than 16 Mpc)43 circles for spectroscopy and unfilled circles for photometry23,37–40). The Firefly are also shown. The Firefly Sparkle star clusters are denser than local star Sparkle is one of the lowest stellar mass systems observed and the only with clusters but are comparable to the higher redshift analogues. spectroscopically confirmed star clusters at z > 8. b, The surface density of the a b Firefly Sparkle 0.6 and JWST sensitivity. JWST observations, combined with those of other SFR (Mത yr –1) NBF Firefly Sparkle distant galaxies12,14,26–29, open a new area of study into the role of massive 0.4 BF star clusters in early galaxy formation. These sites of dense and rapid 0.2 star formation in distant galaxies have an uncertain future. They may BF NBF 0 0 0.2 0.4 0.6 survive as present-day globular clusters30 or be stripped by tidal forces Time (Gyr) in the nascent disk of the galaxy to become nuclear star clusters31–33. c d Some clumps may even survive tidal stripping and loss, as simulations suggest that they can reaccrete gas in the turbulent environment34. 105 106 Future observations by JWST and ALMA will help distinguish these log(M /M ) pc–2 BF possibilities. 1 1 2 2 Online content Any methods, additional references, Nature Portfolio reporting summa- ries, source data, extended data, supplementary information, acknowl- 6 4 6 4 10 8 10 8 edgements, peer review information; details of author contributions 5 5 3 and competing interests; and statements of data and code availability 3 9 9 7 are available at https://doi.org/10.1038/s41586-024-08293-0. 7 1. Willott, C. J. et al. The Near-infrared Imager and Slitless Spectrograph for the James Webb Space Telescope. II. Wide field slitless spectroscopy. Publ. Astron. Soc. Pac. 134, 025002 (2022). Fig. 4 | The smooth component of the arc. a, Image of the Firefly Sparkle, FF-BF 2. Postman, M. et al. The Cluster Lensing and Supernova Survey with Hubble: an overview. and FF-NBF. b, Non-parametric integrated star-formation histories as a function Astrophys. J. Suppl. Ser. 199, 25 (2012). 3. Hoag, A. et al. Spectroscopic confirmation of an ultra-faint galaxy at the epoch of of lookback time of these three galaxies reconstructed from SED fitting. The reionization. Nat. Astron. 1, 0091 (2017). Firefly Sparkle and FF-BF both show a recent burst of star formation in the past 4. Gardner, J. P. et al. The James Webb Space Telescope Mission. Publ. Astron. Soc. Pac. 135, approximately 50 Myr indicative of recent interactions. c, The source plane RGB 068001 (2023). (F444W, F277W and F115W) image of the Firefly Sparkle. d, The source plane 5. Doyon, R. et al. The Near Infrared Imager and Slitless Spectrograph for the James Webb Space Telescope. I. Instrument overview and in-flight performance. Publ. Astron. Soc. Pac. stellar mass surface density map of the Firefly Sparkle reconstructed from the 135, 098001 (2023). individual SED fits of the components and its source plane model; 49–57% of 6. Rieke, M. J. et al. Performance of NIRCam on JWST in flight. Publ. Astron. Soc. Pac. 135, the total stellar mass is concentrated in the star clusters. Scale bars, 1 kpc. 028001 (2023). Nature | Vol 636 | 12 December 2024 | 335 Article 7. Jakobsen, P. et al. The Near-Infrared Spectrograph (NIRSpec) on the James Webb Space 29. Vanzella, E. et al. Early results from GLASS-JWST. VII. Evidence for lensed, gravitationally Telescope. I. Overview of the instrument and its capabilities. Astron. Astrophys. 661, A80 bound protoglobular clusters at z = 4 in the Hubble Frontier Field A2744. Astrophys. J. Lett. (2022). 940, L53 (2022). 8. Kneib, J. P., Ellis, R. S., Santos, M. R. & Richard, J. A probable z ~ 7 galaxy strongly lensed by 30. Belokurov, V. & Kravtsov, A. Nitrogen enrichment and clustered star formation at the dawn the rich cluster A2218: exploring the Dark Ages. Astrophys. J. 607, 697–703 (2004). of the Galaxy. Mon. Not. R. Astron. Soc. 525, 4456–4473 (2023). 9. Langeroodi, D. et al. Evolution of the mass–metallicity relation from redshift z ≈ 8 to the 31. Bournaud F. in Galactic Bulges (eds Laurikainen, E. et al.) 355–390 (Springer, 2016). local universe. Astrophys. J. 957, 39 (2023). 32. Mandelker, N. et al. Giant clumps in simulated high-z galaxies: properties, evolution and 10. Peng, C. Y., Ho, L. C., Impey, C. D. & Rix, H.-W. Detailed decomposition of galaxy images. dependence on feedback. Mon. Not. R. Astron. Soc. 464, 635–665 (2017). II. Beyond axisymmetric models. Astron. J. 139, 2097–2129 (2010). 33. Tamfal, T. et al. The dawn of disk formation in a Milky Way-sized galaxy halo: thin stellar 11. Wilkins, S. M. et al. First Light and Reionization Epoch Simulations (flares) – XIV. The disks at z > 4. Astrophys. J. 928, 106 (2022). Balmer/4000 Å breaks of distant galaxies. Mon. Not. R. Astron. Soc. 527, 7965–7973 (2024). 34. Bournaud, F. et al. The long lives of giant clumps and the birth of outflows in gas-rich 12. Claeyssens, A. et al. Star formation at the smallest scales: a JWST study of the clump galaxies at high redshift. Astrophys. J. 780, 57 (2014). populations in SMACS0723. Mon. Not. R. Astron. Soc. 520, 2180–2203 (2023). 35. Grazian, A. et al. The galaxy stellar mass function at 3.5 ≤ z ≤ 7.5 in the CANDELS/UDS, 13. Harris, W. E. & Racine, R. Globular clusters in galaxies. Annu. Rev. Astron. Astrophys. 17, GOODS-South, and HUDF fields. Astron. Astrophys. 575, A96 (2015). 241–274 (1979). 36. Stefanon, M. et al. The brightest z ≳ 8 galaxies over the COSMOS UltraVISTA field. 14. Adamo, A. et al Bound star clusters observed in a lensed galaxy 460 Myr after the Big Astrophys. J. 883, 99 (2019). Bang. Nature 632, 513–516 (2024). 37. Tang, M. et al. JWST/NIRSpec spectroscopy of z = 7–9 star-forming galaxies with CEERS: 15. Hirschmann, M., Charlot, S. & Somerville, R. S. High-redshift metallicity calibrations for new insight into bright Lyα emitters in ionized bubbles. Mon. Not. R. Astron. Soc. 526, JWST spectra: insights from line emission in cosmological simulations. Mon. Not. R. 1657–1686 (2023). Astron. Soc. 526, 3504–3518 (2023). 38. Austin, D. et al. A large population of faint 8 < z < 16 galaxies found in the first JWST NIRCam 16. Trump, J. R. et al. The physical conditions of emission-line galaxies at cosmic dawn from observations of the NGDEEP survey. Astrophys. J. Lett. 952, L7 (2023). JWST/NIRSpec spectroscopy in the SMACS 0723 early release observations. Astrophys. J. 39. Nakajima, K. et al. JWST census for the mass-metallicity star formation relations at z = 4–10 945, 35 (2023). with self-consistent flux calibration and proper metallicity calibrators. Astrophys. J. 17. Ellison, S. L. et al. The ALMaQUEST Survey XI: a strong but non-linear relationship between Suppl. Ser. 269, 33 (2023). star formation and dynamical equilibrium pressure. Mon. Not. R. Astrophys. Soc. 527, 40. Tacchella, S. et al. On the stellar populations of galaxies at z = 9–11: the growth of metals 10201–10220 (2024). and stellar mass at early times. Astrophys. J. 927, 170 (2022). 18. Kim, C.-G. et al. Metallicity dependence of pressure-regulated feedback-modulated star 41. Vanzella, E. et al. JWST/NIRCam probes young star clusters in the reionization era sunrise formation in the TIGRESS-NCR simulation suite. Astrophys. J. 972, 67 (2024). arc. Astrophys. J. 945, 53 (2023). 19. Ostriker, E. C. & Kim, C.-G. Pressure-regulated, feedback-modulated star formation in disk 42. Gieles, M. & Portegies Zwart, S. F. The distinction between star clusters and associations. galaxies. Astrophys. J. 936, 137 (2022). Mon. Not. R. Astron. Soc. Lett. 410, L6–L7 (2011). 20. Behroozi, P. S. et al. Using cumulative number densities to compare galaxies across 43. Brown, G. & Gnedin, O. Y. Radii of young star clusters in nearby galaxies. Mon. Not. R. cosmic time. Astrophys. J. Lett. 777, L10 (2013). Astron. Soc. 508, 5935–5953 (2021). 21. Papovich, C. et al. ZFOURGE/CANDELS: on the evolution of M* galaxy progenitors from z = 3 to 0.5. Astrophys. J. 803, 26 (2015). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 22. Pillepich, A. et al. Milky Way and Andromeda analogues from the TNG50 simulation. Mon. published maps and institutional affiliations. Not. R. Astron. Soc. 535, 1721–1762 (2024). 23. Feng, Y. et al. The formation of Milky Way-mass disk galaxies in the first 500 million years of a cold dark matter universe. Astrophys. J. Lett. 808, L17 (2015). Open Access This article is licensed under a Creative Commons Attribution 24. Iza, F. G., Nuza, S. E. & Scannapieco, C. The temporal evolution of gas accretion onto the 4.0 International License, which permits use, sharing, adaptation, distribution discs of simulated Milky Way-mass galaxies. Bol. Asoc. Argent. Astron. La Plata Argent. and reproduction in any medium or format, as long as you give appropriate 64, 143–145 (2023). credit to the original author(s) and the source, provide a link to the Creative Commons licence, 25. Nepal, S. et al. Discovery of the local counterpart of disc galaxies at z > 4: the oldest thin and indicate if changes were made. The images or other third party material in this article are disc of the Milky Way using Gaia-RVS. Astron. Astrophys. 688, A167 (2024). included in the article’s Creative Commons licence, unless indicated otherwise in a credit line 26. Fujimoto, S. et al. Primordial rotating disk composed of ≥15 dense star-forming clumps at to the material. If material is not included in the article’s Creative Commons licence and your cosmic dawn. Preprint at https://arxiv.org/abs/2402.18543 (2024). intended use is not permitted by statutory regulation or exceeds the permitted use, you will 27. Mowla, L. et al. The Sparkler: evolved high-redshift globular cluster candidates captured need to obtain permission directly from the copyright holder. To view a copy of this licence, by JWST. Astrophys. J. Lett. 937, L35 (2022). visit http://creativecommons.org/licenses/by/4.0/. 28. Rivera-Thorsen, T. E. et al. The Sunburst Arc: direct Lyman α escape observed in the brightest known lensed galaxy. Astron. Astrophys. 608, L4 (2017). © The Author(s) 2024 336 | Nature | Vol 636 | 12 December 2024 Methods the PSF centre. We also fix the bending mode B2 (2.14), ellipse radius (3.9″), axis ratio (0.08) and position angle (−51.8°) of the arc from the Image preparation F444W fit. We fix the morphology of FF4 also with radius = 0.59″, axis The cluster field MACS J1423.8 + 2404 was observed with JWST/NIRCam ratio = 0.1 and position angle = −53°. imaging using filters F090W, F115W, F150W, F200W, F277W, F356W, We now fit all 11 components in all 10 filters to determine their fluxes. F410M and F444W with exposure times of 6.4 ks each, reaching a The resulting models and residuals are shown in Extended Data Fig. 1. signal-to-noise ratio between 5 and 10 for an mAB = 29 point source. Residuals from the fits are negligible, as shown by χ2/ν ~ 1 in the GALFIT It was also observed with JWST/NIRISS imaging using filters F115W, fits in all filters. This confirms the original visual impression that nine F150W and F200W. of the ten clusters are unresolved and an additional smooth compo- To reduce the imaging data, we use the photometric pipeline that is nent is present. presented in more detail in ref. 44. Briefly, the raw data has been reduced To derive the uncertainty in our flux estimation, we inject the full using the public grism redshift and line analysis software Grizli43, Firefly Sparkle model in 100 random locations in our 10″ × 10″ post- which masks imaging artefacts, provides astrometric calibrations age stamps (avoiding the edge) and refit with the exact same setting based on the Gaia Data Release 3 catalogue13 and shifts images using of GALFIT. We find no significant systematic offset between the fitted Astrodrizzle. The photometric zero-points are applied as described in flux and the injected flux for any of the 11 components, in any of the fil- ref. 34. RGB image created using six filters of NIRCam observation of ters, showing that our photometric technique is robust to background the Firefly Sparkle is shown in Fig. 1. We used images from which bright variations across all filters. The uncertainty in the photometry is calcu- cluster galaxies and intracluster light have been removed, as described lated from the bi-weight scale of the 100 refitted fluxes. The resulting in ref. 25. The methodology for modelling and removing diffuse light photometry and the RGB image of the model and the residual are shown from cluster galaxies and intracluster light (ICL) is presented in ref. 25. in Extended Data Fig. 1. The agreement between NIRISS and NIRCam The NIRCam depths (0.3′ diameter aperture) for F090W, F115W, F150W, fluxes in the three overlapping filters is another confirmation of the F200W, F277W, F356W, F410M and F444W are 7.2, 6.6, 5.2, 4.4, 3.0, 2.9, robustness of photometry. We have used updated zero-points34 and cor- 5.5 and 4.3 nJy, respectively, and the NIRISS depths for F115WN, F150WN rected for Milky Way extinction using the colour excess E(B − V) = 0.0272 and F200WN are 3.6, 4.3 and 4.0 nJy, respectively41. from ref. 6 and assuming the extinction law in ref. 35 using the factor between the extinction coefficient and colour excess RV = 3.1. Photometry of Firefly Sparkle We perform photometry in 10 JWST bands (NIRISS: F115WN, F150WN Spectroscopy extraction and spectral fitting and F200WN; NIRCam: F115W, F150W, F200W, F277W, F356W, F410M NIRSpec spectroscopy has been acquired for MACS J1423.8 + 2404 and F444W) in which the Firefly Sparkle is detected from their mor- and spectra were obtained for the Firefly Sparkle, FF-BF and FF-NBF. phological fit with GALFIT. In other JWST and HST filters, the Firefly The spectra for the FF-BF were part of the sample in ref. 23, with Sparkle is not or barely detected; hence, we place upper limits for the zspec = 8.2953 ± 0.0005. The spectra were observed using the PRISM/ entire source. As the object is resolved into at least 10 distinct clusters CLEAR disperser and filter, through three Micro-Shutter Assembly and a diffuse galaxy component, we perform a morphological fit using (MSA) masks per cluster with a total exposure time of 2.9 ks per MSA Galfit10 to extract the photometric information. configuration. Point spread functions are extracted empirically by median stack- The NIRSpec data were processed using the STScI JWST pipeline ing bright, isolated, non-saturated stars following the methodology (software v.1.8.4 and jwst_1030.pmap) and the msaexp package31. described in ref. 28. Convolution kernels for homogenizing all data We used the standard JWST pipeline for the level 1 processing, in to the F444W resolution are created with photutils.psf.matching which we obtained the rate fits files from the raw data. We enabled using a SplitCosineBellWindow() windowing function to remove the jump step option expand_large_events to mitigate contamina- high-frequency noise, which results from floating-point imprecision tion by snowball residuals and used a custom persistence correction when taking the ratio of Fourier transforms. We optimize the shape that masked out pixels that approach saturation within the following of each window function to minimize the median residual between 1,200 s for any readout groups. We then used msaexp for level 2 pro- convolved stars from each source filter that is convolved and stars cessing, for which we performed the standard wavelength calibration, from the target F444W filter. flat-fielding, path-loss correction and photometric calibration and For the morphological fit, we create 10″ × 10″ postage stamps in all obtained the 2D spectrum before background subtraction. As the 10 filters from the BCG-subtracted images. We determine the priors central and upper shutters contain different clusters (see Fig. 2a to for the centres of the 10 clusters by visual inspection. Although nine find the shutter positions), we need custom background subtraction out of the ten appear as point sources, FF-4 has an elongated shape to avoid self-subtraction. We did this by building the background 2D and appears unresolved. We first determine the central coordinates of spectrum by stacking and smoothing the sky spectrum in the empty the 10 clusters and the arc by fitting (1) an elliptical Gaussian for FF-4; pixels and obtained the background subtracted 2D spectrum of Fire- (2) nine point sources for the other nine clusters; and (3) another ellip- fly Sparkle. We confirmed that this custom background subtraction tical Gaussian with the bending mode turned on for the diffuse arc to method works as well as a standard drizzle background subtraction the F115W image, which has the highest resolution (smallest PSF). The method used in the literature33, using a well-isolated galaxy spectrum free parameters are the centres and total fluxes of all the components, from the CANUCS observation (Asada et al., in prep.). We finally extract the radius and axis ratio of FF-4, and the radius, axis ratio and bend- the 1D spectrum separately in slit 1 and slit 2, by collapsing the 2D ing mode (B2) of the arc. The initial guesses for the coordinates were spectrum using an inverse-variance weighted kernel following the determined by visual inspection of the F115W image. Once we obtain the prescription in ref. 24. We verified that the uncertainty array of the 1D fitted central coordinates of all the components from F115W, we again spectrum has the appropriate normalization by testing the distribu- fit all 11 components in F444W, which has the highest signal-to-noise tion of spectral fluctuations in an empty sky region and finding the ratio for the arc and FF-4, to determine the radius, axis ratio, position fractions of pixels at >1 and >2σ as expected. angles of the ellipses, and the bending mode B2 of the arc. We use the best-fit centre coordinates from F115W as the central Spectral fitting in Firefly Sparkle slit 1. The resulting 1D spectrum of coordinates in all the filters. However, instead of fixing the central Firefly Sparkle in slit 1, dominated by the cluster FF-6, is shown in Fig. 2. coordinates, we allow GALFIT to fit for them in every filter within a very The spectrum exhibits a Balmer jump at λobs ~ 3.5 μm and a turnover narrow range of ±0.5 pixels (0.02″) to account for the uncertainty in at λobs ≲ 1.4 μm, probably because of two-photon emission. These Article features suggest that the nebular continuum should dominate over [Oiii]1661+1666. We assume the electron density to be ne = 103 cm−3, which the stellar continuum in the rest frame UV to optical spectrum within is consistent with recent JWST observations of similarly high-z galaxies7 slit 1 (as found for a z = 5.9 galaxy in ref. 12). We thus model the con- and obtain consistent independent temperature measurements within tinuum of the spectrum with nebular continuum using the photoion- the uncertainties (Te,O++ = 4.0−0.9+2.6 K and 2.9+0.7 4 −0.4 × 10 K , respectively; ization code CLOUDY v.23 (ref. 5). To determine the dust attenuation Extended Data Fig. 3 (right)). Note that because the [Oiii]λλ1661 + 1666 value in the continuum model fitting, we first measure the Hγ/Hβ ratio detection is tentative and potentially blended with Heiiλ1640, we con- by fitting the Gaussian profiles. The ratio agrees well with the case B sider [Oiii]λ4363 to be more reliable. recombination, and no significant dust attenuation is indicated. There- We note that in ref. 16, the authors measured a similar ratio of fore, in the continuum spectral modelling, we use pure hydrogen gas [Oiii]4959+5007/[Oiii]4363 in the z = 6 galaxy RXCJ2248-ID to that of slit 1. irradiated by an ionizing source having black-body SED without dust In ref. 16, medium resolution spectroscopy was used to determine the attenuation. We vary the effective temperature of the black body (Teff ) electron density directly. They found that when using lines with higher and the electron temperature of the (ionized) hydrogen gas (Te,H+), ionization potential than O+, the electron density was higher and search for the best-fitting model continuum by χ2 minimization. (ne ~ 105 cm−3) than is typically found from [Oii]λ3727 (ref. 7). This high In the continuum fitting, we mask out emission line regions and electron density leads to a lower electron temperature for their galaxy all wavelengths λobs < 1.2 μm at which the Lyman break is seen in the of Te,O++ = 2.5 × 104 K. Similarly, if we assume the electron density of slit 2 spectrum, because this region may be affected by a neutral ne = 105 cm−3 instead for our slit 1 spectrum, the electron temperature hydrogen damping wing. The best-fit model has log(Teff/K) = 5.10 and from [Oiii]λ4363 becomes Te,O++ = 3.2 +1.6 −0.96, which is in between the two log(Te,H+ /K ) = 4.34, which is fully consistent with the results in ref. 12. measurements based on [Oiii]λλ1661 + 1666 and [Oiii]λ4363 The result of continuum fitting does not change if we consider a slight when assuming ne ~ 103 cm−3 above. To consider the possibility of a dust attenuation (AV = 0.1 mag) in the fitting. As discussed in ref. 12, the somewhat higher electron density in the highly ionized region, we effective temperature of log(Teff/K) = 5.10 is much hotter than typical adopt the mean value of our two electron temperature measurements massive type O stars and is suggestive of this star-forming cluster (Te,O++ = 3.5 × 104 K) as our fiducial value and propagate the full range having a top-heavy IMF. The IMF of this cluster is further discussed in of the two measurement uncertainties into the following metallicity section ‘SED fitting analysis’. measurement. Note that the UV continuum turnover feature could be because of Based on the electron temperature measurement, we obtained the the absorption from dense neutral hydrogen either in the intergalactic oxygen abundance from [Oiii]4959+5007/Hβ and [Oii]3727/Hβ ratios, fol- medium (IGM) or in the circumgalactic medium (CGM). However, in the lowing the prescription in ref. 8. We assume the electron density to be case of slit 1 spectrum, we expect the effect of IGM and CGM damping ne = 103 cm−3. The total oxygen abundance is calculated from O++/H+ absorption to be negligible or limited at λobs < 1.2 μm based on the blue and O+/H+, and the higher ionizing state oxygen is ignored30. As the continuum and sharp drop-out in the slit 2 spectrum (see section ‘Spec- [Oii]λ3727 emission line is undetected, we can obtain only an upper tral fitting in Firefly Sparkle slit 2’ for details of slit 2 spectrum). Consid- limit for O+/H+, but the upper limit for the abundance of the singly ering the spatial proximity of the slit 1 and slit 2 regions (Fig. 2), we can ionized oxygen is negligibly small as compared with the doubly ionized assume the absorption feature from line-of-sight neutral hydrogen to oxygen. We thus derived the total oxygen abundance from O++/H+, +0.22 be the same in the slit 1 and slit 2 spectra. The slit 2 spectrum is rather yielding 12 + log(O/H) = 7.05−0.37 ( Zgas /Z⊙ = 0.02 +0.04 −0.01 assuming the blue and has a sharp Lyman break starting at λobs = 1.2 μm, whereas solar abundance to be 8.69; ref. 38). the slit 1 spectrum shows the turnover starting at λobs ~ 1.4 μm. Thus, We also derive the ionization parameters using the ionization- the turnover feature should not be because of the neutral hydrogen sensitive emission line ratios: [Oiii]5007/[Oii]3727 and [Neiii]3869/[Oii]3727. absorption, but rather because of the intrinsic continuum shape of the Following the prescription in refs. 45,46, we obtain the lower limit for source. Nevertheless, to avoid the possible effect of the neutral hydro- the ionization parameters (log U) from these two ratios. Both ratios gen absorption, we mask out λobs < 1.2 μm in the nebular continuum provide a similar limit of log U > −2.0. fitting above (corresponding to 1,290 Å in the rest frame). All the emission line flux measurements and the derived physical Having the model continuum, we subtract the underlying model parameters in Firefly Sparkle slit 1 are presented in Extended Data continuum from the observed spectrum and measure the spectro- Table 1. We also compare the diagnostic emission line ratios in Fire- scopic redshift and emission line fluxes by fitting Gaussian profiles. fly Sparkle with those in other galaxy population in Fig. 2d. We use The best-fitting model spectrum with nebular continuum and Gauss- the ionization-sensitive line ratio O32 ([Oiii]5007/[Oii]3727) and the ian profiles is shown in Fig. 2b (red solid curve). We securely detect temperature-sensitive line ratio RO3 ([Oiii]4959+5007/[Oiii]4363) and com- emission lines of [O iii]λλ4959, 5007, Hβ, [Oiii]λ4363, Hγ, Hδ and pare these line ratios with other [Oiii]λ4363-detected galaxies at z = 2–9 [Neiii]λλ3869, 3889. We do not find significant detection of [Oii]λ3727 from previous JWST observations2 and those in the local universe from and obtain an upper limit for the flux of this line. There is a tentative SDSS observations14. Extended Data Fig. 3 (middle) presents a similar detection of the blended line of [Oiii]λλ1661 + 1666, although the comparison but uses another ionization-sensitive line ratio Ne3O2 spectral resolution of the prism is low at this wavelength making this ([Neiii]3869/[Oii]3727) instead of O32. doublet difficult to securely detect and separate from Heiiλ1640. We use these emission line fluxes to estimate the physical parameters in Spectral fitting in Firefly Sparkle slit 2. In contrast to slit 1, the slit 1. We first estimate the dust attenuation based on Balmer decre- extracted 1D spectrum in Firefly Sparkle slit 2 does not show nebular ments. Both the Hγ/Hβ and Hδ/Hβ ratios are consistent with theo- continuum features, and the blue continuum is rather smooth with a retical predictions in case B recombination21 within the uncertainties, sharp drop-out because of the Lyman break at λobs ~ 1.2 μm. We thus suggesting there is no significant dust attenuation (Extended Data derive the emission line fluxes from the slit 2 spectrum by fitting Gauss- Fig. 3, red squares in the left). This result is consistent with the initial ian profiles with the continuum being modelled by a constant offset measurement before the continuum fitting above and supports the around each emission line. We detect [Oiii]λλ4959,5007, Hβ, Hγ, Hδ, validity of the dust-free assumption in the nebular continuum fitting [Neiii]λ3869 and [Oii]λ3727 emission lines in the slit 2 spectrum but process. Therefore, we do not correct for dust attenuation in the fol- do not detect [Oiii]λ4363. lowing measurements of emission line ratios and physical parameters We then derive the physical properties in the same way as done for in this section. Firefly Sparkle slit 1 spectrum. We measure the dust attenuation from We next measure the electron temperature using temperature- Balmer decrement, Hγ/Hβ and Hδ/Hβ, and find both line ratios agree sensitive emission line ratios: [Oiii]4959+5007/[Oiii]4363 and [Oiii]5007/ well with the predicted ratios under case B recombination (blue squares in Extended Data Fig. 3 (left)). This suggests that the dust attenuation that we then use to fit the photometry. We find that the fits are consist- is negligible in the slit 2 spectrum as well, and we do not make a dust ent with negligible dust attenuation, consistent with our estimates attenuation correction. from measuring the Balmer decrement. We also find that our fits rule As we do not detect the temperature-sensitive emission lines of out the part of parameter space consistent with the canonical Chabrier- [Oiii]λ1666 or [Oiii]λ4363 in the slit 2 spectrum, we cannot measure like or Kroupa-like IMF (with the high-mass slope ≈ 2.3) in favour of the electron temperature and the metallicity from the direct- more top-heavy slopes of about 1.5+0.7 −0.6 for slit 1, which contains portions temperature method. We thus obtain only the upper limit for the elec- of clusters 3, 4, 5 and 6. We find weaker constraints from the spectrum tron temperature (Te,O++ ) from the non-detection of [Oiii]λ4363. for slit 2, which still skews towards top-heaviness but with large uncer- The electron temperature in Firefly Sparkle slit 2 is shown to be tainties of about 1.7 +0.9 −0.7. Finally, we find estimates of both stellar and Te,O++ < 1.8 × 104 K (1σ) or 10% of maximum flux. We then computed the best magnification value for all selected pixels and computed the mean and standard deviation values for these to find the magnification of Data availability the arc (μ = 24.4 ± 6.0). All data supporting the findings of this study are publicly available The source plane reconstruction is made using the best GALFIT on the CANUCS website at GitHub (https://niriss.github.io/). Further model to compute the source plane positions and magnification requests for data can be directed to the corresponding authors. for the 10 star clusters. We use Lenstool to generate a source plane image reconstruction of the diffuse light of the galaxy with a smooth PSF-deconvolved model of its light profile. We use GALFIT to add 10 44. Baumgardt, H. & Makino, J. Dynamical evolution of star clusters in tidal fields. Mon. Not. R. point sources convolved with the appropriate PSFs to the diffused Astron. Soc. 340, 227–246 (2003). source plane model at the source plane positions of the star clusters 45. Luridiana, V., Morisset, C. & Shaw, R. A. PyNeb: a new tool for analyzing emission lines. with the demagnified fluxes. This process is repeated to generate source I. Code description and validation of results. Astron. Astrophys. 573, A42 (2015). 46. McLeod, D. J. et al. The evolution of the galaxy stellar-mass function over the last 12 billion plane models in all filters. We also generate a mass map using the same years from a combination of ground-based and HST surveys. Mon. Not. R. Astron. Soc. prescription, replacing the demagnified fluxes with the demagnified 503, 4413–4435 (2021). masses. The resulting source plane RGB image and mass map are shown 47. Jullo, E. et al. A Bayesian approach to strong lensing modelling of galaxy clusters. New J. Phys. 9, 447 (2007). in Fig. 4c,d. 48. Cameron, A. J. et al. Nebular dominated galaxies: insights into the stellar initial mass function at high redshift. Mon. Not. R. Astron. Soc. 534, 523–543 (2024). Size and surface density of star clusters 49. Tan, V. Y. Y. et al. A measurement of the assembly of Milky Way analogs at redshifts 0.5 < z < 2 with resolved stellar mass and star formation rate profiles. Astrophys. J. 964, We now investigate the spatial properties of the star clusters. Nine out 177 (2024). of the ten star clusters are unresolved even in our highest resolution 50. Desprez, G. et al. ΛCDM not dead yet: massive high-z Balmer break galaxies are less common than previously reported. Mon. Not. R. Astron. Soc. 530, 2935–2952 (2024). F115W NIRCam image. FF-4 has a slightly elongated shape visually but has a best-fit major axis size (0.01) smaller than the smallest PSF, making the size estimate unreliable. Hence, we use the half-width half-max of Acknowledgements We thank N. Murray for the discussions about the formation of star the NIRCam F115W PSF (0.02) to set an upper limit on the size of all 10 clusters and their dynamic state, and E. Nelson, J. Antwi-Danso, G. Bryan and R. Somerville star clusters. To determine the upper limits of the sizes of unresolved for their discussions. We thank N. Gosavi for coining the name Firefly and O. Trottier for his help with Fig. 4. This research was funded by grant 18JWST-GTO1 from the Canadian Space sources, we use the tangential eigenvalue of magnification 1/∣λt∣, which Agency and the Natural Sciences and Engineering Research Council of Canada. L.M., K.I. and ranges between 14 and 24. This results in a size upper limit between 4 pc R.A. acknowledge the support of the Dunlap Institute for Astronomy and Astrophysics at the and 7 pc. The central star clusters have the highest magnification and University of Toronto. K.I. was supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51508 awarded by the Space Telescope Science Institute, which is operated by the the smallest upper limits, whereas the ones near the two ends of the arc Association of Universities for Research in Astronomy, for NASA, under contract NAS5-26555. have the lowest. We use the upper limit on sizes and the demagnified Y.A. was supported by a Research Fellowship for Young Scientists from the Japan Society of stellar masses to calculate the lower limit on stellar surface densities the Promotion of Science (JSPS). M.B. and G.R. acknowledge support from the ERC Grant FIRSTLIGHT and from the Slovenian National Research Agency ARRS through grants N1-0238, as shown in Fig. 3b. P1-0188 and the program HST-GO-16667, provided through a grant from the STScI under NASA contract NAS5-26555. This research used the Canadian Advanced Network For Abundance matching for MW and M31 progenitors Astronomy Research (CANFAR) operated in partnership by the Canadian Astronomy Data Centre and The Digital Research Alliance of Canada with support from the National Research To estimate the range of stellar masses of progenitors of both MW-mass Council of Canada, the Canadian Space Agency, CANARIE and the Canadian Foundation for and M31-mass galaxies at higher redshift, we adopt a semi-empirical Innovation. Author contributions L.M. conducted the photometry and size measurements, K.I. performed Additional information the SED fitting, Y.A. conducted the spectral analysis, G.D. conducted the lens modelling and Supplementary information The online version contains supplementary material available at V.Y.Y.T. performed the progenitor mass analysis. L.M., K.I., Y.A., G.D. and V.Y.Y.T. also made the https://doi.org/10.1038/s41586-024-08293-0. figures and wrote the paper. G.B., C.W. and V.S. were responsible for the image processing Correspondence and requests for materials should be addressed to Lamiya Mowla or pipeline and generating image mosaics. N.M. conducted the BCG subtraction and G.S. Kartheik Iyer. performed the PSF modelling. All authors made contributions to the paper and provided Peer review information Nature thanks the anonymous reviewers for their contribution to the assistance in data analysis and interpretation. peer review of this work. Peer reviewer reports are available. Competing interests The authors declare no competing interests. Reprints and permissions information is available at http://www.nature.com/reprints. Article Extended Data Fig. 1 | Morphological fit with GALFIT of the Firefly Sparklein sources in F115W. The full model consists of nine point sources, an elliptical 11 JWST/NIRCam and NIRISS filters for photometric extraction and size Gaussian for cluster (FF-4), and an elliptical Gaussian with a bending mode for determination. The images, their respective models, and residues for 11 filters the diffuse arc. The photometry is derived from the total model flux of the and the RGB image (R: F444W, G: F277W, B: 115W) are shown. The cutouts used 11 components. The error of the photometry is estimated by injecting the full for the fitting have size 10’ × 10’. In this figure we have zoomed in on the central model in random locations in the MACS 1423 field, and refitting them them 7’ × 7’. The FWHM of the point spread functions of the respective filters are with GALFIT. The upper limit on the size of all clusters is determined by the shown as black circles on the lower left corner. The Firefly Sparkleis completely HWHM of the F115W PSF (0.02’), as the deconvolved size of FF-4 is smaller than invisible in the bluest filter (F090W). Based on the reduced χ2 of the fits, nine the F115W PSF. out the ten clusters of the Firefly Sparkleare consistent with being point Extended Data Fig. 2 | NIRSpec Prism spectra for two slits of the Firefly emission lines and a Lyman break in all the spectra unambiguously determine Sparkle, along with those of the nearby FF-BF and FF-NBF companions. the redshifts of all the components. There is a slight oversubtraction of Slit 1 is dominated by FF-6 and contains contributions from FF-5 in Fig. 1, while background at λ > 4 μm for Firefly SparkleSlit 2 and BF due to their locations Slit 2 is dominated by FF-4 and contains contributions from FF-3, FF-5 and FF-6, close to the bar of the NIRSpec MSA shutter. Further analysis of these regions is with both slits getting marginal contributions from the diffuse arc. Strong left for follow-up observations. Article Extended Data Fig. 3 | Inferring properties from the NIRSpec Prism spectrum temperature measurements from [O iii] emission lines in Firefly SparkleSlit 1. of the Firefly Sparkle. Left: Balmer decrements of Hδ/Hβ and Hγ/Hβ in Slit 1 The dashed lines with shaded regions are the measured line ratios. The (red) and Slit 2 (blue) spectra. Black solid line denotes the line ratios under Case B solid lines denotes predicted line ratios as a function of different electron recombination. The line ratios indicates the dust attenuation is not significant temperatures from PyNeb 45. The two different emission line ratios in both spectra. Middle: Similar to bottom middle in Fig. 2, but Ne3O2 ratio is independently suggest a high electron temperature of Te,O++ ~ 40000 K. used instead as an indicator of the ionization parameter. Right: Electron Extended Data Fig. 4 | Inferring physical properties from the covariances, while the diagonal plots show the marginal posteriors for each spectrophotometric fits. The top panel shows the spectrophotometry parameter. In addition to the joint posteriors, the spectra and photometry for Slit 1 (black line and points) along with the best-fit spectrum from Dense posteriors generally agree, with the spectra better able to constrain Basis (orange line and points) and fits using simple stellar populations (green parameters like the gas-phase metallicity. The inset panels on the right line and points). An inset panel shows the region with Hβ + [OIII] where the show the stellar population posteriors from photometry alone (light blue), spectrum has much higher high fluxes. The corner plot shows the posteriors spectroscopy alone (red) and joint (black lines) along with 1σ uncertainties for each parameter being fit with only photometry (light blue contours), only using the MILES+MIST and BPASS templates, again finding that they agree spectroscopy (red contours) and both (black contours) using Dense Basis. between photometric and spectroscopic fits within uncertainties. The contours show the 1-σ and 2σ regions for each posterior, along with their Article Extended Data Fig. 5 | Top: The multiwavelength 11-component model for composite F115W + F277W + F444W images. Bottom: Photometry for the 10 the resolved structure in the Firefly Sparkle consisting of 10 clusters and clusters are shown along with fits using DENSE BASIS (orange) and simple stellar the diffuse arc shown in Extended Data Fig. 1. The three panels show the populations (green), along with estimated stellar masses from the SSP fits. observed image (left), the GALFIT model (middle) and the residuals (right) in Extended Data Fig. 6 | The SFHs of the ten star clusters from the DENSE BASIS-MILES+MIST fits. The solid lines and shaded regions indicate the median and 1σ uncertainties for the SFHs of individual clusters. Article Extended Data Table 1 | Emission line flux measurements and the physical properties from the spectrum of Firefly Sparkle Slit 1 Fluxes are in units of 10−20 erg s−1 cm−2. Extended Data Table 2 | Median stellar mass and age estimates from the various SED modeling configurations described in Section 4 in Methods Article Extended Data Table 3 | Magnification and upper limit on the size of the Firefly Sparklestar clusters, the diffuse arc, BF, and NBF Extended Data Table 4 | Photometry of individual star clusters and the diffuse arc of the Firefly Sparkle Fluxes are in units of nJy.