This site uses cookies. By continuing to use this site you agree to our use of cookies. To find out more, see our Privacy and Cookies policy.

The following article is Open access

Nodal Precession and Tidal Evolution of Two Hot Jupiters: WASP-33 b and KELT-9 b

, , , , , , and

Published 2022 June 1 © 2022. The Author(s). Published by the American Astronomical Society.
, , Citation Alexander P. Stephan et al 2022 ApJ 931 111 DOI 10.3847/1538-4357/ac6b9a

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/931/2/111

Abstract

Hot Jupiters orbiting rapidly rotating stars on inclined orbits undergo tidally induced nodal precession measurable over several years of observations. The Hot Jupiters WASP-33 b and KELT-9 b are particularly interesting targets because they are among the hottest planets found to date, orbiting relatively massive stars. Here, we analyze archival and new data that span 11 and 5 yr for WASP-33 b and KELT-9 b, respectively, in order to model and improve upon their tidal precession parameters. Our work confirms the nodal precession for WASP-33 b and presents the first clear detection of the precession of KELT-9 b. We determine that WASP-33 and KELT-9 have gravitational quadrupole moments $({6.3}_{-0.8}^{+1.2})\times {10}^{-5}$ and $({3.26}_{-0.80}^{+0.93})\times {10}^{-4}$, respectively. We estimate the planets' precession periods to be ${1460}_{-130}^{+170}$ yr and ${890}_{-140}^{+200}$ yr, respectively, and that they will cease to transit their host stars around the years ${2090}_{-10}^{+17}$ CE and ${2074}_{-10}^{+12}$ CE, respectively. Additionally, we investigate both planets' tidal and orbital evolution, suggesting that a high-eccentricity tidal migration scenario is possible to produce both system architectures and that they will most likely not be engulfed by their hosts before the end of their main-sequence lifetimes.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Over recent decades, the number of discovered exoplanets has steadily increased into the thousands, revealing a large variety of exoplanet system orbital architectures around stars of all evolutionary stages and a large range of stellar masses (e.g., Charpinet et al. 2011; Gettel et al. 2012; Howard et al. 2012; Barnes et al. 2013; Nowak et al. 2013; Niedzielski et al. 2015, 2016; Reffert et al. 2015). Of peculiar interest among these exoplanetary systems have been Hot Jupiters (HJs), which orbit their stars with periods of 10 days or shorter, as such planets do not exist in our own solar system. Additionally, many discovered HJs have strong spin–orbit misalignments with their host stars, often even being in retrograde orbits (e.g., Albrecht et al. 2012; Winn & Fabrycky 2015). This, together with planet formation models, supports the idea that HJs are not born at their current orbital configuration, but have migrated post-formation via some dynamical mechanism (e.g., Nagasawa et al. 2008; Naoz et al. 2011, 2013, 2012).

Regardless of the formation or migration mechanisms, spin–orbit misaligned HJs present a unique observational opportunity to measure tidally induced nodal precession of an HJ's orbit, if the host star is oblate due to rapid rotation. As colder, more convective stars tend to "spin down" rapidly due to magnetic braking as they age (e.g., Kawaler 1988), hotter (Teff ≳ 6200 K), stars with radiative envelopes can stay fast-spinning for most of their main-sequence (MS) lifetime (e.g., Kraft 1967; Meynet et al. 2011). Therefore, misaligned HJs around hot stars are ideal targets for measuring nodal precession.

Here, we present our observations and analysis of the nodal precession and tidal dynamics of two HJs orbiting two hot, rapidly spinning stars, WASP-33 b and KELT-9 b. These exoplanets are among the hottest exoplanets that have been discovered over recent decades, orbiting stars with effective temperatures above 7400 K and 10,000 K, respectively (e.g., Collier Cameron et al. 2010; Gaudi et al. 2017). Both planets have been extensively studied by previous works, investigating among other factors their orbits and atmospheres (e.g., Collier Cameron et al. 2010; von Essen et al. 2014; Johnson et al. 2015; Gaudi et al. 2017; Ahlers et al. 2020; Watanabe et al. 2020; Borsa et al. 2021; Pai Asnodkar et al. 2022a, 2022b), and they are known to have large spin–orbit misalignments. In this work, we obtain new data to extend the timeline of nodal precession measurements by several years for both planets compared to previous works (see Sections 2 and 3.3), allowing us to derive more accurate tidal parameters for both systems (see Section 4). Additionally, we estimate the past tidal evolution of both planets, allowing us to set limits on the formation and migration mechanisms that brought them to their current configurations (see Section 5).

2. Observations

2.1. WASP-33 b

We use archival data for WASP-33 b that were obtained by the Harlan J. Smith Telescope (HJST) with the Robert G. Tull Coude Spectrograph (Tull et al. 1995) at McDonald Observatory on UT 2008 November 12th (Collier Cameron et al. 2010) and UT 2014 October 4th (Johnson et al. 2015).

In addition, we include new data from HJST that were obtained on UT 2016 December 11 and data from the PEPSI (Strassmeier et al. 2015) spectrograph on the Large Binocular Telescope (LBT) on UT 2019 November 17 (Cauley et al. 2020).

There are also archival data from the Subaru telescope. However, we do not have access to the Subaru data that were obtained on UT 2011 October 19 (Watanabe et al. 2020). Because the Subaru observation was bracketed by the two ends of the time baseline from 2008 to 2019, the Subaru data are nonessential in the subsequent analyses. For a similar reason, we do not include archival data from HARPS-N that span from 2016 to 2019 (Borsa et al. 2021).

2.2. KELT-9 b

We use archival data in the discovery paper (Gaudi et al. 2017). These data were obtained by the Tillinghast Reflector Echelle Spectrograph (TRES) on the 1.5 m telescope at the Fred Lawrence Whipple Observatory on UT 2014 November 14, 2015 November 06, and 2016 June 12.

