© The Authors 2025

1 Introduction

The clumpy structure, typical of distant galaxies, was first observed in Hubble Space Telescope (HST) images of unlensed sources in deep fields and was attributed to giant star-forming regions with sizes of ~1 kpc (e.g., Elmegreen & Elmegreen 2006; Wisnioski et al. 2012). The increase in intrinsic spatial resolution provided by gravitationally lensed fields allowed these large regions to be broken into smaller sub-components, namely clumps with sizes spanning the range of 10–500 pc (e.g., Johnson et al. 2017; Vanzella et al. 2017; Cava et al. 2018; Vanzella et al. 2019; Meštrić et al. 2022; Messa et al. 2022, 2024). Recently, thanks to the James Webb Space Telescope’s (JWST) resolution and sensitivity, high-z stellar clumps have become observable at ~100 pc scales even in unlensed fields (e.g., Ji et al. 2024; Kalita et al. 2024, 2025a,b); when considering the aid of gravitational lensing, JWST can probe clumps down to <10 pc in size and 10⁵ M_⊙ in mass, reaching the scales of individual star clusters. For instance, “Cosmic Gems” is a gravitational arc at z ~ 10.2 (Adamo et al. 2024), lensed by the galaxy cluster SPTCL J0615-5746 and initially discovered by HST in the RELICS survey (Salmon et al. 2018). Using JWST observations, Adamo et al. (2024) detected five different star clusters located in a region smaller than 70 pc and with intrinsic sizes ~1 pc, with ages younger than 50 Myr. Other notable examples are, for example, the “Sunrise Arc” at z ~ 6 (Vanzella et al. 2023), where six stellar clusters have been detected (characterized by sizes ~1–20 pc and ages in the range 1–30 Myr) and the “Firefly Sparkle” galaxy at z ~ 8.3 (Mowla et al. 2024), containing ten star clusters (with sizes < 10 pc and masses ~ 10⁵ M_⊙) whose star formation histories (SFHs) show a recent starburst, probably triggered by the interaction with a nearby galaxy. The “redder” wavelength range probed by JWST, compared to HST, may also unclose the study of old clumpy systems. Several stellar clumps have been identified in the “Sparkler” galaxy at z ~ 1.37, lensed by the galaxy cluster SMACS 0723; these systems are consistent with being red, old, evolved globular clusters with formation redshift z_form ~ 7–11 and ages between 1 and 4 Gyr (Mowla et al. 2022; Claeyssens et al. 2023; Adamo et al. 2023).

Hence, when the angular resolution increases thanks to more advanced instruments or the magnification granted by gravitational lensing, single stellar clusters are also detectable in the high-redshift Universe. This is a crucial achievement since it allows us to map the star cluster formation in early galaxies. Young star clusters (YSCs) are the natural site of formation for massive stars (e.g., Zinnecker & Yorke 2007; Portegies Zwart et al. 2010); their detection at high-z can unveil their role in the reionization process, which took place in the redshift interval between z = 5 and z = 15 when the first sources started to form and their photons led to the ionization of the neutral hydrogen permeating the Universe. Observing compact star-forming regions is also fundamental for our understanding of galaxy growth and evolution since they alter the properties of the interstellar medium of the hosting galaxy thanks to their radiative and mechanical feedback. This interconnection with the host is closely related to the clump formation scenario, which is still debated and may evolve with the host’s redshift. The first suggested way is a formation in situ, due to the fragmentation of the disk of the host galaxy triggered by instabilities (Agertz et al. 2009; Dekel et al. 2009; Bournaud et al. 2010). The second scenario sees the clumps originating from minor merger events (Mandelker et al. 2017; Zanella et al. 2019).

Large statistics are needed to derive physical properties of the stellar clumps and clusters robustly. In one of the largest attempts, Claeyssens et al. (2025) studied ~2000 individual clumps in 476 galaxies at redshift z ~ 0.7–10 lensed by the galaxy cluster A2744, using JWST NIRCam observations, and recovering the main clump physical properties, such as their sizes, ages, and masses. They show an evolution of the clumps’ properties with redshift (clumps are, on average, smaller and younger at higher redshifts) and correlations between the properties of the clumps and those of their host galaxy (massive galaxies host more massive and older clumps). This sample of high-z clumps is a factor ~10 larger than the one studied in SMACS0723 by Claeyssens et al. (2023), which, in turn, doubled the number of compact (<100 pc) clumps at z > 1 known before JWST (e.g., Meštrić et al. 2022; Messa et al. 2024). The current HST and JWST samples of high-z lensed clusters are largely inhomogeneous in terms of sizes and luminosity probed; the large difference in magnification, from one source to the other within a field, makes the estimate of completeness a very complex problem, even when focusing on subsamples within a small redshift range. This complexity limits the possibility of studying trends among systems and across redshifts.

For this work, we aimed to implement a pipeline to simulate lensing effects by intervening mass distributions such as galaxy clusters on high-redshift star-forming galaxies hosting multiple star clusters. The pipeline also allows observations to be simulated with various instruments, such as HST and JWST. The pipeline uses functions and methods from the multipurpose lensing library PYLENSLIB (Meneghetti 2021), written in Python. In this paper, we illustrate the pipeline functionalities (Sects. 2–5) and discuss some potential usages (Sects. 6–8). In particular, this software can be particularly useful for studying the detectability of lensed stellar clusters at any redshift, especially those at the low-mass end of the star cluster mass function (e.g., Adamo et al. 2020). In this study, we consider a stellar clump to be “detected” only if it is visually identifiable in the images. This choice is motivated by the relatively small sample analyzed here, which allows for careful visual inspection. In future studies involving larger datasets or simulations, we plan to complement visual inspection with automatic detection software, such as SExtractor (Bertin & Arnouts 1996), following the method described in Claeyssens et al. (2025), to efficiently and systematically identify clumps.

In this paper, we use the expression star formation clump when referring to a sub-component of a galaxy (typically on physical scales of 10–1000 pc). In contrast, the term stellar cluster represents the usual gravitationally bound YSC (with typical sizes <10 pc, e.g., Ryon et al. 2015, 2017; Brown & Gnedin 2021). Hence, depending on the angular resolution, a clump can be a complex or a single stellar cluster. We assume a flat cosmology with Ω_M = 0.3, Ω_Λ = 0.7, and H₀= 70 km s⁻¹ Mpc⁻¹.

2 Source and lens planes

When simulating lensing effects by a given gravitational lens, we use the thin screen approximation, that is, we assume that the lens mass distribution can be projected on a plane at redshift z_l, perpendicular to an arbitrary optical axis passing through the observer position. Similarly, the source is assumed to lay on a plane perpendicular to the same axis at redshift z_s. At the moment, we do not consider the combined effects of multiple mass distributions along the line of sight. This scenario can be easily implemented by means of a multiplane ray-tracing algorithm, which is already implemented in PYLENSLIB.

We define the angular coordinates θ = (θ₁, θ₂) and β = (β₁, β₂) on the lens and source planes, respectively. Each point on the lens plane can be mapped onto the source plane using the lens equation,

$$\mathit{\beta }\mathbf{=}\mathit{\theta }\mathbf{-}\mathit{\alpha }\mathbf{(}\mathit{\theta }\mathbf{)},$$(1)

where α(θ) is the lens deflection angle at position θ.

