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

Dozens of planets are now discovered with large orbital obliquity, and have become the proof for the dynamical evolution of planetary orbits. In the current samples, there is an apparent clustering of planets around $90^\circ$, and also an absence of planets around $180^\circ$ although the latter is expected by some theories. Statistical extrapolation using Hierarchical Bayesian Analysis have recently refuted the significant clustering around $90^\circ$ and suggested that the distribution may actually be broader. In this work, the symmetric TESS transit light curve of KELT-19Ab is analyzed using gravity darkening to measure its true obliquity. Its large sky projected obliquity $\lambda = -179.7^{\circ+3.7^\circ}_{\,\,-3.8^\circ}$ makes KELT-19Ab the only currently known planet with obliquity potentially close to $180^\circ$. We apply spectroscopic constraints on $v\mathrm{sin}i$ and $\lambda$ as well as theoretical constraints on the limb-darkening coefficients to find that the KELT-19Ab's obliquity is $\psi = 155^{\circ+17^\circ}_{\,\,-21^\circ}$, in favor of a flipped orbit. The result is consistent with the statistically inferred uniformity of obliquity distribution, and also highlights the applicability of the gravity darkening technique to symmetric light curves.


INTRODUCTION
In the past decade of exoplanetary research, orbits of many exoplanets have been discovered misaligned with respect to the rotation axes of their host stars.The degree of misalignment between the orbital axis of exoplanets and the spin axis of their host stars  (spin-orbit angle or obliquity) has then become an important probe to investigate the dynamical history of these planets.
In the current samples of obliquity , there is an apparent clustering of planets at around  = 90 • .The trend was first mentioned by Albrecht et al. (2021), and also presented in Attia et al. (2023).However, multiple studies have more recently refuted such a preponderance of polar planets via Hierarchical Bayesian Analysis (Siegel et al. 2023;Dong & Foreman-Mackey 2023), arguing that the observed peak is statistically insignificant.Such a context puts the absence of planets with obliquity close to 180 • into question.
Previous measurements of  for planets with  (the sky-projection of ) greatly exceeding 90 • show that their orbits are actually closer to polar than flipped 1 (or at least on the verge of the two classes), as ★ yugo6581@g.ecc.u-tokyo.ac.jp 1 Throughout the paper, 60 • <  < 120 • is referred to as polar, and  > 120 • as flipped.Retrograde orbit ( > 90 • ) or even large sky-projected obliquity ( ≈ 180 • ) does not ensure a flipped orbit, due to the uncertainty in the direction of the stellar rotational axis (Fabrycky & Winn 2009;Xue & Suto 2016).listed in Table 1, which vaguely hints to the rarity of these "flipped orbits".Theoretically, however, excitation of obliquity to such an extreme have been predicted with variations of Kozai-Lidov mechanism including the concurrent evolution of stellar spin (Storch et al. 2014;Storch & Lai 2015;Anderson et al. 2016;Storch et al. 2017).Storch et al. (2014) found that chaotic evolution of stellar spin due to spin-orbit coupling between the precession frequencies of the oblate star and the Kozai frequencies result in extreme obliquities, some of them close to 180 • .The result owes to the fact that obliquity in this context is the actual angle between the stellar spin axis and the planetary orbital axis, instead the conventional angle between the initial and the evolved planetary orbital axis, which is practically limited by the Kozai angle = 140.8• and assumes primordial starplanet alignment.
Since the direction of stellar spin significantly alters the final obliquity after Kozai-Lidov mechanism, primordial misalignment between the stellar spin axis and the planetary orbital axis can also produce  ≈ 180 • (Anderson et al. 2016;Vick et al. 2022).Observational evidence for primordial misalignments exist in forms of and inner disks misaligned with respect to the outer disks (Marino et al. 2015;Ansdell et al. 2020).Misaligned coplanar multi-planet systems (Huber et al. 2013;Hjorth et al. 2021) are also potentially the legacy of primordial misalignment.
Hot Jupiter explored in this paper, has the largest sky projected obliquity  = −179.7 •+3.7 • −3.8 • of currently known exoplanets, and is confirmed to have a retrograde orbit and potentially a flipped one, depending on the stellar inclination.It also belongs to the class of hot Jupiters around host stars above the Kraft break ( eff = 6250 K).Beyond this point, convective layers in the host become absent and planets are unlikely to be tidally realigned (Winn et al. 2010).The outcome of any violent orbital evolution in KELT-19Ab is therfore likely retained.
The host is an Am star with  eff = 7500±110 K and sin = 84.2±2.0 km/s (Siverd et al. 2018), making it an excellent target for analysis with gravity darkening, whose magnitude is more pronounced in hotter and faster rotating stars.The distant companion KELT-19B is a late-G9V/early-K1V star at a projected separation of 160 au (Siverd et al. 2018), providing the necessary grounds for Kozai-Lidov mechanism to take place.
In this paper, we use gravity darkening to estimate the obliquity of hot Jupiter KELT-19Ab, from its symmetric TESS (Ricker et al. 2015) light curve.We apply spectroscopic constraints derived by Siverd et al. (2018) to gauge the extent of gravity darkening apriori and also theoretical constraints on limb-darkening coefficients to test the applicability of the gravity darkening method to symmetric transits.Gravity darkening for KELT-19Ab have previously been explored by Garai et al. (2022) using CHEOPS data, who concluded a non-detection of gravity darkening.They speculated that this was because the degree of anomaly expected for KELT-19Ab's light curve with gravity darkening is about 10 ppm.We argue that the degree of anomaly is 600 ppm instead, and that their non-detection was due to an apparent underestimation of the host's rotational velocity.
The paper is organized as follows.In Section 2, we describe the methods used in this work, including data reduction, gravity darkening, and MCMC including the choice of priors.Section 3 and 4 provides the results and the discussion respectively.

