TOI-1670 b and c: An Inner Sub-Neptune with an Outer Warm Jupiter Unlikely to Have Originated from High-eccentricity Migration

The origin of giant planets interior to the water-ice line remains an open question. A number of theories have been proposed to explain the closest-in (P < 10 days) giant planets, or hot Jupiters (HJs; e.g., Dawson & Johnson 2018; Fortney et al. 2021). These scenarios are primarily divided between dynamically "violent" or "gentle" mechanisms. The former consists of three-body dynamical interactions such as planet–planet scattering or high-eccentricity tidal migration (e.g., Wu & Murray 2003; Fabrycky & Tremaine 2007; Triaud et al. 2010; Naoz et al. 2011; Batygin 2012). The latter refers to disk migration (e.g., Ward 1997; Albrecht et al. 2012; Kley & Nelson 2012) or in situ formation (e.g., Batygin et al. 2016; Boley et al. 2016; Huang et al. 2016; Anderson et al. 2020). These processes have also been used to explain part of the farther-out population of warm Jupiters (WJs; defined here to have 10 days < P < 200 days). However, observed WJ demographics suggest that multiple processes are present in sculpting these more distant giant systems.

WJs can be broadly divided into two classes. The first is a transient population that will likely evolve into HJs. In more disruptive formation mechanisms, such as high-eccentricity tidal migration, giant planets at comparatively wide separations are disturbed onto highly eccentric orbits by a third body via planet–planet scattering or von Zeipel–Lidov–Kozai oscillations and eventually circularize into orbits with shorter periods (Kozai 1962; Lidov 1962; Naoz 2016; Ito & Ohtsuka 2019). Eccentric giant planets undergoing this tidally damped inward migration are caught in a rapid, temporary state (Naef et al. 2001; Dawson & Johnson 2018; Dong et al. 2021; Jackson et al. 2021). They are expected to start their journeys with much higher eccentricities (e ≳ 0.9; Vick et al. 2019), which can decay as rapidly as ∼1 Myr as they settle in near their host star (Patra et al. 2020; Mancini et al. 2022). The majority of WJs are not expected to belong to this transient classification.

Instead, most WJs are a part of a "static" population that will remain stable over long time periods. This group consists of the apparently single systems with low-to-moderate eccentricities, as well as coplanar multiplanet systems containing WJs with low eccentricities. These WJs have periapses larger than what is required for efficient tidal damping of their orbits, which occurs at ≲0.05 au (Anderson et al. 2016; Dong et al. 2021), so these planets cannot be undergoing high-eccentricity migration. If most giant planets form beyond the water-ice line, other migration mechanisms must play a major role in sculpting these WJ orbital properties and demographics (e.g., Veras & Armitage 2005; Fogg & Nelson 2009; Dong et al. 2014; Ortiz et al. 2015; Huang et al. 2016; Anderson & Lai 2017; Anderson et al. 2020; Schlecker et al. 2020). However, the relative importance of these pathways is still unknown. Investigating WJ orbital eccentricities can place additional constraints on the dominant giant planet migration mechanism, since each scenario will produce different observed eccentricity distributions.

WJs have an eccentricity distribution that peaks at e = 0.0 with a tail that extends out to e ∼ 0.8 (Kipping 2013; Dong et al. 2021). In order to produce the population of WJs with moderately eccentric orbits (e ∼ 0.2–0.7), a mechanism is needed that can excite eccentricities. These potential excitation scenarios include interactions involving a disk (e.g., Goldreich & Sari 2003; Petrovich et al. 2019), secular eccentricity oscillations driven by interactions with a distant inclined giant planet (e.g., Anderson & Lai 2017), and planet–planet scattering events (e.g., Mustill et al. 2017; Frelikh et al. 2019; Marzari & Nagasawa 2019; Anderson et al. 2020). An important clue is the observed dependence on metallicity of the giant planet eccentricity distribution, where metal-rich systems (that may more favorably form multiple giant planets) are more likely to host eccentric gas giants (Dawson & Murray-Clay 2013).

Systems hosting WJs with low-mass inner companions on coplanar orbits are especially useful laboratories to test these planet formation and migration theories. Their small orbital eccentricities and low mutual inclinations suggest that disk migration or in situ formation likely helped create this population. WJs have a relatively high close companion rate of nearly 50% (Huang et al. 2016). However, their intrinsically low occurrence rate (∼1%–2%; Cumming et al. 2008) combined with the difficulty of detecting lower-mass inner planets means that only a handful of known multiplanet systems host a WJ (Johnson et al. 2010; Santerne et al. 2016; Fernandes et al. 2019). Increasing the number of systems with this multiplanet architecture may further distinguish this subsample into two WJ populations, each of which likely reflects different formation and migration routes.

Here we present the discovery of the transiting multiplanet system TOI-1670 bc, a WJ (TOI-1670 c) with an inner sub-Neptune (TOI-1670 b) found with the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015). TOI-1670 b and c were originally identified by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016) pipeline as two promising transiting signals that were subsequently promoted to TESS Object of Interest (TOI; Guerrero et al. 2021) status. TOI-1670 (TIC ID 441739020; 2MASS J17160415+7209402; Gaia DR2 1651911084230149248) is a relatively inactive (log${R}_{\mathrm{HK}}^{{\prime} }=-4.93$) old F7 dwarf with a TESS apparent magnitude of 9.5 mag and a moderate projected rotational velocity of ≈9 km s⁻¹ (see Section 3). In this work, we validate both planets and measure the mass of the outer planet, TOI-1670 c. In Section 2, we describe the TESS photometric data and follow-up radial velocity (RV) observations used in the planet validation and mass measurement. Our characterization of the system, including the host star and a global fit to the RVs and light curve, is presented in Section 3. We conclude in Section 4 by contextualizing TOI-1670 in the paradigm of WJs and their formation.

