Main

The Cosmic Gems arc (SPT0615-JD1) was initially discovered in Hubble Space Telescope (HST) images obtained by the Reionization Lensing Cluster Survey (RELICS) of the lensing galaxy cluster SPT-CL J0615–5746 at z = 0.972 and reported as a redshift z = 10 candidate8.

A recent James Webb Space Telescope Near Infrared Camera (JWST/NIRCam) imaging campaign of this field has observed the Cosmic Gems arc with eight bands covering the 0.8–5.0 μm range (Methods). Spectral energy distribution (SED) fitting to the James Webb Space Telescope (JWST) photometry indicates that the Cosmic Gems galaxy has a fairly young stellar population with a mass-weighted age of less than 79 Myr and a lensing-corrected stellar mass in the range of 2.4–5.6 × 107M⊙, with low dust extinction (AV < 0.15 mag) and metallicity (<1% Z⊙) (ref. 1).

The far-ultraviolet (FUV)-to-optical rest frame of Cosmic Gems arc shows bright clumpy structures and extended faint emission over a 5″-long arc (Fig. 1). The symmetry between the south–east (hereafter Img.1) and the north–west (Img.2) part of the arc uncovers two lensed mirror images of the galaxy, implying that the Cosmic Gems arc is observed at very high magnification on the lensing critical curve. Four independent magnification models have been created to account for the galaxy appearance. All the models successfully reproduce the z = 10.2 critical line crossing the Cosmic Gems arc (Methods).

Fig. 1: The Cosmic Gems arc in a JWST colour composite.
figure 1

The filter combination shows the rest-frame UV, blue optical wavelengths (1,200–2,800 Å). The arc is extended over 5″. A foreground galaxy at redshift zphot = 2.6 is visible above and to the right. The field of view is rotated North up. Top right, a magnified inset of the centre of the arc in which the brightest star clusters are located, highlighted by the white square on the left image. Two mirror images are observed because of gravitational lensing. Lensing critical curves based on three models are shown bisecting the arc. Bottom left, each star cluster is labelled in a grey scale FUV rest-frame image of the galaxy. Bottom right, source plane reconstruction of the core of the galaxy in which the star clusters are located showing their relative sizes and positions. The physical distance between A and E and A and D is about 40 pc. Note that there is some uncertainty in the source positions parallel to the lensing caustic. Scale bar, 10 pc (bottom right).

Full size image

Five star cluster candidates are uniquely identified in Img.1. In Img.2, three sources are distinguishable in the F150W filter (Fig. 1), along with a fourth fainter source (E.2). The appearance of Img.2 is probably perturbed by further lensing effects because of the northern galaxy at z ≈ 2.6 (visible in the top right corner) and possibly by an undetected small scale perturber closer to the arc1. Source D.2 is possibly blended with C.2 and is therefore identified for the remaining of the analysis as C.2 + D.2. The source E.2 is only detected at a 2σ level and not included in this analysis (Methods). The observed projected distance between the A.1 and E.1 clusters in Img.1 is about 0.65″. Using forward modelling and Lenstool-A predictions presented in the Methods, we find that their physical distance is \(4{2}_{-5}^{+29}\,{\rm{pc}}\). The star clusters are all located within this compact region (Fig. 1). The size of this region is similar to that reported for individual stellar clumps observed in the moderately lensed galaxies at redshift z ≈ 10 (refs. 2,5), suggesting that they likely harbour star clusters within them.

