© The Authors 2025

1. Introduction

Individual stars in lensed galaxies at z ≳ 1 are now being routinely discovered in deep, repeated Hubble Space Telescope (HST) and James Webb Space Telescope (JWST) imaging of galaxy clusters (e.g., Kelly et al. 2018; Chen et al. 2019; Kelly et al. 2022; Yan et al. 2023; Meena et al. 2023; Fudamoto et al. 2025). These lensed stars must be intrinsically luminous, and be subject to extreme magnifications – from a combination of cluster macrolensing and stellar microlensing – to overcome the cosmological distance modulus and be detected as transient events (Miralda-Escude 1991; Oguri et al. 2018; Diego 2019). The two primary populations that constitute the lensed stars are: (i) the blue supergiants (BSGs, effective temperature of 10 000 − 30 000 K), which can be seen in blue filters that probe the rest-frame ultraviolet (UV) and/or optical (e.g., Kelly et al. 2018); and (ii) the red supergiants (RSGs, effective temperature of 3000 − 4000 K), which can be seen in red filters that cover the rest-frame near infrared (Diego et al. 2023). Bright O/B stars, or other classes of evolved stars such as asymptotic giant branches and Cepheid variables (Diego et al. 2024d) can also be detected as lensed transients, albeit much less frequently, given their lower luminosity.

Massive stars are very short-lived and so are only present in stellar populations at specific ages. For example, BSGs can only exist in stellar populations that have ages ≲10 Myr. By contrast, RSGs are evolved stars that take time (≳10 Myr) to evolve from their progenitors. The relative abundance of massive BSGs and RSGs has been well explored in nearby galaxies (where resolved photometry can be carried out) to characterize properties such as the age and metallicity of the stellar populations therein (Eggenberger et al. 2002; Dohm-Palmer & Skillman 2002). As lensing is achromatic and does not affect the color of lensed stars¹, the ratio between the detection rate of transients in a blue filter to a red filter, B/R, should also directly correlate with the abundance of the corresponding stellar populations in the lensed galaxy, at least, to the first order, as the maximum magnification RSGs can attain is lower owing to their larger radius. Diego et al. (2024c) first investigated how the B/R ratio is related to the relative abundance of RSG and BSG in a lensed galaxy under a simplified scenario, and found that a roughly equivalent abundance of BSG and RSG can reproduce the observed B/R ratio. In light of the strong and different dependencies between BSGs and RSGs with age, the B/R ratio can provide insights into the recent star formation histories (SFHs) of galaxies.

In this paper, we explore how the lensed star B/R ratio depends on star formation parameters and how this ratio can be used in combination with the spectral energy distribution (SED) of the galaxy to better constrain its recent SFH. Throughout this manuscript, we adopt the AB magnitude system (Oke & Gunn 1983), along with standard cosmological parameters: Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. Transient detection rate

As a demonstration, we use a lensed arc in the MACSJ0416 (z = 0.39) field known as “Warhol” (z = 0.94) in which multiple transients have been detected over the last decade (Chen et al. 2019; Kelly et al. 2022; Yan et al. 2023), as is shown in Fig. 1. We adopted the magnification map of Warhol from Palencia et al. (2025), where for each pixel in which light from Warhol was detected, we have a specific tangential magnification, −10^4.7 < μ_t < 10⁴, radial magnification, 1.7 < μ_r < 2.3, and surface mass density of stellar microlenses, 30 < Σ_★ < 115 M_⊙/pc².

Fig. 1. RGB image of the Warhol arc (z = 0.94). We show the transients detected in this arc as coral circles (reported by Yan et al. 2023; Williams et al. 2025) and green triangles (reported by Kelly et al. 2022), respectively. We highlight areas characterized as region 2 by Palencia et al. (2025) in red.

We began by simulating stellar populations that have different ages (t_age), metallicities (Z_★), and dust extinctions (A_V) by using the stellar population synthesis code SPISEA (Hosek et al. 2020) that adopts a MIST isochrone (Choi et al. 2016) and either a Kroupa (2001) (a power-law slope of 1.3 and 2.3 below and above 1.4 M_⊙, respectively) or a top-heavy (power-law slope of 1 for all stellar masses) stellar initial mass function (IMF). A representation of the age isochrones from MIST is shown in Fig. 2 for reference. In our simulations, each stellar population has the same total stellar mass of 10⁶ M_⊙.

Fig. 2. Subsets of MIST age isochrone (Choi et al. 2016) used in our SPISEA (Hosek et al. 2020) stellar population synthesis featuring a range of ages as shown in the legend. The metallicity of each of the subsets is reflected by their line style, as is shown in the bottom legend. None of them have any dust extinction, and the color represents the age of the isochrone, as is shown in the top legend. For each age isochrone, we also annotate some characteristic masses to indicate the mass-age dependency.

Following the methodology in Li et al. (2025) and Palencia et al. (2025), we computed the expected transient detection rate by calculating the expected number of stars in a stellar population that is brighter than the detection threshold when subject to cluster macrolensing and stellar microlensing. The probability of any stars having some magnification, p(μ), is characterized by μ_t, μ_r, and Σ_★ (Palencia et al. 2024). To ensure the unresolved objects in the lensed arcs are individual stars, we require that they be detected in the form of transients, such that our sample would not be contaminated by persistent unresolved objects that vary in brightness, such as star clusters (e.g., Li et al. 2024). We hence adopted the detectable-through-microlensing (DTM) limit introduced by Diego et al. (2024b) in our calculation, which restricts the individual background lensed stars to being undetectable when subjected to cluster macro magnification only, but becoming temporarily detectable when they are additionally boosted in magnification by microlenses. Mathematically, the detection rate for any pixel that has μ_t, μ_r, and Σ_★ is given by

