Formation of a low-mass galaxy from star clusters in a 600-million-year-old Universe

Image preparation

The cluster field MACS J1423.8 + 2404 was observed with JWST/NIRCam imaging using filters F090W, F115W, F150W, F200W, F277W, F356W, F410M and F444W with exposure times of 6.4 ks each, reaching a signal-to-noise ratio between 5 and 10 for an m_AB = 29 point source. It was also observed with JWST/NIRISS imaging using filters F115W, F150W and F200W.

To reduce the imaging data, we use the photometric pipeline that is presented in more detail in ref. ⁴⁴. Briefly, the raw data has been reduced using the public grism redshift and line analysis software Grizli⁴³, which masks imaging artefacts, provides astrometric calibrations based on the Gaia Data Release 3 catalogue¹³ and shifts images using Astrodrizzle. The photometric zero-points are applied as described in ref. ³⁴. RGB image created using six filters of NIRCam observation of the Firefly Sparkle is shown in Fig. 1. We used images from which bright cluster galaxies and intracluster light have been removed, as described in ref. ²⁵. The methodology for modelling and removing diffuse light from cluster galaxies and intracluster light (ICL) is presented in ref. ²⁵. The NIRCam depths (0.3′ diameter aperture) for F090W, F115W, F150W, F200W, F277W, F356W, F410M and F444W are 7.2, 6.6, 5.2, 4.4, 3.0, 2.9, 5.5 and 4.3 nJy, respectively, and the NIRISS depths for F115WN, F150WN and F200WN are 3.6, 4.3 and 4.0 nJy, respectively⁴¹.

Photometry of Firefly Sparkle

We perform photometry in 10 JWST bands (NIRISS: F115WN, F150WN and F200WN; NIRCam: F115W, F150W, F200W, F277W, F356W, F410M and F444W) in which the Firefly Sparkle is detected from their morphological fit with GALFIT. In other JWST and HST filters, the Firefly Sparkle is not or barely detected; hence, we place upper limits for the entire source. As the object is resolved into at least 10 distinct clusters and a diffuse galaxy component, we perform a morphological fit using Galfit¹⁰ to extract the photometric information.

Point spread functions are extracted empirically by median stacking bright, isolated, non-saturated stars following the methodology described in ref. ²⁸. Convolution kernels for homogenizing all data to the F444W resolution are created with photutils.psf.matching using a SplitCosineBellWindow() windowing function to remove high-frequency noise, which results from floating-point imprecision when taking the ratio of Fourier transforms. We optimize the shape of each window function to minimize the median residual between convolved stars from each source filter that is convolved and stars from the target F444W filter.

For the morphological fit, we create 10″ × 10″ postage stamps in all 10 filters from the BCG-subtracted images. We determine the priors for the centres of the 10 clusters by visual inspection. Although nine out of the ten appear as point sources, FF-4 has an elongated shape and appears unresolved. We first determine the central coordinates of the 10 clusters and the arc by fitting (1) an elliptical Gaussian for FF-4; (2) nine point sources for the other nine clusters; and (3) another elliptical Gaussian with the bending mode turned on for the diffuse arc to the F115W image, which has the highest resolution (smallest PSF). The free parameters are the centres and total fluxes of all the components, the radius and axis ratio of FF-4, and the radius, axis ratio and bending mode (B2) of the arc. The initial guesses for the coordinates were determined by visual inspection of the F115W image. Once we obtain the fitted central coordinates of all the components from F115W, we again fit all 11 components in F444W, which has the highest signal-to-noise ratio for the arc and FF-4, to determine the radius, axis ratio, position angles of the ellipses, and the bending mode B2 of the arc.

We use the best-fit centre coordinates from F115W as the central coordinates in all the filters. However, instead of fixing the central coordinates, we allow GALFIT to fit for them in every filter within a very narrow range of ±0.5 pixels (0.02″) to account for the uncertainty in 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 F444W fit. We fix the morphology of FF4 also with radius = 0.59″, axis ratio = 0.1 and position angle = −53°.

We now fit all 11 components in all 10 filters to determine their fluxes. The resulting models and residuals are shown in Extended Data Fig. 1. Residuals from the fits are negligible, as shown by χ²/ν ~ 1 in the GALFIT fits in all filters. This confirms the original visual impression that nine of the ten clusters are unresolved and an additional smooth component is present.