TESS photometry
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 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.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 − 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.

Geometry
The definition for true obliquity  used in this paper is given using orbital inclination   , stellar inclination  * , and sky-projected obliquity  in Figure 1, as well as the equation cos = cos p cos * + sin p sin * cos. ( 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 usually modelled by the von Zeipel theorem (von Zeipel 1924) as The transit is shallower in the former case, where the planet transits through the darker equator of the gravity darkened host star (Top right).The transit gets progressively deeper towards mid-transit for the latter case, as the planet transits through the brighter region near the pole (Bottom right).The difference in transit depth mid-transit is around 600 ppm.
where  () and () are the temperature and surface gravity at a given latitude,  pole and  pole are the temperature and surface gravity at the poles respectively. is the gravity darkening constant, which is equal to 0.25 for a barotropic star in a strict radiative equilibrium.Surface gravity () at at a given latitude can be evaluated with the following equation where  is the universal gravitational constant,  * is the stellar mass, Ω * is the stellar rotation rate,  () is the distance from the stellar center to the point of interest, and  ()sin hence denotes the distance from the stellar rotational axis to the same point.r and r⊥ are both unit vectors pointing to the point of interest, one from the star's center and the other from the rotational axis.These equations relate the magnitude of centrifugal force at a given latitude to temperature, which will be the lowest at the equator, where the star also becomes the darkest.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  * , owing to  = −179.7 •+3.7 • −3.8 • .As  * 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  * = 34 • is given by requiring that KELT-19A rotates slower than the empirical limit  ≈ 150 km/s for an Am star (Siverd et al. 2018).The difference in transit depth mid-transit between  * = 90 • and  * = 34 • is around 600 ppm.
On top of the ways in which  and  * 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 , in order to disentangle the degeneracy between the prograde and retrograde solutions.In the hypothetical case that  * = 90 • and the light curve is symmetric, one cannot tell  ≈ 0 • from  ≈ 180 • by gravity darkening alone without prior knowledge of .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 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).

MCMC parameters and priors
The parameters for the gravity darkening model are listed in Table 2, including the sky-projected obliquity , stellar rotational velocity  sin , stellar radius  * , stellar effective temperature  eff , gravity darkening exponent , stellar inclination  * , as well as regular transit parameters which are radius ratio   / * , impact parameter , stellar density , and the two limb-darkening coefficients  1 ,  2 .Some crucial parameters are also derived from these parameters, which include the stellar rotational period,  rot = 2 * sin / sin , and the stellar oblateness  , approximated in the model to be  = 3Ω 2 * /8 .To sample from the posterior distribution, we conduct a Markov chain Monte Carlo (MCMC) run using Python code 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.