KESPRINT ³⁰ is an international collaboration focused on the discovery, confirmation, and characterization of exoplanet candidates from space-based missions (e.g., Gandolfi et al. 2018; Persson et al. 2018; Livingston et al. 2019; Lam et al. 2020; Šubjak et al. 2020). As part of this consortium, a series of ground-based follow-up observations of TOI-1670 were taken. These data are primarily used to reject the possibility of a false-positive scenario in which the observed transiting signal is caused by something other than a planet. For example, this includes a low-mass eclipsing binary (EB), a grazing transit of an EB, a background EB, or a transiting planet around a background star. Reconnaissance spectra are used to exclude an EB scenario by constraining the maximum amplitude of the RV signal. High-resolution speckle images are taken to exclude binary companions to TOI-1670 and nearby background stars. High-resolution spectra are used to characterize the host star and, when possible, measure the masses of the planets.

2.1. TESS Photometry

TOI-1670 was observed by TESS at 2 minute cadence over 11 sectors (15, 16, 18, 19, 20, 21, 22, 23, 24, 25, and 26) for a total of 323 days. Images were reduced and light curves were analyzed for transit signals with the TESS SPOC pipeline (Jenkins et al. 2016), which identified two potential transit signals (Jenkins 2002; Jenkins et al. 2010, 2020) with periods of 40.7 (TOI-1670.01) and 10.9 (TOI-1670.02) days. The SPOC vetting tests (Twicken et al. 2018; Li et al. 2019) validated both signals as consistent with planets, and they were designated as TOIs (Guerrero et al. 2021) by the TESS Science Office.

We downloaded the SPOC Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) light curve (Smith et al. 2012; Stumpe et al. 2012, 2014) from the MAST data archive ³¹ using the lightkurve (Lightkurve Collaboration et al. 2018) software package. We removed all of the photometric measurements that are flagged as poor quality by the SPOC pipeline (DQUALITY > 0) or where either the flux or flux error is listed as NaN. Outlier rejection was performed at 3σ for positive outliers and 10σ for negative outliers to allow for transit events. The light curve was flattened by removing low-frequency trends using a Savitzky–Golay filter (Savitzky & Golay 1964) after all transit events were masked out. The final light curve for TOI-1670 is shown in Figure 1. The photometric points used in the global model fit are shown in purple and pink. These cover the transit events for TOI-1670 b and c, respectively, and their times of transit are further denoted by the corresponding colored triangles along the time axis.

Zoom In Zoom Out Reset image size

2.2. TRES Reconnaissance Spectroscopy

2.3. OES Reconnaissance Spectroscopy

2.4. Tull Coudé Spectroscopy

2.5. FIES Spectroscopy

2.6. HARPS-N Spectroscopy

2.7. High-resolution Imaging

On the nights of UT 2021 April 5 and June 24, TOI-1670 was observed with the NESSI and 'Alopeke speckle imagers (Scott et al. 2018; Scott 2019), mounted on the 3.5 m WIYN telescope at Kitt Peak and the 8.1 m Gemini North telescope on Maunakea, respectively. Both instruments simultaneously acquire data in two bands centered at 562 and 832 nm using high-speed electron-multiplying CCDs. Observations of TOI-1670 were performed in the 562 and 832 nm bands following the procedures described in Howell et al. (2011). The resulting reconstructed images have a 5σ delta magnitude contrast of 4 to 8 magnitudes at angular separations from 20 mas to 12 in the 832 nm band (Figure 2). No other companion sources are detected in the reconstructed images within these angular limits down to the contrasts obtained. These angular limits correspond to spatial separations of 3.3–200 au at the distance of TOI-1670.

Zoom In Zoom Out Reset image size

3.1. Stellar Parameters

Planetary parameters measured from the global joint fit of the transit and RV data depend on precise stellar mass and radius measurements. In particular, R_p is determined from the transit depth and M_p from the RV semiamplitude, which are dependent on R_* and M_*, respectively. The stellar mass and radius can be inferred with atmospheric and evolutionary models using the spectroscopic parameters (T_eff, log g, [Fe/H], v sin i). We determine both of these spectroscopic and fundamental parameters for TOI-1670 using several approaches described below.

3.1.1. Spectral Analysis

We analyzed the coadded HARPS-N (S/N = 180) spectrum with the spectral analysis package Spectroscopy Made Easy (SME; Valenti & Piskunov 1996; Valenti & Fischer 2005; Piskunov & Valenti 2017). The spectral fitting technique of SME minimizes the χ² value by fitting synthetic spectra of stars based on grids of atmosphere models and observations. We fit the coadded HARPS-N spectrum with the ATLAS12 model spectra (Kurucz 2013) using the non-local thermodynamic equilibrium SME version 5.2.2 following the procedure described in Fridlund et al. (2017) to compute T_eff, log g, v sin i, and chemical abundances. The stellar surface gravity, log g, was estimated using the spectral wings of the Ca i 6102, 6122, 6162 Å triplet and the Ca i 6439 Å line. The microscopic and macroscopic turbulences, V_mic and V_mac, were held fixed to the values determined in the calibration for stars with similar T_eff and log g from Bruntt et al. (2010) and Doyle et al. (2014), respectively.