To derive the uncertainty in our flux estimation, we inject the full Firefly Sparkle model in 100 random locations in our 10″ × 10″ postage stamps (avoiding the edge) and refit with the exact same setting of GALFIT. We find no significant systematic offset between the fitted flux and the injected flux for any of the 11 components, in any of the filters, showing that our photometric technique is robust to background variations across all filters. The uncertainty in the photometry is calculated from the bi-weight scale of the 100 refitted fluxes. The resulting photometry and the RGB image of the model and the residual are shown in Extended Data Fig. 1. The agreement between NIRISS and NIRCam fluxes in the three overlapping filters is another confirmation of the robustness of photometry. We have used updated zero-points³⁴ and corrected for Milky Way extinction using the colour excess E(B − V) = 0.0272 from ref. ⁶ and assuming the extinction law in ref. ³⁵ using the factor between the extinction coefficient and colour excess R_V = 3.1.

Spectroscopy extraction and spectral fitting

NIRSpec spectroscopy has been acquired for MACS J1423.8 + 2404 and spectra were obtained for the Firefly Sparkle, FF-BF and FF-NBF. The spectra for the FF-BF were part of the sample in ref. ²³, with z_spec = 8.2953 ± 0.0005. The spectra were observed using the PRISM/CLEAR disperser and filter, through three Micro-Shutter Assembly (MSA) masks per cluster with a total exposure time of 2.9 ks per MSA configuration.

The NIRSpec data were processed using the STScI JWST pipeline (software v.1.8.4 and jwst_1030.pmap) and the msaexp package³¹. We used the standard JWST pipeline for the level 1 processing, in which we obtained the rate fits files from the raw data. We enabled the jump step option expand_large_events to mitigate contamination by snowball residuals and used a custom persistence correction that masked out pixels that approach saturation within the following 1,200 s for any readout groups. We then used msaexp for level 2 processing, for which we performed the standard wavelength calibration, flat-fielding, path-loss correction and photometric calibration and obtained the 2D spectrum before background subtraction. As the central and upper shutters contain different clusters (see Fig. 2a to find the shutter positions), we need custom background subtraction to avoid self-subtraction. We did this by building the background 2D spectrum by stacking and smoothing the sky spectrum in the empty pixels and obtained the background subtracted 2D spectrum of Firefly Sparkle. We confirmed that this custom background subtraction method works as well as a standard drizzle background subtraction method used in the literature³³, using a well-isolated galaxy spectrum from the CANUCS observation (Asada et al., in prep.). We finally extract the 1D spectrum separately in slit 1 and slit 2, by collapsing the 2D spectrum using an inverse-variance weighted kernel following the prescription in ref. ²⁴. We verified that the uncertainty array of the 1D spectrum has the appropriate normalization by testing the distribution of spectral fluctuations in an empty sky region and finding the fractions of pixels at >1 and >2σ as expected.

Spectral fitting in Firefly Sparkle slit 1

The resulting 1D spectrum of Firefly Sparkle in slit 1, dominated by the cluster FF-6, is shown in Fig. 2.

The spectrum exhibits a Balmer jump at λ_obs ~ 3.5 μm and a turnover at λ_obs ≲ 1.4 μm, probably because of two-photon emission. These features suggest that the nebular continuum should dominate over the stellar continuum in the rest frame UV to optical spectrum within slit 1 (as found for a z = 5.9 galaxy in ref. ¹²). We thus model the continuum of the spectrum with nebular continuum using the photoionization code CLOUDY v.23 (ref. ⁵). To determine the dust attenuation value in the continuum model fitting, we first measure the Hγ/Hβ ratio by fitting the Gaussian profiles. The ratio agrees well with the case B recombination, and no significant dust attenuation is indicated. Therefore, in the continuum spectral modelling, we use pure hydrogen gas irradiated by an ionizing source having black-body SED without dust attenuation. We vary the effective temperature of the black body (T_eff) and the electron temperature of the (ionized) hydrogen gas (\({T}_{{\rm{e}},{{\rm{H}}}^{+}}\)), and search for the best-fitting model continuum by χ² minimization. In the continuum fitting, we mask out emission line regions and all wavelengths λ_obs < 1.2 μm at which the Lyman break is seen in the slit 2 spectrum, because this region may be affected by a neutral hydrogen damping wing. The best-fit model has log(T_eff/K) = 5.10 and \(\log ({T}_{{\rm{e}},{{\rm{H}}}^{+}}/K)=4.34\), which is fully consistent with the results in ref. ¹². The result of continuum fitting does not change if we consider a slight dust attenuation (A_V = 0.1 mag) in the fitting. As discussed in ref. ¹², the effective temperature of log(T_eff/K) = 5.10 is much hotter than typical massive type O stars and is suggestive of this star-forming cluster having a top-heavy IMF. The IMF of this cluster is further discussed in section ‘SED fitting analysis’.

