Giant exoplanets orbiting close to their host stars are unlikely to have formed in their present configurations 1 . These 'hot Jupiter' planets are instead thought to have migrated inward from beyond the ice line and several viable migration channels have been proposed, including eccentricity excitation through angular-momentum exchange with a third body followed by tidally driven orbital circularization 2,3 . The discovery of the extremely eccentric (e = 0.93) giant exoplanet HD 80606 b (ref. 4) provided observational evidence that hot Jupiters may have formed through this high-eccentricity tidal-migration pathway 5 . However, no similar hot-Jupiter progenitors have been found and simulations predict that one factor affecting the efficacy of this mechanism is exoplanet mass, as low-mass planets are more likely to be tidally disrupted during periastron passage [6][7][8] . Here we present spectroscopic and photometric observations of TIC 241249530 b, a high-mass, transiting warm Jupiter with an extreme orbital eccentricity of e = 0.94. The orbit of TIC 241249530 b is consistent with a history of eccentricity oscillations and a future tidal circularization trajectory. Our analysis of the mass and eccentricity distributions of the transitingwarm-Jupiter population further reveals a correlation between high mass and high eccentricity.

The Transiting Exoplanet Survey Satellite (TESS) 9 monitored the apparent brightness of the star TIC 241249530 for 28 days during the second year of its primary mission. These data reveal a transit-like approximately 0.8% dip in brightness on 12 January 2020, the shape and depth of which were consistent with a Jupiter-sized planet passing in front of the star (Fig. 1a). To find out the nature and origin of this signal, we conducted a series of ground-based observations of TIC 241249530. We first used high-spatial-resolution speckle imaging data from NESSI 10 to rule out the presence of contaminating sources and confirm that the signal was not associated with a background eclipsing-binary in the TESS aperture. We then began radial velocity (RV) observations with the NEID spectrograph 11 , which revealed that the TESS transit was probably induced by a giant exoplanet on a highly eccentric (e = 0.94), long-period (P = 167 days) orbit. These measurements were consistent with the absence of a transit detection when TESS re-observed this star for 27 days from December 2022 to January 2023. Further NEID measurements, supplemented by observations with the HPF 12 and HARPS-N 13 spectrographs, were strategically scheduled to be taken when the planet was predicted to be approaching periastron and thus inducing large stellar RV variations. We attempted to detect a second transit using the global Unistellar telescope network 14 in March 2023, but these efforts were unsuccessful as the ephemeris was not yet well constrained. However, RV data collected during the periastron window enabled us to more precisely predict the subsequent transit window. We captured the first half of this transit using the engineered diffuser 15 on the ARCTIC imager 16 on 30 August 2023 (Fig. 1b). We refined the orbit using the ARCTIC data together with further NEID observations, including several concurrent with this transit and the subsequent one on 12 February 2024. Our ensemble of photometric and RV measurements is best explained by a massive exoplanet on a long-period, eccentric orbit.

To characterize the host star TIC 241249530, we separately analysed the NEID and HARPS-N spectra using synthetic spectral fitting techniques and we then fit the spectral energy distribution (SED; see Methods). TIC 241249530 is a main-sequence star that is slightly hotter, larger and more massive than the Sun; the derived parameters, listed in Extended Data Table 1, suggest that the star is 3.2 ± 0.5 Gyr old. The star also has a low-mass binary stellar companion, TIC 241249532, at a projected separation of 4.930 ± 0.104″, or 1,664 ± 11 .

We jointly fit the NEID, HPF, HARPS-N, TESS and ARCTIC measurements, accounting for perturbations to the in-transit RV signal owing to the Rossiter-McLaughlin effect. The transit and RV fits are shown in Figs. 1 and 2, respectively, and the best-fit parameters are given in Table 1. TIC 241249530 b is an exoplanet that is 4.98 -0.18 +0.16 times as massive as Jupiter and it is on a 165.77190 -day -0.00028 +0.00027 orbit around its host star, with an eccentricity of 0.9412 -0.0009 +0.0009

. Our fit to the Rossiter-McLaughlin signal (Fig. 2b) shows that the exoplanet is orbiting in the opposite direction to the projected stellar spin λ ( = 163.5 ) -7.7

+9.4 ∘ and is retrograde to 99.5% confidence. Few exoplanets have orbits as extreme as this; this orbit is more eccentric than that of any other transiting exoplanet, and only a handful of known planets have similarly large projected spin-orbit misalignments 17 .

