Measurement of the photosphere oblateness of $\gamma$ Cassiopeiae via Stellar Intensity Interferometry with the VERITAS Observatory

Spectroscopic measurements of the H$\alpha$ line of the subgiant $\gamma$ Cassiopeiae ($\gamma$ Cas) indicate that it is a Be-type star rotating at nearly its breakup velocity about an axis inclined by $60^{\circ}$ to our line of sight. Thus, it is expected to have a pronounced equatorial bulge. The history of Michelson interferometric measurements of $\gamma$ Cas spans many decades (see §3 of Stee, Ph. et al., 2012, for a brief overview). These measurements reveal an extended decretion disk that suggests a $122^{\circ}$ position angle of the rotational axis, measured east of celestial north (Sigut et al., 2020). However, these measurements, performed at infrared wavelengths, are unable to measure the photosphere itself.

Stellar intensity interferometry (SII) at optical wavelengths has the potential to measure the photosphere itself at sub-milliarcsecond scale. Originally developed in the 1950’s and used through the early 1970’s (Hanbury Brown & Twiss, 1954; Hanbury Brown, 1974), SII has very recently seen a renaissance with the availability of large light collectors in imaging atmospheric Cherenkov telescope (IACT) arrays, high-speed electronics and modern computing at several observatories (Abeysekara et al., 2020; Abe et al., 2024; Zmija et al., 2023; Acharyya et al., 2024). By correlating the intensity of light at separated telescopes, rather than the wave amplitude (as done with Michelson interferometry), SII is robust against effects of optical imperfections, including atmospheric turbulence, and works well at optical wavelengths, providing sub-milliarcsecond resolutions with 100 m baselines. Relative to Michelson interferometry, this comes at a cost of reduced signal-to-noise for the same observing duration.

In this paper, we present the first SII measurement sensitive to the shape, as well as the size, of a star. The work is organized as follows. Section 2 gives a brief description of the VERITAS stellar intensity interferometer (VSII) and the datasets discussed here. Section 3 presents the data analysis and extracted squared visibilities. Section 4 uses simple geometric models of the star to extract size, shape, and orientation information from the measurements. The dependence of the squared visibility on $u$ and $v$ is fitted with uniform disk and uniform ellipse models. Section 5 discusses a complete stellar model of a rapid rotator with limb and gravitational darkening. Comparison of our data with this model provides tight constraints on the size, rotational speed, and orientation of the star. In Section 6, we summarize and provide an outlook for future work.

VERITAS is composed of four 12-m diameter IACTs, usually used for very high-energy ($E>100\,\rm GeV$) gamma-ray astronomy (Holder et al., 2006). During bright-moon periods, the system is used for SII measurements. The four telescopes are individually referred to as T1, T2, T3, and T4, and are unequally spaced with separations ranging between about 80 and 170 meters, and interferometric projected baselines ranging between about 35 and 170 meters. We have presented further details of the VSII stellar intensity interferometer in previous publications (Abeysekara et al., 2020; Acharyya et al., 2024). The optical light from the observed star is detected by a single Hamamatsu super-bialkali R10560 photomultiplier tube (PMT) located at the telescope focal plane, and continuously digitized at a 250 MHz rate with 8 bit resolution. We use a narrow-band filter ($\lambda=416$) with a bandpass of 10 nm to focus on the photosphere of the star. An upgrade associated with the present study is the use of a CPU-based server that simultaneously correlates all six telescope pairs with a frequency-domain algorithm. This correlator produces correlations identical to the FPGA-based correlator using a time-based algorithm that was used in our previous measurements (Scott, 2023). For some observations in the present paper, a 10-MHz clock improved the start-time synchronization.

Table 1 lists the $\gamma$ Cas observations for which data were collected, over 6 bright moon periods between January 2023 and December 2024. Intensity interferometry is generally robust against ambient light for bright stars; VSII operations are reliable for target-to-moon angles from $30^{\circ}$ to $95^{\circ}$ (Kieda et al., 2021). Target-to-moon angles are listed in table 1.

Over 160 pair-hours of data were used, after quality cuts discussed below. The baseline coverage in the star image Fourier reciprocal ($u-v$) plane is shown in figure 1. The relative telescope separation lies in the $u-v$ plane, the celestial east-north plane, following general interferometric convention. (See § 4 of Thompson et al. (2001) and Dravins et al. (2013) for a further description of the $u-v$ plane.)

