Under Einstein’s Microscope: Measuring Properties of Individual Rotating Massive Stars from Extragalactic Microcaustic Crossings

,


INTRODUCTION
When a background galaxy overlaps the lensing caustic cast by a foreground galaxy cluster lens, a portion of that galaxy in the caustic vicinity appears immensely stretched.Stars residing in that region can be magnified by a few hundred to a thousand fold, and the most luminous ones are individually detectable (Miralda-Escudé 1991).A significant advance in our understanding of highly magnified stars is the realization that intracluster stars in the foreground cluster inevitably disrupt the smooth cluster caustic into a network of microcaustics on a fine scale ∼ 10 3 au (Venumadhav et al. 2017;Diego et al. 2018;Oguri et al. 2018).Highly magnified stars therefore appear variable as they traverse the nonuniform magnification pattern on the source plane as a result of the microcaustic network.
Highly magnified stars can be detected based on the variability induced by microlensing.Without microlens-ing, they are expected to be common at z = 1-2 in deep imaging that reaches 28-29 mag in the optical and near-IR (Diego 2019).Microlensing can further boost the flux by 1-3 mag temporarily, which greatly enhances detectability.Several blue supergiant stars have already been discovered by the Hubble Space Telescope (HST) in caustic-crossing giant arcs at z = 0.7-1.5 (Kelly et al. 2018;Chen et al. 2019;Kaurov et al. 2019).Searching for variability using the imaging difference technique, the Flashlight program recently detected a dozen highly magnified stars from this redshift range when the imaging limiting magnitude with HST in the optical improves from 28 to 30 (Kelly et al. 2022).James Webb Space Telescope (JWST) NIRCam imaging with the PEARLS program uncovered that redder magnified stars are commonly detected at infrared wavelengths at a magnitude limit 28.5-30 (Yan et al. 2023).Several intriguing magnified star candidates have been reported in more distant lensed galaxies at z ≃ 2-6 (Welch et al. 2022;Diego et al. 2023a,b;Meena et al. 2023;Welch et al. 2022), including quite bright objects (Vanzella et al. 2020;Diego et al. 2022).The nature of these sources is under active investigation.
Microlensing-induced brightening culminates every time when the source star transits a microcaustic.When the photosphere overlaps the microcaustic, a peak magnification µ pk ∼ 10 3 -10 4 is reached over the stellardiameter crossing time (hours to a few days).Flux evolution around such dramatic times has yet to be measured at high cadence but is observationally feasible owing to the expected immense magnification factor around the peak time.
In this work, we study the potential of caustic crossing as an exquisite gravitational microscope to measure stellar properties.During the caustic crossing, different parts of the stellar surface have significantly different magnification factors, and thus the time-dependent fluxes and colors encode information about the size of the photosphere, as well as its shape and surface temperature profile.
This can be compared to the study of finite source size effects in Galactic binary microlensing events with caustic crossings (Mao & Paczynski 1991), with some observed examples (Lennon et al. 1996;Yoo et al. 2004;Udalski et al. 2005).In extragalatic caustic crossings it is the most massive stars with the highest luminosities that are observationally selected.Therefore, caustic crossings with highly magnified stars provide a unique opportunity to study the properties of individual massive stars in distant star-formation environments that might differ fundamentally from local universe samples (in terms of, e.g., metallicity and interstellar medium density).
In this work, we will study single source stars that have axisymmetric surface properties owing to stellar rotation.We will introduce a parameterized model that accounts for the axisymmetric surface spectral energy distribution (SED) profile due to gravity darkening.By resorting to analytic arguments and by developing a numerical program to fit mock multifilter light curves for caustic-crossing events, we will show that, with all unknown extrinsic geometric and kinematic parameters marginalized over, the ratio of the surface rotation velocity to the breakup velocity (see Eq. ( 3)), denoted as ω, can be measured without a degeneracy with the stellar inclination as in spectroscopic rotation measurements.For massive stars, the evolution of stellar rotation connects to wind mass loss (Maeder & Meynet 2000a), surface element abundance enhancement (e.g.N and He; Maeder & Meynet (2000b)), and the explosive fate of the star (Woosley & Heger 2006;Yoon et al. 2006).Knowing ω for a sample of massive stars at z = 0.7-1.5, combined with the knowledge of the property of each star's host environment, can potentially improve our understanding of massive star evolution in the younger Universe.
Additionally, the equatorial radius R e and the bolometric luminosity L can be constrained to precision levels limited primarily by the unknown microcaustic strength (see Eq. ( 18)) the logarithmic uncertainties (in units of dex) for R e and L are 1/3 and 2/3 of that for the microcaustic strength, respectively.Fortunately, an informative prior for the latter can be derived from microlensing simulations knowing their number density and mass distribution (Venumadhav et al. 2017).
This paper will be organized as follows.In Sec. 2, we discuss an analytical model that quantifies the axisymmetric surface brightness profile of a rotating massive star and introduce intrinsic and extrinsic parameters of the model.In Sec. 3, we summarize how individual microcaustic crossing events can be modeled by a local fold lens model and define the relevant parameters.In Sec. 4, we discuss parameter degeneracy expected in fitting mock multifilter light curves during the caustic-crossing process to our theoretical model.In Sec. 5, we present examples of mock parameter inference under realistic observational conditions, which we carry out using a numerical code we have developed that performs inverse ray shooting and rapidly calculates the light-curve fitting χ 2 as a function of model parameters.Sec.6 is for further discussions, and we will close with concluding remarks in Sec. 7.