$$\begin{array}{c}\hfill {\mathcal{R}}_f({\mu }_t,{\mu }_r,{\mathrm{\Sigma }}_{★};t,Z_{★},A_V)\propto \frac{1}{{\mu }_t{\mu }_r}{\displaystyle {\int }_{-\infty }^{m_{\sigma ,f}}{\int }_{{\mu }_{\text{min}}={10}^{-1}}^{{\mu }_{\text{max}}={10}^4}}·\\ \hfill \mathrm{d}\mu \mathrm{d}m{\prime }_{f}p(\mu ;{\mu }_t,{\mu }_r,{\mathrm{\Sigma }}_{★})N_f(m-2.5{\log }_{10}\mu ;t,Z_{★},A_V),\end{array}$$(1)

where m_σ, f denotes the detection limit in apparent magnitude, m, at the filter, f. By convolving p(μ; μ_t, μ_r, Σ_★) with the stellar luminosity function, N_f(m; t, Z_★, A_V), for a stellar population at a given age, metallicity, and dust extinction, we can compute the expected number of lensed stars and their brightness distribution after considering all possible lensing magnifications they can attain, μ, through the second integral. The choice of μ_min in this integral is motivated by the p(μ) resolution (Palencia et al. 2024), whereas the choice of μ_max is motivated by the maximum size of lensed stars (Li et al. 2025). The first integral determines the expected number of stars that are brighter than the detection threshold, given the combined lensing magnification from cluster and stellar microlenses in any random observation. The whole process is computed over every pixel in the image plane representing the arc, requiring us to divide the pixel area by the macro-magnification, μ_tμ_r, to account for the decrease in source plane area, and thus the decrease in the number of stars. We assume a constant surface brightness in the lensed arc for simplicity, where JWST F200W observation shows that the surface brightness of Warhol only varies by ∼40%, with diminished variations in regions over which transients are usually detected (near the critical curve).

We selected two characteristic filters that best capture the light from BSGs and RSGs at z ≈ 1. For the blue filter, we chose the HST WFC3/UVIS F200LP as it spans the rest-frame UV through the optical (100 nm to 540 nm) and has the highest sensitivity to BSGs; for the red filter, we chose the high-throughput JWST NIRCAM F200W as it covers the rest-frame near infrared (850 nm to 1.2 μm) where RSGs are the brightest. Owing to the high sensitivity of this filter, it is used in many deep programs. The filter response curves and characteristic blackbody spectra of BSGs and RSGs are shown in Fig. 3. We adopted a detection limit of 30 m_AB for F200LP, motivated by the HST Flashlights survey (Kelly et al. 2022), and 29.7 m_AB for F200W, motivated by the JWST PEARLS program (Windhorst et al. 2023; Williams et al. 2025). BSGs can sometimes be detected in F200W and, conversely, RSGs in F200LP. The likelihood of such detections, however, is low given the bandpasses of the respective filters compared with the blackbody spectra of the different stars. An example of the transient detection rate in these two filters is shown in Fig. 4. Those using other combinations of metallicity and dust are shown in Fig. A.1.

Fig. 3. Filter response curve of HST F200LP (blue shade) and JWST F200W (coral shade), and the blackbody spectrum of two stars at 30 000 K (blue curve) and 4500 K (red curve) – characteristic temperature of BSGs and RSGs, respectively – shown at the rest frame and observed frame (z = 0.94) wavelengths. One can see that the two filters best capture the brightest part of the blackbody of the two classes of stars.

Fig. 4. Transient detection rate per unit mass for stellar populations at different ages with a metallicity of 0.5 Z⊙ and AV = 0, calculated with Eq. (1) under the magnification map of Warhol from Palencia et al. (2025). Blue curves show the detection rate in the HST F200LP filter (detection limit of 30 mAB); red curves show the detection rate in the JWST F200W filter (detection limit of 29.7 mAB). Solid and dotted lines show the case of a Kroupa IMF and a top-heavy IMF, respectively.

The transient B/R ratio for any given stellar population is then just the ratio between the detection rate summed over all the n lens model pixels that represent the arc in the two filters:

$$\begin{array}{c}\hfill B/R={\displaystyle \sum _i^n}{\mathcal{R}}_{F200LP}({\mu }_{t,i},{\mu }_{r,i},{\mathrm{\Sigma }}_{★,i})/{\displaystyle \sum _i^n}{\mathcal{R}}_{F200W}({\mu }_{t,i},{\mu }_{r,i},{\mathrm{\Sigma }}_{★,i}).\end{array}$$(2)

We show the dependence of the B/R ratio on age for stellar populations that have different metallicities, Z_★, and for each Z_★ subject to different dust extinctions, A_V, in Fig. 5.

Fig. 5. Expected B/R ratio calculated based on the F200LP and F200W transient detection rate in Warhol arc shown in Fig. A.1, with a Kroupa IMF. Each of the panels shows the case for one of the dust extinctions (upper row) or one of the metallicities (lower row), as is indicated by the subtitles. The colour of each curve indicates the metallicities (upper row) or dust extinction (lower row). The dashed gray line denotes B/R = 1 for reference. The fluctuations in ages earlier than ∼50 Myr are significant and arise from intrinsic variation due to stellar evolution, whereas those later than ∼50 Myr arise from sampling uncertainty given the low transient detection rate at these ages. We do not show the error bars for visualization purposes.

3. Properties of the B/R ratio

As is shown in Fig. 3, the transient detection rate per unit mass at z ≈ 1 in both the blue (HST F200LP) and red (JWST F200W) filters is highest for star formation episodes over the past ∼10 Myr (notice that the peak in transient detection rate is slightly delayed in the red filter, as massive stars take time to evolve into RSGs). Beyond ∼10 Myr, the transient detection rate decreases quickly with increasing age as most BSGs and a significant fraction of RSGs (≳10⁴ L_⊙) only survive up to ∼10 Myr, as is shown in Fig. 2. After this age, better still, the transient detection rate is dominated by less massive evolved stars. These stars, however, contribute much less to the transient detection rate than BSGs or RSGs, given their lower intrinsic brightnesses, as is shown in Fig. 2, requiring a higher (less likely) magnification to be detected as a transient.

