Constraints on Sub-Neptune Planet Candidate KOI-972.01 via Joint Variability/Gravity-darkening Analysis

We analyze Kepler photometry of transiting planet candidate KOI-972.01, accounting for both stellar variability and gravity darkening. KOI-972.01 stands out because of its small radius, less than that of Neptune, and because of its intermediate orbit period at 13.12 days, long enough to avoid significant tidal evolution, and thus it represents an underexplored exoplanet class. The parent star of KOI-972.01 is a rapidly rotating δ-Scuti variable, complicating transit lightcurve interpretation but also offering a potential independent source of stellar parameters. We measure the stellar rotation period (16.2 hr) by identifying the stellar rotation frequency and subsequently place a constraint on the stellar obliquity of no greater than 10, but have difficulty isolating individual oscillation modes in the periodogram owing to time variation of the δ-Scuti oscillations. After subtracting the stellar oscillations, lightcurve fits place the transiting object radius at 3.07 ± 0.09 R ⊕, but the shallow transit prevents useful constraints on the system’s spin–orbit alignment.


Introduction
Observations of exoplanet systems over the past 20 years have revealed a wide range of orbital configurations. The origins of these configurations, including misaligned planets (Winn et al. 2005), Hot Jupiters (Mayor & Queloz 1995), and highly eccentric orbiters (Naef et al. 2001), are not well understood. Many possible explanations for misaligned planets -those whose orbital planes differ significantly from their star's equatorial plane-exist, including disk torques (Batygin 2012), the Kozai mechanisms (Wu & Murray 2003;Libert & Tsiganis 2009), and internal gravity waves (Rogers et al. 2012), but the challenge of evaluating the validity of candidate formation pathways lies with the difficulty of measuring primordial orbital configurations.
Due to the observation bias toward high mass (or high radius), close-in planets, our library of discovered planets may not reflect all planetary systems. Close-in, massive planets experience strong tidal influences from their host stars (Jackson et al. 2008), increasing the likelihood that the orbital configurations that we observe are the result of an evolved system state and not indicative of the system's primordial configuration. Therefore, in order to get a more accurate picture of the planet population, we wish to extend our measurements to planets with orbital configurations outside of this narrow range of orbit periods and masses. Although this limitation is mainly due to the low transit probabilities for smaller, more distant planets, this expansion can also be achieved with the application of new techniques to overcome the inherent signalto-noise ratio problem for planets that are observed in this orbital regime.
One such technique is stellar variability analysis: the method of analyzing the frequency oscillations of variable stars to place constraints on the stellar rotation, and thus the stellar obliquity (in this work we define the stellar obliquity to be the angle between the stellar rotation axis and the plane of the sky). Goupil et al. (2000) showed that the separation between the modes of the dominant frequencies of oscillation can be used to calculate the star's obliquity and rotation rate. This method has since been used to classify star-planet interactions in WASP-33 (Herrero et al. 2011), clean the lightcurve of WASP-33 (von Essen et al. 2014), analyze HD 144277 (Zwintz et al. 2014), and measure the spin-orbit misalignment of KOI-976 (Ahlers et al. 2019).
Another technique is that of gravity-darkening fitting: fitting a value for the stellar obliquity from a transit with an asymmetric shape (Barnes 2009). These asymmetric transits are the product of stars with a pole-to-equator temperature gradient that is due to rapid rotation (von Zeipel 1924). Several previous analyses used the particular lightcurve shapes resulting from transits across gravity-darkened stellar disks to derive transit parameters, including the relative angle between the planet's orbit angular momentum and the stellar spin angular momentum-spin-orbit misalignment (Barnes et al. 2011(Barnes et al. , 2013Zhou & Huang 2013;Ahlers et al. 2014Ahlers et al. , 2015Masuda 2015;Ahlers et al. 2019Ahlers et al. , 2020aAhlers et al. , 2020bZhou et al. 2019). Of particular note, Barnes et al. (2015) constrained Kepler 1115b (KOI-2138) as spin-orbit aligned, one of only a handful of measurements of the spin-orbit alignment of sub-Neptune-radius planets (Bourrier & Hébrard 2014;Dalal et al. 2019;Gaidos et al. 2020;Mann et al. 2020;Hodžić et al. 2021).
Stars that exhibit both of these properties are therefore highly amenable to a joint analysis, where variability analysis is used to both remove the variability from the lightcurve and provide an independent measurement of the stellar obliquity, and a gravity-darkening fitting technique is used to provide an additional measurement of the stellar obliquity. This joint analysis is beneficial in that it allows for two independent measurements of the orbital geometries of a system, using only one data set. Such an approach is useful in that it (1) allows for a tighter constraint based on two different measurements and (2) provides a "back-up" measurement in the case that one method does not yield results. The latter is especially important for placing constraints on the orbital geometry of planets that are neither large nor close in to their star, as data on such systems are often noisier than for large, close-in planets.
Following the example of Ahlers et al. (2019), we apply this joint analysis using both variability analysis and gravity darkening to the system KOI-972, which consists of a transiting planet candidate around a rapidly rotating (Catanzaro et al. 2010) δ-Scuti variable (Uytterhoeven et al. 2011;Tkachenko et al. 2013;Balona 2014;Barceló Forteza et al. 2018) star. The KOI-972 system is both amenable to joint variability/gravitydarkening analysis and has not previously had a successful measurement of its spin-orbit misalignment. (Johnson et al. 2016 was unable to detect the planet in Doppler tomographic observations aimed at measuring spin-orbit misalignment.) In Section 2 we describe our data. In Section 3 we describe the variability analysis component of our technique, and in Section 4 we describe the gravity-darkening fit. We briefly discuss results in Section 5.

