Properties of the brightest young stellar clumps in extremely lensed galaxies at redshifts 4 to 5

We study the populations of stellar clumps in three high-redshift galaxies, at z=4.92, 4.88 and 4.03, gravitationally lensed by the foreground galaxy clusters MS1358, RCS0224 and MACS0940, respectively. The lensed galaxies consist of multiple counter-images with large magnifications, mostly above $\rm \mu>5$ and in some cases reaching $\rm \mu>20$. We use rest-frame UV observations from the HST to extract and analyse their clump populations, counting 10, 3 and 11 unique sources, respectively. Most of the clumps have derived effective radii in the range $\rm R_{eff}=10-100$ pc, with the smallest one down to 6 pc, i.e. consistent with the sizes of individual stellar clusters. Their UV magnitudes correspond to $\rm SFR_{UV}$ mostly in the range $\rm 0.1-1\ M_\odot yr^{-1}$; the most extreme ones, reaching $\rm SFR_{UV}=5\ M_\odot yr^{-1}$ are among the UV-brightest compact ($\rm R_{eff}<100$ pc) star-forming regions observed at any redshift. Clump masses span a broad range, from $10^6$ to $\rm 10^9\ M_\odot$; stellar mass surface densities are comparable, and in many cases larger, than the ones of local stellar clusters, while being typically 10 times larger in size. By compiling published properties of clump populations at similar spatial resolution between redshift 0 and 5, we find a tentative evolution of $\rm \Sigma_{SFR}$ and $\rm \Sigma_{M_\star}$ with redshift, especially when very compact clumps ($\rm R_{eff}\leqslant20$ pc) are considered. We suggest that these trends with redshift reflect the changes in the host galaxy environments where clumps form. Comparisons with the local universe clumps/star clusters shows that, although rare, conditions for elevated clump $\rm \Sigma_{SFR}$ and $\rm \Sigma_{M_\star}$ can be found.


INTRODUCTION
Since the first deep observations with the Hubble space telescope (HST), galaxy morphology was recognised to change from disk-like or elliptical into more irregular appearance at increasing redshifts (e.g.Abraham et al. 1996;Brinchmann et al. 1998).In addition, galaxies around the cosmic noon (z∼1-3) are characterised by the presence of bright stellar clumps dominating their rest-frame ultraviolet (UV) morphology (e.g.Cowie et al. 1995;van den Bergh et al. 1996).The James Webb Space Telescope (JWST) is bringing new insight into the properties of high-z galaxies, especially at the epoch of re-ionisation (z ≳ 7); the first results seem to confirm the overall morphological evolution traced by HST at lower redshifts, with galaxies at redshift 7-12 characterised by irregular (yet compact) structures and in a minor part (∼ 20%) by interaction/mergers (Treu et al. 2023).
Most of clumps at redshift z < 3 may have formed in-situ within their host galaxies.This scenario is supported by several observational evidence, such as: (1) the redshift evolution of clumpy galaxies, closely following evolution of the overall star formation rate (SFR) volume density and inconsistent with the evolutionary trends of minor and major mergers (Lotz et al. 2011;Shibuya et al. 2016); (2) the presence of numerous clumpy galaxies (at least up to z∼ 3) still dominated by disk-like appearance (Shibuya et al. 2016), with comparable disk scale-heights in case of either presence or absence of clumps (Elmegreen & Elmegreen 2006;Elmegreen et al. 2017); (3) the kinematics of the majority of clumpy star-forming galaxies at cosmic noon being dominated by ordered disk rotation (yet with elevated velocity dispersions, Förster Schreiber et al. 2009;Wisnioski et al. 2015;Rodrigues et al. 2016;Swinbank et al. 2017;Turner et al. 2017;Simons et al. 2017;Girard et al. 2018).Simulations show that such turbulent high-z disks fragment because of gravitational instability and can form gas clouds that turn into stellar clumps (Bournaud et al. 2014;Tamburello et al. 2015;Mandelker et al. 2017); (4) the observations of very young clumps (age ≲ 10 Myr), inconsistent with an "external" origin (e.g.Förster Schreiber et al. 2011b;Zanella et al. 2015); (5) finally, the stellar mass function of these clumps follows a power-law with slope −2 (Dessauges-Zavadsky & Adamo 2018), characteristic of nearby star clusters and HII regions (see e.g. the reviews by Krumholz et al. 2019;Adamo et al. 2020a) and expected for stellar structures formed in a hierarchical manner via turbulencedriven fragmentation (e.g.Elmegreen 2010;Guszejnov et al. 2018;Ma et al. 2020).In this scenario, UV-bright clumps are simply star cluster complexes formed in high-z galaxies, and thus likely trace the star-formation process in their host galaxy.
With respect to local stellar clusters and cluster complexes, highredshift clumps are more extreme systems, usually with elevated star formation rate (SFR) and SFR densities (e.g.Livermore et al. 2015;Messa et al. 2022;Claeyssens et al. 2023), and sometimes observed as mini-starburst entities within their host galaxies (Zanella et al. 2015;Iani et al. 2022).Within the in-situ formation scenario, this difference can be explained by high-z disks fragmenting at much larger scales (and possibly densities) than in local MS galaxies because of their gas-rich and turbulent nature, as suggested by models (e.g.Immeli et al. 2004), numerical simulations of turbulent high-redshift galaxies (e.g.Tamburello et al. 2015;Renaud et al. 2021;van Donkelaar et al. 2022), observations of dense giant molecular cloud complexes from CO data in galaxies at z = 1 (Dessauges- Zavadsky et al. 2019Zavadsky et al. , 2023)), as well as by observations in nearby analogs (e.g.Fisher et al. 2017a,b;Messa et al. 2019).
While the in-situ origin seems to explain the formation of ∼ 70% of the clumps at redshifts z < 3 (Zanella et al. 2019), it is still possible that galaxies at earlier times were characterised by clumps formed ex-situ, i.e. by a merger process during which the satellite is stripped and its nucleus appears as a massive clump (e.g.Somerville et al. 2000;Hopkins et al. 2008;Puech 2010;Straughn et al. 2015;Ribeiro et al. 2017).This scenario could be justified by the increasing galaxy (minor and major) merger rate towards higher redshifts (e.g.Lotz et al. 2011), as well as by the increasingly lower number of massive galaxies (Marchesini et al. 2009); simulations show that large-enough galaxy masses are needed for disk fragmentation to happen (e.g.Tamburello et al. 2015).Unfortunately systems before cosmic noon ( ≳ 4), where violent disk instability is thought to have less impact on clump formation, are currently under-represented in clump studies (e.g.Meštrić et al. 2022).The JWST will soon bring a new insight on clump formation at these redshifts, by providing deep, high-resolution rest-frame optical observations (as seen from the first results by e.g.Mowla et al. 2022;Claeyssens et al. 2022;Vanzella et al. 2022a).
In this paper, we analyse the clump populations of three bright and highly magnified gravitationally-lensed galaxies at z > 3.5, in order to pave the way for upcoming JWST studies of larger samples.The selection, main properties and observational data of the sample are presented in Sect.2; the clump analysis methodology is described in Sect.3; results are presented in Sect. 4 and discussed, also in the context of other literature works, in Sect.5; finally, the main details of this analysis are summarised in Sect.6.Throughout this paper, we adopt a flat Λ-CDM cosmology with  0 = 68 km s −1 Mpc −1 and Ω M = 0.31 (Planck Collaboration et al. 2014), and the Kroupa (2001) initial mass function.All quoted magnitudes are on the AB system.

