Spin–Orbit Misalignment and Precession in the Kepler-13Ab Planetary System

For about one-third of measured exoplanets, the orbital axis is tilted relative to the stellar spin axis by more than two standard deviations (Albrecht et al. 2012), a situation referred to as spin–orbit misalignment. In nearly all of these cases the exoplanets are hot Jupiters, which challenges our understanding of how planets form, migrate, and evolve. Studying such misaligned systems is therefore imperative if we wish to understand the origin and dynamical history of systems that differ remarkably from our own.

The Rossiter–McLaughlin (RM) effect (Holt 1893; Schlesinger 1910) has traditionally been the primary method used to probe spin–orbit misalignment in exoplanetary systems. This tool relies on the distortion of stellar line profiles that occurs when a transiting body eclipses part of the rotating star, and it is commonly used in the study of eclipsing binary stars (e.g., McLaughlin 1924; Rossiter 1924). However, the RM approach has two significant drawbacks. First, it can only constrain the sky-projected spin–orbit misalignment λ. It cannot provide information about the stellar inclination i_* or orbital inclination i_orb, so the true spin–orbit misalignment in three-dimensional-space ψ cannot be determined without other means (e.g., Winn et al. 2007). Second, the RM effect is far less precise for hot, rapidly rotating stars. This is because hot stars generally have few spectral lines, and these lines will be significantly broadened by rapid stellar rotation (Johnson et al. 2014).

Yet rapid stellar rotation lends itself to spin–orbit determination through transit photometry alone, as shown by Barnes (2009). Fast rotation results in stellar oblateness, which induces a gradient in the surface gravity g of the star. Von Zeipel (1924) describes the effect this has on the effective temperature (and hence emitted flux) at each point on the stellar surface,

$$\begin{eqnarray}&&\displaystyle \frac{T}{{T}_{\mathrm{pole}}}={\left(\displaystyle \frac{g}{{g}_{\mathrm{pole}}}\right)}^{\beta }\end{eqnarray} \tag{ 1 }$$

where the exponent β characterizes the strength of this gravity-darkening phenomenon and is theoretically 0.25 for a star with a radiative envelope. If the planet's orbit is misaligned with respect to the star's axis of rotation, the planet's transit chord will pass over a nonuniform stellar surface, resulting in an asymmetric transit light curve. This gravity-darkening technique provides an advantage over the RM effect, as it can measure both λ and i_* simultaneously based on the full shape of the transit anomaly.

Kepler-13Ab was the first planetary system in which spin–orbit misalignment was identified using transit photometry alone (Barnes et al. 2011; Szabó et al. 2011). A number of groups have since investigated the transit asymmetries it presents (e.g., Szabó et al. 2014; Masuda 2015; Howarth & Morello 2017). Both Barnes et al. (2011) and Masuda (2015) implement gravity-darkened transit models (based on Barnes 2009) to characterize the spin–orbit misalignment, while Howarth & Morello (2017) employ a model described by Espinosa Lara & Rieutord (2011), approximating the rotational distortion of the host star's surface with a Roche equipotential. There are a number of important factors to consider when comparing these studies. As highlighted by Howarth & Morello (2017), the effective temperature of Kepler-13A has reported values between 7650 and 9107 K (Brown et al. 2011; Szabó et al. 2011; Huber et al. 2014; Shporer et al. 2014), the lower end of which places the star near the boundary between radiative and convective envelope regimes. This makes it unclear how strongly gravity-darkening affects stars like Kepler-13A; for convective envelopes β is theoretically just 0.08 (Lucy 1967). The intensity of gravity darkening is also influenced by stellar mass: low-mass stars show more exaggerated gravity darkening or higher β values (Barnes et al. 2011). However, the reported mass of Kepler-13A is also poorly constrained, with values ranging from 1.72 to 2.47 M_⊙ (Borucki et al. 2011; Szabó et al. 2011; Huber et al. 2014; Shporer et al. 2014). This provides little clarity as to which stellar envelope regime and β value most accurately describe the star. The best-fit parameter values from previous gravity-darkened transit models—and those presented herein—should therefore be interpreted with these factors in mind.