Note that the UV continuum turnover feature could be because of the absorption from dense neutral hydrogen either in the intergalactic medium (IGM) or in the circumgalactic medium (CGM). However, in the case of slit 1 spectrum, we expect the effect of IGM and CGM damping absorption to be negligible or limited at λ_obs < 1.2 μm based on the blue continuum and sharp drop-out in the slit 2 spectrum (see section ‘Spectral fitting in Firefly Sparkle slit 2’ for details of slit 2 spectrum). Considering the spatial proximity of the slit 1 and slit 2 regions (Fig. 2), we can assume the absorption feature from line-of-sight neutral hydrogen to be the same in the slit 1 and slit 2 spectra. The slit 2 spectrum is rather blue and has a sharp Lyman break starting at λ_obs = 1.2 μm, whereas the slit 1 spectrum shows the turnover starting at λ_obs ~ 1.4 μm. Thus, the turnover feature should not be because of the neutral hydrogen absorption, but rather because of the intrinsic continuum shape of the source. Nevertheless, to avoid the possible effect of the neutral hydrogen absorption, we mask out λ_obs < 1.2 μm in the nebular continuum fitting above (corresponding to 1,290 Å in the rest frame).

Having the model continuum, we subtract the underlying model continuum from the observed spectrum and measure the spectroscopic redshift and emission line fluxes by fitting Gaussian profiles. The best-fitting model spectrum with nebular continuum and Gaussian profiles is shown in Fig. 2b (red solid curve). We securely detect emission lines of [O iii]λλ4959, 5007, Hβ, [Oiii]λ4363, Hγ, Hδ and [Neiii]λλ3869, 3889. We do not find significant detection of [Oii]λ3727 and obtain an upper limit for the flux of this line. There is a tentative detection of the blended line of [Oiii]λλ1661 + 1666, although the spectral resolution of the prism is low at this wavelength making this doublet difficult to securely detect and separate from Heiiλ1640. We use these emission line fluxes to estimate the physical parameters in slit 1. We first estimate the dust attenuation based on Balmer decrements. Both the Hγ/Hβ and Hδ/Hβ ratios are consistent with theoretical predictions in case B recombination²¹ within the uncertainties, suggesting there is no significant dust attenuation (Extended Data Fig. 3, red squares in the left). This result is consistent with the initial measurement before the continuum fitting above and supports the validity of the dust-free assumption in the nebular continuum fitting process. Therefore, we do not correct for dust attenuation in the following measurements of emission line ratios and physical parameters in this section.

We next measure the electron temperature using temperature-sensitive emission line ratios: [Oiii]_4959+5007/[Oiii]₄₃₆₃ and [Oiii]₅₀₀₇/[Oiii]_1661+1666. We assume the electron density to be n_e = 10³ cm⁻³, which is consistent with recent JWST observations of similarly high-z galaxies⁷ and obtain consistent independent temperature measurements within the uncertainties (\({T}_{{\rm{e}},{{\rm{O}}}^{++}}={4.0}_{-0.9}^{+2.6}\,{\rm{K}}\) and \({2.9}_{-0.4}^{+0.7}\times 1{0}^{4}\,{\rm{K}}\), respectively; Extended Data Fig. 3 (right)). Note that because the [Oiii]λλ1661 + 1666 detection is tentative and potentially blended with Heiiλ1640, we consider [Oiii]λ4363 to be more reliable.

