The Ambiguous Age and Tidal History for the Ultra-Hot Jupiter TOI-1937Ab

To constrain the TTVs for this system, we performed a one parameter fit for each individual transit to find its center time of transit ($t_{\mathrm{tra}}$) using the Monte Carlo Markov Chain code emcee in conjunction with batman (Foreman-Mackey et al., 2013). The planet parameters (except for $t_{0}$) were set to the best-fit values reported in Yee et al. (2023) and reproduced in Table 1. Each fit was run with 8 walkers for 5000 steps with a uniform prior on $t_{0}$ to within 0.4 days of the predicted time of transit from the Yee et al. (2023) solution. Once the $t_{0}$ values were found, we fit a quadratic model to the values and the remaining residuals are the TTVs. Our fitted best-fit TTV values and uncertainties for each epoch for which data exists are plotted in Figure 1. Our maximum variation in transit time was 6.63 minutes, which occurred during transits that took place during TESS sectors 7 and 9. Due to the 1800 second cadence in these sectors, these transits often had only one or two data points during the transit.

In the full TTV curve, we found no evidence of periodic perturbations in the TTV’s that would indicate the presence of another planet. To evaluate the potential for orbital decay, we fit a quadratic model using a emcee fit with 32 walkers and 10000 steps to our derived center times of transit (as done in Harre et al., 2023):

when $N$ denotes the epoch of the transit event, $t_{0}$ the time of transit at epoch $N=0$, $P$ the orbital period of TOI-1937Ab, and $dP/dN$ the orbital decay rate. The period derivative $\dot{P}=dP/dt$ is computed by $dP/dt=(1/P)\ dP/dN$. Using the posterior generated from this fit, we found a 3$\sigma$ lower limit on the change in orbital period $dP/dt$ of -0.09 sec/yr. In Figure 2, we show 150 random draws from the posterior of the quadratic model overlaid with the TTV data, along with a line corresponding to the 3$\sigma$ limit. In Figure 3, we present a histogram of the posterior of the decay rate $dP/dt$. Further TESS observations of TOI-1937A may extend the baseline and improve this limit, but we find no current evidence of significant decay in TOI-1937Ab’s orbit.

The age of TOI-1937A is currently debated within the literature, with some claiming that the system is a tidal tail member of the 150 Myr NGC 2516 open cluster, while others claim that its membership is ambiguous (Kounkel & Covey, 2019; Yee et al., 2023). However, the low upper limit implied by the Li I $6708\AA$ equivalent width (EW) ($<25\,m\AA$ at 3$\sigma$; Yee et al., 2023) indicates a significantly older age than that of NGC 2516. We investigated the age constraints imposed by this lithium EW upper limit using BAFFLES¹¹1https://github.com/adamstanfordmoore/BAFFLES — a software package capable of computing age posteriors for field stars from measurements of a target’s $B-V$ color and lithium EW (Stanford-Moore et al., 2020). Incorporating $B-V=0.76$ (Henden et al., 2016) and the star’s aforementioned EW upper limit, we calculate a 3$\sigma$ minimum age of 501 Myr (a 1$\sigma$ minimum age of 4.6 Gyr), which is inconsistent with the measured age of the cluster. These findings are illustrated in Figure 4.

A third way to estimate the age of this system is gyrochronology. In order to estimate the photometric rotational period of TOI-1937A, we used a method similar to that laid out in Capistrant et al. (2024). We performed a Lomb-Scargle (LS) periodogram analysis (Lomb, 1976; Scargle, 1982) on the full TESS lightcurve, assisted by the Phase Dispersion Minimization (PDM) algorithm (Stellingwerf, 1978) and an autocorrelation function (ACF) (McQuillan et al., 2013; Shumway & Stoffer, 2005).

All three period-search techniques were implemented on the full light curve, comprised of all seven sectors of TESS data along with the light curves from ground-based facilities. We searched for periodic signatures within a range of 0.04–30 days. We then performed a by eye search of the resulting power spectrum, autocorrelation plot, phase folded light curve, and periodogram to search for periodic structure to determine the rotation period. We determined the stellar rotation period to be 6.5 days, in agreement with the 6.6-day estimate in Yee et al. (2023).

After computing our estimate of the stellar rotation period, we used gyrointerp²²2https://github.com/lgbouma/gyro-interp — a framework that calculates the posterior probability for a star’s age given its observed rotation period and effective temperature (Bouma et al., 2023). We assumed no planet is present and used the stellar temperature found in Yee et al. (2023). We found an age of $510.94^{+1217.51}_{-261.51}$ Myr to $3\sigma$ confidence (see Figure 5).