For a coeval stellar population that has an age > 10 Myr, the transient detection rate decreases more slowly with age in the R filters than in the B filters. This behaviour translates to a decreasing B/R ratio with age, as is shown in Fig. 5. As an example, for a stellar population at the youngest age of ≲3 Myr, the detection rate is significantly higher in the blue filter, yielding a high B/R ratio of ∼10^2 − 3. At such young ages, the DTM stars are mostly BSGs, as is demonstrated by the blue isochrone in Fig. 2, while less massive stars need more time to evolve into the RSG phase. Although BSGs are so bright that they can sometimes also be detected in the red F200W filter, they primarily contribute to the detection rate in the blue filter, as is shown in Fig. 4. On the other hand, red DTM stars (most likely, RSGs) only become appreciable at a slightly older age of ≳5 Myr (the cyan isochrone in Fig. 2), corresponding to the time taken by intermediate-mass stars to evolve away from the main sequence. The B/R ratio thus drops quickly once massive stars start to enter the RSG phase. At even older ages of ≳5 − 40 Myr (depending on the metallicity and dust extinction), the most massive (and thus brightest) BSGs are gone, with few of the dimmer BSGs surviving, as is shown by the green isochrone in Fig. 2. During this time, the B/R ratio falls below unity, as is reflected in Fig. 5, as the only remaining DTM stars are relatively dim RSGs alongside other red evolved stars. This trend continues at even older ages, at which the B/R ratio fluctuates owing to effects related to metallicity and dust extinction, as is explained next.

As can be seen in Fig. 5, the dependence of the B/R ratio on metallicity is weaker than the dependence on age, but still could promote or demote the B/R ratio up to factors of ∼10 − 100 depending on the age of the stellar population. There are a couple of factors that determine how metallicity affects the B/R ratio. The first is the metal absorption in stellar atmospheres, for which there are more metal absorption lines in the B than R filters, such that an increase in metallicity would generally decrease the B/R ratio, as can be seen in Fig. 5. Furthermore, as is shown by the isochrones in Fig. 2, stars with high metallicities (indicated by dotted lines) are cooler (and thus redder) compared with those with a lower metallicity (indicated by dashed lines). Another complicated effect would be the role of metallicity in stellar evolution, and therefore the time that stars of different masses spend at different stages, which is more age-specific. For example, at super-solar metallicity (right-most panel, Fig. 5), the increase in stellar opacity owing to metal absorption inflates the star and leads to increased mass loss via stellar wind (Vink et al. 2001). With the outer envelope stripped away, the massive stars cannot evolve into RSGs and remain as BSGs, leading to a dip in the transient detection rate in the R filter at ∼10 Myr, explaining the bump in the B/R ratio at this age. These effects depend on the isochrone adopted (Choi et al. 2016), and we do not analyze all the trends in further detail in this work.

Compared with metallicity, dust extinction has an even more minor effect on the B/R ratio. Increasing the dust extinction results in more blue light being absorbed, making the observed stellar population redder in general. This increase leads to the detection of more red transients and suppresses the B/R ratio at any age. Such an effect is more appreciable after an age of ∼10 Myr – the DTM stars at such age are much dimmer, and the dust extinction would make them even dimmer in the B filter, further lowering the B/R ratio.

Last but not least, we found that the B/R ratio varies relatively weakly with the two adopted IMFs. We do not show the corresponding results for a top-heavy IMF in Fig. 5, as the changes induced compared with a Kroupa IMF (as is shown in Fig. 5) are much weaker than the observed variation with age. A shallower IMF gives rise to more blue transients at ages of less than 3 Myr, with a fractional change in the B/R ratio of ≲20%. This difference in the B/R ratio for the different stellar IMFs diminishes (≲3%) at older ages, as Kroupa has a similar slope as the top-heavy IMF for low- to intermediate-mass stars.

4. Insights into star formation

In Sect. 3, we showed how the B/R ratio changes dramatically with the age of a stellar population, and to a lesser extent metallicity, and found that this ratio is quite insensitive to the assumed IMF (whether Kroupa or top-heavy) as well as dust extinction. As a consequence, the B/R ratio can be used as a powerful probe of the recent SFH of galaxies that are multiply lensed by galaxy clusters. In this section, we show that the best-fit synthetic stellar population to the SED of the Warhol arc predicts incompatible B/R ratios for this galaxy, and vice versa, no matter whether we adopt a single burst or a more extended duration of star formation. We then show how fitting to both the B/R ratio and SED simultaneously improves our understanding of the recent SFH of this galaxy, predicting a quite different stellar mass formed over the past ∼50 Myr compared to that implied when fitting either constraint alone.

4.1. Single starburst

As a proof of concept, we begin with the simplest possible case, in which there is only one coeval population of young stars. The observed B/R ratio of the Warhol arc is 0.33 ± 0.26 (Kelly et al. 2022; Yan et al. 2023, summarized by Williams et al. 2025). In the case of a single stellar population, we deduce from Fig. 5 that this arc must have experienced a starburst ∼5 − 40 Myr ago, depending on the assumed metallicity. As we show in Fig. 6, however, the observed SED (black curve, from Palencia et al. 2025) is not consistent with a starburst over the aforementioned age range, for which the blue curve shows the predicted SED for an age of 12 Myr (that which minimizes differences with the observed SED, but for which the reduced χ² ≈ 100). Conversely, fitting a single starburst generated by SPISEA to the observed SED of Warhol, we find a best-fit age of ∼350 Myr (requiring 0.3 Z_⊙ and A_V = 0.9) for which the reduced χ² is 1.8 (red curve in Fig. 6). The predicted B/R ratio at this age, however, is on the order of 10⁻⁷, in dramatic contradiction with the observed B/R ratio.