Galaxy Sample
We search for spatially extended lensed galaxies at  > 3.5, with spectroscopically-confirmed redshifts.These criteria are chosen to characterise sub-galactic scales at redshifts where, as outlined in the introduction, we may expect different formation mechanisms than cosmic noon.The search is restricted to cluster lensing, to allow for larger magnifications across the full extent of the arcs.The sample of clusters is taken from the full MAssive Clusters Survey (MACS) survey (e.g.Ebeling et al. (2010); Repp et al. (2016)), as well as other clusters taken from the LoCuSS (Richard et al. 2010), CLASH (Postman et al. 2012), Frontier Fields (Lotz et al. 2017) and RELICS (Coe et al. 2019) cluster samples.In Fig. 1 we observe how, while the typical population of background lensed galaxies behind many massive galaxy clusters is distributed around apparent (and magnified) rest-frame UV magnitudes in the range 26 − 30 mag (see also Richard et al. 2021, Claeyssens et al. 2022), three galaxies clearly stand out in terms of brightness and angular area.Those galaxies are located beyond the lensing clusters MS 1358+62 (at z cl = 0.33), RCS 0224-0002 (at z cl = 0.77) and MACS J0940.9+0744(at z cl = 0.34); the lensed systems are in the redshift range z = 4 − 5 and are all magnified by large factors ( > 5).The selected targets are typical of their redshift, in terms of masses and SFRs (see Sections 2.1.1,2.1.2 and 2.1.3)but appear 'clumpy', i.e. host more than one bright clump (in Fig. 1, the marker size is proportional to the number of observed clumps), due to their large magnification factors, making them optimal candidates for the study of star formation at small subgalactic scales, down to 10 pc (see Section 4.1).In addition, their large magnified observed brightness makes them ideal targets for future JWST observations 1 .
We summarise below the main properties of these galaxies, from previous studies in the literature; for simplicity, we will use the shortened names of the galaxy clusters as the names of the high-z lensed galaxies we analyse in this work.

MS1358
The  = 4.92 lensed galaxy behind MS 1358+62 galaxy cluster was first discovered and studied by Franx et al. (1997).The galaxy rest-frame UV/optical morphology is dominated by several compact star-forming regions.Fitting the broad-band SED from HST and Spitzer imaging, Swinbank et al. (2009) derived a total mass M ★ = (4.7 ± 1.3) × 10 8 M ⊙ 2 , from a young stellar population (14 ± 7 Myr) with sub-solar metallicity (Z = 0.2 Z ⊙ ); the stellar extinction is consistent with being close to zero, E(B − V) = 0.05 ± 0.05 mag.A star formation rate of SFR = 28 ± 5 M ⊙ yr −1 was derived from 1 JWST observations with NIRSpec-IFU and NIRCam imaging have been approved for these targets, GO program 3433. 2 The original values of mass (M ★ = 7 × 10 8 M ⊙ ) and star formation rate (SFR = 42 M ⊙ yr −1 for the entire galaxy and SFR = 18 M ⊙ yr −1 for the two brightest clumps) were derived by Swinbank et al. (2009) assuming a Salpeter (1955) IMF and are here converted to match the assumption of Kroupa (2001) IMF used throughout the current paper.the [OII] emission-line flux (Swinbank et al. 2009).The galaxy is characterised by the presence of two main sub-galactic star-forming regions (IDs: 1 and 2 in Fig. 2) accounting for 12 ± 1 M ⊙ yr −1 in SFR.IRAM PdBI observations in CO  emission suggest that the galaxy has a large gas fraction, f gas = 0.6 (Livermore et al. 2012b), similarly to what is observed in z ≥ 3 galaxies (Dessauges-Zavadsky et al. 2020).

RCS0224
A lensing-magnified arc at  = 4.88 in the RCS 0224-0002 cluster, it was first discovered as a bright Ly halo by Gladders et al. (2002).Swinbank et al. (2007) studied the rest-frame ultraviolet and optical properties of the galaxy by combining HST imaging to VIMOS and SINFONI spectroscopy, deriving a dynamical mass of ∼ 10 10 M ⊙ (within the estimated 2 kpc intrinsic size of the galaxy) and a star formation rate of 8.2 ± 1.4 M ⊙ yr −1 (using [OII] observations3 covering only the central and western, C and W, images in Fig. 3).A further study with MUSE suggested that the detected emission lines are powered by a young (< 5 Myr) and low-metallicity (Z ≲ 0.05 Z ⊙ ) stellar population (Smit et al. 2017).

MACS0940
MACS0940 at z = 4.03, first studied by Leethochawalit et al. (2016), is observed as a strongly stretched arc, made of two mirrored images (see Fig. 4).Two complete counter-images of the same galaxy are seen northern and eastern from the arc.The galaxy is characterised by a bright Ly halo (Claeyssens et al. 2019), suggesting recent episodes of star formation.Bright compact sources are clearly visible along the arc in rest-frame UV observations.

Hubble Space Telescope (HST) observations
HST observations with WFPC2, ACS/WFC and WFC3/IR are available on the HST MAST archive (see Data Availability section at the end of the publication for details about program IDs).The list of filters used is given in Tab. 1.For each galaxy, we choose a reference filter, covering the rest-frame UV wavelengths (see Tab. 2), where we measure the size and UV luminosity of the clumps (as described in Section 3.2); the other filters are used to infer the broad-band SED of the clumps (see Section 3.4).Individual flat-fielded and CTEcorrected exposures were aligned and combined in single images using the AstroDrizzle procedure from the DrizzlePac package (Hoffmann et al. 2021).The astrometry was aligned to the Gaia DR2 (Gaia Collaboration et al. 2018).
We model the instrumental point-spread function (PSF) from isolated bright stars within the field of view of the observations.In each filter, we fit the selected stars with an analytical function described by a combination of Moffat, Gaussian, and 4 th degree power-law profiles.Such a generic combination is chosen to mitigate possible bias introduced by the choice of a specific function.The fits provide good descriptions of the stars up to a radius of ∼ 0.8 ′′ , i.e. more than 10× larger than the half-width at half maximum (HWHM) of PSF in the reference filters (the full width at half maximum, FWHM, of the reference filters are 0.10 ′′ , 0.11 ′′ and 0.12 ′′ for MS1358, RCS0224 and MACS0940, respectively, see Tab. 2).The modeled PSF is used to derive the size and luminosity of the clumps, as described in detail in Section 3.2.

Gravitational lens models
Our study relies on the reconstruction of intrinsic properties of the identified clumps, accounting for the anisotropic lensing magnification which is estimated from a lens model.
As a starting point, we make use of existing mass models of the cluster cores (references in Table 1) which include the lensed arcs as constraints.In summary, these models use a parametric mass distribution describing the cluster-scale and galaxy-scale mass components of the clusters as a combination of double Pseudo Isothermal Elliptical (dPIE) profiles.The models are optimised with the Lenstool Jullo et al. (2007) software 4 based on multiple imaged constraints, similarly to the modelling performed in e.g.Richard et al. (2014).In the case of the older models for MS1358 and RCS0224, we have performed a new optimisation using the latest version (v8) of Lenstool.
In all cases, individual clumps identified as mirrored images in the three lensed arcs are individually used as constraints in the model, ensuring a very good reproduction of the images.The rms between the observed and predicted locations of the images used in the model is 0.06", 0.08" and 0.24" for MS1358, RCS0224 and MACS0940 respectively.
The outcome of the Lenstool optimisation is a statistical sample of mass models sampling the posterior distribution function of the model parameters.We make use of both the best model (achieving the lowest rms) and the range of mass models to estimate the magnification factors (along each direction) at each position across the arc and their 68% percentile error.The range of magnification factors obtained is summarised in Table 1.   1. Galaxy sample and HST data: (1) ID of the galaxy (we remind that we are using the name of the hosting galaxy clusters as names/IDs of the lensed arcs studied); (2) redshift of the galaxies, from Franx et al. (1997), Gladders et al. (2002) and Leethochawalit et al. (2016); (3) number of lensed images analysed (in parenthesis the number of counter-images not included in the clump study); ( 4) range of magnification values covered by the galaxy; (5) references for the lens models; (6) reference filter used to extract the size and UV magnitude of the clumps; (7) other filters, used to derive the broadband SED of the clumps.† All filters are from either WFC3 or ACS, except F606W from WFPC2.‡ F775W is used instead of F814W as the reference filter for MS1358 because of its longer exposure time, providing better signal-to-noise.2, for MACS0940 (r:F160W, g:F125W, b:F814W, zoom-ins:F814W).According to the lens model of the system, the SE and SW arcs (and all the clumps they host) are different lensed images of the main galaxy (labeled as "0" in the N and E fields).The clumps 9, 10 and 11, observed in the N and E fields, do not have counterparts in the SE and SW fields.

