TOI-3686: A Hot Jupiter and Sub-Neptune In A 2:1 Mean-Motion Resonance

Email: [email protected]

Introduction

Most hot Jupiters—giant planets with orbital periods of less than ten days—appear to be isolated, without nearby planetary companions [80, 130, 64, 125]. This observation has been interpreted as evidence that hot Jupiters form farther out in the protoplanetary disk and arrive at their current close-in orbits via a dynamically disruptive, high-eccentricity migration pathway. In this picture, a wide-orbiting giant planet is excited onto a highly eccentric orbit through mechanisms like planet-planet scattering [110, 141] or secular interactions [45] (Y. Wu & N. Murray 2003; Y. Wu & Y. Lithwick 2011), before undergoing orbit circularization due to tidal interactions with the host star [44]. Since any interior planets would likely be disrupted during the high-eccentricity phase [14], such hot Jupiters would end up isolated.

The existence of hot Jupiters with nearby planetary companions implies that not all of them form through high-eccentricity migration. Some may either form in-situ or migrate inward by dynamically quiescent means, such as disk-driven migration (see, e.g., review by R. I. Dawson & J. A. Johnson 2018). The first example to be found was the WASP-47 system [15], which contains a hot Jupiter flanked by two smaller planets that were detected by the K2 mission. The prime Kepler mission [18] also discovered a hot Jupiter with an inner transiting companion, Kepler-730 [27].

However, hot Jupiters are intrinsically rare, occurring around $\sim1\%$ of all Sun-like stars [143], so Kepler observed only a few dozen hot Jupiter systems. NASA’s ongoing TESS mission [114] is observing hundreds of hot Jupiters across the sky [147], making it possible to find more transiting companions to hot Jupiters. Indeed, TESS has detected transiting planets both in previously known hot Jupiter systems (WASP-132, [63]; WASP-84, [89]), as well as in newly discovered systems (e.g., [65, 124, 77, 109, 90]).

Apart from observing their transits directly, close-by companions to hot Jupiters may also be detected through their gravitational influence on the hot Jupiter, causing transit timing variations (TTVs). Strong TTVs might especially be expected if hot Jupiters with companions arrive at their present locations through disk migration, which naturally captures neighboring planets into mean-motion resonances [100, 82, 75], and resonant interactions lead to larger and more detectable TTVs [e.g. 102] (Y. Lithwick et al. 2012). Indeed, a few of the hot Jupiters with known transiting companions show TTVs, with semi-amplitudes ranging from $\sim40$ seconds for WASP-47 b [99] to 8 minutes for TOI-1408 b [77] and 10 minutes for TOI-1130 c [76, 17].

In this manuscript, we report the discovery and characterization of the TOI-3686 system, a new member of the class of hot Jupiters with nearby companions. The system contains a hot Jupiter showing TTVs with a peak-to-peak amplitude of 40 minutes, some of the largest observed for a short-period giant planet. We determined that these TTVs are due to a nearby sub-Neptune in a mean-motion resonance with the hot Jupiter. We describe our observations in § 2, stellar host characterization in § 3, and analysis of light curve data in § 4. We use $N$-body simulations to characterize the system in § 5, before discussing its dynamical state (§ 6) and its relation to the broader sample of hot Jupiters with companions (§ 7). We summarize in § 8.

Observations

TESS Photometry

TESS observed TOI-3686 in Sectors 18, 58, and 85, spanning a five year observation baseline from 2019 through 2024. The TESS data show clear 1.5% dips in the star’s brightness occurring every 5.078 days (Figure 1, Figure 2). The events in Sector 18 were alerted as a transiting planet candidate, TOI-3686.01, by the “faint-star search” [78] carried out for stars observed in the TESS Full Frame Images (FFIs), using light curves produced by the Quick Look Pipeline (QLP; [66 c]). Following the identification of the transit events, the star was selected for higher cadence observations when it was again in TESS’ field of view in Sectors 58 and 85. We downloaded the TESS data for this target using the lightkurve package (Lightkurve Collaboration et al. 2018), selecting the light curves produced by the TESS-SPOC pipeline from the FFIs at 1800s cadence in Sector 18, and from the postage stamps at 120s cadence in Sectors 58 and 85 [70, 26].

The TESS light curves of TOI-3686 also show longer term variability with a peak-to-peak amplitude of 1–2% and which is quasi-periodic on a timescale of $\sim5.55$ days (Figure 2). We attribute this photometric modulation to the rotation of starspots in and out of view, such that the periodicity is at the stellar rotation period, and discuss this signal further in § 3.1.

Ground-Based Photometry

We also obtained time-series transit photometry of TOI-3686.01 from telescopes on the ground. These were made initially for false positive vetting, but once the presence of TTVs were noted, we began monitoring the transits more intensively to measure additional transit times. Most of this follow-up was coordinated through the TESS Follow-up Observing Program Sub Group 1 (TFOP SG1; [32]), which focuses on obtaining seeing-limited time-series photometry of TESS candidates. We used the TESS Transit Finder, a customized version of the Tapir software package ([71]), to schedule the transit observations. We obtained a total of 45 transit light curves from professional facilities, college observatories, as well as private observatories operated by citizen scientists across North America, Europe and Asia. Table 2 lists observation details and all facilities used. Most of the ground-based lightcurves were extracted using AstroImageJ ([33]).