We also include archival data that have been obtained since the discovery. These data include two transits observed with the HARPS-N spectrograph mounted on the Telescopio Nazionale Galileo (TNG) in La Palma, Spain on UT 2017 Aug 01 and 2018 July 21 (Hoeijmakers et al. 2019), and one transit observed with PEPSI at LBT on UT 2018 July 03 (Cauley et al. 2019).

In addition, we add a new data set from PEPSI at LBT that were observed on UT 2019 June 22, therefore extending the time baseline to ∼5 yr from 2014 to 2019. We note that CARMENES (Quirrenbach et al. 2018) has also observed KELT-9 b (Yan et al. 2019). The observation was made on UT 2017 January 5. We do not include this data point because it falls in the time baseline and does not add a significant constraint to the nodal precession measurement.

3. Modeling the Doppler Tomography (DT) Signal

Both WASP-33 b and KELT-9 b have been studied with DT (Johnson et al. 2015; Gaudi et al. 2017; Watanabe et al. 2020). With the new data sets that have been taken in more recent epochs, we perform a uniform analysis to model the DT signal for the two planets. The uniform analysis will eliminate any systematics introduced by using different models.

3.1. Extracting Line Profiles

We follow the least-square deconvolution (LSD; Kochukhov et al. 2010; Pai Asnodkar et al. 2022a) method to extract line profiles (LPs). The LSD method requires a stellar template spectrum, for which we use the PHOENIX BT-Settl model (Allard et al. 2012, 2013, and references therein). We do not interpolate between stellar parameters. Instead, we use the spectra whose effective temperature, surface gravity, and metallicity are the closest to the reported values for WASP-33 b and KELT-9 b, i.e., Teff = 7400 K, log(g) = 4.0, and solar metallicity for WASP-33 b, and Teff = 10,200 K, log(g) = 4.0, and solar metallicity for KELT-9 b.

The LSD method also requires a scalar value for the regularization matrix. The scalar controls the smoothness of the extracted LPs. We experiment with a wide range of possible values that span four orders of magnitude, and we choose a scalar that retains the DT signal and planet absorption signal (if visible) while removing the stochastic noise induced by the inversion process of LSD.

3.2. Removing Stellar Pulsation

WASP-33 is a δ-Scuti variable star that shows a very strong stellar pulsation signal. The stellar pulsation is so strong that it overwhelms the DT signal. As detailed in Cauley et al. (2020), we apply a customized Fourier filter to remove the pulsation signal. Because the pulsation signal and the DT signal are almost orthogonal to each other, the two signals can be easily separated in Fourier space. A similar strategy was adopted in previous works (Johnson et al. 2015; Watanabe et al. 2020).

We do not attempt to remove the stellar pulsation signal for KELT-9 b, although doing so was suggested in Hoeijmakers et al. (2019). This is because the pulsation signal does not pose a challenge for us in the modeling of the DT signal for KELT-9 b.

3.3. Modeling DT

We adopt an analytical framework to model the DT signal as described in Cauley et al. (2020). In short, the DT signal is a Gaussian perturbation of the stellar line profile (Hirano et al. 2011). Therefore, the DT signal can be analytically calculated, and the details of that process are described in Wang et al. (2018). The analytical framework can be applied to modeling not only the DT signal (e.g., WASP-33 b, Cauley et al. 2020) but also (1) the planet absorption signal for high-resolution spectroscopy (as shown for KELT-9 b) and (2) photometric and spectroscopic observations of eclipsing binaries (e.g., EPIC 203868608, Wang et al. 2018).

For WASP-33 b, the modeling parameters for the DT signal are the impact parameter, projected spin–orbit alignment angle, projected rotational velocity, quadratic limb darkening parameters, planet–star radius ratio, systemic velocity, and the mid-transit time. For KELT-9 b, the planet absorption signal during the transit is clearly visible, so we add the following parameters: planet–star mass ratio, planet absorption depth and width, and a net velocity shift in the planet rest frame. Figure 1 and Figure 2 show the modeling results for WASP-33 b and KELT-9 b.

Figure 1.

Figure 1. From top to bottom are our modeling results for WASP-33 b data in 2008 (TULL), 2014 (TULL), 2016 (TULL), and 2019 (PEPSI). Left: Fourier-filtered residual map after subtracting the median line profile. The planet "Doppler shadow" is the diagonal blue track running from the bottom right to the top left. Occasional gaps are due to interpolation onto a time array with a fixed interval. Middle: Modeled DT signal. The overall gradient along the time axis is due to the shift of stellar radial velocity. The enhanced red background during transit is due to the fact that we set the line profile normalization to unity for all line profiles. Right: differences between the maps in the left and middle panels.

Standard image High-resolution image
Figure 2.

Figure 2. From top to bottom are our modeling results for KELT-9 b data in 2014 (TRES), 2015 (TRES), 2016 (TRES), 2017 (HARPS-N), 2018 (HARPS-N), 2018 (PEPSI), and 2019 (PEPSI). Left: Residual map after subtracting the median line profile. The planet "Doppler shadow" is the vertical blue track due to the nearly polar orbit of KELT-9 b. The planet absorption signal is the diagonal red track that runs from bottom left to top right. The planet signal is more apparent in 2017–2019 data taken from HARPS-N and PEPSI, due to a higher signal-to-noise ratio. Occasional gaps are due to interpolation onto a time array with a fixed interval. Middle: Modeled DT and planet absorption signal. Right: differences between the maps in the left and middle panels.

Standard image High-resolution image

3.4. Sampling Posteriors