Clump extraction
Clump candidates have been extracted by running SExtractor (Bertin & Arnouts 1996) requiring 3 detections in at least 4 pixels (px) in background subtracted images (where the background is evaluated on a 30 px scale) in the reference filter.In order to check for 'redder' sources that could be missed in the reference filter, we ran SExtractor with the same set of parameters also in all the other available filters in each galaxy; this test did not produce any further source.The colors of the clumps and the lens models have been used to identify interlopers, mainly bright compact sources at different redshifts, among the extracted sources (some examples are highlighted in Figs. 2 and 4).6) extraction and completeness limits for a pointlike source within the galaxy, as described in Section 3.2.1;note that these values are corrected for the Galactic reddening.
overall.Only 6 of those clumps are seen in the western arc5 (SW), while both the eastern and northern counter-images (E and N) look like a single bright source (consistent with the superposition of the brightest clumps, IDs: 1, 2, 3 and 4) surrounded by diffuse light and by 3 sources (IDs: 9, 10 and 11, not seen in the arcs).

Clump modeling
Clumps are modelled on the image plane, fitting their sizes and fluxes on the observed data of the reference filter; such quantities are later converted into intrinsic values, using the magnification maps produced by the lens model (Section 2.3).The clump modelling follows accurately the methodology applied to the study of stellar clumps in the z=1 galaxy A521-sys1 (Section 3.2 of Messa et al. 2022) of which we summarize here the main features.We assume that clumps have Gaussian profiles in the image plane, and therefore we fit a model consisting of 2D Gaussian functions, convolved to the the instrumental point spread function (i.e. the response of the instrument), to obtain their observed profiles.We performed clump fitting in 9 × 9 px cutouts centered on each of the clumps, where a 1st degree polynomial function is added to account for the galaxy background luminosity: both size (including axis-ratio and orientation, for elliptical sources) and flux of the clumps are left as free parameters of the fit in the reference filter.The best-fit models and residuals are shown in Appendix A. The clump shape is then kept fixed in the other filters, where only the source flux is fitted, under the assumption that the intrinsic clump shape is the same in all bands.In order to reduce the possible contamination from nearby bright sources, the brightest clumps were fitted first, and then their best-fit model subtracted from the data.The fitted fluxes, in units of  − /, are converted into observed AB magnitudes by considering the instrumental zeropoints and by subtracting the reddening introduced by the Milky Way in each of the filters.The (observed) effective radius of the clumps, R eff,obs , here defined as the radius enclosing half of the source's luminosity, is equal to the (circularized) full width at half maximum, in the current assumption of gaussian profiles.In the occurrence of multiple clumps separated by less than 5 px, a single fit is performed in cutouts large enough to include all sources (typically 15 × 15 px).
As already discussed in Messa et al. (2019) and Messa et al. (2022), the oversampling of the HST PSF (FWHM∼ 2 px in the filters considered in this work) allow to resolve sub-pixel clump sizes.By inserting mock clumps of various sizes in the image frames and fitting them as for the real clumps, we derive the minimum resolvable size,  min = 0.40 px for MS1358 and MACS0940, and  min = 0.45 px for RCS0224.These limits, converted into physical sizes at the redshift of the galaxies, are highlighted in the top panels of Fig. 5.We refer to Appendix C of Messa et al. 2022 for more details on the process of estimating the minimum resolvable size.

Completeness of the sample
We separately tested the completeness we reach in extracting and analysing our samples.In order to test the sensitivity of our extraction process we insert, one by one, 1500 synthetic clumps at random positions within the region covered by the lensed systems (in the reference filter used in the original extraction, Section 3.1).The clumps are modelled as symmetric gaussians and divided into 3 'size groups', with  = 0.4 (0.45 for RCS0224), 1.0 and 2.0 px; the first value correspond to the minimum resolvable size (see Section 3.2), while the last encompass the largest clumps observed in the samples (see Fig. 5, top row).Clumps are randomly drawn from a flat distribution in magnitudes.For each realisation we run SExtractor with the same sets of parameters used in Section 3.1 and we check if the synthetic clump has been extracted.For each of the three sizes chosen we consider as extraction limit, lim ext the magnitude below which the fraction of extracted clumps goes above 90%.
In order to test the reliability of the derived photometry, we estimate the completeness of the sample in a second way, following the method described in Appendix D of Messa et al. (2022); the same synthetic clumps used to estimate the extraction limits are photometrically analyzed in the same way as the real sources.We consider as good sources the ones where the relative error on the recovered flux is below 50%, flux rel = |flux in − flux out |/flux in < 0.5.We consider as completeness limits, lim com , the magnitudes below which the fraction of good sources recovered goes above 80%.All the extraction and completeness limits are compared to the photometry of the clumps in the sample in Fig. 5 (top row).The limits in the case of point-like clumps are also reported in Tab 2.

Conversion to intrinsic sizes and magnitudes
In order to recover the intrinsic clumps' properties, we considered the magnification map produced by the best fit lens model of each galaxy; we used as reference amplification,  tot , the median amplification value found in the region centered on each clump coordinates and extending within one FWHM of its size.We use the standard deviation of amplification values found in the same region as an estimate of the magnification uncertainty associated to the clump position.To estimate the uncertainty associated to the magnification map, we consider, for each cluster, the 1 interval (16 th to 84 th percentiles) of the amplification values found in the 500 lens model variations described in Section 2.3.The two magnification uncertainties are combined (by the sum of the squares) into a final uncertainty.We note that the uncertainty from the 500 'variation' maps is the one usually dominating the final value.
The observed magnitudes are converted into absolute ones by subtracting the distance modulus and adding the  correction, a factor 2.5 log(1 + z); the amplification is accounted for by adding 2.5 log( tot ), converting the magnitude into the intrinsic value.
We consider three cases when measuring the intrinsic effective radius, R eff , following the methodology already used and discussed in the literature (e.g.Vanzella et al. 2017a;Claeyssens et al. 2023).If the clump is resolved along both minor and major axis in the image plan, we simply divide the observed R eff,obs by the squareroot of the clump amplification.If the clump is un-resolved along the transversal direction of the magnification, we consider as informative only the size measured in the shear direction, and we divide the latter by the tangential component of the magnification ( tan ).Finally, if the clump is un-resolved in both directions we divide the observed R eff,obs (consistent with the instrumental PSF) by  tan to derive a size upper limit.In the last two cases, the underlying assumption is that the clump has approximately a circular shape in the source plane.The final uncertainties on the intrinsic properties combine both photometric and magnification uncertainties via the root sum squared.