We also derive the stellar parameters using the publicly available SpecMatch-Emp software package (Yee et al. 2017). SpecMatch-Emp compares the HARPS-N template spectrum to a high-resolution (R ∼ 55,000), high-S/N (>100) Keck/HIRES optical spectral library of 404 well-characterized early- to late-type dwarfs (F1 to M5). The empirical spectra are calibrated using interferometry, so SpecMatch-Emp produces estimates for T_eff, [Fe/H], and R_* (instead of log g). Prior to running the code, we convert the HARPS-N spectrum template onto the Keck/HIRES format following the procedure described in Hirano et al. (2018).

The stellar parameters derived from SME and SpecMatch-Emp are in good agreement with each other (Table 1). From SME, we find an effective temperature of T_eff = 6170 ± 61 K and metallicity of [Fe/H] = 0.09 ± 0.07 dex, while SpecMatch-Emp gives T_eff = 6048 ± 110 K and [Fe/H] = 0.05 ± 0.09 dex; these are consistent with each other within 1σ. These results are also in good agreement with the photometrically derived effective temperature from Gaia DR2 (${T}_{\mathrm{eff}}={6162}_{-175}^{+162}$ K) and agree at the 2σ level with the TESS Input Catalog (TIC) v8 (Guerrero et al. 2021) value of 6345 ± 121 K. For this work, we adopt the spectroscopic parameters from SME as it produces all atmospheric parameters. The final adopted stellar parameters are reported in Table 3.

Table 1. Spectroscopic Parameters of TOI-1670 Derived Using SME and SpecMatch-Emp

Method	Teff (K)	log g (g cm−3)	[Fe/H] (dex)	v sin i (km s−1)
SME	6170 ± 61	4.29 ± 0.11	0.09 ± 0.07	9.2 ± 0.6
SpecMatch-Emp	6048 ± 110	⋯	0.05 ± 0.09	⋯

Download table as:  ASCIITypeset image

3.1.2. Stellar Mass and Radius

We infer the stellar radius by fitting the spectral energy distribution (SED) of TOI-1670 using the software package ARIADNE. ³² ARIADNE utilizes a Bayesian model averaging framework that convolves four stellar atmosphere models—Phoenix v2 (Husser et al. 2013), BT-Settl (Allard et al. 2011), Kurucz (1993), and Castelli & Kurucz (2003)—with the response functions of commonly available broadband filters. For our SED fitting, we use the Two Micron All Sky Survey (2MASS) JHK_s , Gaia DR2 (G, G_BP, B_RP), Johnson V and B, and Wide-field Infrared Survey Explorer (WISE; W1 and W2) bandpasses. Synthetic SEDs are created by interpolating in T_eff–$\mathrm{log}g$–[Fe/H] space. Distance, radius, A_V , and excess photometric uncertainty terms are free parameters in the fitting process. We set the priors for T_eff, log g, and [Fe/H] to the values we found in Section 3.1.1, the distance prior to the Bailer-Jones et al. (2021) Bayesian-based value (${165.72}_{-0.38}^{+0.32}$ pc), and the stellar radius prior to the Gaia DR2 value (${R}_{* }={1.38}_{-0.07}^{+0.08}$ R_⊙). The extinction, A_V , has a flat prior limited by the maximum line-of-sight reddening according to the recalibrated SFD galaxy dust map (Schlegel et al. 1998; Schlafly & Finkbeiner 2011). The excess photometric noise parameters all have Gaussian priors centered at zero with a standard deviation equal to 10 times the reported photometric error. The SED of TOI-1670 and best-fitting model are shown in Figure 3.

Zoom In Zoom Out Reset image size

We estimate the mass of TOI-1670 using the stellar isochrone software package isochrones (Morton 2015a) and the MESA Isochrones and Stellar Tracks (Dotter 2016; Choi et al. 2016) evolutionary model grids. The isochrones package infers fundamental stellar parameters by comparing a variety of observational inputs to interpolated model values. We input the Gaia DR2 parallax, broadband photometry (2MASS JHK_s ; Gaia DR2 G, G_BP, and B_RP; Johnson V and B; and WISE W1 and W2), and the SME spectroscopic values (T_eff, log g, and [Fe/H]) as priors. The posteriors are sampled using the MultiNest (Feroz et al. 2009, 2019) sampling algorithm.

All values for the stellar radius and mass are reported in Table 2. We include values from the TIC and Gaia DR2, as well as the typical mass and radius for an F7V dwarf for reference (Cox 2000). The SpecMatch-Emp fit also derives a stellar radius, which we couple to the calibration equations from Torres et al. (2010) to infer a surface gravity of log g = 4.14 ± 0.07 dex and a stellar mass of 1.25 ± 0.09 M_⊙. All values are in good agreement with each other. We adopt the ARIADNE radius (R_* = 1.316 ± 0.019 R_⊙) and the isochrones mass (M_* = 1.21 ± 0.02 M_⊙) as the stellar parameters to be used in the global fit and report all adopted physical, photometric, and kinematic properties of TOI-1670 in Table 3.

Table 2. Stellar Mass and Radius of TOI-1670 Derived from Different Methods

Method	M* (M⊙)	R* (R⊙)
ARIADNE a	1.16 ± 0.16	1.316 ± 0.019
isochrones	1.21 ± 0.02	1.316 ± 0.007
SpecMatch-Emp+Torres b	1.25 ± 0.09	1.57 ± 0.18
TIC c	1.25 ± 0.18	1.312 ± 0.057
Gaia DR2 d	⋯	${1.38}_{-0.07}^{+0.08}$
Typical F7V dwarf e	1.21	1.32
Adopted	1.21 ± 0.02	1.316 ± 0.019

Notes.

^a Mass calculated using derived radius and log g. ^b Mass calculated using SpecMatch-Emp parameters and calibration equations from Torres et al. (2010). ^c Stassun et al. (2019). ^d Gaia Collaboration et al. (2018). ^e Cox (2000).

