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

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 $\sim 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 microlensing, 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\simeq 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 $\mu_{\rm pk}\sim 10^{3}$–$10^{4}$ is reached over the stellar-diameter 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 $\omega$, 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 $\omega$ 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 $\chi^{2}$ as a function of model parameters. Sec. 6 is for further discussions, and we will close with concluding remarks in Sec. 7.

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

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 $\vec{x}$ on the image plane are mapped to the angular coordinates $\vec{y}$ on the source plane through the ray equation

where $\vec{\alpha}(\vec{x})=\vec{\nabla}\psi(\vec{x})$ is the deflecting angle and $\psi(\vec{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,

In addition to the convergence parameter $\kappa_{0}$, the fold is uniquely characterized by the gradient vector

The lensing Jacobian matrix in the caustic vicinity is

The ray deflection $\vec{\alpha}(\vec{x})=(\alpha_{1},\alpha_{2})$ is given by

The parameter $\kappa_{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 $\vec{d}_{\star}$ 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, $\vec{d}_{\star}$ has a characteristic magnitude (Venumadhav et al., 2017)

where $\kappa_{\star}$ is the contribution of intracluster microlenses to the coarse-grained convergence and $\theta_{\star}$ 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 $\vec{d}_{\star}$ based on the knowledge of $\kappa_{\star}$ and $\theta_{\star}$, 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 $\vec{v}_{s}$ and $\vec{v}_{l}$ are the proper velocities of the source and the lens (relative to the Earth), $\hat{\vec{s}}$ 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_{\rm cr}$ is unknown and has to be another free parameter.

The surface temperature $T_{\rm eff}$ is mainly constrained by the multiband colors. The maximal observed fluxes constrain the luminosity $L$ multiplied by the peak magnification $\mu_{\rm pk}$, i.e. $\mu_{\rm pk}\,L\simeq\mu_{\rm pk}\,4\pi\,R^{2}_{e}\,T^{4}_{\rm eff}$. The value of $\mu_{\rm pk}$ is unknown; it is approximately given by $\mu_{\rm pk}\sim(d_{\star}\,R_{e}/D_{s})^{-1/2}$. Therefore, the combination $R^{3/2}_{e}/d^{1/2}_{\star}$ is well measured. It follows that the logarithmic uncertainty in $R_{e}$ is $1/3$ of the uncertainty in $d_{\star}$ and the logarithmic uncertainty in $L$ is $2/3$ of the uncertainty in $d_{\star}$. 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_{\star}$. 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}$, $\omega$ 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_{\rm 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 $\kappa_{0}$ as a known parameter. In reality, its inference from galaxy cluster lens modeling is subject to uncertainty. Since $\kappa_{0}$ directly impacts only the magnification perpendicular to the direction of image elongation, we expect that if uncertainty in $\kappa_{0}$ is included it will add to the degeneracy involving $L$, $R_{e}$, $\vec{d}_{\star}$ and $v_{t}$, for which uncertainty is in any case dominated by the significant uncertainty in $\vec{d}_{\star}$. We therefore expect that the uncertainty in $\kappa_{0}$ does not impact the inference of $\omega$.

We note that for the most massive stars, especially the evolved supergiants, the bolometric luminosity $L$ approaches the Eddington limit $L_{\rm Edd}=4\pi\,G\,M\,c/\kappa_{T}$, where $\kappa_{T}$ is the opacity of Thomson scattering. The Eddington ratio $L/L_{\rm 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_{\star}$, which is not directly measurable, it is possible to derive a reliable and informative prior on $d_{\star}$. 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 distribution 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_{\star}$ can be derived, for instance, by numerically simulating random microlensing realizations. We crudely estimate that the fractional RMS in $d_{\star}$ will be $\simeq 0.35\,$dex, a number we use to set our priors for $d_{\star}$ in the next section. For this value, we expect $\sim 0.1\,$dex uncertainty in $R_{e}$ and $v_{t}$ and $\sim 0.2\,$dex uncertainty in $L$, provided that errors of photometry do not dominate the error budget.

In this section, we present mock parameter inference examples. We consider two scenarios: one star with a dynamically significant rotation speed $\omega=0.4$, and another star with a slow rotation $\omega=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_{\odot}$ blue supergiant with a surface temperature $T_{\rm eff}\approx 20,000\,$K and a bolometric luminosity $L=10^{6}\,L_{\odot}$ at a metallicity $Z=0.2\,Z_{\odot}$ 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$\sigma$ 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 $\lesssim 4\,$hr for three filters). JWST will be more sensitive but can only observe $\lambda>0.7\,\mu\mathrm{m}$. The logarithmic likelihood function is taken to be minus one-half of the overall fitting $\chi^{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 $\chi^{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 $\omega=0.4$. The true value of $\omega$ 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_{\star 1}$, the component of $\vec{d}_{\star}$ along the degenerate direction, which is the limiting factor and is set by the prior we adopt for $\vec{d}_{\star}$. This is reflected in the degeneracy between $d_{\star 1}$ and any of $R_{e}$, $L$ and $v_{t}$ as explained in Section 4. No significant degeneracy is seen between $d_{\star 2}$, the other component of $\vec{d}_{\star}$, 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 $\Phi$.

For a comparison, Figure 7 shows the posterior distributions derived from observing the identical source star but with a dynamically unimportant rotation $\omega=0.1$. A similar absolute uncertainty is inferred for $\omega$. 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 $\omega=0.4$ case as is set by the prior used for $\vec{d}_{\star}$.

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.

We have studied a method to measure stellar parameters of individual highly magnified stars in caustic-overlapping 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$\sigma$ magnitude limit per epoch $\sim 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 $\sigma_{\omega}=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$\sigma$), respectively. Mass inference is less constraining 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 higher-$z$ 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).

The authors would like to express thanks to the Berkeley Physics International Education (BPIE) program through which X.H. visited UC Berkeley as an undergraduate student and during which this research project was initiated. The authors thank Maude Gull for useful discussions. L.D. acknowledges research grant support from the Alfred P. Sloan Foundation (award No. FG-2021-16495) and the support of the Frank and Karen Dabby STEM Fund in the Society of Hellman Fellows.