Fig. 6. SED of Warhol arc (black curve, from Palencia et al. 2025). The best fit (red curve) corresponds to a single star burst model at an age of ∼350 Myr with Z★ = 0.32 ± 0.12 Z⊙ and AV = 0.91 ± 0.04. This model predicts B/R = 10−7, orders of magnitude from the observed B/R = 0.33 ± 0.26. Among all the models that can reproduce the observed B/R ratio well, we show the best fit SED as a blue curve. This model corresponds to a single star burst at an age of ∼80 Myr with Z★ = 0.8 ± 0.4 Z⊙, and AV = 0.5 ± 0.3. We adopted a Kroupa IMF in both models.

The corresponding Markov chain Monte Carlo (MCMC, see details in Appendix B) sampling of the constraints imposed by either the B/R ratio or SED alone on the properties of a single starburst is shown in Fig. B.1. As can be seen, regardless of whether we adopt a Kroupa or top-heavy IMF, there is a ∼4σ tension between the age of the starburst inferred from the B/R ratio (blue contours) or SED (red contours) alone. This contradiction clearly shows that the Warhol arc must have a more complicated SFH so as to reproduce the SED and B/R ratio simultaneously.

4.2. Continuous star formation

We now consider a more complicated SFH, in which the star formation rate (SFR) is allowed to be freely fit in five different age bins spanning 1 − 5 Myr, 5 − 10 Myr, 10 − 50 Myr, 50 − 250 Myr, and 250 Myr − 1 Gyr. We adopted a Kroupa IMF, having shown earlier that a top-heavy IMF makes little change to the predicted B/R ratio.

We show the best-fit SED in Fig. 7, demonstrating a 2σ tension between the resulting predicted and observed B/R ratios ($1.0_{-0.1}^{+0.2}$ versus 0.33 ± 0.26) when the SED is used as the sole constraint on the SFH. Similarly, we also demonstrate the resultant tension between the predicted (blue curve) and observed (black curve) SED for the best-fit SFH model when constrained by the B/R ratio alone. As the B/R ratio does not depend on the stellar mass formed, the model SED plotted in blue is that which best reproduces the measured SED given the best-fit SFH to the B/R ratio. Evidently, increasing the complexity of the SFH does not resolve the tension between best-fit model SFHs as constrained by the SED or B/R ratio individually.

Fig. 7. Same as Figure 6, but we use a nonparametric SFH model instead of a single star burst to fit the data. Again, the SFH that is constrained by the SED (red curve) does not predict the B/R ratio well (predicted $1.0_{-0.1}^{+0.2}$, versus observed 0.33 ± 0.26), while that is constrained by the B/R poorly predicts the SED (blue curves). Notice that the B/R ratio is not sensitive to the underlying stellar mass; thus, the SED inferred from B/R deviates significantly from the observed SED (dotted blue curve). We enforced a normalized version of the same SED (solid blue curve), whereby the relative SFR in each age bin of the underlying SFH is invariant, to acquire the best fit to the observed SED. Nevertheless, imposing both the SED and B/R as constraints infers a solution that fits both of the constraints to a sensible degree, which signifies the importance of including the B/R ratio as an extra constraint to better capture the true SFH. The best-fit solution has a metallicity of 0.44 ± 0.21 Z⊙ and a dust extinction of 0.91 ± 0.07, respectively.

In light of the persistent tension, we conducted another MCMC analysis to test whether there is any model SFH, parameterized in the manner described above, that can provide an acceptable fit to both the B/R ratio and SED simultaneously. We were able to find such a solution, for which the model SFH is shown in Fig. 8, to be compared with the best-fit model SFHs when fitting for either the SED or B/R ratio alone. The corresponding predicted SED is shown by the gray curve in Fig. 7 with the MCMC result shown in Fig. B.2. The metallicity inferred in this solution is 0.44 ± 0.21 Z_⊙ with dust extinction of 0.91 ± 0.07. As can be seen, the model SED closely resembles the observed SED (black curve). Furthermore, despite the different implied SFHs as shown in Fig. 7, the model-predicted SED constrained by the SED and B/R ratio simultaneously is very similar to that constrained solely by the SED. At the same time, when constrained by both the SED and B/R ratio simultaneously, the B/R ratio is now also reproduced to well within 1σ. When constrained by both the SED and B/R ratio simultaneously, the stellar mass formed over the past ∼50 Myr is $\sim 12{300}_{-4500}^{+7200}{M}_{\odot }$ ($0.9_{-0.3}^{+0.5}\%$ of the total stellar mass), two orders of magnitude higher than that implied by the corresponding SFH models constrained by the SED alone ($\sim {260}_{-150}^{+580}{M}_{\odot }$, $0.{02}_{-0.01}^{+0.04}\%$ of the total stellar mass) with a ∼3σ tension.

Fig. 8. Corresponding best-fit SFH models of the Warhol arc that generate the SED and B/R ratio shown in Fig. 7 with the band representing the 1σ tolerance from MCMC analysis. One can see that the inclusion of the B/R ratio as a constraint on top of the SED gives rise to a different SFH, with more star formation over the past 10 − 50 Myr (factor of ∼60 with ≳3σ significance) and less star formation over the last 50 − 250 Myr (factor of ∼60 with 2σ).

We also compared the Bayesian information criterion (BIC) to see which of the two aforementioned models is preferred. The BIC is given as

$$\begin{array}{c}\hfill BIC=k\text{ln}(N)-2\ln (L_{\max }),\end{array}$$(3)

where k is the number of model parameters, N the number of data points, and L_max the maximum likelihood of the given model. The best-fitting model that is constrained only by SED has a BIC of 45, and the best-fitting model that is constrained by both SED and B/R ratio has a BIC of 30. This means that there is very strong evidence (difference in BIC > 10) that the latter model is better at reproducing the observations.