The intrinsic physical properties of these five star clusters are particularly meaningful for probing proto-globular cluster formation mechanisms as well as their potential evolution. As described in the Methods, we find that A.1 is only marginally resolved, with an observed effective half-light radius Reff,obs = 0.6 px, whereas C.1 and B.2 are consistent with being unresolved. For the latter sources, we assumed an upper limit to their radii coincident with the half-width half maximum (HWHM) of the stellar PSF in the F150W (0.025″ = 1.25 px). To derive lensing-corrected Reff, we assumed the predicted Lenstool-A tangential magnifications at the location of the star clusters. The five star clusters have intrinsic Reff close to 1 pc (within uncertainties; Table 1). Using the other lensing models produces similar size ranges (0.3–0.9 pc for Lenstool-B and 0.3–1.2 pc for Glafic). Independent intrinsic sizes (Reff,FM in Table 1) have been derived by projecting the star cluster shapes from the source plane into the image plane. The latter method recovered intrinsic sizes in excellent agreement, within the uncertainties, with those measured in the image plane, strengthening the reliability of the derived values.

Table 1 Estimated physical properties of the Cosmic Gems arc star clusters
Full size table

The star clusters have been fitted with BAGPIPES9 and PROSPECTOR10. We tested different star formation history (SFH) assumptions that simulate a single burst (inherent to the small sizes of the stellar systems analysed), different high-mass limits of the initial mass function (IMF) and models with stellar binaries (Methods). Despite the assumptions, the resulting physical properties of the clusters (ages, masses, extinction and metallicities) are in reasonable agreement. In the analysis presented here, we use the physical values derived with an SFH based on a single exponential decline with τ = 1 Myr (referred to as BAGPIPES-exp, Extended Data Table 2).

The recovered ages of the star cluster candidates are between 9 Myr and 36 Myr. The age range suggests that star formation has been propagating within this compact area of the galaxy for a few tens of Myr. The measured rest-frame UV slopes of the star clusters (β between −1.8 and −2.5, with Fλ ∝ λβ; Extended Data Table 1) are similar to those found for more evolved star clusters in the Sunburst arc at redshift 2.37 (ref. 11). Although the Cosmic Gems clusters are not extremely young, they have likely delivered large amounts of energy and momentum to their host galaxy.

The lensing-corrected stellar masses range between 1.0 × 106M⊙ and 2.6 × 106M⊙, for a total combined stellar mass of 8.3 × 106M⊙. The total mass of the clusters is close to 30% of the total stellar mass of the host. As the mass-weighted age of the galaxy and those of the star clusters are comparable, we can extrapolate the cluster formation efficiency (CFE)12 to be around 30%. A caution note is necessary because the mass estimates (both for the galaxy and star clusters) are subjected to magnification values and SED fit uncertainties, making the quoted CFE uncertain. A more direct way to establish the CFE is to use the fraction of observed FUV light in star clusters with respect to the host. This quantity is not affected by the same degeneracy as the mass estimates and, thus, is a more reliable indicator of the CFE, under the assumption that the FUV light is produced by stellar populations formed during a similar timescale (as we find here). The analysed star clusters account for about 60% of the total F150W flux of the host extracted within an elliptical Kron aperture (0.51 ± 0.01 μJy, corresponding to an intrinsic FUV ABmag of −17.8 after lensing correction1), thus reinforcing the conclusion that star formation in star clusters is the main mode for the Cosmic Gems arc and high-redshift galaxies with similar physical properties. This observationally driven conclusion is supported by high-resolution numerical simulations13 and analytical models14 that find that compact star clusters with sizes of 0.5–2 pc are the dominant star formation mode in the first low-metallicity dwarf galaxies. The compactness of the star clusters also seems to drive the leakage of hydrogen ionizing radiation from their natal molecular cloud15, making the star clusters observed here potential contributors to cosmic reionization. Massive star clusters similar to those observed in the Cosmic Gems arc are predicted by the feedback-free starburst model in ref. 16 and could be at the root of the super-Eddington conditions necessary to launch strong outflows in short timescales17, both models aimed to explain the bright UV luminosity reported for z > 9 galaxies.