Download table as:  ASCIITypeset image

Table 3. Adopted Physical, Photometric, and Kinematic Properties of TOI-1670

Parameter	Value	Source
TIC ID	441739020	1
TOI ID	1670	1
Gaia ID	1651911084230149248	2
2MASS ID	J17160415 + 7209402	3
Gaia α (J2000.0)	17:16:04.16	2
Gaia δ (J2000.0)	+72:09:40.17	2
Gaia epoch	2015.5	2
Gaia parallax (mas)	5.92 ± 0.02	2
Distance (pc)	${165.72}_{-0.38}^{+0.32}$	4
Gaia μα cos δ (mas yr−1)	−6.09 ± 0.05	2
Gaia μδ (mas yr−1)	5.154 ± 0.05	2
B (mag)	10.43 ± 0.03	5
V (mag)	9.89 ± 0.03	5
T (mag)	9.423 ± 0.0061	1
G (mag)	9.8232 ± 0.0004	2
GRP (mag)	9.4145 ± 0.0014	2
GBP (mag)	10.0747 ± 0.0010	2
J (mag)	8.97 ± 0.02	3
H (mag)	8.75 ± 0.03	3
Ks (mag)	8.724 ± 0.02	3
W1 (mag)	8.689 ± 0.023	6
W2 (mag)	8.702 ± 0.020	6
Teff (K)	6170 ± 61	This work
log g (g cm−3)	4.29 ± 0.11	This work
[Fe/H] (dex)	0.09 ± 0.07	This work
v sin i (km s−1)	9.2 ± 0.6	This work
M* (M⊙)	1.21 ± 0.02	This work
R* (R⊙)	1.316 ± 0.019	This work
ρ* (g cm−3)	0.752 ± 0.036	This work
Age (Gyr)	2.53 ± 0.43	This work
AV (mag)	0.010 ± 0.006	This work

References. (1) Stassun et al. (2019), (2) Gaia Collaboration et al. (2018), (3) Cutri et al. (2003), (4) Bailer-Jones et al. (2021), (5) Høg et al. (2000), (6) Cutri et al. (2021).

Download table as:  ASCIITypeset image

3.2. Stellar Activity

Stellar activity in the form of rotationally modulated starspots and granulation can both mimic and mask the signals of planets in light curves (Llama & Shkolnik 2015, 2016) and RVs (e.g., Figueira et al. 2013). Thus, prior to running the global model fit, we first examine whether stellar activity significantly influences the light curve and RV time series of TOI-1670. We measured a low average value of log${R}_{\mathrm{HK}}^{{\prime} }=-4.93\pm 0.01$ from the HARPS-N spectra, which suggests that TOI-1670 is a quiet star not dominated by stellar activity (Mamajek & Hillenbrand 2008). The TESS light curve prior to detrending also does not exhibit any significant rotation or activity-induced variability.

A common statistical tool used to detect periodic signals in unevenly sampled time series data is the Lomb–Scargle periodogram (Lomb 1976; Scargle 1982). We utilize this algorithm to search for periodicity in both the TESS photometry and RV activity indicators to distinguish stellar activity–based signals from those induced by planetary motion. We compute the generalized Lomb–Scargle (GLS) periodograms (Zechmeister & Kürster 2009) for the "undetrended" PDCSAP light curve with the transit events removed, the RVs, the aforementioned activity indicators, and the spectral window function over the frequency range 0.0005–0.5 day⁻¹ (2–2000 days) in Figure 4. The GLS power thresholds corresponding to false-alarm probability (FAP) levels of 1% and 0.1% computed via a bootstrap approach are shown as blue dotted lines (Kuerster et al. 1997). The GLS periodogram for the RVs was computed for the combined HARPS-N and Tull Coudé data after subtracting the systematic velocity offsets as reported in Table 4. The periodogram of the TESS photometry has very low power with no peaks that have significance higher than the 1% FAP level. This is consistent with the flat nature of the undetrended PDCSAP light curve and indicates that TOI-1670 does not have a large starspot coverage fraction. The strongest signal in the periodogram of the RVs is at the ∼40.7 day orbital period of TOI-1670 c, which has an FAP < 0.1%. This peak has no counterparts in the periodograms of the activity indices, which would be the case if that signal originated from stellar activity. Activity signals can also appear at the frequency of the stellar rotation period. Using the stellar radius and v sin i, we can place a lower limit of P_rot ≥ 7.2 days (f ≤ 0.138 day⁻¹). No significant peaks are visible in the GLS periodogram of the S-indices at this frequency.

Zoom In Zoom Out Reset image size

Table 4. Priors and Posteriors on the Global System Parameters of TOI-1670 b and c