We note that in ref. ¹⁶, the authors measured a similar ratio of [Oiii]_4959+5007/[Oiii]₄₃₆₃ in the z = 6 galaxy RXCJ2248-ID to that of slit 1. In ref. ¹⁶, medium resolution spectroscopy was used to determine the electron density directly. They found that when using lines with higher ionization potential than O+, the electron density was higher (n_e ~ 10⁵ cm⁻³) than is typically found from [Oii]λ3727 (ref. ⁷). This high electron density leads to a lower electron temperature for their galaxy of \({T}_{{\rm{e}},{{\rm{O}}}^{++}}=2.5\times 1{0}^{4}\,{\rm{K}}\). Similarly, if we assume the electron density of n_e = 10⁵ cm⁻³ instead for our slit 1 spectrum, the electron temperature from [Oiii]λ4363 becomes \({T}_{{\rm{e}},{{\rm{O}}}^{++}}={3.2}_{-0.96}^{+1.6}\), which is in between the two measurements based on [Oiii]λλ1661 + 1666 and [Oiii]λ4363 when assuming n_e ~ 10³ cm⁻³ above. To consider the possibility of a somewhat higher electron density in the highly ionized region, we adopt the mean value of our two electron temperature measurements (\({T}_{{\rm{e}},{{\rm{O}}}^{++}}=3.5\times 1{0}^{4}\,{\rm{K}}\)) as our fiducial value and propagate the full range of the two measurement uncertainties into the following metallicity measurement.

Based on the electron temperature measurement, we obtained the oxygen abundance from [Oiii]_4959+5007/Hβ and [Oii]₃₇₂₇/Hβ ratios, following the prescription in ref. ⁸. We assume the electron density to be n_e = 10³ cm⁻³. The total oxygen abundance is calculated from O⁺⁺/H⁺ and O⁺/H⁺, and the higher ionizing state oxygen is ignored³⁰. As the [Oii]λ3727 emission line is undetected, we can obtain only an upper limit for O⁺/H⁺, but the upper limit for the abundance of the singly ionized oxygen is negligibly small as compared with the doubly ionized oxygen. We thus derived the total oxygen abundance from O⁺⁺/H⁺, yielding \(12+\log ({\rm{O/H}})=7.0{5}_{-0.37}^{+0.22}\) (\({Z}_{{\rm{gas}}}/{Z}_{\odot }=0.0{2}_{-0.01}^{+0.04}\) assuming the solar abundance to be 8.69; ref. ³⁸).

We also derive the ionization parameters using the ionization-sensitive emission line ratios: [Oiii]₅₀₀₇/[Oii]₃₇₂₇ and [Neiii]₃₈₆₉/[Oii]₃₇₂₇. Following the prescription in refs. ^45,46, we obtain the lower limit for the ionization parameters (log  U) from these two ratios. Both ratios provide a similar limit of log U > −2.0.

All the emission line flux measurements and the derived physical parameters in Firefly Sparkle slit 1 are presented in Extended Data Table 1. We also compare the diagnostic emission line ratios in Firefly Sparkle with those in other galaxy population in Fig. 2d. We use the ionization-sensitive line ratio O32 ([Oiii]₅₀₀₇/[Oii]₃₇₂₇) and the temperature-sensitive line ratio RO3 ([Oiii]_4959+5007/[Oiii]₄₃₆₃) and compare these line ratios with other [Oiii]λ4363-detected galaxies at z = 2–9 from previous JWST observations² and those in the local universe from SDSS observations¹⁴. Extended Data Fig. 3 (middle) presents a similar comparison but uses another ionization-sensitive line ratio Ne3O2 ([Neiii]₃₈₆₉/[Oii]₃₇₂₇) instead of O32.

Spectral fitting in Firefly Sparkle slit 2

In contrast to slit 1, the extracted 1D spectrum in Firefly Sparkle slit 2 does not show nebular continuum features, and the blue continuum is rather smooth with a sharp drop-out because of the Lyman break at λ_obs ~ 1.2 μm. We thus derive the emission line fluxes from the slit 2 spectrum by fitting Gaussian profiles with the continuum being modelled by a constant offset around each emission line. We detect [Oiii]λλ4959,5007, Hβ, Hγ, Hδ, [Neiii]λ3869 and [Oii]λ3727 emission lines in the slit 2 spectrum but do not detect [Oiii]λ4363.