We used PyMultiNest (Buchner et al. 2014) for posterior sampling in a Bayesian framework. We fit each data set independently and obtain posterior samples for all parameters as described in Section 3.3. In this work, we focus on measuring the nodal precession as manifested by the change of projected orbital obliquity (λ) and impact parameter (b). Priors are uniform distributions between −180° and 180° for λ and between −1 and 1 for b. The likelihood function is $\exp [-{({ \mathcal D }-{ \mathcal M })}^{2}/{{ \mathcal E }}^{2}]$, where ${ \mathcal D }$ is data, ${ \mathcal M }$ is model, and ${ \mathcal E }$ is the error term, which we use as the standard deviation of the difference between the LPs and the median LP over time. The posterior distributions of λ and b are summarized in Table 1.

Table 1. Measured Values for λ and b and Derived Values for ip , Ω, and I for WASP-33 b and KELT-9 b

WASP-33 b
SourceHJSTHJSTHJSTPEPSI
Mid-transit Time [UT]2008.11.12 10:022014.10.04 06:292016.12.11 06:522019.11.17 07:59
λ [°] $-{109.82}_{-0.55}^{+0.54}$ $-{110.27}_{-0.41}^{+0.40}$ $-{110.29}_{-0.60}^{+0.55}$ $-{109.97}_{-0.41}^{+0.41}$
b ${0.1502}_{-0.0075}^{+0.0199}$ ${0.0821}_{-0.0059}^{+0.0090}$ ${0.0548}_{-0.0143}^{+0.0084}$ $-{0.0074}_{-0.0080}^{+0.0077}$
ip [°] ${87.67}_{-0.31}^{+0.12}$ ${88.72}_{-0.14}^{+0.09}$ ${89.15}_{-0.13}^{+0.22}$ ${90.11}_{-0.12}^{+0.12}$
Ω [°] ${87.52}_{-0.33}^{+0.13}$ ${88.64}_{-0.15}^{+0.10}$ ${89.09}_{-0.14}^{+0.23}$ ${90.12}_{-0.13}^{+0.13}$
I [°] ${109.80}_{-0.54}^{+0.55}$ ${110.27}_{-0.40}^{+0.41}$ ${110.28}_{-0.55}^{+0.60}$ ${109.97}_{-0.41}^{+0.41}$
KELT-9 b (First Analysis Procedure)
SourceTRESTRESTRESHARPS-NPEPSIHARPS-NPEPSI
Mid-transit Time [UT]2014.11.14 05:092015.11.06 03:582016.06.12 08:552017.08.01 02:042018.07.03 07:142018.07.21 01:482019.06.22 06:57
λ [°] $-{86.23}_{-0.23}^{+0.26}$ $-{85.06}_{-0.84}^{+0.70}$ $-{87.36}_{-0.34}^{+0.41}$ $-{84.03}_{-0.33}^{+0.31}$ $-{85.67}_{-0.17}^{+0.16}$ $-{84.78}_{-0.38}^{+0.37}$ $-{87.34}_{-0.08}^{+0.07}$
b ${0.1176}_{-0.0037}^{+0.0057}$ ${0.1297}_{-0.0059}^{+0.0085}$ ${0.1766}_{-0.0052}^{+0.0047}$ ${0.1912}_{-0.0033}^{+0.0039}$ ${0.1873}_{-0.0010}^{+0.0009}$ ${0.1954}_{-0.0033}^{+0.0038}$ ${0.2004}_{-0.0045}^{+0.0006}$
ip [°] ${87.86}_{-0.10}^{+0.07}$ ${87.64}_{-0.15}^{+0.11}$ ${86.79}_{-0.09}^{+0.10}$ ${86.52}_{-0.07}^{+0.06}$ ${86.59}_{-0.02}^{+0.02}$ ${86.45}_{-0.07}^{+0.06}$ ${86.36}_{-0.01}^{+0.08}$
Ω[°] ${87.86}_{-0.10}^{+0.07}$ ${87.63}_{-0.16}^{+0.11}$ ${86.78}_{-0.09}^{+0.10}$ ${86.50}_{-0.07}^{+0.06}$ ${86.59}_{-0.02}^{+0.02}$ ${86.43}_{-0.07}^{+0.06}$ ${86.35}_{-0.01}^{+0.08}$
I [°] ${86.23}_{-0.26}^{+0.23}$ ${85.06}_{-0.70}^{+0.84}$ ${87.37}_{-0.41}^{+0.34}$ ${84.04}_{-0.31}^{+0.33}$ ${85.68}_{-0.16}^{+0.17}$ ${84.79}_{-0.37}^{+0.38}$ ${87.34}_{-0.07}^{+0.08}$
KELT-9 b (Second Analysis Procedure)
λ [°] $-{85.20}_{-0.29}^{+0.26}$ $-{85.26}_{-0.72}^{+0.57}$ $-{86.69}_{-0.47}^{+0.43}$ $-{84.58}_{-0.32}^{+0.30}$ $-{85.91}_{-0.09}^{+0.09}$ $-{85.18}_{-0.33}^{+0.33}$ $-{86.53}_{-0.09}^{+0.08}$
b ${0.1171}_{-0.0037}^{+0.0056}$ ${0.1363}_{-0.0054}^{+0.0073}$ ${0.1725}_{-0.0051}^{+0.0042}$ ${0.1905}_{-0.0029}^{+0.0036}$ ${0.1884}_{-0.0014}^{+0.0007}$ ${0.1994}_{-0.0031}^{+0.0031}$ ${0.2016}_{-0.0010}^{+0.0006}$
ip [°] ${87.87}_{-0.10}^{+0.07}$ ${87.52}_{-0.13}^{+0.10}$ ${86.86}_{-0.08}^{+0.09}$ ${86.54}_{-0.07}^{+0.05}$ ${86.57}_{-0.01}^{+0.02}$ ${86.37}_{-0.06}^{+0.06}$ ${86.33}_{-0.01}^{+0.02}$
Ω [°] ${87.86}_{-0.10}^{+0.07}$ ${87.51}_{-0.13}^{+0.10}$ ${86.86}_{-0.08}^{+0.09}$ ${86.52}_{-0.07}^{+0.05}$ ${86.57}_{-0.01}^{+0.02}$ ${86.36}_{-0.06}^{+0.06}$ ${86.33}_{-0.01}^{+0.02}$
I [°] ${85.20}_{-0.26}^{+0.29}$ ${85.27}_{-0.57}^{+0.72}$ ${86.69}_{-0.43}^{+0.47}$ ${84.59}_{-0.30}^{+0.32}$ ${85.92}_{-0.09}^{+0.09}$ ${85.19}_{-0.33}^{+0.33}$ ${86.54}_{-0.08}^{+0.09}$