The lens equation also allows us to relate the surface brightness of the lensed images, I(θ), to the intrinsic surface brightness of the source, I_s(β), using the ray-tracing method. This technique consists of tracing bundles of light rays through the lens plane at positions (θ₁, θ₂), where the deflection angles can be read off the deflection angle maps, α₁(θ₁, θ₂) and α₂(θ₁, θ₂). Then, the arrival position of each ray on the source plane can be derived from the lens equation, and the image surface brightness is finally computed as

$$I(\mathit{\theta })=I_s[\mathit{\theta }-\mathit{\alpha }(\mathit{\theta })].$$(2)

3 Modeling the sources

We aim to generate clumpy galaxies for image simulations, including lensing. We model these sources using multiple luminous components, describing the hosts and the stellar clumps. Analytical functions describe the brightness profiles of each component. Specifically, we use the Sérsic profile (Sérsic 1963), given by

$$I_s(\beta )=I_e\exp \{-b(n)[{\left(\frac{\beta }{R_e}\right)}^{1/n}-1]\},$$(3)

where β is the angular distance from the source center (β_s,1, β_s,2),

$$\beta =\sqrt{{({\beta }_1-{\beta }_{s,1})}^2+{({\beta }_2-{\beta }_{s,2})}^2},$$(4)

I_e is the surface brightness at the effective radius R_e, n is the Sérsic index and b is a function of the Sérsic index well approximated by $b(n)=2n-\frac{1}{3}+\frac{4}{405n}$ (Sérsic 1963).

On average, bulges and elliptical galaxies have the steepest central profiles, with 2 < n < 10. Galactic disks have exponential profiles with n ~ 1. Bars have flatter central profiles with n ≤ 0.5, while the peculiar case with n = 0.5 corresponds to a Gaussian brightness profile, often used to describe stellar clumps (e.g., Cava et al. 2018).

3.1 Host galaxy

Our simplest description of galaxies hosting stellar clumps assumes that they can be well approximated by dual brightness distributions consisting of a spherical bulge and an elliptical, exponential disk. Each of these components has its own effective surface brightness and radius, I_e,bulge, R_e,bulge and I_e,disk, R_e,disk. For the examples in this work, we use Sérsic profiles with n = 1 for the disk and n = 4 for the central bulge. These values can be adapted to the sources one wishes to simulate.

Thus, for such a galaxy, the total surface brightness distribution on the source plane is

$$\begin{array}{rl}{\mathit{I}}_s(\mathit{\beta })=& I_{e,\text{bulge}}\exp \{-b(4)[{\left(\frac{\beta }{R_{e,\text{bulge}}}\right)}^{1/4}-1]\}\\ & +I_0\exp (-\frac{\tilde{\beta }}{R_h}),\end{array}$$(5)

where I₀ = I_e,disk exp (1.678), and R_h = R_e,disk/1.678. To include ellipticity in the disk surface brightness, we define

$$\tilde{\beta }=\sqrt{{({\beta }_1-{\beta }_{s,1})}^2/q^2+{({\beta }_2-{\beta }_{s,2})}^2},$$(6)

where q is the disk axis ratio.

We can simulate arbitrary position angles φ for the source galaxy by applying a rotation of the reference frame:

$$\begin{array}{c}{\beta }_1^{\mathrm{\prime }}=({\beta }_1-{\beta }_{s,1})\cos \phi +({\beta }_2-{\beta }_{s,2})\sin \phi ,\\ {\beta }_2^{\mathrm{\prime }}=-({\beta }_1-{\beta }_{s,1})\sin \phi +({\beta }_2-{\beta }_{s,2})\cos \phi ,\end{array}$$(7)

where (β₁, β₂) and (${\beta }_1^{\mathrm{\prime }},{\beta }_2^{\mathrm{\prime }}$) are the original and new coordinates, respectively.

Fig. 1 (a) Example of a host galaxy consisting of an exponential disk with two spiral arms and a spherical bulge. Details on the adopted parameters are reported in the text. (b) Model of a host galaxy, made of an exponential disk, a spherical bulge, and two spiral arms, populated with Sérsic stellar clumps. Details on the adopted parameters are reported in the text. The galaxy center coincides with the center of the field of view. (c) RGB color image of the same source. Here, we assume the bulge is 1 Gyr old, while the disk is 500 Myr old. For the stellar clumps, their ages span in the range 1–500 Myr. In this case, the red, green, and blue channels correspond to effective wavelengths λe = 5, 3, and 2 μm, respectively (0.70, 0.42 and 0.28 μm, rest-frame, taking into account an assumed redshift z = 6.145). These images have an angular resolution of 0.25 mas/pix, which at a z ~ 6 corresponds to a physical scale of approximately 1 pc per pixel. As a result, our clumps, all of which have sizes between 1 and 10 pc, are resolved in the high-resolution images.

3.2 Spiral arms

The disk in the source model outlined in Sect. 3.1 has a smooth brightness distribution. Our model allows the addition of a spiral structure to the smooth disk.

To achieve this goal, we follow the method of Metcalf et al. (2019), where the surface brightness of a circular disk including N_a spiral arms is described by

$$\begin{array}{c}I(\beta ,\vartheta {)}_{disk}=I_0\exp (-\frac{\beta }{R_h})[1+A\cos (N_a\vartheta +{\varphi }_r)],\\ {\varphi }_r=\alpha \log \left(\frac{2\beta }{R_h}\right)+{\varphi }_d,\end{array}$$(8)

where ϑ is a polar angle, ϕ_d is the phase angle of the arms, and α describes how the arms are wrapped around the bulge.

The equations above correspond to a circular disk observed face-on. Therefore, in this specific case Eq. (6) simplifies to:

$$\beta =\sqrt{{({\beta }_1-{\beta }_{s,1})}^2+{({\beta }_2-{\beta }_{s,2})}^2}.$$(9)

To simulate different viewing angles of the same galaxy, we introduce a new parameter, $\tilde{\phi }$, which accounts for the inclination of the disk with respect to the optical axis. By incorporating this parameter, Eq. (7) becomes:

$$\begin{array}{c}{\beta }_1^{\mathrm{\prime }}=({\beta }_1-{\beta }_{s,1})\cos \phi +({\beta }_2-{\beta }_{s,2})\sin \phi ,\\ {\beta }_2^{\mathrm{\prime }}=-({\beta }_1-{\beta }_{s,1})\sin \phi +({\beta }_2-{\beta }_{s,2})\cos \phi /\cos \tilde{\phi }.\end{array}$$(10)

Finally, to make the spiral arms more irregular, we perturb the disk simulating a Perlin noise¹.

In the first panel of Fig. 1, we show an example of a face-on spiral galaxy $(\tilde{\phi }=0)$ on the source plane. We adopt an axis ratio q = 1 for the disk, while the effective radius is R_e,disk = 0.18″. For the arms we use A = 1.0, N_a = 2, ϕ_d = 1, and α = 5.0. The bulge effective radius is R_e,bulge = 0.1″. We use Eqs. (5) and (8) to compute the source brightness distribution on a regular grid with 2000 × 2000 pixels, with an angular resolution of 0.25 mas/pix. The resolution of this grid can be properly chosen according to the purpose of the simulation, as it will be discussed below (see Sect. 5).