Epoch	Facility/Instrument	Location	Aperture (m)	Filter	Date (UT)	Cadence (s)	Transit Time$^{a}$ (BJD$_{\mathrm{TDB}}$)
0	TESS Sector 18	Space	0.1	TESS	2019 Nov 05	1800	2458792.55799(67)(78)
1	TESS Sector 18	Space	0.1	TESS	2019 Nov 10	1800	2458797.63638(71)(68)
3	TESS Sector 18	Space	0.1	TESS	2019 Nov 20	1800	2458807.78781(70)(66)
4	TESS Sector 18	Space	0.1	TESS	2019 Nov 25	1800	2458812.86795(84)(59)
138	OAUV/T50	Valencia, Spain	0.5	$R$	2021 Oct 05	75	2459493.35915(94)(90)
143	Waffelow Creek Observatory	Nacogdoches, TX, USA	0.36	$r'$	2021 Oct 31	60	2459518.74545(59)(48)
152	Observatory de Ca l’Ou	Sant Martí Sesgueioles, Barcelona, Spain	0.4	$R$	2021 Dec 15	100	2459564.4406(10)(11)
153	SUTO/OTIVAR	Pyskowice, Poland	0.3	$R$	2021 Dec 20	90	2459569.51586(91)(81)
215	TESS Sector 58	Space	0.1	TESS	2022 Oct 31	120	2459884.34713(61)(50)
216	TESS Sector 58	Space	0.1	TESS	2022 Nov 05	120	2459889.42431(82)(69)
217	TESS Sector 58	Space	0.1	TESS	2022 Nov 11	120	2459894.50254(55)(45)
218	TESS Sector 58	Space	0.1	TESS	2022 Nov 16	120	2459899.58002(53)(52)
219	TESS Sector 58	Space	0.1	TESS	2022 Nov 21	120	2459904.65817(52)(51)
219	Acton Sky Portal	Acton, MA, USA	0.36	$g'$	2022 Nov 21	32	2459904.65646(40)
220	TESS Sector 58	Space	0.1	TESS	2022 Nov 26	120	2459909.73720(51)(57)
278	Kourovka Observatory	Yekaterinburg, Russia	0.4	$R$	2023 Sep 16	240	2460204.26530(88)(92)
278	Kourovka Observatory	Yekaterinburg, Russia	0.4	$V$	2023 Sep 16	240	2460204.2694(13)
337	Pine Mountain Observatory/CDK14	Bend, OR, USA	0.35	$r'$	2024 Jul 12	110	2460503.88499(49)(56)
344	OAUV/TURIA2	Valencia, Spain	0.3	$R$	2024 Aug 16	195	2460539.44256(51)
345	Kotizarovci Observatory	Viskovo, Croatia	0.3	$V$	2024 Aug 21	120	2460544.51810(98)(72)
348	Peters Observatory	Clinton, NY, USA	0.5	$r'$	2024 Sep 05	60	2460559.75519(25)(28)
348	Nube Blanca Observatory	Berthoud, CO, USA	0.4	$R$	2024 Sep 06	160	2460559.75356(72)(81)
349	FLWO/KeplerCam	Amado, AZ, USA	1.2	$i'$	2024 Sep 11	40	2460564.83296(18)(17)
349	LCO McDonald/Sinistro	Fort Davis, TX, USA	1.0	$z'$	2024 Sep 11	67	2460564.83224(23)(24)
349	Lookout Observatory	Colorado Springs, CO, USA	0.28	$V$	2024 Sep 11	20	2460564.8311(14)(12)
349	Union College Observatory	Schenectady, NY, USA	0.5	$V$	2024 Sep 11	60	2460564.83450(25)
358	Acton Sky Portal	Acton, MA, USA	0.36	$r'$	2024 Oct 26	20	2460610.52944(30)(31)
358	Privat Observatory Herges-Hallenberg	Steinbach-Hallenberg, Germany	0.28	$BB$	2024 Oct 26	120	2460610.5280(37)(30)
359	TESS Sector 85	Space	0.1	TESS	2024 Nov 01	120	2460615.60747(57)
360	TESS Sector 85	Space	0.1	TESS	2024 Nov 06	120	2460620.68311(59)(54)
360	Mason Observatory	Oklahoma City, OK, USA	0.28	$V$	2024 Nov 05	90	2460620.6848(11)(13)
360	LCO McDonald/Sinistro	Fort Davis, TX, USA	1.0	$z'$	2024 Nov 06	95	2460620.68322(24)(26)
Epoch	Facility/Instrument	Location	Aperture (m)	Filter	Date (UT)	Cadence (s)	Transit Time$^{a}$ (BJDTDB)
361	LCO McDonald/Sinistro	Fort Davis, TX, USA	1.0	$z'$	2024 Nov 11	67	2460625.75964(23)(25)
362	TESS Sector 85	Space	0.1	TESS	2024 Nov 16	120	2460630.83783(57)(61)
362	LCO McDonald/Sinistro	Fort Davis, TX, USA	1.0	$z'$	2024 Nov 16	67	2460630.83758(47)(43)
363	TESS Sector 85	Space	0.1	TESS	2024 Nov 21	120	2460635.91428(51)(49)
369	Kotizarovci Observatory	Viskovo, Croatia	0.3	$V$	2024 Dec 21	120	2460666.3719(16)(17)
371	LCO Teide/Sinistro	La Laguna, Tenerife, Spain	1.0	$i'$	2024 Dec 31	17	2460676.52554(32)(28)
373	FLWO/KeplerCam	Amado, AZ, USA	1.2	$i'$	2025 Jan 11	60	2460686.68456(24)(29)
374	BYU OPO	Provo, UT, USA	0.3	$V$	2025 Jan 16	75	2460691.7617(18)(24)
411	OAUV/TURIA2	Valencia, Spain	0.3	$R$	2025 Jul 22	195	2460879.66949(77)(72)
421	OAUV/TURIA2	Valencia, Spain	0.3	$R$	2025 Sep 11	195	2460930.43795(43)(50)
422	OAUV/TURIA2	Valencia, Spain	0.3	$R$	2025 Sep 16	195	2460935.51959(40)(36)
423	LCO Teide/Sinistro	La Laguna, Tenerife, Spain	1.0	$i'$	2025 Sep 21	42	2460940.59470(20)(19)
424	LCO Teide/Sinistro	La Laguna, Tenerife, Spain	1.0	$i'$	2025 Sep 27	42	2460945.67077(20)(18)
424	Union College Observatory	Schenectady, NY, USA	0.5	$R$	2025 Sep 27	250	2460945.6704(13)(12)
425	Pritchard Observatory	York, PA, USA	0.2	$V$	2025 Oct 01	60	2460950.74219(46)(39)
425	Union College Observatory	Schenectady, NY, USA	0.5	$R$	2025 Oct 02	69	2460950.7495(13)(15)
426	BYU OPO	Provo, UT, USA	0.6	$R$	2025 Oct 07	150	2460955.82511(48)(46)
433	Maidanak Astronomical Observatory	Qashqadaryo, Uzbekistan	0.6	$R$	2025 Nov 11	60	2460991.37734(66)(56)
437	LCO McDonald/Sinistro	Fort Davis, TX, USA	1.0	$i'$	2025 Dec 02	17	2461011.69166(20)(21)
437	Nube Blanca Observatory	Berthoud, CO, USA	0.4	$R$	2025 Dec 02	60	2461011.6902(12)(11)
439	Unistellar evScope	Placitas, NM, USA	0.1	$V$	2025 Dec 12	120	–
445	Privat Observatory Herges-Hallenberg	Steinbach-Hallenberg, Germany	0.28	$V$	2026 Jan 11	120	2461052.3283(11)(09)
446	Salerno University Observatory	Salerno, Italy	0.6	$R$	2026 Jan 16	120	2461057.4076(10)(11)
448	Marvel Telescope	La Palma, Spain	0.8	$g'$	2026 Jan 26	45	2461067.56733(49)(58)
486	SUTO/OTIVAR	Pyskowice, Poland	0.3	$R$	2026 Aug 07	90	2461260.5154(16)(15)
488	Two-meter Twin Telescope	La Laguna, Tenerife, Spain	0.8	$i'$	2026 Aug 18	40	2461270.67621(61)(78)
489	Two-meter Twin Telescope	La Laguna, Tenerife, Spain	0.8	$i'$	2026 Aug 23	40	2461275.7549(15)

These observations resolved TOI-3686 from all nearby bright stars and confirmed that the transit events occur on-target with achromatic depths, favoring a planetary interpretation for the transits. The geographic dispersion of the observers greatly extended the TTV baseline and improved phase coverage of the TTV curve compared with the TESS data alone.

Speckle Imaging

To check for even closer stellar companions to TOI-3686, we also observed it with the Speckle Polarimeter [115, 131] on the 2.5 m telescope at the Caucasian Mountain Observatory (CMO) of Sternberg Astronomical Institute (SAI) of Lomonosov Moscow State University. These observations achieved a limiting contrast ratio of $\Delta I = 5.2$ mag in $I$-band at a separation of 1″, with no nearby sources detected down to instrumental limits.

We note that Gaia DR3 (Gaia Collaboration et al. 2022) lists a faint source (Gaia DR3 465296935576885376, $G = 18.0$ mag) separated by 1″.66 from TOI-3686. No astrometric data apart from position were available, so it is unclear whether this is a bound or background object. Given

this source’s faintness ($\Delta G = 6.24$ mag), it is impossible for it to be the source of either the transit events or long-term photometric modulation; thus, we do not consider this object further.

Spectroscopic Observations

We observed TOI-3686 with several high-resolution spectrographs, which we used to characterize the host star as well as measure precise radial velocities (RVs). All RV data are provided in a machine-readable table accompanying this paper.

**TRES—**We obtained 11 observations between Sep 2020 and Nov 2024 with the Tillinghast Reflector Echelle Spectrograph (TRES; [46]) on the 1.5m telescope at the Fred Lawrence Whipple Observatory (FLWO). The spectra were reduced according to the pipeline described in L. A. Buchhave et al. [24] and S. N. Quinn et al. [108]. We derived relative RVs by first constructing a median-filtered template from all the observations and then cross-correlating this template against each individual observation. Long-term changes in the zero-point offset due to instrumental effects are corrected using standard star observations.

**NEID—**We obtained 41 spectra with the NEID spectrograph ([119, 56]) on the WIYN 3.5m telescope at Kitt Peak National Observatory (KPNO) between Oct 2021 and Mar 2025. We observed in High Resolution (HR) mode ($R \sim110{,}000$) with exposure times of 300–900s. The data were reduced and RVs extracted using v1.4.2 of the NEID Data Reduction Pipeline.² The pipeline measures RVs by performing a cross-correlation between each observed spectrum and a weighted line mask derived for the star’s spectral type ([7, 105]). These observations spanned several “eras” of NEID operations,³ between which offsets of a few $\mathrm{m\,s^{-1}}$ may occur, so we treat the different eras as having independent instrumental offsets. In our subsequent analysis, we used only the 31 observations taken during the third RV era, as the earlier data were more sparsely sampled and had lower S/N.