Kepler-13Ab has also been the subject of some disagreement between literature values of the sky-projected spin–orbit misalignment λ. Doppler tomography performed by Johnson et al. (2014) was inconsistent with the results of Barnes et al. (2011), but Masuda (2015) reconciled the two by adjusting the quadratic limb-darkening parameters and fixing λ to the later spectroscopic measurement. Howarth & Morello (2017) find a similar solution with their model.

Moreover, not only does the transit show asymmetries characteristic of a spin–orbit misaligned system, but it is also known to vary in duration over the length of the Kepler observations. This is indicative of spin–orbit precession caused by the rapidly rotating star's quadrupole moment (Szabó et al. 2012, 2014; Masuda 2015). The precession can very clearly be seen in the temporal variation of the impact parameter between transits.

To date, no evidence of this variation has been seen in the phase curve and secondary eclipse of the planet, though given the amplitude of the variability and the sublime quality of the Kepler data, this possibility is worth exploring. In a system not unlike Kepler-13Ab, Armstrong et al. (2016) recently presented the first detection of variations in the peak offset of a planet's phase curve. HAT-P-7b, a hot Jupiter, showed signs of atmospheric variability indicative of changing cloud coverage—something only seen previously in the atmospheres of brown dwarfs (e.g., Radigan et al. 2012; Wilson & Rajan 2014; Rajan et al. 2015).

In this paper, we perform an entirely independent analysis of the transit data of Kepler-13Ab by modeling the binned and phase-folded transit light curves for 12 short-cadence quarters, determining both the spin–orbit misalignment and precession rate. We also analyze the secondary eclipse in an effort to further confirm the change in impact parameter seen in the transit data.

We use the Kepler simple aperture photometry (SAP) data outlined in Table 1 for our analysis of the light curve. For the transit, we use the short-cadence (60 s exposure) data, while for the phase curve we use both long- and short-cadence data, binning the latter into 30 minute bins to match the cadence of the former. We do not include the long-cadence data in our transit analysis, as it tends to smear out important features and would consequently reduce the precision of our fitted transit model. We discard any bad data points flagged by Kepler in the SAP files and normalize the data to the out-of-transit median. We exclude quarter 13 from our analysis because the standard deviation of its raw out-of-transit data was much higher than other quarters.

2.1. Companion Flux Contribution

2.2. Removal of Phase Variations and Systematics

We then remove the following out-of-transit variations from the data in each SAP file.

• 1.  
  Variations due to Doppler boosting and ellipsoidal signals, F_m, both of which are effects due to the mass of the planet. Doppler boosting is a combined result of brightness increasing along the radial velocity vector and a red/blueshift of the stellar spectrum that shifts parts of the spectrum in and out of the bandpass, influencing the observed brightness. Ellipsoidal variations are due to the tidal bulge raised on the stellar surface as the planet orbits, changing the observed surface area of the star.
• 2.  
  The brightness contribution of the planet, F_p, which varies as the illuminated side of the planet moves in and out of view.
• 3.  
  The third harmonic signal in the planet's period, F₃, identified in Esteves et al. (2013) and (2015, from here on E15) but not yet well understood in terms of the underlying cause. They suggest that it could be due to the spin–orbit misalignment presenting as a perturbation of the ellipsoidal variations.

Notably, all of these variations depend on orbital phase ϕ, and we refer the reader to E15 for the detailed equations describing them. We remove the variations from the normalized out-of-transit flux such that

$$\begin{eqnarray}&&F^{\prime} (\phi )=\displaystyle \frac{F(\phi )-{F}_{p}(\phi )}{1+{F}_{m}(\phi )+{F}_{3}(\phi -{\theta }_{3})}\end{eqnarray} \tag{ 2 }$$

where ϕ runs from 0 to 1 with mid-transit at $\phi =0$, and ${\theta }_{3}$ is the cosine third harmonic phase offset. For ${\theta }_{3}$ and other constants used to calculate the phase variations, we use the values from Table 8, Model 1 in E15.