The resulting stellar surface densities of the Cosmic Gems are around 105M⊙ pc−2 (Table 1). Consistent physical properties have been reported in star clusters detected in the Sunrise arc at z ≈ 6 (ref. 18) and the Sunburst arc at z = 2.37 (ref. 19) (Fig. 2). Using the derived ages, masses and intrinsic sizes, we also determine whether these stellar systems are gravitationally bound. According to the framework introduced in ref. 20, a star cluster is considered bound if its age is greater than the crossing time of the system (where \({t}_{{\rm{cross}}}=10\sqrt{{R}_{{\rm{eff}}}^{3}/GM}\)), or in other words, under the assumption of virial equilibrium, a cluster is gravitationally bound if Π = Age/tcross > 1. The log(Π) values reported in Table 1 are all significantly larger than unity, indicating that we are detecting gravitationally bound star clusters in an early galaxy, 460 Myr after the Big Bang. This conclusion is valid in spite of the uncertainties inherent to physical quantity estimates as well as lensing models.

Fig. 2: Cluster stellar surface density versus half-light radius Reff.
figure 2

The star clusters in the zphot ≈ 10.2 Cosmic Gems arc are shown colour-coded by their dynamical age, log(Π). The plotted values have been derived using the reference lensing model. Predictions by other models do not change the observed trends. The error bars do not account for magnification uncertainties. However, refer to Extended Data Fig. 4, in which we show that magnification uncertainties do not affect these results. Other gravitationally bound star clusters detected in lensed galaxies at z ≈ 6 (ref. 18) and z = 2.38 (ref. 19) are included, along with the z = 0 Milky Way globular clusters 20 and young star cluster properties of star-forming spiral galaxies in the Local Volume (distance < 16 Mpc, ref. 7). Lines of equal mass show the change in density as a function of Reff. GCs, globular clusters; YMCs LVL, young massive clusters in the Local Volume.

Full size image

The Cosmic Gems arc clusters (Fig. 2) have substantially higher stellar densities and smaller sizes than typical young star clusters observed in the local Universe7 as well as global clusters in the Milky Way21. The offset with respect to young star clusters in the local Universe is expected because the conditions under which star formation operates in reionization-era galaxies are more extreme (for example, galaxies are more compact, harbour harder ionizing radiation fields, and reach higher electron densities and temperatures6,22). The offset with respect to local global clusters could be explained in terms of dynamical evolution. Global clusters are hot stellar systems in which stars continuously exchange energy and momentum. Three different internal mechanisms contribute to their dynamical evolution over a Hubble time: (1) mass loss due to stellar evolution; (2) relaxation due to N-body interactions; and (3) formation and dynamics of stellar black holes (SBHs)23,24. Mass loss due to stellar evolution drives the adiabatic expansion of global clusters under the condition of virial equilibrium. A typical mass loss of 50% will expand the initial radius of the Cosmic Gems proto-globular clusters by a factor of 2, whereas the density will decrease correspondingly by a factor of 8 (ref. 23). Relaxation time scales (shortened by the presence of SBHs25) will also contribute to their expansion. Finally, external tidal fields will further affect the dynamical evolution of these bound stellar systems, which appear to be bona fide proto-globular clusters.

Very dense stellar clusters (Σ* ≈ 105M⊙ pc−2, which for Reff = 1 pc correspond to ρh ≈ 105M⊙ pc−3), similar to those detected in the Cosmic Gems arc, are predicted to form in low metallicity and highly dense gas26, in which radiative pressure cannot counteract the collapse, resulting in extremely high star-formation efficiencies (about 80%; ref. 27). The high stellar densities found in these proto-globular clusters imply a notable increase in stellar black hole mergers in their interiors28,29 and therefore pave the way to intermediate-mass black hole seeds30. With stellar masses greater than 105M⊙, these star clusters naturally harbour Wolf–Rayet and very massive stars31, and because of their elevated stellar densities, satisfy the necessary condition to form supermassive stars in runaway collisions within their cores32. These different classes of stars are among the potential polluters that could explain the observed nitrogen enrichment in the ionized gas of high-redshift galaxies33, possibly linked to the formation of the chemically enriched stellar populations ubiquitously found in Milky Way global clusters34.