In the following Sections, we will use the source galaxy in the left panel of Fig. 1 as the host for the stellar clumps. In addition, we will calibrate the source flux to mimic the observations of the Cosmic Archipelago arc in galaxy cluster MACS J0416.1-2403; this is a highly magnified z=6 Lyman-α emitter (Vanzella et al. 2017; Vanzella et al. 2019), hosting massive star clusters, as revealed by recent JWST observations (Messa et al. 2025). Since we do not have a proper estimate of the host galaxy stellar mass from SED fitting, we assume a value M_*,tot = 10⁸ M_⊙. This value is consistent with stellar masses measured for some high-z galaxies (e.g., Vanzella et al. 2022), and is therefore a plausible choice for our simulated source.

3.3 The stellar clumps

This Section describes how we populate the host galaxy with stellar clumps. We make the following assumptions:

1. the clump spatial distribution follows the underlying brightness distribution of the host;

2. the total mass in clumps is a predetermined fraction of the total mass of the host;

3. the distribution of clump masses obeys a Schechter-like function (e.g., Schechter 1976);

4. the clump size is a growing function of the clump mass.

While our modelization is extremely flexible, the model parameters discussed below should be set to match observations. Unfortunately, the statistical properties of stellar clumps in distant, high-z galaxies are not well constrained (see, e.g., Meštrić et al. 2022; Claeyssens et al. 2025). For this reason, we use observations of local galaxies as mild references. The similarity of these sources to their high-z equivalents is not guaranteed. In the future, we plan to investigate the statistical properties of stellar clumps in high-z sources by comparing simulations realized with our tool to real observations. Such an application can help us understand the level of similarity between clumpy sources in the high- and low-z Universe.

Firstly, we constrain the clump’s overall number, considering that only a fraction Q of the host galaxy’s total stellar mass goes into star clumps. This fraction varies with the morphological classification of a galaxy: it is reasonable to assume Q = [0–0.1] for early-type, lenticular galaxies and Q = [0.2–0.4] for late-type galaxies (Plazas et al. 2019), while even higher fractions can be reached at high-z, where extreme gas physical conditions are not rare in star-forming galaxies (Adamo et al. 2020) which can be characterized by high stellar and gas densities (Walker et al. 2015). We expect late-type galaxies to be characterized by stronger star formation and host more star-forming regions than early-type galaxies. For the simulation of the source at z = 6.145, we can assume extreme star formation conditions by using Q = 0.8.

We assume that the clump stellar masses follow a distribution which is a generalization of a Schechter mass function,

$$f(x)=x^{\beta }\exp [-\delta {\left(\frac{x}{x_{cut}}\right)}^{\gamma }],$$(11)

which combines a power-law behavior with slope β at low masses with an exponential cut-off above a characteristic scale x_cut. In the previous equation, x is the clump mass in solar masses and δ and γ are parameters defining the shape of the mass function in the high-mass tail. This parameterization of the clump mass function is borrowed from cosmological hydrodynamical simulations of galaxy formation and evolution (e.g., Larsen 2010). It can be adapted to observations of local galaxies, which are consistent with setting β = −2.0, δ = 1.0, and γ = 1.0 (e.g., Krumholz et al. 2019; Adamo et al. 2020).

We sample the mass function in the range between x_min = 10⁴ M_⊙ and x_max = 10⁷ M_⊙. Furthermore, we define x_cut = 5 × 10⁶ M_⊙, to limit the number of very massive stellar clumps, considering that the total galaxy mass is assumed to be 10⁸ M_⊙. The sampling continues until we saturate the mass budget in stellar clumps, as defined by the parameter Q. To conserve the total galaxy mass, we subtract the sum of the clump masses from the host mass. The upper panel of Fig. 2 shows the stellar mass function of the star clumps resulting from this sampling.

As anticipated, we assume that the clump spatial distribution follows the surface brightness distribution of the host galaxy. This assumption is reasonably well supported by observations of local galaxies (see e.g., Adamo et al. 2020). Thus, we draw the clump positions in the host by sampling the source pixel surface brightness distribution.

We model the stellar clumps using Sérsic profiles with n = 0.5, assuming they have elliptical shapes with axis ratios uniformly distributed in the range [0.3,1). Their position angles are also randomly chosen in the range [0, π). To determine their sizes, we use the relation between effective radius and clump mass measured from YSCs in the Local Universe by Brown & Gnedin (2021):

$$R_e=2.55{\left(\frac{M}{{10}^4{M}_{\odot }}\right)}^{0.24}.$$(12)

Claeyssens et al. (2023) proposed a different mass-size relation for clumps with sizes ≥10 pc in high-z (z > 1) galaxies². Since we want to simulate objects with sizes down to ~ pc scales, we prefer to rely on the Brown & Gnedin (2021) relation. Hence, as shown in the lower panel of Fig. 2, our clumps can span a range of effective radii from ~10 pc down to ~1 pc. Moreover, at a given stellar mass, we assume that the size distribution follows a lognormal distribution with a standard deviation of σ = 0.3.

The second panel of Fig. 1 shows the same galaxy in the first panel but after including the stellar clumps as discussed in this Section. Further examples of possible morphological variations of the source galaxy, as well as variation in the cluster mass fraction, are presented in Appendix A, highlighting the flexibility of our pipeline.

Fig. 2 Upper panel: magenta histogram showing the stellar mass function of the star clumps populating our simulated galaxy. Lower panel: clump sizes (in units of parsec) as a function of the clump stellar mass. Clumps are represented by magenta points. The underlying average size-mass relation given by Eq. (12) (Brown & Gnedin 2021) is shown in a blue solid line. At a given stellar mass, we assume the size distribution follows a lognormal distribution with a standard deviation of σ = 0.3.

Fig. 3 Upper panel: histogram of the clump’s ages. Lower panel: clump stellar masses as a function of their ages in years.

3.4 Simulating colors

Since we want to produce multiband observations of the clumpy galaxy, we need to assign Spectral Energy Distributions (SEDs) to the host and clumps. To derive the SEDs of each component and, therefore, to produce the expected magnitudes for different instruments and filters, we use stellar population synthesis models. In particular, we use the Yggdrasil population synthesis code by Zackrisson et al. (2011), designed for modeling the SEDs of high-z galaxies³. Yggdrasil predicts the SEDs of integrated stellar populations across a wide range of metallicities, from zero to supersolar, accounting also for nebular emission and dust extinction. The software allows the SEDs to be computed directly or the magnitudes to be computed in any photometric band for a chosen instrument. Among the input parameters required by Yggdrasil to compute the SEDs are the metallicity (Z) and the Initial Mass Function (IMF) for the stellar population; for the examples of the current work, we always assume Z = 0.004, with a Kroupa IMF (e.g., Kroupa 2001). One further assumption is the Star Formation History (SFH) of the system. In particular, the Star Formation Rate (SFR) can be instantaneous, which produces a single-age stellar population, or constant over a specific period of time. For our simulation, we assume for both components of the host galaxy a constant SFR for 100 Myr, while, for the stellar clumps, we assume a constant SFR for 10 Myr. We assume the bulge to be 1 Gyr old, while the spiral arms are assumed to be younger, with an age of 500 Myr. Regarding the stellar clumps, we define their ages by sampling a power-law distribution with slope α = −0.5 and spanning between 1 Myr and 500 Myr. This is in line with cluster studies in nearby galaxies (see e.g., Krumholz et al. 2019). The upper panel of Fig. 3 shows the distribution of the clump’s ages, while the lower panel shows the clump masses as a function of their ages. We point out that different assumptions can be made on the age and SFHs of the clusters and their host; also the use of other stellar population synthesis models can be easily implemented within the pipeline.