Parameter	Adopted Prior		Posterior Values
Fitted Parameters	b	c	b	c
T0 (BJDTDB−2,457,000)	${ \mathcal U }(1721.92,1721.99)$	${ \mathcal U }(1750.82,1750.92)$	${1721.9423}_{-0.0062}^{+0.0071}$	${1750.88286}_{-0.00083}^{+0.00085}$
P (days)	${ \mathcal U }(10.980,10.988)$	${ \mathcal U }(40.7485,40.7505)$	${10.98462}_{-0.00051}^{+0.00046}$	${40.74976}_{-0.00021}^{+0.000022}$
K (m s−1)	${ \mathcal U }(0.0,10.0)$	${ \mathcal U }(10.0,100)$	${4.6}_{-3.0}^{+3.3}$	${32.7}_{-4.3}^{+4.7}$
b	${ \mathcal U }(0.0,1.0)$	${ \mathcal U }(0.0,1.0)$	${0.61}_{-0.37}^{+0.22}$	${0.76}_{-0.04}^{+0.02}$
a/R*	${ \mathcal U }(1.1,20.0)$	${ \mathcal U }(1.1,50.0)$	${16.88}_{-0.27}^{+0.27}$	${40.68}_{-0.66}^{+0.66}$
Rp/R*	${ \mathcal U }(0.0,0.05)$	${ \mathcal U }(0.0,0.15)$	${0.014}_{-0.001}^{+0.001}$	${0.077}_{-0.002}^{+0.002}$
$\sqrt{e}$ sin ω	${ \mathcal U }(-1.0,1.0)$	${ \mathcal U }(-1.0,1.0)$	${0.18}_{-0.44}^{+0.36}$	${0.27}_{-0.10}^{+0.08}$
$\sqrt{e}$ cos ω	${ \mathcal U }(-1.0,1.0)$	${ \mathcal U }(-1.0,1.0)$	$-{0.63}_{-0.21}^{+0.62}$	$-{0.07}_{-0.13}^{+0.14}$
Derived Parameters
Mp	⋯	⋯	${13.8}_{-8.7}^{+9.5}$ M⊕	${0.63}_{-0.08}^{+0.09}$ MJup
Rp	⋯	⋯	${2.06}_{-0.15}^{+0.19}$ R⊕	${0.987}_{-0.025}^{+0.025}$ RJup
ρp (g cm−3)	⋯	⋯	${8.6}_{-5.6}^{+6.9}$	${0.81}_{-0.11}^{+0.13}$
e	⋯	⋯	${0.59}_{-0.26}^{+0.17}$	${0.09}_{-0.04}^{+0.05}$
ω (deg)	⋯	⋯	${163.6}_{-53.7}^{+41.7}$	${105.5}_{-29.4}^{+28.6}$
i (deg)	⋯	⋯	${86.87}_{-1.07}^{+1.16}$	${88.84}_{-0.04}^{+0.04}$
a (au)	⋯	⋯	${0.103}_{-0.002}^{+0.002}$	${0.249}_{-0.005}^{+0.005}$
T14 (hr)	⋯	⋯	${2.80}_{-0.19}^{+0.16}$	${5.40}_{-0.05}^{+0.06}$
Teq (K)	⋯	⋯	${1062}_{-13}^{+14}$	${684}_{-9}^{+9}$
Additional Parameters
M* (from scaled parameters) (M⊙)	${ \mathcal N }(1.215,0.023)$		${1.218}_{-0.076}^{+0.081}$	${1.239}_{-0.079}^{+0.084}$
ρ* (from transit) (g cm−3)	⋯		${0.754}_{-0.035}^{+0.037}$	${0.767}_{-0.037}^{+0.038}$
q1	${ \mathcal U }(0.0,1.0)$		${0.35}_{-0.11}^{+0.19}$
q2	${ \mathcal U }(0.0,1.0)$		${0.32}_{-0.23}^{+0.39}$
γTull (km s−1)	${ \mathcal U }(-3.6849,-3.3299)$		$-{3.5132}_{-0.0043}^{+0.0042}$
γFIES (km s−1)	${ \mathcal U }(-0.1217,0.1139)$		$-{0.0108}_{-0.0061}^{+0.0061}$
γHARPS−N (km s−1)	${ \mathcal U }(-11.6134,-11.3602)$		$-{11.4805}_{-0.0027}^{+0.0025}$
RV jitter (Tull) (m s−1)	⋯		${22.5}_{-3.6}^{+3.9}$
RV jitter (FIES) (m s−1)	⋯		${5.3}_{-3.8}^{+7.2}$
RV jitter (HARPS-N) (m s−1)	⋯		${4.3}_{-2.8}^{+4.5}$

Download table as:  ASCIITypeset image

3.3. Statistical Validation of TOI-1670 b

3.4. Joint Modeling of RVs and Photometry

We perform a multiplanet global fit to the available RV and transit observations of TOI-1670 using the pyaneti modeling suite (Barragán et al. 2019). As the planetary Doppler signals are not recovered at a significant level in the TRES and Ondřejov spectra, we only use the 49 Tull Coudé and eight HARPS-N RVs in the modeling. We limit the light-curve data to photometry spanning four full transit durations before and after all transit events of TOI-1670 b and c to improve computation efficiency; this results in a total of 24,772 photometric points. These regions are shown in purple and pink in Figure 1 for TOI-1670 b and c, respectively.

We simultaneously fit the Keplerian orbit and TESS light curve for eight parameters: orbital period (P), central time of transit (T₀), RV semiamplitude (K), transit impact parameter (b), planetary-to-stellar radius (R_p/R_*), scaled semimajor axis (a/R_*), and parameterized forms of eccentricity and argument of periastron ($\sqrt{e}\ \sin \ \omega$ and $\sqrt{e}\ \cos \ \omega$). This last parameterization by Anderson et al. (2011) is used because the eccentricity posterior distribution for orbits with low e and broad ω is poorly sampled by Markov chains (e.g., Lucy & Sweeney 1971; Ford 2006; Wang & Ford 2011). By defining e and ω in a polar form, we avoid truncating the posterior distribution at zero and impose a uniform prior on e. We also adopt the parameterization of b as defined by Winn (2010),

$$\begin{eqnarray}&&b=\displaystyle \frac{a\ \cos {i}_{* }}{{R}_{* }}\left(\displaystyle \frac{1-{e}^{2}}{1+e\sin {\omega }_{* }}\right),\end{eqnarray} \tag{ 1 }$$

where i_* is the stellar inclination, in order to impose priors that exclude nontransiting orbits $(b\gt 1+\tfrac{{R}_{p}}{{R}_{* }})$.