In this section, we fit the visibilities with simple models, to extract an effective size, shape, and orientation of the star. It is common to model a star as a uniformly illuminated disk of angular radius $\theta_{UD}$. In this case, the squared visibility for a projected baseline $b$ is

where $\lambda=416$ nm and $J_{1}$ is the first-order Bessel function. $C_{\rm UD}$ is the proportionality constant between actual squared visibilities and peak integrals measured from correlation functions, which we treat as a fit parameter. The red curve in figure 7 is the best fit of equation 4 to the data. The effective uniform disk angular radius is $0.44\pm 0.01{\rm(stat)}\pm 0.02{\rm(syst)}$ mas, consistent with the result reported by the MAGIC collaboration (Abe et al., 2024) of $0.515\pm 0.038{\rm(stat)}\pm 0.023{\rm(syst)}$. Our systematic uncertainties were calculated by varying the peak quality cuts and threshold and not removing tracking corrections or accounting for the background light contributions. As discussed below, this star is not expected to be circular, so some tension between the values of $\theta_{UD}$ measured by different observatories may be expected.

A rapid rotator like $\gamma$ Cas is expected have an equatorial bulge and anisotropic gravitational darkening, which would lead to an oblate source with spatial extension perpendicular to the axis of rotation. Therefore, fitting the data with a uniform circular disk may not be appropriate, as the length scale may depend on the orientation of the baseline. An intuitive way to explore this possibility is to measure the uniform disk effective angular size as a function of angle of the projected baseline in the $u-v$ plane, $\phi_{b}$. The angle with the smallest effective angular diameter would correspond to the position angle of the star rotation axis.

In figure 8, the visibility measurements are sorted according to ranges in $\phi_{b}$. In the center panel, the baseline for each measurement is shown in the $u-v$ plane, and $22.5^{\circ}$-wide slices are defined and shown as differently colored regions. Measurements that fall within each slice are shown in the surrounding panels and fitted with the uniform disk model of equation 4 to extract an effective one-dimensional angular diameter for a small range in $\phi_{b}$. While the visibility measurement from T1T4 with baseline of 145 m passes our threshold cut (see figures 7 and 10), we suspect it may be a fluctuation, like those discussed in section 3.6. Therefore, we report fit parameters obtained both including and excluding this point, and it is drawn as a red square in Figures 8, 10, and 11.

Figure 9 shows the effective angular diameter versus the average position angle of the slice. The data are well fit with the function

The average effective diameter is $\overline{\theta}_{UD}=0.48\pm 0.01$ mas. The minimum diameter, the minor axis, is at $\phi_{b}=\phi^{*}=117^{\circ}\pm 4^{\circ}$, which we associate with the position angle of the stellar rotation axis and is consistent with previous measurements of the accretion disk of this star (Sigut et al., 2020). The extracted magnitude of the modulation, $\delta=0.13\pm 0.02$ corresponds¹¹1If an oscillation $A\sin(2\phi)$ is uniformly integrated over bins of width $\Delta$, the resulting histogram will exhibit oscillations of magnitude $A\tfrac{\sin\Delta}{\Delta}$. In our case, this is small effect: the finite-binning-corrected oscillation magnitude is $0.13\tfrac{\pi/8}{\sin(\pi/8)}=0.133$. to a maximum-to-minimum ratio of $\tfrac{1+\delta}{1-\delta}=1.30\pm 0.05$.

While figures 8 and 9 provide an intuitive way to estimate the anisotropy and orientation of the rapid rotator, it is more accurate to perform a simultaneous fit of all visibilities as a function of $u$ and $v$. The hour angle, $H$, for a given source and observatory, uniquely defines the $u-v$ coordinates for each telescope pair, so we fit the data as a function of hour angle (Dravins et al., 2013). The data is fitted with a uniformly illuminated ellipse model (Tycner et al., 2006)

where

Here, $r$ is the ratio of the major to minor axes ($r\geq 1$), $\theta_{\rm min}$ is angular diameter of the minor axis, and $\phi^{*}$ is the position angle of the minor axis of the ellipse.²²2Note that Tycner et al. (2006) describes the model for $r\leq 1$. For a rapid rotator, then, $\phi^{*}$ is the position angle of the axis of rotation.

The proportionality between the actual squared visibility and the peak integrals is $C_{\rm norm}$, which we treat as a fit parameter. Figure 10 shows the peak integrals (divided by $C_{\rm norm}$) for each telescope pair, as a function of local hour angle.