The main difference between the SFH constrained by the SED alone and by the SED and the B/R ratio simultaneously is that the latter one has a significantly higher SFR (factor of ∼60, ≳2σ) in the third age bin of 10 − 50 Myr, and a lower SFR in the fourth age bin of 50 − 250 Myr (2σ). Since a certain number of UV bright stars (only formed in the last ∼10 Myr) are required to fit the rest-frame UV part of the SED, the only way to reduce the overpredicting B/R ratio is to form more red stars that boost the detection rate in the R filter. The age bin 10 − 50 Myr contains a relatively high abundance of RSGs compared to BSGs, whereas an increased SFR in this bin suppresses the B/R ratio (cf. Fig. 5, middle panel) to fit the observation without sacrificing the quality of fits towards the SED. On the other hand, the SFR is lower in the fourth bin, possibly a result of compensating for the increased total fluxes due to the increased SFR in the third age bin. The transient detection rate between ages of 50 − 250 Myr is expected to be very low (cf. Fig. 4), such that the variation in the SFR in this bin would not affect the inferred B/R ratio at all. This result also hints at the recent star formation in the Warhol galaxy perhaps being rather bursty in nature, with the B/R ratio possibly constraining the burstiness.

4.3. Spatial variation in the star formation history

One of the benefits of using the B/R as a constraint on the SFH is that it has extremely high spatial sensitivity down to the pixel level. The Warhol arc has varying colors over different regions, which is a strong indicator of varying SFHs in recent times over different parts of the galaxy. Among the six regions of Warhol, as characterized by their photometric color (Palencia et al. 2025), only region 2 (as is highlighted in Fig. 1) has a transient detection in both the B and R filters with a B/R ratio of 0.3 ± 0.4. Since we only have a measurement of B/R at region 2, we cannot investigate how the inclusion of the B/R ratio as a constraint imposes different SFHs compared with those constrained only by the SED in different parts of the arc; instead, we can compare the local SFH confined to region 2 with the globally inferred SFH.

We applied the same calculation described earlier in Sect. 4.2 to region 2 and show the results in Fig. 9 and Fig. 10, where the MCMC result can be found in Fig. B.3. When SED is the only constraint, the B/R ratio is poorly predicted ($5.7_{-1.3}^{+0.3}$ versus 0.3 ± 0.4) with worse tension (> 3σ) than the case for the entire Warhol arc. Again, when B/R is incorporated as an additional constraint, the tension is resolved. One can tell from Fig. 10 that the SFH inferred is very different when the SED alone (red) is used as a constraint, or the SED and B/R ratio are used as constraints together (gray). The recent SFR (1 − 10 Myr) is suppressed by a factor of ∼10 with ∼1 − 2σ, and the SFR at 10 − 50 Myr ago is boosted by a factor of ∼20 with > 3σ upon the inclusion of B/R ratio as constraints. The total stellar mass formed thus increases from ${450}_{-200}^{+1300}{M}_{\odot }$ ($\sim 0.1_{-0.1}^{+0.3}\%$ of the total stellar mass) to ${6600}_{-2200}^{+3300}{M}_{\odot }$ ($\sim 1.4_{-0.4}^{+0.7}\%$ of the total stellar mass). We again compared the BICs of the two models, finding that the former model has a BIC of 26 and the latter has a BIC of 23. The difference between the BICs decreases; however, the model that is constrained by both SED and B/R is still moderately preferred over the one that is only constrained by the SED.

Fig. 9. Same analysis in Fig. 7, but for region 2 of Warhol (highlighted in Fig. 1), as it is characterized by Palencia et al. (2025). The best fit model that is constrained by both SED and B/R ratio has a metallicity of 0.12 Z⊙ and a dust extinction of AV = 0.53, agreeing with the SED fitting result in Palencia et al. (2025) to within 1σ.

Fig. 10. SFH model corresponding to the analysis in Fig. 9, but for region 2 of warhol, as it is characterized by Palencia et al. (2025). Again, the inclusion of B/R ratio as a model constraint on top of the SED predicts a different SFH. The disagreement between the two solutions of SFH when the B/R ratio is incorporated as a constraint, however, is more drastic in the case confined to region 2 than in the case for the whole Warhol arc. This demonstrates the spatial variation of the SFH in the lensed galaxy, and the constraining power of B/R at different regions.

Region 2 is a subset of Warhol – the fact that the inferred SFH of region 2 forms more stellar mass in the last ∼50 Myr than that inferred for the whole Warhol arc indicates that there are severe degeneracies in the SED fitting (e.g., long-known stellar mass and metallicity degeneracy, Worthey 1994), especially when the differences between the different regions of galaxies are smoothed out during the computation of the SFH of the galaxy as a whole, as was previously found in literature (e.g., Sorba & Sawicki 2018; Giménez-Arteaga et al. 2024). This stresses the importance of having spatially varying information and highlights the possible role of B/R in improving the constraints on the SFH for lensed galaxies that are spatially resolved.

Also, our inferred SFH of region 2 is different from the one inferred in Palencia et al. (2025), even when the exact same SED is used as the constraint. The main driver behind this discrepancy is the model complexity, as Palencia et al. (2025) has only explored an exponential decay model that does not freely allow for recent star formation. Although all the SFHs have the highest SFR at ∼250 Myr – 1 Gyr, our solution prefers some elevated SFR during the last ∼50 Myr when a nonparametric SFH model is constrained by the SED, where the exact time of the elevated star formation can be better constrained when the B/R ratio is used as a constraint. As Palencia et al. (2025) did not compute the detection rate in HST F200LP, we cannot further comment on the difference between the different underlying SFH models adopted and how they predict the B/R ratio when only the SED is used as the constraint.

5. Discussion

5.1. Merits of B/R as constraints on SFH

Further constraints on the SFH can be obtained from spectroscopy targeting emission lines. Williams et al. (2025) conducted SED fitting on Warhol, which also included the emission line fluxes (Hα, Hγ, Hδ, [O III], [O II], and [Ne III]) of Warhol as constraints. Similarly to us, they used a nonparametric SFH with almost the same set of age bins (except that they have two bins of 1 − 2 Myr and 2 − 5 Myr, instead of a single bin of 1 − 5 Myr as we do; and they have an extra bin for ≳1 Gyr). Based on the inferred SFH, they predict a B/R ratio of $0.{07}_{-0.04}^{+0.07}$ when adopting an IMF slope of 2.3 (same as Kroupa IMF at M_⊙ ≥ 1.4). This is within ∼1σ of the observed value of 0.33 ± 0.26. Since the SFH shown in Williams et al. (2025) is not corrected for lensing magnification, we can only compare the relative SFR at different age bins between our SFHs by normalizing the SFR to the first star formation bin. We find their inferred relative SFR to be consistent with ours in each age bin to within 2σ, except for the last age bin of ∼250 Myr−1 Gyr, which is not associated with either emission line strength or transient detection (and thus the B/R ratio).