SOURCE STAR MODEL
In this section, we discuss how we model the surface brightness profile of axisymmetric rotating stars.

Surface Property of Rotating Stars
We consider general stars with surface rotation.The surface rotation angular velocity may not equal that of the stellar interior.The surface of a rotating star has an axisymmetric shape and has a latitudinal effective temperature profile T eff (θ) (where θ is the polar angle and is equal to π/2 minus the latitude).The von Zeipel theorem (von Zeipel 1924) states that the bolometric radiation flux F emerging from the surface of the rotating star is proportional to the local effective surface gravity g eff , which is centrifugally reduced on the equator.As a result, the polar caps are hotter than the equatorial region and have a higher surface brightness.
We adopt the analytic model of Espinosa Lara & Rieutord (2011) for the axisymmetric distribution of surface radiation flux F (r, θ, ω) and surface gravity g eff (r, θ, ω): More sophisticated than the simple power-law model T eff ∝ g β eff with β ≈ 1/4 for massive stars with a radiative envelope, this model agrees well with 2D numerical rotating star models.In this model, the reduced radius r is defined as the dimensionless ratio of the radial position of the stellar surface to the equatorial radius R e , i.e. r ≡ r/R e .This quantifies the axisymmetric deformation of the surface due to the centrifugal force.The quantities L and M are the bolometric luminosity and the stellar mass, respectively, and σ is the Stefan-Boltzmann constant.Following Espinosa Lara & Rieutord (2011), we have introduced an auxiliary variable ϑ, which is related to the polar angle θ via Eq.( 4).
The dimensionless rotation parameter ω is a measure of how dynamically important the surface rotation is.It is defined as Here Ω is the angular velocity of the surface, and Ω k is the Keplerian angular velocity on the equatorial surface, which is the critical rotating rate for material on the equatorial surface to be bound in a circular orbit.
The model requires 0 ⩽ ω < 1.If the value is close to breakup ω ≳ 0.7-0.8, a decretion disk is likely to form on the equator to form a Be star (Ekström et al. 2011).
The following implicit equation gives an analytic relation between the polar angle θ and ϑ: (4) The variable ϑ reduces to θ in the absence of rotation.
The shape of the star r(θ) can be found from a constant effective potential on the surface, whose value can be determined by considering the equator: It is assumed here that the gravitational potential at the stellar surface is well approximated by GM/r (i.e. the Roche potential), as the mass distribution of a massive star is rather centralized (Espinosa Lara & Rieutord 2011).Recasting the equation into the dimensionless form, we derive the shape equation The sequential procedure to compute the rotating star model is the following.We find r(θ, ω) by solving Eq. ( 6) and then find ϑ(θ, ω) = ϑ(r(θ, ω), θ, ω), and then we determine the effective surface gravity g eff (θ, ω) = g eff (r(θ, ω), θ, ω) according to Eq. (2) as well as the surface radiation flux F (θ, ω) = F (g eff (θ, ω), ϑ(θ, ω), θ) according to Eq. (2).Knowing F (θ, ω), the surface temperature T eff (θ, ω) is found from F (θ, ω) = σ T 4 eff (θ, ω).Defining characteristic quantities T eff,0 ≡ L/4πσR 2 e 1/4 and g eff,0 ≡ GM /R 2 e , we can express the surface temperature and the effective surface gravity in dimensionless forms In Figure 1, we show these reduced quantities as a function of the polar angle θ for a few different values of ω.