This more complete model describes the data better than the uniform disk model discussed above; see table 2. The angular size of the minor axis is $\theta_{\rm min}=0.43\pm 0.02{\rm(stat)}\pm 0.02{\rm(syst)}$ mas, and the axis ratio is $r=1.28\pm 0.04\pm 0.02$. Fitting an elliptical source with a uniform circular model (equation 4) would result in a uniform disk radius ($\theta_{D}$) somewhere between 0.43 mas and 0.55 mas, depending on the orientation of the baseline. Therefore, the value of $\theta_{UD}$ extracted from any given measurement will depend on the orientation of the baselines used during observation. This complicates comparisons between, say VSII and MAGIC (Abe et al., 2024), when using a circular source model on a rapid rotator.

It also means that $\overline{\theta}_{UD}$ from equation 5 and $\theta_{UD}$ from equation 4 represent slightly different quantities. While figures 8 and 9 are an intuitive way to display the non-circular nature of the star, the characteristics of the elliptical appearance of the star are most accurately extracted by a fit to an elliptical image model.

The extracted position angle of the rapid rotator star axis is estimated to be $\phi^{*}=116^{\circ}\pm 5^{\circ}{\rm(stat)}\pm 7^{\circ}{\rm(syst)}$, consistent with measurements of the decretion disk at longer wavelengths (Sigut et al., 2020), indicating that the decretion disc is generally aligned along the equatorial bulge of the photosphere itself.

We have computed synthetic squared visibilities, high-resolution spectra and spectral energy distributions for the photosphere of $\gamma$ Cas using PHOENIX (Hauschildt & Baron, 1999, version 20.01.02B), a code for nova, supernova, stellar and planetary model atmospheres. We use a Roche–von Zeipel model for a rapidly rotating star (Aufdenberg et al., 2006; Sackrider & Aufdenberg, 2023), modeled as an infinitely concentrated central mass under uniform angular rotation.

Each axially symmetric model is parameterized by a trigonometric parallax, $\varpi$, equatorial angular diameter, $\theta_{\rm eq}$, polar effective temperature, $T_{\rm pole}$, polar gravity, $\log(g)_{\rm pole}$, fraction of the critical (breakup) angular rotation rate, $\Omega/\Omega_{\rm c}$, and the von Zeipel gravity darkening parameter, $\beta$. As seen from Earth, the model is further parameterized by the inclination $i$, and the position angle, $\phi^{*}$, of the visible stellar rotation axis.

The intensity of each point on the surface of the model star, as viewed from Earth, is interpolated from a grid of 106 PHOENIX stellar atmosphere radiation fields: specific intensities at 78 direction cosines between surface normal and the direction of an observer and 106 ($T_{\rm eff},\log(g)$) pairs (with step sizes of 500 K in $T_{\rm eff}$ and 0.1 in $\log(g)$) between the hot, high gravity pole and the cooler, lower gravity equator. A model rotating very close to the critical limit ($\Omega/\Omega_{\rm c}$ = 0.9999) has a surface temperature difference of $\Delta T_{\rm eff}\simeq$ 12000 K and a gravity difference of $\Delta\log(g)\simeq 1.2$.

For our model of $\gamma$ Cas we fixed $\varpi=5.94$ mas (van Leeuwen, 2007, $5.94\pm 0.12$ mas), $i=60^{\circ}$ as constrained by Lailey & Sigut (2024) (see also Sigut et al. (2020) and Tycner et al. (2006)), and $\beta=0.20$, as empirically determined for rapid rotating stars imaged with CHARA/MIRC (Che et al., 2011).

We computed 7750 models for $\gamma$ Cas with 31 $\theta_{\rm eq}$ values from 0.40 mas to 0.70 mas, 25 $\phi^{*}$ values from 90^∘ to 138^∘ and 10 $\Omega/\Omega_{\rm c}$ values from 0.9000 to 0.9999. With six baselines and 118 hour angles, this corresponds to a table of 5,487,000 synthetic visibility values.

Figure 11 shows the VERITAS data (the same as in Figure 10) with a Roche-von Zeipel stellar model. Figure 12 shows the distributions of the best-fit parameters. The $\chi^{2}/ndf$ fit to the data confronts the upper physical bound of unity, therefore we can only define a $1\sigma$ lower limit on the rotation rate at 97.7%

Figure 13 is a rendering of the calculated light intensity distribution corresponding to the Roche-von Zeipel model that fits the data best.