The SEDs obtained from all these ingredients must be red-shifted according to the source redshift z_s. Then, they need to be normalized to the correct flux. Along with the wavelength and flux F_Ygg data points, the software provides also the stellar mass (in solar masses) M_Ygg of the stellar population associated with that SED. The normalization is thus achieved by assuming a constant mass-to-light ratio:

$$F_{norm}=F_{Ygg}\frac{M}{M_{Ygg}},$$(13)

where M is the mass of our object.

In the right panel of Fig. 1, we show the RGB colored image of our source galaxy, built using the monochromatic flux from simulated images at the three observed wavelengths: 5 μm (red channel), 3 μm (green channel) and 2 μm (blue channel). The intrinsic, rest-frame wavelengths are approximately 0.70 μm, 0.42 μm, and 0.28 μm, respectively, taking into account the source redshift of z = 6.145.

Fig. 4 Mass map of galaxy cluster MACS J0416.1-2403 (Bergamini et al. 2023) at z = 0.396 (RA: 64.038142 deg, Dec: −24.067472 deg). This model is generated with the software Lenstool. The yellow lines represent the caustics of a source at z = 6.145, while the white lines are the corresponding critical lines on the lens plane. The magenta marker indicates the position of our source (β1 = −13.72 arcsec and β2 = 3.62 arcsec), and the three cyan markers show the corresponding multiple images on the lens plane: a stretched arc at the position θ1 = −35.57 and θ2 = 18.85, and two counter-images in θ1 = −42.66 and θ2 = 8.46 and θ1 = −20.16 and θ2 = 31.53, respectively. In this work, we focus on the gravitational arc.

4 Simulating lensing effects

Our pipeline allows us to include lensing effects by providing deflection angle maps. These maps can be derived differently (e.g., using analytical lens models or ray tracing through numerical simulations). In this work, we use maps obtained by modeling existing galaxy clusters with the public software Lenstool (see, e.g., Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009). More precisely, the examples in this work use the reconstruction of galaxy cluster MACS J0416.1-2403 at redshift z = 0.396. The details of the mass modeling can be found in several papers (Bergamini et al. 2019, 2021, 2023). The mass map for this cluster is shown in Fig. 4, where the white lines represent the cluster critical lines, that is, the contours corresponding to an infinite magnification on the lens plane, for a source at redshift z_s = 6.145. By mapping the critical lines onto the source plane, we compute the lens caustics, represented in yellow in the same Figure. From the position of the caustics, we can identify the points on the source plane where sources can be placed to maximize their magnification and distortion. We select the position on the source plane at z_s = 6.145 for our simulated galaxy that reproduces the same configuration and distortion of the Cosmic Archipelago in MACS J0416.1-2403 (Messa et al. 2025). This configuration features an elongated arc and two counter-images, which, however, are not later included in our simulations. The chosen source position lies inside the cusp of the cluster caustic, ensuring that the lensing effect produces three multiple images with high tangential magnification, as they are located near the tangential critical lines. In Fig. 4, the magenta point marks the selected source position, while the cyan crosses indicate the corresponding multiple images positions. Note that 1) the lensed images are stretched and have a spatial extension much larger than the intrinsic source size due to lens magnification; and 2) the stellar clumps in the image are also lensed and multiply imaged inside the arcs, as they follow the same distortion pattern as the overall galaxy.

5 Inclusion of observational noises

This Section illustrates the main steps to mimic observations with specific instruments and add the appropriate observational noise to the previously generated noise-free images.

For example, we simulate observations of the clumpy galaxies with the Hubble Space Telescope (HST) and James Webb Space Telescope (JWST). The procedure we follow is the same for both of them and consists of the following steps:

1. we draw the proper surface brightness distribution on a pixel grid resembling the resolution of the detector;

2. we convert the source surface brightness into units of counts s⁻¹ pixel⁻¹;

3. we convolve the resulting image with the instrument Point-Spread-Function (PSF);

4. we add the sky background;

5. we add photon noise.

We do not consider other instrumental effects such as Charge Transfer Efficiency or Brighter Fatter effects, nor do we simulate cosmic rays, radiation damage, hot pixels, etc. We assume that these effects have already been corrected in the images we produce.

We begin by assigning a stellar mass to the source we want to simulate, based on observations of the Cosmic Archipelago, in order to ensure that the simulated image is comparable to the observed one. As mentioned in Sect. 3.2, we adopt a total stellar mass of M_*,tot = 10⁸ M_⊙. While we currently lack a robust estimate of the Cosmic Archipelago’s stellar mass from SED fitting, this value is consistent with the stellar masses measured for some high-z galaxies (e.g., Vanzella et al. 2022), and is therefore a plausible assumption for our simulation. In this work, we split the total stellar mass into the disk and the bulge. Moreover, the star clump mass is a fraction of the disk’s. In particular, we refer to a disk of M_*,disk = 9.0 × 10⁷ M_⊙, a bulge of M_*,bul = 1.0 × 10⁷ M_⊙, and 329 stellar clumps with a total mass of M_*,cl = 7.2 × 10⁷ M_⊙, corresponding to a fraction Q = 80% of M_*,disk. We point out that the disk mass includes the mass in clumps.

At first, our tool constructs the source with an angular resolution significantly higher than the telescope detector’s chosen for the simulation. In this example, we adopt 0.25 mas/pix instead of the resolution of the Wide Field Camera 3 (WFC3) onboard HST (0.13 arcsec/pix) or the Near Infrared Camera (NIRCam) onboard JWST (short wavelength channel: 0.031 arcsec/pix, long wavelength channel: 0.063 arcsec/pix). This choice is driven by the need to ensure that the clumps within the source are well sampled before simulating lensing effects and the observation with a specific instrument. If the source were not constructed at such a high resolution, the lensing effect would be applied to intrinsically unresolved clumps, causing them to appear as point-like sources rather than revealing their actual extended structure. Consequently, unresolved clumps would not exhibit the stretched features characteristic of gravitational arcs, but would instead appear as compact, magnified points in the lens plane. While a fine resolution of 0.25 mas/pix is necessary on the source plane to resolve such high-z clumps, we use a coarser resolution of 2.5 mas/pix on the lens plane, which still preserves the morphological details of lensed features in the final simulated images.

Using the instrument Zero-Point (ZP), we convert the magnitudes obtained in Sect. 3.4 into a total number of counts (or data numbers, DN) s⁻¹ registered by the detector. The ZP of an instrument, by definition, is the magnitude of an object that produces one DN per second. The magnitude of an object producing DN counts in an observation of length t_exp is, therefore:

$$m_{\mathrm{A}\mathrm{B}}=-2.5{\log }_{10}\left(\frac{\mathrm{D}\mathrm{N}}{t_{exp}}\right)+{\mathrm{Z}\mathrm{P}}_{\mathrm{A}\mathrm{B}}.$$(14)

To obtain DN s⁻¹, we invert the equation above, using the exposure time t_exp from the real observation that serves as a reference for our simulations. We then scale the brightness distribution of the source so that the total number of DN s⁻¹ matches the value corresponding to the desired input magnitude.

In the next step, we consider the PSF of the chosen instrument. We apply it to the image of the clumpy galaxy in units of DN s⁻¹ pixel⁻¹ via a two-dimensional convolution:

$$\tilde{I}({\theta }_1,{\theta }_2)=\int \int d{\theta }_1^{\mathrm{\prime }}d{\theta }_2^{\mathrm{\prime }}I({\theta }_1^{\mathrm{\prime }},{\theta }_2^{\mathrm{\prime }})\mathrm{PSF}({\theta }_1-{\theta }_1^{\mathrm{\prime }},{\theta }_2-{\theta }_2^{\mathrm{\prime }}).$$(15)

The sky background (bkg hereafter) in an astronomical observation depends on several factors, including the direction in which the observation is carried out. In the examples shown in this paper, we set the background level according to real observations of the sources we wish to simulate. More precisely, we place 50 circular regions over the darkest areas of real images, particularly avoiding every source. The mean bkg flux associated with each pixel can be obtained from these regions.

Finally, the photon noise is computed on the pixel grid, assuming it follows a Poisson distribution with variance given by the sum of the source and background counts. The simulated raw observations are obtained by summing the noise maps and the PSF-convolved source images. More details about this procedure can be found in Meneghetti et al. (2008); Plazas et al. (2019), for example.

6 Preparing the virtual observations

As anticipated earlier, to illustrate how our simulation software works, we chose to focus on a source inspired by the gravitationally lensed source Cosmic Archipelago. This source is particularly interesting for the stellar clumps it contains. We intend not to reproduce the Cosmic Archipelago exactly, but to simulate a similar source size and lensing configuration. Hereafter, we detail the parameters required to model the source as a Sérsic galaxy, specifying those related to the bulge and the disk components. The full list of keywords defined in our pipeline to model the source and to set up the virtual observations are reported in Appendix B.

• Bulge: Sérsic index n = 4, effective radius R_e = 0.1″, axis ratio q = 1, position angle ϕ = π/8;

• Disk: Sérsic index n = 1, effective radius R_e = 0.18″, axis ratio q = 1.0, position angle ϕ = π/8.

Both components are centered at source positions β_s,1 = −13.716581″ and β_s,2 = 3.6176634″, and set at redshift z_s = 6.145. Regarding the spiral arms, we use these parameters: scale height R_h = 0.1″, normalization A = 1.0, number of arms N_a = 4, linking angle ϕ_d = 1.0, wrapping factor α = 10.0, tilting angle $\varphi =\frac{3}{4}\pi$. At this stage, the simulation produces an ideal high-resolution image of the source galaxy, without any instrumental effects such as PSF convolution, exposure time, or noise. This image is shown in the central panel of Fig. 1.

In this example, we simulate observations with two cameras onboard the HST and JWST, namely WFC3 in the IR channel and NIRCam. Different pixel scales and PSFs characterize these instruments. Some details are given below.

HST WFC3/IR. We simulate observations through the IR Channel of the imager WFC3, onboard HST. The native pixel size of this detector is 0.13″/pix. We consider the bands F105W, F125W, and F160W. The PSF models used in our simulations are obtained with the public PSF modeling tool Tiny Tim (Krist et al. 2011).

To define the sky background and exposure time, we use the observations of MACS J0416.1-2403 from the CLASH program. While t_exp is explicitly stated in the data header, the bkg level is estimated following the same approach described in Sect. 5. Specifically, we analyze the CLASH files by placing 50 circular regions over the darkest areas, carefully avoiding all luminous sources. The mean bkg flux associated with each pixel can be obtained from these regions. The t_exp values are listed in the first row of Table 1, while the corresponding bkg are reported in the second row of the same Table.

JWST NIRCam. We simulate observations in the imaging mode of the JWST NIRCam through the short-wavelength channel (0.6–2.3 μm) in the bands F115W, F150W, and F200W and through the long-wavelength channel (2.4–5.0 μm) in the bands F277W, F356W and F444W. In this case, the pixel scale is 0.031″/pix for the short wavelength channel and 0.063″/pix for the long wavelength one. The PSF models for the short-wavelength channel are taken from Plazas et al. (2019)⁴, who generated them using the tool WebbPSF, while those for the long-wavelength channel are directly generated from the same tool. We analyze the data of MACS J0416.1-2403 from the Cycle 1 GTO program (PI Windhorst, Windhorst et al. (2022)) and define t_exp and bkg following the same method used for HST WFC3/IR. These values are listed in Table 2.

7 Results

Figure 5 shows the observations simulated with our pipeline. In particular, in the first row, we can see the images obtained simulating HST WFC3/IR observations, and in the second and third rows, those coming from JWST NIRCam. A preliminary visual inspection of our simulations reveals the detection of only a few clumps, consistently with what is observed in real images of the Cosmic Archipelago (Messa et al. 2025). The comparison between HST WFC3/IR and the short-wavelength (SW) filters of JWST NIRCam highlights how, despite the shorter exposure times, the better resolution and sensitivity of JWST allow the detection of some (observationally) faint clusters. Nevertheless, we note how sources that are at the detection limit in the SW NIRCam filters may be lost again in the long-wavelength (LW) ones. As a demonstration, the clump highlighted in Fig. 5 is undetectable in the HST simulations (with estimated signal-to-noise ratios S/N ~0.7 in F105W, ~0.4 in F125W, and ~0.3 in F160W), while it remains clearly visible in all short-wavelength NIRCam images (reaching S/N ~14 in F115W, ~8 in F150W, and ~10 in F200W). In contrast, its significance decreases in the long-wavelength NIRCam bands (S/N ~4 in F277W, ~5 in F356W, and ~1.5 in F444W). If we consider S/N < 3 as very low signal-to-noise ratios (Messa et al. 2025), it is clear that in the HST simulations the clump is undetectable, in the JWST LW bands it lies close to the detection limit, while it is clearly detectable in the JWST SW bands.

To assess the validity of our tool and enable a meaningful comparison with real observations, we compared the limit magnitude of our simulated images with that of the real images. By placing 100 circular apertures over dark regions of the simulated observation, we can define a numerical distribution of the background counts associated to each region. In the case of JWST/NIRCam F200W band, the standard deviation of those apertures corresponds to an AB limit magnitude mag_lim = 27.38. By repeating the same steps on the real data, we get the comparable limiting magnitude mag_lim = 27.42, supporting the validity of our methodology. From a qualitative comparison to the compact rest-UV sources seen in the Cosmic Archipelago (e.g., Fig. 2 in Messa et al. 2025), it is clear that the clusters we are simulating are fainter than their ‘real’ counterpart; indeed, while all our clusters have masses ≤5.8 · 10⁵ M_⊙ (see Fig. 2), the lowest mass reported in Messa et al. (2025) is 10⁶ M_⊙.

In Fig. 6, we show two panchromatic views of the simulated Cosmic Archipelago. In particular, in the two bottom panels, we show a color composite image created by combining three JWST NIRCam/SW channels (red F200W, green F150W, and blue F115W); and three JWST NIRCam/LW channels (red F444W, green F356W, and blue F277W). It is worth noticing that these images reproduce the same zoom-in over the Cosmic Archipelago seen with high-resolution in the first panel of Fig. 6. We use the red circle to trace back the position of the observed brightest clump to the full-resolution lensed image (top-left panel). Several individual clusters are imaged within the circle, as further revealed by the convolution to the F115W PSF (the top-right panel), providing an intermediate view of the lensed galaxy, prior to the inclusion of instrumental noise and sky background. In our simulated system, the fact that only a ‘single’ clump is detected within the red circle, is caused by one of the clumps being much brighter than the others. This highlights how faint clumps, or regions with low signal, can conceal a much more complex underlying structure. A more quantitative analysis will be carried out in future works.