The planet that most closely resembles TIC 241249530 b is HD 80606 b (ref. 4), which has a mass 4.1 times that of Jupiter and is also on a misaligned orbit with a period of 111 days and an eccentricity of 0.93. HD 80606 b is an archetypal example of an exoplanet destined to become a hot Jupiter with an eventual orbital period of less than 10 days. The eccentric orbit carries the planet close enough to its host star at periastron that tides raised on the planet and star will sap energy from the orbit, causing it to gradually shrink and circularize. Also, simulations 3,5 show that the present orbit of HD 80606 b is consistent with a history of von Zeipel-Lidov-Kozai (vZLK) eccentricity oscillations [18][19][20] driven by angular momentum exchange with HD 80607, the stellar companion to the host star. Our own simulations of the dynamical history and trajectory of TIC 241249530 b (see Methods) show that the orbit is consistent with this same type of perturber-coupled, high-eccentricity tidal migration. Eccentricity oscillations would have continued until the most recent few hundred million years, at which point general relativistic precession overtook the torque exerted by the companion, locking cadence and binned 15-min cadence measurements, along with the residuals to the fitted transit signal. d, Residuals to the best-fit model for the ARCTIC transit detection. All brightnesses are given in parts per thousand (ppt). Error bars on individual data points indicate the 1σ measurement uncertainties.

Table 1 | TIC 241249530 b system parameters Parameter Value Description T 0 2458860.8007 0.0016 0.0015 -+ Time of mid-transit (BJD) P 165.77190 0.00028 0.00027 -+ Orbital period (days) e -+ 0.9412 0.0009 0.0009 Orbital eccentricity ω -+ 42.32 0.36 0.40 Argument of periastron (°) i 85.17 0.51 0.57 -+ Orbital inclination (°) K 463.3 4.0 4.1 -+ RV semi-amplitude (m s -1 ) M p 4.98 0.18 0.16 -+ Exoplanet mass (M J ) R p 1.19 0.04 0.04 -+ Exoplanet radius (R J ) λ 163.5 7.7 9.4 -+ Projected spin3orbit obliquity (°) M ⋆ -+ 1.271 0.068 0.061 Stellar mass (M ⊙ ) R ⋆ -+ 1.397 0.028 0.025 Stellar radius (R ⊙ ) vsini ⋆ 4.60 0.63 0.56 -+ Projected rotational velocity (km s -1 )

We report the median values of the posterior distributions from our joint fit to the observed transits and RVs. The uncertainties represent the 68% confidence intervals (±1σ) for each parameter.

the exoplanet on an eccentric orbit that is now gradually circularizing. The architectures of the HD 80606 and TIC 241249530 systems lend support to this process as a plausible hot-Jupiter-formation mechanism. However, although other giant exoplanets on tidal migration tracks have been discovered 21 , including two that probably have vZLK-driven dynamical histories 22,23 , no previous examples have eccentricities >0.9 and none have formation scenarios as clear as that of HD 80606 b.

The observed occurrence rate of super-eccentric progenitors to hot Jupiters 24,25 falls well short of predictions from simulations 26 , suggesting that giant-planet migration is dominated by other channels. With the discovery of TIC 241249530 b, a second super-eccentric exoplanet in a hierarchical triple system has been added to the sparse sample, providing a new lens through which we can explore the formation of these planets. Not only do the TIC 241249530 b and HD 80606 b systems share similar orbital architectures but these exoplanets also have similar masses. The masses and eccentricities of all transiting warm Jupiters, which we define as giant planets with intermediate periods (10 days < P < 365 days), are shown in Fig. 3b. These two planets, which are the only members of the sample with super-eccentric orbits (e > 0.9), are also among the most massive. A correlation between exoplanet mass and eccentricity has been identified in previous works 27330 , each of which found that higher-mass planets are more likely to have larger orbital eccentricities. We find that our narrower sample of transiting giant planets conforms to this known trend (see Methods); the eccentricity distributions of high-mass (M p > 2 M J ) and low-mass (0.3 M J ≤ M p ≤ 2 M J ) members of this population are statistically distinct (Fig. 3a). Although lower-mass planets are more likely to be found on low-eccentricity orbits, high-mass planets exhibit a broad, nearly flat distribution from circular to highly eccentric orbits.

Although the observed mass3eccentricity correlation may be shaped by several processes, such as collisional eccentricity growth 30 or resonant interactions with the protoplanetary disk 31333 , the high masses of TIC 241249530 b and HD 80606 b may offer a clue as to the dearth of super-eccentric giant planets. During the high-eccentricity phase of vZLK oscillations, orbital eccentricities can be driven so close to unity that the exoplanet will approach, or even breach, the tidal radius of the host star. Because the tidal radius is inversely proportional to the planet3star mass ratio, lower-mass planets more easily cross this threshold and experience tidal disruption. A relative dearth of low-mass, eccentric progenitors to hot Jupiters is a consistent outcome of simulations of high-eccentricity tidal migration following vZLK oscillations under an equilibrium tide assumption 638 . For planets susceptible to chaotic, or diffusive, dynamical tidal evolution, whereby oscillations excited in the planet accelerate orbital decay, this mass dependence is largely erased, as low-mass planets can become decoupled from the perturber before being disrupted 34 . However, chaotic tides facilitate circularization on a much shorter timescale (<100 Myr); these planets will spend very little time with intermediate-period orbits 34336 . It is possible that only the most massive eccentric giant planets last long enough in this period regime to be represented in the observed sample.