For the best-fitting median values, $\theta_{\rm eq}$ = 0.60 mas and $\Omega/\Omega_{c}$ = 0.99, $T_{\rm pole}$ = 26500 K provides a reasonable match to available absolute archival spectrophotometry in the UV optical range (see Figure 14). A $\log(g)_{\rm pole}$ = 3.82 yields a mass of 15.1 M_⊙, at the top of the range 13-15 M_⊙ adopted by Smith (2019). Table 3 shows the best fit and the derived parameters.

Our adopted polar gravity and inclination yield a projected rotation velocity of $v\sin i=389\pm 20~{\rm km\,s^{-1}}$. Chauville et al. (2001) found $v\sin i$ values 420-425 ${\rm km\,s^{-1}}$ by comparing models to the He I $\lambda$ 4471 line. Figure 15 shows a comparison between the best-fit model to the squared visibility data and a high-resolution spectrum in the VSII band around 416 nm. Although the H$\delta$ line is filled with emission from the disk, the model roughly matches the wings of the line profile. The other recognizable absorption line is He I $\lambda$ 4026, which is stronger than the model with roughly the same width.

We have used the VERITAS telescope array to perform intensity interferometry measurements on the rapid rotator $\gamma$ Cassiopeiae, combining data from more than 160 pair-hours of observation. While the size and elongation of the decretion disk of this star have been well measured at infrared wavelengths, measurements of the photosphere at visible wavelengths have previously only been sensitive to its size. The six telescope pairs of the VSII observatory provide good coverage of the Fourier $u-v$ plane. By varying the orientation as well as the magnitude of the baseline between telescopes, we were able to measure the size, the equatorial bulge, and the orientation of the axis of rotation of the photosphere. This is the first time that intensity interferometry has been used to measure the oblateness of a star.

A uniform ellipse fit of our data indicates that the angular size of the minor axis of the photosphere is $0.43\pm 0.02\pm 0.02$ mas, with a major-to-minor axis ratio of $1.28\pm 0.04\pm 0.02$. The position angle of the minor axis, which corresponds to the rotator axis, is $116^{\circ}\pm 5^{\circ}\pm 7^{\circ}$, in agreement with earlier measurements of the accretion disk.

We have also compared our data to a Roche-von Zeipel model of a rapid rotator, which includes effects of limb and gravitational darkening as well as temperature gradients. The temperature profile was adjusted to match published spectrophotometry measurements. Our data is well fit within this model with an equatorial diameter of $0.604^{+0.041}_{-0.034}$ mas, position angle of $114.7^{+6.4}_{-5.7}$ degrees, and a lower limit on the angular velocity very near the critical value, $\Omega/\Omega_{c}=0.977$. The equatorial radius is 10.9${}^{+0.8}_{-0.6}\,\rm R_{\odot}$ and the mass of the star in the model is $15\pm 2\ M_{\odot}$.

Previous stellar intensity interferometry measurements have been limited to extracting the overall size of the photosphere and resolving the extended decretion disk in its H$\alpha$ light (Matthews et al., 2023). In order to take advantage of our coverage in the $u-v$ plane and gain sensitivity to the stellar shape, we found it important to remove artifacts from external radio interference and manual tracking corrections. Our remaining systematic uncertainties are dominated by uncertainties in the clock start synchronization, the peak-finding algorithm, and by residual correlated noise in the background. The overall systematic error is on the order of the statistical error, suggesting that only adding hours at previously unobserved hour angles would increase the precision of our measurement.

Having established the capability of stellar intensity interferometry to measure not only the size, but also the shape and orientation of a photosphere, we are continuing analysis of other rapid rotators, to expand the scope of data in the stellar catalog of bright stars. This will provide valuable input and constraints to stellar models that must match both the spectral and the geometrical consequences of extreme rotation.

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, by NSERC in Canada, and by the Helmholtz Association in Germany. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and the collaborating institutions in the construction and operation of the VERITAS and VSII instruments. The authors gratefully acknowledge support under NSF grants #AST 1806262, #PHY 2117641, and #PHY 2111531 for the construction and operation of the VSII instrumentation. We acknowledge the support of the Ohio Supercomputer Center and of an Ohio State University College of Arts and Sciences Exploration Grant. We acknowledge the financial support of the French National Research Agency (project I2C, ANR-20-CE31-0003). This work made use of the High Performance Cluster “Vega” at Embry-Riddle Aeronautical University. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original description of the VizieR service was published in 2000, A&AS 143, 23. Part of this research is based on INES data from the IUE satellite.