Note. Given are the 68th percentile confidence intervals for each parameter determined from the posterior distributions discussed in Section 3.4. As mentioned in Section 3.4, we applied two analysis procedures to the KELT-9 b data. As explained in Section 3.5, for our precession model, we inflated the uncertainties of λ and b to minimum values of 0fdg75 and 0.01, respectively, where necessary. This increase also implicitly accounts for systematic errors introduced by stellar oblateness.

Download table as:  ASCIITypeset image

3.5. Estimating Systematics

We eliminate systematics by uniformly analyzing all data sets using the analytical modeling framework as described in previous sections. However, our model has limitations in describing the physical process of a planet occulting the stellar disk. For example, the perturbation of the planet shadow to the LP may not be sufficiently described by our analytical model, which is a convolution of a rotational-broadening kernel and an Gaussian instrument profile. This could explain some of the residuals that we see in the high signal-to-noise ratio (S/N) case, e.g., the lower four panels in Figure 2 for data from HARPS-N and PEPSI. For the low-S/N cases, e.g., 2014 February 16 data from TRES for KELT-9 b (top three panels in Figure 2), we also notice structured residual patterns after subtracting the model from the data.

To estimate the systematics in the nodal precession measurement, we take the following two approaches. First, we have two independent data sets for KELT-9 b that are taken close in time: HARPS-N data on 2018 July 20 and PEPSI data on 2018 July 03. The differences in the results for λ and b between these two data sets can be used to estimate the systematics, under the assumption that the nodal precession is smooth over timescales of a year. The differences of λ and b between the two data sets are 0fdg73 and 0fdg01, respectively. These are at least ∼4 and ∼6 times larger, respectively, than the uncertainties based on the posteriors from the best-fit models for λ and b (See Table 1).

The second approach we take is to slightly change our data analysis procedure and check how much the results vary accordingly. For example, instead of subtracting the out-of-transit median LP, we subtract the overall median LP. We find that changes in the analysis procedure can result in changes of λ and b by ∼2–3 times the value of the reported uncertainty.

Given the above exercises in understanding the systematics in the analyses, we inflate the measurement uncertainties for λ and b to 0fdg75 and 0fdg01 for cases that return values lower than these uncertainties based on posterior distributions, in order to account for unknown sources of systematics.

4. Modeling Planetary Precession

We model the observed changes of the projected obliquity λ and impact parameter b over time, assuming precession of the nodes due to tidal interactions between the planets and their host stars. Our analysis is similar to that discussed in previous works (e.g., Watanabe et al. 2020; Borsa et al. 2021). We translate the observed values of λ and b into values for the ascending node versus the line of sight, Ω, and the inclination of the planetary orbital plane versus the plane of the sky, I, via the equations

Equation (1)

and

Equation (2)

where ip is the angle between the observer's line of sight and the planet's angular momentum. The latter can be derived given the impact parameter b and the ratio of the planetary orbit's semimajor axis to the star's equatorial radius, (a/R*), via

Equation (3)

The values for these parameters for the different observation dates are shown in Table 1. The geometries of these various angles are shown in Figure 3.

Figure 3.

Figure 3. Geometries of the orbital elements and observables as defined in this study. The observer's line of sight is defined as the y-axis, with the sky defined along the xz plane. The stellar spin axis is defined to be on the yz plane (see solid red arrow), forming an angle i* with the line of sight. If the star is viewed equator-on, i* = 90°. The orbital angular momentum vector of the planet (see green solid arrow) forms an angle ip with the line of sight and an angle I with the z-axis (I is also the inclination between the orbital plane and the xy plane). The projection of the orbital angular momentum vector onto the xz plane defines the planet's projected obliquity vs. the z axis, λ, while the projection of the vector onto the xy plane defines the ascending node vs. the line of sight, Ω. See also Figure 4 of Watanabe et al. (2020), on which this graphic is based.

Standard image High-resolution image

The long-term precession of the longitude of the ascending node, Ω, and inclination, I, of any planet with significantly smaller orbital angular momentum than its host star's stellar rotational angular momentum can be described by the equations (e.g., Iorio 2016)

Equation (4)

and

Equation (5)

with i* being the angle of the star's rotation axis to the line of sight of the observer, P being the planet's orbital period, a being the planet's orbital semimajor axis, R* being the stellar equatorial radius, and J2 being the star's quadrupole gravitational moment. The values of most of these parameters are taken from the previous literature and are listed in Table 2.

Table 2. Stellar and Planetary Parameters of WASP-33 and KELT-9 from the Literature Used for the Nodal Precession and Tidal Evolution Models

WASP-33
ParameterValueReference
$V\sin {i}_{* }$ (km s−1) ${86.63}_{-0.32}^{+0.37}$ Johnson et al. (2015)
a/Rs 3.69 ± 0.01Collier Cameron et al. (2010)
Rp /Rs 0.1143 ± 0.0002Collier Cameron et al. (2010)
Rs (R)1.444 ± 0.034Collier Cameron et al. (2010)
Mp (MJup)2.81 ± 0.53von Essen et al. (2014)
Ms (M)1.495 ± 0.031Collier Cameron et al. (2010)
P (days)1.2198675 ± 0.0000011von Essen et al. (2014)
KELT-9
ParameterValueReference
$V\sin {i}_{* }$ (km s−1)111.4 ± 1.3Gaudi et al. (2017)
a (au) ${0.03462}_{-0.00093}^{+0.00110}$ Gaudi et al. (2017)
Rp (RJup) ${1.891}_{-0.053}^{+0.061}$ Gaudi et al. (2017)
Rs (R) ${2.362}_{-0.063}^{+0.075}$ Gaudi et al. (2017)
Mp (MJup)2.88 ± 0.84Gaudi et al. (2017)
Ms (M) ${2.52}_{-0.20}^{+0.25}$ Gaudi et al. (2017)
P (days)1.4811235 ± 0.0000011Gaudi et al. (2017)