Cosmological simulations that focus on Milky Way disk-like assembly find that most of its global cluster population form at redshift z < 7 (refs. 35,36,37), suggesting that these star clusters forming at z ≈ 10 might build up the global cluster populations of more massive early-type galaxies in the local Universe. It is difficult to predict whether the proto-globular clusters of the Cosmic Gems arc will survive a Hubble time. Their chances would be highly enhanced if they were ejected into their host halo during dynamical interactions12.

Methods

Size and flux measurements

The JWST/NIRCam38 observations of the SPT-CL J0615–5746 galaxy cluster were obtained in 2023 September (GO 4212: principal investigator Bradley) using four short-wavelength (SW) filters (F090W, F115W, F150W and F200W) and four long-wavelength (LW) filters (F277W, F356W, F410M and F444W) spanning 0.8–5.0 μm. Each filter had an exposure time of 2920.4 s. The data were reduced with the GRIZLI (v.1.9.5) reduction package39. They suffered from strong wisps40 and required special background subtraction as described in ref. 1. The final data are in units of 10 nJy. The NIRCam SW images were drizzled to a pixel scale of 0.02″ per pixel and the LW images were drizzled to a pixel scale of 0.04″ per pixel. For further details on the observations and image reduction, please see ref. 1. We assume throughout the analysis a cosmology with H0 = 67.7 km s−1 Mpc−1 and Ωm = 0.31 (ref. 41). Under these assumptions, 1 SW pixel (0.02″) corresponds to 83.7 pc at z = 10.2.

We derived the star cluster radii and multiband photometry by applying the method published and tested in refs. 42,43,44. We simultaneously fitted for the shape of the light distribution, the flux and the local background (including the galaxy diffuse light) of each of the identified clusters in the reference filter, F150W, which offers the sharpest view of the clusters. Empirical PSFs in all filters used in this work were built by selecting stars in each band. We fitted the PSF out to 0.4″ with an analytical expression, which was then convolved with varying two-dimensional (2D) Gaussians to create a grid of models. Intrinsic sizes (de-convolved by PSF) have then been derived by fitting the observed light distribution of each cluster with this grid of 2D Gaussian models in the F150W. Owing to the large shear effects, compact sources might appear resolved in the shear direction (ystd) (ref. 45). We assumed that the measured major-axis ystd of the 2D Gaussian ellipse is the standard deviation of a 2D circular Gaussian that we translated into the observed PSF-deconvolved effective radius, \({R}_{{\rm{eff}},{\rm{obs}}}={y}_{{\rm{std}}}\times \sqrt{2\,{\rm{ln}}(2)}\). The derived Reff,obs are reported in Table 1. The de-lensed intrinsic effective half-light radii, Reff, have been determined by dividing Reff,obs by the tangential magnification μtan (Extended Data Table 3).

For each star cluster, the flux in the reference filter has been determined by integrating the fitted shape and subtracting the local background. We then measured the fluxes in the other bands by convolving the derived intrinsic shape in F150W with the empirical PSF of the respective bands. Fitting this model to the other bands, we let the centre, normalization and local background (including the galaxy diffuse light) as free parameters.