Data and Observations
We set out to measure the characteristics of the system Kepler Object of Interest (KOI)-972 (KIC 11013201). KOI-972 has one planet candidate that has a false positive probability of 0.01 (Morton & Johnson 2011) and shows no evidence of any transit timing variations (Ford et al. 2011Steffen et al. 2012). The 13.12 day orbital period (Borucki et al. 2011) is inside of the 10 to 100 day window that we classify as an intermediate orbital period within which we expect the measured orbital configuration to be primordial in nature. The star itself is a seismically active, 1.684 ± 0.004 solar mass star with polar effective temperature 7800 ± 300 K (Niemczura et al. 2015). (Note that this value and its error bars are set to account for all three temperature models of Niemczura et al. 2015 and their entire range of outputs.) Furthermore, KOI-972 has a reported v i sin of 117 ± 2 km s −1 (Niemczura et al. 2015) and analysis of the short cadence data and a periodogram of the data (Section 3) leads us to confirm KOI-972 as a δ-Scuti variable (Figure 1, 2).
Due to this stellar variability, the transits are obscured from view in the short cadence PDC_SAP data (all lightcurves are obtained from the Mikulski Archive for Space Telescopes: https://archive.stsci.edu/hst/). However, the δ-Scuti oscillations are not visible in the long cadence PDC_SAP data set, as its 30 min integration time exceeds most oscillation periods. We therefore work with the short cadence data set and employ the Linear Algorithm for Significance Reduction (LASR) technique described in Ahlers et al. (2018) to identify and fit frequencies in the KOI-972 lightcurve. Short cadence data are available for three quarters, Q2, Q5, and Q8, with large data gaps in Q2 and Q8 at times when Kepler did not make short cadence observations. We initially work with all three quarters, but include only data from Q5 in our final analysis (Section 3). Our final analysis therefore includes nine transit events.
In addition to being variable, we expect KOI-972 to be gravity-darkened as a result of its rapid rotation and size. We calculate a theoretical gravity-darkening parameter for KOI-972 using Claret & Bloemen (2011; this calculation is described in Section 4). We also use theoretical values of the firstand second-order limb-darkening parameters as described in Sing (2010) to parameterize the brightness profile of the stellar disk. In Table 1 we list these and all other premeasured or precalculated values for KOI-972 that we use in our transit analysis.
To conduct our transit analysis, we use the χ 2 -minimization fitting routine described in Barnes et al. (2011) that accounts for gravity-darkening across the stellar disk. However, this technique only works for nonvariable lightcurves. Analysis of the system is therefore broken up into two phases: (1) using the LASR technique described in Ahlers et al. (2018), we remove stellar variability from the short cadence lightcurve and perform the variability analysis, and (2) using the cleaned lightcurve, we apply the gravity-darkened fitting technique described in Barnes et al. (2011).

Variability Removal
Based on the classification of KOI-972 as a δ-Scuti variable, we are able to use the Linear Algorithm for Significance Reduction (LASR) technique outlined in Ahlers et al. (2018) to subtract off sinusoidal stellar variability in preparation for an orbital parameter fit. LASR operates under the assumption that the total stellar variability of a δ-Scuti is the result of a linear we also show the residual, after the subtraction of the variability fit from the original data. A transit is clearly visible at the center of the residual that was previously obscured by the high level of photometric variability. The gray box marks the area that was removed from the data set during the LASR variability fit. Final fits are done with multiple folded and binned transit events. combination of individual frequencies of sinusoidal stellar oscillation. Taking one oscillation mode at a time, consisting of the frequency, amplitude, and phase of the oscillation, LASR employs a least-squares fitting technique to fit for the values of the oscillation mode. This is similar to traditional prewhitening, however, LASR operates with a goodness-of-fit parameter in frequency space, rather than in time space. Initial guesses for the frequencies are obtained by generating a Lomb-Scargle periodogram of the entire short cadence baseline and identifying frequencies that are statistically significant in frequency space ( Figure 2).
Before generating our periodogram, we remove long-term secular variations, such as those arising from instrument anomalies and spacecraft orientation, using a boxcar filter with width equal to three times the transit duration (as measured in the long cadence data). High error, outlier data points are also identified and removed by hand. We then use the long cadence data to identify the time of transit center and remove the transits from the short cadence data. This is done to ensure the oscillation removal does not affect the transit profile.
We create the Lomb-Scargle periodogram from the short cadence data with transits removed, and we identify initial guesses for the significant frequencies of oscillation. Working from the most to least statistically significant peak, we fit the exact parameters for each frequency of oscillation sequentially and remove it from the baseline. This process is stopped when the peaks being removed match the level of noise in the periodogram. When all peaks have been iteratively fit and removed in this manner, we sum each fit together to produce a total fit of the stellar variability. The list of all oscillations Figure 2. Lomb-Scargle periodogram generated from one 30 day chunk of Q5 short cadence data of KOI-972. Here, we plot Lomb-Scargle power (VanderPlas 2018) as a function of frequency, showing the false "triplet" with central peak located around 390 μHz, which in reality consists of three independent oscillation modes. Left: we plot on a linear scale to highlight the large power of the three central oscillation modes. The identified rotation peak, ν å , is indicated. Right: we plot on a logarithmic scale to highlight the large number of oscillations present in the photometry. The horizontal dashed line indicates the lower Lomb-Scargle power limit at which we stopped fitting peaks. The identified rotation peak is indicated. identified in this manner is given in the Appendix. Finally, we remove the variability signal from the original lightcurve with the transits, interpolating through the previously removed transit points, which produces both a cleaned lightcurve and the variability data we use for variability analysis of the system (Figure 1). Further analysis of the short cadence data shows that quarters Q2 and Q8 have a standard deviation of the data points higher than that of Q5, likely as a result of the large data gaps within the quarters. We therefore elect to work only with the data from Q5. Initial work fitting the frequencies of oscillation in this quarter additionally shows that the frequencies, amplitudes, and phases of the oscillations vary with time. These time dependencies are further discussed in Section 5. As a result of these variations, we break the quarter into three sections, each roughly 30 days in length. We choose break points at these intervals because they (1) correspond to already occurring breaks in the data, and (2) balance the need to reduce the influence on the oscillation modes of the time variations with the requirement that LASR have access to enough data points to obtain an accurate fit.
This time variation of the oscillation modes also makes it difficult to identify the primary rotational splitting triplet commonly used for variability analysis. Rotational splitting peaks are identified by the presence of a central peak of highest significance, with two lower significance peaks roughly equally spaced on either side of the central peak. Furthermore, the peaks making up the primary triplet tend to be the highest significance peaks in the periodogram. The three most significant peaks in our data, all of which are stable across the three sections of data, visually appear to make up such a triplet in the Lomb-Scargle periodogram, with the central peak located at roughly 390 μHz (Figure 2). Using the rotational splitting equations from Ahlers et al. (2019), a first-order calculation of the stellar rotation frequency suggested by these three peaks gives an estimated stellar rotational period between 3.16 and 3.95 hr. Both of these values indicate rotation faster than the 4.13 hr break-up speed calculated using the Gaia Collaboration et al. (2018) values. Incorporating second-order effects, we revise this estimate to approximately 4.67 hr, corresponding to a rotation frequency of approximately 67.6 μHz. This rotation period is still on the edge of the star's break-up speed and we find no frequency around this value. Therefore, our calculated rotation rates are inconsistent with both our measured rotation period and the v i sin reported by Niemczura et al. (2015). We therefore conclude that the configuration of the three primary frequencies is inconsistent with rotational splitting from an L = 1 mode.
Furthermore, we identify a significant frequency at 17.11 μ Hz that is constant across all three data sections and has the visual properties associated with stellar rotation frequencies ( Figure 2; Ahlers et al. 2019). A peak at this value corresponds to a rotation with a period of 16.2 ± 0.3 hr, consistent with the v i sin reported by Niemczura et al. (2015). We therefore conclude 17.11 μHz to be the true stellar rotation frequency, ν å , and take this peak as further evidence that the appearance of the three primary frequencies of oscillation as a rotational splitting triplet is a coincidence.
Although both this lack of a rotational splitting triplet and the time variation in the frequencies prohibits in-depth variability analysis, the identification of the stellar rotation frequency allows us to place a constraint on the stellar obliquity. Ahlers et al.   . Note that we cannot break the degeneracy between a positive or negative impact parameter. KOI-972 is also shown to be gravity-darkened. The darker horizontal band on the star is the stellar equator, with the north pole located at the top of the image, tilted slightly into the page. The darker band is indicative of the cooler surface temperature at the equator, which results from the gravity darkening. Top right: we show the final transit (dots) along with our best fit (line) plotted around time at transit center. The transit shown is the folded and binned transit used for the gravity-darkening fit. The small transit depth makes it difficult to properly characterize, but we would expect the transit to be slightly asymmetric due to the differences in the stellar surface temperatures the planet passes over during its transit. Bottom right: we show the residual after subtracting the transit fit from the transit itself. , and the v i sin listed in Table 1 to place a 2σ upper limit on the stellar obliquity of no greater than 10°.
We are also able to remove the variability signature from the original baseline. Once all frequencies of oscillation are identified in each of the three data chunks, the linear combination of all the frequencies is phase fit to the original short cadence data. These fit values are then subtracted from the original data, producing a lightcurve with constant out-oftransit flux and visible transits (Figure 1). The resulting δ-Scuti oscillation-subtracted lightcurve can then be used with our gravity-darkening fitting technique.

Transit Fitting
Gravity-darkening is a pole-to-equator temperature gradient that arises due to rapid stellar rotation. The rapid rotation results in a smaller effective surface gravity at the stellar equator and thus a lower effective surface temperature. The von Zeipel theorem (von Zeipel 1924) relates this difference in surface gravity to a difference in temperature by the equation where β is referred to as the gravity-darkening exponent. Using the model outputs of Claret & Bloemen (2011), we assign a theoretical value for the gravity-darkening exponent of KOI-972 to be 0.13, based on the stellar surface gravity, polar temperature, and metallicity. This value indicates that KOI-972 is at least weakly gravity-darkened. Furthermore, taking the star's q-value (defined as the ratio of the star's centrifugal force to its equatorial surface gravity) to be 0.06, we calculate an effective (equatorial) temperature of 7700 K. The difference between the effective and polar temperatures results in a large stellar flux disparity between the polar and equatorial regions, as flux is proportional to T 4 . In combination with the star's reported v i sin of 117 ± 2 km s −1 (Niemczura et al. 2015), these results indicate rapid rotation. We therefore predict KOI-972 to be gravity-darkened and apply the gravity-darkened fitting technique described in Barnes (2009) to the cleaned data we get from LASR's variability removal.
First, for each of the three chunks of data used in our LASR analysis, we subtract the variability from the original short cadence data chunk, which includes the transits. (We mask all transit events from the original oscillation fitting process.) The three resulting data chunks without variability are coadded together to create one time series. We then fold this time series across the length of one orbital period, producing one composite transit with data from nine individual transit events. The time series is then recentered around time t = 0 and binned to 120 seconds to reduce fit computation time, where 120 seconds is chosen to match the original cadence of the data. We then use the resulting transit as the base on which we perform our gravity-darkened fit (Figure 3). The known system parameters listed in Table 1-with the exception of surface gravity and metallicity, which are not fit parameters-are used as fixed inputs for our fits.
The fits themselves are done with the transitfitter program described in Barnes et al. (2011). We do polar fits in which transitfitter numerically integrates the stellar flux based on the input parameters and fits the resulting lightcurve to the data using a Leavenburg-Marquardt χ 2 minimization technique (Press et al. 2007). The program transitfitter can account both for oblate star shape and a transit asymmetry resulting from gravity-darkening. Ultimately, we wish to constrain the true spin-orbit misalignment of the system, which we define as the angle between the star's rotation axis and the planet's orbital angular momentum vector. This quantity cannot be fit for directly, but can be calculated with a combination of the orbital inclination, the stellar obliquity, and the longitude of the ascending node. These three parameters, and others, can be fit for directly with our lightcurve fitting method.
Our initial fits are unable to constrain the stellar obliquity or longitude of the ascending node due to the small transit depth resulting from the planet's small radius. The small transit depth makes it difficult to disambiguate any observed effects from noise, preventing a useful gravity-darkening based constraint of the stellar obliquity. We attempt to break this degeneracy by fixing the value of the longitude of the ascending node, thereby reducing the number of parameters to fit. The stellar radius is left fixed at the Gaia Collaboration et al. (2018) value for the same reason. In a series of fits, the longitude of the ascending node is fixed at 10°intervals across its range, and the resulting fit behavior of the stellar obliquity is tracked. In these fits we fit for planetary radius, orbital inclination, transit time at center, normal star flux, and stellar obliquity. These fits show a slight preference for low stellar obliquities but are ultimately not significant, as the χ 2 difference between the individual fits falls outside of the required 2σ confidence level.
The best fit in this set-defined as the fit with the lowest χ 2 value-corresponds to an ascending node of 150°. We therefore take this to be our most plausible value for the longitude of the ascending node, and thus leave it fixed in subsequent fits. Since the stellar radius plays an important role in determining the transit profile, we repeat our best fit allowing transitfitter to fit for the stellar radius. Therefore, we refit the time series, measuring the stellar radius, planet radius, orbital inclination, transit time at center, normal star flux, and stellar obliquity, while still keeping the stellar rotation period and v i sin constant. Using these results we then calculate the impact parameter and stellar oblateness. Best-fit results and corresponding χ 2 values are shown in Table 2. The final fit, overlayed on the binned transit, and the fit residual, are shown in Figure 3.