Download table as:  ASCIITypeset image

We use Equations (4) and (5) to calculate the evolution of Ω and I over time and fit them to our observed data in Table 1 via maximum likelihood estimation, with i*, J2, Ω0, and I0 as free parameters to be determined by the best-fit solution. To determine the confidence intervals for these parameters, we use the MCMC code emcee (Foreman-Mackey et al. 2013).

4.1. WASP-33 b

We apply the model discussed in the previous section to the WASP-33 b data, as shown in Table 1, with inflated uncertainties as explained in Section 3.5. We use the earliest measured values and uncertainties for Ω and I as priors for Ω0 and I0. From the posterior, we determined 1σ confidence intervals (determined by the 16th, 50th, and 84th percentiles) and best-fit model values of J2, i*, Ω0, and I0, which we list in Table 3 (see Figure 4 for a visualization of the probability distributions).

Figure 4.

Figure 4. Corner plot of MCMC-determined values for i*, ${\mathrm{log}}_{10}({J}_{2})$, Ω0, and I0 of WASP-33 b. The black, dashed lines mark the 1σ confidence interval for the values of each parameter, while the blue lines mark the parameter values of the best-fit model solution. See Table 3 for parameter values.

Standard image High-resolution image

Table 3. Modeling Results for the Precessions of WASP-33 b and KELT-9 b

WASP-33
Parameter68th-percentile Confidence IntervalBest-fit SolutionLiterature Values (References)
J2 $({6.3}_{-0.8}^{+1.2})\times {10}^{-5}$ 5.92 × 10−5 (9.14 ± 0.51) × 10−5 (Watanabe et al. 2020), (6.73 ± 0.22) × 10−5 (Borsa et al. 2021)
i* [°] ${85}_{-19}^{+20}$ 84.63 ${96}_{-14}^{+10}$ (Watanabe et al. 2020), 90.11 ± 0.12 (Borsa et al. 2021)
Ω0 [°] ${87.4}_{-0.2}^{+0.2}$ 87.34 
I0 [°] ${109.97}_{-0.64}^{+0.63}$ 109.91 
Precession period [years] ${1460}_{-130}^{+170}$ 1493∼840 (Watanabe et al. 2020), 1108 ± 19 (Borsa et al. 2021)
Predicted year of last transit [CE] ${2090}_{-10}^{+17}$ 2090.3∼2068 (Borsa et al. 2021)
b0 0.16 ± 0.010.1609 
db/dt [yr−1]−0.0143 ± 0.0013−0.0145 
λ0 [°]−110.0 ± 0.5−109.931 
d λ/dtyr−1]+0.02 ± 0.07 + 0.019 
KELT-9
J2 $({3.26}_{-0.80}^{+0.93})\times {10}^{-4}$ 3.38 × 10−4 5.6 × 10−5 < J2 < 2.5 × 10−4 (Ahlers et al. 2020, , estimated)
i* [°] ${139}_{-7}^{+7}$ 139.35 ${142}_{-7}^{+8}$ (Ahlers et al. 2020)
Ω0 [°] ${87.47}_{-0.19}^{+0.18}$ 87.51 
I0 [°] ${84.83}_{-0.43}^{+0.41}$ 84.78 
Precession period [years] ${890}_{-140}^{+200}$ 870 First detection by this work
Predicted year of last transit [CE] ${2074}_{-10}^{+12}$ 2074.3 
b0 0.138 ± 0.0080.1366 
db/dt [yr−1]−0.0144 ± 0.00250.0148 
λ0 [°]−84.8 ± 0.3−84.784 
d λ/dt [°yr−1]−0.31 ± 0.08−0.315 

Note. Given are the determined values of the stellar quadrupole moment, J2, the angle of the stellar spin axis versus the line of sight, i*, the longitude of the ascending node, Ω0, and the inclination, I0, at the time of the initial observations, the nodal precession period, and the estimated year when transits will cease. We also estimate the linear ephemerides for the observable parameters b and λ, with the starting points being the times of the first observations listed in Table 1 for each system. We give the 68th-percentile confidence interval for each parameter, as well as the best-fit solutions. Where possible, we compare our results to previous results from the literature.

Download table as:  ASCIITypeset image

These values translate to the precession model shown in Figure 5, with a precession period of roughly ${1460}_{-130}^{+170}$ yr for the ascending node. In Figure 5, we show the evolution of the ascending node Ω, inclination of the orbital plane to the sky I, inclination of the planet's orbit to the line of sight ip and impact parameter b, respectively. Note that the planet can only transit in respect to an Earth-bound observer if b is between −1 and 1 (excluding "grazing" transits). Based on our model estimates, WASP-33 b's impact parameter will dip below −1 around the year ${2090}_{-10}^{+17}$ and cease to transit for approximately 500 ± 200 yr (see last frame of Figure 5). Overall, WASP-33 b only transits for about 20% of its precession period.

Figure 5.