We then calculate the out-of-transit median and remove any outliers above 3σ. To remove instrumental signals, we use a least-squares fit of a third-order polynomial to ±0.2 days around the transit and divide the variation-subtracted SAP flux in this phase range by the resulting fit. We elected to use this method rather than fit and divide out the Kepler co-trending basis vectors (CBVs) in order to avoid complications introduced by the chosen number of CBVs fitted and their relative contributions to the systematic signal.

2.3. Transit Phase-folding and Binning

Before performing a quantitative analysis of the phase curve and transit data, we consider whether variations are evident from a purely qualitative perspective. Figure 1 shows the difference between each phase curve and the mean overall phase curve. Unlike the transit data, we use the reduced phase curve data from E15, which includes both short- and long-cadence data binned into 30-minute bins. The secondary eclipse occurs at $\phi =0.5$. It is difficult to discern any noticeable gradient in the phase curves by eye; however, the boundary of the secondary eclipse suggests the possibility for variability; the entirety of its ingress and egress fit into single bins, so these points are more strongly affected by discrepancies between phase curves. We therefore elect to investigate the secondary eclipse (though with a less rigorous model than that used for the transit data) and we outline this analysis in Section 5.4.

Zoom In Zoom Out Reset image size

Figure 2 shows the difference between each corrected short-cadence transit and the mean overall transit. Here, we find qualitative confirmation of the change in transit shape caused by decreasing impact parameter, as reported by Szabó et al. (2012) and Masuda (2015). The vertical gradient seen near ingress and egress indicates that the transit duration is increasing with time; an effect of decreasing impact parameter (and thus orbital inclination). The gradient at egress is noticeably less pronounced than its ingress counterpart. This could be due to the slight variation in the time of mid-transit T₀ we identify in Section 5.1. In the following section, we describe our transit model used to quantify the variation we see in the transit data.

Zoom In Zoom Out Reset image size

We use the modeling code simuTrans,⁶ a numerical integrator created to model transit light curves affected by gravity-darkened stars. It is based on the methods of Barnes (2009). The star is defined on a grid with the brightness at each point modeled as blackbody emission for a temperature described by Equation (1). After setting a series of fixed parameters, simuTrans employs emcee, a Markov Chain Monte Carlo (MCMC) ensemble sampler (Foreman-Mackey et al. 2013) to explore the posterior probability distribution for the remaining fitted parameters and determine their best-fit values.

The fixed parameters are as follows,

• 1.  
  Orbital period, P,
• 2.  
  Stellar mass, M_*,
• 3.  
  Stellar radius, R_*,
• 4.  
  Stellar rotation period, P_rot,
• 5.  
  Gravity-darkening exponent, β,
• 6.  
  Stellar oblateness, ${f}_{* }=1-\tfrac{{R}_{\mathrm{pole}}}{{R}_{\mathrm{eq}}}$,
• 7.  
  Limb-darkening coefficient, ${q}_{1}={({u}_{1}+{u}_{2})}^{2}$,
• 8.  
  Limb-darkening coefficient, ${q}_{2}=\tfrac{{u}_{1}}{2({u}_{1}+{u}_{2})}$,
• 9.  
  Stellar radius scaled by semimajor axis, ${R}_{* }/a$, and
• 10.  
  Planetary radius scaled by stellar radius, ${R}_{p}/{R}_{* }$,

where q₁ and q₂ are the re-parameterized quadratic limb-darkening coefficients outlined by Kipping (2013). We also assume a circular orbit, as the eccentricity has been shown to be very small: on the order of 10⁻⁴ (Benomar et al. 2014; Shporer et al. 2014; Esteves et al. 2015).

The fitted parameters are as follows,

• 1.  
  Time of mid-transit, T_o,
• 2.  
  Impact parameter, b,
• 3.  
  Angle between sky plane and stellar spin axis, φ, and
• 4.  
  Sky-projected angle between stellar spin axis and planetary orbital axis, λ.