Although in this case the most recent relative SFR inferred by either emission line strength or the B/R ratio on top of the SED is similar, it is nevertheless worth noticing that the SFH constrained by emission lines, such as the Balmer lines, could be biased toward the most recent SFR within the last ∼10 Myr (Flores Velázquez et al. 2021), which is contributed primarily by star-forming regions that emit strongly in UV. On the other hand, the B/R ratio can trace star formation up to ∼50 Myr, until the detection rate per stellar mass begins to decrease significantly (cf. Fig. 4).

One of the remaining tests for the merits of the B/R ratio as a constraint of the SFH is whether our solution can reproduce the measured emission line fluxes, or, ultimately, the extent to which we can constrain the most recent SFH when the SED, emission line strength, and B/R ratio are all used as constraints together. The different sensitivities of all these three constraints could allow one to better constrain the recent SFR, especially regarding the burstiness of star formation. Limited to our customized code that constrains the SFH with the B/R ratio alongside the SED, we cannot predict the emission line strength and see how well our best-fitting solution reproduces the emission line strength reduced in Williams et al. (2025) under the current framework. A full exploration is deferred to future work.

5.2. Differential dust extinction

One could argue that differential dust extinction in galaxies (e.g., Mawatari et al. 2025) could also result in a change in B/R across regions as is shown in Fig. 5, without altering the SFH. A_V is always a free parameter in our MCMC analysis with a flat prior, whereas the two solutions (constrained by SED only, and both SED and B/R ratio) have a slightly different A_V of 0.85 ± 0.07 and 0.73 ± 0.12, respectively. Although the degree of freedom is well given, the fact that both solutions still prefer a similar amount of dust extinction means that such a dusty solution is required to explain the SED. Thus, the solution varies the SFH in order to resolve the B/R ratio without sacrificing the quality of fits towards the SED.

Also, despite the fact that one cannot rule out extremely high dust extinction, the current solution of A_V ≈ 0.8 is already very dusty. And, even with such a level of dust extinction, the solution that is only constrained by the SED also overpredicts the B/R ratio. It is not impossible that the detected transient stars are more heavily dust-obscured (for example, from their own stellar wind, e.g., Tambovtseva & Grinin 2008) compared with the other stellar population, and therefore way redder than what is expected from our stellar population synthesis. However, such a scenario can only be confirmed upon a spectroscopic analysis of these lensed stars in future detections; hence, the variation in SFH might be the simplest, most viable solution under the principle of Occam’s razor.

5.3. Uncertainties

Throughout the whole discussion, we have been assuming that the magnification map we adopted is a truthful representation of the real lensing magnification. The two primary factors that could have affected the lensing magnification are the macro lens model itself (e.g., Perera 2025), as well as perturbations in the magnification distribution arising from, for example, substructures in dark matter haloes (e.g., Moore et al. 1999; Dai et al. 2018).

The Warhol arc contains many unresolved features, presumably star clusters that are multiply lensed. These features form anchor points that allow one to pinpoint where the macro model’s critical curve should cut through (Broadhurst et al. 2025), with little room for uncertainty. The lens model adopted by Palencia et al. (2025) is the Diego et al. (2024a) lens model, in which one can see that the critical curve closely resembles the expectation demonstrated in the last figure in Broadhurst et al. (2025). In such a case, the region where the macro-magnification can be wrongly estimated is rather small. Even though the critical curve is shifted by a bit, the high magnification region (μ_m ≳ 10²) that contributes the most to the detection rate (see Palencia et al. 2025; Li et al. 2025) remains inside the Warhol arc, as is predicted by the Diego et al. (2024a) lens model. Far regions at the two sides of the critical curve of the Warhol arc are sufficiently far away from the critical curve and have low magnification (∼1 − 10)², such that the predicted magnification therein would not change a lot even if there were a small shift in the position of the critical curve. These low macro-magnification regions barely contribute to the detection rate (Li et al. 2025; Palencia et al. 2025), and hence the B/R inferred should not change a lot with the choice of lens model.

For substructures, Palencia et al. (2025) studied the role of nonluminous cold dark matter subhaloes (or millilenses) with masses of ∼10^6 − 8 M_⊙ in the detection rate of transients in the Warhol arc. The main finding is that although the existence of substructures redistributes the magnification, they primarily affect the spatial distribution of transient events (see also Williams et al. 2024; Broadhurst et al. 2025), and only affect the detection rate of transients down to a ∼1% level. Propagating this subtle uncertainty barely affects the resulting B/R ratio, which is currently dominated by Poisson noise. In such a case, our inference of SFH based on the B/R ratio should be invariant under the effect of substructures. Nevertheless, there are more possible cases of substructures to be considered (e.g., density modulations from wave-like dark matter, Amruth et al. 2023), and a full exploration is required in future works.

6. Conclusion

In this paper, we have found that the ratio of transient detection in two characteristic filters, HST F200LP and JWST F200W, which best capture the light of blue and red supergiants at z ≈ 1, is a very sensitive probe of the age of the underlying stellar population in lensed galaxies. Based on this fact, we have demonstrated that the B/R ratio, the transient detection rate in the two aforementioned filters, poses challenges for the SFH inferred solely by SED fitting, regardless of the model’s complexity. We hence demonstrate for the first time that including the B/R ratio itself as a constraint on top of the SED resolves such challenges and offers an alternative solution for the SFH that can explain both the B/R and SED for relevant galaxies. The solution constrained by the B/R ratio as well also has a lower BIC than the solution constrained by the SED only. With the higher spatial resolution of the B/R ratio measurement and the higher sensitivity to slightly older SFH up to ∼50 Myr ago, the B/R ratio could be a powerful method of studying the spatially resolved stellar population in relevant galaxies.