SED of Rotating Stars
To calculate the observed SED of a lensed star, we must integrate the spectral surface brightness over the lensed image.For OB stars, the blackbody SED that depends on a single parameter T eff would be too crude an approximation.We adopt the more realistic TLUSTY stellar spectrum templates computed for nonrotating OB stars (Lanz & Hubeny 2003, 2007).Each model outputs the emergent SED, which is characterized by three parameters: effective temperature, surface gravity, and metallicity.
The TLUSTY grid we use covers a range of effective temperature from T eff = 15, 000 to 30, 000 K and a range of surface gravity from log g = 1.75 to 4.75.To model rotating stars, we use the approximation (Collins & Harrington 1966;Maeder & Peytremann 1970) that the SED of the local emergent radiation at any point on the surface of the rotating star is that of a nonrotating star that has an effective temperature T eff , a surface gravity g eff and a metallicity Z.By assuming that radiation transfer through the stellar atmosphere is a local process, we can model the varying SED across the rotating star surface by interpolating a grid of nonrotating TLUSTY models.Limb darkening is not modeled in this work, but we will comment on that in Section 6.2.Our approximation is imperfect, as certain spectral features depend on the radiation-hydrodynamic solution describing the line-driving wind of a massive star.The wind solution can differ between nonrotating and rotating stars.Since the relevant observable for our purpose is wide-filter photometry, wind features are not expected to be important.
For examples shown in the later part of the paper, we will choose a TLUSTY grid with a metallicity Z = 0.2 Z ⊙ and a microturbulence v turb = 2 km s −1 .We will suppress explicit dependence on these two parameters in equations.
Note that highly magnified stars may be cool stars for which the TLUSTY models are not applicable.They may be red supergiants with T eff = 3500-4500 K (Diego et al. 2023a).Even the first detected star has a lower effective temperature T eff = 11, 000-14, 000 K (Kelly et al. 2018).Nevertheless, our framework is sufficiently general that any other stellar atmosphere models appropriate for different stellar types can be incorporated if needed.For a proof of concept, we will consider in this work relatively hot stars for which the TLUSTY templates are valid.
Let Sλ (λ, T eff , g) be the surface brightness of the TLUSTY model at wavelength λ.The surface brightness for a source star at a redshift z s observed at wavelength λ is The observed flux, with the lensing magnification accounted for, is then given by an angular integration over the lensed image: The observed flux convolved with a given photometric filter is where f (λ) is the filter throughput.The AB magnitude (Oke & Gunn 1983) is then computed as where λ p is the pivot wavelength of the filter.

Viewing Angles
For a rotating star, the observed surface brightness profile depends on the orientation of the star with respect to the line of sight, which is random.We parameterize this using two of the three Euler angles: Θ is the angle between the line of sight and the star's rotation axis, and Φ is the azimuthal position of the rotation axis when projected onto the plane of the sky.On the source plane, the surface brightness profile is simply related by rotations in the plane of the sky for situations with the same Θ but different Φ values.Figure 2 shows an example of how the surface brightness of a rotating star varies with the observer's perspective.The angle Φ is relevant when the star is lensed because the multifilter photometric outcome of the time-varying differential magnification effect during a caustic-crossing event depends on the orientation of the star in the sky plane relative to the direction of the lensing caustic.
In Figure 3, we consider the same example star as in Figure 2. We vary the dimensionless surface rotation ω while keeping other intrinsic parameters fixed.Broadband colors, corresponding to three promising HST filters (F435W, F555W, and F814W) with which the star may be observed, are shown as a function of ω and along different viewing angles.In this example, there are only up to ∼ 0.1 mag subtle color differences between a nonrotating star and a rapidly rotating star, in the absence of differential magnification.As we will see later, the detailed shape of the flux light curve during the caustic crossing (as a result of time-dependent differential magnification of the photosphere), rather than the color variation, provides the key information to distinguish between different ω values.

Summary of Stellar Parameters
In summary, we characterize a rotating source star with four intrinsic parameters: equatorial radius R e , bolometric luminosity L, mass M , and the dimensionless rotation rate 0 ⩽ ω < 1; two external parameters describe the viewing angle: the inclination 0 ⩽ Θ ⩽ π and the azimuthal angle 0 ⩽ Φ < 2π.