While the values for fixed parameters 1–5 are taken directly from previous studies (see Table 2), the values for fixed parameters 6–10 were determined through an iterative MCMC fit. All of the short-cadence data for Kepler-13Ab were divided into four groups based on their year of observation. These were then phase-folded into four transits following the same method described in Section 2, and then fitted with simuTrans. Including the four previously listed fitted parameters, we fit for nine parameters total. We expected parameters 6–10 to remain constant, as the limb-darkening coefficients describe an intrinsic property of the star, and ${R}_{* }/a$ and ${R}_{p}/{R}_{* }$ are intrinsic to the planetary system. The parameter f_* is related to the star's rotational velocity, which should not change over the timescales considered. Therefore, after discerning that these parameters remained constant throughout, we calculated their average values and set these as fixed parameters when fitting the quarter-folded light curves. The values reported in Table 2 for the four fitted parameters are the mean values calculated from all of the quarter-folded transit fits.

Table 2.  Transit Model Parameters

Parameter	Value
P (days)a	1.763588
M* (M${}_{\odot }$)b	1.72
R* (R${}_{\odot }$)b	1.74
Prot (hr)c	22.5
β	0.25
f*	0.02 ± 0.01
q1	0.261 ± 0.017
q2	0.350 ± 0.040
R${}_{* }/a$	0.2205 ± 0.0013
Rp/R*	0.0905 ± 0.0004
To (days)d	120.56578 ± 0.000011
b	0.1960 ± 0.0013
φ (o)	13.0 ± 0.8
λ (o)	−27.9 ± 1.1

Notes.

^aFrom E15. ^bFrom Shporer et al. (2014). ^cFrom Barnes et al. (2011). ^dT_o is given in BJD—2454833 days.

Download table as:  ASCIITypeset image

Figure 3 shows the definitions of the various angles used to describe the configuration of the system. The angle of the stellar spin axis φ is measured from the sky plane down along the observer's line of sight (LOS) and is defined to be in the range $[-90^\circ$, 90°].⁷ The angle from the sky-projected planetary orbital axis to the sky-projected stellar spin axis, λ, is measured on the sky plane perpendicular to the LOS. It is defined to be in the range $[-180^\circ$, 180°]. To determine the true spin–orbit misalignment ψ in three-dimensional space (see Figure 3), the planet's orbital inclination must also be taken into account:

$$\begin{eqnarray}&&\cos {i}_{\mathrm{orb}}=\left(\displaystyle \frac{{R}_{* }}{a}\right)b.\end{eqnarray} \tag{ 3 }$$

The true spin–orbit misalignment can then be found by modifying Equation (1) of Benomar et al. (2014) to match our parameter definitions such that

$$\begin{eqnarray}&&\cos \psi =\sin \varphi \cos {i}_{\mathrm{orb}}+\cos \varphi \sin {i}_{\mathrm{orb}}\cos \lambda .\end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

In this section, we examine the results of the models fitted to the phase-folded and binned transits for the 12 quarters considered, as well as the implications on the future evolution of Kepler-13Ab. The fitted models are shown in Figure 4. There is a slight asymmetry in the residuals during egress in some of the models, which could be due to the fixed value of the gravity-darkening parameter β. As discussed in Section 1, it is unclear whether the host star is best described by a convective or radiative stellar envelope, and exploring the resultant influence on β lies outside the scope of this paper. We leave this work for future consideration and instead adopt the value $\beta =0.25$, which is consistent with previous studies of this system (e.g., Barnes et al. 2011; Masuda 2015).

Zoom In Zoom Out Reset image size

5.1. Spin–Orbit Precession

The evolution of the fitted parameters b, λ, φ, and T₀ are plotted in Figure 5. The top panel shows a clear decrease in the impact parameter over the time span of the short-cadence data quarters, which indicates that the inclination of Kepler-13Ab's orbit is undergoing temporal variation. We fit a linear model to this time series using a least-squares method and determine a rate of change of ${db}/{dt}\,=(-4.1\pm 0.2)\times {10}^{-5}$ day⁻¹. This is consistent with previously reported values; Szabó et al. (2012) find ${db}/{dt}=(-4.4\pm 1.2)\times {10}^{-5}$ day⁻¹ by individually fitting each transit in Kepler quarters 2 and 3, while Masuda (2015) fit individual transits in all short-cadence quarters (including quarter 13) and find $d| \cos ({i}_{\mathrm{orb}})| /{dt}=(-7.0\pm 0.4)\,\times {10}^{-6}$ day⁻¹. Equation (3) gives the relationship between $\cos {i}_{\mathrm{orb}}$ and b.