We then derive the physical properties in the same way as done for Firefly Sparkle slit 1 spectrum. We measure the dust attenuation from Balmer decrement, Hγ/Hβ and Hδ/Hβ, and find both line ratios agree well with the predicted ratios under case B recombination (blue squares in Extended Data Fig. 3 (left)). This suggests that the dust attenuation is negligible in the slit 2 spectrum as well, and we do not make a dust attenuation correction.

As we do not detect the temperature-sensitive emission lines of [Oiii]λ1666 or [Oiii]λ4363 in the slit 2 spectrum, we cannot measure the electron temperature and the metallicity from the direct-temperature method. We thus obtain only the upper limit for the electron temperature (\({T}_{{\rm{e}},{{\rm{O}}}^{++}}\)) from the non-detection of [Oiii]λ4363. The electron temperature in Firefly Sparkle slit 2 is shown to be \({T}_{{\rm{e}},{{\rm{O}}}^{++}} < 1.8\times 1{0}^{4}\,{\rm{K}}\) (1σ) or <4.5 × 10⁴ K (3σ). To visualize the difference in physical properties in slit 1 and slit 2, we show the diagnostic emission line ratios of Firefly Sparkle slit 2 in Fig. 2d and Extended Data Fig. 3 (middle) as well.

SED fitting analysis

SEDs derived from our photometry were analysed using a slightly modified version of the Dense Basis method^18,47 to determine non-parametric SFHs, masses and ages for our sources in Firefly Sparkle. We adopt the Calzetti attenuation law⁴⁸ and a Kroupa IMF³² with a flat prior for the high-mass slope α ∈ [1., 4.]. We run fits using both the MILES stellar libraries²⁹ and MESA Isochrones and Stellar Tracks (MIST; ref.  ¹⁷), as well as the Binary Population and Spectral Synthesis (BPASS; refs. ^26,36) models to consider for the presence of binary populations. As the latest BPASS version in FSPS (-bin-imf135all 100) assumes a Salpeter IMF with an upper mass cutoff of 100M_⊙ and does not allow for a varying IMF, we only vary the top-heavy slope of the Kroupa IMF in the MILES + MIST runs with an upper mass cutoff of 120M_⊙. This should be considered while comparing the physical properties from the two runs, as allowing for a varying IMF based on the MILES + MIST configuration results in lower stellar masses for those runs because they are preferentially fit with top-heavy SSPs with a greater fraction of light coming from more massive stars. We fix the redshift to that found from the NIRSpec Prism spectroscopy by the [Oiii] λ4959 line at z_spec = 8.296 ± 0.001. All other parameters are left free. We run the SED fits in two configurations to account for different possibilities of the nature of the individual clusters:

1. 1.

   SSP fits: to account for the possibility that the individual clumps are star clusters, which is likely given the physical scales of the clusters and the extreme emission lines in the spectra, we modify the code to fit for instantaneous bursts of star formation, described by SSPs. In this case, we assume a flat prior in the log age of the SSP from 10⁵ years to the age of the universe at z_spec = 8.296 ± 0.001 instead of the non-parametric defaults for the SFH in Dense Basis.

2. 2.

   Non-parametric SFH fits (Dense Basis): to fit the diffuse body of the galaxy and to account for the possibility that the clusters are nuclear star clusters or remnants of previous mergers, we also run fits with non-parametric SFHs with a Dirichlet prior. The main advantage of using Dense Basis with non-parametric SFHs is that they allow us to account for flexible stellar populations, which is important at these redshifts⁴⁹ because star formation is expected to be stochastic and may be underestimated if fit using traditional parametric assumptions^39,50.