Broadband SED fitting
We use the broadband photometry to estimate the clumps' masses.Following the methodology of Messa et al. (2022) and Claeyssens et al. (2023), photometry in the filters other than the reference one is performed by assuming that each source has the same observed size (derived in Section 3.2) in all the bands, and thus fitting only the flux and the background.Due to the elevated redshift of the sources, the available HST filters cover only the rest-frame wavelength range ∼ 1000−3000 Å.As a consequence our SED-derived values are UVweighted quantities.In particular, ages and extinctions are only poorly constrained; on the other hand, mass estimated are more robust, and within a factor of few correct (see also Messa et al. 2022 for the robustness of mass estimates over different model assumptions).
To mitigate the effects of degeneracies between parameters of the fit (in particular ages, extinctions, metallicities and star formation histories, SFHs), we limit the number of free parameters.First of all, we fixed the metallicity of the stellar models; following Stark et al. (2009) we use a sub-solar metallicity, Z = 0.2 Z ⊙ .Smit et al. (2017) suggest that the nebular CIV lines observed in RCS0224 must come from a young stellar population with metallicity Z = 0.05 Z ⊙ or lower.For this reason we perform, for all galaxies, a second 'control' fit using stellar models with metallicity Z = 0.02 Z ⊙ .As second strong assumption, we consider SFHs described by a 10 Myr continuous star formation.This choice is driven by the sizes of the clumps (Section 4.1), larger than typical stellar clusters (for which an instantaneous burst is usually assumed); assuming longer histories, like a 100 Myr continuous star formation, would lead to larger masses on average (as already found by, e.g.Adamo et al. 2013;Messa et al. 2022;Claeyssens et al. 2023).We use the stellar models from the Yggdrasil stellar population synthesis code6 (Zackrisson et al. 2011), based on Starburst99 Padova-AGB tracks (Leitherer et al. 1999;Vázquez & Leitherer 2005) with a universal Kroupa (2001) IMF (in the mass interval 0.1 − 100 M ⊙ ), processed through Cloudy (Ferland et al. 2013), assuming a 50% nebular covering fraction, to obtain the evolution of the nebular continuum and line emission produced by the ionized gas.The model spectra, at each age, are attenuated with a color excess ranging between E(B − V) = 0 − 1 mag, using the Calzetti et al. (2000) law, before being convoluted with the filter throughput.We do not use the differential expression formulated by Calzetti et al. (2000), but we apply the same reddening to both stellar and nebular emission, assuming that stars and gas are well mixed.We check a-posteriori that the accepted solutions for the SED fit have low extinctions, E(B − V) ⩽ 0.4 mag in most cases.Age, mass and extinction are left as free-parameters in the SED fitting process.Best-fit parameters are given by the model with the lowest reduced  2 ( 2 red.,min ).Their uncertainties are given by the entire range of solutions whose reduced  2 satisfies  2 red.≤ 2 ×  2 red.,min ; this value was chosen, by inspecting the fit results, as it encompasses similarly-good solutions.Like it was the case for magnitudes in Section 3.3, de-lensed masses are derived by dividing the observed mass of each clump by its total magnification.

Clumps sizes and luminosities
The observed (i.e.not de-lensed, and therefore not 'intrinsic') clump sizes and magnitudes are shown in the top row of Fig. 5, where they are compared to the extraction and completeness limits described in Section 3.2.1.The magnitude ranges are similar among the different galaxies (∼ 28 − 25 AB mag) and, in the great majority of cases, the analyzed clumps are above the completeness limit, suggesting that their photometry is robust.The shallower extraction and completeness limits in RCS0224 (consistent with the exposure time of this galaxy being the shortest, see Tab 2) causes the average clump magnitude in this galaxy to be bright (∼ 26 − 25 mag).This result may suggest that we are missing clumps at lower surface brightness.A shaded-gray area in the plot marks the region below the size lower limits discussed in Section 3.2; in the absence of gravitational lensing we wouldn't be able to study clumps below ∼ 100 − 200 pc scales.
The intrinsic (i.e.de-lensed) clump sizes and magnitudes are shown in the bottom row of Fig. 5 and are reported in Tab. 3. Intrinsic sizes are smaller than 200 pc and in many cases reach values close to 10 pc; in particular, due to the magnification factors associated to our galaxies, we are studying clumps at scales comparable to large individual star clusters in RCS0224 (R eff = 5 − 25 pc) and slightly larger scales in the other two cases, consistent with the sizes of compact star-forming regions (R eff = 13 − 200 pc in MACS0940 and R eff = 30 − 150 pc in MS1358).The rest-frame UV magnitudes, on the y-axis of Fig. 5 (bottom panels), have also been converted into a SFR value (using the conversion factors from Kennicutt & Evans 2012), as this is a commonly used parameter in high-z clump studies.In the same panels we show lines of uniform SFR surface densities (Σ SFR = 1, 10, 10 2 , 10 3 M ⊙ yr −1 kpc −2 ).The extraction limits described above, when translated into surface-brightness limits (dashed lines), correspond to Σ SFR ∼ 10 M ⊙ yr −1 kpc −2 .In all cases the completeness is above the typical Σ SFR values of clumps in local main sequence galaxies ∼ 0.6 M ⊙ yr −1 kpc −2 (Kennicutt et al. 2003;Livermore et al. 2015).In some cases, we are not able to robustly constrain the clump properties, mainly due to very large uncertainties on the magnification.Those cases are discussed for each galaxy individually in Section 5.1 and are then removed from following analyses.A detailed discussion of clump intrinsic properties for each of the galaxies is given in Section 5.1, while a comparison of clump sizes and SFRs with other literature samples is given in Section 5.2.

Clumps masses
We report in Tab. 3 the properties derived from the broadband SED fitting of the clumps.Due to the large uncertainties on ages and extinctions, we choose to provide only a range of allowed values for those two degenerate properties.The majority of clumps are consistent with being as young as ∼ 1 Myr (this is consistent with the observations of bright nebular emission in Ly, CIV and [OII] reported in literature, e.g.Smit et al. 2017;Claeyssens et al. 2019) but with equally-probable solutions at ages as old as ∼ 100 Myr.
In the latter case, a star formation longer than what is considered should be accounted for, in order to explain the nebular emission.
The derived color excesses mainly span the range E(B-V)=0.0-0.3 mag (up to 0.6 mag in MS1358), suggesting low overall extinctions in these galaxies.
Model degeneracies have lower effect on clump masses, which we provide as best-fit values; intrinsic masses are affected by uncertainties on the lens model, as seen for absolute magnitudes in the previous section.Derived clump masses span the range M = 10 6 − 10 9 M ⊙ .Using models with lower metallicity (Z = 0.02 Z ⊙ ) than the reference one (Z = 0.2 Z ⊙ ) we get, on average, the same masses (within 0.1 dex) as the ones reported in Tab. 3 for MS1358 and MACS0940.In the case of RCS0224, we get 0.3 dex larger masses associated to older best-fit ages, on average.We also expect that, using models with longer SFHs, we would derive, on average, larger masses (e.g. by ∼ 0.1 dex for 100 Myr continuous star formation, see Messa et al. 2022).Both these alternative assumptions would cause mass differences that are within the mass uncertainty ranges reported in Tab. 3.
The comparison of clump masses to their sizes, in Fig. 6, reveals that we are looking at very dense systems; they have, on average, mass surface densities Σ M ★ ∼ 10 3 M ⊙ pc −2 , similar to those of stellar clusters in local galaxies (e.g. Brown & Gnedin 2021) but on scales which are up to ∼ 10 times larger, i.e. up to ∼ 100 pc in the case of MS1358.The densest systems observed in these three high-z galaxies reach values > 10 4 M ⊙ pc −2 , matching the most extreme stellar ensembles observed at any redshift (see the discussion in Section 5.2).
The combination of mass (M ★ ) and size (R eff ) of the clumps provides an estimate of their crossing time, defined as (Gieles & Portegies Zwart 2011): In local studies of stellar clusters, crossing times are compared to cluster ages to derive the so-called 'dynamical age' of clusters (Π ≡ Age/T cr ); a value Π > 1 indicates that the stars in the system remained clustered together, without freely expanding into their surroundings, for a time longer than their crossing time, implying that the system is likely gravitationally bound (e.g.Gieles & Portegies Zwart 2011;Ryon et al. 2015Ryon et al. , 2017;;Krumholz et al. 2019;Brown & Gnedin 2021).The same kind of analysis has been recently applied to the study of high-z stellar clumps (e.g.Vanzella et al. 2021;Messa et al. 2022;Claeyssens et al. 2023;Vanzella et al. 2022a).Crossing times for the clumps in the three galaxies of the current sample are reported in Tab. 3. Given the large age uncertainties, we do not attempt to calculate the respective dynamical ages; however, we notice that T cr ≲ 10 Myr for ∼ 55% of the cases, suggesting that (even if  their young ages are confirmed) a large fraction of the clumps we are observing can be of gravitationally bound systems.

Stellar clump populations
We discuss in this section the results presented in Fig. 5 and Fig. 6, and collected in Tab. 3, by putting the clump properties in the context of their host, for each of the galaxies studied.

