Orbital Alignment of the Eccentric Warm Jupiter TOI-677 b

Warm giant planets, those with radii comparable to that of Jupiter and orbital periods in the range of ∼10–200 days, are well suited for advancing our understanding of close-in giant-planet formation. In contrast to their hotter counterparts—the so-called hot Jupiters (period ≲ 10 days)—warm giants are not expected to be subject to significant tidal friction (e.g., Alexander 1973; Zahn 1977; Hut 1981), thus better preserving their primordial orbital configurations. Consequently, characterization of warm giant orbits, albeit a significant observational challenge, can help better constrain planet formation models.

Mechanisms through which close-in giant planets form are hotly debated, but generally speaking, there are two families of models: (i) in situ formation and (ii) planetary migration. In situ scenarios rely on the core accretion model (Pollack et al. 1996) to work at small stello-centric distances, provided there is enough gas, and that critical cores can form from a sufficiently dense distribution of solids (e.g., Batygin et al. 2016) or from the consolidation of smaller cores (Boley et al. 2016). Planetary migration, on the other hand, relies on the significant reduction of planet's semimajor axis from initial separations beyond the ice line. Migration can be mediated by the tidal interaction with a gaseous, Keplerian disk (Lin & Papaloizou 1979; Goldreich & Tremaine 1980; Ward 1997), or mediated by extreme eccentricity growth followed by circularization and orbital decay (e.g., Mazeh & Shaham 1979), which result naturally from tidal friction (e.g., Goldreich 1963; Goldreich & Soter 1966; Hut 1981). This "high-eccentricity migration" can be triggered by planet–planet scattering (e.g., Rasio & Ford 1996) or by different types of secular perturbations (e.g., Eggleton & Kiseleva-Eggleton 2001; Wu & Murray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2011; Wu & Lithwick 2011; Petrovich 2015).

One may also choose to categorize these different formation mechanisms as either "dynamically cold" or "dynamically hot" (e.g., Tremaine 2015). In dynamically cold channels, the eccentricities, inclinations, and obliquities remain low; in dynamically hot evolution, on the other hand, the orbital elements can vary widely. For instance, in situ formation and disk-driven migration do not typically involve growth in inclination nor eccentricity, and can be deemed dynamically cold. High-eccentricity migration, on the other hand, is by definition, dynamically hot. Thus, measuring a warm giant's eccentricity and/or inclination relative to the stellar spin axis can serve as a discriminant between "hot" and "cold" dynamical histories, and consequently, serve as a crucial diagnostic of planet formation theories.

In principle, a sufficiently large number of spin–orbit measurements could prove extremely powerful for discerning between different planet migration models (e.g., Morton & Johnson 2011). Nonetheless, measurement of the spin–orbit angle (or projected stellar obliquity) λ is more difficult for warm Jupiters than for hot Jupiters, due to the rarity and longer duration of their transits. The angle between the stellar rotation axis and the planet's angular momentum vector, projected onto the plane of the sky, is measured through the observations of the Rossiter–McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924) with spectroscopic observations during the exoplanet transit, which has thus far limited these observations to close-in planets around bright stars. Recently, however, high resolution spectroscopic observations at large aperture telescopes have made RM measurements of warm Jupiter systems possible, suggesting that these planets might represent a population significantly different from their hotter counterparts (Rice et al. 2022).

In this work we present the sky-projected obliquity measurement for the warm Jupiter planet TOI-677 b (Jordán et al. 2020), through the analysis of the RM effect, observed with high resolution, spectroscopic observations of a single primary transit of the exoplanet. This study is structured as follows: in Section 1 an introduction to the analysis is presented; in Section 2 we briefly present the observations of the target with ESPRESSO and the subsequent data reduction process; in Section 3 the underlying analytical model, as well as the non-parametric noise model are presented, as well as the determination of the orbital obliquity angle from the modeling of the ESPRESSO transit data, while additional Fiber-fed Extended Range Optical Spectrograph (FEROS) radial velocity observations are analyzed together with previous data to infer a possible presence of an outer companion in the system; in Section 4 we discuss the possible implications of our results in the greater context of giant-planet formation theories; and finally in Section 5 we summarize this work and present the final conclusions of the study.

We observed a single primary transit of TOI-677 b on 2021 December 9, with the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observation (ESPRESSO) spectrograph (Pepe et al. 2021), installed at the Incoherent Combined Coudé Focus of ESO's Very Large Telescope (VLT) at Paranal Observatory, Chile. TOI-677 b is a 1.24 ± 0.07 M_J, 1.17 ± 0.03 R_J planet on an eccentric (e = 0.44 ± 0.02) 11.2366 ± 0.0001 day orbit around a late F-type star. This host star has an effective temperature of 6295 ± 77 K, with $\nu \sin i$ = 7.80 ± 0.19 km s⁻¹ (Jordán et al. 2020), estimated with the zaspe code (Brahm et al. 2017b). The stellar parameters determined in the detection study have been summarized in Table 1, in addition to some of those that have been determined from the spectral synthesis analysis of the out of transit ESPRESSO spectra obtained in this study, using zaspe.