We perform our fitting in two stages—we initially perform a joint spectrophotometric fit to the NIRSpec Prism spectrum along with the HST + NIRISS + NIRCam photometry in the slits in which both exist (Extended Data Fig. 4). We correct for slit loss considering two factors—the amount of light lost due to the changing PSF as a function of wavelength and an overall multiplicative correction to match the spectrum against the photometric measurements. We modify the default Dense Basis method in this stage to additionally fit for the slope at the massive end of the IMF, the gas-phase metallicity and the ionization parameter, using the relevant parameters from FSPS (imf3, gas_logz and gas_logu). Doing so allows us to substantially constrain priors on star formation rate, IMF, dust, ionization parameter and metallicity that we then use to fit the photometry. We find that the fits are consistent with negligible dust attenuation, consistent with our estimates from measuring the Balmer decrement. We also find that our fits rule out the part of parameter space consistent with the canonical Chabrier-like or Kroupa-like IMF (with the high-mass slope ≈ 2.3) in favour of more top-heavy slopes of about \({1.5}_{-0.6}^{+0.7}\) for slit 1, which contains portions of clusters 3, 4, 5 and 6. We find weaker constraints from the spectrum for slit 2, which still skews towards top-heaviness but with large uncertainties of about \({1.7}_{-0.7}^{+0.9}\). Finally, we find estimates of both stellar and gas-phase metallicities to be sub-solar, consistent with estimates from the line ratios.

Using our photometry (Extended Data Table 4), we now determine the stellar properties of each individual component by running a second set of fits using the same set of parameters that are used to fit the spectrophotometry. Although parameters such as the metallicity and ionization parameter are only loosely constrained by these fits, we obtain parameter estimates for the stellar masses, star formation rates and ages of the individual star clusters with uncertainties that marginalize over the variations in the other parameters.

Both photometry and corresponding fits to the SED fit are shown in Extended Data Fig. 5, with variations in the stellar mass, age and reduced χ² of the fits for each of the four scenarios (SSPs fits with MILES + MIST and BPASS, and Dense Basis fits with MILES + MIST and BPASS) shown in Extended Data Table 2. All 10 components have intrinsic (corrected for magnification) stellar masses of about 10⁵–10⁶M_⊙ and sSFR of 10⁻⁷ yr⁻¹. Although the error bars are large, the distinct colours of the clusters hint at different formation times. Although the smooth component contains a large fraction of the stellar mass, the bulk (about 57%) lies in the clusters. Extended Data Table 3 lists the physical properties of the individual components as well as the full Firefly Sparkle, BF and NBF galaxies.

We find that the SSP fits are generally less massive compared with the Dense Basis fits, because the light from the SED is modelled by a single epoch of star formation instead of an extended episode. As light from the massive stars responsible for young star formation are much brighter than older stellar populations, they can describe the observed SED with a lower mass. However, the SSP fits often cannot capture both the UV slopes and the nebular emission in the rest-optical, as seen for clusters 1, 3, 7 and 8 in Extended Data Fig. 5 and often approximate it using a Balmer break, leading to posteriors consistent with much older ages than the median values.

Although the t_age from SSP and t₅₀ from Dense Basis fits (Extended Data Table 2) may seem inconsistent, it is important to note that the Dense Basis fits for most star clusters indicate a sharp burst of star formation within the past 10 million years (Extended Data Fig. 6). By design, an SSP is biased towards this recent burst, whereas a non-parametric SFH can accommodate extended episodes of star formation. However, with our current data, we cannot distinguish between extended SFH in the star clusters and the contribution of light from the diffused arc.

The masses of the clusters also scale with the top-heaviness of the high-mass end of the IMF in the MILES + MIST fits, with lower masses for more top-heavy IMF values as that scale the amount of light from massive stars. In comparison, the BPASS fits in the current setup are done at the canonical Kroupa IMF, leading to higher masses for those fits. At the same IMF slope, the masses are comparable within uncertainties for the different SPS models, and the sSFR and age/t₅₀ values are consistent even marginalizing over the IMF posteriors. Given the observational constraints and the χ² from the fits in Extended Data Table 2, it is not currently possible to definitively rule out any of the current fitting approaches.

Lens modelling

We use Lenstool⁹ to build a strong lensing model of the MACS 1423 cluster, to be fully presented in Desprez et al. (manuscript in preparation). This model is constrained with the three multiple image systems that were leveraged in ref. ³, for which we provide additional information obtained from the CANUCS data. The first two systems are those presented in ref. ²⁷, one at z = 2.84 for which we account for the two clusters visible in the four images of the objects, and the second one with three images at z = 1.779 for which we identify another cluster in the two northernmost images for improved constraints. The last system is composed of five images¹¹ for which we provide a new spectroscopic redshift measurement of z = 1.781 that is in agreement with photometric and geometric redshifts previously measured.