The clump population of MS1358
MS1358 is characterised by the presence of several clumps; 10 were extracted in the western image (the one with the largest overall magnification), 8 of which are seen also in the eastern image.The analyses on the two images of the galaxy, produce comparable clump properties (see Tab. 3; Figs. 5 and 6).For clump 8, visible only in the western part of the galaxy, we were not able to derive a robust size estimate.
As suggested by previous studies, and as clearly visible from Fig. 2, MS1358 is dominated, in the rest-frame UV, by two main clumps (ID: 1 and 2); those account for ∼ 25% of the flux in F775W, rest-frame ∼ 1300 Å (clump 1 alone accounts for ∼ 20%, while if all clumps are considered, their contribution to the rest-frame UV emission of the galaxy is ∼ 40%), suggesting that they are major contributors to the recent star formation of the galaxy.Their rest-frame UV fluxes, converted to the SFR via the Kennicutt & Evans (2012) relation, correspond to 5 and 1 M ⊙ yr −1 .Nebular tracers of the current star formation indicate that their contribution is even larger (∼ 40% of the SFR of the galaxy, Swinbank et al. 2009).The [OII]-derived SFR, 7.5 and 3.7 M ⊙ yr −1 for the two clumps 7 ), suggests that the SFR UV values we derived for our clumps may be underestimated; this could either be due to the SFR episode in the galaxy being more recent than what is assumed by the Kennicutt & Evans (2012) conversion and/or due to the presence of some extinction.The derived masses span the range 10 7 − 10 9 M ⊙ and are therefore comparable to the mass of the entire galaxy as estimated by Swinbank et al. (2009).The derived young ages for some of the clumps (≲ 10 Myr), combined to the nebular emission observed in the galaxy (Swinbank et al. 2009) suggest that a considerable fraction of the stellar mass in the galaxy is being formed during the current star formation episode, and that the clumps are an important contributor to that process.

The clump population of RCS0224
RCS0224 is characterised by 3 bright clumps, accounting for 45% of the rest-frame UV emission.The brightest one (ID: 1, SFR UV = 0.7 M ⊙ yr −1 ) is not multiply imaged and therefore is only visible in the eastern image of the arc.The other two clumps do not have robust intrinsic size (and flux) estimates in the eastern image, due to the large uncertainties related to the magnification values; we can however rely on their properties derived in the western image.They are seen down to very small scales; both their sizes (below 10 pc for clump 2) and masses (10 6 − 10 7 M ⊙ ) are close to the ones of super star clusters in nearby galaxies (e.g.Leitherer et al. 2018;Vanzella et al. 2019;Adamo et al. 2020b).The mass surface density of clump 3 (Σ M ★ ∼ 10 3 M ⊙ pc −2 ) is also consistent with nearby star clusters (e.g. Brown & Gnedin 2021), while clump 2 has a density close to Σ M ★ ∼ 10 5 M ⊙ pc −2 , one of the largest observed for star-forming regions (see also the discussion in Section 5.2).A single source is visible in the central image of the galaxy (see Fig 3) but the lensing model cannot distinguish which source it is.The intrinsic size of the source is 6 pc and both the luminosity and mass suggest that this is another image of clump 2. Summing the SFR UV of the clumps in the central and western images and comparing it to the value estimated from [OII] emission in the same region (8.2 M ⊙ yr −1 , Swinbank et al. 2007), we derive that only < 10% of the galaxy SF is in the observed clumps.This result is consistent with bright [OII] emission being observed along the entire arc seen in F814W (a long part of which is devoid of compact sources, see Fig 3); it may also indicate, as pointed out for MS1358, that SFR UV values are underestimates and therefore that the bulk of star formation in the galaxy is very recent.This galaxy is also the one with the shallowest data (as discussed in Section 4.1), and therefore there is also the possibility that current observations are missing a population of low surface brightness clumps, currently undetected along the arc.

The clump population of MACS0940
Compact sources are seen all along the lensed arcs of MACS0940; the galaxy is dominated by 4 bright clumps (IDs: 1 to 4) contributing to ∼ 40% of the rest-frame UV emission of the galaxy (the contribution of all the clumps in the arcs is ∼ 50%).For some of the clumps we were not able to derive robust size estimates.SE_5 and SE_6 have large size uncertainties; contrary to the rest of the sample, they both have de-lensed magnitudes that differ largely from their SW counterpart (SW_5,6), hinting to a possible mis-association.We point out that there are no detected SW clumps with similar photometric properties to the former.Many of the clumps in the SW region (e.g.IDs: 2, 3, 4 and 7) also have large size uncertainties, mainly due to uncertainties in the lens model in that region of the arc; on the other hand, their SE counterparts have robust measurement that help inferring their intrinsic properties and will be used in the further analyses of this work, Section 5.2.For some clumps we could not perform SED fitting due to the lack of signal in some of the filters (clumps SE_5, SE_7, N_10, N_11, E_11 and E_12).Sizes and masses of the bright clumps in the galaxy, distributed mainly between 15 − 60 pc and 10 6 − 10 8 M ⊙ , lead to large mass surface densities, Σ M ★ = 10 3 − 10 4 M ⊙ pc −2 .While little is known about the properties of MACS0940, the clump contribution to the restframe UV emission and the spatial superposition of the peaks of Ly emission (Claeyssens et al. 2019) on the location of the brightest clumps suggest that the observed clumps are strongly contributing to the host galaxy recent star formation.

Size and SFR across redshifts
The sizes and luminosities of the clumps in these 3 galaxies are compared to samples from literature in Fig. 7 (top-left panel); in case of clumps with multiple images we will consider only the values recovered from the one with highest magnification.The samples are color-coded according to their redshift.The figure shows that, at any scale between ∼ 1 and ∼ 10 3 pc, clumps become on average increasingly brighter moving to higher redshifts; a similar trend was already suggested by Livermore et al. (2015) using samples studied in H emission.We point out that the SFR scale in Fig. 7 (derived from the Kennicutt & Evans 2012 conversion, as already described in Section 4.1), is used to directly compare samples studied in H (SINGS, L12, L15, DYNAMO, see figure caption) to the others, studied in rest-frame UV.We are not including, in this study, SFR values derived from the analysis of the clump spectral energy distributions 8 .The average clump SFR H surface densities for z = 0, 1, 3 and 5 proposed by Livermore et al. (2015) are shown as dashed lines in Fig. 7 (top-left panel), and are consistent with the overall densities of the plotted samples.In the same figure we show, as black-grey contours (enclosing 1 and 2 of the sample), the sizes and SFRs of HII regions in the SINGS sample of local MS galaxies (Kennicutt et al. 2003); already at z = 1, clumps are clearly detached from the region covered by the SINGS sample, as pointed out by previous studies (e.g.Livermore et al. 2012aLivermore et al. , 2015;;Messa et al. 2022;Claeyssens et al. 2023).
The redshift evolution of clump Σ SFR is made explicit in the bottom-left panel of Fig. 7 via violin distributions.For the whole sample we consider only clumps with R eff < 100 pc, as we want to focus this analysis on compact star-forming regions; however, similar results would be derived considering clumps at all scales (see Appendix B).A clear redshift evolution of Σ SFR is observed comparing local clumps (z = 0) to samples at cosmic 'afternoon' (0 < z < 1.5) and at cosmic noon (1.5 ⩽ z < 3.5).A shallower evolution is seen at earlier cosmic times: while the distributions are on average shifted towards denser Σ SFR at higher redshifts, extreme values can be observed already at  ∼ 2, as it is the case for the clump population of 8 The published sample from Meštrić et al. ( 2022) includes SED-derived SFR values; however, for consistency with the other samples included in Fig. 7, we chose to plot only their rest-frame UV magnitudes.
the Sunburst galaxy (purple triangles in the top-left panel of Fig. 7), studied at scales < 10 pc (Vanzella et al. 2022b).
A necessary caveat is that the surface brightness completeness is different for each of the samples considered.Completeness limits are unavailable for many of the samples; it can however be assumed that they become brighter moving to higher redshift.For example, the completeness limits derived in the current work, 1 − 10 M ⊙ yr −1 kpc −2 , are larger than the densest clumps at z < 1.5.We argue that the completeness limits, biasing the study of high-z systems towards the most extreme (i.e.densest) sources, may be the main driver of the Σ SFR redshift evolution overall, especially of the median values of their distributions.As mentioned in the discussion of individual galaxies (Section 5.1) other factors possibly affecting the overall Σ SFR distributions are the ages of the clumps and their extinctions.On the other hand, we expect that the high-density end of each distribution is not affected by the completeness; in this case the evolution of the densest clumps would be real and may reflect the redshift evolution of the properties of their host galaxies.We also remind that there is a great in-homogeneity in how the considered literature samples were selected and analysed.Systematic statistical studies of clumps across redshift ranges are needed to truly prove and map the evolution of clumps; forthcoming studies with the JWST will indeed provide this statistics.
We can interpret the clump Σ SFR redshift evolution discussed in this section as an evolution of the host galaxy conditions where the clumps form.Main sequence galaxies are characterised by higher SFR densities at higher redshift, with the most extreme change happening during the cosmic noon (z = 1 − 3.5, e.g.Schreiber et al. 2015).In order to test this possibility we consider, in the next section, the properties of clumps in local starburst galaxies, i.e. galaxies which do not fall into the local main sequence but rather present typical properties of higher redshift systems.