**MAROON-X—**We also observed TOI-3686 with the MAROON-X spectrograph ([120, 121, 122, 123]) on the Gemini-North telescope (Prog. ID GN-2024B-FT-108; PI: Yee) on Maunakea, Hawaii, obtaining 23 epochs with 600s exposure times. The data were reduced using the standard instrument pipeline.⁴ Precise RVs were measured using a modified version of the SERVAL package ([149]), which constructs an empirical co-added template from all observations, and then computes order-by-order shifts between each observation and the template through least squares fitting. A combined RV for each of the blue and red channels of MAROON-X are derived from a weighted mean of the RVs across all orders in each arm. To account for potential offsets between the two arms and between observing runs ([9]), we treated each combination of the above as separate instruments.

Preliminary RV modeling

We performed a preliminary fit of a single-planet Keplerian to the RV data, constrained to the best-fit linear ephemeris from the transit data. This produced a reasonable fit to the data, with an RV semi-amplitude of $K = 55.5^{+3.6}_{-3.3}\ \mathrm{m\,s^{-1}}$, corresponding to a planetary mass of $M_{\mathrm{p}}\sin i = 0.48 \pm0.04 M_{\mathrm{J}}$, although this measurement does not account for other sources of RV variation. Together with the ground-based photometry and imaging, this confirms TOI-3686.01 as a bona fide planet with a mass roughly half that of Jupiter.

However, the RV residuals show significant scatter, and their Lomb-Scargle periodogram ([87, 116]) peaks close to the stellar rotation period of $P_{\mathrm{rot}} \sim5.55$ days, motivating us to include a Gaussian Process (GP) model to account for activity-induced RV variations in our later modeling (§ 5).

Stellar Characterization

We used our high-resolution spectroscopic observations of TOI-3686 to measure stellar atmospheric parameters. We applied the Stellar Parameter Classification code (SPC; [25]) to match the observed TRES spectra against the R. L. Kurucz [79] grid of synthetic model spectra. We took the mean and standard deviation as the final value and uncertainty of each parameter, adopting error floors of 50 K in $T_{\mathrm{eff}}$, 0.1 dex in $\log g$, 0.08 dex in [Fe/H], and 0.5 $\mathrm{km\,s^{-1}}$ in $v\sin i$ to account for systematic model uncertainties. This analysis yields $T_{\mathrm{eff}} = 6047 \pm75$ K, $\log g = 4.35 \pm0.13$ dex, [Fe/H] $= 0.09 \pm0.08$ dex, and $v\sin i = 11.6 \pm0.5\ \mathrm{km\,s^{-1}}$.

We also analyzed a co-added NEID spectrum with SpecMatch-Empirical ([146]), which compares the target spectrum to an empirical library References— Catalog astrometry and photometry are drawn from Gaia DR3 (Gaia Collaboration et al. 2022), 2MASS ([34, 127]), WISE ([35]), and the TESS Input Catalog ([128, 129]).

of observed spectra of well-characterized stars. We obtained stellar parameters consistent with the results from SPC: $T_{\mathrm{eff}} = 6050 \pm110$ K, $R_{\star} = 1.12 \pm0.20 R_{\odot}$, and [Fe/H] $= 0.04 \pm0.09$ dex. We adopt the SPC parameters, which have smaller uncertainties.

To complement the stellar atmospheric parameters, we used the EXOFASTv2 code (J. Eastman et al. [42]; J. D. Eastman et al. [43]) to constrain the physical parameters of TOI-3686 by modeling its broadband photometry with MIST stellar evolutionary models (A. Dotter [41]; J. Choi et al. [30]). We fitted photometry from the Gaia DR3, 2MASS, and WISE catalogs, imposing error floors of 0.02, 0.02, and 0.03 mag on the photometric uncertainties respectively. We also used as input the zero-point corrected parallax from Gaia (Gaia Collaboration et al. (2022); L. Lindegren et al. [84]). We imposed a Gaussian prior on [Fe/H] from our SPC analysis and an upper limit on the visual extinction $A_V$ from galactic dust maps from D. J. Schlegel et al. [118]; ); E. F. Schlafly & D. P. Finkbeiner [117]. We sampled the posterior using a differential evolution Markov Chain Monte Carlo (DE-MCMC) algorithm, requiring a Gelman-Rubin statistic below 1.01 for convergence (A. Gelman & D. B. Rubin [53]).

Our final stellar properties are presented in Table 1, indicating that TOI-3686 is a star slightly more massive and metal-rich than the Sun. The MIST modeling provides only a wide constraint on the age of TOI-3686; however, a few indicators of youth suggest that it may have an age $\lesssim1$ Gyr, which we discuss below.

Parameter	Value	Source
TIC ID	13533330	TICv8
Gaia DR3 ID	465296935580218368	Gaia DR3
2MASS ID	02273616+5957098	2MASS
WISE ID	J022736.14+595709.8	WISE
Astrometric Measurements
R.A. (J2016.0)	02:27:36.13	Gaia DR3
Dec (J2016.0)	+59:57:09.79	Gaia DR3
$\mu_{\alpha*}$ (mas yr$^{-1}$)	$-14.012 \pm0.012$	Gaia DR3
$\mu_{\delta}$ (mas yr$^{-1}$)	$-3.428 \pm0.013$	Gaia DR3
$\varpi$ (mas)	$3.920 \pm0.013$	Gaia DR3
Catalog Photometry
$T$ (mag)	$11.2380 \pm0.0060$	TICv8
$G$ (mag)	$11.7659 \pm0.0028$	Gaia DR3
$G_{\mathrm{BP}}$ (mag)	$12.1898 \pm0.0031$	Gaia DR3
$G_{\mathrm{RP}}$ (mag)	$11.1681 \pm0.0039$	Gaia DR3
$J$ (mag)	$10.427 \pm0.020$	2MASS
$H$ (mag)	$10.148 \pm0.018$	2MASS
$K_{\mathrm{s}}$ (mag)	$10.038 \pm0.016$	2MASS
W1 (mag)	$9.943 \pm0.041$	WISE
W2 (mag)	$9.988 \pm0.020$	WISE
W3 (mag)	$9.754 \pm0.061$	WISE
Spectroscopic Parameters
$T_{\mathrm{eff}}$ (K)	$6047 \pm75$	SPC
$\log g$ (cgs)	$4.35 \pm0.13$	SPC
[Fe/H] (dex)	$0.086 \pm0.080$	SPC
$v \sin i_{\star}$ (km s$^{-1}$)	$11.64 \pm0.50$	SPC
Physical Stellar Parameters
$M_{\star}$ ($M_{\odot}$)	$1.049^{+0.092}_{-0.085}$	EXOFASTv2
$R_{\star}$ ($R_{\odot}$)	$1.184^{+0.047}_{-0.044}$	EXOFASTv2
$T_{\mathrm{eff}}$ (K)	$5990^{+280}_{-270}$	EXOFASTv2
$\log g$ (cgs)	$4.313^{+0.050}_{-0.052}$	EXOFASTv2
[Fe/H] (dex)	$0.042^{+0.084}_{-0.076}$	EXOFASTv2
Age (Gyr)	$5.9^{+4.4}_{-3.5}$	EXOFASTv2
Rot. Age (Gyr)	$0.73^{+0.62}_{-0.23}$	gyro-interp
Li. Age (Gyr)	$1.27^{+2.2}_{-1.1}$	EAGLES

Indications of youth for TOI-3686

The TESS light curves of TOI-3686 show significant, 1–2% peak-to-peak quasi-periodic variability that is likely due to rotational modulation. Modeling this with a GP (§ 4.2), we measure a mean $P_{\mathrm{rot}} = 5.55 \pm0.06$ days. Assuming that its rotation has only been modified by magnetic braking, this moderately rapid rotation indicates that TOI-3686 is somewhat youthful compared to the typical field star. We used the gyro-interp code (L. G. Bouma et al. [20]) to infer a gyrochronological age of $730^{+620}_{-260}$ Myr. One caveat is that gyrochronological ages may not be reliable for stars that host hot Jupiters due to tidal spin-up (e.g., D. J. A. Brown [22]; P. F. L. Maxted et al. [93]; R. A. Tejada Arevalo et al. [134]), although the relatively low mass ($M_{\mathrm{p}} \approx0.5 M_{\mathrm{J}}$) and long orbital period ($P = 5.078$ days) of TOI-3686.01 would weigh against significant modification of its host star’s spin period.