We set narrow uniform priors on both orbital period and time of transit based on visual inspection of the light curve and the SPOC preliminary parameters. For the inner sub-Neptune, the ranges are ${T}_{0,b}={ \mathcal U }(1721.92,1721.99)$ in units of (BJD_TDB – 2,457,000) days, ${P}_{b}={ \mathcal U }(10.980,10.988)$ days, and ${K}_{b}={ \mathcal U }(0.0,10.0)$ m s⁻¹. For the outer Jupiter, the ranges are ${T}_{0,c}={ \mathcal U }(1750.82,1750.92)$ in units of (BJD_TDB – 2,457,000) days, ${P}_{c}={ \mathcal U }(40.7485,40.7505)$ days, and ${K}_{c}={ \mathcal U }(10.0,100.0)$ m s⁻¹. The stellar mass and radius are also free parameters with Gaussian priors of ${M}_{* }={ \mathcal N }(1.215,0.023)\ {M}_{\odot }$ and ${R}_{* }={ \mathcal N }(1.316,0.019)\ {R}_{\odot }$. ³⁴ These parameters are further constrained by the stellar mean density, which is affected by P and a/R_* (Seager & Mallén-Ornelas 2003; Winn 2010). We assumed a quadratic limb-darkening law following the equations from Mandel & Agol (2002), who defined the linear and quadratic coefficients as u₁ and u₂, respectively. The parameterization of ${q}_{1}={({u}_{1}+{u}_{2})}^{2}$ and ${q}_{2}=0.5{u}_{1}{({u}_{1}+{u}_{2})}^{-1}$ from Kipping (2013) is adopted. We set broad uniform priors for all other parameters and report them in Table 4. A "jitter" term is added to the RVs to account for any systematic and astrophysical variance not reported in the observational uncertainties. ³⁵

Posterior distributions of fitted and derived parameters were sampled using a Markov Chain Monte Carlo Metropolis–Hasting algorithm following the description by Sharma (2017) as implemented by pyaneti. The distributions were sampled using 50 chains for 10,000 iterations with a thinning factor of 10. The convergence of each chain was determined with the Gelman–Rubin diagnostic test (Gelman & Rubin 1992).

Using the TESS photometry and RVs, we jointly model the transits of TOI-1670 b and c and the RV curve of TOI-1670 c using the priors as previously described. ³⁶ The posterior values of the fitted and derived system parameters from pyaneti for TOI-1670 are given in Table 4. The best-fitting phased TESS light curves for TOI-1670 b and c and RV model for TOI-1670 c are plotted in Figures 5–7. Figure C1 in Appendix C displays the posterior distributions for the fitted parameters. We find a mass, radius, and density of TOI-1670 c of ${M}_{{\rm{c}}}={0.63}_{-0.08}^{+0.09}$ M_Jup, ${R}_{{\rm{c}}}={0.987}_{-0.025}^{+0.025}$ R_Jup, and ${\rho }_{{\rm{c}}}={0.81}_{-0.11}^{+0.13}$ g cm⁻³, respectively. For TOI-1670 b, we find a radius of ${2.06}_{-0.15}^{+0.19}$ R_⊕ and a 3σ mass upper limit of M_b < 0.13 M_Jup.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The existence of WJ systems hosting one or more smaller coplanar inner companions such as TOI-1670 is inconsistent with dynamical migration routes. During high-eccentricity tidal migration, multiple close-in planets would likely interact with each other and potentially lead to ejections or collisions (e.g., Rasio & Ford 1996; Chatterjee et al. 2008; Mustill et al. 2015). Similarly, planet–planet scattering and von Zeipel–Lidov–Kozai interactions require an outer companion (Veras & Armitage 2005; Anderson & Lai 2017). Multiplanet systems hosting a WJ with low eccentricity represent another type of system that experienced comparatively gentle dynamical histories, such as inward disk migration or in situ formation.

TOI-1670 joins a handful of confirmed systems with an outer warm giant exoplanet (M_p > 0.25 M_Jup, 10 days < P < 200 days) and at least one inner smaller companion (see Table 5). Figure 8 shows the eccentricity versus semimajor axis of the confirmed WJs. The WJs in all 11 systems (including TOI-1670) with similar configurations have low eccentricities, whereas the eccentricities of other WJs without known inner companions are widely distributed. This divergence further suggests that this group of TOI-1670-like systems may have formed and migrated along a similar evolutionary pathway.

Zoom In Zoom Out Reset image size

Table 5. Confirmed Systems with an Outer Warm Giant Planet and at Least One Smaller Inner Companion

System	Min	Rin	Pin	Mout	Rout	Pout	eout	σe b	Pout/Pin	Reference
	(MJup)	(RJup)	(days)	(MJup)	(RJup)	(days)
Kepler-89	<0.05	0.385	10.42	0.33	1.005	22.34	0.022	0.038	2.14	1, 2
TOI-216	0.06	0.714	17.16	0.56	0.901	34.53	0.0046	${}_{-0.0012}^{+0.0027}$	2.01	3, 4
Kepler-117	0.09	0.719	18.80	1.84	1.101	50.79	0.0323	0.0033	2.70	5, 6
Kepler-30	0.03	0.348	29.22	1.69	1.097	60.32	0.011	0.001	2.06	7, 8
HIP 57274	0.04 b	⋯	8.14	0.41 b	⋯	32.03	0.05	0.03	3.94	9
GJ 876	0.76 c	⋯	30.13	2.39 c	⋯	61.08	0.027	0.002	2.03	10
K2-290	<0.07	0.273	9.21	0.77	1.006	48.37	0 (fixed)	<0.241	5.25	11
Kepler-56	0.07	0.581	10.50	0.57	0.874	21.40	0.00	0.01	2.04	12, 13
Kepler-88	0.03	0.307	10.92	0.67	⋯	22.26	0.0572	0.0005	2.04	14
Kepler-289	0.01	0.239	66.06	0.42	1.034	125.85	0.005	0.015	1.91	15
TOI-1670	<0.13	0.184	10.98	0.63	0.987	40.75	0.09	${}_{-0.04}^{+0.05}$	3.71	This work