Our solution of Warhol’s SFH after incorporating the B/R ratio as a constraint on top of the SED increases the estimated stellar mass formed in the last ∼50 Myr by a factor of ∼60 for the whole Warhol arc, and a factor of ten for region 2, a subset of Warhol, with 3σ confidence. Although this change is subtle as it nevertheless only contributes to a small fraction of the stellar mass of the entire Warhol, this kind of study could drastically change our understanding of galaxies in the early Universe, when we are witnessing them grow in stellar mass (e.g., Weaver et al. 2023; Weibel et al. 2024). Also, the current signal-to-noise in the B/R ratio is low and has inferred a low significance in the change in inferred stellar mass. More transient observations will further test the significance of the B/R ratio as a constraint on the SFH.

The current record holder of lensed star candidates has a redshift up to ∼6 (Welch et al. 2022). With sufficient observations, and thus transient detections of galaxies across cosmic history, the B/R ratio can be a powerful tool in better assessing the star formation of galaxies in the Universe as early as ≲1 Gyr, and therefore understanding of the evolution of galaxies over cosmic time. Another interesting line of investigation is whether this change in stellar mass formed over the last ∼50 Myr in relevant galaxies is a systematic bias revealed by the B/R ratio. Highly resolved lensed galaxies with many transients detected by both HST and JWST, such as the “Dragon” (z = 0.73 with 40+ transients Kelly et al. 2022; Fudamoto et al. 2025), would serve as an excellent testing ground.

Acknowledgments

We sincerely thank the referee for the comments, which have eventually led to a more comprehensive analysis. We also thank Yoshinobu Fudamoto, Juno Li, and Arsen Levitskiy for the insightful discussion. S.K.L., J.L., and J.N. acknowledge support from the Research Grants Council (RGC) of Hong Kong through the General Research Fund (GRF) 17302023. J.M.P. received 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 support from project PID2022-138896NB-C51 (MCIU/AEI/MINECO/FEDER, UE) Ministerio de Ciencia, Investigación y Universidades. J.N. also acknowledges the support of the Dissertation Year Fellowship issued by the University of Hong Kong. P.L.K. and L.L.R.W. acknowledge the support by NASA/HST grants GO-15936 and GO-16278 from STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. AZ acknowledges support by grant 2020750 from the United States-Israel Binational Science Foundation (BSF) and grant 2109066 from the United States National Science Foundation (NSF), and by the Israel Science Foundation Grant No. 864/23. We made use of the following software: Python, NumPy, Matplotlib, SciPy, Astropy, SPISEA, emcee, and corner.

References

Appendix A: Selected simulated cases

Fig. A.1. Same as Fig. 4, but for the cases with different metallicities and dust extinction as denoted in the subtitles of each of the panels. We do not show the error bar for clarity issues.

Appendix B: MCMC analysis

We use emcee to carry out our MCMC analysis throughout this paper. The model likelihood, L, is given as:

$$\begin{array}{c}\hfill L\propto \text{exp}[-\frac{1}{2}{\displaystyle \sum _i^N}{(\frac{x_i(\theta )-{\mu }_i}{{\sigma }_i})}^2],\end{array}$$(B.1)

where x_i(θ) is the model prediction given a set of parameters to be optimized upon, θ, whereas μ_i is the observed value with uncertainty σ_i over N data points, i.

Fig. B.1. Markov chain Monte Carlo analysis regarding the constraining power of B/R (blue) and SED (red) of Warhol on the parameter space, given that we assume there is only one single burst of star formation at age t with metallicity Z★ and dust extinction AV. Since we only care about the inference of these three parameters, all the SED templates generated by SPISEA are pre-normalized to best fit the observed SED of Warhol. The contour levels correspond to 1 − 4σ level confidences. Histograms/contours in solid lines represent the case for a Kroupa IMF, where those in dotted lines (with a slight color tweak to improve visibility) represent the case for a top-heavy IMF. There is no significant difference in preference for star formation parameters constrained by either B/R or SED when a different IMF is adopted.

Fig. B.2. Same as Fig. B.1, but with a nonparametric SFH instead of two single star bursts. Since we concluded that the IMF is not important for our analysis, here we only show the case where a Kroupa IMF is adopted. The gray contour represents the case where both B/R and SED are used as constraints. Mass 1-5 corresponds to the mass formed in the 5 age bins of 1 − 5 Myr, 5 − 10 Myr, 10 − 50 Myr, 50 − 250 Myr, and 250 Myr−1 Gyr, respectively. The MCMC here corresponds to the results shown in Fig. 7 and Fig. 8.

Fig. B.3. Same as Fig. B.2, but confined the analysis to region 2 of Warhol. The MCMC here corresponds to the results shown in Fig. 9 and Fig. 10.

All Figures

	Fig. 1. RGB image of the Warhol arc (z = 0.94). We show the transients detected in this arc as coral circles (reported by Yan et al. 2023; Williams et al. 2025) and green triangles (reported by Kelly et al. 2022), respectively. We highlight areas characterized as region 2 by Palencia et al. (2025) in red.
In the text

Fig. 2. Subsets of MIST age isochrone (Choi et al. 2016) used in our SPISEA (Hosek et al. 2020) stellar population synthesis featuring a range of ages as shown in the legend. The metallicity of each of the subsets is reflected by their line style, as is shown in the bottom legend. None of them have any dust extinction, and the color represents the age of the isochrone, as is shown in the top legend. For each age isochrone, we also annotate some characteristic masses to indicate the mass-age dependency.

In the text

Fig. 3. Filter response curve of HST F200LP (blue shade) and JWST F200W (coral shade), and the blackbody spectrum of two stars at 30 000 K (blue curve) and 4500 K (red curve) – characteristic temperature of BSGs and RSGs, respectively – shown at the rest frame and observed frame (z = 0.94) wavelengths. One can see that the two filters best capture the brightest part of the blackbody of the two classes of stars.