Another indicator of youth is the presence of lithium in the stellar atmosphere, as stars with significant convective envelopes should transport Li into the interior where it is destroyed, depleting photospheric Li over time. Using the high S/N stacked NEID spectrum (§ 3), we detected significant absorption in the Li I 6707.8 Å line, with an equivalent width (EW) of $51.8^{+1.6}_{-1.9}$ mÅ. This was measured following the procedure from L. G. Bouma et al. [19], which corrects for blending of the Li line with a nearby Fe line by subtracting a calibrated mean offset. We then used the EAGLES code (R. D. Jeffries et al. [69]) to estimate the stellar age given the stellar $T_{\mathrm{eff}}$ and Li EW measurement, finding an age of $1270^{+2200}_{-1100}$ Myr. The weaker constraining power of the Li method for stars with this $T_{\mathrm{eff}}$ is expected, given the limited depth of their convection zones. Still, both the Li absorption and gyrochronology appear to point consistently to a relatively young, $\lesssim1$ Gyr age for TOI-3686.

Light Curve Analysis

Transit Time Extraction

To characterize the TTVs of TOI-3686.01, we measured individual transit times from the space- and ground-based observations. We used the exoplanet code (D. Foreman-Mackey et al. 2021) to model the photometry, which implements the differentiable transit model of E. Agol et al. [2] for rapid evaluation of transit light curves and their derivatives. We modeled each observation with its own mid-transit time $T_{\mathrm{c},i}$ and derived a mean orbital period $P_{\mathrm{orb}}$ from a least-squares fit to these transit times. The other transit parameters $b$, $R_{\mathrm{p}}/R_{\star}$, and stellar pseudo-density $\tilde{\rho}_{\star}$, were shared across all transits, and with broad uniform priors on these parameters. We assumed a quadratic limb-darkening law, with coefficients shared across observations in the same filter, sampling in the $(q_1, q_2)$ basis from D. M. Kipping [74].

We simultaneously detrended each ground-based light curve against a linear combination of instrumental systematics such as pixel position, airmass, and point spread function (PSF) width. For the TESS observations, we first masked the transit events and fitted a basis spline to the out-of-transit flux with the Keplerspline code to remove the long-term variability ([137, 126]). We then split each sector into separate light curves for each transit, retaining $\pm2$ transit durations of data around the transit center. To account for remaining uncorrected flux variations in the TESS data, we allowed for simultaneous detrending against a quadratic function of time.

We report the best-fit transit times and $1$-$\sigma$ uncertainties in Table 2. The TTVs appear well-described by a sinusoidal model, with a peak-to-peak amplitude of 40 minutes, some of the largest ever seen for a hot Jupiter, and a period of $\sim260$ days (Figure 3). These TTVs clearly indicate the presence of a nearby planet in or close to a mean-motion resonance (MMR; e.g. [102] (Y. Lithwick et al. 2012)), motivating us to search for additional transits in the TESS data.

Detection of a Second Transiting Planet

Previous searches of the TESS data had not detected any other transits, but shallow events could have been missed given the stellar variability and potentially large TTVs. We thus decided to perform our own search of the TESS light curves.

We again used the exoplanet code (D. Foreman-Mackey et al. 2021) to model the TESS data, this time using the SPOC Simple Aperture Photometry (SAP) light curves, as these retained more of the stellar variability signal than those that have undergone Pre-search Data Conditioning (PDCSAP light curves). We fitted the full light curve from each sector separately to allow for changes in stellar variability and the orbital ephemeris due to TTVs. We first included only the hot Jupiter’s transits and a GP model for the stellar variability, where we used the RotationTerm GP kernel implemented in celerite2 ([49, 48]), which models stellar rotation as a mixture of two stochastically-driven and damped harmonic oscillators. We assumed a constant linear ephemeris within each sector, since we found no evidence for high-frequency “chopping” signals, while the $\sim25$ day sector baseline is only $10\%$ of the TTV super-period such that the long-term TTV signal would be absorbed by changes in the mean orbital period and time of conjunction. We then optimized over the hot Jupiter’s transit and orbital parameters and the GP hyperparameters.

We then subtracted the MAP single-planet and GP model from the data, and searched the residuals for additional transit signals using the transit least squares (TLS) algorithm (M. Hippke & R. Heller [61]). A peak at $\sim2.54$ days, roughly half the hot Jupiter’s orbital period, is the most significant periodogram peak in all three sectors, with Signal Detection Efficiencies (SDEs) of 6.1, 7.4, and 5.3 respectively (Figure 2). No statistically significant peaks were seen at either the hot Jupiter’s orbital period or the stellar rotation period, indicating that our model removed those signals effectively.

