The flipped orbit of KELT-19Ab inferred from the symmetric TESS transit light curves

We extract the light curve for KELT-19Ab from the TESS full-frame images created at 2-minute cadence in Sector 7, which returned four full transits, using the python package $\mathrm{\tt eleanor}$ (Feinstein et al., 2019). From the options provided with the package, we select PCAFLUX as the flux values, which has the systematics removed via the subtraction of cotrending basis vectors (CBVs) provided by the Science Payload Operations Center (SPOC) pipeline (Jenkins et al., 2016). CBVs are obtained through Principal Component Analysis and represent systematic trends unique to each operational sector.

To correct for flux contamination in the light curve, we use information from the TESS target pixel file to remove excess flux from sources other than the target. Specifically, we use the CROWDSAP and FLFRCSAP attributes, which represent the fraction of flux within the aperture that comes from the target and the fraction of target flux captured within the aperture, respectively. We first subtract the excess flux, calculated as $(1-\mathrm{CROWDSAP})*$ median flux, from the raw flux values. Then, we divide the corrected flux values by FLRCSAP to obtain a final flux that accounts for the effects of crowding.

The definition for true obliquity $\psi$ used in this paper is given using orbital inclination $i_{p}$, stellar inclination $i_{*}$, and sky-projected obliquity $\lambda$ in Figure 1, as well as the equation

Here, orbital inclination and stellar inclination are defined as the angle between the planetary orbital axis and stellar rotational axis with the line of sight respectively.

Gravity darkening is the latitude dependent flux variation of a star, due to its rapid rotation. This effect causes a transiting planet to occult areas of differing surface brightness at different points in its transit, allowing the obliquity to be constrained by fitting the anomaly in the light curve. The technique has been used to constrain obliquity for planets including KELT-9b (Ahlers et al., 2020a), KOI-89.01 and KOI-89.02 (Ahlers et al., 2015), KOI-368.01(Zhou & Huang, 2013), Kepler-13Ab (Barnes et al., 2011; Masuda, 2015), HAT-P-7b (Masuda, 2015), HAT-P-70b (Zhou et al., 2019), PTFO 8-8695 (Barnes et al., 2013), MASCARA-1b (Hooton et al., 2022), MASCARA-4b (Ahlers et al., 2020b), WASP-33b (Dholakia et al., 2022), and WASP-189b (Lendl et al., 2020; Deline et al., 2022).

Gravity darkening is usually modelled by the von Zeipel theorem (von Zeipel, 1924) as

The degree of anomaly is predominantly governed by the stellar effective temperature and rotational velocity from Equation 2, and the general shape of the gravity darkened light curve is governed mainly by the sky-projected obliquity and the stellar inclination. For KELT-19Ab, the light curve retains its symmetry regardless of $i_{*}$, owing to $\lambda=-179.7^{\circ+3.7^{\circ}}_{\,\,-3.8^{\circ}}$. As $i_{*}$ decreases (approaching pole on), the light curve gets progressively deeper towards mid-transit, as the planet transits through the brighter region near the pole as illustrated in Figure 2. The lower limit of $i_{*}=34^{\circ}$ is given by requiring that KELT-19A rotates slower than the empirical limit $v\approx 150$ km/s for an Am star (Siverd et al., 2018). The difference in transit depth mid-transit between $i_{*}=90^{\circ}$ and $i_{*}=34^{\circ}$ is around 600 ppm.

On top of the ways in which $\lambda$ and $i_{*}$ affect the shape of the gravity darkened light curve, it is also worth noting that gravity darkening is usually used in tandem with other methods of constraining $\lambda$, in order to disentangle the degeneracy between the prograde and retrograde solutions. In the hypothetical case that $i_{*}=90^{\circ}$ and the light curve is symmetric, one cannot tell $\psi\approx 0^{\circ}$ from $\psi\approx 180^{\circ}$ by gravity darkening alone without prior knowledge of $\lambda$. This is because gravity darkening is a latitude dependent effect which can determine the stellar inclination but not the direction of the host star’s spin with respect to the planetary orbit. To get the sense of stellar spin, observations of longitude dependent effects such as the Rossiter-Mclaughlin effect are required.

We use the gravity darkening model implemented in the python package $\mathrm{\tt PyTransit}$ (Parviainen, 2015) to estimate the model parameters. For each discretized point on the stellar surface, the model evaluates the temperature based on Eq. 2 and derives the flux given an emission spectrum and a response function of the instrument. We use the synthetic spectra from the PHOENIX library (Husser et al., 2013) and the TESS instrument response function provided in Ricker et al. (2015).