Fig. 5 Simulated observations using HST WFC3 (the top three panels, for F105W, F125W, and F160W) and JWST NIRCam (the remaining panels for F115W, F150W, F200W, F277W, F356W, and F444W). In each subfigure, we show the same source at z = 6.145 lensed by MACS J0416.1-2403. Each panel uses its own flux range to optimize visibility of the source structure. The typical background level and exposure time for each filter are listed in Tables 1 and 2. Further details are reported in the text. The black squares in the bottom right corner of every subfigure show how the source galaxy appears when the gravitational lensing effect is not included. The unlensed and lensed sources are represented using the same pixel scale. A blue circle highlights the same clump in both the lensed and unlensed images, showing its corresponding position in the lens and source planes.

Fig. 6 RGB simulation of the Cosmic Archipelago (z = 6.145), which is in MACS J0416.1-2403. This is the same source shown in the right panel of Fig. 1, but including strong lensing effects by MACS J0416.1-2403. The unlensed source is assumed to be at z = 6.145 and distorted into a gravitational arc if placed near the cluster tangential caustic. Top left: a high-resolution RGB simulation (2.5 mas/pix), where the red, green, and blue channels correspond to wavelengths of λe = 5, 3, and 2 μm, respectively; top right: the same high-resolution image, but monochromatic, constructed using the magnitudes of the host and the clumps in the JWST/F115W filter and convolved with the JWST/F115W PSF to illustrate the effect of instrumental resolution; bottom left: an RGB image created by combining three JWST NIRCam short wavelength channels (red F200W, green F150W, and blue F115W), including realistic noise and sky background; bottom right: an RGB image created by combining three JWST NIRCam long wavelength channels (red F444W, green F356W, and blue F277W), including realistic noise and sky background. The background level and effective exposure time are based on real JWST observations of MACS J0416.1-2403 and are listed in Table 2. The red circle in each image highlights an example of a clump, illustrating how the instrument’s resolution is a limiting factor: while clumps appear as distinct structures in the HR image (see the zoom-in region shown at the bottom right of the figure), they become blended together in the JWST observations. For each image, the origin of the coordinate system coincides with the center of the cluster, located at (RA: 64.038142 deg, Dec: −24.067472 deg).

8 Abell 2744

This Section presents another example of a gravitational arc simulated with our tool. In particular, we focus on the so-called System 3 (Vanzella et al. 2022), a source at z = 3.98 lensed by Abell 2744. The mass model of this cluster is taken from Bergamini et al. (2023). The galaxy is assumed to have a total stellar mass M_*,tot = 2 × 10⁸ M_⊙, divided into a bulge with M_bul = 5 × 10⁷ M_⊙ and a disk with M_disk = 1.5 × 10⁸ M_⊙. Here, we assume the 50% of M_disk goes into stellar clumps, for an overall mass of 7.5 × 10⁷ M_⊙. Moreover, to model the bulge and the disk as Sérsic sources, we adopt the following parameters:

• Bulge: Sérsic index n = 4, effective radius R_e = 0.1″, position angle ϕ = π/8;

• Disk: Sérsic index n = 1, effective radius R_e = 0.3″, position angle ϕ = π/8.

Both components are centered at source position β_s,1 = 28.395195″ and β_s,2 = 40.2729467″, and at redshift z_s = 3.98; In order to build the spiral arms, we assume R_h = 0.3″, A = 1.0, N_a = 4, ϕ_d = 1.0, α = 10.0, $\varphi =\frac{3}{4}\pi$.

Regarding the stellar clumps, we use the same Schechter mass function (β = −2.0, δ = 1.0, γ = 1.0, x_min = 10⁴ M_⊙, x_max = 10⁷ M_⊙, and x_cut = 5 × 10⁶ M_⊙) and the same mass-size relation (Brown & Gnedin 2021) used in the previous example of the Cosmic Archipelago. We also consider the same parameters, such as the Sérsic index, position angle, and axis ratio, while the clumps’ positions have been modified to suit the new host galaxy.

To use the Yggdrasil models, we consider the same input parameters (Z = 0.004, Kroupa IMF), SFH and ages for each component considered: the bulge has a SFR which is constant over 100 Myr and is assumed to be 1 Gyr old, the disk has the same SFR, but it is 500 Myr old, while for the clumps the SFR is constant over 10 Myr and their ages span between 1 Myr and 500 Myr.

Here we simulate observations with the Advanced Camera for Surveys (ACS) onboard HST (angular resolution: 0.05 arcsec/pix), and the usual HST WFC3/IR and JWST NIRCam. Since HST WFC3/IR and JWST NIRCam have already been described in the previous example, we give some details only for HST ACS. The native pixel size of the ACS WFC camera onboard the HST is 0.0495″/pix. We generate the PSF models using the public PSF modeling tool Tiny Tim (Krist et al. 2011), simulating observations in the bands F435W, F606W, and F814W.

To define the sky background, we rely on the data of Abell 2744 from the HST proposals 11689 (PI: Dupke) for ACS simulations and 13386 (PI: Rodney) for WFC3/IR, while for JWST NIRCam, we consider the GTO program (PI Windhorst). Following the method already described in the previous Section, we list all the t_exp and bkg values in Table 3.

In Fig. 7, we show the observations simulated using the aforementioned combinations of filters and instruments. In this case we are able to detect many more clumps (>10 in F115W) than in the previous example, despite the initial conditions of the two simulations being quite similar. We point out that in this case the simulated galaxy is made of two mirrored images, thus each of the individual clusters is imaged twice. The observed difference with the Cosmic Archipelago is mainly driven by two factors: the difference in redshift and the one in magnification. The redshifts considered (z = 3.980 compared to 6.145 of the Cosmic Archipelago) correspond to a difference in observed magnitude of 1.1 mag. The lensing magnification of the system in A2744 ranges from a factor ~30 in the external regions (i.e., the left and right ends, in the lensed galaxy) to >100 near the critical line (toward the center of the lensed image, see Vanzella et al. 2022); on the other end, the largest magnifications achieved in the Cosmic Archipelago reach only a factor ~25 (Messa et al. 2025). For instance, the clump highlighted in Fig. 7 is detectable in all the simulated observations, both with HST (with estimated S/N ~11 in F435W, ~27 in F606W, ~15 in F814W, ~13 in F105W, ~4 in F125W, and ~8 in F160W) and with JWST (reaching S/N ~55 in F115W, ~45 in F150W, ~40 in F200W, ~22 in F277W, ~22 in F356W, and ~7 in F444W). In contrast to the case discussed in Sect. 7, here the significance never falls below the commonly adopted threshold of S/N = 3 (Messa et al. 2025), implying that the clump remains consistently detectable.

The top left panel of Fig. 8 shows the high-resolution simulation of System 3 at z = 3.98 in Abell 2744. In this Figure the red, green and blue channels correspond to rest-frame wavelengths λ_rf = 1.0, 0.6 and 0.4 μm, respectively. As done previously for Fig. 6, in addition to the lensed image at high-resolution, we provide the image of the convolution with the F115W PSF, without the inclusion of noise and background (top-right panel in Fig. 8). The middle row presents RGB images constructed from HST data: the left panel combines three ACS bands (red F814W, green F606W, and blue F435W), while the right one uses three WFC3/IR bands (red F160W, green F125W, and blue F105W). The bottom row shows RGB reconstructions from JWST NIRCam observations, using SW channels (red F200W, green F150W, and blue F115W) on the left and LW channels (red F444W, green F356W, and blue F277W) on the right. As already pointed out in the simulations of the Cosmic Archipelago system in Section 7, while the simulated high-resolution image clearly resolves numerous small-scale structures within the system, the lower resolution of the HST and JWST observations leads to significant blending, particularly among the faintest clumps. As a result, the observed emission is often dominated by the brightest components, potentially hiding a more complex distribution of stellar clusters.