Figure 5. Evolution of the ascending node Ω (first row), inclination I (second row), the planet's orbital inclination vs. the observer's line of sight ip (third row), and impact parameter b (fourth row) over time for WASP-33 b. The colored symbols and error bars mark the measured values and uncertainties of the parameters, with light blue triangles showing the HJST data and green dots showing PEPSI data, while the black line shows the best-fit model solution based on the data. The red lines show 2000 example model solutions from our MCMC with model likelihoods within the 1σ confidence level of the best-fit model. The left frames focus on the time span of the measurements, while the right frames extrapolate the model for the next 2000 yr. The precession period for Ω is approximately ${1460}_{-130}^{+170}$ yr. The value of b must be between −1 and 1 for observable transits to occur (see black dotted lines). Based on our model estimate, WASP-33 b will cease to transit around the year ${2090}_{-10}^{+17}$ CE.

Standard image High-resolution image

4.2. KELT-9 b

Similar to our analysis of the WASP-33 system discussed above, we apply our model to the KELT-9 b data from the 1st analysis procedure, as shown in Table 1, with inflated uncertainties as explained in Section 3.5. We use additional orbital constraints from Ahlers et al. (2020) in our modeling analysis. They observed KELT-9 b's transit light curve to be clearly asymmetric, from which they were able to constrain the value of i* to be ${i}_{* ,\mathrm{Ahlers}}={142}_{-7}^{+8}$ deg. 5 We thus performed our model fit using i*,Ahlers as the prior. From the posterior, we determined 1σ confidence intervals (determined by the 16th, 50th, and 84th percentiles) and best-fit model values of J2, i*, Ω0, and I0, which we list in Table 3 (see Figure 6 for a visualization of the probability distributions).

Figure 6.

Figure 6. Corner plot of MCMC-determined values for i*, ${\mathrm{log}}_{10}({J}_{2})$, Ω0, and I0 of KELT-9 b, using the prior on i* derived from Ahlers et al. (2020) for the fitting model. The black, dashed lines mark the 1σ confidence interval for the values of each parameter, while the blue lines mark the parameter values of the best-fit model solution. See Table 3 for parameter values.

Standard image High-resolution image

Figure 7 shows the precession evolution based on these fit values, with a precession period of roughly ${890}_{-140}^{+200}$ yr for the ascending node. Like WASP-33 b, KELT-9 b will eventually precess such that it will no longer transit its star from the perspective of an Earth-bound observer, which we estimate will occur around the year ${2074}_{-10}^{+12}$ CE. KELT-9 b will reappear as a transiting planet approximately 300 ± 100 yr afterward.

Figure 7.

Figure 7. Evolution of the ascending node Ω (first row), inclination I (second row), the planet's orbital inclination vs. the observer's line of sight ip (third row), and impact parameter b (fourth row) over time for KELT-9 b, using the prior on i* derived from Ahlers et al. (2020) for the fitting model. Dark blue downward-pointing triangles and error bars show TRES data, yellow squares and error bars show HARPS-N data, and green dots and error bars show PEPSI data. The black line shows the best-fit model, with red lines showing 2000 example models with model likelihoods within the 1σ confidence interval around the best-fit model. The precession period for Ω is approximately ${890}_{-140}^{+200}$ yr.

Standard image High-resolution image

We also fit the data with no priors. The resulting posterior distributions are completely consistent with those shown in Figure 6. The model fit using no priors is provided in the Appendix A for comparison.

5. Discussion

5.1. Improvements on Previous Results and New Findings

One of the goals of this study was to improve upon the established results for the orbital parameters of WASP-33 b and KELT-9 b published in previous works (e.g., Johnson et al. 2015; Gaudi et al. 2017; Watanabe et al. 2020), given the longer baseline for observations and making use of the measured nodal precession of the planets' orbits. As such, we confirm for WASP-33 b that i* is close to 90°, and we determine a best-fit value for J2 that is significantly lower compared to previous works with shorter baselines (e.g., Watanabe et al. 2020) and consistent with recent work with similar extended baseline coverage (Borsa et al. 2021) (see also Table 3).

For KELT-9 b, our work presents the first clear detection of the precession of KELT-9 b. We determine a very large J2 value for KELT-9. We caution that the measurement of i* was complicated due to the variety of observational instruments used for observations. By using the prior based on previous work by Ahlers et al. (2020), we confirm that its best-fit value must be close to 140°. Furthermore, The measured precession rate of approximately −0fdg3 per year is consistent with previous predictions by Ahlers et al. (2020).

A crucial contribution of our study to the understanding of WASP-33 b and KELT-9 b is in the exploration of both planets' tidal evolution, which we explore more in the following section.

5.2. Tidal and Orbital Evolution

Our observations have clearly shown the tidal influence on the precession of the ascending nodes of the orbits for both WASP-33 b and KELT-9 b, with best-fit values for J2, the stellar gravitational quadrupole moments, of 5.92 × 10−5 and 3.38 × 10−4, respectively. These values are in line with expected J2 values for fast-spinning stars, which create their quadrupole moments via their oblateness as a response to their rapid rotation (e.g., Kraft 1967; Ward et al. 1976). The interplay of rotation and internal mass redistribution is a highly complex phenomenon that depends on factors such as the size of a star's radiative and convective regions, differential rotation rates, magnetic field strengths, and many more. Based on our best-fit parameter results derived from the orbital precession, and defining the stellar oblateness f via the flattening formula (Murray & Dermott 2000),

Equation (6)

with ΩS being the stellar angular spin frequency, we estimate the stellar oblatenesses of WASP-33 and KELT-9 to be approximately 0.0185 ± 0.0025 and 0.07 ± 0.01, respectively. Our estimate for KELT-9 is consistent with the previous estimate by Ahlers et al. (2020) of 0.087 ± 0.017. We note here that the systematic errors introduced in b due to the measured oblatenesses is consistent with the minimum adopted error value of 0.01, discussed in Section 3.5.

In our further analysis, we are mostly interested in the tidal interactions between the stars and their planets, in particular in terms of the planets' past and future dynamical evolution. Our data clearly show very strong spin–orbit misalignments. This puts limits on the strengths of the tidal dissipation within the stars that naturally would work toward realignment.