Results
Using LASR, we are able to fit and remove the variability signal from three short cadence data sections roughly 30 days in length. First-and second-order calculations of the stellar rotation speed using these variability fits predict a rotation speed inconsistent with the observed rotation peak at 17.11 μHz. Due to the agreement between this observed rotation peak and the v i sin reported by Niemczura et al. (2015), we conclude 17.11 μHz is the star's true rotation frequency. The discrepancy between the observed rotation frequency and the theoretical value calculated using the three most significant oscillation modes indicates that the three primary modes cannot actually result from rotational splitting and thus do not represent a rotational splitting triplet. Despite the lack of a clear rotational splitting triplet, we are able to place a constraint on the stellar obliquity using the measured rotation frequency, indicating the stellar obliquity is no greater than 10°.. Analysis of the variability fits across the three data chunks also shows that the oscillation modes are time dependent, with the frequencies, amplitudes, and phases of these oscillation modes changing across the three data chunks; however, our initial comparison across the three data chunks shows no easily discernible pattern in these changes. Ahlers et al. (2019) previously used the LASR technique on Kepler short cadence data using PDC_SAP flux and did not observe similar time variations in the oscillation modes. Furthermore, we do not observe time variations in the three most significant modes, nor in our identified stellar rotation frequency. We take this as evidence that these variations are therefore astrophysical in origin, and not instrumental or systematic. Further analysis of this time dependence is beyond the scope of this work.
The results of our gravity-darkened transit fit show that KOI-972.01 is a 0.274 ± 0.008 R J planet. We also corroborate the 1.55 R e stellar radius measurement made by the Gaia Collaboration et al. (2018); however, we are not able to place a constraint on the stellar obliquity, and thus on the orbital alignment of the system, using gravity-darkening lightcurve fitting due to the low transit depth resulting from the planet's small radius. All fit results are shown in Table 2.

