© The Authors 2025

1. Introduction

Galaxy clusters are the most powerful gravitational lenses in the Universe. Comprising hundreds to thousands of galaxies, hot gas, and dark matter, these clusters distort space-time, bending the paths and amplifying the light of distant galaxies behind them. This well-known gravitational lensing effect results in a range of phenomena for background sources, such as altering their apparent position and shape, amplifying their brightness, and creating multiple images with different arrival times (Schneider et al. 1992).

Galaxies near the critical curves (hereafter CCs) of galaxy clusters form elongated images, commonly known as lensed arcs, which are typically magnified by a factor of $\mathcal{O}(10)$. CCs map onto caustics in the source plane with a well-defined dependency of the magnification, μ, with respect to the distance, d, to the caustic, $\mu ={\mu }_o/\sqrt{d}$, where for clusters, μ_o usually takes values of ∼10–30 when d is expressed in arcseconds. Sources with radius, R, (such as compact star-forming regions where R≲0.1″) within a lensed galaxy crossing a caustic, can achieve maximum magnifications of ${\mu }_{\max }={\mu }_o/\sqrt{R}\approx 50$. Even smaller objects, such as stars R≲10⁻⁸ arscec, can theoretically achieve extreme magnification factors, μ_max∼10⁶ (Miralda-Escude 1991). Such large magnification is possible if the distribution of mass is smooth and the lensed image lies extremely close to the CC. However, microlenses (such as stars from the intracluster medium) and small-scale structures within the cluster (dwarf galaxies, globular clusters, sub-haloes predicted in abundance by massive particle cold dark matter (CDM), or pervasive density modulations predicted in ultra-light particle CDM or Wave DM) perturb the gravitational potential, which reduces this maximum achievable magnification of the order of tens of thousands (Venumadhav et al. 2017; Diego et al. 2018; Weisenbach et al. 2024). This reduction in magnification near the CC is compensated (as demanded by flux conservation) by larger magnification factors around the microlenses and small-scale substructures farther away from the CC. These small-scale structures create a web of micro- and milli-caustics near the cluster caustic.

When a background star moves through this web of microcaustics the combined lensing effect from different microlenses produces many unresolved images (microimages) of the same star, which results in a combined magnification typically of the order of thousands, hence the term ‘extremely magnified stars’. The angular separation between these images is typically on the scale of microarcseconds (Venumadhav et al. 2017), hence far smaller than the resolution capabilities of current telescopes that operate with resolutions of tens of milliarcseconds for the most powerful optical and infrared telescopes. As a result, we observe a single image of the background star, which represents the combined magnification (or flux) of all the microimages.

Due to the relative motion between the source and the web of microcaustics, this magnification varies over time, and peaks when the source crosses one of the multiple microcaustics and then gradually decreases until the star is no longer detectable. These events can last from a few days to several years, depending primarily on the redshift of the source, its luminosity, the mass of the microlens, and more importantly the relative velocity and direction of motion with respect to the web of microcaustics.

Microlenses can also be made of compact dark matter objects, such as primordial black holes (PBHs) (Venumadhav et al. 2017; Diego et al. 2018; Carr & Kühnel 2020; Green & Kavanagh 2021; Palencia et al. 2024). However, current studies are consistent with a fully stellar microlens scenario (Oguri et al. 2018; Müller & Miralda-Escudé 2025). The discovery of Icarus (Kelly et al. 2018), a lensed star at redshift z = 1.49, marked the beginning of this field. Since then, O(100) of these lensed high-redshift stars have been discovered with both the Hubble Space Telescope (HST) and the James Webb Space Telescope (JWST).

Notable discoveries include Earendel, the farthest known star with a redshift of z≈6 (Welch et al. 2022a, b) (see Ji & Dai (2025) for caveat on the source size constraint), and the Dragon Arc in Abell 370, with 46 detected events in a single galaxy (Fudamoto et al. 2025). These discoveries open a new frontier in gravitational lensing, and enable the study of stars in the early Universe, shortly after reionization, their evolution, and potentially even the first generation of Population III stars (Zackrisson et al. 2024). They also offer new ways to investigate the nature of dark matter, either through its presence as compact objects within galaxy clusters (Diego et al. 2018; Oguri et al. 2018), or via small-scale fluctuations in the dark matter density, as predicted by wave dark matter models (Amruth et al. 2023; Broadhurst et al. 2025). High-cadence light curves of caustic-crossing events can serve as sensitive probes of dark matter microstructures formed in the early Universe (Dai & Miralda-Escudé 2020).

MACS J0416.1-2403 (Ebeling et al. 2001; Mann & Ebeling 2012), commonly referred to as MACS 0416, at redshift z = 0.396, is one of six galaxy clusters studied by HST as part of the Hubble Frontier Fields programme (HFF, Lotz et al. 2017). Later observed by JWST, MACS 0416 holds the record for the lens with the largest number of spectroscopically confirmed lensed galaxies (Rihtaršič et al. 2025). Several of the lensed galaxies cross the CC, hence with portions of these galaxies reaching maximum magnification.

This cluster has been extensively studied in the optical and infrared (IR) bands, with data from both telescopes being used by experts to derive best-fit lens models. Prior to JWST observations of this cluster, some of the earliest lensed stars at cosmological distances (z≳1) were discovered in this cluster by HST (Chen et al. 2019; Kaurov et al. 2019). Later on, Flashlights, a HST programme for taking extremely deep images of lensing galaxy clusters, identified about a dozen additional lensed stars (Kelly et al. 2022) in two background galaxies, named Spock and Warhol, both at redshift z≈1. Future deeper observations of these galaxies are expected to yield yet more microlensing events.

The modest redshift of z≈1 for Spock and Warhol, with a distance modulus of 44.10, increases the detection rates, since we can observe fainter stars. Typically, observations are biased towards the brightest stars, as they require a smaller (and more likely) magnification, in order to be detectable. This makes these two galaxies prime targets for studying the abundance of different star types at redshift z = 1, specifically probing the high-mass end of the initial mass function (IMF), as we are limited to the most luminous and hence brightest stars.

In this paper, we focus on Warhol and provide a detailed forecast for the number of detectable events across eight NIRCam filters. We also present the first spatial distribution forecast for the number of events. Previous studies have made similar predictions for Spock (Diego et al. 2024b). Warhol and Spock are not the only galaxies in MACS 0416 that host microlensing events. Mothra (Diego et al. 2023b), a lensed star at redshift z = 2.09, was also discovered in a galaxy behind MACS 0416; it provides interesting constraints on dark matter, as it demands the presence of a small millilens near the lensed star in order to explain its peculiar magnification.

Future dedicated observations of this galaxy cluster are expected to lead to many more detections, which would offer valuable insights into the nature of dark matter within galaxy clusters, as well as the formation and evolution of stars at high redshift. A larger number of events will allow us to test the results derived from this paper.

This paper is structured as follows. Section 2 details the data used in our analysis. In Section 3, we outline the methodology used to transition from data analysis to forecasting the rate of lensed events. The results, which cover both spatial and integrated event counts across eight NIRCam filters, are presented in Section 4. In Section 5, we discuss the prospects for detecting lensed high-redshift stars in MACS 0416 and similar arcs. Finally, our conclusions are summarised in Section 6.