To estimate these tidal strengths, we assume the equilibrium tide model (following equations by Hut (1980), Eggleton et al. (1998), and Eggleton & Kiseleva-Eggleton (2001); see Appendix B). In this model, an orbiting body raises a tidal bulge on its companion that either lags or leads the bulge-raising object by an angle that is determined by the ratio of the spin period to the orbital period. In such a situation, dissipative forces then act both to synchronize the spin and orbital periods as well as to realign the orbital inclination and spin axis. The nature of these dissipative forces, however, are not well-understood for many objects; in particular, the equation of state and internal structure can majorly impact the strengths of such forces. Previous works have shown that a strong distinction exists between mostly radiative versus mostly convective stars (for an overview, see Ogilvie 2014), with radiative stars generally having much weaker dissipative forces.

We test both the orbits of WASP-33 b and KELT-9 b using the equilibrium tide model. We determine that only the weaker dissipative forces for mostly radiative stars can allow their planets to exist at their current orbital configurations for an extended period of time. This result is not surprising, as both WASP-33 and KELT-9 are rather large and hot stars, and should both have radiative envelopes. Additionally, we estimate that the tidal dissipation forces acting within the planets vastly out-scale (by at least five orders of magnitude during peak migration) the forces within the stars. As such, we determine that (1) the planets' spin periods must be synchronized to their orbital periods, (2) the planet's spin and orbital axes are aligned, and (3) the planetary orbits have no residual eccentricities.

Moreover, we can estimate the timescale of migration. We can assume that the planets were originally formed further away from their hosts and migrated inward via high-eccentricity migration, which can be caused by a scattering event or secular mechanisms such as the Kozai–Lidov effect by currently unseen companions (Kozai 1962; Lidov 1962; Naoz 2016). This scenario is consistent with the observed large spin–orbit misalignment.

Starting with original orbital semimajor axes between 10 and 20 au (assuming the formation of the gas giants to occur at or close to the water ice lines within the protoplanetary disks of the respective systems), migration and circularization at their currently observed orbital separations would take on the order of 100 to 200 Myr for WASP-33 b and KELT-9 b, assuming the equilibrium tide model. 6 These timescales are not unrealistic, given the estimated ages of roughly 100 Myr for WASP-33 and 300 Myr for KELT-9 (e.g., Collier Cameron et al. 2010; Gaudi et al. 2017). However, this would imply that migration began very soon after planet formation. An important caveat here is that, if tidal migration indeed began shortly after planet formation, the host stars were likely still in their pre-main-sequence stages, during which their internal convection and tidal energy dissipation rates were significantly increased (e.g., Zanazzi & Wu 2021). Therefore, the migration times we mention above can be seen as conservative upper limits.

These migration scenarios would require a significant dissipation of energy by the planets as heat, reaching on the order of Lp ∼ 1030 erg s−1 during the fastest migration episodes, or ∼0.1% of L, the solar energy output. This energy dissipation is well within realistic levels and would not heat the planets to more than about 2500 K during peak migration. This calculation ignores additional heat sources and assumes blackbody radiation. As the current dayside temperatures of both planets caused by their host stars' incoming irradiation is higher than that, the tidal energy would not significantly alter the planets' chemistry or internal structure. However, if the planets still have any residual internal energy from migration, it could show up as higher-than-expected nightside temperatures. If such a discrepancy exists, it would be a telltale sign of migration and allow us to estimate how recently migration must have occurred, with the caveat that winds and other energy transport mechanisms in the planetary atmospheres could dilute such a signal. High nightside temperatures on the order of 2500 to 3000 K have indeed been observed for exoplanets TOI-1431 b/MASCARA-5 b (Addison et al. 2021) and KELT-9 b (Wong et al. 2020), though in both cases current models explain these temperatures with atmospheric energy transport.

6. Summary

Multi-epoch observations of transiting planetary systems can (1) predict the long-term dynamical evolution and (2) infer stellar and orbital properties that are responsible for the tidal interactions within such systems. We use multi-epoch observations to study two ultra-hot Jupiter systems: WASP-33 b and KELT-9 b. The data span 11 and 5 yr, respectively. Using the Doppler Tomography technique, we adopt an analytical framework to model the high-resolution spectroscopic data sets during the transit. Through the framework, we infer the posterior distribution of the projected obliquity λ and impact parameter b over time. The time-varying λ and b data reveal that both WASP-33 b and KELT-9 b are undergoing orbital precession.

By applying a tidal interaction model to the time-varying λ and b data, we can constrain the J2 values for the two planet host stars (see Table 3) and predict the dynamical fate of the two systems. We determine the nodal precession periods for the orbits of WASP-33 b and KELT-9 b to be ${1460}_{-130}^{+170}$ and ${890}_{-140}^{+200}$ yr, respectively. As such, WASP-33 and KELT-9 b are expected to cease transiting their stars around ${2090}_{-10}^{+17}$ CE and ${2074}_{-10}^{+12}$ CE, respectively. We furthermore find that either planet can only exist at its current orbit for an extended period of time given weak stellar tidal dissipation consistent with stars with radiative envelopes. The weak tidal dissipation also allows for the observed high stellar spin–orbit misalignments. As such, most tidal dissipation within these systems occurs within the planets, implying planetary spin–orbit alignment and synchronization. If the planets reached their current orbital configurations due to high-eccentricity migration, such migration would have occurred on timescales of 100 to 200 Myr. We estimate that the planets can exist at approximately their current orbital configurations for the remaining main-sequence lifetimes of their host stars, to be eventually engulfed during post-MS inflation.