The parameters for the gravity darkening model are listed in Table 2, including the sky-projected obliquity $\lambda$, stellar rotational velocity $v\sin i$, stellar radius $R_{*}$, stellar effective temperature $T_{\mathrm{eff}}$, gravity darkening exponent $\beta$, stellar inclination $i_{*}$, as well as regular transit parameters which are radius ratio $R_{p}/R_{*}$, impact parameter $b$, stellar density $\rho$, and the two limb-darkening coefficients $u_{1},u_{2}$. Some crucial parameters are also derived from these parameters, which include the stellar rotational period, $P_{\mathrm{rot}}=2\pi R_{*}\sin i/v\sin i$, and the stellar oblateness $f$, approximated in the model to be $f=3\Omega_{*}^{2}/8\pi G\rho$. To sample from the posterior distribution, we conduct a Markov chain Monte Carlo (MCMC) run using Python code $\mathrm{\tt emcee}$, which implements an affine-invariant ensemble sampler (Foreman-Mackey et al., 2013). We chose 50 walkers, thinning of 10 steps for the entirety of 500,000 steps, and the first 10,000 steps are discarded as burn-in.

From the analysis, we obtain that the stellar inclination $i_{*}=97^{\circ+29^{\circ}}_{\,\,-34^{\circ}}$ and $i_{*}=96^{\circ+33^{\circ}}_{\,\,-37^{\circ}}$ for the two gravity darkening models (GD7400 and GD7600) respectively. They are fully consistent with $33.5^{\circ}<i_{*}<146.5^{\circ}$ suggested by Siverd et al. (2018), assuming that KELT-19A rotates slower than the empirical limit, $v_{\mathrm{eq}}<150\,\mathrm{km/s}$. Fit and derived parameters are listed in Table 3 with the corner plot in Appendix A, and the best-fit models to the data is presented in Figure 3.

The obtained stellar inclination results in $\psi=157^{\circ+15^{\circ}}_{\,\,-20^{\circ}}$ with GD7400 and $\psi=155^{\circ+17^{\circ}}_{\,\,-21^{\circ}}$, with GD7600, both of which place stronger constraints than the $1\sigma$ constraint of $119^{\circ}<\psi<180^{\circ}$ in Siverd et al. (2018). At 2 $\sigma$, the constraint with the two models becomes $121^{\circ}<\psi<177^{\circ}$ and $119^{\circ}<\psi<177^{\circ}$ respectively. The posteriors for both the obliquity and stellar inclination are shown in Figure 4. The result is in favor of a flipped orbit, and KELT-19Ab likely has the largest obliquity among currently known exoplanets, when compared to $124^{\circ}\pm 6.0^{\circ}$ of K-290b and c (Hjorth et al., 2021).

For KELT-19Ab, there are discrepancies in the reported values of the radius ratio, impact parameter and semi-major axis between observations by KELT (Siverd et al., 2018), CHEOPS (Garai et al., 2022) and TESS (Yang et al., 2022) as listed in Table 3. The values obtained in this work for these parameters is consistent with the analysis of Yang et al. (2022) within $1\sigma$. Yang et al. (2022), who explored the role of atmosphere in this discrepancy but reached an inconclusive result, mentioned that a transit fit with inclination and semi-major axis fixed to KELT values provides a fit with slightly larger reduced $\chi^{2}$ but still consistent with TESS data, although not the preferred solution. It is therefore likely that the discrepancy is due to some degeneracy.

Barnes (2009) showed that when $\psi$ and $|\lambda|=0^{\circ}$ or $180^{\circ}$ (i.e. the light curve is symmetric and U-shaped), the radius ratio of a Jupiter sized planet around an gravity darkened Altair (A7V star) with $b=0.9$ can be overestimated by almost $30\%$ with a regular transit model. This is because in such a gravity darkened system, larger impact parameter causes the planet to transit nearer to the brighter pole, resulting in larger transit depth. Overestimation is less pronounced at lower impact parameters, where increased projected area around the oblate equator makes up for the area being darker. Such an overestimation might be able to explain the trend of larger impact parameter coupled with larger radius ratio seen in Table 3. If we then assumed larger estimates of impact parameter in previous literature as the ground truth and perform a fit with gravity darkening, radius ratio should become smaller than in literature.

When impact parameter is fixed to value in Siverd et al. (2018), we obtain $0.09870^{+0.00073}_{-0.00084}$ with GD7400 and $0.09863^{+0.00075}_{-0.00087}$ with GD7600, which are both smaller than the values in Siverd et al. (2018) by 6.7 $\sigma$. The values, however, is marginally smaller but consistent with the $0.09940\pm 0.00048$ obtained by Yang et al. (2022) who did the same analysis using TESS light curve and a quadratic model. When fixed to values from Garai et al. (2022), we obtain $0.09743^{+0.00067}_{-0.00076}$ and $0.09741^{+0.00070}_{-0.00080}$ respectively, which are slightly smaller than Garai et al. (2022), but still consistent within 1$\sigma$. It is therefore difficult to conclude that gravity darkening is the lone culprit of the reported discrepancies, although it must be playing a role in making the problem degenerate.