Unless stated otherwise, we assume a flat cosmological model with Ω_m = 0.3, Λ = 0.7, and h = 0.7 (100 km s⁻¹ Mpc⁻¹). All magnitudes are given in the AB system. The Warhol arc is at redshift z = 0.94, and within the adopted cosmological model, 1 arcsec corresponds to ∼15.3 kpc.

2. Data

In this section, we describe the high-resolution HST and JWST mosaics of MACS 0416 that we used, as well as the scientific programmes from which the data were obtained. Finally we present the calibrated lens model used in this work.

To perform all photometric spectral energy distributions (SED) fitting, we used a combination of 11 photometric mosaic images of MACS 0416, using public data from HST (described further below) and the JWST Prime Extragalactic Areas for Reionisation and Lensing Science (PEARLS) program (Windhorst et al. 2023). We focus on a 24-arcsecond region centred on Warhol, using mosaics at a pixel scale of 30 mas per pixel.

For the lensing model of MACS 0416, we used the free-form lens model from Diego et al. (2024a) based on the latest JWST data. For the contribution of microlenses at the position of Warhol, we assumed a mean surface mass density of microlenses, Σ_* = 59.39 M_⊙ pc⁻² (Kaurov et al. 2019), which is subsequently modulated by the intra-cluster light (ICL). The rapidly declining brightness near the brightest cluster galaxy (BCG) suggests that the surface mass density varies along the arc, rather than remaining constant.

2.1. HST

We used HST mosaics that were constructed by PEARLS team members following the techniques first described by Koekemoer et al. (2011, 2013). These mosaics include public archival HST data on MACS 0416, specifically data from the Hubble Frontier Fields (HFF, Lotz et al. 2017), and the Beyond Ultradeep Frontier Fields and Legacy Observations program (BUFFALO, Steinhardt et al. 2020), as well as other public archival HST programs available in the STScI MAST Archive¹. For this paper, the mosaics that we used include all the HST Advanced Camera for Surveys (ACS) observations on this cluster in the F435W, F606W, and F814W filters, which cover a wavelength range of approximately 0.35 to 0.96 microns.

2.2. JWST

We used JWST mosaics constructed by PEARLS team members, using data from the PEARLS program obtained in eight filters from the NIRCam instrument: F090W, F115W, F150W, F200W, F277W, F356W, F410M, and F444W, which cover wavelengths from 0.795 to 4.981 microns. The corresponding magnitude depths, at which a point-source signal-to-noise ratio of 5 is achieved for fake sources injected (Williams, in preparation), are 29.7, 29.5, 29.5, 29.5, 29.5, 29.6, 29.0, and 29.3 m_AB, respectively. The HST and JWST mosaics that we used were all produced in a way that aligns them to a common astrometric grid, at a pixel scale of 30 mas per pixel. See Windhorst et al. (2023) for further details about the construction of the mosaics that we use here.

2.3. Lens model

We utilised the free-form code WSLAP+ (Diego et al. 2005, 2007) to derive the mass distribution and magnification of the macromodel. WSLAP+ is a hybrid method that allows the combination of weak and strong lensing constraints. For MACS0416, we only used strong lensing constraints. The method exploits the linearity (with mass) of the deflection field, which can be expressed as a linear combination of functions of the mass. It describes the mass distribution as the combination of two components: (i) a superposition of Gaussians located either on a regular grid or a adaptive grid, with additional Gaussians or increased resolution in specific regions of the cluster where higher mass is located, effectively acting as a prior for the mass distribution. This represents the smooth global component of the cluster lens. (ii) A compact component that accounts for the individual masses of cluster members, assumed to be proportional to their flux in a certain reference band. For this component, there is only one free parameter: the mass scaling factor or mass-to-light ratio of the cluster members. Different layers can be introduced to account for variations in the mass-to-light ratio among different types of galaxies within the cluster. The strong- and weak-lensing problem is then formulated as a system of linear equations:

$$\mathbf{\Theta }=\mathbf{\Gamma }\text{X},$$(1)

where the array X contains the free parameters: the grid points for the Gaussian approximation of the smooth component, the mass-to-light scalings, and the identified sources positions. The lensing observables, including the positions of strongly lensed sources (and when available the two shear components), are arranged in the Θ array. The Γ matrix is known and depends on the fiducial grid and redshifts of the lensed galaxies. The solution X is obtained by inverting a scalar function constructed from the system of linear equations. The lens model for MACS0416 is derived from a very large number of constraints (343) and is described in detail in Diego et al. (2024a).

3. Methodology

The objective of this paper is to derive a realistic forecast for the expected number of microlensing events on the Warhol arc. These transient events originate from stars within the lensed galaxy whose apparent magnitude, m, is typically well below the detection threshold, which makes them undetectable in most observational epochs. However, a temporary magnification boost, μ_micro, arises when the source approaches a microcaustic, which brings their apparent magnitude, m_obs, above the detection threshold:

$$m_{\mathrm{obs}}=m-2.5{\text{log}}_{10}({\mu }_{\mathrm{micro}})\lesssim m_{\varepsilon },$$(2)

where m_ɛ represents the detection threshold. The average number of detected stellar transient events depends on two main factors: the micro-magnification probability distribution function (PDF) and the background stellar population. The former describes the micro-magnification as a random variable that depends on the tangential component of the macrolens, μ_t, the radial component, μ_r, and the surface mass density of microlenses, Σ_* (see Palencia et al. 2024, for further details). Furthermore, local perturbations in the macrolens caused by millilenses can reshape the local macro-magnification map, thereby altering the magnification PDFs. The latter is influenced by various parameters, such as the age of the population, its metallicity, environmental constraints like dust attenuation, and the total stellar mass, among others. These parameters collectively shape the SED of the population, though some degeneracies exist, such as the dust and star formation history (SFH) degeneracy, where dusty and older populations both produce similar features with attenuated UV-rest-frame emission. Consequently, the number and properties of the source stars can be inferred from the total SED measured from the galaxy.

In this section, we present our methodology. We used a Markov chain Monte Carlo (MCMC) for the estimation of the stellar population parameters from SEDs across different regions of Warhol, which allowed us to construct a set of source stellar populations. We then integrated local magnification PDFs that reflect the peculiarities of the lens model, and applied them to the source stars to forecast the average expected number of events in various NIRCam filters at multiple detection thresholds. The results are presented and further discussed in Sections 4 and 5, respectively.

3.1. SED fitting

The first step in our analysis involves retrieving the SED of the Warhol arc. To achieve this, we carefully removed the contaminating foreground light from the galaxy cluster, specifically the ICL. We utilised SExtractor (Bertin & Arnouts 1996), a publicly available code that enables source detection, deblending of overlapping or nearby sources, and image background estimation. Additionally, we used SExtractor to delineate the region of pixels that corresponds to the Warhol arc. The next step involved subtracting the remaining part of the BCG emission, which SExtractor identifies as a source and treats independently from the image background. In our specific case, the BCG is almost directly in front of Warhol, so its subtraction is crucial for recovering an accurate representation of Warhol's SED. For this task, we relied on PetroFit (Geda et al. 2022), a Python package designed for tasks such as galaxy light profile fitting. We modelled the main BCG using a Sérsic profile (Sérsic 1963), which was then subtracted from the images. This removal of both foreground contaminants was performed separately for each filter image, and resulted in a set of contamination-free images across different filters for the Warhol arc.