We then added the transits of a second planet to our model, placing uniform priors on its orbital period and time of transit centered on the periodogram peaks. As before, we assumed a constant period within each sector. We used the Hamiltonian Monte Carlo (HMC) No U-Turn Sampler (M. D. Hoffmann & A. Gel(man 2014) implemented in pymc (O. Abril-Pla et al. [1]) to explore the posterior probability distributions of the fitted parameters in this model. Each sector’s fit was conducted independently, yet we found properties for the second planet that were consistent between them; notably the transit depths yielded $R_p/R_\star= 0.024 \pm0.006$, $0.023 \pm0.003$, and $0.024 \pm0.007$ respectively, strongly suggesting that they are indeed due to a planet and not spurious noise. The transit ephemerides of this second planet vary between sectors, but we show in § 5 that these variations are consistent with a TTV signal 180° out of phase with that of the hot Jupiter. The best-fit transit models along with random draws from the posterior, are plotted in Figure 2.

Given that the transit depths are measured at 4-, 7- and 3.5-$\sigma$ significance in individual sectors with consistent shapes and depths, we regard this as a confident detection of a second transiting planet in the TOI-3686 system. The transits correspond to a sub-Neptune-sized ($R_p \approx3 R_\oplus$) planet in a 1:2 period commensurability

with the hot Jupiter. In the remainder of this paper, we focus our analysis on this system architecture, referring to the inner sub-Neptune as TOI-3686 b and the outer hot Jupiter as TOI-3686 c.

**Note—**Acronyms for facility names are provided here. OAUV: Observatori Astronòmic de la Universitat de València; SUTO: Silesian University of Technology Observatories; FLWO: Fred Lawrence Whipple Observatory; LCO: Las Cumbres Observatory ([23]); BYU OPO: Brigham Young University Orson Pratt Observatory.

Photodynamical Analysis

To leverage the rich information encoded by the TTVs about the planetary masses, eccentricities and dynamical states, we modeled the TOI-3686 planets with 3-body integrations. We used the jnkepler code (K. Masuda et al. [92]; K. Masuda [91]), an $N$-body code written in JAX, which enables automatic differentiation of the model with respect to the fitted parameters and efficient sampling of their posterior distributions with HMC algorithms. Given a set of initial conditions for the two planets, we integrate the evolution of the system for 2500 days, covering full observational baseline.

Because the transits of the inner sub-Neptune are too shallow to measure individual transit times, we implemented a “photodynamical” model, directly fitting the TESS flux time series. In jnkepler, the flux loss due to each planet’s transits is computed using the modeled sky-plane positions and velocities at the transit mid-point, with the E. Agol et al. [2] transit model implemented in jaxoplanet (S. Hattori et al. [59]). The model for Sector 18 was evaluated at 2-minute cadence and integrated over the actual 30-minute observation cadence to account for the finite integration time (D. M. Kipping [73]). We also allowed for potential dilution of the transit depth from nearby stars by fitting for a flux dilution term $D_s$ (e.g., T. Han et al. [58]), as well as a flux jitter term to account for underestimated flux uncertainties. As in § 4.2, we modeled the stellar variability in each TESS sector with a RotationTerm GP kernel. Thus, the likelihood for the data from TESS sector $s$ is given by:

where $\mathcal{N}(\mathbf{r}_s \mid\boldsymbol{\mu}, \sigma^{2})$ denotes the multivariate Gaussian probability density evaluated at $\mathbf{r}_s$; $\mathbf{K}_{\mathrm{rot},s}$ is the RotationTerm kernel; and $\mathbf{r}_s$ is the residual vector taken by subtracting the 2-planet transit model from the observed fluxes.

To save on computational cost, we did not model the ground-based photometry directly but fitted the mid-transit times for those observations extracted in § 4.1. We assumed a standard Gaussian likelihood, including a global timing jitter term added in quadrature with each timing measurement’s uncertainty.

Finally, our preliminary analysis (§ 2.4.1) revealed that the RV data also contain a significant activity-induced signal with periodic modulation at the stellar rotation period. Thus, our model for the RV data includes both the stellar velocities computed from the $N$-body integration at the actual observation times, as well as a GP model for stellar variability, again using the RotationTerm kernel.

Because the majority of the RV data are contemporaneous with the TESS Sector 85 observations, we shared the rotation period hyperparameter $P_{\mathrm{GP}}$ between the RV GP model and that for the Sector 85 photometry, but treated the remaining GP hyperparameters independently. We also fitted for a per-instrument offset $\gamma_i$ and jitter, such that each instrument $i$ contributes the likelihood:

and $\mathbf{r}_i = \mathbf{RV}_{\mathrm{obs},i} - \mathbf{RV}_{N\text{-}\mathrm{body}} - \gamma_i$ is the residual vector for the RV observations.

The total log-likelihood is then

where the first sum is over the three TESS sectors and the second is over the five RV instruments.

We varied all orbital elements, masses, and $R_p/R_\star$ for both planets. Instead of mean anomaly at the initial epoch, we sampled the osculating time of conjunction at the start of integration, which allowed more intuitive bounds to be placed on this parameter. We parameterized the planetary inclinations in terms of impact parameters and bounded them between 0 and 1. To convert from mass and radius ratios with the star to physical planet properties, we also fitted for the stellar mass and radius, imposing Gaussian priors centered at the values from § 3. We note that the stellar mass and radius are linked in the photodynamical model because the transit durations depend on the stellar density for a given orbital period. Because the fitted dilution terms are degenerate with the planet-to-star radius ratios in affecting the measured TESS transit depths, we imposed a Gaussian prior on $R_p/R_\star$ for the hot Jupiter, derived from a fit similar to that described in § 4.1, but using the ground-based data only. Since the ground-based observations resolved TOI-3686 from all neighboring bright stars, we treat the transit depth derived from those observations as the undiluted ground truth. The remaining fitted parameters received broad uniform or log-uniform priors, listed in Table 3.

We first investigated the solution landscape by conducting a grid search of possible initial conditions: we randomly sampled system parameters within their prior ranges, then performed a local optimization over the model posterior probability with the L-BFGS-B algorithm. We found multiple modes corresponding to different orbital periods of TOI-3686 b, arising from a different integer number of elapsed orbits between TESS

sectors. We then performed a fine sweep in $P_b$ to more closely sample each alias. Ultimately, we found that the mode corresponding to 722 orbits of TOI-3686 b between the reference $T_0$ in Sectors 18 and 85 was most favored, with the next best mode having significantly lower posterior probability ($\Delta\ln p = -33.8$).

We used the maximum a posteriori (MAP) result from this mode to initialize a numpyro HMC sampler ([16, 106]), running four parallel HMC chains with 1500 warm-up steps and 1000 sampling steps. We checked for convergence by ensuring the rank-normalized $\hat{R}$ statistic ([138]) was below 1.01 for all parameters, and an effective sample size ESS > 500.

Figure 3 shows the data and predictions from the best-fitting photodynamical model, while the parameter inferences are reported in Table 3. The GP hyperparameters for the stellar flux variability model are similar to those obtained in our earlier light curve analysis (§ 4.2), while for the RVs, the activity GP model has a standard deviation $\sigma= 30\ {\rm m\ s^{-1}}$. The predicted TTV model for TOI-3686 c is essentially sinusoidal, with a peak-to-peak amplitude of 40 minutes, while the TTVs of the inner planet have a peak-to-peak amplitude of 10 hours, $\sim16\%$ of its orbital period. The period ratio oscillates around a mean value of $P_c/P_b = 2.0025$, corresponding to a deviation from perfect period commensurability of

This small value of $\Delta$ indicates that the system is likely in, not merely close to, an MMR ([54]), which we explore further in § 6.

The modeling of the RV and TTV data allow us to constrain the masses of both planets. We find that TOI-3686 c has a mass of $M_p = 0.56 \pm0.04\ M_J$, consistent with that measured in our preliminary RV fit (§ 2.4.1). With an equilibrium temperature of $T_{\rm eq} = 1300\ {\rm K}$, the planet is subject to the “hot Jupiter radius anomaly” that is active above $T_{\rm eq} \gtrsim1000\ {\rm K}$ ([40, 97]). The measured $R_p = 1.27 \pm0.03\ R_J$ is consistent with empirical relationships for this anomalous inflation ([135]) and is significantly inflated from the $R_p \sim0.9\ R_J$ that would be expected for a planet of this mass at lower insolation levels ([136]). Meanwhile, TOI-3686 b has a mass and radius typical of a sub-Neptune, ($M_p = 10.5^{+1.3}_{-1.2}\ M_\oplus$, $R_p = 2.8 \pm0.2\ R_\oplus$). We note that this mass constraint comes entirely from the TTVs of TOI-3686 c, as its expected RV semi-amplitude of $K \approx5\ {\rm m\ s^{-1}}$ is too small to be detected by the RV data alone.

Table 3 continued

Parameter	Prior	Posterior
Planet b
$P_b$ (days)	$\mathcal{U}(2.4886, 2.5389)$	$2.54320^{+0.00066}_{-0.00086}$
$T_{c,b}$ (BJD $-T_{\rm ref}$)	$\mathcal{U}(0.8, 1.7)$	$1.384^{+0.015}_{-0.013}$
$e_b \cos\omega_b$	$\mathcal{U}(-0.4, 0.4)^a$	$0.019^{+0.013}_{-0.009}$
$e_b \sin\omega_b$	$\mathcal{U}(-0.4, 0.4)^a$	$-0.1584 \pm0.0026$
$e_b$	Derived $^a$	$0.1595^{+0.0013}_{-0.0016}$
$M_b$ ($M_\odot$)	$\mathcal{LU}(10^{-6}, 10^{-3})$	$0.0000315^{+0.0000040}_{-0.0000036}$
$M_b$ ($M_\oplus$)	Derived	$10.5^{+1.3}_{-1.2}$
$R_b/R_\star$	$\mathcal{U}(0.001, 0.2)$	$0.0231^{+0.0015}_{-0.0014}$
$R_b$ ($R_\oplus$)	Derived	$2.84^{+0.19}_{-0.20}$
$b_b$	$\mathcal{U}(0, 1)$	$0.41^{+0.18}_{-0.24}$
$a_b/R_\star$	Derived	$7.23 \pm0.14$
Planet c
$P_c$ (days)	$\mathcal{U}(5.0270, 5.1286)$	$5.076790^{+0.000074}_{-0.000068}$
$T_{c,c}$ (BJD $-T_{\rm ref}$)	$\mathcal{U}(-0.4, 0.4)$	$0.00053^{+0.00041}_{-0.00046}$

Parameter	Prior	Posterior
$e_c\cos\omega_c$	$\mathcal{U}(-0.4,0.4)$a	$0.00024^{+0.00084}_{-0.00074}$
$e_c\sin\omega_c$	$\mathcal{U}(-0.4,0.4)$a	$-0.0006^{+0.0010}_{-0.0011}$
$e_c$	Deriveda	$0.0012^{+0.0005}_{-0.0011}$
$M_c\ (M_\odot)$	$\mathcal{LU}(10^{-6},10^{-3})$	$0.000538^{+0.000034}_{-0.000038}$
$M_c\ (M_{\mathrm{Jup}})$	Derived	$0.564^{+0.036}_{-0.040}$
$R_c/R_\star$	$\mathcal{N}(0.1158,0.0005)$	$0.11575^{+0.00050}_{-0.00051}$
$R_c\ (R_{\mathrm{Jup}})$	Derived	$1.272^{+0.038}_{-0.028}$
$b_c$	$\mathcal{U}(0,1)$	$0.287^{+0.066}_{-0.056}$
$a_c/R_\star$	Derived	$11.46^{+0.23}_{-0.22}$
$\Delta\Omega_{cb}$ (deg)	$\mathcal{U}(-180,180)$	$-7.1^{+2.7}_{-2.9}$
$i_{\mathrm{mut}}$ (deg)	Derived	$7.5^{+2.3}_{-2.8}$
Stellar Parameters
$M_\star\ (M_\odot)$	$\mathcal{N}(1.049,0.09)$	$1.130^{+0.075}_{-0.074}$
$R_\star\ (R_\odot)$	$\mathcal{N}(1.184,0.045)$	$1.130^{+0.030}_{-0.027}$
$u_{1,\mathrm{TESS}}$	$\mathcal{U}(0,1)$b	$0.358^{+0.088}_{-0.069}$
$u_{2,\mathrm{TESS}}$	$\mathcal{U}(0,1)$b	$-0.05^{+0.10}_{-0.17}$
TTV Jitter
$\sigma_{\mathrm{TTV}}$ (days)	$\mathcal{LU}(10^{-5},10^{-1})$	$0.00139^{+0.00017}_{-0.00022}$
TESS Sector 18 Photometry
$F_{0,\mathrm{S18}}$	$\mathcal{U}(-0.1,0.1)$	$0.00170^{+0.00064}_{-0.00059}$
$\sigma_{\mathrm{GP,S18}}$ (ppt)	$\mathcal{LU}(0.1,50)$	$9.3^{+2.7}_{-4.3}$
$P_{\mathrm{GP,S18}}$ (days)	$\mathcal{U}(4,20)$	$5.60^{+0.12}_{-0.15}$
$\log Q_{0,\mathrm{S18}}$	$\mathcal{U}(0,5)$	$2.89^{+0.74}_{-0.85}$
$\log\Delta Q_{\mathrm{S18}}$	$\mathcal{U}(-10,10)$	$-4.0^{+3.3}_{-4.8}$
$f_{\mathrm{GP,S18}}$	$\mathcal{U}(0,1)$	$0.035^{+0.019}_{-0.032}$
$\sigma_{\mathrm{jit,S18}}$ (ppt)	$\mathcal{LU}(10^{-3},10)$	$0.015^{+0.024}_{-0.014}$
$D_{\mathrm{S18}}$	$\mathcal{U}(0.5,1.5)$	$1.060^{+0.017}_{-0.018}$
TESS Sector 58 Photometry
$F_{0,\mathrm{S58}}$	$\mathcal{U}(-0.1,0.1)$	$-0.00192^{+0.00031}_{-0.00034}$
$\sigma_{\mathrm{GP,S58}}$ (ppt)	$\mathcal{LU}(0.1,50)$	$6.6^{+1.5}_{-2.2}$
$P_{\mathrm{GP,S58}}$ (days)	$\mathcal{U}(4,20)$	$5.367^{+0.093}_{-0.089}$
$\log Q_{0,\mathrm{S58}}$	$\mathcal{U}(0,5)$	$3.08^{+0.66}_{-0.66}$
$\log\Delta Q_{\mathrm{S58}}$	$\mathcal{U}(-10,10)$	$-3.6^{+3.1}_{-5.6}$
$f_{\mathrm{GP,S58}}$	$\mathcal{U}(0,1)$	$0.17^{+0.08}_{-0.11}$
$\sigma_{\mathrm{jit,S58}}$ (ppt)	$\mathcal{LU}(10^{-3},10)$	$0.010^{+0.012}_{-0.009}$
$D_{\mathrm{S58}}$	$\mathcal{U}(0.5,1.5)$	$1.001^{+0.017}_{-0.014}$
TESS Sector 85 Photometry
$F_{0,\mathrm{S85}}$	$\mathcal{U}(-0.1,0.1)$	$0.00028^{+0.00081}_{-0.00093}$
$\sigma_{\mathrm{GP,S85}}$ (ppt)	$\mathcal{LU}(0.1,50)$	$8.2^{+1.2}_{-1.7}$
$P_{\mathrm{GP,S85}}$ (days)	$\mathcal{U}(4,20)$	$5.55^{+0.22}_{-0.20}$
$\log Q_{0,\mathrm{S85}}$	$\mathcal{U}(0,5)$	$0.099^{+0.064}_{-0.099}$
$\log\Delta Q_{\mathrm{S85}}$	$\mathcal{U}(-10,10)$	$6.8^{+3.0}_{-1.4}$
$f_{\mathrm{GP,S85}}$	$\mathcal{U}(0,1)$	$0.73^{+0.27}_{-0.14}$
$\sigma_{\mathrm{jit,S85}}$ (ppt)	$\mathcal{LU}(10^{-3},10)$	$0.008^{+0.011}_{-0.007}$
$D_{\mathrm{S85}}$	$\mathcal{U}(0.5,1.5)$	$1.052^{+0.018}_{-0.020}$
RV Instrumental Parameters
Parameter	Prior	Posterior
$\sigma_{\mathrm{jit},\mathrm{MX-Blue-1}}$ (m s$^{-1}$)	$\mathcal{LU}(0.1,100)$	$7.6^{+3.0}_{-7.5}$
$\gamma_{\mathrm{MX-Blue-1}}$ (m s$^{-1}$)	$\mathcal{U}(-100,100)$	$12.4^{+3.5}_{-2.3}$
$\sigma_{\mathrm{jit},\mathrm{MX-Blue-2}}$ (m s$^{-1}$)	$\mathcal{LU}(0.1,100)$	$1.8^{+3.5}_{-1.7}$
$\gamma_{\mathrm{MX-Blue-2}}$ (m s$^{-1}$)	$\mathcal{U}(-100,100)$	$35.0^{+3.7}_{-3.2}$
$\sigma_{\mathrm{jit},\mathrm{MX-Red-1}}$ (m s$^{-1}$)	$\mathcal{LU}(0.1,100)$	$1.9^{+2.9}_{-1.8}$
$\gamma_{\mathrm{MX-Red-1}}$ (m s$^{-1}$)	$\mathcal{U}(-100,100)$	$10.0^{+3.4}_{-2.4}$
$\sigma_{\mathrm{jit},\mathrm{MX-Red-2}}$ (m s$^{-1}$)	$\mathcal{LU}(0.1,100)$	$1.7^{+2.6}_{-1.6}$
$\gamma_{\mathrm{MX-Red-2}}$ (m s$^{-1}$)	$\mathcal{U}(-100,100)$	$15.7^{+3.8}_{-3.3}$
$\sigma_{\mathrm{jit},\mathrm{NEID-3}}$ (m s$^{-1}$)	$\mathcal{LU}(0.1,100)$	$11^{+10}_{-10}$
$\gamma_{\mathrm{NEID-3}}$ (m s$^{-1}$)	$\mathcal{U}(-12000,-11900)$	$-11954.1^{+1.5}_{-1.1}$
RV Activity GP
$\sigma_{\mathrm{RV,GP}}$ (m s$^{-1}$)	$\mathcal{U}(0.05,50)$	$30.7^{+5.8}_{-6.7}$
$\log Q_{0,\mathrm{RV}}$	$\mathcal{U}(0,20)$	$2.0^{+1.2}_{-1.4}$
$\log\Delta Q_{\mathrm{RV}}$	$\mathcal{U}(-10,20)$	$-2.9^{+5.0}_{-3.8}$
$f_{\mathrm{GP,RV}}$	$\mathcal{U}(0,1)$	$0.36^{+0.17}_{-0.30}$

Note—$\mathcal{U}(a,b)$, $\mathcal{LU}(a,b)$ and $\mathcal{N}(\mu,\sigma)$ denote uniform, log-uniform and Gaussian priors. The posterior column reports the median and 68% highest-density intervals for each parameter. All orbital elements are osculating elements at the beginning of the integration, taken to be 0.1 days before the first observed transit of TOI-3686 c, $T_{\mathrm{ref}} = 2458792.557987$.

Dynamics & Formation of TOI-3686

We investigate the resonant dynamics of the two planets in the TOI-3686 system. The top panel of Figure 4 shows the time evolution of the mixed resonant angle, $\psi= 2\lambda_c - \lambda_b - \hat{\varpi}_{bc}$. Here, $\lambda_b$ and $\lambda_c$ are the planets’ mean longitudes, while $\hat{\varpi}_{bc} = \operatorname{Arg}(f e_b e^{i\varpi_b} + g e_c e^{i\varpi_c})$ is the generalized longitude of pericenter (e.g., S. Hadden [55]), and $f$, $g$ are functions of Laplace coefficients (e.g., Table 3 in Y. Lithwick et al. 2012). We find that $\psi$ librates around 180°, as expected for a system in an MMR, with an amplitude $A_\psi\sim20^\circ$ and a period $P_\psi\sim260$ days; it is the libration of this angle that drives the primary TTV signal for systems in resonance.

We transformed the orbital elements to the 1-dimensional resonance model derived by D. Nesvorný & D. Vokrouhlický [102], with the following Hamiltonian:

Here, $\Psi$ is the action conjugate to the mixed resonant angle $\psi$ and is a function of the planets’ masses, eccentricities, and semimajor axes. The parameter $\delta$ denotes the system’s proximity to resonance and differs from the distance from period commensurability $\Delta$ (Equation (4)), as $\delta$ depends on the planets’ period ratios, masses, and eccentricities. Resonance requires $\delta\gtrsim0.945$, as no separatrix exists for lower values of $\delta$. For the TOI-3686 system, we find $\delta= 3.64$, indicating that the system is deep in resonance (lower panels of Figure 4). Capture into MMR is a natural outcome of migration within a protoplanetary disk ([100, 82, 75]), thus one plausible formation history of the TOI-3686 planets is that they arrived at their present locations after migrating through the disk.

This resonant configuration differs from other known hot Jupiter systems, even those with nearby transiting companions (Figure 5). While the TOI-1130 ([65, 76, 17]) and TOI-1408 ([77]) systems contain hot Jupiters and inner companions with near-integer period ratios and large TTVs, both have $\delta< 0.945$ and are formally outside the resonant regime. Instead, this configuration calls to mind the TOI-216 system (R. I. Dawson et al. 2019, 2021; D. Nesvorný et al. 2022), which contains an inner Neptune-mass super-puff ($P_b = 17.2$ days) and an outer warm Jupiter ($P_c = 34.5$ days) that are also in a resonant state ($\delta= 2.1$). Interestingly, while being at longer orbital periods, the TOI-216 planets are similar to TOI-3686 in mass ($M_b = 18 M_{\oplus}$, $M_c = 0.53 M_{\mathrm{J}}$; B. J. McKee & B. T. Montet [95]) and eccentricity ($e_b = 0.16$, $e_c \sim0$). Capture into MMR requires that the two planets undergo convergent migration, such that their period ratio approaches integer commensurability from above (M. H. Lee & S. J. Peale [82]; K. Batygin [10]). However, the outer planet TOI-3686 c is likely massive enough to have opened a gap in the disk leading to slower type II migration, so the two planets may have undergone divergent migration, depending on disk properties (K. D. Kanagawa & E. Szuszkiewicz [72]). In the case of TOI-216, D. Nesvorný et al. [101] argued that the $P_{\mathrm{orb}} = 17$ day inner Neptune was trapped in the transition region between the disk’s inner magnetorotational instability (MRI)-active zone and the outer dead zone, allowing the outer planet’s migration to catch up and capture the planet pair into resonance. Given the smaller or-

bital separation of TOI-3686 b ($a/R_{\star} = 7.2 \pm0.1$), its inward migration may instead have stalled at the inner disk edge. For a young, solar-mass star, the magnetospheric truncation radius is expected to be $\sim7R_{\star}$ (J. Bouvier et al. [21]), precisely where we observe TOI-3686 b today. The planets in the TOI-3686 system are thus consistent with having formed further out in the disk and migrated inward, with trapping at the disk edge enabling resonant capture and setting the innermost orbital period at a location similar to planets in other exoplanetary systems (e.g., E. J. Lee & E. Chiang [81]; K. Batygin et al. [11]).

Given its short orbital period, planet b’s moderate eccentricity ($e_b = 0.16$) is surprising as tides raised on the planet by its host star should act to circularize its orbit. Assuming the constant phase lag model of equilibrium tides, its eccentricity should be damped on a timescale of (e.g., B. Jackson et al. [68]):

This eccentricity damping should lead to resonant repulsion, with the planet’s semimajor axis shrinking at the rate of $\dot{a}/a = -2e^2/\tau_e$ (J. C. B. Papaloizou & C. Terquem [104]; Y. Lithwick & Y. Wu 2012; K. Batygin & A. Morbidelli [12]). We ran a long-term forward integration of the system using the rebound code (H. Rein & S.-F. Liu [113]), incorporating a dissipative eccentricity damping force through reboundx (D. Tamayo et al. [133]). These integrations showed that within a time of $0.3 \times\tau_e$, the orbit of TOI-3686 b would have decayed sufficiently that the system is no longer in the resonant regime, with $\delta< 0.945$.

The resonance could be preserved if the hot Jupiter’s orbit is shrinking more quickly than that of the inner planet, counteracting the resonant repulsion. Because the host star is presently spinning more slowly than the orbital period of TOI-3686 c, tides raised on the host star would shrink the planet’s orbit (F. A. Rasio et al. [111]), on a timescale of

Estimates for the stellar reduced tidal quality factor $Q'_{\star}$ vary by several orders of magnitude (e.g., S. Meibom & R. D. Mathieu [96]), but even using the relatively high dissipation rate of $Q'_{\star} = 2 \times10^5$ measured from the orbital decay of WASP-12b (S. W. Yee et al. [148]), we find an orbital decay timescale for TOI-3686 c of $\tau_{a,c} = 37$ Gyr, too slow to counterbalance the shrinking of planet b’s orbit due to resonant repulsion. In fact, stellar tides alone would also cause the inner planet’s orbit to decay faster than the outer planet, with $\tau_{a,b}/\tau_{a,c} = 0.85$. Thus, tidally-driven orbital decay cannot preserve the system’s resonant state over timescales much longer than $\tau_e$.

A more plausible explanation is that the current dynamical state of TOI-3686 is only a temporary phase of the system’s evolution, and it will eventually relax toward more circular, nonresonant orbits. As Equation (6) indicates, a planetary reduced tidal quality factor of $Q'_b \sim10^5$ implies an eccentricity damping timescale of 0.6 Gyr, comparable to the system’s gyrochronological age (§ 3.1). This value is somewhat less dissipative than the $10^2$–$10^4$ typically assumed for sub-Neptunes, but is consistent with recent estimates of this quantity—for example, M.-T. Wang et al. [140] recently suggested $Q' > 10^5$ to explain the remnant eccentricities in the TOI-1449 system, which contains two near-resonant sub-Neptunes.

If TOI-3686 is indeed young and still undergoing tidal circularization and resonant repulsion, this could explain why it is the only one of the known hot Jupiter systems with transiting companions to be in a resonant state. The other such systems all have older ages, and could simply have had more time to evolve out of resonance. Indeed, L. Borsato et al. [17] attributed the location of the TOI-1130 planets just outside the formal resonant domain to past tidal evolution. TOI-3686 may therefore be an early snapshot of an evolutionary sequence in which hot Jupiters and their inner companions are captured into resonance during disk migration and then gradually pushed out of it by tides, leaving behind the near-resonant or non-resonant pairs seen around older stars. A similar trend has been observed for compact multi-planet systems containing smaller planets, which appear to be more commonly found to be in or close to MMR around young stars, with their frequency declining over 100 Myr to Gyr timescales (J. H. Hamer & K. C. Schlaufman [57]; F. Dai et al. [36]).

A Size Dichotomy Between Hot Jupiters and their Close Friends

Figure 5 shows TOI-3686 alongside the other known systems containing a hot Jupiter ($R_p > 7 R_{\oplus}$, $P_{\mathrm{orb}} < 10$ days) and a transiting companion.⁷ An interesting pattern that can be observed from the growing number of such systems is that they exhibit a strong size dichotomy, with the companions to the hot Jupiter being significantly smaller in size. This disparity stands in contrast to the “peas-in-a-pod” pattern observed for compact multi-planet systems of super-Earths and sub-Neptunes ([98, 142]), where planets within the same system tend to have similar sizes. We performed a Kolmogorov-Smirnov test comparing the distribution of radius ratios between hot Jupiters and their companions with the same distribution for the compact multi-planet systems in L. M. Weiss et al. [142], finding that they are clearly drawn from different distributions ($p = 8 \times10^{-8}$).

The radius contrast for hot Jupiters with companions also stands apart from systems containing a longer-period warm Jupiter ($10$ days $< P_{\rm orb} < 100$ days) and another transiting planet, as seen in Figure 5. While there are systems containing one warm Jupiter and a much smaller planet, there are also many examples of systems containing warm Jupiters with additional transiting giant planets. In contrast, there are no known hot Jupiter systems containing multiple transiting giant planets.⁸ A Kolmogorov-Smirnov test comparing the radius ratios between hot and warm Jupiters with companions is marginally significant, with $p = 0.061$.

A connected observation is that companions to hot Jupiters are typically among the largest and most massive of the super-Earth/sub-Neptune population for their orbital period, with a median companion mass of $11 M_{\oplus}$, close to the $\sim10 M_{\oplus}$ threshold thought to be necessary for initiating runaway gas accretion ([107]). One possible interpretation then is that these systems may have initially formed as pairs of planets with similar masses, but only one subsequently accreted a substantial gas envelope, leading to a stark size dichotomy. Further enhancement of the size and mass contrast for systems containing the shortest-period planets may then subsequently occur due to photoevaporation ([103, 88]), as some of these lie on the lower edge of the “Neptunian desert” ([132, 13, 94, 29]). The difference between the hot and warm Jupiters with companions may be a consequence of the fact that the warm Jupiter systems are less tightly-packed, potentially allowing multiple cores to grow into giant planets.

One caveat to this interpretation is that only eight of the twelve known transiting companions to hot Jupiters have measured masses, and both the RV and TTV techniques are more sensitive to more massive planets. Detection with the transit method is also biased toward larger radii, and most of the known companions, including TOI-3686 b, lie close to TESS’ detection limits. Mass measurements for the remaining systems, as well as more sensitive transit surveys like PLATO (H. Rauer et al. [112]) that could detect smaller planets, would test whether companions to hot Jupiters are systematically larger and more massive than typical sub-Neptunes. If so, these systems may be the high-mass extremes of the processes that form other compact multi-planet systems, albeit where one planet crossed the threshold for undergoing runaway gas accretion.

Conclusions

We have presented the discovery and characterization of the TOI-3686 system. The 40-minute TTVs of TOI-3686 c are amongst the largest ever observed for a hot Jupiter, while the 10 hour TTVs of the inner sub-Neptune, have the third largest TTV amplitude of any planet measured relative to its orbital period, after TOI-1408 c (J. Korth et al. [77]) and TOI-4504 d (M. Ví tková et al. (2025); J. M. Almenara et al. [3]).

Dynamical modeling of TESS and extensive ground-based data reveal that the two planets are deep in a 2:1 MMR. TOI-3686 joins the handful of hot Jupiters known to have nearby transiting companions, and is the only such system known to be in a truly resonant state. Because tidal dissipation is expected to eventually drive the system out of resonance, that we are observing the system in this resonant configuration implies a reduced tidal quality factor for the inner sub-Neptune of $Q^\prime_{\mathrm{p}} \gtrsim10^{5}$ for its eccentricity to have been maintained over the $\lesssim1$ Gyr system age.

TESS’ all-sky transit survey and growing observation baseline is enabling it to discover new, rare examples of non-isolated hot Jupiters, by detecting the transits of or TTVs induced by the companions to hot Jupiters. Such systems trace a dynamically quiescent pathway for producing close-in giant planets, likely involving early migration within the disk, and may represent the high mass end of planet formation processes that produce compact, multi-planet systems of smaller planets. More examples of hot Jupiters with significant TTVs continue to be discovered (e.g., [77, 17 S. W. Yee et al. accepted; S. Stefanov et al. submitted]), adding to the benchmark sample of hot Jupiters that are in or close to MMRs; these will be valuable systems that may help shed light on how giant planets form and interact with their protoplanetary disks.

Acknowledgments

S.W.Y. gratefully acknowledges support from the Heising-Simons Foundation. T.I.D. acknowledges support from the Harvard Origins of Life initiative. K.A.C. acknowledges support from the TESS mission via subaward s3449 from MIT and NASA grants 80NSSC24K1889 and 80NSSC26K0081. E.N. acknowledges support from the TESS mission via NASA grants 80NSSC24K1889 and 80NSSC26K0081. Funding for K.B. was provided by the European Union (ERC AdG SUBSTELLAR, GA 101054354).

The authors thank the observers at the WIYN 3.5m and Gemini-North telescopes for their skillful execution of the observations, and the NEID and MAROON-X instrument teams for providing the reduced data products.

AI tools (ChatGPT 5.6 Sol, Claude Sonnet 5.0, Claude Opus 5.5) were used to improve the visualizations in this manuscript. No text in this manuscript was generated by LLMs. Responsibility for the final manuscript lies entirely with the authors.

This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST): doi: 10.17909/8ps9-1y38. Funding for the TESS mission is provided by the NASA’s Science Mission Directorate. We acknowledge the use of public TOI Release data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. This research has made use of the Exoplanet Follow-Up Observation Program (ExoFOP; DOI: 10.26134/ExoFOP5) website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Part of the data presented herein were obtained with the WIYN 3.5m from telescope time allocated under the NN-EXPLORE program, with support from the National Aeronautics and Space Administration and the U.S. National Science Foundation. Some of the data were obtained with the NEID spectrograph built by Penn State University and operated at the WIYN Observatory by NSF NOIRLab, under the NN-EXPLORE partnership of the National Aeronautics and Space Administration and the U.S. National Science Foundation. The authors are honored to be permitted to conduct astronomical research on I’oligam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham.