Size and SFR of nearby samples
In order to test the effects of galaxy environments on their stellar clump population, we consider local samples of galaxies characterised by properties typical of high-z systems.Galaxies in the DYNAMO sample, at z = 0.2, were selected to contain gas-rich galaxies with turbulent and marginally stable disks, well representing the star formation conditions at cosmic noon, z ∼ 1 − 3 (Green et al. 2014); their star-forming regions were studied in H down to ∼ 50 pc resolution (Fisher et al. 2017a).In a similar way, the Lyman-Alpha Reference Sample (LARS), at redshift z = 0.03−0.20,contains galaxies selected to be analogs of z ∼ 3 (Lyman-break and Lyman-emitter) galaxies, characterised by elevated UV luminosities and SFRs (Östlin et al. 2014); their stellar clumps were studied in far-UV emission down to 10 pc scales (Messa et al. 2019).Finally, blue compact galaxies (BCGs) are (usually low-mass) star-forming systems with high specific star formation rate (e.g.Östlin et al. 2001); among the best studied ones, SBS 0335-052E, ESO 338-IG04 and Haro 11 are known to host populations of bright young stellar clusters/clumps (e.g.Östlin & Kunth 2001;Östlin et al. 2003;Adamo et al. 2010Adamo et al. , 2011;;Sirressi et al. 2022).Sizes and UV magnitudes of clumps in ESO 338-IG04 were derived in Messa et al. (2019).No size measurements are available in the literature for clumps in SBS 0335-052E and Haro 11; we therefore perform, for the clumps in these galaxies, the same size-luminosity analysis described in Section 3.2 in the F140LP filter, tracing FUV emission.
The sizes, UV magnitudes (and SFRs) and Σ SFR distributions of the clumps in these nearby samples are shown in Fig. 7 (top-right and bottom-left panels).In the case of SBS 0335-052E we plot the .The sample at z = 0 from SINGS (Kennicutt et al. 2003) is shown as black density contours (enclosing 1 and 2 of the distribution).On the y-axis we are showing either H luminosities converted to SFR values (for L12 and L15) or UV magnitudes (for the rest of the samples).Values of typical SFR surface densities for samples at  = 0, 1, 3 and 5, as derived by L15, are shown as dashed lines.Individual redshifts (or ranges) for each of the samples is given in the bottom-right panel.The references for the samples considered are: the Cosmic snake arc (Cava et al. 2018), A521-sys1 (Messa et al. 2022), L12 Livermore et al. (2012a), the Sunburst arc (Vanzella et al. 2022b), SDSSJ1110+6459 (Johnson et al. 2017), Abell2744-arc (Vanzella et al. 2022c), the Sunrise arc (Vanzella et al. 2022a), SMACS0723 (considering only clumps at z > 5, i.e.where JWST-NIRCam traces their rest-frame UV emission, see Claeyssens et al. 2023), L15 (Livermore et al. 2015), V17&V19 Vanzella et al. (2017a,b, 2019), MACSJ0416 (Meštrić et al. 2022).In case of multiply-imaged clumps, we consider only the one with the largest magnification.(Top right:) comparison to nearby clump samples of: blue compact galaxies (ESO338-IG04, Messa et al. 2019;SBS0335-052E and Haro11, this work), LARS (Messa et al. 2019) and DYNAMO (Fisher et al. 2017a).On When considering only the compact (R eff < 100 pc) clumps, in Fig. 7 (bottom-left), their median SFR surface densities are similar to the ones of clumps observed in galaxies at 0 <  < 3.5; the most extreme cases are as dense as the clumps in the z ∼ 5 of the sample studied in this work.These distributions suggest that the environmental conditions setting the starburst nature of these galaxies (interactions, mergers, elevated gas fraction and turbulence) can indeed drive the formation of clumps with elevated SFR densities.

Mass surface densities
We investigate if the SFR redshift evolution discussed in Section 5.2.1 (and Fig. 7) is also reflected in the clump masses.We plot in Fig. 8 the clumps mass surface densities, Σ M ★ , as a function of the clump sizes (left panel) and as violin distributions (right panel); we use the same clump samples as in Fig. 7, when mass values are available.In more detail, masses are unavailable for L12, L15, SDSSJ1110, the nearby DYNAMO sample and the SINGS sample of local galaxies; instead, we plot the LEGUS sample of local clusters (from the analysis of Brown & Gnedin 2021), to help the comparison of compact clumps to single stellar clusters.For comparison we provide also the stellar surface density of a typical 10 5.2 M ⊙ globular cluster with R eff = 3 pc (Brodie & Strader 2006).In this analysis we consider the entire SMACS0723 sample by Claeyssens et al. (2023) and not only the one at z > 5 as done in Fig. 7, because the selection here is based on clump mass and not on FUV luminosity.The references for the samples are the same reported in the caption of Fig. 7, except for masses of clumps in: (i) SBS0335-052E, from Adamo et al. ( 2010); (ii) Haro 11, from Sirressi et al. (2022); (iii) ESO338-IG04 and LARS, from a forthcoming paper (Messa et al., in prep).Mass surface densities span almost 5 orders of magnitudes (Σ M ★ ∼ 10 0 − 10 5 M ⊙ pc −2 ).Contrary to the Σ SFR values (Fig. 7), no clear redshift trend is visible in the left panel of Fig. 8.Some of the clumps in the current study, from MS1358 and MACS0940, stands out, together with other few z > 1.5 sources as very dense (Σ M ★ > 10 4 M ⊙ pc −2 ) large (20 − 100 pc) clumps; these are the densities of the densest gravitationally bound local stellar systems (young and globular clusters), yet on much larger scales.The absence of a clear redshift evolution of Σ M ★ can be deduced also from redshift-binned violin distributions (see Appendix B).
A tentative redshift evolution of Σ M ★ can be seen when considering the clumps at R eff ⩽ 20 pc i.e. at the scales of star clusters, Fig. 8, right panel.An apparent feature of the figure is the increase in the upper-end of the Σ M ★ distribution at z > 5, but this is driven by a single massive (M ∼ 10 7 M ⊙ ) compact (R eff = 1.4 pc) source observed in the Sunrise arc (Vanzella et al. 2022a).The high stellar densities of high-z clumps become evident when local clusters are considered (black contours and black distribution in Fig. 8); there is only little overlap between local clusters and their counterparts at higher redshifts, the latter being on average denser.As already suggested in the previous section, this evolution probably reflects the ambient pressure where clusters form; galaxies at higher redshift are characterised by densest environment, which in turn are able to form denser clusters.This environmental effect was observed locally for stellar clusters (e.g.Johnson et al. 2017;Messa et al. 2018;Adamo et al. 2020b).This trend seems to be confirmed by the clumps in the local starburst BCGs, reaching densities comparable to their z > 1.5 counterparts.We would like to point out that the current sample of high-z clumps with R eff ⩽ 20 is very limited; new insight will come from samples observed with JWST, able to combine extreme spatial resolution with a better age and mass characterisation of the clumps (see e.g.Claeyssens et al. 2022;Vanzella et al. 2022a,c).We also notice that at larger scales (R eff > 20 pc, dashed violin distributions in Fig. 8) only clumps from the z > 1.5 samples are observed to reach the most extreme densities (Σ M ★ ≳ 10 4 M ⊙ pc −2 ), again suggesting a redshift evolution even for large star-forming complexes.
As final remark, we remind that UV luminosities, and consequently SFR UV values, are strongly affected by the age and extinction of the clumps; this could be the cause of the different strengths in the redshift evolution of Σ M ★ and Σ SFR UV .