The different mass components are parameterised as double Pseudo-Isothermal Elliptical (dPIE) profiles⁴. The model is composed of a large cluster scale mass halo, an independent galaxy scale centred on the brightest cluster galaxy and small galaxy scale mass components to account for the contribution of all cluster members that follow a mass–luminosity scaling relation²². For all galaxies, their positions, ellipticities and orientations have been fixed to these measured from the images. The final best model manages to reproduce the position of the input multiple images with a distance rms of 0.46″.

Magnifications are obtained by generating convergence and shear maps around the Firefly Sparkle with a size of 20″ and a resolution of 10 milli-arcsec per pixel. Uncertainties in the magnifications are computed from 100 randomly selected models from the optimization of Lenstool after its convergence around the minimum χ². The numbers provided in Extended Data Table 3 are the median and ±1σ limits on the distribution of the 100 values obtained at the position of each cluster. We measured the average magnification of the FF-arc by using the GALFIT model of the arc (in F200W) and selecting all pixels with flux >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 the arc (μ = 24.4 ± 6.0).

The source plane reconstruction is made using the best GALFIT model to compute the source plane positions and magnification 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 point sources convolved with the appropriate PSFs to the diffused source plane model at the source plane positions of the star clusters with the demagnified fluxes. This process is repeated to generate source plane models in all filters. We also generate a mass map using the same prescription, replacing the demagnified fluxes with the demagnified masses. The resulting source plane RGB image and mass map are shown in Fig. 4c,d.

Size and surface density of star clusters

We now investigate the spatial properties of the star clusters. Nine out of the ten star clusters are unresolved even in our highest resolution 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 the NIRCam F115W PSF (0.02) to set an upper limit on the size of all 10 star clusters. To determine the upper limits of the sizes of unresolved sources, we use the tangential eigenvalue of magnification 1/∣λ_t∣, which ranges between 14 and 24. This results in a size upper limit between 4 pc and 7 pc. The central star clusters have the highest magnification and the smallest upper limits, whereas the ones near the two ends of the arc have the lowest. We use the upper limit on sizes and the demagnified stellar masses to calculate the lower limit on stellar surface densities as shown in Fig. 3b.

Abundance matching for MW and M31 progenitors

To estimate the range of stellar masses of progenitors of both MW-mass and M31-mass galaxies at higher redshift, we adopt a semi-empirical approach combining both simulations and observations. We assume an evolving co-moving number density with redshift, as determined by the abundance matching code from ref. ²⁰, with z = 0 number densities of \(\log (n/{{\rm{Mpc}}}^{3})=-\,2.95\) and \(\log (n/{{\rm{Mpc}}}^{3})=-\,3.4\), respectively, for MW and M31 mass analogues. The code calculates a past median galaxy number density at z₂, given an initial number density at z₁, using peak halo mass functions. As the merger rate per unit halo per unit Δz is roughly constant, the evolution of the cumulative number density of progenitors of any given galaxy is a power law, with the change described by (0.16Δz) dex.

In ref. ²⁰, peak halo mass functions are used because the resultant median number densities are less affected by the scatter in stellar mass and luminosity. However, this scatter does affect the 1σ errors in cumulative number density. The 1σ or 68 percentile range grows with increasing redshift, but this growth is also higher for more massive galaxies.

As the code from ref. ²⁰ does not calculate stellar masses, we obtain the stellar mass ranges of the progenitor populations of MW and M31 analogues using stellar mass functions (SMFs) from various surveys^15,19,40. We take the median cumulative number densities at each Δz to find the stellar mass associated with that number density from the corresponding SMF. Moreover, the 1σ errors on the given number density for each redshift are then used to determine the 1σ errors on the stellar mass of the progenitors. At z = 8.3, the median stellar mass of MW progenitors is \(\log ({M}_{\star }/{M}_{\odot })=6.4\pm 0.7\) and the median stellar mass of M31 progenitors is \(\log ({M}_{\star }/{M}_{\odot })=6.9\pm 0.8\). The Firefly Sparkle with a stellar mass of \(\log ({M}_{\star }/{M}_{\odot })={7.0}_{-0.3}^{+1.0}\) is definitely within 1σ stellar mass range of both Milky Way and M31 progenitors. More details on the progenitor matching technique can be found in ref. ³⁷.