This work makes use of observations from the LC OGT network. Part of the LCOGT telescope time was granted by NOIRLab through the Mid-Scale Innovations Program (MSIP). MSIP is funded by NSF. This article is based on observations made in the Two-meter Twin Telescope (TTT)⁹ sited at the Teide Observatory of the Instituto de Astrofísica de Canarias (IAC), that Light Bridges operates in Tenerife, Canary Islands (Spain). The observation time rights (DTO) used for this research were consumed in the PEI ”SB-STLLAR25”. This research used storage and computing capacity in ASTRO POC’s EDGE computing center at Tenerife under the form of Indefeasible Computer Rights (ICR). The ICR were consumed in the PEI “SBSTL-LAR25” with the collaboration of Betchle AG.

This work is based on observations obtained at the international Gemini Observatory, a program of NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation on behalf of the Gemini Observatory partnership: the U.S. National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This work was enabled by observations made from the Gemini North telescope, located within the Maunakea Science Reserve and adjacent to the summit of Maunakea. We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance.

We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through projects PID2021-125627OB-C32 and PID2024-158486OB-C32. This work is supported by the ERC Grant (ERC Advanced Grant SPEAR, GA 101200674). Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency (ERCEA). Neither the European Union nor the granting authority can be held responsible for them.

We acknowledge the Marvel consortium (Institute of Astronomy, KU Leuven (PI); Stockholm University; University of Vienna; Instituto de Astrofísica de Canarias; Institute of Space Studies of Catalonia (IEEC), Barcelona; School of Mathematical and Physical Sciences, Macquarie University; University of Warwick; Ruprecht Karl University of Heidelberg) for the telescope time granted on the Marvel Telescope Array in La Palma, Spain.

Facilities: TESS, Gemini:Gillett (MAROON-X), WIYN (NEID), FLWO:1.5m (TRES), FLWO:1.2m (KeplerCam), LCOGT 1.0m (Sinistro)

Software: AstroImageJ ([33]), astropy ([4] (Astropy Collaboration et al. 2018; Astropy Collaboration et al. 2022)), EXOFASTv2 ([42, 43]), exoplanet ((D. Foreman-Mackey et al. 2021; D. Foreman-Mackey et al. 2021)), jaxoplanet ([59]), jupyter ([92, 91]), numpy ([16, 106]), Tapir ([71])