The intrinsic sizes and observed fluxes in the reference filter F150W were derived using a cutout box centred on the source with a size of 11 × 11 px (about four times the FWHM) for A.1 and A.2. Larger box sizes did not produce noticeable differences in the output sizes and fluxes. Owing to their proximity, B.1, C.1, D.1 and E.1 have been fitted simultaneously within a box of 15 × 15 px. A larger box size produces consistent values within the uncertainties of measurements for B.1 and C.1, whereas E.1 gets increasingly elongated, affecting the fit of D.1. To avoid this degeneracy, we fix the box size to 15 × 15, which would correspond to fix the source ellipticity of the faint E.1 to 2. Similarly, B.2, C.2 and D.2 were fitted simultaneously within a box of 11 × 11 (changing the box size does not produce noticeable effects on the recovered parameters). We did not detect two maxima at the location of C.2 and D.2, so we allowed the fit to optimize the centre of a second hidden source. We also repeated the fit of this region by assuming only one source. Both approaches produced similar residuals. The flux extracted by assuming only one source is comparable within uncertainties to the flux extracted by fitting for C.2 and D.2. Owing to the degeneration in identifying the position of D.2, we extracted the physical properties by fitting only one source that we refer to as C.2 + D.2. Owing to the faintness of E.2 (2σ), our method did not produce meaningful constraints. We, therefore, excluded E.2 from our analysis. In the BAGPIPES-exp fit, we find that C.1 and D.1 have similar ages and a combined total mass of about 2 × 106M⊙. C.2 + D.2 has a slightly older age (but in agreement within 1σ) than C.1 and D.1. The total mass and size of C.2 + D.2 is a factor of two higher than their counterparts, corroborating the idea that the two star clusters are blended in Img.2.

The size uncertainties were derived by bootstrapping the fit of the source taking into account the RMS of the local background. The photometric errors include the latter uncertainties as well as the sum in quadrature of the local background variance estimated within the box in which the sources have been fitted. Aperture corrections have been extrapolated up to 0.4″ in all bands.

In Extended Data Fig. 1, we show the best model of the star clusters and the residual image in the reference filter and two more bands. The extraction of the sources does not produce significant artificial residuals above the RMS of the image.

Independent measurements of intrinsic Reff have been obtained following the forward modelling method in ref. 46. Briefly, this method creates a model of the galaxy in the source plane and then projects that model into the image plane. After convolving with the measured empirical PSF, the image plane model is compared with the observed data. The source plane model parameters are first optimized using a downhill simplex algorithm, then sampled using an MCMC with the Python package emcee47.

For SPT0615-JD1, the source plane model consisted of five Sersic profiles centred on the five identified clumps A.1–E.1 (Extended Data Fig. 2). No diffuse component of the arc has been included in this analysis because of the faintness of this component with respect to the clusters. Separately, we modelled clumps A.2–C.2 on the other side of the lensing critical curve. Uncertainties in the Lenstool-A lens model resulted in slight offsets between the source plane positions of clumps on the opposite sides of the critical curve, which prevents simultaneous fitting of the two images of the arc. We found similar results for clump sizes on both sides of the critical curve, with clump radii ranging from 0.7 pc to 1.1 pc (Table 1).

SED fitting analysis

We performed SED fitting with BAGPIPES9 and tested the derived physical properties against different assumptions, as well as with a different software PROSPECTOR10. For all the runs, we fixed the redshift at z = 10.2, as measured by ref. 1. The standard stellar population templates were reprocessed with CLOUDY to generate nebular continuum and line information (see ref. 48 for comparisons of the two code implementations). In both codes, we assumed a Kroupa IMF, unless otherwise specified. We constrained SFHs to prescriptions that reproduce a short burst in all tests except the one in which τ is set free to vary. The short burst assumption is in agreement with the studies of stellar cluster and global cluster populations in the local Universe (refs. 12,34). The recovered median of the posterior distributions of age, mass, AV, metallicity and associated 68% uncertainties are reported in Extended Data Table 2. We let the ionization parameter, U, to change between −2 and −3.5. We assumed a Calzetti attenuation49 but tested also the SMC extinction. In the reference set, used to produce results reported in Fig. 2 and Table 1 and referred to as BAGPIPES-exp, we assumed an exponential decline with a very short τ = 1 Myr and Calzetti attenuation. Extended Data Fig. 3 shows the observed SEDs of the five star clusters identified in Img.1 (black dots with uncertainties). When available, we include the observed SED of the corresponding clusters in Img.2 (orange stars with associated errors). The latter have been normalized by the median flux ratio in the six bands of the corresponding source in Img.1 to match the flux level while preserving the intrinsic SED shape. The best spectral and integrated photometry model obtained for the BAGPIPES-exp fit is included. The overall shapes of the observed SEDs of mirrored clusters in both images are similar within uncertainties, confirming the symmetry.