We thank the anonymous referee for providing comments and suggestions. A.P.S. acknowledges partial support from a President's Postdoctoral Scholarship from the Ohio State University and the Ohio Eminent Scholar Endowment. A.P.S. and B.S.G. acknowledge partial support by the Thomas Jefferson Chair Endowment for Discovery and Space Exploration. J.W. acknowledges the academic gift from Two Sigma Investments, LP, which partially support this research. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflects the views of Two Sigma Investments, LP. This paper includes data taken at The McDonald Observatory of The University of Texas at Austin, as well as data obtained through the LBT and HARPS-N. The LBT is an international collaboration among institutions in the United States, Italy and Germany. LBT Corporation partners are: The Ohio State University, representing OSU, University of Notre Dame, University of Minnesota, and University of Virginia; The University of Arizona on behalf of the Arizona Board of Regents; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society, The Leibniz Institute for Astrophysics Potsdam, and Heidelberg University. The HARPS-N project was funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origin of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), and the Italian National Astrophysical Institute (INAF), University of St. Andrews, Queen's University Belfast, and University of Edinburgh.

Appendix A: Precession Model Fit for KELT-9 b without Priors

As mentioned in Section 4.2, our observations for KELT-9 b most likely contained some unknown uncertainties that made fitting the model to the data without additional priors difficult. Here, we provide our model fit to the data using no priors, as a comparison to our cited results in Section 4.2 using the Ahlers et al. (2020) prior. Using no priors for the MCMC results in 1σ confidence intervals for the values of the tidal parameter, ${J}_{2}=({3}_{-0.8}^{+1.8})\times {10}^{-4}$, of the stellar spin axis angle versus the line-of-sight, ${i}_{* }={104}_{-34}^{+30}$ deg, of the initial longitude of the ascending node, ${{\rm{\Omega }}}_{0}={87.50}_{-0.19}^{+0.19}$ deg, and of the initial inclination, ${I}_{0}={85.53}_{-0.64}^{+0.57}$ deg (see Figure 8). These values correspond to the precession evolution shown in Figure 9. However, the large uncertainty in the value of i* highlights the limited information content of the data. In particular, the scatter in the observed inclination (I) values (see Figure 9, second row, left frame), on the same order as the uncertainties, makes finding an appropriate model fit challenging.

Figure 8.

Figure 8. Corner plot of MCMC-determined values for i*, ${\mathrm{log}}_{10}({J}_{2})$, Ω0, and I0 of KELT-9 b, using no priors for the fitting model. The black, dashed lines mark the 1σ confidence interval for the values of each parameter, while the blue lines mark the parameter values of the best-fit model solution. See text of Appendix A for parameter values.

Standard image High-resolution image
Figure 9.

Figure 9. Evolution of the ascending node Ω (first row), inclination I (second row), the planet's orbital inclination vs. the observer's line of sight ip (third row), and impact parameter b (fourth row) over time for KELT-9 b, using no priors for the fitting model. Dark blue downward-facing triangles and error bars show TRES data, yellow squares and error bars show HARPS-N data, and green dots and error bars show PEPSI data. The black line shows the best-fit model, with red lines showing 2000 example models with model likelihoods within the 1σ confidence interval around the best-fit model.

Standard image High-resolution image

Appendix B: Basic Equations of the Equilibrium Tide Model

As discussed in Section 5.2, we use the equilibrium tide equations, following Hut (1980), Eggleton et al. (1998), and Eggleton & Kiseleva-Eggleton (2001), to estimate the tidal circularization and migration timescales for WASP-33b and KELT-9b, assuming high-eccentricity migration with original orbital semimajor axes of tens of au. The most important equations in this context describe the evolution of the inner orbit eccentricity vector, e 1 , and specific angular momentum vector, J 1, which form the ($\hat{{\boldsymbol{q}}}$, $\hat{{{\boldsymbol{J}}}_{1}}$, $\hat{{{\boldsymbol{e}}}_{1}}$) reference frame with the unit vector $\hat{{\boldsymbol{q}}}$ perpendicular to both, via

Equation (B1)

Equation (B2)

with the V, W, X, Y, and Z terms following the functional forms

Equation (B3)

Equation (B4)

Equation (B5)

Equation (B6)

Equation (B7)

Here, n is the inner orbital mean motion, μ is the reduced mass of the inner star and planet companion (individual masses m1 and m2, respectively), R1 and R2 are the inner star and planet radii, a1 and e1 are the inner orbit semimajor axis and eccentricity, k1 and k2 are the respective classical apsidal motion constants (about 0.014 for an n = 3 polytrope and 0.25 for an n = 1 polytrope; see, for example, Kiseleva et al. 1998) of the inner star and the planet, G is the gravitational constant, and Ωs1 and Ωs2 are the spin vectors of the inner star and the planet in the ($\hat{{\boldsymbol{q}}}$, $\hat{{{\boldsymbol{J}}}_{1}}$, $\hat{{{\boldsymbol{e}}}_{1}}$) frame. The complete set of terms can easily be derived by swapping the subscripts for the body-dependent factors, as well as for the tidal friction timescale, tF1. This timescale is described by the equation

Equation (B8)

where tV1 is a viscous timescale on the order of a year (see also Eggleton & Kiseleva-Eggleton 2001).

By numerically integrating Equations (B1) and (B2), we estimate the migration and circularization timescales of WASP-33b and KELT-9b, assuming initial orbital semimajor axes between 10 and 20 au, to be on the order of 100 to 200 Myr, with some uncertainties due to possible variations in initial planetary spin rates and orientations.

Footnotes

  • 5  

    Ahlers et al. (2020) defined i* in their work as 0° when the star is viewed equator-on, while in our work i* is defined as 0° when the "northern" pole is directly pointing away from the observer, as shown in Figure 3, which is simply an offset of 90°. We thus translated their published measurement of ${i}_{* }={52}_{-7}^{+8}$ deg in accordance with our definition in this paper.

  • 6  

    We note here that these timescales depend on the assumed value of the viscous timescale tV , which is usually estimated to be on the order of years or decades for stars and gas giants (see Eggleton & Kiseleva-Eggleton 2001), as well as the exact nature of the eccentricity-inducing mechanism.

Please wait… references are loading.
10.3847/1538-4357/ac6b9a