Conclusion
Although we are not able to measure the orbital geometry of the system comprehensively, we are able to positively identify a stellar rotation period of 16.2 ± 0.3 hr and place an upper limit of 10°. on the stellar obliquity. This rotation speed is high enough to predict the presence of rotational splitting. Future work to better understand the time variation of the oscillation modes may allow for the identification of the rotational splitting triplet. Measurements of the system's orbital geometry could then be made by further constraining the stellar obliquity.
Based on the planet's size, and the star's unique variability signature, we identify KOI-972 as a prime target for future observations and analysis. KOI-972.01ʼs size places it in an interesting range, as there is no corresponding solar system planet in this radius range. Furthermore, only a handful of planets in this radius range have measured spin-orbit misalingments (Bourrier & Hébrard 2014;Barnes et al. 2015;Dalal et al. 2019;Gaidos et al. 2020;Mann et al. 2020;Hodžić et al. 2021).
Finally, we determine that the joint variability analysis and gravity-darkening fitting technique described in Ahlers et al. (2019) requires stellar variability that is not time dependent to be successful.   Note. Each oscillation is subtracted from its respective chunk to remove the variability from the lightcurve. The time variations of the oscillations can be seen when comparing oscillations between the chunks. The three chunks are presented from first to last as they appeared in the original time series. The first three oscillations in each chunk are the three oscillations that make up the false "triplet." The identified stellar rotation peak is bolded in each chunk. All oscillations are presented in order of decreasing significance, as determined by the periodogram.