Fig. 7 Simulated observations of the System 3 at z = 3.98 lensed by Abell 2744, using the HST ACS filters F435W, F606W, and F814W; the HST WFC3/IR filters F105W, F125W, and F160W; and the JWST NIRCam filters F115W, F150W, F200W, F277W, F356W, and F444W. Each panel uses its own flux range to optimize visibility of the source structure. The typical background level and exposure time for each filter are listed in Table 3. Further details, including the S/N analysis of the clump marked with a blue circle, are reported in the text. are reported in the text. The white squares in the bottom left corner of every panel show how the source galaxy appears when the gravitational lensing effect is not included. The unlensed and lensed sources are represented using the same pixel scale.

Fig. 8 RGB simulations of the System 3 at z = 3.98 in Abell 2744. Top left: a high-resolution RGB simulation (2.5 mas/pix), where the red, green, and blue channels correspond to wavelengths of λe = 5, 3, and 2 μm, respectively; top right: the same high-resolution image, but monocromathic, constructed using the magnitudes of the host and the clumps in the JWST/F115W filter and convolved with the JWST/F115W PSF to illustrate the effect of instrumental resolution; middle left: an RGB image created by combining three HST ACS bands (red F814W, green F606W, and blue F435W), including realistic noise and sky background; middle right: an RGB image created by combining three HST WFC3/IR bands (red F160W, green F125W, and blue F105W), including realistic noise and sky background; bottom left: an RGB image created by combining three JWST NIRCam/SW channels (red F200W, green F150W, and blue F115W), including realistic noise and sky background; bottom right: an RGB image created by combining three JWST NIRCam/LW channels (red F444W, green F356W, and blue F277W), including realistic noise and sky background. The background level and effective exposure time are based on the data of Abell 2744 from real HST and JWST observations and are listed in Table 3. The red circle in each image highlights an example of a clump, illustrating how the instrument’s resolution is a limiting factor: while clumps appear as distinct structures in the HR image (see the zoom-in region shown at the bottom right of the figure), they become blended together in the JWST and HST observations. For each image, the origin of the coordinate system coincides with the center of the cluster, located at (RA: 3.586257 deg, Dec: −30.400172 deg).

9 Conclusions

In this work, we have described the main phases and the fundamental ingredients for implementing a Python code to create realistic simulated observations of lensed, high-redshift, and star-forming galaxies. As examples, we present two simulations. The first one, based on the same source size and lensing configuration of the real arc Cosmic Archipelago (Messa et al. 2025), a source at z = 6.145 lensed by MACS J0416.1-2403, is used to describe the complete procedure of the tool. The second simulation follows a similar approach; however, in this case, we only present the final images and list the corresponding key parameters. The pipeline can be summarized in the following fundamental steps:

• Modeling of the host galaxy: the tool allows the user to choose which components to add to the source, including a disk, a bulge, and spiral arms. The brightness profiles of each component are described by Sérsic brightness distributions. As input parameters, we have to define each component’s stellar masses and sizes. Moreover, we have to define the source position on the source plane;

• Adding the stellar clumps: since we aim to simulate star-forming galaxies, our sources are designed to be populated by young stellar populations. The user selects which fraction of the host galaxy’s stellar mass goes into stellar clumps. Then, the tool constructs them based on the following steps: it samples their masses following a Schechter distribution, places them according to the underlying host’s brightness profile, and assigns them a radius that depends on their stellar mass;

• Simulating colors: the software can take as input parameters also the magnitudes in the bands that the user wishes to simulate. This step allows the simulation of different combinations of ages and SFHs for both the host galaxy and the stellar clumps. In the example described in this paper, we adopted a bulge of 1 Gyr and a disk of 500 Myr, while the clumps’ ages span between 1 and 500 Myr. Regarding the SFH, we assumed an SFR that is constant over 100 Myr for both components of the host galaxy, while for the clumps it is constant for 10 Myr. The tool relies on the Yggdrasil population synthesis code (Zackrisson et al. 2011);

• Simulating the lensing effects: by providing deflection angle maps, the pipeline can also include lensing effects. In the example presented here, we relied on the maps of MACS J0416.1-2403 from Lenstool (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009) by Bergamini et al. (2023);

• Including observational noises: observational noise and instrumental contamination can affect the previously generated noise-free images by providing the appropriate instrumental PSF, the flux from the sky background, and the proper exposure time, required to define the photon noise.

Following this framework, the tool can simulate realistic observations of lensed sources characterized by the presence of stellar clumps. In Figs. 5 and 6, we present examples of multiband simulations for HST and JWST. The reliability of our software has been tested and confirmed through a direct comparison with real images of the same sources, with a particular focus on the limiting magnitude, as discussed in Sect. 7. Sect. 8 describes and shows the simulations of the other example we considered in this work: the System 3 (z = 3.98) in A2744 (Vanzella et al. 2022).

The software presented in this work is a powerful and flexible tool designed to study the detectability of stellar clumps, particularly the smallest ones, in terms of mass and radius. It allows us to explore the fraction of clumps that appear blended in single and larger structures when, in reality, they consist of multiple smaller components. This approach enables us to probe the low-mass end of the stellar cluster IMF and to statistically investigate the physical properties of these systems.

One of the main strengths of the code lies in its modular design: all assumptions and parameters used to define the source galaxy, the gravitational lens configuration, the instrumental response, and the noise characteristics can be easily customized by the user. This makes the tool highly adaptable to a wide range of scientific goals and observational setups.

Beyond its immediate application in lensing simulations, the code can be used to generate realistic training datasets for machine learning algorithms aimed at detecting high-redshift clumps, or to benchmark and validate forward modeling frameworks such as GravityFM (Bergamini et al., in prep). Its flexibility and broad applicability make it a valuable resource for the community working on galaxy formation and gravitational lensing.

Data availability

The code used in this work is publicly available at https://github.com/irenemini/clumPyLen-master.

Acknowledgements

The research activities described in this paper have been co-funded by the European Union – NextGenerationEU within PRIN 2022 project no. 20229YBSAN – Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch. MM and EV acknowledge support from PRIN-MUR 2020SKSTHZ. MM was supported by INAF Grant “The Big-Data era of cluster lensing”. LM acknowledges the financial contribution from the PRIN-MUR 2022 20227RNLY3 grant “The concordance cosmological model: stress-tests with galaxy clusters” supported by Next Generation EU and from the grant ASI no. 2024-10-HH. 0 “Attività scientifiche per la missione Euclid – fase E.”

References

Appendix A Alternative Morphologies

Fig. A.1 Four examples of alternative morphologies for the source galaxy. Top left: a spiral galaxy composed exclusively of spiral arms, with a clump fraction Q = 80%; top right: an elliptical-like galaxy with a single ellipsoidal component, and Q = 80%; bottom left: a spiral galaxy including both spiral arms and a central spheroidal bulge, with Q = 1%; bottom right: a spiral galaxy composed of an exponential disk and a central spheroidal bulge, with Q = 90%. All models are shown at an angular resolution of 0.25 mas/pixel