Specstroscopic constraints
We apply spectroscopic constraints on the stellar radius, rotational velocity and sky-projected obliquity with normal priors, taking the mean values and uncertainties from Siverd et al. (2018) as listed in Table 2.
A normal prior on rotational velocity constrains the degree of gravity darkening and hence the degree of anomaly.As discussed by Barnes (2009), it then helps distinguish a gravity darkened and symmetric light curve like the one of KELT-19Ab from a regular non gravity darkened light curve.In the exploration of gravity darkening for KELT-19Ab by Garai et al. (2022), they have used Ω/Ω crit = 0.23, the ratio of rotational velocity to the break-up velocity, as a constraint instead of sin.This seems to be an underestimation, given that ratio of sin ≈ 85 km/s to the break-up velocity for A-type stars ≈ 250 km/s already gives 0.34, which might have resulted in their nondetection of gravity darkening.
A normal prior of  = −179.7±3.8 • on the sky-projected obliquity constrains the expected shape of the gravity darkened light curve for the reasons discussed in the previous section.It also serves to avoid any nonphysical or discrepant values with the results from Doppler tomography as mentioned in Masuda (2015).

Limb darkening and gravity darkening coefficients
As obliquity is estimated from subtle anomalies in transit shape, the result of gravity darkening analysis is also sensitive to the choice of limb darkening coefficients (Masuda 2015;Ahlers et al. 2015).This is especially true for planets with || ≈ 0 • or 180 • , because the transit will always be symmetric and gravity darkening will only be imprinted as the difference in transit depth mid-transit, which is exactly the case for .When a gravity darkened light curve retains such U-shaped symmetry, limb-darkening coefficients, if let free, can alter the shape of ingress and egress, causing the problem to be fully degenerate.Siegel et al. (2023) ran a Mote Carlo simulation to show that the aforementioned degeneracy is the cause of an observational bias in gravity darkening.This was done by showing that a regular transit model with free limb-darkening coefficients can provide as good of a fit to any gravity darkened light curves that are U-shaped and symmetric, without the need for a gravity darkening model.Their conclusion was that the method is heavily sensitive to polar or nearly polar planets, which are the limited configurations that potentially produce asymmetric light curves depending on  (See Barnes (2009) for possible types of asymmetry with gravity darkening).
However, if gravity darkening is expected to occur for a certain system, choice of the gravity darkening model over a regular transit model is not out of the quality of fit, but rather out of the necessity to account for an existing physical phenomenon.In a case a regular transit model is used for gravity darkened light curves that are U-shaped and symmetric, parameters such as radius ratio can be overestimated at the expense of setting limb-darkening coefficients free, because transit depth changes with stellar inclination as seen in Figure 2. Therefore, prior constraints on limb-darkening coefficients are necessary.
A common choice is to rely on theoretical values, be it fixing to catalog values as has been done in Ahlers et al. (2015), or calculating theoretical priors based on the theoretical stellar intensity profile as in Deline et al. (2022) 2 .In this work, we decided to apply normal priors with standard deviation of 0.03 on the limb-darkening coefficients with theoretical values given by Claret (2017) for TESS bandpass as the mean, given the stellar effective temperature of 7400 K and 7600 K, as well as surface gravity of 4.0 and metallicity of 0. This is based on  eff = 7500 ± 110 K, log * = 4.127 ± 0.029 g/cm 3 and [Fe/H] = −0.12± 0.51 of KELT-19A from Siverd et al. (2018), and on the intent to test the robustness of our results with different choice of limb-darkening coefficients.We hence run two sets of analysis on the light curve, first fixing the stellar effective temperature to 7400 K and then 7600 K.
The gravity darkening constant is fixed to 0.25.Previous studies have shown that the choice of this constant can lead to conflicting results (e.g.Ahlers et al. 2020b).To address this, we examined the impact of the gravity darkening exponent on the obliquity of KELT-19Ab by also performing a fit using the theoretical value of around 0.2 proposed by Claret (2016) for stars with effective temperature, metallicity and surface gravity similar to that of KELT-19A.We find that varying this parameter does not affect the conclusion that KELT-19Ab's orbit could be flipped, and hence leave it to 0.25.