MICROCAUSTIC CROSSING
The situation of the highest observational promise is when a source star approaches a portion of a microcaustic cast by intracluster microlenses.Around this time, the star appears the brightest owing to a pair of highly magnified microimages.This lensing situation can be well described by a universal fold model applied to the microcaustic vicinity (Blandford & Narayan 1986;Schneider et al. 1992).In this section, we discuss this fold model in the context of this work.
In geometrical optics and under the thin-lens approximation, the angular coordinates x on the image plane are mapped to the angular coordinates y on the source plane through the ray equation where α(x) = ∇ψ(x) is the deflecting angle and ψ(x) is the lensing potential.
The lensing Jacobian matrix takes the following form: Following the notation in Venumadhav et al. (2017), we can Taylor-expand the quantities, 16) In addition to the convergence parameter κ 0 , the fold is uniquely characterized by the gradient vector The lensing Jacobian matrix in the caustic vicinity is  The ray deflection α(x) = (α 1 , α 2 ) is given by The parameter κ 0 is set by the macroscopic lens surface mass density over roughly arcsecond scale around the magnified star and can be determined to a decent precision with galaxy cluster mass modeling (Venumadhav et al. 2017;Kelly et al. 2018;Dai et al. 2020).By contrast, the relevant vector d ⋆ at each microcaustic crossing is set by the microlenses and has a random value.In the region where the macro critical curve is fully disrupted by microlenses, d ⋆ has a characteristic magnitude (Venumadhav et al. 2017) where κ ⋆ is the contribution of intracluster microlenses to the coarse-grained convergence and θ ⋆ is the typical angular Einstein radius of microlenses (Venumadhav et al. 2017;Diego et al. 2018;Oguri et al. 2018).In practice, an informative prior can be supplied for d ⋆ based on the knowledge of κ ⋆ and θ ⋆ , i.e. the abundance and typical masses of microlenses.
The above fold model enables us to calculate the lensed flux using the inverse ray-shooting method given a source star and its source-plane location.
Figure 4 shows snapshots of a lensed rotating star on the source plane and on the image plane at four different epochs during the microcaustic crossing.As the star approaches the caustic, a pair of images become increasingly bright and stretched.The image pair merge when the source star touches the caustic.The merged image shrinks and eventually disappears as the source star hides on the other side of the caustic.
The motion of the star sets the timescale of the caustic-crossing light curve.To the extent that the curvature radius of the caustic is large compared to the stellar radius, the relevant velocity component is the one perpendicular to the caustic, corrected for time dilation in the expanding Universe (Miralda-Escudé 1991;Venumadhav et al. 2017), Here v s and v l are the proper velocities of the source and the lens (relative to the Earth), ŝ is a unit vector in the sky perpendicular to the local caustic, z s and z l are the source redshift and lens redshift, and D s and D l are the angular diameter distances to the source and to the lens, respectively.For a given highly magnified star, z l and z s are often both known, and so are D s and D l .For each microcaustic crossing event, the value of v t is not known.We therefore include it as a free parameter in light-curve fitting.Finally, the exact time of caustic crossing t cr is unknown and has to be another free parameter.

PARAMETER DEGENERACY
The surface temperature T eff is mainly constrained by the multiband colors.The maximal observed fluxes constrain the luminosity L multiplied by the peak magnification µ pk , i.e. µ pk L ≃ µ pk 4π R 2 e T 4 eff .The value of µ pk is unknown; it is approximately given by µ pk ∼ (d ⋆ R e /D s ) −1/2 .Therefore, the combination R 3/2 e /d 1/2 ⋆ is well measured.It follows that the logarithmic uncertainty in R e is 1/3 of the uncertainty in d ⋆ and the logarithmic uncertainty in L is 2/3 of the uncertainty in d ⋆ .The combination R e /v t can be inferred from the timescale of peak magnification, and hence the logarithmic uncertainty in v t is 1/3 that in d ⋆ .These conclusions are confirmed by our mock parameter inferences (see Figure 6 and Figure 7).
Combining Eq. ( 1) and Eq. ( 2) and taking into account the dimensionless equations Eq. ( 4) and Eq. ( 6), we see that if two stars have identical L, R e , ω and Z but different M values, they have the same surface shape but only exhibit subtle differences in their SEDs that reflect the effect of g eff .This provides the unique information to infer M .Since this information cannot be efficiently extracted through wide-filter photometry, we expect the inference precision for M to be poor (see Figure 7 and Figure 6) and perhaps prone to uncertainty in the modeling of spectral line features in the SED.
In this study, we treat the macro convergence κ 0 as a known parameter.In reality, its inference from galaxy cluster lens modeling is subject to uncertainty.Since κ 0 directly impacts only the magnification perpendicular to the direction of image elongation, we expect that if uncertainty in κ 0 is included it will add to the degeneracy involving L, R e , d ⋆ and v t , for which uncertainty is in any case dominated by the significant uncertainty in d ⋆ .We therefore expect that the uncertainty in κ 0 does not impact the inference of ω.
We note that for the most massive stars, especially the evolved supergiants, the bolometric luminosity L approaches the Eddington limit L Edd = 4π G M c/κ T , where κ T is the opacity of Thomson scattering.The Eddington ratio L/L Edd is not precisely unity, and the precise value depends on stellar structure and wind models.Combining constraints from direct observation of caustic crossing and massive star models may enable a significantly improved mass measurement.
While the inference precisions for multiple intrinsic parameters are fundamentally limited by d ⋆ , which is not directly measurable, it is possible to derive a reliable and informative prior on d ⋆ .For a given highly magnified star on a given caustic-crossing giant arc, the parameters of the local macro caustic can be obtained through lens modeling, while both the abundance and mass dis-tribution of the microlenses can be obtained from photometric analysis of the intracluster light (e.g., Kelly et al. (2018); Kaurov et al. (2019); Dai et al. (2020)).With such knowledge, an accurate prior for d ⋆ can be derived, for instance, by numerically simulating random microlensing realizations.We crudely estimate that the fractional RMS in d ⋆ will be ≃ 0.35 dex, a number we use to set our priors for d ⋆ in the next section.For this value, we expect ∼ 0.1 dex uncertainty in R e and v t and ∼ 0.2 dex uncertainty in L, provided that errors of photometry do not dominate the error budget.