To check the consistency of the derived physical properties of the cluster, we made different assumptions. The outputs are summarized in Extended Data Table 2 in which we list the median and 68% values produced by the different fits. Changing the attenuation prescription from Calzetti to SMC produces noticeably smaller Av, but all the recovered parameters are still within the 68% uncertainties associated with the recovered values. In BAGPIPES-burst, we assumed a single burst. We recovered slightly older ages and larger masses (but noticed uncertainties) that would prefer higher stellar surface densities and older dynamical ages, confirming that we were looking at dense and bound star clusters. In a third SED fitting set, BAGPIPES-BPASS, we used BPASS v.2.2.1 SED templates50 and the fiducial BPASS IMF with the maximum stellar mass of 300M⊙ and a high-mass slope similar to that in ref. 51. Also, this model reproduces values that are very close to the reference value, suggesting that the clusters are compatible with being slightly older and therefore less sensitive to the presence of very massive stars and binary systems in their light (and the limitation of fitting only six broad and medium bands covering FUV-blue optical). Letting τ = free (we report mass-weighted parameters in Extended Data Table 2) produces significantly older ages, but similar masses, thus not affecting the results presented in this study.

PROSPECTOR allows us to test single stellar population SFH. In this case, we find that the age of A.1 is slightly younger (but within uncertainties) than those produced by the BAGPIPES, resulting in lower masses. This would result in slightly lower intrinsic mass M*,int = 0.39 × 106M⊙, log(Σ*) = 4.5M⊙ pc−2, and log(Π) = 1.2, but leaving unchanged any of the conclusions of this study.

Finally, given that some of the knots in the Cosmic Gems arc are unresolved (C.1 and B.2) or only marginally resolved (A.1) in our current images, we have also explored scenarios in which these sources are individual, highly magnified stars. Using SED models for stars at high redshifts52 we find that, although the slopes of the SEDs of these sources would be broadly consistent with individual stars at effective temperatures greater than about 20,000 K, these scenarios would require magnifications well in excess of what our macrolens models predict at the positions of these sources. Even the most massive and luminous stars (initial mass 560–575M⊙) described by the stellar evolutionary tracks of ref. 53 would require magnifications μ > 1,000 to explain the observed fluxes of C.1, B.2 or A.1.

Lens models and uncertainties on the derived star cluster physical properties

Four different lensing models have been created for the SPT-CL J0615−5746 cosmological field. The models are presented in detail in ref. 1. We include here below a short description.

Lenstool-A, here used as a reference model for the analysis presented in this Article, is based on the software LENSTOOL54, which uses a parametric approach and MCMC sampling of the parameter space to identify the best-fit model and uncertainties. In Lenstool-A, we model the cluster lens as a combination of three main halos and cluster member galaxies, all parameterized as pseudo-isothermal mass distributions. The model uses as constraints the positions of 43 multiple images of 14 clumps, belonging to 9 unique source galaxies. The redshifts of three sources are used as constraints (the z = 10.2 arc, and sources at z = 1.358 and z = 4.013; ref. 55), whereas the rest of the redshifts are treated as free parameters. Three clumps on each side of the main arc were used as constraints, A, B and C, assumed to be at z = 10.2. The model predicts a counterimage at (right ascension (R.A.), declination (Decl.)) = (93.9490607, −57.7701814). A possible candidate of this counterimage, observed near this location (about 1.8″), was not used as a constraint. The image plane RMS of the best-fit model is 0.36″. All the observed lensed features are well reproduced by this model.