TIC 241249530 b passes through periastron just six hours before each transit, presenting a unique opportunity to observe how an exoplanet atmosphere responds to a rapid, extreme heating event. Temporal variations in exoplanet atmospheres are best explored through studies of planets on eccentric orbits, for which we may see signatures of time-varying irradiation and changing pressure3temperature profiles, such as turbulent surface flows 37 and disequilibrium chemistry 38 , depending on the heat-redistribution timescales. The atmospheres of several eccentric giant planets have been studied 39341 , but the periastron phase has not been captured in transit for these systems. The orbital geometry of TIC 241249530 b will make such measurements possible for the first time. The planetary atmosphere can also be studied by means of emission measurements during other orbital phases, but the orientation precludes a secondary eclipse at a 6.3σ confidence.

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-024-07688-3.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article9s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article9s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2024 1 U.S. National Science Foundation National Optical-Infrared Astronomy Research Laboratory (NSF NOIRLab), Tucson, AZ, USA. 2 Department of Astronomy and Astrophysics, The Pennsylvania State University, University Park, PA, USA. 3 Center for Exoplanets and Habitable Worlds, The Pennsylvania State University, University Park, PA, USA. 4 Department of Physics,

TIC 241249530 was observed with TESS 9 from 24 December 2019 to 21 January 2020 (Sector 20) at a 30-min cadence and from 21 December 2022 to 18 January 2023 (Sector 60) at a 2-min cadence. A single transit-like dip (flux depth about 8 parts per thousand) was identified by the TESS Single Transit Planet Candidate Working Group (TSTPC WG) in the Sector 20 Quick Look Pipeline 42,43 light curve using a box least-squares search. The TSTPC WG focuses on searching full-frame TESS light curves for isolated transit events and validating and confirming those that are true planets, with the aim of increasing the yield of TESS planets with period >30 days (for example, refs. 44347). There is no flux centroid motion during the transit event for TIC 241249530 and we identify no other sources brighter than Δm G = 5 in the target aperture. Although there is flux contamination from two nearby stars with 6 > Δm G > 5, TIC 241249532 and TIC 241249533, both of which were centred on the same pixel as TIC 241249530 in Sector 20, these are too faint to have been responsible for the observed change in brightness. No notable brightness variations were detected in the Sector 60 light curve. For all subsequent analysis in this work, we rely on the pre-search data conditioned simple aperture photometry 48350 (PDCSAP) light curve from the Science Processing Operations Center 51 (SPOC) for Sector 60 and the TESS-SPOC 52 light curve for Sector 20 (Extended Data Fig. 1).

To verify that the transit signature detected by TESS was indeed associated with TIC 241249530 and not with a nearby star or binary system that was blended in the TESS aperture, we used the NN-EXPLORE Exoplanet Stellar Speckle Imager (NESSI) 10 on the WIYN 3.5-m telescope at Kitt Peak National Observatory to conduct high-spatial-resolution observations of the target on 21 April 2021. A sequence of 1,000 40-ms exposures was taken in the 832-nm and 562-nm narrow-band filters simultaneously with the red and blue NESSI cameras, respectively. These diffraction-limited exposures were used to reconstruct high-contrast images (Extended Data Fig. 2) following the steps outlined in ref. 53. The achieved 5σ contrast limits are sufficient to rule out the presence of faint stellar companions and background sources with Δmag 562 < 3.3 and Δmag 832 < 3.7 at a separation of 0.2″ and Δmag 562 < 3.9 and Δmag 832 < 4.8 at a separation of 1″.

We used the Unistellar Network, a collaboration of citizen scientists using Unistellar telescopes 54 in support of astronomical research, to observe TIC 241249530 from locations in Japan, Europe and the United States in search for transit signatures in March 2023. Observations were taken at various times from 7 to 19 March 2023, when the companion orbital period and transit ephemeris were still highly uncertain. After removing off-target and saturated frames, we calibrated the remaining images, binned them in sets of 15330 to amplify the signal-to-noise ratio (S/N) and performed differential photometry 55,56 . No signatures of statistical significance were found in the Unistellar data and, based on our subsequent orbit fit, we confirm that none of these observations were taken during the transit.