In the text

Fig. 4. Transient detection rate per unit mass for stellar populations at different ages with a metallicity of 0.5 Z⊙ and AV = 0, calculated with Eq. (1) under the magnification map of Warhol from Palencia et al. (2025). Blue curves show the detection rate in the HST F200LP filter (detection limit of 30 mAB); red curves show the detection rate in the JWST F200W filter (detection limit of 29.7 mAB). Solid and dotted lines show the case of a Kroupa IMF and a top-heavy IMF, respectively.

In the text

Fig. 5. Expected B/R ratio calculated based on the F200LP and F200W transient detection rate in Warhol arc shown in Fig. A.1, with a Kroupa IMF. Each of the panels shows the case for one of the dust extinctions (upper row) or one of the metallicities (lower row), as is indicated by the subtitles. The colour of each curve indicates the metallicities (upper row) or dust extinction (lower row). The dashed gray line denotes B/R = 1 for reference. The fluctuations in ages earlier than ∼50 Myr are significant and arise from intrinsic variation due to stellar evolution, whereas those later than ∼50 Myr arise from sampling uncertainty given the low transient detection rate at these ages. We do not show the error bars for visualization purposes.

In the text

Fig. 6. SED of Warhol arc (black curve, from Palencia et al. 2025). The best fit (red curve) corresponds to a single star burst model at an age of ∼350 Myr with Z★ = 0.32 ± 0.12 Z⊙ and AV = 0.91 ± 0.04. This model predicts B/R = 10−7, orders of magnitude from the observed B/R = 0.33 ± 0.26. Among all the models that can reproduce the observed B/R ratio well, we show the best fit SED as a blue curve. This model corresponds to a single star burst at an age of ∼80 Myr with Z★ = 0.8 ± 0.4 Z⊙, and AV = 0.5 ± 0.3. We adopted a Kroupa IMF in both models.

In the text

Fig. 7. Same as Figure 6, but we use a nonparametric SFH model instead of a single star burst to fit the data. Again, the SFH that is constrained by the SED (red curve) does not predict the B/R ratio well (predicted $1.0_{-0.1}^{+0.2}$, versus observed 0.33 ± 0.26), while that is constrained by the B/R poorly predicts the SED (blue curves). Notice that the B/R ratio is not sensitive to the underlying stellar mass; thus, the SED inferred from B/R deviates significantly from the observed SED (dotted blue curve). We enforced a normalized version of the same SED (solid blue curve), whereby the relative SFR in each age bin of the underlying SFH is invariant, to acquire the best fit to the observed SED. Nevertheless, imposing both the SED and B/R as constraints infers a solution that fits both of the constraints to a sensible degree, which signifies the importance of including the B/R ratio as an extra constraint to better capture the true SFH. The best-fit solution has a metallicity of 0.44 ± 0.21 Z⊙ and a dust extinction of 0.91 ± 0.07, respectively.

In the text

Fig. 8. Corresponding best-fit SFH models of the Warhol arc that generate the SED and B/R ratio shown in Fig. 7 with the band representing the 1σ tolerance from MCMC analysis. One can see that the inclusion of the B/R ratio as a constraint on top of the SED gives rise to a different SFH, with more star formation over the past 10 − 50 Myr (factor of ∼60 with ≳3σ significance) and less star formation over the last 50 − 250 Myr (factor of ∼60 with 2σ).

In the text

	Fig. 9. Same analysis in Fig. 7, but for region 2 of Warhol (highlighted in Fig. 1), as it is characterized by Palencia et al. (2025). The best fit model that is constrained by both SED and B/R ratio has a metallicity of 0.12 Z⊙ and a dust extinction of AV = 0.53, agreeing with the SED fitting result in Palencia et al. (2025) to within 1σ.
In the text

Fig. 10. SFH model corresponding to the analysis in Fig. 9, but for region 2 of warhol, as it is characterized by Palencia et al. (2025). Again, the inclusion of B/R ratio as a model constraint on top of the SED predicts a different SFH. The disagreement between the two solutions of SFH when the B/R ratio is incorporated as a constraint, however, is more drastic in the case confined to region 2 than in the case for the whole Warhol arc. This demonstrates the spatial variation of the SFH in the lensed galaxy, and the constraining power of B/R at different regions.

In the text

	Fig. A.1. Same as Fig. 4, but for the cases with different metallicities and dust extinction as denoted in the subtitles of each of the panels. We do not show the error bar for clarity issues.
In the text

Fig. B.1. Markov chain Monte Carlo analysis regarding the constraining power of B/R (blue) and SED (red) of Warhol on the parameter space, given that we assume there is only one single burst of star formation at age t with metallicity Z★ and dust extinction AV. Since we only care about the inference of these three parameters, all the SED templates generated by SPISEA are pre-normalized to best fit the observed SED of Warhol. The contour levels correspond to 1 − 4σ level confidences. Histograms/contours in solid lines represent the case for a Kroupa IMF, where those in dotted lines (with a slight color tweak to improve visibility) represent the case for a top-heavy IMF. There is no significant difference in preference for star formation parameters constrained by either B/R or SED when a different IMF is adopted.

In the text

Fig. B.2. Same as Fig. B.1, but with a nonparametric SFH instead of two single star bursts. Since we concluded that the IMF is not important for our analysis, here we only show the case where a Kroupa IMF is adopted. The gray contour represents the case where both B/R and SED are used as constraints. Mass 1-5 corresponds to the mass formed in the 5 age bins of 1 − 5 Myr, 5 − 10 Myr, 10 − 50 Myr, 50 − 250 Myr, and 250 Myr−1 Gyr, respectively. The MCMC here corresponds to the results shown in Fig. 7 and Fig. 8.

In the text

	Fig. B.3. Same as Fig. B.2, but confined the analysis to region 2 of Warhol. The MCMC here corresponds to the results shown in Fig. 9 and Fig. 10.
In the text