Notes.

^a 1σ uncertainties or 3σ upper limit on e_out. ^b Minimum mass, M_p sin i. ^c Mass determined assuming coplanar model with fixed inclinations.

References. (1) Hirano et al. (2012), (2) Weiss et al. (2013), (3) Dawson et al. (2019), (4) Dawson et al. (2021), (5) Rowe et al. (2014), (6) Bruno et al. (2015), (7) Sanchis-Ojeda et al. (2012), (8) Panichi et al. (2018), (9) Fischer et al. (2012), (10) Trifonov et al. (2018), (11) Hjorth et al. (2019), (12) Huber et al. (2013), (13) Otor et al. (2016), (14) Weiss et al. (2020), (15) Schmitt et al. (2014).

Download table as:  ASCIITypeset image

One way to disentangle whether disk migration or in situ formation plays the dominant role in sculpting these multiplanet systems is by examining their period ratios in search of near mean-motion resonances (MMRs). Disk migration is expected to efficiently capture giant planets into MMRs close to small integer period ratios such as 2:1, 3:1, 3:2, and 4:3 (e.g., Goldreich & Tremaine 1980; Lee & Peale 2001; Armitage 2010; Winn & Fabrycky 2015). In situ formation can also create planets in orbital resonances, either coincidentally or by eccentricity damping via interactions with the protoplanetary or planetesimal disk (Dawson et al. 2016; Morrison et al. 2020). In this formation scenario, there should be a population of systems that congregate at or near these different integer ratios.

Within the sample of 11 known systems that have a giant planet with a small inner companion, eight WJs (73%) are in or near a 2:1 or 3:1 resonance with the inner planet (Table 5). With a period ratio of 3.7, TOI-1670 joins two other systems, K2-290 and HIP 57274, that have non-MMR orbital period ratios greater than 3. The planets in these systems may have formed in situ or migrated inward together. Alternatively, that the planets in these systems are not locked in an MMR could also indicate that they formed independently and did not migrate together or became unstable over time (Petit et al. 2020; Pichierri & Morbidelli 2020; Izidoro et al. 2021). This may hint at a division within this small class of WJs in multiplanet systems in which some migrate into place via disk migration (those with integer period ratios), while others formed where we see them today or experienced further dynamical interaction later in their lifetime. This hypothesis can be further investigated by increasing the number of warm giant planets with smaller inner companions and examining population trends within this sample.

We thank Benjamin Tofflemire, Daniel Krolikowski, Michael Gully-Santiago, and Erik Petigura for insightful discussions on the generalized Lomb–Scargle periodogram, gas giant occurrence rate, and light-curve analysis.

Q.H.T. and B.P.B. acknowledge support from a NASA FINESST grant (80NSSC20K1554). This work benefited from involvement in ExoExplorers, which is sponsored by the Exoplanets Program Analysis Group (ExoPAG) and NASA's Exoplanet Exploration Program Office (ExEP). B.P.B. acknowledges support from National Science Foundation grant AST-1909209 and NASA Exoplanet Research Program grant 20-XRP20_2-0119.

C.M.P. and M.F. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19 and 177/19). J.K. gratefully acknowledges the support of the Swedish National Space Agency (SNSA; DNR 2020-00104). P.K., M.S., J.S., and R.K. acknowledge the financial support of Inter-transfer grant No. LTT-20015. M.K. acknowledges support from ESAs PEA4000127913. M.E. acknowledges the support of the DFG priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (HA 3279/12-1). Funding for the Stellar Astrophysics Centre is provided by the Danish National Research Foundation (grant agreement No. DNRF106). D.G. and L.M.S. gratefully acknowledge financial support from the Cassa di Risparmio di Torino (CRT) foundation under grant No. 2018.2323 "Gaseous or rocky? Unveiling the nature of small worlds." This work is partly supported by JSPS KAKENHI grant No. JP20K14518 and SATELLITE Research from Astrobiology Center (AB022006).

Funding for the TESS mission is provided by NASA's Science Mission Directorate. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. This paper includes data collected by the TESS mission that are publicly available from the Mikulski Archive for Space Telescopes (MAST).

Observations in the paper made use of the NN-EXPLORE Exoplanet and Stellar Speckle Imager (NESSI). NESSI was funded by the NASA Exoplanet Exploration Program and the NASA Ames Research Center. NESSI was built at the Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. The authors are honored to be permitted to conduct observations on Iolkam Du'ag (Kitt Peak), a mountain within the Tohono O'odham Nation with particular significance to the Tohono O'odham people.

Observations in the paper made use of the High-Resolution Imaging instrument 'Alopeke. 'Alopeke was funded by the NASA Exoplanet Exploration Program and built at the NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. 'Alopeke was mounted on the Gemini North telescope of the international Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini partnership: the 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.

Based on observations made with the Nordic Optical Telescope, operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias under program 59-210.