MOCK PARAMETER INFERENCE
In this section, we present mock parameter inference examples.We consider two scenarios: one star with a dynamically significant rotation speed ω = 0.4, and another star with a slow rotation ω = 0.1.In both cases, we set the same values for the other intrinsic and extrinsic parameters (true values and assumed prior distributions summarized in Table 1), corresponding to an M = 60 M ⊙ blue supergiant with a surface temperature T eff ≈ 20, 000 K and a bolometric luminosity L = 10 6 L ⊙ at a metallicity Z = 0.2 Z ⊙ and from z s = 1.
The mock light curves are the sum of injected light curves and random photometric noise.To collect crucial color information, we consider monitoring in three HST wide filters, F435W, F555W, and F814W, with the added random photometric noise corresponding to the 1σ AB magnitude limit 29.0, 29.0, and 29.3 per epoch, respectively.The photometric sensitivities assumed here are realistic for the HST, as they are comparable to the magnitude limits of the Hubble Frontier Fields (HFF) program (Lotz et al. 2017) when the exposure time per observation is rescaled to be about 1 hr per filter (compatible with the cadence ≲ 4 hr for three filters).JWST will be more sensitive but can only observe λ > 0.7 µm.The logarithmic likelihood function is taken to be minus one-half of the overall fitting χ 2 , under the assumption that individual photometric measurements at different epochs or in different filters have uncorrelated measurement errors.
In Figure 5 we show the example mock light curves, color variability, and the best model fit.The peak magnitudes in the optical are 25.2-25.5.These are comparable to the brightest magnitudes observed in the light curve of the first discovered highly magnified star in Kelly et al. (2018).Given that that star is from z s = 1.5, the observability represented by our mock event is not exaggerated.
In Figure 5, we also show model light curves (dashed lines) corresponding to a different parameter set that yields a fitting χ 2 that is approximately 55 units worse than the best-fit value.From both panels of Figure 5, we conclude that it is the detailed shapes of the light curves during caustic crossing, more than the color changes, that carry the dominant information for constraining the surface rotation rate.
In the Appendix A, we show posterior distributions that result from the mock parameter inferences we perform.Figure 6 shows the case of the fast-rotating star with ω = 0.4.The true value of ω is recovered with a fractional error of about 35%.The fractional uncertainties of R e , L and v t scale with the uncertainty in d ⋆1 , the component of d ⋆ along the degenerate direction, which is the limiting factor and is set by the prior we adopt for d ⋆ .This is reflected in the degeneracy between d ⋆1 and any of R e , L and v t as explained in Section 4. No significant degeneracy is seen between d ⋆2 , the other component of d ⋆ , and intrinsic stellar parameters.We are able to distinguish between the nearly pole-on or moderately inclined perspective and a perspective in the equatorial plane, and we derive a decent constraint on the position angle Φ.
For a comparison, Figure 7 shows the posterior distributions derived from observing the identical source star but with a dynamically unimportant rotation ω = 0.1.A similar absolute uncertainty is inferred for ω.The constraints on the viewing perspective are significantly weaker, as expected for a star with much weaker variation of temperature across its surface.The fractional uncertainties in R e , L and v t are similar to the ω = 0.4 case as is set by the prior used for d ⋆ .
It is worth mentioning the possible issue of nonaxisymmetric features on the surfaces of rotating massive stars.For instance, early-type B stars have been shown to display rotational variations with amplitudes of up to 0.01 mag (Balona 2022).Given our assumed photometric precision 0.03 mag under realistic observing conditions, this effect is not likely to have a significant impact on the parameter inference results.