Table 1. Stellar Parameters of TOI-677

Parameter	Jordán et al. (2020)	This Work
Age [Gyr]	${2.92}_{-0.73}^{+0.80}$	3.1 ± 0.7
J-band magnitude, mJ	8.722 ± 0.020	⋯
Mass, M⋆ [M⊙]	1.17 ± 0.06	${1.158}_{-0.027}^{+0.029}$
Radius, R⋆ [R⊙]	1.28 ± 0.03	1.281 ± 0.012
Temperature, Teff [K]	6295 ± 77	6295 ± 80
$\mathrm{log}g$ [dex]	4.291 ± 0.025	${4.286}_{-0.015}^{+0.016}$
Metallicity, [Fe/H] [dex]	0.00 ± 0.05	−0.02 ± 0.05
$\nu \sin {i}_{\star }$ [km s−1]	7.80 ± 0.19	7.42 ± 0.5

Download table as:  ASCIITypeset image

The observations were performed in the single-UT, HR mode (i.e., using the ${1}^{{\prime} }$ entrance fiber) with Unit Telescope 1 (UT1). The spectrograph records cross-dispersed echelle spectra through two blue- and red-optimized cameras, at a median resolving power ${ \mathcal R }$ ≈ 140,000. The detectors were read in the unbinned readout mode at an average spectral sampling of 4.5 pixels per resolution element, where each spectral order is recorded onto two slices owing to the anamorphic pupil slicing unit of the spectrograph. Starting at 04:49 UT, a total of 62 spectra were recorded (6 before, 40 during and 16 after transit) at exposure times of 180 s, with S/N values of ≈70 at 550 nm, across the two slices. A more detailed view of the observations is presented in the top panel of Figure 1. The observations were performed with the principal fiber (A) on the target and calibration fiber (B), which is at ${7}^{{\prime} }$ from A, on sky.

Zoom In Zoom Out Reset image size

The spectra were reduced using the dedicated data reduction pipeline (version 2.3.3), provided by the ESPRESSO consortium and ESO, and run on the esoreflex environment. Briefly, the reduction cascade includes bias and dark subtraction, flat-field correction, slice identification and wavelength calibration. For the purpose of solving the dispersion solution, day-time calibration frames taken with the Thorium–Argon lamp are used. We chose not to use the sky-subtracted spectra as lunar contamination in the science spectra is expected to be negligible due to its phase and angular distance (41% at 114 deg) and given the magnitude of the target (m_V = 9.82), thereby avoiding an additional source of noise in the final reduced spectra.

The pipeline also calculates the cross-correlation function (CCF) of the spectra with a binary mask for the stellar type matching closest the spectral type of the observed target (F9 in our case). We calculated the CCF at steps of 0.5 km s⁻¹, for ±40 km s⁻¹ centered on the expected systemic velocity of the star. The CCF from individual slices are summed (excluding those slices heavily contaminated by telluric absorption lines) and a Gaussian fit to this final CCF determines the central position of the profile and therefore the radial velocity. These calculated radial velocities together with their respective uncertainties, are presented in Table 4 (Appendix) and demonstrated in Figure 2, where the RM anomaly is clearly evident. Additionally, the pipeline provides S/N calculations for the individual spectral orders (middle panel of Figure 1), as well as a series of diagnostics determined from the CCF, which we used to search for correlations with the residuals of our eventual model. Further to the data reduction pipeline, we also used the dedicated Data Analysis Software (DAS, version 1.3.3) to determine activity indices from the spectra, such as the S-index and $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$, the latter of which is shown in the bottom panel of Figure 1.

Zoom In Zoom Out Reset image size

It has been shown that depending on the methodology through which the radial velocities are extracted from the observed spectra, one obtains different shapes for the RM effect (Boué et al. 2013). As we obtained RV values for TOI-677 through fitting of a Gaussian function to the CCF, we use the publicly available code ARoME (Boué et al. 2013) to model the RM effect, as this code provides instantaneous RM function definitions for RVs estimated through the cross-correlation and iodine cell techniques, as well as the weighted mean method. In the function definition, for the treatment of the stellar limb darkening, we use the quadratic law (Kopal 1950). We calculate the Cartesian coordinates of the planet at a given observation time as:

$$\begin{eqnarray}&&x(t)=r(t)(\cos \lambda \cos u(t)-\sin \lambda \sin u(t)\cos i)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&y(t)=r(t)(\sin \lambda \cos u(t)+\cos \lambda \sin u(t)\cos i)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&z(t)=r(t)(\sin u(t)\sin i)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\mathrm{with}\,r(t)=\displaystyle \frac{a(1-{e}^{2})}{1+e\cos \nu (t)}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{and}\,u(t)=\omega +\nu (t)\end{eqnarray} \tag{ 5 }$$

