TESS photometry

TIC 241249530 was observed with TESS⁹ 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. ^44,45,46,47). 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^48,49,50 (PDCSAP) light curve from the Science Processing Operations Center⁵¹ (SPOC) for Sector 60 and the TESS-SPOC⁵² light curve for Sector 20 (Extended Data Fig. 1).

High-contrast imaging

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)¹⁰ 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. ⁵³. The achieved 5σ contrast limits are sufficient to rule out the presence of faint stellar companions and background sources with Δmag₅₆₂ < 3.3 and Δmag₈₃₂ < 3.7 at a separation of 0.2″ and Δmag₅₆₂ < 3.9 and Δmag₈₃₂ < 4.8 at a separation of 1″.

Ground-based photometric observations

We used the Unistellar Network, a collaboration of citizen scientists using Unistellar telescopes⁵⁴ 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 15–30 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)¹⁶ 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¹⁵. We used the Semrock narrow-band filter (838–876 nm) to avoid atmospheric absorption bands⁵⁷. 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⁵⁸ 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.

Spectroscopic observations

We monitored the RV signal of TIC 241249530 with the NEID spectrograph¹¹ 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 ARCTIC, 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⁶⁰. 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. ⁶³. 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⁻¹ and σ_RV,DRP = 7.9 m s⁻¹, 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¹², which is on the Hobby–Eberly 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⁶⁶ and barycorrpy⁶⁷ 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⁻¹. 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⁷⁰ 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⁻¹, with an average precision of 0.1 km s⁻¹. The median resulting RV measurement precision is 3.4 m s⁻¹. We show the complete RV time series from NEID, HPF and HARPS-N in Extended Data Fig. 3.

Stellar characterization of TIC 241249530

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⁷³, MARCS atmospheric models⁷⁴, solar abundances from 3D hydrodynamic models⁷⁵ and the sixth version of the Gaia ESO survey (GES) atomic line list⁷⁶. 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⁷⁷. 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 Levenberg–Marquardt fitting algorithm, using constraints from the aforementioned models and line lists.

The HARPS-N spectra were independently analysed with BACCHUS⁷⁸, using MARCS atmospheric models, the GES atomic line list and the TURBOSPECTRUM 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⁷². 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. ^81,82,83. We use JHK_S magnitudes from 2MASS⁸⁴, W1–W3 magnitudes from WISE⁸⁵, G_BPG_RP magnitudes from Gaia⁸⁶, BVgri magnitudes from APASS⁸⁷ and the NUV magnitude from GALEX⁸⁸. We also used the Gaia spectrophotometry spanning 0.4–1.0 µm. Altogether, the available photometry spans the full stellar SED over the wavelength range 0.2–10.0 µm. We fit the SED using PHOENIX stellar atmosphere models⁸⁹, 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⁹⁰. 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⁻¹⁰ erg s⁻¹ cm⁻². 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⁹¹, 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 Yonsei–Yale isochrone models⁹². 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⁹², 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}_{-2.6}^{+3.8}\,{\rm{days}}\). 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⁹³ 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⁹⁴. 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 Lomb–Scargle 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.

Stellar characterization of TIC 241249532

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. ⁸⁶). The probability of a chance alignment between TIC 241249530 and TIC 241249532 is R = 9.73 × 10⁻⁵ (ref. ⁹⁵), suggesting that the pair is indeed gravitationally bound. Gaia’s 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 mass–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⁻¹ year⁻¹, with 99% of our orbits returning slopes <30 cm s⁻¹ year⁻¹. The amplitude of this signal is small compared with the km-s⁻¹-level variations induced by TIC 241249530 b and we therefore do not include it as an extra body in our joint fitting.

Joint transit + RV analysis

We use the exoplanet software package⁹⁸ to fit a transit model and a Keplerian orbit with Rossiter–McLaughlin 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. ¹⁰¹ to fit the transits, and the orbital parameter posteriors are sampled using the PyMC3 Hamiltonian Monte Carlo package¹⁰². The orbit model consists of a full Keplerian with tight Gaussian priors on the orbital period, P, and time of conjunction, T₀, 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 \(\sqrt{e}\sin \omega \) and \(\sqrt{e}\cos \omega \) 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¹⁰³. Separate quadratic limb-darkening coefficients, reparameterized as q₁ and q₂ as in ref. ¹⁰⁴, 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 Rossiter–McLaughlin signal, we adopt the formalism of ref. ¹⁰⁷ 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 spin–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.

Stellar obliquity

The stellar obliquity, ψ, is related to the projected obliquity, λ, by

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 \(\psi ={141}_{-24}^{{+15}^{\circ }}\). 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⁷. This is not definitive proof of the formation history, however, as retrograde orbits such as this can also be produced through planet–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.

Dynamical history—analytic constraints

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¹¹².

The planet is at present close enough to the primary star such that short-range forces—general relativity (GR), tides and rotational distortions—have 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¹¹³

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₂, a₂ and e₂ 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₂ = 0.453 M_⊙, a₂ = 1,664 au and e₂ = 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}_{{\rm{f}}}\equiv a\left(1-{e}_{\max }^{2}\right)\) 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. ¹¹²:

Here \({j}_{{\rm{i}},\min }=\sqrt{1-{e}_{{\rm{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 ‘strength’ 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.

Dynamical history—simulations

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 equations 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. ³. 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¹¹⁴. 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₂ = 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(1-{e}_{\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(1-{e}_{\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 system’s 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.

Modelling the transiting-warm-Jupiter eccentricity distribution

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. ⁴¹, which draws from the RV planet sample and thus uses projected planet mass (M_psini) 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¹¹⁷ 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¹¹⁸. The best-fit distributions are shown in Fig. 3 and the resulting hyperparameters are \({\mu }_{{\rm{low}}}={0.18}_{-0.03}^{+0.04}\), \(\log {\kappa }_{{\rm{low}}}={1.23}_{-0.29}^{+0.27}\), \({\mu }_{{\rm{high}}}={0.44}_{-0.05}^{+0.05}\) and \(\log {\kappa }_{{\rm{high}}}={0.91}_{-0.24}^{+0.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 3–5σ.