Surface Rotation
An important finding from our study is the ability to infer dimensionless surface rotation.The most massive stars, relevant for extragalactic caustic crossings, are expected to have significant rotation when they start on the main sequence (Maeder & Meynet 2012).For an M = 50 M ⊙ zero-age main-sequence (ZAMS) star, the breakup velocity can be v eq = 800-900 km s −1 at Z = 0.1-1 Z ⊙ and higher at lower Z (Ekström et al. 2008).Based on observed examples, v i = 200-500 km s −1 are considered plausible values for the ZAMS surface ro-tation velocity.Dynamically significant initial values ω = O(1) are hence expected for massive stars.
However, the most massive stars on the main sequence (O stars) are currently not the most promising target highly magnified stars observationally (from z s = 0.7-1.5).Because of their high T eff ≳ 25, 000 K, they do not appear the brightest in HST or JWST bands (Kaurov et al. 2019).Indeed, the first detected highly magnified star is a cool supergiant with T eff = 11, 000-14, 000 K (Kelly et al. 2018).The magnified star reported by Chen et al. (2019) appears to be photometrically fit by a cool supergiant with T eff = 13600 K, while a dust-reddened, hot O supergiant with T eff ≃ 40, 000 K also fits photometry decently but is deemed much less probable.
Except under metal-poor conditions, angular momentum loss from stellar wind is expected to significantly reduce surface rotation through the main-sequence lifetime.It raises the question whether warm or cool supergiants (B supergiants) as evolved massive stars may have dynamically significant surface rotation.In the Galaxy, some examples of rotating blue supergiants are known.Sher 25 in the H II region NGC 3603 has M ≃ 40 M ⊙ , T eff = 22, 000 K and ω ≃ 0.2 (Hendry et al. 2008).Blue supergiant SBW1 has T eff = 21, 000 K and ω ≃ 0.2 (Smith et al. 2017).Another blue supergiant HD 168625 has R = 105 R ⊙ and T eff = 14000 K (perhaps most similar to the highly magnified star reported in Kelly et al. (2018)) and has a surface rotation ω ≃ 0.27 (Mahy et al. 2016).The significant surface rotations of these example blue supergiants are also evidenced from their associated ring-shaped nebulae.For this reason, SBW1 is considered an analog of the SN 1987A progenitor (Smith et al. 2017).Recently, it has been found in the IACOB spectroscopic survey that Galactic blue supergiants mostly have v sin Θ < 70 km/s for T eff < 21, 000 K, except for only a handful of stars with v sin Θ = 100-150 km s −1 (de Burgos et al. 2023).
It might be that there are rotating blue supergiants with higher ω in lower-metallicity galaxies given that wind-driven spin-down from initial rotation is less efficient (Meynet & Maeder 2002).Additionally, critical surface rotation can be reached for a blue supergiant that has recently contracted from a red supergiant phase (Heger & Langer 1998;Meynet & Maeder 2000).It is also plausible that tides in tight massive star binaries can make fast-rotating supergiant stars.In any case, the example rotating stars we consider here with ω = 0.1 and 0.4 bracket the likely range of ω for rotating blue supergiants with T eff < 21, 000 K. , and the temporal unit t0 is the time it takes for the star to move a distance of its diameter 2Re at vt = 8.2 × 10 2 km s −1 , which is t0 = 39.25 hr in this case.The observation cadence in each filter is about 0.1 t0 ≈ 4 hr.Injected model parameters can be found in Table 1.The deviated light curves and the associated color variability are also shown (dashed lines), in which the parameter ω = 0.7.In this example, the surface rotation rate ω is constrained mainly by the detailed shapes of the flux light curves (top panel) than by the subtle color variations (bottom panel).