Both age estimates, whether obtained through lithium abundance methods or through gyrochronological methods, are at odds (to greater than $3\sigma$) with the 150 Myr reported age of NGC 2516. Additionally, while the $3\sigma$ minimum age estimate for the lithium abundance does allow for a 510 Myr age estimate, it is near the minimum of the 3$\sigma$ age range. All three age estimates differ wildly from each other, with the 1$\sigma$ minimum age of 4.6 Gyr and the 510 Myr gyrochronological age estimate both indicating that TOI-1937A is a field star, and thus not belonging to NGC 2516. Though the rotation period of the star could reasonably fit within the upper edge of the spread of the rotation sequence of NGC 2516 found in Bouma et al. (2021), the potential for tidal spin up (discussed further in Section 3) leads to gyrochronological methods being an unreliable measure of the star’s age on its own.

POET provides a self-consistent framework for modeling the evolution of planetary orbits and stellar rotation, taking into account the effects of tidal interactions and stellar wind. It also incorporates the impact of core-envelope differential rotation within the stellar interior. In this section, we will outline the main equations POET uses to compute the orbital evolution of the planet and the rotational evolution of its host star, which underpin our numerical explorations of the tidally-altered stellar rotation rate of TOI-1937A.

First, the orbital evolution of the planet’s orbit due to the (static) equilibrium tide can be written as a pair of coupled differential equations (Goldreich, 1963; Barnes et al., 2008),

and

Here, $Q^{\prime}=3Q/(2k_{2})$ is the modified tidal dissipation factor, where $Q$ represents the tidal dissipation efficiency and $k_{2}$ the Love number, $m_{p}$ and $r_{p}$ denote the planetary mass and radius, respectively, $G$ is the gravitational constant, $M_{*}$ and $R_{*}$ are the stellar mass and radius, $a$ is the planetary semi-major axis, $e$ is its eccentricity, $n$ denotes the planetary mean motion, $\omega_{\text{surf}}$ the angular velocity of the stellar surface, and $\sigma=\text{sign}(2\omega_{\text{surf}}-3n)$ attains the value 1 when the planet’s orbital frequency is lower than the stellar rotation frequency and -1 in the opposite case.

To model stellar rotation evolution over time, we must write an equation for how the stellar angular momentum changes over a star’s lifetime due to two factors: (a) in response to the planet’s tidally-induced orbital evolution, and (b) due to angular momentum loss caused by stellar winds (Kawaler, 1988). Following the formulation and notation of Penev et al. (2014), we can write an expression for the change in angular momentum due to these two effects as follows:

where $K$ is a parameter describing the strength of the stellar magnetic wind, $\omega_{\text{sat}}$ the star’s angular velocity at which the magnetic wind saturates. Solving Equation 4 simultaneously with Equations 2 and 3 yields the evolution over time for the planetary orbital ($a$ and $e$) and the stellar rotation rate (set by $L_{*}$).

POET also allows, for low-mass stars like TOI-1937A, the rotation rates of the stellar core and envelope to be decoupled. POET computes the orbital evolution of the planet simultaneously with the evolution of the stellar angular momenta and moments of inertia for the convective and radiative zones inside the star and the stellar mass and radius. These evolutions are based on pre-computed MESA tracks. Finally, POET also models disk-locked stellar rotation while the protoplanetary disk is present, and the transition between non-tidally locked and tidally locked planet rotation. The details of how these calculations are implemented are given in Penev et al. (2014). POET has been used successfully in the past to assess how stellar parameters such as $Q_{*}$ affect the orbital evolution of planets (Penev et al., 2018), evaluate tidal dissipation in stellar binaries (Penev & Schussler, 2022), and derive parameter constraints for $Q_{p}$ and $Q_{*}$ (Mahmud et al., 2023; Patel et al., 2023). For the work in the remainder of this section, we use the publicly available version of POET³³3https://github.com/kpenev/poet.

In attempting to run a tidal evolution simulation that represents the true evolution of TOI-1937A and TOI-1937Ab, we must choose initial orbital parameters for the planet and initial stellar parameters (such as its initial angular momentum $L_{*}$ and tidal quality factor $Q_{*}$) to initialize the simulation. None of these parameters are exactly known, as all we can measure from the observational data is current-day values. Because of this, choosing our initial simulation parameters is complicated by the fact that our age estimate of TOI-1937A is uncertain, with different methods giving very different age estimates. (ranging between 150 Myr and 4.5 Gyr or more). This means that the inferred $L_{*}$ needed to match the current-day observed rotation rate will be different for the different stellar age estimates.