The second model, here referred to as Lenstool-B, uses the same algorithm, but with noticeably different assumptions. This model uses 43 multiple images from 11 unique sources. A secondary galaxy cluster scale halo is placed around the location of dusty galaxies nearly 50 arcsecs north of the bright centre galaxy and allowed to move within a 20″ box around this position. As with the previous model, the z = 10.2 arc has a predicted counterimage near the possible candidate and is only about 2″ away from the Lenstool-A model. The main differences between these models are the assumptions about the mass distribution of the lens and the addition of constraints. The image plane RMS of the best-fit model is 0.68″.

The third model used in this analysis has been created with Glafic. The Glafic56,57 mass model is constructed with three elliptical NFW58 halos, external shear and cluster member galaxies modelled by pseudo-Jaffe profile. The model parameters are fitted to reproduce the position of 44 multiple images generated from 15 background sources. Spectroscopic redshifts are available for 7 of the 15 sources. We include positions of A.1/A.2 and B.1/B.2 in the Cosmic Gems arc as constraints, with small positional errors of 0.04″ to accurately predict the magnifications of each star cluster image. For the other multiple images, we adopt the positional error of 0.4″. Our best-fitting model reproduces all the multiple image positions with the RMS of image positions of 0.41″.

As a consistency check, we excluded the positional constraints from the Cosmic Gems arc to construct the mass model and confirmed that the critical curve of this mass model still passes through the arc. Our Glafic best-fitting mass model also predicts a counterimage of the Cosmic Gems arc at around (R.A., Dec.) = (93.9504865, −57.7696559). We find that there is a candidate counterimage at around 2″ from the predicted position, (R.A., Dec.) = (93.9500002, −57.7702197). Both the consistency check and the presence of the candidate counterimage confirm the validity of this mass model.

A fourth model has also been produced with WSLAP+59,60. The WSLAP+ lens models offer an alternative to parametric models and are free of assumptions made about the distribution of dark matter. When the z = 10.2 arc is not included as a constraint, the WSLAP+ model predicts the critical curve passing at about 1″ from the z = 10.2 arc. This solution predicts a mirrored image of the arc that is not observed, reinforcing the expectation that the Cosmic Gems is a double image with the critical curve passing through the middle. When the arc is included as a constraint, the predicted critical curve passes between C.1 and D.1, just 0.3″ from the alleged symmetry point in the arc and within the uncertainties typical of WSLAP+ models. Moreover, this model predicts the position of a third counterimage consistent with the previous models. This model is currently under development with the goal of explaining the perturbation seen in Img.2 and is therefore not included in this analysis.

The photometric redshift of the candidate counter image is \({z}_{{\rm{phot}}}=10\,.\,{8}_{-1.4}^{+0.6}\) (95% confidence)1, in agreement with the expectation.

The total and tangential magnifications at the position of the star clusters, μtot = μtang × μrad, are reported in Extended Data Table 3. For the two Lenstool-based models, we estimated uncertainties following the method presented in ref. 44 based on magnification maps produced from the Lenstool MCMC posterior distributions of the lens model. Uncertainties are omitted for the Glafic model.

In Extended Data Fig. 4, we show the impact that magnification predictions have in the recovered physical properties (intrinsic half-light radius and mass, black and blue solid lines) and derived quantities (dynamical age and stellar surface density in magenta and orange solid lines). We use logarithmic scales so that all quantities can be included. The coloured bands show the level of uncertainties recovered from the analysis. We also include the upper limits on the Reff (assuming that the source is unresolved and has a size smaller than the stellar PSF) and what type of lower limits it will translate for the physical quantities that depend on the size estimates as dashed lines (notice that these are the reference quantities for C.1 that is unresolved). The magnifications (total in the bottom image and tangential in the top image) are reported on the x-axis. As we move to lower magnifications, the derived masses and radii become larger, consequently predicting lower stellar surface densities and dynamical ages. However, even in the unlikely case that the magnifications are wrong by one order of magnitude, the stellar surface density will remain above 104M⊙ pc−2 and dynamical ages will still be significantly larger than 1 (log(Π) > 0), leaving the main conclusion of this analysis unchanged: we are detecting bound proto-globular clusters within the first 500 Myr of our Universe.