CONCLUSIONS
We presented the analysis of a small sample of stellar clumps in three galaxies at redshift between 4 and 5, namely the lensed arcs beyond the galaxy clusters MS1358, RCS0224 and MACS0940.These galaxies were chosen to be among the most highly magnified systems hosting multiple stellar clumps in a redshift range currently under-represented in clump studies.Each galaxy is multiply-imaged; many of the images are amplified by factors  > 5 and reaching in some cases  > 20.The clumps were studied using multi-band photometry from the HST, providing a maximum angular resolution, for the clump radii, of ∼ 0.02 ′′ ; combined to the large amplifications of these systems, it allows the study of physical scales down to ∼ 10 pc, i.e. comparable to the sizes of individual stellar clusters.
Clump populations in the three galaxies were extracted from a reference rest-frame UV filter (ACS-WFC-F775W for MS1358, ACS-WFC-F814W for RCS0224 and MACS0940); we further checked that our extraction did not miss clumps of different (redder) colors by testing the clump extraction also on the other available HST filters.Clump colors, in combination with the lens models, were used to recognise (and discard) foreground sources in the field; we find in total 10 unique clumps in MS1358, 3 in RCS0224 and 11 in MACS0940.
Clumps sizes and magnitudes were derived in the reference restframe UV filter for each of the galaxies; the other filters were used to fit a broad-band SEDs and derive the clump masses.Intrinsic sizes, magnitudes and masses were derived using the lens models; due to the large amplifications involved, the uncertainties associated to the amplification factors dominate these intrinsic quantities.The derived effective radii range from ∼ 10 to ∼ 200 pc; the smallest sizes are reached in RCS0224, the galaxy with the largest amplification, where we observe sources down to R eff = 6 pc.UV magnitudes are also converted to SFR UV following the conversion by Kennicutt & Evans (2012); they range from ∼ 10 −2 to the most extreme values of 5 M ⊙ yr −1 in MS1358.The completeness limits of the samples, when converted into SFR UV surface densities, Σ SFR , are typically ∼ 10 M ⊙ yr −1 kpc −2 , and in all cases above 1 M ⊙ yr −1 kpc −2 , i.e. above the typical value for clumps in local main sequence galaxies (Kennicutt et al. 2003;Livermore et al. 2015); we deduce that an eventual population of clumps with low surface brightnesses would be missed in this study.
By focusing individually on the galaxies of this study we find that: • the morphology of MS1358 is dominated by two bright clumps, accounting for 20% and 5% of its rest-frame UV emission; when all the 10 observed clumps are considered together, this fraction raises to ∼ 40%.The UV-derived SFR, SFR UV , of the brightest clumps are lower than the literature values derived from nebular emission (Swinbank et al. 2009), suggesting either the presence of dust extinction or that the current star-formation episode is much younger than the assumption used in the Kennicutt & Evans (2012) FUV-to-SFR conversion.Clump masses range between 10 7 and 10 9  ⊙ ; they are consistent with the entire galaxy mass as derived by Swinbank et al. (2009).We deduce that clumps in MS1358 are major contributors to the recent star formation episode(s) in the galaxy and to the build-up of its mass.The two main clumps of MS1358 are the UV-brightest among the galaxies studied in the current work and among all the compact (R eff < 100 pc) clumps known in literature.Despite most of the clumps in this galaxy have sizes in the range of 50 − 100 pc, their mass surface densities are comparable (and in many cases higher) to the ones of the densest local stellar clusters.
• RCS0224 is characterised by 3 bright and compact clumps, accounting for ∼ 45% of the rest-frame UV emission.The galaxy also shows a larger region of diffuse UV emission, appearing as a very elongated arc devoid of clumps, but actively forming stars, as derived from nebular [OII] emission by Swinbank et al. (2007); deeper observation would be needed to test the presence of lowsurface brightness clumps along the arc.Despite RCS0224 has the shallowest data in the studied sample, its large magnification allows to reach very small intrinsic scales; the clumps sizes and masses, R eff = 6 − 25 pc and M ★ = 10 6 − 10 7.5 M ⊙ respectively, are close to the ones of the most massive stellar clusters in local galaxies.As is the case for clumps MS1358, also in RCS0224 the clump densities reach higher values than what typically observed in local samples.
• Four bright clumps characterise the rest-frame UV morphology of the lensed arc MACS0940, accounting for ∼ 40% of the emission, with other 4 sources contributing to another ∼ 10%.Their derived intrinsic sizes (R eff = 10−100 pc) and masses (M ★ = 10 6 −10 8 M ⊙ ) suggest also for these clumps extreme stellar densities.
Finally, we compare the SFR and stellar mass surface densities of the clumps to the ones of known samples in the literature.We find overall an increase of Σ SFR with redshift, particularly when comparing clumps in local main sequence galaxies to their counterparts at cosmic noon z ∼ 1 − 3.5; the evolution is less prominent at higher redshifts.A weaker evolution is suggested for Σ M ★ (more evident for very compact sources, R eff ⩽ 20 pc).We can interpret the evolution of both clump quantities in the context of the evolution of the properties of their host galaxies; the latter, at higher redshifts, are characterised by increasingly denser environments, which produce denser galactic (and sub-galactic) star-forming regions and, in turn, also denser stellar products.This interpretation is supported by the study of nearby starburst galaxies, whose clump properties resemble high-z samples more than the local average ones.With current data we do not find any discontinuity in the redshift evolution of SFR and M ★ surface densities when considering redshift earlier than cosmic noon (z > 3.5); this may imply that the clump formation conditions at z ∼ 5 are not different compared to later times.This result is in line with recent ALMA findings suggesting the presence of turbulent disk galaxies even at z > 5 (Lelli et al. 2021;Jones et al. 2021;Rizzo et al. 2021;Herrera-Camus et al. 2022;Parlanti et al. 2023).
We remark that these redshift trends and relative interpretations are still tentative, for two main reasons.First, populations of clumps below ∼ 100 pc are still limited, leading to a small sample statistics when a binning in redshifts is considered.Second, the samples considered are in-homogeneous in how they have been extracted and analysed; this is true especially for what concerns SED-derived quan-tities, e.g.clump masses, ages and extinctions.On going observations with the JWST will go in the direction of tackling both problems by observing large galaxy samples and extending the rest-frame optical coverage also to high-redshift sources.

Figure 1 .
Figure 1.Observed (lensed)rest-frame UV magnitudes and angular area covered by spectroscopically confirmed  > 3.5 galaxies gravitationally lensed by galaxy clusters taken from the literature (see text for details).The sizes of the markers is proportional to the number of clumps observed in the galaxy.The final selection of 3 galaxies analysed in the current work are highlighted by the following markers: a circle (MS1358), a diamond (RCS0224) and a square (MACS0940).

Figure 2 .Figure 3 .
Figure 2. RGB composite (r:F160W, g:F110W, b:F775W) of the galaxy cluster field MS 1358+62, containing the  = 4.93 lensed western (W) and eastern (E) arcs and their counter-image (CI).The zoomed-in inset, showing the observations in the reference F775W filter, highlights the positions and names of the extracted stellar clumps in all the different multiple images.A partial northern image (N) of the lensed galaxy contains only one clump and is labeled as 'N_1'.Dashed lines show regions of constant magnifications at  = 5, 10 and 50 (green, red and cyan, respectively), at the redshift of the lensed arcs.Reference angular scale are given; the image is aligned with north-up and east-left.