In this subsection, our goal is to determine how much the tidal migration of a planet like TOI-1937Ab could have been expected to alter the rotation rate of its host star. To make our exploration computationally feasible, we choose one of our computed possibilities of the stellar age to use for our simulations ($\sim 500$ Myr) and choose $L_{*}$ such that the star has its observed rotation rate ($P_{rot}$ = 6.5 days) at this age.

With the stellar initial parameters chosen, we used POET to compute the evolution for a set of 300 combinations of initial semi-major axis and eccentricity of the planet, spanning $10<a/R_{*}<100$ and $0<e<1$. For all simulations, stellar mass and radius and the planet mass were set to the best-fit values given in Table 1, and the stellar quality factor set to an assumed value of $Q_{*}=10^{5}$. The stellar companion TOI-1937B was not modeled in the simulations, as its large projected separation of 1030 AU found in Yee et al. (2023) suggests it is dynamically decoupled. The simulations for TOI-1937A begun with disk-locked rotation with a stellar rotation period of 5 days, and the disk set to dissipate at 4 Myr. At that point, the planet is assumed to form and the stellar rotation is allowed to evolve under the influence of the stellar wind and planetary tidal evolution. For each simulation, we computed the stellar rotation period at 500 Myr.

In addition to the 300 simulations evaluating the effect of the companion planet, we also ran one control simulation under the same conditions outlined above, but with no planet. This integration shows the star’s rotational evolution in the absence of planet-induced tidal evolution, yielding a rotation rate of 6.5 days after 500 Myr. However, for each individual simulation, the behavior of the planet affects the observed rotation rates. Two examples of this effect can be seen in Figure 6. In both panels, the control simulation is given as a dashed line, and the red line shows the evolution of the stellar rotation as computed in the presence of the planet. A second line (right y-axis, tan line) shows the semi-major axis evolution of the planet. These examples were chosen as representative of each of their behaviors, as described in the caption on the figure. Two two panels show ways in which the planet can increase or decrease the stellar rotation rate. In the upper panel, the planet’s orbit expands as it extracted angular momentum from the star and the stellar rotation period slows compared to the no-planet integration. In the bottom panel, the planet experiences quick tidal in-spiral immediately after the disk dissipates and it is subsequently engulfed by the star. This results in it depositing its angular momentum onto the star and increasing the stellar rotation rate as compared to the no-planet case.

The results of the full suite of 300 simulations are shown in Figure 7. Our simulated planets followed a variety of tracks in their orbital evolution, including being engulfed by their host. We note that due to the coarseness of our sampled grid, none of our simulations had planets that ended the simulation at a semi-major axis matching that of that of TOI-1937Ab ($a/R_{*}\approx 3.8)$. As such, our simulations do not recreate the exactly geometry of TOI-1937A, but rather compute the range of alterations to stellar rotational period that could be possible due to interactions with a planet with TOI-1937Ab’s parameters. For simulated systems where the planet avoided being engulfed by the star, the stellar rotation period ranged between a minimum of 5 - 5.75 days (for high-$e$, low-$a$ initial conditions) and a maximum of 7 days (for low-$e$, low-$a$ initial conditions). Low-$e$ and high-$a$ planet parameters resulted in no alteration to the stellar rotation period as compared to the 6.5 days seen in the no-planet case.

Our initial, limited set of simulations does not encompass the entire range of possible initial parameters for this system. We have used only one value of initial spin angular momentum $L_{*}$ for the star and simplified the problem by assuming constant tidal quality factors for the star and the planet. However, for this particular set of simulations, it is evident that the effect of the planet-induced tidal evolution over the first 500 Myr of evolution is to change the inferred stellar rotation by (at maximum) a couple of days. Depending on the exact initial parameters, this change can either spin up the star, resulting in a faster rotation period, or, for a small range of initial parameters, it can slow the stellar rotation. Given the observed orbit of this planet, it is unlikely that it has slowed the stellar rotation int the past, as it would have had to in the past have resided on an even closer orbit as compared to its current orbit. It is more likely that the effect of the planet has been to increase the star’s rotation compared to what it would have been if no planet were present.