Remaining parameters
We let the stellar inclination to vary freely between 0 • and 180 • , where 0 • denotes a pole-on configuration.We also apply uniform priors to non-gravity darkening parameters   / * ,  and  * in the range specified in Table 2.

Obliquity
From the analysis, we obtain that the stellar inclination  * = 97 •+29 • −34 •   and  * = 96 •+33 • −37 • for the two gravity darkening models (GD7400 and GD7600) respectively.They are fully consistent with 33.5 • <  * < 146.5 • suggested by Siverd et al. (2018), assuming that KELT-19A rotates slower than the empirical limit,  eq < 150 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.

Radius ratio, impact parameter and semi-major axis
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.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  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  and || = 0 • or 180 • (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  = 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 .The values, however, is marginally smaller but consistent with the 0.09940±0.00048obtained 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.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 ( = 90 • ), 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.

Possibility for stronger constraints
As evident from the posterior distribution, the obtained precision in  is limited by the precision in  * , which is one of the primary parameters fitted with the gravity darkening model.Improved observational precision will therefore result in stronger constraints on  * and thus .With double the number of data points, we obtain the photometric precision per two-minute-binned data of ∼ 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  or sky-projected obliquity  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,  () is the nodal angle, which expresses how much the orbital axis has precessed around the stellar spin axis, and is expressed, where  2 is the stellar gravitational quadruple moment.The formulation assumes that the angular momentum due to stellar rotation is much larger than that due to the planetary orbit, causing the orbital axis to precess around the stellar rotational axis.Such assumption may not always hold for hot Jupiters with sufficient orbital angular momentum (Barnes et al. 2013), although this does not significantly affect the precession rate.We calculate () and () for KELT-19Ab, assuming  2 = (1.36 +0.15  −0.12 )×10 −4 taken from Watanabe et al. (2022) in their analysis of WASP-33b, in which the host WASP-33 has a similar  eff ≈ 7430 K and sin * ≈ 85.6 km/s as KELT-19A.Under this assumption, impact parameter for the polar ( * ≈ 34 • ) and flipped ( * ≈ 97 • ) case will only deviate by around 0.04 and by about a degree for skyprojected obliquity in the next three decades, neither of which are deviations large enough to be detected.

Kozai-Lidov mechanism
Considering that a suitable companion exists in KELT-19B, which is a late-G9V/early-K1V star located at a projected separation of 160 au (Siverd et al. 2018), Kozai-Lidov mechanism is a plausible hypothesis to explain the orbital evolution of KELT-19Ab.The timescale for Kozai-Lidov mechanism can be calculated using the following approximation.
where   ,   , and   are the mass, the orbital period and the eccentricity of the companion.To estimate the this timescale for KELT-19Ab, we assume 0.9 ⊙ for the late-G9V/early-K1V companion KELT-19B from Pecaut & Mamajek (2013), and take the projected separation of the companion of 160 au from Siverd et al. (2018) to calculate   .Considering a full range of orbital eccentricity of KELT-19B and accounting for the underestimation of the actual separation with respect to the projected separation, we obtain that the timescale does not exceed 0.2 Gyr, which is well below the estimated system age of 1.1 Gyr (Siverd et al. 2018).Hence, the estimated timescale is consistent with the hypothesis that Kozai-Lidov mechanism played a role in KELT-19Ab's orbital evolution.
If KELT-19Ab's orbit is flipped and is beyond 140.8 • (known as the Kozai angle), however, traditional application of Kozai-Lidov mechanism (Fabrycky & Tremaine 2007, e.g.) fails to explain such large obliquity.This is because the libration of pericenter essential in inducing the oscillation of eccentricity and inclination only occurs for systems below this critical angle.This angle is therefore the upper limit of obliquity resulting from Kozai-Lidov mechanism (lower limit exists similarly for prograde planets at 39.2 • ), where obliquity in this context is defined to be the angle between the initial (pre-Kozai) and the evolved (post-Kozai) planetary orbital axis.Numerical simulations show that octupole level effects allow for some dispersion around these angles (Petrovich 2015;Petrovich & Tremaine 2016), but not to the extent that planets with  ≈ 180 • are expected.
Meanwhile, obliquites beyond the Kozai angle can be achieved when the concurrent evolution of stellar spin is taken into account during Kozai-Lidov mechanism.Storch et al. (2014) found that chaotic evolution of stellar spin due to spin-orbit coupling between the precession frequencies of the oblate star and the Kozai frequencies result in extreme obliquities, some of them close to 180 • .The result owes to the fact that obliquity in this context is no longer the angle between the pre and post-Kozai planetary orbital axis assuming primordial alignment, but rather the angle between the stellar spin axis and the planetary orbital axis.The likeliness of a star experiencing this chaotic spin evolution is quantified in Storch & Lai (2015) by the adiabacity parameter , where Ω pl and Ω ps are the precession rate of the planet and star around the companion respectively, and  on the order of unity (i.e.Ω pl ≈ Ω ps ) is the required criterion for a chaotic evolution.Relevant parameters are the stellar moment of inertia constant  * , 'rotational distortion constant'   which is one third the Love number  2 , stellar raidus  * , stellar rotation rate Ω★ = Ω ★ /(  * / 3 * ) 1/2 , companion and planetary mass   and   as well as semi-major axis   and , and the initial mutual inclination between them  0 lc , which we use 140.8 • to calculate the upper bound of .We also follow Storch & Lai (2015) and adopt the canonical values  * ≈ 0.1 and   ≈ 0.05 (Claret & Gimenez 1992).We obtain  ≈ 1.7 for KELT-19Ab when semi-major axis of 1 au is assumed, but a gas giant like KELT-19Ab is unlikely to have formed at such a short distance considering the location of snowline (Kennedy & Kenyon 2008).The parameter is very sensitive to the assumed semi-major axis and quickly shoots up ( ≈ 5 × 10 4 when  = 10 au), likely indicating that a chaotic evolution of stellar spin is not a suitable mechanism to explain the obliquity of KELT-19Ab.
Even if the direction of stellar spin remains unchanged for the duration of Kozai-Lidov mechanism, flipped orbits beyond the Kozai angle can also result from primordial misalignment between the stellar spin axis and the planetary orbital axis (Anderson et al. 2016;Vick et al. 2022).Misaligned inner disks (Marino et al. 2015;Ansdell et al. 2020) provide strong observational evidence for primordial misalignments, and misaligned coplanar multi-planet systems (Huber Theoretically, gravitational interactions between host stars, protoplanetary disks, and inclined binary companions are a way to form primordial misalignments (Batygin & Adams 2013;Lai 2014;Spalding & Batygin 2014;Zanazzi & Lai 2018).Zanazzi & Lai (2018) showed, however, that when taken into account spin-orbit coupling between the host star and the planet as well as planet-disc interactions, excitation of primordial misalignment can either be completely suppressed or at best significantly reduced for hot Jupiters whether formed in situ or migrated inwards via type-II migration.There are nonetheless several other ways to form primordially misaligned planets including the nonaxisymmetric collapse of the molecular cloud core (Bate et al. 2010;Takaishi et al. 2020), magnetic star-disc interactions (Lai et al. 2011), or dissipation of internal-gravity waves in the host star (Rogers et al. 2013), which when combined with Kozai-Lidov mechanism may explain the obliquity of KELT-19Ab.

Other mechanisms
Coplanar high-eccentricity migration (Li et al. 2014;Petrovich 2015) is another theory predicting the formation of planets with  ≈ 180 • , where secular gravitational interactions between two coplanar and eccentric planets cause their orbits to potentially flip.In such a case, even when the orbit of a planet initially begins prograde with respect to the orbit of the distant companion, the secular interaction causes the orbit to completely flip.For such an orbital flip to occur, the semimajor axis of the inner planet must be large enough to avoid being tidally disrupted from the extreme eccentricity growth preceding the flip, which is hardly achieved (Petrovich 2015;Xue & Suto 2016). 3ue et al. ( 2014) have also suggested that tidal dissipation of inertial waves in the convective envelope of the star can temporarily realign the planet's orbit to  ≈ 180 • .However, KELT-19A with the  eff = 7500 ± 110 K is expected to have a radiative envelope rather than a convective one, and tidal realignment is improbable.While limited to systems in a very dense stellar cluster, another proposed mechanism fly-bys of a prograde and equal-mass perturber (Breslau & Pfalzner 2019).

Implications on the inferred obliquity distribution
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 cos < 0.75 or  > 41 • ) is significantly different from a uniform distribution between 41 • <  < 180 • (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  was already sufficient to establish the uniformity of the inferred distribution. 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 .The ability of  in the context of Hierarchical Bayesian Analysis, nevertheless, does not undermine the importance of measuring  for individual systems.Especially when flipped orbits are in interest, large  does not guarantee a flipped orbit (Fabrycky & Winn 2009;Xue & Suto 2016).Then, measuring  is insufficient, unlike when we have  ≈  regardless of  * , where  ≈ 90 • .Therefore, measurements of  is indispensable when searching for flipped orbits.

CONCLUSION
In this work, KELT-19Ab's symmetric light curve from TESS was analyzed with gravity darkening to estimate its obliquity.We find that with constraints on stellar rotational velocity, sky-projected obliquity as well as limb-darkening coefficients, KELT-19Ab's obliquity is  = 155 •+17 • −21 • (119 • <  < 177 • at 2), in favor of a flipped orbit.The finding is consistent with the uniformity of the inferred obliquity distribution presented in both Siegel et al. (2023) and Dong & Foreman-Mackey (2023) using Hierarchical Bayesian Analysis.
From theoretical standpoints, formations of flipped orbits are rare but possible with Kozai-Lidov mechanism even for the most extreme cases, when accounting for the concurrent evolution of stellar spin (Storch et al. 2014;Storch & Lai 2015, e.g.) or primordial misalignments between the planetary orbital axis and stellar spin axis (Anderson et al. 2016;Vick et al. 2022).For KELT-19Ab, the latter case of Kozai-Lidov mechanism coupled with primordial misalignments is a plausible hypothesis.
We also investigate the previously discussed but unresolved discrepancy in transit depth, impact parameter and semi-major axis which might be explained by the gravity darkened transit of  Although the trend of larger impact parameters coupled with larger transit depth precisely matches the characteristic of degeneracy caused by gravity darkened transits of planets with  and  both nearing 180 • (Barnes 2009), the results were inconclusive as to whether gravity darkening was the lone culprit of the discrepancies in .Similar discrepancies, however, if more pronounced might be another way to tell a flipped orbit from a polar one in other systems.
The above results highlight the applicability of the gravity darkening method to symmetric light curves, which was previously set aside the visibly asymmetric ones.As discussed by Barnes (2009), it is nonetheless crucial that a constraint on stellar rotational velocity is given apriori to gauge the degree of gravity darkening.Constraints on sky-projected obliquity is also vital to limit the possible shapes of gravity darkened light curves, and on limb-darkening coefficients to reduce the degeneracy in gravity darkening.MNRAS 000, 1-10 (2015)

Figure 2 .
Figure 2. The gravity darkened light curve simulated for KELT-19Ab.The light curve retains its symmetry regardless of  * , owing to sky-projected obliquity  = −179.7 •+3.7 • −3.8 • .The blue and orange lines represent the light curve for when  * = 90 • and  * = 34 • respectively.The transit is shallower in the former case, where the planet transits through the darker equator of the gravity darkened host star (Top right).The transit gets progressively deeper towards mid-transit for the latter case, as the planet transits through the brighter region near the pole (Bottom right).The difference in transit depth mid-transit is around 600 ppm.

Figure 3 .
Figure3.Fit of the two gravity darkening models (GD7400,GD7600) to the phase folded transit light curve from TESS.Both models return very similar fits, resulting in the overlapping lines in the middle plot.

Figure 4 .
Figure 4. Violin plot showing the posterior distribution for obliquity and stellar inclination derived from the two gravity darkening models.The analysis favors a flipped orbit, and KELT-19Ab likely has the largest obliquity among currently known exoplanets.The over-plotted histogram shows the obliquity of currently known exoplanets taken from Albrecht et al. (2022) with the count on the bottom x-axis.

Table 1 .
and  of retrograde planets

Table 2 .
List of parameters and their chosen priors

Table 3 .
Planetary and stellar parameters obtained from gravity darkening analysis of TESS data assuming  eff,pole = 7400 K and  eff,pole = 7600 K as well as from existing literature.