where λ is the sky-projected obliquity angle, i is the orbital inclination angle, a is the orbital semimajor axis scaled to the stellar radius, e is the orbital eccentricity, ν(t) is the true anomaly, ω is the argument of the periapsis, u is the argument of latitude (not to be confused with the limb darkening coefficients), and r is the radius from true anomaly. The x and y axes point along the plane of the sky (pointing arbitrarily to the right and up, respectively) and the z axis toward the observer.

Once the position of the planet is defined at each time of observation, the anomalous radial velocity value is calculated using the RM model, introduced above. The final model ${ \mathcal M }$ is subsequently the sum of the underlying RV trend of the star, which has a systemic and a planetary component, and this RM anomaly:

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal M }(t) & = & {\gamma }_{0}+{\mathrm{RV}}_{\mathrm{orb}}(t,P,K,e,{T}_{0},\omega )\\ & & \,\,\,+{\mathrm{RV}}_{\mathrm{RM}}(t,{T}_{0},a,P,e,i,\nu \sin {i}_{\star },\\ & & \ \ \omega ,\lambda ,{u}_{1},{u}_{2},{R}_{p},{\sigma }_{0},{\beta }_{0},\zeta \end{array}\end{eqnarray} \tag{ 6 }$$

where γ₀ is the systemic velocity, P is the orbital period, K is the RV semiamplitude, T₀ is the time of mid-transit, u₁ and u₂ are the limb-darkening coefficients, R_p is the planetary radius scaled to the stellar radius, σ₀ is the width of the CCF (FWHM of a Gaussian fit to the CCF), β₀ is the line width of the nonrotating star, and ζ is the stellar macroturbulence velocity.

3.1. Noise Consideration

We initially fitted the data with this analytical model, assuming only uncorrelated noise. However, the residuals of this fit, not presented in this manuscript, presented a distribution clearly deviating from the expected Gaussian. This points to the presence of correlated noise in the observations, caused by astrophysical and/or instrumental sources. The presence of active regions on the stellar surface has been shown to introduce anomalies in the observed photometric light curves both in and out of transit (Rackham et al. 2018), as well as in the radial velocity measurements (Huerta et al. 2008). We therefore measured a series of activity indices from our observed spectra, including the FWHM, bisector span and contrast of the CCF, $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ and the S-index, as well as line indices for Hα_0.6, Hα_1.6, ⁸ He i and Na i lines (Gomes da Silva et al. 2011). Two of these indices are presented in the bottom panel of Figure 1. Furthermore, the Ca i activity-insensitive line index was also calculated as control. The estimation of these line indices was made using the ACTIN python package (Gomes da Silva et al. 2018). Of all these indices, the variations in the Hα_1.6 index present the only clear sign of correlation with the residuals. This was searched for visually, as well as with a simple correlation analysis, this index showing a significant correlation (Pearson's correlation coefficient of 0.86). Namely, the two sharp decreases in this index at approximately 0.75 and 0.5 hr before mid-transit (indicated with red arrows in the bottom panel of Figure 1), coincide with RV deviations from the noise-free model. We checked for possible contamination of the Hα line with mirco-telluric absorption lines, whose variation could mimic line index variability. To this end we modeled the telluric absorption in the spectral series with ESO's molecfit (v. 4.2.3; Smette et al. 2015; Kausch et al. 2015), and note that those telluric lines, included in the region used for the calculation of the Hα_1.6 index, are not responsible for the variability observed. The calculated indices together with their uncertainties are presented in Table 4.

We incorporate these index measurements into our model via a Gaussian process (GP). The covariance matrix Σ of the GP is modeled with a squared exponential kernel:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{ij}}={ \mathcal A }\exp \left(\displaystyle \frac{-{({\rm{H}}{\alpha }_{i}-{\rm{H}}{\alpha }_{j})}^{2}}{2{{\ell }}^{2}}\right)+{\delta }_{{ij}}{\sigma }^{2},\end{eqnarray} \tag{ 7 }$$

with Hα being the line indices measured for the 1.6 Å bandpass, ${ \mathcal A }$ and ℓ the kernel amplitude and the length scale, respectively, and σ the uncorrelated or white noise in the data. The implementation of this GP noise model is performed with the GeePea python module (Gibson et al. 2012).

To sample the posterior distributions we ran 3 independent Markov chain Monte Carlo (MCMC) simulations of 120,000 steps each, using an Affine invariant ensemble sampler (Goodman & Weare 2010), assuming restrictive Gaussian prior distribution for the stellar macroturbulence velocity ⁹ (ζ ), as well as the orbital period (P), the eccentricity (e), the argument of periastron (ω), the RV semiamplitude (K), and the relative planetary radius R_p /R_⋆, whose values were taken from Jordán et al. (2020) through the analysis of TESS light curves and RV monitoring data. This approach was taken to ensure that the uncertainties on those parameters are correctly propagated. However, one caveat to note is that such restrictive Gaussian priors do not account for the impact of the existing correlations between the scaled semimajor axis, $\nu \sin {i}_{\star }$ and the eccentricity, and therefore the quoted uncertainties could be slightly underestimated. The priors are detailed in Table 2. The two coefficients of the quadratic limb-darkening law are fixed to those calculated from PHOENIX stellar spectrum model library of Husser et al. (2013), for the ESPRESSO bandpass using PyLDTK (Parviainen & Aigrain 2015). For all other model parameters, we assumed flat, uninformative prior distributions details of which are presented in Table 2. Additionally, for the kernel parameters we assumed very restrictive gamma priors with the shape parameter equal to 1 to maximize the distribution at 0 and the scale parameter to 0.1 in order to encourage the probability distributions to converge toward 0. This approach ensures that the included GP regressor contributes to the covariance only when there exists a significant correlation with the systematic noise present.

Table 2. Best-fit Parameter Values from the MCMC Simulations

Parameter	Prior a	Jordán et al. (2020)	This Work
Mid-transit time, T0 [−2459558 BJDTDB]	${ \mathcal U }(0.7191,0.8191)$	⋯	${0.769306}_{-0.000655}^{+0.000700}$
Orbital period, P [days]	${ \mathcal N }(11.2366,{0.00011}^{2})$	11.23660 ± 0.00011	11.23660 ± 0.00011
Orbital eccentricity, e	${ \mathcal N }(0.435,{0.024}^{2})$	0.435 ± 0.024	0.443 ± 0.021
Argument of periastron, ω [rad]	${ \mathcal N }(1.23,{0.06}^{2})$	1.23 ± 0.06	1.23 ± 0.06
RV semiamplitude, K [m s−1]	${ \mathcal N }(111.6,{0.5}^{2})$	⋯	111.6 ± 0.5
Systemic velocity, γ0 [−37 km s−1]	${ \mathcal U }(0.935,0.945)$	⋯	0.94068${}_{-0.00535}^{+0.00497}$
Scaled semimajor axis, a/R⋆	${ \mathcal U }(12,25)$	17.44 ± 0.69	15.86 ${}_{-1.32}^{+1.58}$
Relative planetary radius, Rp /R⋆	${ \mathcal N }(0.0942,{0.0012}^{2})$	${0.0942}_{-0.0012}^{+0.0010}$	0.0942 ± 0.0012
Orbital inclination, i [deg]	${ \mathcal U }(80,90)$	${85.86}_{-0.10}^{+0.11}$ b	84.80${}_{-0.79}^{+0.80}$
Orbital impact parameter, b (derived)	⋯	${0.723}_{-0.024}^{+0.018}$	${0.858}_{-0.220}^{+0.272}$
Sky-projected obliquity, λ [deg]	${ \mathcal U }(-180,180)$	⋯	0.3 ± 1.3
Equatorial stellar rotation, $\nu \sin {i}_{\star }$ [km s−1]	${ \mathcal U }(0,15)$	7.80 ± 0.19	6.91${}_{-1.20}^{+1.32}$
Stellar macroturbulence velocity, ζ [km s−1]	${ \mathcal N }(5.5,{0.5}^{2})$	⋯	5.53 ${}_{-0.51}^{+0.50}$
Linear limb-darkening coefficient, u1	⋯	0.50TESS	0.5153 (fixed)
Quadratic limb-darkening coefficient, u2	⋯	−0.06TESS	0.1518 (fixed)
GP kernel amplitude, ${ \mathcal A }$ [m s−1]	Γ(1, 0.1)	⋯	${6.4}_{-1.7}^{+2.3}$
GP kernel regressor length scale, ℓ	Γ(1, 0.1)	⋯	${0.0067}_{-0.0028}^{+0.0022}$
White noise, σw [m s−1]	${ \mathcal U }(0,\infty )$	⋯	4.2${}_{-0.4}^{+0.5}$

Notes.

^a The distributions in the prior column are defined as: ${ \mathcal U }(l,u)$ is a uniform distribution between l and u, ${ \mathcal N }(\mu ,{\sigma }^{2})$ is a normal distribution with a mean of μ and a variance of σ², and Γ(k, θ) is a gamma distribution with a shape parameter k and a scale parameter θ. ^b Inclination was reported incorrectly in Jordán et al. (2020) as i = 87.63 deg. The fitted parameter in that work was actually b, whose reported value is correct. The reported inclination was derived incorrectly using the relation between i and b valid only for a circular orbit.

Download table as:  ASCIITypeset image

The best-fit analytical model, as well as the noise model, together with their residuals, are presented in Figure 2. The posterior probability distributions and the joint posteriors are presented in Figure 3, with the independent chains over-plotted (the initial 20,000 steps of which are burnt in). The best-fit results for the fitted parameters, given in Table 2, are derived from the median and the 16th and 84th percentiles of those distributions. We obtain a perfectly aligned sky-projected orbit of TOI-677 b with respect to the spin orbit of its host star, with λ = 0.3 ± 1.3 deg. All other estimated parameters are in general agreement with previously obtained results, whereby our parameters result in a slightly more inclined orbit.

Zoom In Zoom Out Reset image size

3.1.1. True Obliquity

As the angle measured from this analysis of the RM effect is the sky-projected (λ) portion of the true obliquity angle (ψ), we attempted to estimate this true value. One can potentially de-project this measurement if the stellar line-of-sight inclination (i_⋆) can be measured via ${\sin }^{-1}\left(\nu \sin {i}_{\star }/(2\pi {R}_{\star }/{P}_{\mathrm{rot}})\right)$, although this approach suffers from biases due to the fact that $\nu \sin {i}_{\star }$ and 2π R_⋆/P_rot (≡ν) are not statistically independent measurements (Morton & Winn 2014; Masuda & Winn 2020). ψ would subsequently be estimated through its geometrical relation to the planetary orbital plane inclination (i) and i_⋆:

$$\begin{eqnarray}&&\cos \psi =\cos {i}_{\star }\cos i+\sin {i}_{\star }\sin i\cos \lambda .\end{eqnarray} \tag{ 8 }$$

To this effect, we attempted to measure the stellar rotation period (P_rot) through modulations in the TESS light curves, both the simple aperture photometry (SAP) and the presearch data conditioning (PDC) LCs, imprinted by the rotation of active regions on the star. However, a Lomb–Scargle periodogram search of all available observations of TOI-677 from TESS sectors 9, 10, 35 and 36 did not result in a viable detection.

3.2. Possible Outer Companion

In the analysis of TOI-677 RV measurements, Jordán et al. (2020) detected an underlying slope of 1.58 ± 0.19 m s⁻¹day⁻¹. In order to investigate possible roots of such trend in the data, in addition to the ESPRESSO data presented previously, we also observed TOI-677 with FEROS, mounted at the MPG/ESO 2.2 m telescope at La Silla observatory, in nine distinct epochs. The stellar RVs were subsequently derived from the spectra via processing with the CERES pipeline (Brahm et al. 2017a), similar to what was performed in Jordán et al. (2020). These additional radial velocities are given in Table 5 in Appendix, which together with the initial RV data of Jordán et al. (2020) and the ESPRESSO data presented in this work, ¹⁰ are used to search for possible outer companions to TOI-677 b.

We analyze this newly assembled RV data using the juliet package (Espinoza et al. 2019), which utilizes Keplerian orbital radial velocity perturbation formalism via the radvel package (Fulton et al. 2018). In contrast to the model fit performed by Jordán et al. (2020), instead of an underlying linear trend, we include a second body inducing the long period trend observed in the data (Figure 4). We performed two separate fits to the data, whereby the outer component is assumed to be on either a circular or eccentric orbit. In both scenarios, we fitted for orbital parameters of both bodies, as well as instrumental parameters. We found the instrumental dependent systemic velocity (γ), as well as instrumental jitter (η) values consistent with those reported by Jordán et al. (2020), in both sets of analyses. We sampled the Bayesian posterior distributions using the importance nested sampling and MultiNest algorithms (Feroz et al. 2019), implemented by juliet via the PyMultiNest python package (Buchner et al. 2014), using 2000 live points. The two sets of posterior codistributions and probability distribution functions are presented in Figure 5, where orbital parameters for only the possible outer companion are presented. It must be noted that we do not observe any other significant peaks in the posterior distributions of any of the parameters in either fit. All derived orbital parameters for both bodies in the system, from both modeling approaches, have been presented in Table 3, with the best-fit models and their respective residuals shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Best-fit Parameter Values from the Nested Sampling Simulations for the Two-planet Fit of the RV data, Assuming both Circular and Eccentric Orbits for the Possible Outer Companion

Parameter	Prior	Circular Orbit	Eccentric Orbit
Mid-transit, T0,1 −2458547 [BJD]	${ \mathcal N }(0.4743,0.0012)$	${0.474303}_{-0.001133}^{+0.001118}$	${0.474312}_{-0.001045}^{+0.001060}$
Period, P1 [days]	${ \mathcal U }(11.1,11.5)$	${11.236518}_{-0.000609}^{+0.000589}$	11.236168 ± 0.000537
Eccentricity, e1	${ \mathcal N }(0.434,0.05)$	0.442 ± 0.018	0.436 ± 0.019
Argument of periastron, ω1 [deg]	${ \mathcal N }(70.47,1.0)$	70.47 ± 0.09	70.47 ± 0.44
RV semiamplitude (planet 1), K⋆,1 [m s−1]	${ \mathcal N }(111.6,5.0)$	111.59 ± 2.84	${113.48}_{-2.64}^{+2.74}$
Mass,mp,1 [Mjup]	⋯	1.24 ± 0.13	1.26 ± 0.14
Mid-transit, T0,2 −2460000 [BJD]	${ \mathcal U }(0,{10}^{4})$	${579.329132}_{-61.530429}^{+53.885799}$	${1755.376054}_{-549.826780}^{+392.459980}$
Period, P2 [days]	${ \mathcal U }(11.5,{10}^{4})$	$\gtrsim {3953.751212}_{-127.364290}^{+111.367231}$	$\gtrsim {4901.138760}_{-644.759743}^{+523.775478}$
Eccentricity, e2	${ \mathcal U }(0,1)$	0 (fixed)	${0.436}_{-0.058}^{+0.067}$
Argument of periastron, ω2 [deg]	${ \mathcal U }(-180,180)$	90 (fixed)	$-{123.68}_{-6.99}^{+6.96}$
RV semiamplitude (planet 2), K⋆,2 [m s−1]	${ \mathcal U }(0,1000)$	${450.98}_{-13.46}^{+12.39}$	${594.91}_{-35.34}^{+35.01}$
Mass lower limit,${m}_{p,2}\sin i$ [Mjup]	⋯	≳39.20 ± 2.81	≳49.99 ± 14.07
Log evidence, $\mathrm{ln}{ \mathcal Z }$	⋯	−683.72 ± 0.18	−675.04 ± 0.04

Notes.

^a Subscript 1 refers to the inner planet and 2 to the outer companion. The period, and consequently the lower mass limit, of the possible outer companion are presented only as a lower limits since the fitted orbit is not closed.

Download table as:  ASCIITypeset image

In order to evaluate the statistical significance of the two-body model as compared to the single-planet case, we also fitted all the RV data assuming just one planet in the system, which results in ${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }\sim 1500$, as compared to the two-body scenarios, i.e., pointing to a significant preference of the data for the two-body model and the possible presence of an outer companion. Furthermore, there is also strong preference for an eccentric orbit of this possible outer companion, as compared to a circular orbit, from the ratio of the likelihoods of the two models, with ${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }=8.68$ in its favor. This points to a very strong ($2{\rm{\Delta }}\mathrm{ln}{ \mathcal Z }\gt 10$) preference for the two-body model with an outer companion on an eccentric (e = 0.44 ± 0.07) and very wide orbit of ≳4901 days (≳13.4 yr) period. These estimated period values are taken only as lower limits since the orbit is not closed. The lower limit for the mass of this possible outer companion is estimated as ∼39 and ∼50 M_jup for the circular and eccentric cases, respectively, putting it in the brown dwarf regime in either case. Assuming the outer companion is on a relatively coplanar orbit to the inner companion, the true mass of this outer companion is likely close to this lower limit.

However, it must be stressed that such analysis only points to the possible presence of the outer companion, since the data simply do not cover a long enough baseline for any definitive conclusions to be made. This fact is reflected in the relatively large uncertainties in the determination of the orbital parameters for this potential outer companion, as presented in Table 3.

Additionally, we attempted to fit the RV data with a three-planet model, keeping e and ω as free parameters; however, convergence was not achieved as the stopping criterion for the nested sampling algorithm could not be reached. The final $\mathrm{ln}{ \mathcal Z }$ at the moment of stopping the algorithm still pointed to a very strong preference for the two-planet model.

TOI-677 b is now one of ∼ 200 exoplanets for which the projected stellar spin–orbit misalignment λ has been measured. ¹¹ Nearly 85% of these systems correspond to close-in gas giants (R_p > 0.7R_J, P ≤ 200 days), of which ∼90% are "hot" (P < 10 days) and ∼10% are "warm." We show the distribution of ∣λ∣ in Figure 6, plotted as a function of planet semimajor axis, where the symbol sizes scale with planet radii, and the color scale represents orbital eccentricity. From the figure, one can distinguish the hot population from the warm population: in the former case, obliquities are distributed broadly (e.g., Fabrycky & Winn 2009; Morton & Winn 2014; Muñoz & Perets 2018), whereas in the latter case, obliquities are distributed rather narrowly (Rice et al. 2022). In addition, as it is well known, the most compact orbits have zero eccentricity, an indication of circularization owing to tidal dissipation in the planet (e.g., Goldreich & Soter 1966). Since the circularization rate is a steep function of separation (Goldreich & Soter 1966; Hut 1981), wider orbits may allow for nonzero eccentricity. Indeed, several warm Jupiters have eccentricities above 0.4: e.g., HD 80606b (e ≈ 0.93; Naef et al. 2001), Corot-10 b (e ≈ 0.53; Bonomo et al. 2010), Kepler-419 b (e ≈ 0.85; Dawson et al. 2014), Kepler-420 b (e ≈ 0.77; Santerne et al. 2014), or TOI-2179b (e ≈ 0.58; Schlecker et al. 2020).

Zoom In Zoom Out Reset image size

We do not expect TOI-677 b to have been fully circularized over the lifetime of its host star (≈3 Gyr; Jordán et al. 2020). Indeed, assuming that tidal dissipation takes place primarily within the planet, the characteristic circularization timescale is given by:

$$\begin{eqnarray}&&{\tau }_{\mathrm{circ}}\equiv -\displaystyle \frac{e}{\dot{e}}=\displaystyle \frac{2F(e)}{7}{\tau }_{\mathrm{dec}}\,\,\,\mathrm{with}\,\,\,{\tau }_{\mathrm{dec}}\equiv \displaystyle \frac{P}{9\pi }{Q}_{{\rm{p}}}^{\prime} \displaystyle \frac{{m}_{{\rm{p}}}}{{M}_{\star }}{\left(\displaystyle \frac{a}{{R}_{{\rm{p}}}}\right)}^{5}\end{eqnarray} \tag{ 9 }$$

(Goldreich & Soter 1966), where τ_dec is the characteristic orbital decay timescale, ${Q}_{{\rm{p}}}^{\prime}$ is the planet's modified tidal quality factor (e.g., Goldreich & Soter 1966; Ogilvie & Lin 2007) and F(e) is an eccentricity-dependent correction factor (Hut 1981). For e = 0.435, and assuming that the planet is in pseudosynchronous rotation, ¹² we have F(e) ≈ 0.2. Further assuming that ${Q}_{{\rm{p}}}^{\prime} ={10}^{5}\mbox{--}{10}^{6}$ (e.g., Goldreich & Soter 1966; Yoder & Peale 1981), we have τ_circ ∼ 1–10 Gyr for TOI-677 b. Had we not ignored tidal dissipation in the star, ¹³ these timescales would be negligibly shorter for any value of ${Q}_{\star }^{\prime}$ greater than 10⁷, which is to be expected in this type of system (e.g., Barker & Ogilvie 2010, 2011; Penev & Sasselov 2011). Similarly, for such values of ${Q}_{\star }^{\prime}$, tidal realignment of the star itself would take hundreds of times longer than the age of the system. Moreover, even if stellar realignment did take place, it would come at the expense of planetary engulfment (Barker & Ogilvie 2009). Thus, only a tidal theory that goes well beyond weak friction (e.g., Lai 2012) could possibly permit the realignment of the stellar spin while sparing the planet's orbit.

TOI-677 b belongs to an intriguing, emerging group of eccentric, spin–orbit-aligned systems bracketed between the hot and warm populations. At these orbital separations, dissipation of energy within the star—responsible for obliquity damping (Hut 1981)—is extremely weak, which rebuffs the hypothesis of tidal reprocessing of the spin–orbit angle over long timescales (e.g., see Winn et al. 2010; Albrecht et al. 2012). Instead, planets such as TOI-677 b are likely to have attained their unusual orbital configurations soon after the planet formation and migration process was finalized.

Given these properties, systems like TOI-677 b present a significant challenge to standard theories of planet migration, both for the dynamically hot and dynamically cold scenarios. On the one hand, planets like TOI-677 b are unlikely to have attained their eccentricities during dynamically cold, disk-driven planet migration. On the other hand, low spin–orbit alignment would also disfavor dynamically hot, eccentricity excitation mechanisms, which are usually accompanied by large changes in inclination. Moreover, the unlikeliness of large-amplitude eccentricity oscillations being responsible for these elongated orbits would in turn reject the notion that TOI-677 b-like systems are "failed" or "proto" hot Jupiters (e.g., Dong et al. 2014; Petrovich & Tremaine 2016). In fact, while there are indeed two planetary systems–HD 80606 and Kepler-420 –that appear to be quintessential examples of ongoing (or failed) high-eccentricity migration driven by Lidov–Kozai oscillations (Wu & Murray 2003), these planets do not appear to be representative of the warm giant population. They have, in addition to high eccentricity, known binary companions and high obliquities (see Figure 6).

Thus, we must consider the option of in situ excitation of high eccentricity while at low inclination. This is indeed possible due to an exterior, nearly coplanar perturber of mass ≳m_p and of moderate-to-high eccentricity (Lee & Peale 2003). If sufficiently eccentric, the candidate substellar mass companion discovered in this study presents an ideal potential explanation for TOI-677 b's peculiar orbit. In order to excite the planet's eccentricity from a circular orbit to its current value of e₁ = 0.435, the outer perturber must satisfy the approximate condition:

$$\begin{eqnarray}&&\frac{{e}_{2}}{(1-{e}_{2}^{2})}\geqslant \frac{4}{15}\frac{{e}_{1}}{1+\frac{3}{4}{e}_{1}^{2}}{\left(\frac{{P}_{2}}{{P}_{1}}\right)}^{2/3}\approx 5,\end{eqnarray} \tag{ 10 }$$

proposed by Petrovich (2015). Being a necessary yet not a sufficient condition (it ignores suppressing effects such as general relativistic precession; e.g., Liu et al. 2015), Equation (10) can only provide a lower limit on the required value of e₂. The currently estimated eccentricity and period for the outer companion, given in Table 3, fail to satisfy this condition by one order of magnitude.

Consequently, we may conclude that no known mechanism of planet migration can explain the current orbital eccentricity and alignment of TOI-677 b with the currently known objects in the system. This puzzle highlights the importance of obtaining RM observations of a wider class of exoplanets, pushing the boundaries of high-precision spectroscopy.

In this study we presented single transit observations of the warm Jupiter TOI-677 b with the ESPRESSO spectrograph, obtaining the Rossiter–McLaughlin effect in order to measure the sky-projected obliquity angle λ. This angle was determined to be 0.3 ± 1.3 deg, putting the planet on a perfectly aligned orbit with the stellar spin axis. In modeling the effect together with the correlated noise we uncovered a strong correlation with the Hα activity index, which was used as the regressor in the calculation of the covariance matrix in the Gaussian Process model. In the analysis, MCMC methods were used to determine parameter uncertainties, while evaluating posterior codistributions. An attempt was made to measure the true obliquity angle of the system, which was unsuccessful due to the inability to measure the stellar rotation period from TESS photometry, owing to the absence of activity-induced light-curve modulations. Follow-up radial velocity monitoring revealed a long-term periodic signal, which together with the initial data from Jordán et al. (2020) was modeled with a two-component Keplerian model. The analysis revealed a significant preference for a companion on an eccentric orbit, as opposed to a circular one. This solution pointed to the possible presence of a companion with a lower mass limit in the brown dwarf regime (M_p ≈ 50 M_J), on a wide (P ≈ 13.4 yr) and moderately eccentric (e ≈ 0.44) orbit. Posteriors obtained from a nested sampling approach revealed relatively well-constrained distributions, although no definitive conclusion was made about the presence of this outer companion, due to the lack of sufficient coverage of this long orbital period. We finally discussed the orbital architecture of this system in the context of currently known planet migration mechanisms, and the challenges it poses to them. Namely, while it is likely the system attained its eccentricity through disk migration, the aligned orbit disfavors eccentricity excitation mechanisms. Furthermore, we argued that it is also highly unlikely that the system is a failed or proto hot Jupiter. We finally discussed the possibility of an in situ excitation of the eccentricity by the substellar outer companion. However, with the current and limited analysis, it was concluded that this outer companion does not possess high enough eccentricity to cause the elevated eccentricity in the inner planetary companion. This result, subsequently, highlights the need and the importance of obtaining RM measurements for planets in the warm giant regime, to better test and refine planet migration theories.

E.S., A.J., R.B., and C.P. acknowledge support from ANID—Millennium Science Initiative—ICN12_009. A.J. acknowledges additional support from FONDECYT project 1210718. R.B. acknowledges support from FONDECYT project 11200751. C.P. acknowledges support from ANID Millennium Science Initiative-ICN12_009, CATA-Basal AFB-170002, ANID BASAL project FB210003, FONDECYT Regular grant 1210425, CASSACA grant CCJRF2105, and ANID+REC Convocatoria Nacional subvencion a la instalacion en la Academia convocatoria 2020 PAI77200076.

Facilities: VLT/UT1(ESPRESSO) - , MPG/ESO: 2.2m(FEROS) - .

Software: ACTIN (Gomes da Silva et al. 2018), ARoME (Boué et al. 2013), astropy (Astropy Collaboration et al. 2022), CERES (Brahm et al. 2017b), corner (Foreman-Mackey 2016), emcee (Foreman-Mackey et al. 2013), GeePea (Gibson et al. 2012), juliet (Espinoza et al. 2019), matplotlib (Hunter 2007), molecfit (Smette et al. 2015), MultiNest (Feroz et al. 2019), NumPy (Harris et al. 2020), PyLDTK (Parviainen & Aigrain 2015), PyMultiNest (Buchner et al. 2014), radvel (Fulton et al. 2018), SciPy (Virtanen et al. 2020).

In this appendix we present the radial velocities, as well as the Hα_1.6 index, measured for the TOI-677 spectra obtained with ESPRESSO in Table 4, as well as the additional RVs obtained with FEROS in Table 5.