Figure 4 .
Figure 4. Same as Fig.2, for MACS0940 (r:F160W, g:F125W, b:F814W, zoom-ins:F814W).According to the lens model of the system, the SE and SW arcs (and all the clumps they host) are different lensed images of the main galaxy (labeled as "0" in the N and E fields).The clumps 9, 10 and 11, observed in the N and E fields, do not have counterparts in the SE and SW fields.
The extracted clumps are shown, together with their assigned IDs, in Figs.2,3 and 4. The lens models predictions were used to match the same clump over different images.Ten and 8 clumps are observed in the western (W) and eastern (E) arcs of MS1358, respectively; only 5 of them are clearly visible in the counter-image (CI).The brightest clump (ID:1) is the only one visible in the northern image (N_1).RCS0224 appears in the eastern image (E) as composed of two bright regions, one of which is composed by two sub-clumps.Only one of such regions (containing clumps 2 and 3), is visible in the western image (W); the central image (C) is instead characterised by a single source (either clump 2 or 3), with a very high magnification.The counter-image (CI) of RCS0224 appears as a single source, where clumps are indistinguishable.MACS0940 counts a total of 8 clumps in the eastern arc (SE), the one with the largest magnification

Figure 5 .
Figure 5. (Top panels:) observed sizes and magnitudes (in the reference filter) for all the clumps in this study.Estimates of the observed sizes consider the angular diameter distance of the host galaxy but do not take into account their lensing magnification.The gray shaded areas highlight the size regions below our resolution limits.Black dot markers (connected by dashed lines) indicate the extraction limits for each of the simulated sizes considered (see main text in Section 3.2.1 for details); similarly, markers connected by dotted lines are used to indicate completeness limits.(Bottom panels:) Intrinsic sizes and magnitudes of the clumps.UV magnitudes have been converted into SFR UV values using the Kennicutt & Evans (2012) relation on the right-hand side of each panel.Uncertainty bars combine both photometric and lensing uncertainties.The extraction limits shown in the top panels have been converted here to surface brightness limits (dashed lines, each line corresponds to the limit derived for one simulated size).Solid gray lines are of constant SFR surface density for 1, 10, 10 2 and 10 3 M ⊙ yr −1 kpc −2 .In all panels, open symbols are used for upper limits and same markers are used for the same clumps in different counter-images.

Figure 6 .
Figure 6.Intrinsic sizes and masses of the clumps.Solid gray lines are of constant mass surface density for 10 2 , 10 3 and 10 4 M ⊙ pc −2 .The same color and marker notation as in Fig. 5 is used.

Figure 7 .
Figure7.(Top left:)  sizes and magnitudes of the clumps in the current study compared to clump samples from literature, color-coded by their redshift (blue: z < 1.5, purple: 1.5 ≤ z < 3.5, red: 3.5 ≤ z < 5, orange: z ≥ 5).The sample at z = 0 from SINGS(Kennicutt et al. 2003) is shown as black density contours (enclosing 1 and 2 of the distribution).On the y-axis we are showing either H luminosities converted to SFR values (for L12 and L15) or UV magnitudes (for the rest of the samples).Values of typical SFR surface densities for samples at  = 0, 1, 3 and 5, as derived by L15, are shown as dashed lines.Individual redshifts (or ranges) for each of the samples is given in the bottom-right panel.The references for the samples considered are: the Cosmic snake arc(Cava et al. 2018), A521-sys1(Messa et al. 2022),L12 Livermore et al. (2012a), the Sunburst arc(Vanzella et al. 2022b), SDSSJ1110+6459(Johnson et al. 2017), Abell2744-arc(Vanzella et al. 2022c), the Sunrise arc(Vanzella et al. 2022a), SMACS0723 (considering only clumps at z > 5, i.e.where JWST-NIRCam traces their rest-frame UV emission, seeClaeyssens et al. 2023), L15(Livermore et al. 2015),V17&V19 Vanzella et al. (2017a,b, 2019), MACSJ0416 (Meštrić et al.  2022).In case of multiply-imaged clumps, we consider only the one with the largest magnification.(Top right:) comparison to nearby clump samples of: blue compact galaxies (ESO338-IG04,Messa et al. 2019; SBS0335-052E and Haro11, this work), LARS(Messa et al. 2019) and DYNAMO(Fisher et al. 2017a).On the y-axis we are showing H luminosities converted to SFR values for DYNAMO, and UV magnitudes for the rest of the samples.(Bottom left:) distribution of SFR densities (shown as violin plots, with median and extreme values of the distribution marked) of the clump samples, binned by redshift; only clumps with R eff < 100 pc are considered.The clumps from the current work are shown also as a separate sample.
Figure7.(Top left:)  sizes and magnitudes of the clumps in the current study compared to clump samples from literature, color-coded by their redshift (blue: z < 1.5, purple: 1.5 ≤ z < 3.5, red: 3.5 ≤ z < 5, orange: z ≥ 5).The sample at z = 0 from SINGS(Kennicutt et al. 2003) is shown as black density contours (enclosing 1 and 2 of the distribution).On the y-axis we are showing either H luminosities converted to SFR values (for L12 and L15) or UV magnitudes (for the rest of the samples).Values of typical SFR surface densities for samples at  = 0, 1, 3 and 5, as derived by L15, are shown as dashed lines.Individual redshifts (or ranges) for each of the samples is given in the bottom-right panel.The references for the samples considered are: the Cosmic snake arc(Cava et al. 2018), A521-sys1(Messa et al. 2022),L12 Livermore et al. (2012a), the Sunburst arc(Vanzella et al. 2022b), SDSSJ1110+6459(Johnson et al. 2017), Abell2744-arc(Vanzella et al. 2022c), the Sunrise arc(Vanzella et al. 2022a), SMACS0723 (considering only clumps at z > 5, i.e.where JWST-NIRCam traces their rest-frame UV emission, seeClaeyssens et al. 2023), L15(Livermore et al. 2015),V17&V19 Vanzella et al. (2017a,b, 2019), MACSJ0416 (Meštrić et al.  2022).In case of multiply-imaged clumps, we consider only the one with the largest magnification.(Top right:) comparison to nearby clump samples of: blue compact galaxies (ESO338-IG04,Messa et al. 2019; SBS0335-052E and Haro11, this work), LARS(Messa et al. 2019) and DYNAMO(Fisher et al. 2017a).On the y-axis we are showing H luminosities converted to SFR values for DYNAMO, and UV magnitudes for the rest of the samples.(Bottom left:) distribution of SFR densities (shown as violin plots, with median and extreme values of the distribution marked) of the clump samples, binned by redshift; only clumps with R eff < 100 pc are considered.The clumps from the current work are shown also as a separate sample.

Figure 8 .
Figure 8. (left panel): mass surface densities of clumps, in function of their intrinsic sizes; the color and marker coding of the plotted samples are the same presented in Fig. 7; the only differences are the LARS sample, plotted as yellow contours instead of square markers, and the inclusion of the LEGUS sample of local stellar clusters (from the analysis of Brown & Gnedin 2021), instead of SINGS, as black contours, due to the lack of available masses for the latter.Both yellow and black contours enclose 1 and 2 of the relative distribution.The location of a typical globular cluster with M = 10 5 M ⊙ and R eff = 3pc (Brodie & Strader 2006) is shown with a black star marker.(right panel): Σ M★ distributions, shown as violin plots; the distributions of clumps with R eff ⩽ 20 pc are shown as filled violins, while clumps in range 20 < R eff ⩽ 100 pc are shown with empty violins with dashed edges (a dashed line marks also the median value of each distribution).

Figure A1 .
Figure A1.Data (left column), best-fit models (central column) and residuals (right column) for the main lensed arcs, in the reference filters of each galaxy.Clumps IDs are reported in the central panels.Foreground sources are marked with red arrows in the right-hand panels.Grids of 1 × 1 arcsec are plotted to facilitate the comparison between panels.For the position of the galaxy images within the cluster fields, see Fig. 2, 3 and 4.