We observed TIC 241249530 again on 30 August 2023 with the Astrophysical Research Consortium Telescope Imaging Camera (ARCTIC) 16 on the ARC 3.5-m telescope at Apache Point Observatory (APO). Observations were conducted using a beam-shaping diffuser, which creates a stable top-hat point spread function of the star to improve photometric precision 15 . We used the Semrock narrow-band filter (8383876 nm) to avoid atmospheric absorption bands 57 . We began observing when the target rose above an air mass of 4 (altitude ≈ 10°) and continued until 12° morning twilight, collecting a continuous 4.3-h baseline of consecutive 30-s exposures. As the star rose above air mass approximately 1.5, about 2.5 h after the start of the observing sequence, a transit-like decrease in brightness was observed.

To reduce the ARCTIC data, a median-combined master bias image was constructed and subtracted from the individual science frames, which were flat-fielded using dome flat exposures taken at the start of the night. We performed differential aperture photometry on the reduced data using AstroImageJ 58 with a 17-pixel (7.7″) aperture and four reference stars that were carefully selected to minimize the scatter of the out-of-transit flux. Flux uncertainties were calculated following the procedures in refs. 15,59, which account for photon noise from the star and background, detector read noise and air-mass-dependent scintillation noise. We removed exposures flagged by AstroImageJ for approaching the detector saturation limit, as well as exposures taken during intermittent cloud cover that introduced further scatter.

The diffused point spread function of TIC 241249530 overlapped with that of TIC 241249532. Before initiating our ARCTIC observing sequence, we collected several individual exposures without the diffuser in the optical path. We used these data to calculate the relative brightness contributions of the two stars. TIC 241249532 contributes just 0.53% of the total flux in the Semrock bandpass.