Extensions of Model
Our model can be extended to account for other realworld effects.In practice, dust reddening has to be accounted for to correctly interpret the observed colors (Kelly et al. 2018).The dust reddening effect is not included in the model we have presented here, but it would be straightforward to introduce one more free parameter to the model provided that the reddening curve is known.Such a new parameter would not be degenerate with ω because flux variations due to differential magnification are time-dependent signals, while dust reddening is time independent.
There can be other microimages than the two that are highly magnified around the time of caustic crossing (Venumadhav et al. 2017).Those additional images are significantly fainter and have nearly constant fluxes over the short diameter crossing time t 0 = 2 R e /v t .Therefore, one additive flux constant per filter can be introduced as additional free parameters, and we do not expect them to be degenerate with the intrinsic parameters.
We have not accounted for limb darkening.This effect has been studied for nonrotating OB stars (Wade & Rucinski 1985;Claret 2004;Howarth 2011;Reeve & Howarth 2016) and is expected for rotating stars too.While it is not the intention of this work to incorporate the limb-darkening effect in addition to gravity darkening, our modeling framework can be generalized to account for it.In that case, the emergent SED at any given point on the stellar surface not only is a function of Z and the local (T eff , g eff ) but also depends on µ, the cosine of the oblique angle (the angle between the line of sight and the local surface normal).The dependence on µ can be derived from radiation transfer calculations for stellar atmosphere models.Reeve & Howarth (2016) pointed out that for OB stars limb darkening is moderately sensitive to surface gravity, which may provide new information for improving mass inference.In terms of observable effects, limb darkening is partially degenerate with gravity darkening for | cos Θ| ≈ 1 but not so for | cos Θ| ≈ 0. However, if limb darkening can be predicted by radiation transfer calculations as a function of (ω, L, M, Z), then it merely complicates the SED profile model but is not a free parameter.
Another obvious direction to extend this work would be to construct a model for binary systems.Massive stars are commonly found to be in binary systems, and the fraction of interacting binaries can be higher than 70% (Sana et al. 2012).If the binary separation is large enough (say, larger than several au), the surface brightness profiles of the two stars may each be well described by single-star models with appropriate parameters.In this case, the orbital period is expected to be much longer than the diameter crossing time, and the light curves would be trivial superpositions of two caustic crossings, with a time lag and a doubling of parameters.The more complicated situations could involve tidally induced surface deformations or ongoing mass transfer.The caustic-crossing phenomenon in principle would enable us to "see" the details.
If multiple microcaustic crossings of the same source star are observed, it would be possible to tighten up parameter inference by jointly analyzing light curves of all crossings.All crossings share the same model parameters except for v t and d ⋆ , and a reasonable prior distribution for d ⋆ specific to each highly magnified star can be theoretically derived.An in-depth study of this interesting method may be pursued in future work.

CONCLUSION
We have studied a method to measure stellar parameters of individual highly magnified stars in causticoverlapping lensed galaxies at cosmological redshifts.The method exploits light variability as a result of the differential magnification effect when the source stars transit microlensing caustics cast by intracluster stars.
We have constructed an analytic SED model that quantifies the axisymmetric surface brightness profile induced by the gravity-darkening effect for rotating O/B stars (Section 2).The SED model specifically builds on the TLUSTY SED templates for nonrotating O/B stars, although our framework is applicable to general stellar SED templates.Gravitational lensing in the vicinity of the microcaustic is modeled, at least near the epoch of caustic crossing, by a simple fold (Section 3).Combining the model SED profile and the lensing effect, we have introduced a light-curve model parameterized by intrinsic stellar parameters and extrinsic geometric and kinematic parameters.Based on the parameterized light-curve model, we have developed a general code for performing parameter inference with multifilter light curves.
We have demonstrated the method by performing mock parameter inference under observational conditions that are realistic for the HST (three wide filters, 1σ magnitude limit per epoch ∼ 29, and at a cadence of a few hours over a few days), for a blue supergiant source star at z s = 1 (Section 5).Our results suggest that the surface rotation velocity relative to the breakup value can be measured to an error σ ω = 0.1-0.2,without a geometric degeneracy, while other intrinsic parameters such as equatorial radius R e and bolometric luminosity L are limited by uncertainty about the microcaustic strength and are likely to be determined to 0.1 and 0.2 dex (1σ), respectively.Mass inference is less con-  1. Model parameters and the prior distributions we use in mock parameter inference.The physically allowed range for the dimensionless rotation parameter ω is 0 ⩽ ω < 1.We consider a prior range between 0 and 0.8 because numerical difficulty arises for ω values near breakup (ω ≈ 1) and the mock stars we choose are not rotating that fast in any case.Due to reflection symmetries it is sufficient to sample one quadrant of the position angle range, 0 ⩽ Φ ⩽ π/2.For d⋆, the same lognormal prior applies to both Cartesian components d⋆1 and d⋆2.The true values of the model parameters we choose are not fine-tuned.
straining but could be improved with a careful treatment of limb darkening.While there is much room for further improvement in the modeling of the surface brightness profile, this first study suggests a new opportunity to survey the individual properties of the most massive and luminous stars in distant high-z galaxies, which may advance our understanding of massive star formation and evolution.While we have not examined particularly lensed stars in higherz lensed galaxies (say, z > 2), that case could be even more interesting, since there would be increased chance of seeing stars in metal-poor systems (Welch et al. 2022) and cosmological redshifting should make hotter (and more common) massive stars observationally more accessible.Observational prospects for this might arise in the future.
In the situation of recurring microlensing brightening events in the same lensed galaxy, the feasibility to significantly constrain intrinsic and extrinsic parameters can help us determine which caustic-crossing events are associated with the same source star.This knowledge would be useful in realizing highly magnified stars as a novel dark matter probe, such as uncovering star-free subgalactic cold dark matter halos through astrometry (Dai et al. 2018;Abe et al. 2023;Williams et al. 2023) or dense minihalos composed of QCD axion particles as the dark matter (Dai & Miralda-Escudé 2020;Xiao et al. 2021).1D and 2D marginalized posterior distributions for the case of rapid surface rotation ω = 0.4.The blue and yellow lines mark the true parameters and the maximal likelihood parameters, respectively.In the panels showing 1D distributions, the prior distribution is shown by the shaded gray histogram.Among the three vertical dashed lines, the middle one marks the median of the posterior distribution, and the left and right ones mark the 16th and 84th percentiles of the distribution, respectively.In the panels showing 2D distributions, the three contours, from the innermost to the outermost, enclose the 68.27th, 95.45th, and 99.73th percentiles of the distribution, respectively.We refer readers to Table 1 for the units of the dimensionful parameters.