While our simulations in the previous subsection demonstrated the range of possible effects of tidal evolution due to a massive planet like TOI-1937Ab, our simulations were not sufficiently focused to answer the question: what is the future tidal evolution of TOI-1937Ab expected to look like? In order to probe the potential dynamical future of TOI-1937Ab, we ran seven simulations initialized to the current state of TOI-1937A, including all the best-fit parameters reported in Table 1 ($a/R_{*}=3.8$, $e=0$, $R_{p}=1.247R_{J}$, $M_{p}=2.0M_{J}$, $M_{*}=1.072M_{\odot}$, and a stellar rotation rate of 6.5 days at 500 Myr, as assumed in the previous section). The simulation was initialized at 500 Myr with those parameters, and the only value changed between simulations was the stellar tidal quality factor $Q_{*}$. $Q_{*}$ is a parameter expected to vary significantly across physical and orbital (Patel et al., 2023) regimes. It cannot be measured directly, and must instead be inferred from dynamical constraints, such as planet survival (Hamer & Schlaufman, 2020).

We tested values of $Q_{*}$ logarithmically spaced between $10^{5}$ and $10^{11}$ and ran integrations for 4.5 Gyr (from a starting time of 500 Myr to a end time of 5 Gyr). The results of these simulations are plotted in Figure 8. All seven integrations show that TOI-1937Ab is currently in a state of tidal decay. However, the different values of $Q_{*}$ indicate different likely futures for this planet: for $Q_{*}\leq 10^{7}$, the planet is expected to be engulfed by the star within the next 500 Myr; for $10^{8}\approx Q_{*}$, the planet’s orbit will decay and it will be engulfed in slightly over 1 Gyr from now; for $10^{9}\leq Q_{*}$, the planet will not be engulfed within 5 Gyr.

Next, we computed the instantaneous decay rate of the orbital period, with a goal of comparing how these theoretically-derived decay rates compared to the upper limit on the decay rate computed in Section 2. The results of these calculations are plotted in Figure 9. However, none of our simulated planets had decay rates that were greater than our derived upper limit on the decay rate for TOI-1937Ab. For this reason, a concrete constraint on the $Q_{*}$ of TOI-1937A cannot be made with currently available data.

Our simulations provide partial evidence suggesting that this planet attained its observed orbit recently, given that, unless the host star’s $Q_{*}$ is anomalously high, it would have been engulfed by its host star relatively quickly. However, this cannot be taken as evidence on the host star’s age, as it is possible that the dynamical interaction that placed TOI-1937Ab on its high-periastron orbit (which allowed for the start of the tidal evolution process) happened recently, even if the star is older. Additionally, our POET simulations to asses the rotational period alterations caused by the planet’s tidal evolution suggest that the change in rotational period for the host star is marginal (or the order of 1-2 days, which is commensurate with the $3\sigma$ errors on the rotation rate) unless the planet is engulfed by its host star.

From our TTV analysis, we find no evidence of a non-linear ephemeris, and find a $3\sigma$ limit on the orbital period decay rate of -0.09 sec/yr. Additionally, our TTV analysis shows no evidence of additional planets near TOI-1937Ab. Since it is expected that this planet migrated via tidal migration—a process that excites higher orbital eccentricity and could lead this planet to cross the orbits of any other companion planets (Dawson & Johnson, 2018)—it is not surprising that we find no evidence of nearby companions. The majority of hot Jupiters are observed not to have close planetary companions. However, the absence of detected companions does not necessarily mean they are not present. Due to the large measured mass of this planet, any nearby planetary companions would have to be relatively large to cause observable TTVs.

We also considered the tidal evolution of TOI-1937Ab. TOI-1937Ab is a hot Jupiter in an ultra-short period orbit, one of only roughly eleven similar planets in the exoplanet census. Most observed USP planets are less than 2 $R_{\odot}$ (Winn et al., 2018), and population-wide evidence suggests that hot Jupiters are subject to significant tidally-driven orbital decay on the main sequence (Hamer & Schlaufman, 2019). While the age of the systems remains amibiguous, if TOI-1937A were a young ($<1$ Gyr) host star, it would place a useful constraint on the earlier arrival of hot Jupiters to their orbits. We also note that our POET simulations, presented in Section 3, demonstrate that it is likely that TOI-1937Ab is in actively decaying orbit, with its final time of engulfment set by the stellar $Q_{*}$. For smaller ($\approx 2.6M_{Earth}$) USP planets, Hamer & Schlaufman (2020) finds that stars need to have $Q_{*}>10^{7}$ for the planets to be stable against tidal inspiral as a population. For TOI-1937Ab, a $Q_{*}\approx 10^{10-11}$, a value far above the median $Q_{*}$ for smaller USP planets found in Hamer & Schlaufman (2020), would protect the planet against engulfment on the main sequence.