We monitored the RV signal of TIC 241249530 with the NEID spectrograph 11 on the WIYN 3.5-m telescope, collecting measurements on 40 separate nights between 2 September 2021 and 1 March 2024. On all but three of these nights, single exposures were taken, with exposure times ranging from 500 to 1,800 s, depending on the observing conditions. On the night of 30 August 2023, four consecutive 20-min exposures were taken simultaneously with the partial transit as observed with ARC-TIC, and on the subsequent night, we secured a pair of measurements separated by an hour. We also obtained a sequence of 15 consecutive 20-min exposures on the night of 12 February 2024; this sequence covered a full transit as well as several measurements before ingress and after egress. We discard two spectra that were taken on nights for which the wavelength calibration was identified to be unreliable, leaving us with 56 high-quality measurements with a median S/N per extracted pixel of 25 at 550 nm. The raw echelle spectra were processed with version 1.3 of the NEID Data Reduction Pipeline (DRP; https://neid.ipac. caltech.edu/docs/NEID-DRP/), which produces wavelength-calibrated 1D spectra and then calculates RVs using the cross-correlation function (CCF) method 60 . We also independently calculated the RVs from the calibrated 1D spectra using a modified version of the SpEctrum Radial Velocity AnaLyser (SERVAL) template-matching algorithm 61,62 that has been optimized for NEID spectra as described by ref. 63. The SERVAL RVs were calculated using the central 7,000 pixels of 79 orders centred between 4,010 and 8,400 Å (order indices 20 to 100, corresponding to echelle orders 153 to 73). The template-matching results outperform the CCF-based RVs from the NEID DRP, with median single measurement precisions of σ RV,SERVAL = 6.3 m s -1 and σ RV,DRP = 7.9 m s -1 , so we chose to use the SERVAL RVs for the analysis performed in this work.

Further RV measurements were taken with the Habitable-zone Planet Finder (HPF) spectrograph 12 , which is on the Hobby3Eberly Telescope (HET) 64,65 at McDonald Observatory, and the HARPS-N spectrograph, mounted on Telescopio Nazionale Galileo (TNG) in La Palma, as TIC 241249530 approached periastron in March 2023. Six HPF observations were made between 6 and 31 March 2023, for which each observation consisted of two consecutive 945-s exposures with a median nightly binned S/N per extracted pixel of 137 at 1,000 nm. These data were processed using the HxRGproc 66 and barycorrpy 67 packages and the RVs were calculated using a version of SERVAL that has been modified for HPF 68,69 . We achieve a median RV measurement precision of 15.0 m s -1 . We also observed the target five times with HARPS-N between 7 and 18 March 2023, with an exposure time of 3,300 s and a mean (min, max) S/N of 55 (37, 75). We reduced the data with the offline version of the HARPS-N data-reduction software through the Yabi web interface 70 installed at the Italian Center for Astronomical Archives Data Center. To extract the RVs, we used a G2 mask template and obtained a CCF width of 9.9 km s -1 , with an average precision of 0.1 km s -1 . The median resulting RV measurement precision is 3.4 m s -1 . We show the complete RV time series from NEID, HPF and HARPS-N in Extended Data Fig. 3.

To determine the stellar atmospheric parameters of TIC 241249530, we analysed the out-of-transit NEID spectra collected before September 2023 (cumulative S/N ≈ 100 at 550 nm) using the iSpec 71,72 Python package to perform synthetic spectral fitting. We used the SPECTRUM radiative transfer code 73 , MARCS atmospheric models 74 , solar abundances from 3D hydrodynamic models 75 and the sixth version of the Gaia ESO survey (GES) atomic line list 76 . The microturbulence velocity was treated as a free parameter to allow for flexibility in accounting for small-scale motions in the stellar atmosphere. Macroturbulence was determined using an empirical relation, making use of established correlations with other stellar properties 77 . To streamline the fitting, we restricted the analysis to specific spectral regions from 480 to 680 nm, encompassing the wing segments of the Hα, Hβ and Mg I triplet lines, which are sensitive to T eff and logg, and the Fe I and Fe II lines, which provide precise constraints on [Fe/H] and vsini ⋆ . We minimize the difference between the synthetic and input spectra by applying the nonlinear least-squares Levenberg3Marquardt fitting algorithm, using constraints from the aforementioned models and line lists.

The HARPS-N spectra were independently analysed with BACCHUS 78 , using MARCS atmospheric models, the GES atomic line list and the TUR-BOSPECTRUM radiative transfer code 79,80 . For our fit, we constrained T eff by requiring Fe I line abundances to be uncorrelated with their respective excitation potentials in the synthetic spectrum and we constrained logg by requiring ionization balance for the Fe I and Fe II lines. We also required the Fe I line abundances to be uncorrelated with their equivalent widths and the stellar metallicity ([Fe/H]) was calculated as the average of these abundances. The projected rotational velocity was estimated by fitting the broadening of the Fe I lines, accounting for the best-fit microturbulence and assuming the same macroturbulence contribution as in the iSpec analysis. The stellar parameters derived from the NEID and HARPS-N spectra are largely in good agreement (<1σ). Discrepancies between the [Fe/H] values and vsini ⋆ values at the 1.2σ level probably result from differences between the fitted microturbulence, which is known to exhibit small variations for different fitting methods 72 . We adopt the iSpec T eff , logg, [Fe/H] and vsini ⋆ for the rest of the analysis in this work.

We performed an analysis of the broadband SED of TIC 241249530 together with the Gaia DR3 parallax following the procedures described in refs. 81383. We use JHK S magnitudes from 2MASS 84 , W13W3 magnitudes from WISE 85 , G BP G RP magnitudes from Gaia 86 , BVgri magnitudes from APASS 87 and the NUV magnitude from GALEX 88 . We also used the Gaia spectrophotometry spanning 0.431.0 µm. Altogether, the available photometry spans the full stellar SED over the wavelength range 0.2310.0 µm. We fit the SED using PHOENIX stellar atmosphere models 89 , with the effective temperature, surface gravity and metallicity set to the spectroscopically determined values. The remaining free parameter is the extinction (A V ), which we limited to the maximum line-of-sight value of A V = 0.44 mag from galactic dust maps 90 . The resulting fit is shown in Extended Data Fig. 4. Integrating the unreddened model SED yields the bolometric flux at Earth, F bol = 7.19 ± 0. 20 × 10 -10 erg s -1 cm -2 . Taking the F bol and T eff together with the Gaia parallax, we calculate the stellar radius to be R ⋆ = 1.404 ± 0.028 R ⊙ . Also, the stellar mass is inferred using empirical relations 91 , giving M ⋆ = 1.24 ± 0.07 M ⊙ , and we estimate the age to be 3.2 ± 0.5 Gyr by fitting the evolutionary state with the Yonsei3Yale isochrone models 92 . Our reported 0.5-Gyr uncertainty accounts for the uncertainties on each of the inputs to the isochrone fit: effective temperature, surface gravity, metallicity and stellar mass. However, this does not account for systematic uncertainties arising from our choice of stellar models, which can be on the order of 1 Gyr.

The best-fit extinction for our SED model is A V = 0.31 ± 0.02. This large value is supported by a clear detection of interstellar absorption in the Na D doublet and the K I 770 nm lines in the NEID spectra. Both spectroscopic analyses yield T eff values that are substantially hotter than the literature value from Gaia DR3 spectrophotometric analysis 92 , which is consistent with the effect of reddening from dust along the line of sight to the star.

Using the projected rotational velocity and stellar radius, we place an upper limit on the rotation period of 16.9 days -2.6 +3.8

. We attempt to make a more precise measurement of the rotation period to determine the stellar inclination, but the existing data are insufficient. An analysis of the TESS light curves using the TESS Systematics-Insensitive Periodogram package 93 shows no notable photometric modulation on timescales shorter than the length of each individual sector. We also examine archival photometry of the star from the WASP survey 94 . These data consist of 2,178 measurements on 33 nights, with two isolated epochs in April 2006 and March 2007, and the remaining data covering October 2007 to March 2008. In spite of the substantially longer baseline than TESS, a Lomb3Scargle periodogram analysis of the WASP measurements reveals no notable peaks besides the half-day, one-day and two-day sampling aliases. The lack of photometric modulation is reflected in the spectroscopic data as well; we do not detect periodic variation in the activity-sensitive Ca II H & K, Na I or Hα spectral lines as measured by the NEID DRP. Also, there is no emission in the Ca II H & K line cores in the NEID and HARPS-N spectra, suggesting that the star is chromospherically quiet.

TIC 241249530 shares a common parallax and proper motion with TIC 241249532 as measured by Gaia, and the two stars are separated on the sky by 4.930 ± 0.104″ (ref. 86). The probability of a chance alignment between TIC 241249530 and TIC 241249532 is R = 9.73 × 10 -5 (ref. 95), suggesting that the pair is indeed gravitationally bound. Gaia9s photometric measurements of TIC 241249532 place it firmly along the main sequence. We do not perform an independent SED analysis on this star but instead estimate its mass using empirical mass3 luminosity relations 96,97 . We calculate the mass to be 0.453 ± 0.012 times that of the Sun based on the 2MASS K s -band magnitude and 0.400 ± 0.016 times that of the Sun based on the Gaia G RP -band magnitude. The stellar mass, coordinates and broadband photometry are given in Extended Data Table 1.

On the basis of the weighted mean of the Gaia parallax measurements for the system, the on-sky separation corresponds to a projected physical separation of 1,664.00 ± 10.85 au. The relative motions of these two stars are not constrained well enough by Gaia to meaningfully estimate an orbital solution. However, as an effort to quantify the dynamical impact of TIC 241249532 in our analysis, we simulated 10 million orbits sampled randomly in phase, uniformly in cosi and thermally (f(e) = 2e) in eccentricity. We determine the orbital period of the system to be >10,000 years, with a peak in frequency at 35,000 years. Long-period stellar companions such as this can directly bias RV analyses of exoplanets in the form of a linear RV slope. However, our simulations show that TIC 241249532 probably induces a linear trend in the observed RVs of TIC 241249530 at the level of just 1 cm s -1 year -1 , with 99% of our orbits returning slopes <30 cm s -1 year -1 . The amplitude of this signal is small compared with the km-s -1 -level variations induced by TIC 241249530 b and we therefore do not include it as an extra body in our joint fitting.

We use the exoplanet software package 98 to fit a transit model and a Keplerian orbit with Rossiter3McLaughlin perturbations to the observed photometric and RV signals for TIC 241249530. The exoplanet package relies on starry 99,100 and the underlying analytic models from ref. 101 to fit the transits, and the orbital parameter posteriors are sampled using the PyMC3 Hamiltonian Monte Carlo package 102 . The orbit model consists of a full Keplerian with tight Gaussian priors on the orbital period, P, and time of conjunction, T 0 , broad uniform priors on the exoplanet mass, M p , transit impact parameter, b, and transit depth, δ, and Gaussian priors on the stellar mass, M ⋆ , and radius, R ⋆ . We reparameterize the eccentricity, e, and argument of periastron, ω, as e ω sin and e ω cos and we sample these on the unit disk. We do not impose an extra eccentricity prior, as the global warm-Jupiter eccentricity distribution is not well constrained 103 . Separate quadratic limb-darkening coefficients, reparameterized as q 1 and q 2 as in ref. 104, are used for each instrument for which in-transit observations were taken (TESS, ARCTIC, NEID). For the RV data, we fit individual zero-point offset terms (σ) and jitter terms (γ) for each instrument, splitting the NEID RVs into two separate datasets before and after the instrument restart. Dilution terms are included for both TESS and ARCTIC, as both transit measurements suffered from flux contamination. The TESS data products already account for dilution, but previous works have demonstrated that these results are susceptible to overcorrection 105,106 , so we allow the TESS dilution term to float uniformly from 0.1 to 1.2. For ARCTIC, we fix the dilution to be 0.9947 based on the out-of-transit data for which the target was well resolved from its companion. To model the Rossiter3McLaughlin signal, we adopt the formalism of ref. 107 along with their prior distributions for the Gaussian line dispersion parameter, β, the Lorentzian line dispersion parameter, γ, and macroturbulence, ζ. We place a Gaussian prior on the projected stellar rotational velocity, vsini ⋆ , and a uniform prior on the projected spin3 orbit misalignment, λ. The prior distributions and posterior results for all of these fit parameters, as well as for some derived values, are given in Extended Data Table 2.

The stellar obliquity, ψ, is related to the projected obliquity, λ, by ψ i i λ i i cos = sin sin cos + cos cos .

Here i ⋆ is the inclination of the stellar spin axis and i is the inclination of the exoplanet orbit. We cannot directly calculate ψ because the stellar inclination is not known. Instead, assuming that the stellar inclination is drawn from an isotropic distribution, uniform in cosi ⋆ , we use the above equation to determine the possible values of ψ and their relative probabilities. For our derived posteriors on λ and i, we find that the orbit is indeed retrograde (that is, ψ > 90°) at 99.5% confidence in this scenario, and we calculate the obliquity to be ∘ ψ = 141 -24 +15 . This value is consistent with expectations for vZLK-driven migration; simulations show that the final obliquity can be as large as 180° for systems such as this 7 . This is not definitive proof of the formation history, however, as retrograde orbits such as this can also be produced through planet3 planet interactions 108,109 . Regardless, we warn that our result is strongly dependent on the naive assumption of an isotropic stellar inclination distribution, which is not always valid 110,111 .

The high eccentricity and tight orbit of TIC 241249530 b and the presence of the distant stellar binary companion indicate a likely history of high-eccentricity migration driven by vZLK oscillations and tidal dissipation. To determine how this formation channel could have delivered the exoplanet to its current orbit, we first identify a set of initial conditions consistent with the present-day architecture of the system. We work in the context of the secular approximation for the evolution of hierarchical triple configurations 112 .

The planet is at present close enough to the primary star such that short-range forces4general relativity (GR), tides and rotational distortions4have quenched any vZLK oscillations driven by the companion. We calculate the semimajor axis at which this quenching occurs by assuming that GR dominates the short-range forces and examining the ratio between the timescale of GR precession of the inner orbit and the timescale for vZLK oscillations. To leading order (quadrupole limit), this ratio is 113

⋆ t t a a e M c e G M M = 3 (1 -) (1 -) ( + ) GR quad 4 2 3 2 2 2 2 3/2 p 2

in which M ⋆ is the mass of the primary star; M p , a and e are the mass, semimajor axis and eccentricity of the planet; and M 2 , a 2 and e 2 are the mass, semimajor axis and eccentricity of the binary companion. For the exoplanet and primary-star parameters, we adopt the median posteriors from our joint fit. For the binary companion, we set M 2 = 0.453 M ⊙ , a 2 = 1,664 au and e 2 = 0.5. If the planet started with a low initial eccentricity of e = 0.1, vZLK oscillations would have started only if the initial semimajor axis of the planet was a i > 4.2 au, for which this value is calculated by setting t GR /t quad = 1.

We can now constrain the initial eccentricity by requiring that the periastron distance of the first vZLK oscillation was sufficiently small to trigger efficient tidal dissipation. In particular, in the quadrupole limit, the quantity a a e ≡ (1 -)

is approximately conserved throughout the tidal migration, as the orbital angular momentum is conserved both during episodes of maximum eccentricity, as well as after vZLK oscillations have been quenched. Here e max indicates the maximum eccentricity reached during a vZLK oscillation and a f is equal to the final semimajor axis once the orbit has fully circularized. If a f is taken to be conserved, the maximum eccentricity of the initial vZLK oscillation must have been e i,max > 0.9947.

Exciting an eccentricity this high on the initial vZLK oscillation must have required a substantial initial inclination, I i , between the orbit of the planet and that of the binary companion. We derive a lower bound on I i using the following equation from ref.

112:               ε j e j j I 1 -1 = 9 8 -5 3 cos . GR i,min i,max 2 i,min 2 i,min 2 2 i

Here j e = 1i,min i,max 2

and we have assumed GR perturbations to be dominant over those from tidal and rotational distortion. The dimensionless quantity ε GR measures the 8strength9 of perturbations from GR relative to those of the binary companion and it is defined as

Extended Data Fig. 5 shows the required initial inclination between the planetary and binary orbital planes and the resulting initial maximum eccentricity with respect to the initial semimajor axis of the planet. Although vZLK oscillations are present when a i > 4.2 au, not all values above this threshold yield defined values for I i because the short-range forces are too strong for the planet to reach the required high initial eccentricity unless the initial semimajor axis exceeds a i > 7.0 au. The maximum eccentricity of the initial vZLK cycle must have been e i,max > 0.9947 to generate the present-day semimajor axis and eccentricity. Attaining a maximum eccentricity this large is only possible with a nearly polar initial inclination between the orbit planes of the planet and the binary companion. We find that the initial inclination I i > 86.8° for a i > 7.0 au. Altogether, these results indicate that it is possible to reach the present-day parameters of the system if the planet started beyond a i > 7.0 au and the binary companion started on an orbit nearly perpendicular to that of the planet.

We now use our derived constraints on the initial orbital conditions to explore the planetary orbital evolution through numerical simulations.

We conduct integrations of the secular of motion for TIC 241429530 through KozaiPy, a publicly available software package that simulates hierarchical three-body systems (https://github.com/ djmunoz/kozaipy). The equations of motion are provided in ref. 3. We adopt initial values of a i = 10 au and e i = 0.1, considering the minimum semimajor axis necessary for vZLK oscillations to be present. We consider perturbations to the octupole order and also account for tidal evolution in the constant-time-lag model of equilibrium tides 114 . Tidal parameters are adjusted so that the system reaches its present-day orbital parameters at an age of 3 Gyr, approximately equal to the derived age of the system. Specifically, the Love number of the planet is set to k 2 = 0.25 and its viscous timescale is set to t v = 0.01 days.

The simulation results are presented in Extended Data Fig. 6. They indicate the presence of vZLK oscillations that trigger periods of very large eccentricities. At the times that the periapse distance is minimized, tidal dissipation is strong and the semimajor axis shrinks. Eventually, the semimajor axis becomes small enough that the vZLK oscillations are suppressed owing to short-range forces and the planet decouples from the binary companion. After the vZLK oscillations are quenched, the mutual inclination is approximately conserved and the eccentricity of the planet slowly damps owing to continued tidal dissipation. We observe that there is an instant in time at which the eccentricity and semimajor axis of the planet are very close to the present-day values. We also note that the value of a e (1 -) max 2

is conserved during episodes of maximum eccentricity of each individual vZLK cycle, ranging within only a few percent of the average value of a e (1 -) max 2

. According to this simulation, continued dissipation will cause the planet to reach a circular orbit in about a billion years. Altogether, this simulation provides a plausible proof of concept of the system9s dynamical history of coupled vZLK oscillations and tidal migration. We suggest future work on the system to explore the role that dynamical tides might have played in its formation.

To explore the relation between exoplanet mass and eccentricity for warm Jupiters, or intermediate-period giant planets, we start with the sample of all transiting exoplanets with masses between 0.3 and 15 times that of Jupiter and orbital periods between 10 and 365 days. For each system in this sample, we adopt the most up-to-date mass and eccentricity constraints for which the eccentricity was fit as a free parameter when solving for the orbit. We discard four planets for which all literature solutions assumed a circular orbit with the eccentricity fixed at 0. All of these planets are less than 1.3 Jupiter masses. We also remove two planets in P-type circumbinary orbits, as the dynamical environments of these systems are expected to differ from those of planets orbiting single stars 115,116 . Our sample differs from that analysed in ref. 41, which draws from the RV planet sample and thus uses projected planet mass (M p sini) instead of true mass. By restricting our analysis to transiting exoplanets, we ensure that the measured masses are not degenerate with orbital inclination. This approach also mitigates the susceptibility of our results to detection biases, as the completeness fractions of transit surveys should be largely insensitive to exoplanet mass in the Jupiter-sized-planet regime.

The median mass of our sample is 1.935 Jupiter masses. We divide the sample into two groups of equal size, placing planets more massive than the median into one group and planets less massive than the median into the other group. The population-level eccentricity distribution of each group is then modelled in PyMC 117 using a hierarchical Bayesian framework. For our model, we adopt a beta distribution with two hyperparameters, θ = {µ, κ}, in which µ describes the mean of the distribution and 1/κ describes its variance. These hyperparameters represent a reparameterization of the standard beta distribution parameters α and β, in which α = µκ and β = (1 - µ)κ. A beta distribution is chosen for its flexibility in shape and because it is naturally bounded between 0 and 1. We adopt a uniform hyperprior for µ ~ U(0, 1) and a log-normal hyperprior for log κ ~ N(0, 3). These choices reduce the impact of hyperprior choices on the inference results, especially when the sample size is small, as is the case here 118 . The best-fit distributions are shown in Fig. 3 .22 . The mean values, µ, of the eccentricity distributions of low-mass and high-mass transiting warm Jupiters differ by 4.2σ. To assess the robustness of this result, we repeated the process for mass cutoffs between 1 Jupiter mass and 2.7 Jupiter masses. These bounds were chosen such that the size ratio of the two groups does not exceed 2:1. At each cutoff, we ran 1,000 trials, drawing the planet masses from asymmetric Gaussian distributions with means and widths determined by their literature values and uncertainties. For all mass-cutoff values over this range, the mean values of the two eccentricity distributions differ by 335σ.