Zoom In Zoom Out Reset image size

The consistency between these results is reassuring considering the differences in our methods. We use phase-folded transit light curves, which considerably improve the signal-to-noise ratio of the data compared to individual transits, but do not wash out the clear trend in impact parameter over the timescale considered. We also utilize more than six times as much data as Szabó et al. (2012) by including Kepler quarters 7–12 and 14–17 in our analysis.

The temporal variation of the impact parameter has important implications, as first established by Szabó et al. (2012) and confirmed by both Masuda (2015) and this work. Namely, the variations are indicative of spin–orbit precession caused by the host star's quadrupole moment. Based on the rate of change in impact parameter, we expect Kepler-13Ab to become a non-transiting exoplanet in 75–85 years.

We also find a small rate of change in the time of mid-transit T₀ of $(-1.7\pm 0.1)\times {10}^{-2}$ s day⁻¹, as shown in the bottom panel of Figure 5. This is most likely not caused by spin–orbit precession and instead likely results from the value at which we fix the orbital period (see Table 2). The decrease in T₀ suggests that this period is too large by $3.6\times {10}^{-7}$ days, or about 0.03 s. Beyond this small correction we do not find any evidence for a temporal change in the period at the significance level of our measurements.

5.2. Spin–Orbit Misalignment

Because the angular momenta of the stellar spin and orbital motion have comparable magnitudes in the Kepler-13A system (due to the star's rapid rotation and the planet's small semimajor axis), they will both precess about the total angular momentum vector. As a result the angles i_orb, λ, and φ are all expected to undergo variation. Masuda (2015) computes the future evolution of these parameters and shows that neither λ nor φ should change by more than a few degrees on the relatively short timescale we consider. Indeed, the middle and bottom panels of Figure 5 exhibit no discernible temporal variation within error.

We therefore use the average values of b, λ, and φ to calculate the true spin–orbit misalignment ψ using Equations (3) and (4). We report a misalignment of $29^\circ \pm 1^\circ$ for this system. Our value is lower than that given in Masuda (2015), likely due to the difference in our fitted value of λ. In Figure 6 we use the average values of our fitted parameters to plot the mean transit model on top of a single phase-folded transit encompassing all transits in the 12 short-cadence quarters investigated.

Zoom In Zoom Out Reset image size

5.3. Seasonal Variations

The top panel of Figure 5 contains a few systematically offset values which are from quarters 8, 12, and 16. These three quarters correspond to the same "season" in each year of the Kepler observations meaning the telescope is rotated onto the same side and Kepler-13A is located on the same CCD each time. Because of this, the offset has previously been attributed to systematic rather than physical variation, both for Kepler-13Ab (Masuda 2015) and other planetary systems such as HAT-P-7b (Van Eylen et al. 2013). Van Eylen et al. (2013) suggest instrumental artifacts or unknown field crowding could be responsible, which would translate to an increase in flux that could contaminate the star. Much like the flux contributed by Kepler-13A's stellar companion, this has the effect of making the transit shallower, and if unaccounted for, will result in a larger impact parameter when fitting the gravity-darkened transit model.

To investigate the effect of flux dilution on the fitted impact parameter, we produce and compare quadratically limb-darkened transit models (without gravity-darkening) for various percentages of contaminating light and a range of impact parameter values. We use the quadratic limb-darkening model from Mandel & Agol (2002), as implemented with Python in the occultquad function⁸ from EXOFAST (Eastman et al. 2012). We use the same values of q₁, q₂, ${R}_{p}/{R}_{* }$, and ${R}_{* }/a$ as those in our gravity-darkened model (see Table 2), changing only the impact parameter and percentage of contaminating light. Excluding gravity darkening from this analysis provides a less robust but sufficiently accurate transit model for our investigation.