Appendix B Key words

All Tables

All Figures

Fig. 1 (a) Example of a host galaxy consisting of an exponential disk with two spiral arms and a spherical bulge. Details on the adopted parameters are reported in the text. (b) Model of a host galaxy, made of an exponential disk, a spherical bulge, and two spiral arms, populated with Sérsic stellar clumps. Details on the adopted parameters are reported in the text. The galaxy center coincides with the center of the field of view. (c) RGB color image of the same source. Here, we assume the bulge is 1 Gyr old, while the disk is 500 Myr old. For the stellar clumps, their ages span in the range 1–500 Myr. In this case, the red, green, and blue channels correspond to effective wavelengths λe = 5, 3, and 2 μm, respectively (0.70, 0.42 and 0.28 μm, rest-frame, taking into account an assumed redshift z = 6.145). These images have an angular resolution of 0.25 mas/pix, which at a z ~ 6 corresponds to a physical scale of approximately 1 pc per pixel. As a result, our clumps, all of which have sizes between 1 and 10 pc, are resolved in the high-resolution images.

In the text

Fig. 2 Upper panel: magenta histogram showing the stellar mass function of the star clumps populating our simulated galaxy. Lower panel: clump sizes (in units of parsec) as a function of the clump stellar mass. Clumps are represented by magenta points. The underlying average size-mass relation given by Eq. (12) (Brown & Gnedin 2021) is shown in a blue solid line. At a given stellar mass, we assume the size distribution follows a lognormal distribution with a standard deviation of σ = 0.3.

In the text

	Fig. 3 Upper panel: histogram of the clump’s ages. Lower panel: clump stellar masses as a function of their ages in years.
In the text

Fig. 4 Mass map of galaxy cluster MACS J0416.1-2403 (Bergamini et al. 2023) at z = 0.396 (RA: 64.038142 deg, Dec: −24.067472 deg). This model is generated with the software Lenstool. The yellow lines represent the caustics of a source at z = 6.145, while the white lines are the corresponding critical lines on the lens plane. The magenta marker indicates the position of our source (β1 = −13.72 arcsec and β2 = 3.62 arcsec), and the three cyan markers show the corresponding multiple images on the lens plane: a stretched arc at the position θ1 = −35.57 and θ2 = 18.85, and two counter-images in θ1 = −42.66 and θ2 = 8.46 and θ1 = −20.16 and θ2 = 31.53, respectively. In this work, we focus on the gravitational arc.

In the text

Fig. 5 Simulated observations using HST WFC3 (the top three panels, for F105W, F125W, and F160W) and JWST NIRCam (the remaining panels for F115W, F150W, F200W, F277W, F356W, and F444W). In each subfigure, we show the same source at z = 6.145 lensed by MACS J0416.1-2403. Each panel uses its own flux range to optimize visibility of the source structure. The typical background level and exposure time for each filter are listed in Tables 1 and 2. Further details are reported in the text. The black squares in the bottom right corner of every subfigure show how the source galaxy appears when the gravitational lensing effect is not included. The unlensed and lensed sources are represented using the same pixel scale. A blue circle highlights the same clump in both the lensed and unlensed images, showing its corresponding position in the lens and source planes.

In the text

Fig. 6 RGB simulation of the Cosmic Archipelago (z = 6.145), which is in MACS J0416.1-2403. This is the same source shown in the right panel of Fig. 1, but including strong lensing effects by MACS J0416.1-2403. The unlensed source is assumed to be at z = 6.145 and distorted into a gravitational arc if placed near the cluster tangential caustic. Top left: a high-resolution RGB simulation (2.5 mas/pix), where the red, green, and blue channels correspond to wavelengths of λe = 5, 3, and 2 μm, respectively; top right: the same high-resolution image, but monochromatic, constructed using the magnitudes of the host and the clumps in the JWST/F115W filter and convolved with the JWST/F115W PSF to illustrate the effect of instrumental resolution; bottom left: an RGB image created by combining three JWST NIRCam short wavelength channels (red F200W, green F150W, and blue F115W), including realistic noise and sky background; bottom right: an RGB image created by combining three JWST NIRCam long wavelength channels (red F444W, green F356W, and blue F277W), including realistic noise and sky background. The background level and effective exposure time are based on real JWST observations of MACS J0416.1-2403 and are listed in Table 2. The red circle in each image highlights an example of a clump, illustrating how the instrument’s resolution is a limiting factor: while clumps appear as distinct structures in the HR image (see the zoom-in region shown at the bottom right of the figure), they become blended together in the JWST observations. For each image, the origin of the coordinate system coincides with the center of the cluster, located at (RA: 64.038142 deg, Dec: −24.067472 deg).

In the text

Fig. 7 Simulated observations of the System 3 at z = 3.98 lensed by Abell 2744, using the HST ACS filters F435W, F606W, and F814W; the HST WFC3/IR filters F105W, F125W, and F160W; and the JWST NIRCam filters F115W, F150W, F200W, F277W, F356W, and F444W. Each panel uses its own flux range to optimize visibility of the source structure. The typical background level and exposure time for each filter are listed in Table 3. Further details, including the S/N analysis of the clump marked with a blue circle, are reported in the text. are reported in the text. The white squares in the bottom left corner of every panel show how the source galaxy appears when the gravitational lensing effect is not included. The unlensed and lensed sources are represented using the same pixel scale.

In the text

Fig. 8 RGB simulations of the System 3 at z = 3.98 in Abell 2744. Top left: a high-resolution RGB simulation (2.5 mas/pix), where the red, green, and blue channels correspond to wavelengths of λe = 5, 3, and 2 μm, respectively; top right: the same high-resolution image, but monocromathic, constructed using the magnitudes of the host and the clumps in the JWST/F115W filter and convolved with the JWST/F115W PSF to illustrate the effect of instrumental resolution; middle left: an RGB image created by combining three HST ACS bands (red F814W, green F606W, and blue F435W), including realistic noise and sky background; middle right: an RGB image created by combining three HST WFC3/IR bands (red F160W, green F125W, and blue F105W), including realistic noise and sky background; bottom left: an RGB image created by combining three JWST NIRCam/SW channels (red F200W, green F150W, and blue F115W), including realistic noise and sky background; bottom right: an RGB image created by combining three JWST NIRCam/LW channels (red F444W, green F356W, and blue F277W), including realistic noise and sky background. The background level and effective exposure time are based on the data of Abell 2744 from real HST and JWST observations and are listed in Table 3. The red circle in each image highlights an example of a clump, illustrating how the instrument’s resolution is a limiting factor: while clumps appear as distinct structures in the HR image (see the zoom-in region shown at the bottom right of the figure), they become blended together in the JWST and HST observations. For each image, the origin of the coordinate system coincides with the center of the cluster, located at (RA: 3.586257 deg, Dec: −30.400172 deg).

In the text

Fig. A.1 Four examples of alternative morphologies for the source galaxy. Top left: a spiral galaxy composed exclusively of spiral arms, with a clump fraction Q = 80%; top right: an elliptical-like galaxy with a single ellipsoidal component, and Q = 80%; bottom left: a spiral galaxy including both spiral arms and a central spheroidal bulge, with Q = 1%; bottom right: a spiral galaxy composed of an exponential disk and a central spheroidal bulge, with Q = 90%. All models are shown at an angular resolution of 0.25 mas/pixel

In the text