It should also be noted, however, that the said overestimation of radius ratio happens only if the larger impact parameters results in occultation of the stellar surface at higher latitude. If the orbit was completely polar ($\psi=90^{\circ}$), larger impact parameter will cause the occultation at lower latitude, which must instead lead to an underestimation of radius ratio (Barnes, 2009). This means that the discrepancies can in principle be used to bolster the conclusion about the orbit of planets with symmetric gravity darkened light curves.

As evident from the posterior distribution, the obtained precision in $\psi$ is limited by the precision in $i_{*}$, which is one of the primary parameters fitted with the gravity darkening model. Improved observational precision will therefore result in stronger constraints on $i_{*}$ and thus $\psi$. With double the number of data points, we obtain the photometric precision per two-minute-binned data of $\sim~{}600$ ppm, which is tantamount to the expected difference in transit depth for the polar and flipped orbit scenario. KELT-19Ab hence remains as an important follow up target for TESS.

Comparing the transit depth over different wavelengths can in principle constrain obliquity, although the required precision is much higher. The degree of gravity darkening varies over different wavelengths, and stronger gravity darkening is expected at shorter wavelengths, resulting in more pronounced pole to equator contrast (Barnes, 2009). When comparing the TESS and CHEOPS passband for example, KELT-19Ab’s transit is deeper on the order of $10^{1}$ ppm with TESS than with CHEOPS, if the planet transits through the darker equator. The transit will instead be deeper on the order of $10^{2}$ ppm if the planet transits through the pole. Observations at longer wavelengths where limb-darkening is reduced might also be helpful, but this comes at the expense of reduced gravity darkening altogether.

We also consider the chance that changes in impact parameter $b$ or sky-projected obliquity $\lambda$ over time can tell a polar orbit from a flipped one, as the rate of nodal precession is faster for planets on polar orbits than a flipped one. (e.g. Watanabe et al., 2022). The two parameters are expressed as equations of time as follows.

Here, $\theta(t)$ is the nodal angle, which expresses how much the orbital axis has precessed around the stellar spin axis, and is expressed,

We calculate $b(t)$ and $\lambda(t)$ for KELT-19Ab, assuming $J_{2}=(1.36^{\mathrm{+0.15}}_{\mathrm{-0.12}})\times 10^{-4}$ taken from Watanabe et al. (2022) in their analysis of WASP-33b, in which the host WASP-33 has a similar $T_{\mathrm{eff}}\approx 7430$ K and $v\mathrm{sin}i_{*}\approx 85.6$ km/s as KELT-19A. Under this assumption, impact parameter for the polar ($i_{*}\approx 34^{\circ}$) and flipped ($i_{*}\approx 97^{\circ}$) case will only deviate by around 0.04 and by about a degree for sky-projected obliquity in the next three decades, neither of which are deviations large enough to be detected.

The large obliquity of KELT-19Ab is consistent with the recent finding due to Hierarchical Bayesian Analysis that there is no statistically significant preference towards polar planets, and that obliquity can be distributed more uniformly (Siegel et al., 2023; Dong & Foreman-Mackey, 2023). These studies refuted the frequentist test that showed that the observed distribution of misaligned planets (defined as $\mathrm{cos}\psi<0.75\,\,\mathrm{or}\,\,\psi>41^{\circ}$) is significantly different from a uniform distribution between $41^{\circ}<\psi<180^{\circ}$ (Albrecht et al., 2021).

For the purpose of Hierarchical Bayesian Analysis, both Siegel et al. (2023) and Dong & Foreman-Mackey (2023) concluded that measurements of $\lambda$ was already sufficient to establish the uniformity of the inferred distribution. $\lambda$ measured in Siverd et al. (2018) was already incorporated in both studies, and we therefore decide not to replicate the the analysis with our constraints on $\psi$. The ability of $\lambda$ in the context of Hierarchical Bayesian Analysis, nevertheless, does not undermine the importance of measuring $\psi$ for individual systems. Especially when flipped orbits are in interest, large $\lambda$ does not guarantee a flipped orbit (Fabrycky & Winn, 2009; Xue & Suto, 2016). Then, measuring $\lambda$ is insufficient, unlike when we have $\lambda\approx\psi$ regardless of $i_{*}$, where $\lambda\approx 90^{\circ}$. Therefore, measurements of $\psi$ is indispensable when searching for flipped orbits.