Figure 1 .
Figure 1.Reduced surface temperature Teff (θ, ω) and reduced effective surface gravity geff (θ, ω) as a function of the polar angle θ for three different values of ω.

Figure 2 .
Figure 2. Variation of the effective temperature T eff across the surface of a rapidly rotating star with ω = 0.8 on the source plane.From top to bottom we show three perspectives: along the rotation axis (Θ = 0, Φ = 0), in the equatorial plane (Θ = π/2, Φ = 0), and at an inclined angle (Θ = π/4, Φ = π/4).The star has an equatorial radius Re = 83.25 R⊙, a bolometric luminosity L = 10 6 L⊙, and a stellar mass M = 60 M⊙, corresponding to an evolved blue supergiant at Z = 0.2 Z⊙.We use an angular diameter distance Ds = 1.65 Gpc corresponding to a source redshift zs = 1.

Figure 4 .
Figure4.A series of snapshots of a caustic crossing event for the same fast rotating star as in Figure2, as seen on the image plane (left panels) or on the source plane (right panels), with the red line being the critical curve or the caustic.From top to bottom are four epochs during the process.As the star transits the caustic, the lensed image pair get closer, merge on the critical curve and finally disappear.For the caustic, we set κ0 = 0.83, d⋆ = (5.59,5.59) × 10 3 arcsec −1 .Note that the panels on the left (showing the image plane) cover a vastly different angular scale along the x1-direction compared to along the x2-direction, reflecting that the lensed images are highly elongated along the x1-direction.

Figure 5 .
Figure 5. Mock multiband light curves (top panel) and the associated subtle color variability (bottom panel) for a lensed rotating star with ω = 0.4 crossing a microcaustic.The center of the star passes the caustic at t = tcr (vertical yellow dashed lines), and the temporal unit t0 is the time it takes for the star to move a distance of its diameter 2Re at vt = 8.2 × 10 2 km s −1 , which is t0 = 39.25 hr in this case.The observation cadence in each filter is about 0.1 t0 ≈ 4 hr.Injected model parameters can be found in Table1.The deviated light curves and the associated color variability are also shown (dashed lines), in which the parameter ω = 0.7.In this example, the surface rotation rate ω is constrained mainly by the detailed shapes of the flux light curves (top panel) than by the subtle color variations (bottom panel).

Figure 6 .
Figure6.1D and 2D marginalized posterior distributions for the case of rapid surface rotation ω = 0.4.The blue and yellow lines mark the true parameters and the maximal likelihood parameters, respectively.In the panels showing 1D distributions, the prior distribution is shown by the shaded gray histogram.Among the three vertical dashed lines, the middle one marks the median of the posterior distribution, and the left and right ones mark the 16th and 84th percentiles of the distribution, respectively.In the panels showing 2D distributions, the three contours, from the innermost to the outermost, enclose the 68.27th, 95.45th, and 99.73th percentiles of the distribution, respectively.We refer readers to Table1for the units of the dimensionful parameters.

Figure 7 .
Figure 7. Same as Figure 6, but for the case of slower surface rotation ω = 0.1.