The next step is to match the point spread function (PSF) of our images to prevent biases and ensure consistent flux measurements at the pixel level, thereby guaranteeing accurate photometry. Since broader PSFs spread the light of sources more, matching the images to the broadest PSF among all the filters ensures that we achieve bias-free photometry in each pixel of the arc. To accomplish this, we first need to determine the PSF for each filter. Foreground stars, being point-like sources (significantly smaller than the instrument's response), directly reflect the shape of the instrument's PSF in their images. We used two stars located north-west of Warhol (see Fig. 1) to construct an empirical PSF using Photutils (Bradley et al. 2024). Next, we calculated a convolution kernel for each image. This kernel is the inverse Fourier transform of the ratio between the Fourier transforms of the PSFs, specifically the narrower PSF over the broader one. We then convolved the images with narrower PSFs using the corresponding kernels, which resulted in a set of PSF-matched, contaminant-free images that were ready for accurate photometric estimation.

Fig. 1. MACS 0416 colour image showing a 24″×24″ region centred on the Warhol arc. Other prominent caustic crossing galaxies where microlensing transients have been identified are also marked. The colour image was created by combining filters from two instruments: NIRCam (F444W and F356W for the red channel; F277W, F200W, and F150W for the green channel; F115W and F090W for the blue channel) and WFC3 (F814W, F606W, and F435W for the blue channel). The yellow square outlines the central region analysed in this study.

The most accurate estimation of the background stellar population would involve fitting the SED at each pixel independently. However, as our method relies on MCMC to find the set of stellar population parameters that best reproduces the SEDs, it would be extremely computationally demanding. Instead, we grouped pixels with similar SEDs. To refine this process, we first masked the known brightest globular clusters (GCs) and potential transient events in our images; these appear as black dots in Fig. 2.

Fig. 2. Similar SED regions after clustering. Top-left panel: Colour image of Warhol. Bottom-left panel: regions of similar SED pixels on Warhol. Black dots indicate masked regions, globular clusters from the intracluster medium, or transient events whose SEDs do not represent the underlying stellar population. Second, third, and fourth columns: SEDs for each pixel in the different regions. The solid black line represents the median SED, while the dashed black lines indicate the 1σ contour. Violet crosses and red squares indicate the photometric measurements obtained using HST and JWST filters, respectively, all introduced in Section 2.

Globular clusters in MACS 0416 present challenges for observing microlensing transients behind them, but their local modifications to the macromodel and enhanced microlens density can significantly affect the rate of events around such millilenses. As for pixels containing the brightest transients in Warhol, we masked them because their SEDs do not represent the underlying stellar population; the transient's flux outshines the background, which prevents accurate SED determination. Both GCs and bright transients appear as local maxima, which makes their identification straightforward. We removed an area equivalent to twice the σ of an effective Gaussian centred at each maximum. Since the masked areas are relatively small compared to the size of the arc, we do not anticipate a significant underestimation of the total number of transient events expected in the arc, though recovering the full spatial distribution of events within the arc remains challenging.

Once transients and GCs are masked, we applied various methods to group pixels with similar SEDs. First, we performed a principal component analysis (PCA) to reduce the dimensionality of our dataset from 11 filters to three principal components, which capture 96.6% of the colour variability in our data. Additionally, we conducted a colour–colour analysis by plotting the differences between pairs of filters. Both PCA and colour–colour analysis can be combined with clustering techniques such as K-means, DBSCAN, or by applying linear cuts. We clustered the pixels based on two criteria: cuts performed in the reduced data spaces from PCA and colour–colour analysis, and spatial correlation to ensure that regions are spatially connected. The final regions identified through clustering, along with their median SEDs and 1σ contours, are shown in Fig. 2.

To estimate the photometric error, we used aperture photometry with different apertures around the Warhol arc. By analysing the observed variability in flux over area measurements around the arc, we can infer an approximate photometric error per square arcsecond. This error can then be extrapolated to the Warhol regions with known areas.

With accurate photometric SEDs for specific regions, calculated as the total contribution from all pixels within that region, and estimates of the flux errors for each filter, we can now fit these SEDs to reference models. For this purpose, we used the Flexible Stellar Population Synthesis (FSPS) library (Conroy et al. 2009, 2010; Conroy & Gunn 2010; Johnson et al. 2024). FSPS is a sophisticated tool for modelling the SED of stellar populations and offers flexible input options that allow users to specify various parameters, including IMFs, metallicities, masses, redshifts, ages, and more. It also allows for the incorporation of environmental conditions, such as dust attenuation, dust and nebular emission, or even emission from active galactic nuclei (AGNs) with different torus optical depths, which shape the AGN's contribution to the total SED.

One of the many outputs that FSPS generates is the photometric fluxes measured in the filters we are using. This enables us to compare model predictions to our data and explore the parameter space to find the best-fit parameters that produce an SED that minimises a χ² with our observations. We performed this fitting for all our regions, assuming different SFH models, which influence the quantity and types of stars in our selected regions.

The fitting process was carried out using modified packages from the PixedFit (Abdurro’uf et al. 2021) library, which utilises FSPS to generate the models. Additionally, we employed emcee (Foreman-Mackey et al. 2013), a Python implementation of the Affine Invariant MCMC Ensemble sampler (Goodman & Weare 2010), to explore the parameter space and find the best-fit solutions. Once the MCMC fitting was completed, we reconstructed the SFH of the stellar population based on the assumed model and best-fit parameters. This provides the star formation rate, i.e., the total amount of stellar mass formed per year from the beginning of the SFH to the current age of the galaxy at its redshift, z∼1. By integrating the SFH, we can determine the total amount of stars within various age bins. Since we used a region-integrated SED and assumed a homogeneous population of stars throughout the region, the mass distribution across the region was derived by weighting the total mass per age bin for that specific region by the flux distribution.

Once this process is completed for all six regions, we can produce a map that contains the mass distribution for each age bin (see Fig. 3 for the integrated value over all SFH). As the total mass distribution is weighted by the flux distribution, the mass obtained is biased by lensing magnification. To correct for this, we divided the mass in each pixel by its magnification. Since higher magnification corresponds to a smaller area in the source plane, this correction naturally reduces the mass in regions of high magnification. Conversely, higher magnifications make it easier to detect events in these regions. The corrected version of the mass distribution is shown in Fig. 3.

Fig. 3. Distribution of total stellar mass formed per pixel in Warhol throughout the entire SFH shown before and after accounting for magnification bias correction, in the left and right panels respectively. Brighter pixels may suggest a larger stellar mass, but magnification can make smaller stellar populations appear brighter, which would introduce a bias. To correct for this effect, we divide the stellar mass formed in each pixel by μm, which accounts for the change in source-plane area relative to the lens-plane pixel area due to magnification. Near the critical curve, the number of stars in the source plane imaged in this region is lower than in areas with lower magnification, but the magnification is extremely large, which enables the detection of individual fainter stars but within a reduced source area. One pixel corresponds to 32 mas, or 490 pc on the source plane if μ = 1.

3.2. Microlensing magnification

By using FSPS, we can not only reconstruct the SFH based on the best-fit parameters from the MCMC, but also generate a set of isochrones using different stellar libraries. In our case, we employed the MIST library (Dotter 2016; Choi et al. 2016; Paxton et al. 2011, 2013, 2015). FSPS outputs a file containing various stellar properties, such as age, type, luminosity, temperature, and surface gravity. All stars with the same age form an isochrone. Each star has a weight parameter (fraction by mass), which represents the fraction of stars with similar characteristics per solar mass formed. Since we know the total stellar mass formed per pixel and age, we can now determine the total number of stars in each pixel along with their luminosity, temperature, and other required stellar parameters. See Fig. 4 for an output example of FSPS isochrones.

Fig. 4. Set of isochrones outputted by FSPS from the MIST stellar library at different stellar ages in region 1 and best-fit metallicity Z = 0.31 Z⊙. Each colour represents the same population as it evolves with time, with older populations growing dimer and redder. The solid black line delimits the regions of the diagram in which a magnification larger than a factor of 104 would be necessary to observe the star in at least one of the eight NIRCam filters used in this work at a limiting magnitude of 30. The shaded region is thus the stars that cannot be observed in any filter with these conditions. The arrows mark the turn-off points, where stars leave the main sequence. The corresponding masses indicate the stellar mass at which this transition occurs for a subset of the isochrones shown in the plot.

The final step before incorporating microlensing is to determine the apparent magnitude of each star in the different NIRCam filters. This was achieved using the luminosity and temperature of each star. The apparent magnitude of a source is given by

$$m=-2.5{\log }_{10}\left(F_{\nu }\right)-48.6,$$(3)

where F_ν is the flux received in the selected filter, corrected for the filter response S(λ), and redshift dependence (Bessell & Murphy 2012):

$$F_{\nu }=(1+z)\frac{\int f_{\lambda }\left(\frac{\lambda }{1+z}\right)S(\lambda )\lambda \mathrm{d}\lambda }{\int S(\lambda )(c/\lambda )\mathrm{d}\lambda }.$$(4)

The specific flux per unit frequency interval is given by

$$f_{\lambda }(\lambda )=\frac{L}{4\pi D_{\mathrm{l}}(z{)}^2}\frac{\mathrm{SED}(\lambda )}{\int \mathrm{SED}(\lambda )\mathrm{d}\lambda },$$(5)

where L is the stellar bolometric luminosity, D_l(z) is the luminosity distance, and the stellar SED per wavelength is represented by a blackbody spectrum without loss of generality. Using Eqs. (3), (4), and (5), we can calculate the apparent magnitude of each star in the isochrones. For each star in the isochrone, we have their luminosity as shown in Fig. 4. We repeated this for each of the NIRCam filters.

Now that we have all the necessary information from the perspective of the sources (i.e. the stars), we can proceed to the second half of the process: microlensing. For this task, we gathered information from the macromodel and microlenses, and calculated the magnification statistics. Depending on these elements, a source will experience a flux boost, μ_micro, such that its perceived flux is μ_micro×F_f(ν). Following Eq. (3), its new apparent magnitude becomes

$$m_{\mathrm{micro}}=m-2.5{\log }_{10}\left({\mu }_{\mathrm{micro}}\right).$$(6)

The value of μ_micro is drawn from a magnification PDF estimated using M_SMiLe (Palencia et al. 2024). The PDF varies from pixel to pixel as the macrolensing and microlensing models change. M_SMiLe computes a semi-analytical magnification PDF based on three input elements: tangential macro-magnification, μ_t, radial macro-magnification, μ_r, and microlens surface mass density, Σ_*. The output at each pixel, i, is

$$p_i({\log }_{10}({\mu }_{\mathrm{micro}});{\mu }_{\mathrm{t},i},{\mu }_{\mathrm{r},i},{\mathrm{\Sigma }}_{*,i}),$$(7)

which we carefully mapped to p_i(μ) and normalised to unity. At first order, this probability density function peaks at μ_m≈μ_t × μ_r, while Σ_eff=μ_t×Σ_* controls the width of the PDF. As Σ_eff increases, higher μ_micro values become more likely, until Σ_eff≫Σ_crit, where the ‘more is less’ effect begins (Diego et al. 2018; Dai & Pascale 2021; Welch et al. 2022a; Palencia et al. 2024; Kawai & Oguri 2024). At this point, extreme values become less likely until the PDF converges to an attractor universal form, and no further differences are observed even at higher Σ_eff or μ_m. For a more detailed discussion, see Diego et al. (2018), Palencia et al. (2024).

We used the WSLAP+ mass distribution, κ, and the shear, γ, to estimate the macrolens model parameters. Both macro-magnification components (see Fig. 5) can be easily obtained as

$$\begin{array}{c}\hfill {\mu }_{\mathrm{t}}=1/(1-\kappa -\gamma ),\hfill \\ \hfill {\mu }_{\mathrm{r}}=1/(1-\kappa +\gamma ).\hfill \end{array}$$(8)

The Σ_* value was derived from ICL measurements from (Kaurov et al. 2019). We assumed a pivot value of 59.39 M_⊙/pc², which was modulated according to the ICL variation over the arc, as show in Fig. 6.

Fig. 5. Magnification model from WSLAP+ around Warhol. The left panel shows the smooth radial component. The middle panel shows the rapidly varying tangential magnification, with the divergent critical curve near the diagonal. The right panel shows the combined macromodel value, μm=μt×μr. Both μm and μt are shown in log scale for better visualisation. The region shown is the central 200 × 200 pixels in Fig. 1, enclosed in the yellow square around the Warhol arc. The grey line marks the position of Warhol.

Fig. 6. Map of Σ* modulated by the ICL, relative to a pivotal value of ∼60 M⊙/pc2. Open black circles indicate the positions of a random realisation of DM sub-haloes, with their location probability following the ICL distribution. The sub-haloes exhibit a mass distribution derived from the small-scale results of N-body simulations of dark matter within galaxy clusters at the MACS 0416 redshift. Larger masses are represented by larger circles. The masses of the sub-haloes that intersect the arc are indicated in the figure. This region shows the central 200 × 200 pixels in Fig. 1.

To obtain a more realistic picture of the lensing scheme, we also included GCs as millilenses. As seen in Diego et al. (2024a), one can expect ∼70 GCs in the 6″×6″ area around Warhol. Only a fifth of them will be detectable through photometry, but all of them will locally perturb the macromodel shown in Fig. 5. To achieve this, we randomly placed a single realisation of 72 sub-haloes in the 200 × 200 pixels around Warhol with masses (10⁶ M_⊙≲M_milli≲10⁷ M_⊙) following the distributions obtained by the MOKA simulations (Giocoli et al. 2012, 2016) in a galaxy cluster at the MACS 0416 redshift. For each position and mass, we placed a Gaussian mass with a full width at half maximum equal to one pixel (corresponding to 239 pc at the redshift of the lens), computed the deflection angle of each lens in the pixels of the field of view, added them linearly to those obtained from WSLAP+, and recalculated κ, and γ. Following Eq. (8), we obtained μ_t, and μ_r. In order to investigate the possible impact of compact dark matter (such as primordial black holes), we also assumed two additional scenarios, where we included an extra 3% and 10% of the total κ as microlenses that do not contribute to the ICL, i.e. compact dark matter. For reference, the critical density in the MACS 0416 and Warhol lensing system is 2970 M_⊙/pc², and κ=Σ/Σ_crit, where Σ is the surface mass density in the lens plane.

By inputting μ_t,i, μ_r,i, and Σ_*,i into M_SMiLe, we obtained the magnification probability density function in each pixel that we combined with the estimated number of stars, along with their apparent magnitudes, which we retrieved from our MCMC fitting. Integrating the PDF above a given minimum magnification, set by the star luminosity, dust attenuation, filter of choice, its redshift, and limiting magnitude of the observations, provides the probability of observing a star, j, in the corresponding pixel, i.

$$P_{i,j}={\displaystyle {{\displaystyle \int }}_{{\mu }_{\min ,j}}^{\infty }}{p}_i(\mu ;{\mu }_{\mathrm{t},i},{\mu }_{\mathrm{r},i},{\mathrm{\Sigma }}_{*,i})\mathrm{d}\mu .$$(9)

Summing over all stars, weighted by their expected abundance in that pixel, w_i,j, yields the expected number of events in that pixel:

$$n_i={\displaystyle \sum _{j=1}^{N_{\mathrm{star}}}}{P}_{i,j}{\mathit{w}}_{i,j}.$$(10)

The quantity μ_min represents the minimum magnification required for the observed magnitude, following Eq. (6), to be above a given detection threshold, such that m_micro<m_ɛ. Naturally, μ_min depends on the stellar apparent magnitude and the adopted threshold, μ_min=μ_min(m, m_ɛ). Summing over all pixels yields the average number of stars detected across the entire arc above the specified threshold:

$$\langle N\rangle ={\displaystyle \sum _{i=1}^{N_{\mathrm{pix}}}}{n}_i.$$(11)

The primary source of uncertainty in this calculation arises from the SFH. The uncertainties in the total stellar mass formed affect the number of stars susceptible to microlensing. By repeating the calculation performed with the median SFH, but now using the 16th and 84th percentiles obtained from the MCMC chains, we derived a pixel-level uncertainty in the number of events, σ_ni. The variance in the expected number of events is then given by

$${\sigma }^2\langle N\rangle =\sqrt{{\displaystyle \sum _{i=1}^{N_{\mathrm{pix}}}}{\sigma }_{n_i}^2}.$$(12)

For a given telescope pointing towards the arc with a limiting depth or magnitude threshold, there is an associated Poissonian uncertainty in the expected number of events, which can be approximated as

$${\sigma }_{\mathrm{Poisson}}=\sqrt{\langle N\rangle }.$$(13)

Finally, the total uncertainty in the number of events in the arc was obtained by summing all sources of error in quadrature:

$${\sigma }_{\mathrm{total}}=\sqrt{{\sigma }_{\langle N\rangle }^2+{\sigma }_{\mathrm{Poisson}}^2}.$$(14)

4. Results

After removing contaminating foregrounds and applying linear cuts in both the PCA dimensions and colour–colour spaces, we identified six distinct regions, as shown in Fig. 2. The second to fourth columns display the normalised SEDs of the pixels, alongside the median SED with the corresponding 1σ contours.

Fig. 7 presents the best SED fitting results obtained through MCMC analysis for each region, alongside the residuals for the 11 filters used in the photometric fitting (see Appendix A for median values and 1σ limits). The extracted SFHs, modelled as exponentially decaying SFHs, are shown in Fig. 8. By integrating the SFH over cosmic time, we derive the stellar mass formed in each region at different age bins. This mass is distributed across the region's pixels by normalising it according to the flux in each pixel, as illustrated in the left panel of Fig. 3. However, this flux is biased by lensing magnification, which necessitates a correction by dividing each pixel's mass by the macro-magnification predicted at its position within our lens model. The right panel of Fig. 5 displays this magnification, and the corrected stellar mass distribution, free from magnification bias, is shown in the right panel of Fig. 3.

Fig. 7. SED of Warhol regions and best-fit model from the MCMC realizations. The observed data points are showed in red and the solid black line represents the best model from the Flexible Stellar Population Synthesis code. The grey shaded area represents the 1σ contours. The residuals are showed in the bottom sub-panels for each region.

Fig. 8. SFHs of Warhol regions from the best fit in our MCMC realizations. The solid black lines represent the median SFR at each lookback time while the grey shaded contour shows the 1σ contour.

Summing across all regions, we find that the total stellar mass formed in the multiply lensed region of Warhol over its full SFH is $1.27{\pm }_{0.03}^{0.41}\times {10}^6$ M_⊙, or half that value if we consider that we have two connected images of the same background galaxy. The oldest stellar population, at 794 Myr (whose SFR is constant for the most recent times), is found in region 4, while the youngest, at 60 Myr, corresponds to region 3. All six regions exhibit sub-solar metallicities, with regions 1 and 2 having the highest values at 0.31 and 0.23 Z_⊙, respectively, while the remaining regions show metallicities at approximately 10% of the solar value.

After generating the isochrones that correspond to the stellar populations that best reproduce the observed photometric spectra in Warhol, we applied Eq. (9) to each star in the isochrones, and weighted them by the total stellar mass formed in that pixel at the corresponding isochrone age. By combining this information via Eq. (10), we obtained the probability density of microlensing events across the Warhol arc at different fiducial observational thresholds, as depicted in Fig. 9. This figure illustrates the distribution of events at various magnitude thresholds in the JWST NIRCam F090W filter. By integrating the expected event rates across all pixels using Eq. (11), we derive the expected number of events per pointing in F090W as a function of limiting magnitude, shown in Fig. 10. Figs. 9 and 10 incorporate a lens model that includes millilensing by sub-haloes, whose mass distribution follows the results obtained from the MOKA simulations (Giocoli et al. 2012, 2016) in a galaxy cluster at MACS 0416 redshift. The positions of these sub-haloes are marked by open black circles in Fig. 6, and their localised impact on the event distribution is evident in Fig. 9, particularly at shallower detection thresholds where their relative influence is greater. Fig. 11 presents the integrated number of expected events per pointing in F090W under the assumption of no millilenses. The overall impact of millilenses on the total event count is at the level of 1%, well below our uncertainties, which are dominated by Poisson statistics and SFH uncertainties from the SED fitting. However, their effect on the spatial distribution of events is more pronounced, as seen in Fig. 12. This figure shows the ratio of pixel-wise event rates for a 27-magnitude detection threshold in F090W. At the location of a millilens, the event rate is enhanced, whereas the surrounding regions exhibit a localised suppression. This results in a total integrated expectation that remains largely unchanged, as previously seen in Fig. 10. In the southern part of the arc, the combined influence of two neighbouring millilenses produces a similar effect to that observed in the northern region, albeit with a more pronounced suppression in the intermediate region between them.

Fig. 9. Spatial distribution of event probability predicted for different thresholds in F090W. Each pixel value represents the expected average number of detections above a certain threshold per pixel. The total number of pixels included in the galaxy mask and after masking out the globular clusters and possible events is roughly 5000. One pixel length corresponds to 32 mas, or 490 pc on the source plane if μ = 1. In shallow observations we only expect to see events at the critical curve. In deep observations, we expect to see events more uniformly distributed across the arc. Lower thresholds predict fewer events, which are asymmetrically distributed, favouring a negative parity. This is due to the increased Σ* closer to the BCG and the statistical properties of the negative parity regime. As the detection threshold increases, this asymmetry diminishes.

Fig. 10. Total number of events observed in Warhol for different thresholds in F090W, adding millilenses. Each value is the integrated value of all the pixels in the arc as shown in Fig. 9.

Fig. 11. Total number of events predicted in Warhol for different thresholds in F090W. Each value is the integrated value of all the pixels in the arc similar to Fig. 9.

Fig. 12. Ratio of expected stars magnified above 27.5 mag, assuming millilenses at the positions marked by black circles in Fig. 6, compared to the model without millilenses.

The different colours in Figs. 10 and 11 correspond to varying fractions of microlenses added to the constrained value from the ICL fitting. These additional microlenses, assumed to be compact dark matter objects such as PBHs, are parametrised as a fraction of the smooth surface mass density at the Warhol location. We find that the total expected number of events across the arc remains largely unchanged. Higher concentrations of compact dark matter result in a steeper slope, which indicates an increased sensitivity of event rates to limiting magnitude. This can be understood as due to the reduction in maximum magnification near the critical curves as more microlenses (PBHs) are added (Venumadhav et al. 2017). The effect should be compensated in the faint end of the luminosity function, where more events are expected farther away from the critical curve and associated to the increase in number of compact dark matter microlenses. This scenario is confirmed when we look at the spatial distribution of events when PBHs are added, which varies with the addition of extra microlenses, as shown in Fig. 13, where we plot the event ratio for a scenario with 10% additional κ in microlenses compared to a purely baryonic microlens model. We observe a reduction in the number of events near the critical curve due to the ‘more is less’ effect (Diego et al. 2018; Palencia et al. 2024), while the event rate is boosted at larger distances. The turnover point in this behaviour depends on the added fraction of microlenses: for an additional 10% in κ, we find turnover radii of 1.98″ and 0.33″ on the positive and negative parity sides of the arc, respectively, while for a 3% increase in κ, these values shift to 0.96″ and 0.15″, approximately half the previous distances.

Fig. 13. Ratio of expected stars magnified above 27.5 magnitudes. Expected number of events with the Σ* derived from the ICL over the events expected with a 10% increased Σ*. This mimics the effect of a hypothetical PBH (or other compact dark matter candidate) in the cluster medium.

We conclude that, after introducing millilenses into the macromodel and varying the fraction of microlenses in our model, both modifications imprint distinct patterns on the event rate distribution, while their effect on the integrated transient rate remains low or nearly undetectable. Millilenses cluster the event rate probabilities at their location in the lens plane, which simultaneously reduces the rate in the surrounding regions. In contrast, increasing the microlens abundance in Warhol, which is already large due to the baseline value extracted from ICL analysis, further decreases the event rate probability near high magnification regions such as the cluster CC, while enhancing the rates in areas with lower macromagnifications.

We predict an expected number of events per pointing at a limiting magnitude of 29.5 in eight NIRCam filters of order unity (See Table 1). This result is consistent with Yan et al. (2023), who report seven events in Warhol across four epochs, which yields a mean of 1.75 events per pointing at similar limiting magnitudes. The predicted event rates for all filters and magnitude thresholds are provided in Table 1.

5. Discussion

In this section we discuss the results obtained on the predicted microlensing event rates on eight NIRCam filters on Warhol, a lensed galaxy at redshift z = 0.94. First we explain the importance of the chosen models, both lensing models and stellar population models. We then evaluate the strengths and weaknesses of our forecasting pipeline.

Our results depend on two key factors: the lens modelling, both macro and micro, and the assumed stellar properties, including priors on stellar population parameters, dust attenuation models, and IMFs. The choice of these models directly impacts the predicted event rates. Here, we discuss our selections and their consequent effects on our results.

5.1. Macrolens

The choice of macrolens affects our predictions in two primary ways. Firstly, the local values of shear and convergence and their derivatives determine the macromagnification, which serves as an input for our code M_SMiLe, and influence the statistical properties of microlensing magnification. Secondly, different models predict slightly varying positions for the critical curve, thereby altering the spatial extent of the parity regimes, which in turn affects the statistical behaviour of micromagnification.

In the case of Warhol, its distinct morphology and bright knots, as shown in Fig. 2 top left panel, define critical points mirrored in both parity images. The critical curve precisely bisects the arc (Broadhurst et al. 2025), which leaves little margin for error in its location at the resolution of our telescopes. Thus, only minor variations in local macromagnification are expected at the arc's location, which implies a minimal impact on the total expected event rate. Moreover, MACS 0416 is one of the best-studied galaxy cluster lenses and hosts the largest number of spectroscopically confirmed sources and hundreds of multiply lensed images. This results in some of the most accurate lens models available for any galaxy cluster. However, different lens models may produce the critical curve at the right position yet have different slopes on the potential, which affects the scaling of the magnification with distance to the critical curve, μ∝d⁻¹. Overall, the entire problem is degenerate in Σ_eff=Σ_*×μ, so uncertainties in the macromodel and the ICL contribution can be reduced to uncertainties in Σ_eff.

However, this level of precision regarding the location of the critical curve in Warhol does not extend to all arcs. For example, Spock, another arc in MACS 0416, is located near a cluster member galaxy and lacks multiple critical points, which leads to greater uncertainty in its macrolens model and, consequently, in its predicted event rates (Li et al. 2025). Other arcs, such as the Dragon arc in Abell 370, which holds the record for the largest number of detected microlensed stars at cosmological distances (Fudamoto et al. 2025), are also more challenging to model than Warhol. These examples highlight the necessity of high-quality spectroscopic and photometric data in galaxy clusters to refine mass distribution models and improve constraints on lensing efficiency.

5.2. Micro and Millilens

The second major component of the lens model are the small-scale perturbers, including both millilenses and microlenses. As illustrated in Fig. 12, the spatial distribution of events is influenced by the presence of millilenses. These structures locally enhance the cluster's lensing efficiency by increasing both the macromagnification and, if the millilens is a globular cluster, the number of microlenses. This results in a higher event rate at the millilens position, albeit at the cost of a slightly reduced rate in the surrounding regions.

Given that tens of globular clusters and possibly DM sub-haloes are expected to intersect arcs like Warhol, accurately modelling their effects is crucial. Additionally, they can alter the local parity of the macro+milli model, which could potentially lead to long-duration events in negative parity regions, as is suspected in the case of Mothra (Diego et al. 2023a). While our study primarily focuses on the detection of events, light curves also play an essential role in a comprehensive analysis of microlensed stars in galaxy clusters.

Besides globular clusters, the actual CDM paradigm predicts the existence of non-luminous millilens-scale haloes (Dai et al. 2018, 2020; Williams et al. 2024; Ji & Dai 2025), which act as millilenses of 10⁶−10⁸ M_⊙. The spatial distribution of microlensing events can serve as a probe for these small-scale structures, as has been demonstrated for the case of wave DM (Broadhurst et al. 2025), particularly when combined with event duration data, although the latter is partially degenerate with the microlens mass.

Microlenses, a key component in our statistical framework M_SMiLe, also play a critical role in shaping the total event rate. A larger microlens population allows for events to occur farther from the critical curve, while simultaneously suppressing events closer to it due to the ‘more is less’ effect. This results in an overall reduction of event rates at shallower thresholds, whereas deeper pointings exhibit an increase in total expected events. These effects are illustrated in Figs. 10 and 13. We propose leveraging these effects to infer the dark matter abundances at microlens scales in galaxy clusters. This method is expected to be particularly effective in less crowded cluster regions, at higher redshifts where the critical curve moves outwards, or in merging clusters such as El Gordo (Diego et al. 2023a; Caminha et al. 2023; Frye et al. 2023), where the intracluster medium remains relatively pristine.

5.3. Star formation history

The SFH of a galaxy describes how much stellar mass has formed over time. Different galaxies exhibit a variety of SFHs, ranging from continuous star formation to single bursts or multiple episodes of star formation. To model SFHs, we generally adopt one of two approaches: a non-parametric method, where the galaxy's lifetime is divided into age bins of chosen widths and the total stellar mass formed in each bin is fitted independently, or a parametric approach, where the SFH is assumed to follow a predefined functional form governed by a set of parameters.

In general, non-parametric SFH modelling is more flexible and offers a better chance of capturing the true SFH due to a larger number of free parameters. However, having more free parameters (one per age bin) increases the computational complexity, particularly in MCMC-based models, where exploring high-dimensional parameter spaces can hinder convergence and can also increase the risk of over-fitting. Parametric models, on the other hand, are less flexible since they impose a predefined SFH shape, but they require fewer parameters, which significantly improves the efficiency of MCMC sampling. Nevertheless, incorrect assumptions about the SFH shape can introduce systematic biases in the inferred SFHs.

In this work, we adopt a simple, exponentially decaying SFH, also known as a τ-model, where the star formation rate evolves as

$$\mathrm{\Psi }(t)=A{e}^{-\frac{t-t_0}{\tau }}\mathrm{\Theta }(t-t_0),$$(15)

where t₀ represents the onset of star formation, τ controls the decay rate, and A is a normalisation factor.

Other commonly used parametric SFHs include the delayed exponential decay, double power-law models, constant star formation, multi-burst scenarios, and combinations of these models. The exponentially decaying SFH is a reasonable assumption, as it successfully reproduces the observed SED of many galaxies. Similarly, Li et al. (2025) conducted an extensive model comparison for Spock, and found that a non-parametric fit also favoured an exponentially decaying SFH. We tested different parametric SFH models in each of the six regions of Warhol: constant SFH, delayed exponentially decaying, double power law, and exponentially decaying, and we found a preference for the latter model in each region. However, since we did not perform non-parametric modelling or test a combination of parametric models, we may be missing a more complex SFH. It is possible that the true SFH is primarily characterised by an exponential decline but includes short bursts of star formation that are not captured in our analysis. Despite this limitation, analysing discrete regions within Warhol allows us to better account for variations in SFH across the arc.

Disentangling more complex SFHs would require higher quality data. We conclude that 11 photometric filters provide a reasonable understanding of the SFH in these arcs, but a more detailed approach could be achieved, particularly for the most recent SFH epochs, through spectroscopic analysis. Emission lines such as H_α and Pa_α are well-known tracers of recent star formation, and future studies should incorporate them to improve the accuracy of SFH reconstruction.

5.4. IMFs

The stellar IMF describes the number of main sequence stars formed as a function of stellar mass (Salpeter 1955). Many fundamental properties of stellar populations, including their SFH, metallicity, luminosity function, and stellar evolution, are directly influenced by the IMF. Measuring the IMF through resolved photometry has only been possible in nearby galaxies such as the Small and Large Magellanic Clouds (SMC and LMC) (Massey et al. 1995; Hunter et al. 1995, 1997; Sirianni et al. 2000, 2002; Da Rio et al. 2009), M31 (Weisz et al. 2015), and the Milky Way (Chabrier 2003). In all these cases, the observed IMF is typically modelled as a broken power law, with a commonly accepted slope of 2.3 for high-mass stars (>1.4 M_⊙) and a shallower slope for lower-mass stars. While the 2.3 slope for massive stars is widely agreed upon, different models vary in their treatment of the low-mass end of the IMF. The most commonly adopted IMFs include those of Kroupa (Kroupa 2001), Salpeter (Salpeter 1955), and Chabrier (Chabrier 2003), among others.

Transient events typically arise from the brightest members of the stellar populations, including O and B-type main sequence stars, red super-giants (RSGs), and possibly the tip of the red giant branch if the detection threshold is at a sufficiently high magnitude and the distance modulus of the arc is not too large. These stars are not only the most luminous but also the most massive, meaning that transient detections are primarily sensitive to the high-mass tail of the IMF. While lower-mass stars are not directly observed in these events, their presence influences the evolution of nearby massive stars, and due to their shallower IMF slope, they are significantly more abundant.

Whether the IMF evolves with redshift remains an open question (Davé 2011; Li et al. 2023). Complementary studies, such as Li et al. (2025), or Williams (in preparation), have investigated this issue, finding that at z∼1, a power-law slope of 2.3 is favoured. However, since our focus is on implementing M_SMiLe to assess its ability to spatially resolve transient event rate distributions rather than directly constraining the IMF, we adopt a Kroupa IMF (Kroupa 2001) for our analysis. Additionally, introducing the high-mass end slope of the IMF as a free parameter would significantly increase the computational time required for MCMC fitting.

5.5. Source plane modelling

In this work, we have modelled the Warhol arc in the lens plane. This approach was chosen over source-plane modelling to remain independent of any specific lens model. A lens-plane SED fitting yields a stellar population that directly represents the observed data in the arc, whereas source-plane modelling would require delensing the arc based on a particular lens model. Since gravitational lensing is achromatic, the SED of a lensed portion of the arc corresponds to that of the background source, with magnification only affecting the total brightness. This amplification would impact the inferred stellar mass, which could later be associated with a given lens model. In contrast, lens-plane modelling remains tied to the initially assumed model throughout the analysis.

One advantage of source-plane modelling is that it allows mirror-imaged pixels to be combined, effectively averaging out foreground contamination due to imperfect subtraction of the BCG, ICL, or unmodelled GCs in the lens plane. However, for an arc like Warhol, where the critical curve is well constrained and there is little dispersion among available lens models, the lens-plane modelling advantage has a limited impact. The dominant source of uncertainty in our estimates of the event rate remains the degeneracy in the SED fitting, which arises from the intrinsic complexity of stellar population synthesis models. Nevertheless, we argue that since our goal is to establish a general pipeline to estimate microlensing rates in lensed arcs, a lens-plane approach is more robust, particularly in cases like the Spock arc, where the lens models differ significantly, unlike in Warhol.

6. Conclusions

In this work, we presented a detailed forecast for the detection of highly magnified stars in the Warhol arc at redshift z = 0.94 lensed by the galaxy cluster MACS 0416, utilising a combination of strong lensing models, stellar population synthesis, and microlensing statistical techniques. Our analysis, based on photometric data from JWST and HST, provides insights into the spatial distribution and expected number of transient events in different filters and magnitude thresholds.

Our results indicate that the predicted microlensing event rates in the JWST NIRCam filters are in good agreement with current observations (Yan et al. 2023). We have characterised the spatial distribution of microlensing events in Warhol and shown that the highest event rates are not necessarily located at the positions of peak macromagnification, i.e., directly on the CC. Instead, we find that regions with milder macromagnification values report the highest number of expected events. This is consistent with the balance between the increased density of sources at lower magnifications or larger areas in the source plane and the statistical effects of microlensing in highly magnified regions, where the ‘more is less’ effect suppresses the number of observable transients.

Additionally, we have explored the impact of small-scale structures, including globular clusters, primordial black holes (PBHs), and millilenses, on the spatial distribution of microlensing events. Our results demonstrate that, while the total number of detected events remains largely unchanged, the presence of millilenses modifies the clustering of events, which enhances the microlensing efficiency in their immediate vicinity while slightly suppressing it in surrounding regions. This effect is particularly evident at shallow detection thresholds. The presence of dark matter in the form of PBHs and other millilens-scale sub-structures further alters the spatial pattern of events, creating distinct signatures that could be used to probe the small-scale distribution of dark matter within galaxy clusters.

Despite the limitations inherent to our modelling assumptions, our work demonstrates the power of microlensing as a statistical tool for studying lensed stellar populations at high redshift. Future spectroscopic follow-up observations, particularly targeting emission lines such as H_α and Pa_α, will be crucial in refining SFH models and breaking degeneracies between dust and age effects.

With ongoing and upcoming JWST observations, we expect a growing sample of detected lensed stars across multiple galaxy clusters. The methodology developed in this study provides a robust framework for interpreting these discoveries and advancing our understanding of both stellar populations at cosmological distances and the nature of dark matter on sub-galactic scales.

Acknowledgments

J.M.P. acknowledges financial support from the Formación de Personal Investigador (FPI) programme, ref. PRE2020-096261, associated with the Spanish Agencia Estatal de Investigación project MDM-2017-0765-20-2. J.M.D. acknowledges the support of project PGC2018-101814-B-100 (MCIU/AEI/MINECO/FEDER, UE) Ministerio de Ciencia, Investigación y Universidades. B.J.K. acknowledges funding from the Ramón y Cajal Grant RYC2021-034757-I, financed by MCIN/AEI/10.13039/501100011033 and by the European Union “NextGenerationEU”/PRTR. S.K.L. acknowledges support from the Research Grants Council (RGC) of Hong Kong through the General Research Fund (GRF) 17312122. This research is based on observations made with the NASA/ESA Hubble Space Telescope and James-Webb Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with programs GO-13496, GO-15936, and GO-16278, as well as by GTO-1176 and GTO-1208. The data are available at MAST: https://doi.org/10.17909/gvxn-w422 We acknowledge Santander Supercomputacion support group at the University of Cantabria who provided access to the supercomputer Altamira Supercomputer at the Institute of Physics of Cantabria (IFCA-CSIC), member of the Spanish Supercomputing Network, for performing simulations. We made use of the following software packages: Python, NumPy, SciPy, FSPS, PixedFit, Matplotlib, SExtractor, Petrofit, Astropy, Photutils, and DS9.

References

Appendix A: SED Fitting results

All Tables

All Figures

Fig. 1. MACS 0416 colour image showing a 24″×24″ region centred on the Warhol arc. Other prominent caustic crossing galaxies where microlensing transients have been identified are also marked. The colour image was created by combining filters from two instruments: NIRCam (F444W and F356W for the red channel; F277W, F200W, and F150W for the green channel; F115W and F090W for the blue channel) and WFC3 (F814W, F606W, and F435W for the blue channel). The yellow square outlines the central region analysed in this study.

In the text

Fig. 2. Similar SED regions after clustering. Top-left panel: Colour image of Warhol. Bottom-left panel: regions of similar SED pixels on Warhol. Black dots indicate masked regions, globular clusters from the intracluster medium, or transient events whose SEDs do not represent the underlying stellar population. Second, third, and fourth columns: SEDs for each pixel in the different regions. The solid black line represents the median SED, while the dashed black lines indicate the 1σ contour. Violet crosses and red squares indicate the photometric measurements obtained using HST and JWST filters, respectively, all introduced in Section 2.

In the text

Fig. 3. Distribution of total stellar mass formed per pixel in Warhol throughout the entire SFH shown before and after accounting for magnification bias correction, in the left and right panels respectively. Brighter pixels may suggest a larger stellar mass, but magnification can make smaller stellar populations appear brighter, which would introduce a bias. To correct for this effect, we divide the stellar mass formed in each pixel by μm, which accounts for the change in source-plane area relative to the lens-plane pixel area due to magnification. Near the critical curve, the number of stars in the source plane imaged in this region is lower than in areas with lower magnification, but the magnification is extremely large, which enables the detection of individual fainter stars but within a reduced source area. One pixel corresponds to 32 mas, or 490 pc on the source plane if μ = 1.

In the text

Fig. 4. Set of isochrones outputted by FSPS from the MIST stellar library at different stellar ages in region 1 and best-fit metallicity Z = 0.31 Z⊙. Each colour represents the same population as it evolves with time, with older populations growing dimer and redder. The solid black line delimits the regions of the diagram in which a magnification larger than a factor of 104 would be necessary to observe the star in at least one of the eight NIRCam filters used in this work at a limiting magnitude of 30. The shaded region is thus the stars that cannot be observed in any filter with these conditions. The arrows mark the turn-off points, where stars leave the main sequence. The corresponding masses indicate the stellar mass at which this transition occurs for a subset of the isochrones shown in the plot.

In the text

Fig. 5. Magnification model from WSLAP+ around Warhol. The left panel shows the smooth radial component. The middle panel shows the rapidly varying tangential magnification, with the divergent critical curve near the diagonal. The right panel shows the combined macromodel value, μm=μt×μr. Both μm and μt are shown in log scale for better visualisation. The region shown is the central 200 × 200 pixels in Fig. 1, enclosed in the yellow square around the Warhol arc. The grey line marks the position of Warhol.

In the text

Fig. 6. Map of Σ* modulated by the ICL, relative to a pivotal value of ∼60 M⊙/pc2. Open black circles indicate the positions of a random realisation of DM sub-haloes, with their location probability following the ICL distribution. The sub-haloes exhibit a mass distribution derived from the small-scale results of N-body simulations of dark matter within galaxy clusters at the MACS 0416 redshift. Larger masses are represented by larger circles. The masses of the sub-haloes that intersect the arc are indicated in the figure. This region shows the central 200 × 200 pixels in Fig. 1.

In the text

	Fig. 7. SED of Warhol regions and best-fit model from the MCMC realizations. The observed data points are showed in red and the solid black line represents the best model from the Flexible Stellar Population Synthesis code. The grey shaded area represents the 1σ contours. The residuals are showed in the bottom sub-panels for each region.
In the text

	Fig. 8. SFHs of Warhol regions from the best fit in our MCMC realizations. The solid black lines represent the median SFR at each lookback time while the grey shaded contour shows the 1σ contour.
In the text

Fig. 9. Spatial distribution of event probability predicted for different thresholds in F090W. Each pixel value represents the expected average number of detections above a certain threshold per pixel. The total number of pixels included in the galaxy mask and after masking out the globular clusters and possible events is roughly 5000. One pixel length corresponds to 32 mas, or 490 pc on the source plane if μ = 1. In shallow observations we only expect to see events at the critical curve. In deep observations, we expect to see events more uniformly distributed across the arc. Lower thresholds predict fewer events, which are asymmetrically distributed, favouring a negative parity. This is due to the increased Σ* closer to the BCG and the statistical properties of the negative parity regime. As the detection threshold increases, this asymmetry diminishes.

In the text

	Fig. 10. Total number of events observed in Warhol for different thresholds in F090W, adding millilenses. Each value is the integrated value of all the pixels in the arc as shown in Fig. 9.
In the text

	Fig. 11. Total number of events predicted in Warhol for different thresholds in F090W. Each value is the integrated value of all the pixels in the arc similar to Fig. 9.
In the text

	Fig. 12. Ratio of expected stars magnified above 27.5 mag, assuming millilenses at the positions marked by black circles in Fig. 6, compared to the model without millilenses.
In the text

	Fig. 13. Ratio of expected stars magnified above 27.5 magnitudes. Expected number of events with the Σ* derived from the ICL over the events expected with a 10% increased Σ*. This mimics the effect of a hypothetical PBH (or other compact dark matter candidate) in the cluster medium.
In the text