We first produce a quadratically limb-darkened transit model with zero additional light and a fixed impact parameter value of 0.205, which is approximately the value we would expect for quarter 8 based on the fit shown in Figure 5. We then calculate new diluted transit models with additional light between 0 and 1% and vary the impact parameter in the range [0.18, 0.23] for each. We introduce random noise to the data, with standard deviation equal to the median out-of-transit standard deviation of the transits phase-folded by quarter. We then calculate the ${\chi }^{2}$ value with respect to the initial model with fixed impact parameter and zero dilution.

Our results are presented in Figure 7, with black dots indicating the minimum ${\chi }^{2}$ for each dilution percentage. From this, we find that a percentage of contaminating light less than 1% is enough to cause the jump in impact parameter seen in quarters 8, 12, and 16, supporting the idea that systematic effects are indeed the cause of such seasonal variations. Determining the precise percentage of contamination in these quarters exceeds the scope of this paper; however, we treat the offset values as scatter while noting their likely origin.

Zoom In Zoom Out Reset image size

5.4. Secondary Eclipse Analysis

Despite strong evidence of spin–orbit precession from the Kepler transit data, it has not yet been ascertained whether similar variations are perceptible in the phase curve data. We therefore analyze the secondary eclipse as a function of time and describe our findings herein.

Taking into account the appropriate time binning, we model the secondary eclipses of all short- and long-cadence data. We use the same occultquad code described in Section 5.3, with no limb darkening or gravity darkening. We implement a simple MCMC process to fit for the secondary eclipse depth F_ecl for all phase curves, then perform a joint fit assuming constant depth and allowing the impact parameter to vary linearly with time. We also allow for a phase offset to account for the slight deviation in the time of mid eclipse identified in the transit data (see Section 5.1). Our results are presented in Table 3. Though at a much lower significance, the change in impact parameter we measure here is consistent with the value we report from the transit data.

Table 3.  Phase Curve Results

Parameter	Value
Fecl (ppm)	172.4 ± 0.7
db/dt	$(-3.2\pm 1.3)\times {10}^{-5}$ day−1
dT${}_{0,\mathrm{ecl}}$/dt	$(-1.7\pm 0.7)\times {10}^{-1}$ s day−1

Download table as:  ASCIITypeset image

This is the first time confirmation of spin–orbit precession has been identified in the secondary eclipse data. This new result serves as further evidence of the unique capabilities of Kepler owing to its long baseline and high photometric precision.

The prospect of asymmetric transit light curves resulting from spin–orbit misaligned planets orbiting gravity-darkened stars was initially discussed by Barnes (2009). Szabó et al. (2011) were the first to show that Kepler-13Ab displays signs of an asymmetric light curve, and Barnes et al. (2011), Masuda (2015), and Howarth & Morello (2017) each determined the net spin–orbit misalignment of the system using independent transit models. Transit duration variations due to precession of the planetary orbital axis were identified by Szabó et al. (2012) and confirmed by Masuda (2015).

Through analysis of the phase-folded transit light curves for 12 short-cadence data quarters, as well as the individual secondary eclipses, we similarly conclude that the detected change in impact parameter is indicative of a precessing planetary orbit and determine the net spin–orbit misalignment of the system. The calculated rates of change we report based on the transit and phase curve are self-consistent and agree with previously published values. We posit that the systematic offset in quarters 8, 12, and 16 identified in both our light curve analysis and that of Masuda (2015) can be explained by the addition of contaminating light from either instrumental artifacts or field crowding, as described by Van Eylen et al. (2013).

Though we are able to use a gravity-darkened transit model to characterize this system, the origin of its misalignment is not yet fully understood, and discerning its complex evolution will require additional short-cadence, high-precision observations like those that will be possible with the upcoming TESS mission. Similar observations of other planetary systems like Kepler-13Ab will also provide insight into the prevalence of such fascinating dynamical phenomena.

We thank Lisa Esteves for her guidance during the early stages of this work, as well as Dan Tamayo for helpful discussions. This work was supported in part by grants from the Natural Sciences and Engineering Research Council (NSERC) of Canada to R.J. Also, E.d.M. was in part funded by the Michael West Fellowship.