Facilities: McDonald Observatory: 2.7 m Harlan J. Smith Telescope (Tull Coudé); Roque de los Muchachos Observatory: 3.58 m Telescopio Nazionale Galileo (HARPS-N), 2.56 m Nordic Optical Telescope (FIES), TESS; FLWO: 1.5 m (TRES); Ondřejov Observatory: 2 m Perek Telescope (OES); Gemini Observatory: 8.1 m Gemini North telescope ('Alopeke).

Software: pyaneti (Barragán et al. 2019), lightkurve (Lightkurve Collaboration et al. 2018), Austral (Endl et al. 2000), astropy (Astropy Collaboration et al. 2018), matplotlib (Hunter 2007), SME (Valenti & Piskunov 1996; Valenti & Fischer 2005; Piskunov & Valenti 2017), SpecMatch-Emp (Yee et al. 2017), forecaster (Chen & Kipping 2017), isochrones (Morton 2015a), TRICERATOPS (Giacalone et al. 2021).

Table A1 records the reconnaissance RV measurements for TOI-1670. Table A2 lists the relative Tull Coudé and HARPS-N precise RVs used in the joint global orbit fit of TOI-1670. See Section 2 for details.

Model selection balances the quality of the model fit and the model complexity, or number of parameters. Including extraneous parameters can lead to overfitting of the data, while excluding physically motivated aspects of the model can result in underfitting of the data and introduce bias. Our model fit is susceptible to the former scenario when we consider parameters that cannot be robustly estimated from the data, such as the mass of the sub-Neptune or the eccentricities of either planet.

Here we assess whether two alternative, less complex models are more justified by the data using different model selection criteria. We first compare model fits using only the RV data (excluding the light curve) for a one- and two-planet model, whereby we exclude the inner sub-Neptune in the one-planet fit. In the second comparison, we examine two joint model fits, one where the WJ eccentricity is fixed to zero and another in which the WJ eccentricity is a free parameter. In both cases, the inner sub-Neptune is not modeled with the RVs, and its eccentricity is fixed to zero.

Several metrics can be used to establish whether a model is justified by the data. For each fit, pyaneti reports the Bayesian information criterion (BIC; Schwarz 1978; Raftery 1986), defined as

$$\begin{eqnarray}&&\mathrm{BIC}\equiv -2\ \mathrm{ln}\ { \mathcal L }+k\ \mathrm{ln}\ N,\end{eqnarray} \tag{ B1 }$$

and the Akaike information criterion (AIC; Akaike 1998), defined as

$$\begin{eqnarray}&&\mathrm{AIC}\equiv -2\ \mathrm{ln}\ { \mathcal L }+2k.\end{eqnarray} \tag{ B2 }$$

Here ${ \mathcal L }$ is the model likelihood, k is the number of model parameters, and N is the number of data points used in the fit. We further calculate the AIC corrected for small sample sizes (AIC_c; Sugiura 1978; Burnham & Anderson 2004):

$$\begin{eqnarray}&&{\mathrm{AIC}}_{{\rm{c}}}=\mathrm{AIC}+\displaystyle \frac{2k(k+1)}{N-k-1}.\end{eqnarray} \tag{ B3 }$$

This metric is preferred over the AIC, as it can be understood as a relative model likelihood using Akaike weights (Akaike 1981; Burnham & Anderson 2004; Liddle 2007), where the weight w for each model i is

$$\begin{eqnarray}&&{w}_{i}=\displaystyle \frac{{e}^{-{\rm{\Delta }}{\mathrm{AIC}}_{c,i}/2}}{{\sum }_{r=1}^{R}{e}^{-{\rm{\Delta }}{\mathrm{AIC}}_{c,r}/2}}.\end{eqnarray} \tag{ B4 }$$

These metrics take into account the model likelihood while penalizing each additional free parameter. The BIC can be interpreted as a model evidence ratio, such that according to the Jeffreys scale, a BIC difference between two models of ≥5 is strong evidence and ≥10 is decisive evidence against the model with a higher BIC (Jeffreys 1935; Kass & Raftery 1995; Liddle 2007).

We find that both the BIC and Akaike weight criteria favor the one-planet model for the RV data and the WJ circular orbit model in the global joint fits. For the RV-only fits, we find BIC = −269.3, AIC = −297.3, AIC_c = −290.1, and an Akaike weight of >0.99 for the one-planet model and BIC = −254.0, AIC = −288.5, AIC_c = −277.0, and an Akaike weight of <0.01 for the two-planet model. The BIC comparison (ΔBIC > 15) and the Akaike weight (99% relative likelihood) strongly favor the one-planet model fit. In the global joint fits, we find BIC = −283,106.1, AIC = −283,260.4, AIC_c = −283,260.3, and an Akaike weight of 0.85 for the zero-eccentricity global model and BIC = −283,086.4, AIC = −283,256.9, AIC_c = −283,256.9, and an Akaike weight of 0.15 for the free eccentricity global model. Both the BIC comparison (ΔBIC > 15) and the Akaike weights (85% relative likelihood) suggest that there is strong evidence in favor of the circular orbit model.

Despite these model selection preferences, we find that the final parameter uncertainties are slightly underestimated when compared to the full model fit, which suggests that the less complex models are underfitting the data. Orbital parameter uncertainties are larger by approximately 10%–50% in the full global model (which includes the smaller inner planet and the eccentricities of both planets) compared to the simpler fixed zero-eccentricity joint model. For example, the uncertainties of the orbital period, planet mass, and planet radius of the WJ increased by 10%, 30%, and 40%, respectively. By neglecting the smaller inner planet and the eccentricity of either planet in order to avoid overfitting the data, the simpler model artificially underestimates the uncertainties of other model parameters. Ultimately, to more accurately determine the parameter uncertainties, we adopt the global simultaneous joint model fit and report upper limits on parameter posteriors when robust detections are not possible.

Figure C1 displays the posterior distributions of the fitted parameters from the global joint fit. See Section 3.4 for more details.

Zoom In Zoom Out Reset image size