The presence of the stellar companion TOI-1937 B at a projected separation of approximately 1030 AU (Yee et al., 2023) raises the possibility of dynamical interactions that could influence the orbit of TOI-1937Ab. One potentially relevant interaction is the Kozai-Lidov mechanism (Lidov, 1962; Kozai, 1962), which can induce large oscillations in a planet’s orbital eccentricity and inclination, potentially setting the stage its subsequent tidal migration into its final observed short-period orbit. While the full orbital parameters of the stellar companion (TOI-1937B) - including its mass, eccentricity, and semi-major axis - are not well constrained, we performed an analysis to evaluate whether the Kozai-Lidov mechanism could have contributed to the formation of TOI-1937Ab’s current orbit.

As the orbital and physical parameters of companion TOI-1937 B are not known, for this analysis we assumed a companion mass of 0.3 $M_{\odot}$, approximated based on the brightness difference between the primary and secondary and consistent with the results derived in Christian et al. (2024), and assumed that its measured projected separation of 1030 AU is its semi-major axis. The Kozai-Lidov timescale can be computed by (see Equation 42 of Antognini 2015; see also Holman et al. 1997; Christian et al. 2022):

where $m$ denotes masses, $P$ orbital periods, $e$ orbital eccentricity, and subscripts $1,2,3$ denote the primary star, planet, and stellar companion, respectively. We computed the Kozai-Lidov timescale over a range of plausible initial orbital periods and eccentricities for TOI-1937Ab, prior to any scattering or migration event. The results are presented in Figure 10, which shows Kozai-Lidov timescales for different initial conditions of the planet’s orbit. To contextualize these timescales, we overlaid contour lines corresponding to the estimated age of the system.

Only orbital configurations to the right of these lines have sufficiently short Kozai-Lidov timescales to operate within the system’s lifetime. In Figure 10, we over-plot grey points corresponding to the cold Jupiters ($M_{p}>0.5M_{jup}$) discovered via radial velocity observations.⁴⁴4Data obtained from the IPAC Exoplanet Archive, 2/17/2025. These points populate the parameter space where the Kozai-Lidov timescale would be short enough to be expected to operate for all feasible system ages. As a result, we consider it dynamically feasible that if TOI-1937Ab had started as a typical cold Jupiter from the RV sample (i.e., Bryan et al., 2016; Wittenmyer et al., 2020; Rosenthal et al., 2021), TOI-1937B could have played a role in setting TOI-1937Ab on an initial eccentric orbit.

However, we note that for the mechanism to result in the observed hot Jupiter orbit, there are additional conditions beyond the Kozai timescale being short enough to operate within the stellar lifetime: the planet’s periastron distance would need to become sufficiently small to allow for significant tidal dissipation and circularization within the system’s lifetime. Additionally, a Rossiter-McLaughlin measurement by Yee et al. (2023) found that the orbit of the hot Jupiter TOI-1937Ab is roughly aligned with the spin axis of the host star TOI-1937A. This would suggest that TOI-1937Ab either migrated via coplanar high-eccentricity migration (Petrovich, 2015) or had its obliquity damped post-migration (Lai, 2012; Spalding & Winn, 2022). While these additional factors require further investigation, our analysis supports the plausibility of Kozai-Lidov cycles as a contributing mechanism in the evolutionary history of TOI-1937Ab.

The age estimates computed in this work for TOI-1937A are presented in Table 2. Different pieces of dynamical and physical evidence give age estimates. The gyrochronological age estimate is around 510 Myr, the Li abundance gives a $1\sigma$ lower limit of roughly 4.5 Gyr, and were TOI-1937A a member of the NCG 2516 cluster it would be roughly 150 Myr old. The non-detection of lithium, described in Yee et al. (2023), remains the most confounding piece in age dating this star. While this star is likely not part of the NGC 2516 open cluster, as previously suggested in Kounkel & Covey (2019), much of the evidence, with the exception of lithium abundance, points to this star being younger than 1 Gyr. In Sevilla et al. (2022), it is suggested that a planetary engulfment could decrease surface lithium abundance in more massive stars ($1.3-1.4$$M_{\odot}$) as a result of an enhancement of internal mixing and diffusion processes, however, TOI-1937A is not in this mass regime. Thus, we believe that the lithium-derived age estimate is inconclusive and further investigation and observations are required to confidently assign an age for this system.