A Highly Eccentric Warm Jupiter Orbiting TIC 237913194

For identifying warm Jupiter candidates, we generated light curves for all bright stars of the TICv8 catalog from the full -frame images (FFIs; Jenkins et al. 2016) of TESS using the tesseract ¹⁶ package (F. I. Rojas 2020, in preparation). Briefly, tesseract receives any TIC ID or coordinate as input and performs simple aperture photometry on the FFIs via the TESSCut (Brasseur et al. 2019) and lightkurve (Lightkurve Collaboration et al. 2018) packages. Aperture selection was done following Lund et al. (2015). Specifically, 293,253 and 479,184 light curves of bright objects (T < 14 mag) have been generated from Sectors 1 and 2, respectively. For identifying warm Jupiter candidates, we apply a simple algorithm that goes through the light curve searching for zones that significantly deviate in the negative direction from the median flux around a given region. Then we check by visual inspection if these zones are consistent with a transit-like feature. This procedure allows us in principle to identify also single transiters in a given TESS Sector (e.g., Gill et al. 2020). By using this algorithm we found that the star TIC 237913194 presented transit-like periodic features in the two first sectors of TESS. An initial fit to the photometric data indicated a period of P ≈ 15.17 and a transit depth of δ_TESS ≈ 0.8%, consistent with the properties of a warm giant candidate given the parameters of the star according to the TICv8 catalog. The TESS light curve of TIC 237913194 is shown in the upper panel of Figure 1.

Zoom In Zoom Out Reset image size

Due to the limited angular resolution of TESS, ground-based photometry is required to reject false-positive scenarios like blended eclipsing binaries. TIC 237913194 was observed on the night of 2019 December 12 with the 0.7 m CHAT ¹⁷ installed at the Las Campanas Observatory. Observations were performed with a Sloan i^' photometric filter using a mild defocus and exposure times of 110 s. We processed the data with a dedicated pipeline developed to produce high-precision light curves using differential photometry (e.g., Espinoza et al. 2019a; Jones et al. 2019) with the LCOGT 1.0 m telescopes (Brown et al. 2013). The optimal photometric precision was obtained with an aperture of 14 pixels (83). We plot the obtained light curve in the bottom right panel of Figure 1. We recovered a full transit, which confirms that the transit-like features identified in the TESS data occur in TIC 237913194. The transit depth of δ_CHAT = 0.0087 ± 0.0004 is consistent with the signal identified in the TESS photometry.

Because of a rather grazing transit configuration, the posterior probability densities from our initial fits contained a strong degeneracy between the scaled planetary radius R_P/R_⋆ and the impact parameter b. To lift this degeneracy and to improve the constraint on the planet radius, we obtained additional transit photometry with the LCOGT Network on 2020 July 16 (see Figure 1). The measurements were taken in the i^' band and cover all phases of the transit. To maximize photometric precision, we chose an aperture of 24 pixels (94). We recover a transit depth of δ_LCOGT = 0.0083 ± 0.0002. Within the uncertainties, this is consistent with the values obtained for the other instruments. Including the additional data in the fit leaves only little residual correlation between R_P/R_⋆and b and strongly improved the posterior on R_P (see Section 3.2.4).

We make all our follow-up light curves available on exoFOP. ¹⁸

We obtained high-resolution (R = 48,000) spectra with FEROS (Kaufer et al. 1999), mounted at the 2.2m MPG telescope at La Silla Observatory. In total, 25 exposures of 1200 s were taken between 2019 June 19 and 2020 March 9. From these, we extracted RVs using the CERES pipeline (Brahm et al. 2017a), which performs all steps from bias, dark, and flat-field calibration to cross-correlation matching of the resulting spectrum with a G2-type binary mask. The observations were performed in the simultaneous calibration mode for tracking the instrumental velocity drift produced by changes in the spectrograph environment. This procedure involves the monitoring of a ThAr spectrum with a second fiber. The typical signal-to-noise ratio of these spectra was about 70. The time series of the FEROS RV measurements are shown in Figure 2 and listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1. FEROS RVs and Accompanying Data for TIC 237913194 Used in This Paper

BJD	RV (km s−1)	${\sigma }_{\mathrm{RV}}$ (km s−1)	texp (s)	BIS (km s−1)	${\sigma }_{\mathrm{BIS}}$ (km s−1)
2458669.798	29.536	0.009	1200	−0.016	0.012
2458670.816	29.519	0.010	1200	−0.014	0.013
2458672.879	29.533	0.012	1200	0.014	0.016
2458718.916	29.505	0.007	1200	−0.006	0.011
2458721.768	29.528	0.007	1200	−0.000	0.011
2458722.674	29.580	0.008	1200	−0.001	0.011
2458723.758	29.595	0.007	1200	−0.019	0.011
2458724.793	29.611	0.009	1200	−0.019	0.012
2458742.740	29.892	0.009	1200	−0.002	0.013
2458783.690	29.545	0.007	1200	−0.033	0.011
2458785.610	29.646	0.010	1200	−0.026	0.013
2458787.671	29.802	0.008	1200	0.012	0.012
2458791.610	29.521	0.009	1200	−0.000	0.012
2458792.571	29.508	0.008	1200	0.007	0.011
2458800.668	29.614	0.008	1200	−0.017	0.011
2458801.688	29.673	0.008	1200	0.014	0.011
2458802.619	29.739	0.007	1200	−0.014	0.010
2458805.675	29.543	0.007	1200	−0.001	0.010
2458847.592	29.701	0.008	1200	−0.024	0.011
2458848.591	29.858	0.008	1200	0.006	0.011
2458849.592	29.805	0.008	1200	0.005	0.011
2458850.541	29.619	0.007	1200	−0.010	0.011
2458911.524	29.560	0.008	1200	−0.038	0.011
2458915.515	29.480	0.010	1200	0.023	0.013
2458917.514	29.510	0.010	1200	0.039	0.013

Download table as:  ASCIITypeset image

We checked for possible close sources that could contaminate our photometric aperture with their light. Any sources within  ∼10'', which is the photometric aperture we used for our LCOGT photometry, could cause such contamination. After querying the GAIA DR2 catalog (Gaia Collaboration et al. 2018) we found the closest source to TIC 237913194 at an angular separation of ∼46''. We thus find no evidence for significant contamination of our photometry.

For characterizing the host star, we first determined its atmospheric parameters from the co-added FEROS spectra. Specifically, we used the ZASPE code (Brahm et al. 2017b) which compares the observed spectrum against a grid of synthetic ones generated from the ATLAS9 model atmospheres (Castelli & Kurucz 2004). We then used the PARSEC evolutionary models (Bressan et al. 2012), as described in Brahm et al. (2019), to determine the physical parameters of the star. Briefly, we compared the observed broadband photometric magnitudes of the star with those generated with models having different physical parameters by taking into account the distance determined from the Gaia DR2 parallax and assuming an extinction law of Cardelli et al. (1989) dependent on the A_V parameter. The parameter space was explored using the emcee (Foreman-Mackey et al. 2013) package. The obtained atmospheric and physical parameters of TIC 237913194 are listed in Table 2. TIC 237913194 is a main-sequence G-type star with a mass of ${M}_{\star }$ = ${1.026}_{-0.055}^{+0.057}$ ${M}_{\odot }$, a radius of ${R}_{\star }$ = ${1.088}_{-0.012}^{+0.012}$ ${R}_{\odot }$, and an age of $5.7\pm 1.7$ Gyr. TIC 237913194 is slightly metal-rich ([Fe/H] = $+0.14\pm 0.05$ dex) and has an effective temperature of ${T}_{\mathrm{eff}}$ = $5788\pm 80$ K. We note that the quoted uncertainties do not account for possible systematic errors in the stellar evolutionary models.

Table 2. Stellar Properties of TIC 237913194

Parameter	Value	Reference
Names	TIC 237913194	TESS
	2MASS J01294694-6044238	2MASS
	UCAC4 147-001388	UCAC 4
R.A. (J2015.5)	01h29m4699	Gaia
Decl. (J2015.5)	−60d44m2367	Gaia
pmR.A. (mas yr−1)	18.053 ± 0.036	Gaia
pmdecl. (mas yr−1)	10.523 ± 0.034	Gaia
π (mas)	3.23 ± 0.02	Gaia
T (mag)	11.486 ± 0.006	TESS
B (mag)	12.746 ± 0.015	APASS
V (mag)	12.144 ± 0.069	APASS
J (mag)	10.858 ± 0.023	2MASS
H (mag)	10.571 ± 0.024	2MASS
Ks (mag)	10.485 ± 0.021	2MASS
WISE1 (mag)	10.463 ± 0.023	WISE
WISE2 (mag)	10.518 ± 0.021	WISE
WISE3 (mag)	10.408 ± 0.059	WISE
${T}_{\mathrm{eff}}$(K)	$5788\pm 80$	This work
$\mathrm{log}g$(dex)	$4.376\pm 0.021$	This work
$[\mathrm{Fe}/{\rm{H}}]$(dex)	$+0.14\pm 0.05$	This work
$v\sin i$(km s−1)	$2.18\pm 0.41$	This work
${M}_{\star }$(${M}_{\odot }$)	${1.026}_{-0.055}^{+0.057}$	This work
${R}_{\star }$(${R}_{\odot }$)	${1.088}_{-0.012}^{+0.012}$	This work
${L}_{\star }$ (L⊙)	1.196 ± 0.050	This work
Age (Gyr)	$5.7\pm 1.7$	This work
AV (mag)	${0.117}_{-0.063}^{+0.068}$	This work
${\rho }_{\star }$ (gcm−3)	$1.12\pm 0.11$	This work

Download table as:  ASCIITypeset image

For the joint photometry and RV modeling, we employed the python package juliet (Espinoza et al. 2019b). This tool uses existing codes to model transit photometry (batman; Kreidberg 2015) and RV time series (radvel; Fulton et al. 2018). It further allows us to incorporate Gaussian process regression (GP) via the celerite package (Foreman-Mackey et al. 2017), which we employ to model systematic nuisance signals. To explore the parameter space, it uses the MultiNest nested sampling technique (Feroz et al. 2009), implemented in the pyMultiNest software package (Buchner et al. 2014). juliet further calculates evidences ${Z}_{i}=P({M}_{i}(\theta )| { \mathcal D })$ for models M_i with sets of parameters θ given the data ${ \mathcal D }$. To compare two models M_i , M_j , we compute the differences of their log-evidences,

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{ln}{Z}_{i,j}=\mathrm{ln}{Z}_{i}/{Z}_{j}=\mathrm{ln}\left[P({M}_{i}(\theta )| { \mathcal D })/P({M}_{j}(\theta )| { \mathcal D })\right].\end{eqnarray} \tag{ 1 }$$

Here, we adopted a general rule of thumb that if ${\rm{\Delta }}\mathrm{ln}{Z}_{i,j}\geqslant 3$, the model with the larger log-evidence is favored. If ${\rm{\Delta }}\mathrm{ln}{Z}_{i,j}\lesssim 3$, we consider the models to be indistinguishable and prefer the simpler one. As the MultiNest algorithm is known for showing scatter in $\mathrm{ln}Z$ that exceeds the reported uncertainties (e.g., Nelson et al. 2020), we always repeated the calculations several times. The variations among such runs were always smaller than one and therefore negligible for our purposes.

3.2.1. Model Parameters

3.2.2. Limb Darkening

3.2.3. RV Analysis

The RV time series show a strong signal with a period corresponding to the candidate transiting planet (P = 15.17 days; see Figures 2 and 3). To assess the evidence of this signal being of planetary origin, we compared models with and without a planet based on only the FEROS RV data set. We further evaluated models including more than one planet and compared the log-evidences of the different cases:

• 1.  
  No planet. We assumed that all RV variations are due to astronomical and instrumental jitter. The only free parameters were μ_FEROS and σ_FEROS using the same prior distributions as in Table 3. This flat model resulted in a log-evidence $\mathrm{ln}(Z)=-161$.
• 2.  
  Single planet. We assumed there is a planetary signal in the RV data and widened the orbital period prior to a uniform distribution of 1–30 days. The RV semi-amplitude K was free to vary between zero and 1000 m s⁻¹ (uniform prior). For T₀, we chose a uniform prior ranging from the first photometric observation to 124 days later, which corresponds to half the RV baseline. All other free parameters had the same priors as in our nominal model. This fit converged to a similar solution as our final model with a period distribution consistent with the intervals between the observed transit events, $\mathrm{ln}(Z)=-135$.
• 3.  
  Single planet, circular orbit. Same as number 2, but fixing the eccentricity to zero (and ω to an arbitrary 90°). The fit converged to a solution with P similar to the period distribution in our final joint fit, but the jitter term is strongly increased to account for the large RV variations. With $\mathrm{ln}(Z)=-160$, the evidence of this model is similar to the one belonging to the no-planet hypothesis.
• 4.  
  Two planets, circular orbits. Same as number 3, but assuming a second planet in the system. For this hypothetical additional planet, we let the orbital period vary within 1–30 days and used the same uniform prior ${ \mathcal U }(0,1000){\rm{m}}\,{{\rm{s}}}^{-1}$ for its RV semi-amplitude K₂. The ∼15 day periodicity is recovered, but no stable solution in favor of a two-planet scenario is evident, $\mathrm{ln}(Z)=-164$.
• 5.  
  Two planets, circular and eccentric orbits. Same as number 4, but one planet with freely varying eccentricity. The eccentric, 15 day candidate signal is recovered. The period and RV semi-amplitude of the second planet are poorly constrained, $\mathrm{ln}(Z)=-151$.
• 6.  
  Two planets, both on eccentric orbits. Same as number 4, but with free eccentricity for both planets. Again, the 15 day signal is strongly recovered, while the weak signal of an additional planet is poorly constrained, $\mathrm{ln}(Z)=-144$.

We list all model evidences in Table 4. The log-evidence difference between the preferred model (number 2, single planet) and the runner-up (number 6, two planets, both on eccentric orbits) ${\rm{\Delta }}\mathrm{ln}{Z}_{\mathrm{6,2}}\approx 9$, which corresponds to a Bayes factor of $\sim {10}^{4}$. The difference to the flat model is as large as ${\rm{\Delta }}\mathrm{ln}{Z}_{\mathrm{1,2}}\approx 26$, implying a Bayes factor of $\sim {10}^{11}$. Thus, the planetary model is strongly favored above the flat model and an eccentric single-planet solution is preferred.

Zoom In Zoom Out Reset image size

Table 3. Prior Parameter Distributions

Parameter name	Prior	Units	Description
Stellar Parameters
ρ⋆	${ \mathcal N }(1120,{110}^{2})$	kg m−3	Stellar density
Planetary Parameters
P	${ \mathcal N }(15.16,{0.2}^{2})$	day	Period
t0	${ \mathcal N }(2458319.17,{0.2}^{2})$	day	Time of transit center
${R}_{{\rm{P}}}/{R}_{\star }$	${ \mathcal U }(0.0,1.5)$	⋯	Impact factor
$b=(a/{R}_{\star })\cos (i)$	${ \mathcal U }(0.0,0.5)$	⋯	Planet-to-star ratio
K	${ \mathcal U }(140.0,260.0)$	m s−1	RV semi-amplitude
${S}_{1}=\sqrt{e}\sin \omega$	${ \mathcal U }(-1,1)$	⋯	Parameterization for e and ω
${S}_{2}=\sqrt{e}\cos \omega$	${ \mathcal U }(-1,1)$	⋯	Parameterization for e and ω
RV Instrumental Parameters
${\mu }_{\mathrm{FEROS}}$	${ \mathcal U }(-30,30)$	m s−1	Systemic velocity for FEROS
${\sigma }_{\mathrm{FEROS}}$	${ \mathcal J }(1.0,{100.0}^{2})$	ppm	Extra jitter term for FEROS
${\mathrm{RV}}_{\mathrm{linear}}$	${ \mathcal N }(0.0,{1.0}^{2})$	m s−1 day−1	Linear term for the RVs a
${\mathrm{RV}}_{\mathrm{intercept}}$	${ \mathcal N }(0.0,{10000}^{2})$	m s−1	Intercept term for the RVs a
Photometry Instrumental Parameters
${D}_{\mathrm{TESS}}$	1.0 (fixed)	⋯	Dilution factor for TESS
${M}_{\mathrm{TESS}}$	${ \mathcal N }(0.0,{0.1}^{2})$	ppm	Relative flux offset for TESS
${\sigma }_{\mathrm{TESS}}$	${ \mathcal J }({10}^{-5},{({10}^{5})}^{2})$	ppm	Extra jitter term for TESS
${q}_{1,\mathrm{TESS}}$	${ \mathcal U }(0.0,1.0)$	⋯	Linear limb-darkening parameterization
${q}_{2,\mathrm{TESS}}$	${ \mathcal U }(0.0,1.0)$	⋯	Quadratic limb-darkening parameterization
${D}_{\mathrm{CHAT}}$	1.0 (fixed)	⋯	Dilution factor for CHAT
${M}_{\mathrm{CHAT}}$	${ \mathcal N }(0.0,{0.1}^{2})$	ppm	Relative flux offset for CHAT
${\sigma }_{\mathrm{CHAT}}$	${ \mathcal J }({10}^{-5},{({10}^{5})}^{2})$	ppm	Extra jitter term for CHAT
${q}_{1,\mathrm{CHAT}}$	${ \mathcal U }(0.0,1.0)$	⋯	Linear limb-darkening parameterization
${D}_{\mathrm{LCOGT}}$	1.0 (fixed)	⋯	Dilution factor for LCOGT
${M}_{\mathrm{LCOGT}}$	${ \mathcal N }(0.0,{0.1}^{2})$	ppm	Relative flux offset for LCOGT
${\sigma }_{\mathrm{LCOGT}}$	${ \mathcal J }({10}^{-5},{({10}^{5})}^{2})$	ppm	Extra jitter term for LCOGT
${q}_{1,\mathrm{LCOGT}}$	${ \mathcal U }(0.0,1.0)$	⋯	Linear limb-darkening parameterization
Additional Parameters
${\sigma }_{\mathrm{TESS}}^{\mathrm{GP}}$	${ \mathcal J }({10}^{-8},{0.0005}^{2})$	⋯	Amplitude of the GP component
${\tau }_{\mathrm{TESS}}^{\mathrm{GP}}$	${ \mathcal J }(0.0001,{2}^{2})$	⋯	Timescale of the GP component

Notes. N(μ, σ) stands for a normal distribution with mean μ and standard deviation σ, $U(a,b)$ stands for a uniform distribution between a and b, and $J(a,b)$ stands for a Jeffrey's prior (that is, a log-uniform distribution) defined between a and b.

^a These parameters were only used to search for an additional linear RV trend and are not included in our nominal joint fit.

Download table as:  ASCIITypeset image

Table 4. Model Evidences from RV Fits for Different Models

Model	log-evidence ln Z	Δ ln Z
0 planets	−160.94 ± 0.14	−26.41
1 planet, eccentric orbit	−134.53 ± 0.02	0
1 planet, circular orbits	−160.46 ± 0.04	−25.93
2 planets, circular orbits	−164.34 ± 0.02	−29.81
2 planets, circular and eccentric orbits	−150.76 ± 0.01	−16.23
2 planets, eccentric orbits	−144.46 ± 0.01	−9.93

Note. Δ ln Z states the difference in log-evidence compared to our best model "one planet, eccentric orbit."

Download table as:  ASCIITypeset image

To test if the candidate signal could potentially be associated with stellar activity, we produced generalized Lomb–Scargle periodograms (Zechmeister & Kürster 2009) for the RV time series, as well as for common activity indicators based on FEROS data (Figure 3). In particular, we obtained the H_α , log(${R}_{\mathrm{HK}}^{{\prime} }$), Na II, and He I activity indices, which trace chromospheric activity. We computed ${{\rm{H}}}_{\alpha }$ following Boisse et al. (2009). As TIC 237913194 is a G-type star, we used the regions defined by Duncan et al. (1991) and the calibrations of Noyes et al. (1984) for log(${R}_{{HK}}^{{\prime} }$). For Na II and He I we followed Gomes da Silva et al. (2011). For each time series (RV, H_α , log(${R}_{{HK}}^{{\prime} }$), Na II, and He I), we computed the power levels for 1%, 0.5%, and 0.01% false-alarm probabilities (FAP) by a bootstrap method and plot them as solid, dashed, and dotted lines, respectively. There is a strong signal in the periodogram of the RVs at the 15.17 days period, below 0.5% FAP. Meanwhile, there are no significant signals visible in the periodograms of any of the activity indices, indicating that the RV signal is unlikely to come from quasi-periodic stellar activity.

To further show beyond doubt that the measured RV variations represent reflex motions of the star, we also tested if they can be caused by variations in the stellar photosphere. A well-established method to do this is the inspection of atmospheric line profiles, which should be constant in time for actual stellar velocity changes. Specifically, the bisector span (BIS) can serve as a diagnostic to search for possible false-positive scenarios (e.g., Queloz et al. 2001). We are interested in its correlation with the RV time series and orbital phase, and confront these variables in the right panel of Figure 2. Here, we plot the FEROS RV measurements against their bisector spans. The points are color coded by the orbital phase of the measurements, where zero phase is at T₀ + n P using median values from our nominal fit. While bisector span and RV show no evidence for correlation, this may not be true for bisector span and orbital phase. A Spearman's rank coefficient of ${0.45}_{-0.39}^{+0.31}$, where we quote 95% confidence intervals from a bootstrapped sampling method, ¹⁹ permits the suspicion of a positive correlation. However, due to the small number of data points we cannot reject the null hypothesis that there is no monotonic association between the two variables. In addition, the bisector span variations are on the order of 10 m s⁻¹ and cannot account for the observed RV semi-amplitude of $K\approx 191\,{\rm{m}}\,{{\rm{s}}}^{-1}$. The line profiles thus provide further evidence that the observed RV variations are indeed due to velocity changes of the target star and not caused by atmospheric variations.

In addition to the tests described above, the non-sinusoidal pattern of the RVs is a strong indication for orbital motion as opposed to stellar activity. We conclude that the RV time series independently confirms the planet hypothesis. In the following, we refer to the confirmed exoplanet as TIC 237913194b. ²⁰

3.2.4. Joint Photometry and RV Fit

We performed a simultaneous fit on photometric and spectroscopic data sets (TESS, CHAT, LCOGT, and FEROS) to jointly constrain all planetary parameters of TIC 237913194b. To account for the long cadence in the TESS light curve, we modeled the transits with twenty-fold supersampling using the exposure time of the actual observations. The initial photometric fit (see Section 2.1) provided narrow constraints on the orbital period P and time of mid-transit T₀. We used the median values obtained there to construct Gaussian priors for these parameters, but we enlarged the dispersions (compare Table 3). For the instrument-specific flux offsets ${M}_{\mathrm{TESS}}$, ${M}_{\mathrm{CHAT}}$, and ${M}_{\mathrm{LCOGT}}$ we assumed Gaussian priors based on our photometric fit. A Gaussian prior on the stellar density was motivated by the analysis of the stellar parameters presented in Section 3.1. Additional confidence for this prior stems from a separate joint fit where we used an uninformative prior (${ \mathcal J }({10}^{2},{({10}^{4})}^{2})$) on ρ_⋆. The resulting posterior probability, ${\rho }_{\star }={996}_{-421}^{+257}$, and the result from our stellar analysis agree within the uncertainties. For all other parameters, we chose uninformative priors to sample the whole physically plausible parameter space.

There are potentially time-correlated processes such as instrumental red noise, stellar variability, or blended sources that are not covered by our astrophysical model. We account for this red noise by adding a GP component to the TESS photometry with an exponential kernel as implemented in the celerite software package (Foreman-Mackey et al. 2017). This adds two additional hyperparameters: an amplitude ${\sigma }_{\mathrm{TESS}}^{\mathrm{GP}}$ and a timescale ${\tau }_{\mathrm{TESS}}^{\mathrm{GP}}$. For comparison, we performed the same fit with and without a GP component. The variant including GP performed significantly better than the white-noise model and we thus consider it our nominal model.

In the same manner, we tested adding a GP component to the LCOGT photometry, which shows possible systematic effects in the residuals (compare Figure 1). Here, we chose a Matérn 3/2 kernel, which again adds two hyperparameters for timescale and amplitude to the model. The model including the GP component consistently performed comparable ($| {\rm{\Delta }}\mathrm{ln}Z| \lt 1$) or worse than the model without, which is why we chose to continue with the less complex noise model without GP.

Figures 1 and 2 show the photometric and RV time series resulting from this model using its median parameters (solid lines). Dark (light) blue bands show the 68% (95%) credibility bands of the model. The residuals below the time series show the measured data with the median posterior model subtracted; both RV and photometry residuals appear inconspicuous. Figure 4 shows the posterior distributions for the planet's main parameters as sampled in our nominal fit. All distributions are approximately Gaussian and barely correlated, except for the planet-to-star ratio ${R}_{{\rm{P}}}/{R}_{\star }$ (see Appendix A for a discussion). We present the posterior distributions of the model parameters alongside the derived physical parameters in Table 5, where we state for each parameter distribution the 16th, 50th, and 84th percentile. Notably, with a planetary bulk density of ${\rho }_{{\rm{P}}}\,\approx \,1847\,\,\,\mathrm{kg}\ \,\ {{\rm{m}}}^{-3}$, TIC 237913194b's average density is comparable to Neptune's. By sheer coincidence, the planet's period and eccentricity resemble that of the TESS spacecraft (13.7 days, e = 0.55; Ricker et al. 2014).

Zoom In Zoom Out Reset image size

Table 5. Posterior Parameters

Parameter	Value
Stellar Parameters
${\rho }_{\star }\,(\mathrm{kg}\,{{\rm{m}}}^{-3})$	${1076}_{-93}^{+95}$
Planetary Parameters
P (day)	${15.168865}_{-0.000018}^{+0.000018}$
t0 (BJD UTC)	${2458319.15055}_{-0.00077}^{+0.00077}$
a/R⋆	${23.85}_{-0.69}^{+0.67}$
$b=(a/{R}_{\star })\cos (i)$	${0.900}_{-0.017}^{+0.017}$
${R}_{{\rm{P}}}/{R}_{\star }$	${0.1031}_{-0.0042}^{+0.0048}$
K (m s−1)	${191.3}_{-6.2}^{+6.4}$
e	${0.575}_{-0.011}^{+0.011}$
ω	${24.1}_{-2.4}^{+2.4}$
${S}_{1}=\sqrt{e}\sin \omega$	${0.309}_{-0.030}^{+0.029}$
${S}_{2}=\sqrt{e}\cos \omega$	${0.692}_{-0.015}^{+0.014}$
RV Instrumental Parameters
${\mu }_{\mathrm{FEROS}}$ (m s−1)	${14.0}_{-3.3}^{+3.4}$
${\sigma }_{\mathrm{FEROS}}$ (ms−1)	${13.6}_{-2.6}^{+3.2}$
Photometry Instrumental Parameters
${M}_{\mathrm{TESS}}$ (ppm)	${0.00017}_{-0.00039}^{+0.00041}$
${\sigma }_{\mathrm{TESS}}$ (ppm)	${221}_{-31}^{+28}$
${q}_{1,\mathrm{TESS}}$	${0.33}_{-0.24}^{+0.39}$
${q}_{2,\mathrm{TESS}}$	${0.27}_{-0.20}^{+0.35}$
${\sigma }_{\mathrm{TESS}}^{\mathrm{GP}}$	${0.00000152}_{-0.00000041}^{+0.00000110}$
${\tau }_{\mathrm{TESS}}^{\mathrm{GP}}$	${0.40}_{-0.18}^{+0.16}$
${M}_{\mathrm{CHAT}}$ (ppm)	${0.00149}_{-0.00025}^{+0.00025}$
${\sigma }_{\mathrm{CHAT}}$ (ppm)	${1625}_{-140}^{+150}$
${q}_{1,\mathrm{CHAT}}$	${0.55}_{-0.24}^{+0.20}$
${M}_{\mathrm{LCOGT}}$ (ppm)	$-{0.00038}_{-0.00016}^{+0.00016}$
${\sigma }_{\mathrm{LCOGT}}$ (ppm)	${480}_{-80}^{+85}$
${q}_{1,\mathrm{LCOGT}}$	${0.67}_{-0.20}^{+0.17}$
Derived Parameters
i (deg)	${87.0}_{-1.7}^{+1.5}$
${R}_{{\rm{P}}}\,[{R}_{\mathrm{Jup}}]$	${1.117}_{-0.047}^{+0.054}$
${M}_{{\rm{P}}}\,[{M}_{\mathrm{Jup}}]$	${1.942}_{-0.091}^{+0.092}$
$a\,[\mathrm{au}]$	${0.1207}_{-0.0037}^{+0.0037}$
${\rho }_{{\rm{P}}}\,[\mathrm{kg}\ \,\ {{\rm{m}}}^{-3}]$	${1847}_{-260}^{+280}$
${T}_{\mathrm{eq}}\,[{\rm{K}}]$ a	$974$

Note.

^a Time-averaged equilibrium temperature computed according to Equation (16) of Méndez & Rivera-Valentín (2017). We assumed zero albedo, a unity broadband thermal emissivity, and $\beta =0.5$, i.e., only half of the planetary surface re-radiates the absorbed flux.

Download table as:  ASCIITypeset image

3.2.5. Search for Additional Planets

The equilibrium temperature T_eq that the planet maintains if it is in energy balance with the radiation input from the host star is a determining factor for the physical properties of its atmosphere. Due to TIC 237913194b's considerable eccentricity, this input is not constant over its orbit. To give a first-order estimate on the temperature range that the planet can assume, we investigated two extreme cases of planetary heat adjustment:

• 1.  
  Instantaneous heat adjustment (${T}_{\mathrm{eq},\mathrm{inst}}$). Here, we assumed that the planetary atmosphere adjusts to the changing irradiation without any time delay. For this situation, we used the approximation in Kaltenegger & Sasselov (2011, Equation (3)).
• 2.  
  Orbitally averaged heat adjustment (${T}_{\mathrm{eq},\mathrm{avg}}$). In this case, the planetary temperature remains in equilibrium with the incoming stellar energy, i.e., ${T}_{\mathrm{eq}}=\mathrm{const}.$ over one orbit. To approximate this temperature, we used a temporal average for elliptic orbits (Méndez & Rivera-Valentín 2017, Equation (16)).

For both extremes, we assumed that the heat flux from the planet's interior is negligible compared to the stellar irradiation and ignored any internal energy sources. The infrared emissivity ${\varepsilon }_{\mathrm{IR}}$ was fixed to unity. We further assumed two cases for how atmospheric flow distributes the incoming stellar energy over the planetary surface and parameterized this property with the fraction of planetary surface that re-radiates stellar flux β. We distinguished between β = 0.5, i.e., emission only from one hot hemisphere, and β = 1 where the whole globe emits (Seager et al. 2005; Kaltenegger & Sasselov 2011; Carone et al. 2014; Méndez & Rivera-Valentín 2017).

Planets colder than 1000 K are expected to be relatively cloudy (e.g., Parmentier et al. 2016; Stevenson 2016). Here, we parameterized different cloudiness with albedos α = 0, 0.3, and 0.6 following Parmentier et al. (2016). With the above assumptions and in the case of instantaneous heat adjustment, we derived a range of ${T}_{\mathrm{eq},\mathrm{inst}}\approx 900\mbox{--}1300\,{\rm{K}}$ at secondary eclipse and ${T}_{\mathrm{eq},\mathrm{inst}}\approx 700\mbox{--}1100\,{\rm{K}}$ at transit (Figure 5). We list ${T}_{\mathrm{eq},\mathrm{inst}}$ at critical times in Table 6 together with the values for orbitally averaged heat adjustment ${T}_{\mathrm{eq},\mathrm{avg}}$. The latter is constant over one orbit and covers a temperature range of ${T}_{\mathrm{eq},\mathrm{avg}}\approx 650\mbox{--}975\,{\rm{K}}$.

Zoom In Zoom Out Reset image size

Table 6. Theoretical Temperature Constraints of TIC 237913194b in the Course of One Orbit

Time	Orbital Distance	${T}_{\mathrm{eq},\mathrm{inst}}\,$(K)			${T}_{\mathrm{eq},\mathrm{avg}}\,$(K)
	(au)	α = 0	α = 0.3	α = 0.6	α = 0	α = 0.3	α = 0.6
β = 1
Apoastron	0.1900	668	611	531	819	749	651
Transit	0.1047	900	823	716
Periastron	0.0520	1277	1168	1016
Secondary eclipse	0.0668	1127	1031	896
β = 0.5
Apoastron	0.1900	795	727	632	974	891	775
Transit	0.1047	1070	979	851
Periastron	0.0520	1519	1389	1208
Secondary eclipse	0.0668	1340	1226	1066

Download table as:  ASCIITypeset image

Due to the orientation of the orbit (compare Figure 6), the temperature ${T}_{\mathrm{eq},\mathrm{inst}}$ during transit is assumed to be about 200 K colder compared to the temperature at secondary eclipse. In reality, however, some delay in heat adjustment based on radiative and dynamical timescales is expected. Therefore, the temperature during a secondary eclipse, which would occur before passage of periastron, could be colder than in our estimate. Likewise, the temperature at transit, occurring after periastron passage, would be warmer than expected (see, e.g., Lewis et al. 2013, for a qualitative discussion of the thermal evolution of an exoplanet on an eccentric orbit). We emphasize again that our goal is to estimate to the first order the possible temperature ranges for TIC 237913194b for which these simplified assumptions are sufficient.

Zoom In Zoom Out Reset image size

No secondary eclipses or phase curve signals are evident in the photometric time series. We used the starry software package (Luger et al. 2019) to estimate the planetary phase curve based on our derived orbital parameters and a simple toy model for the planetary brightness distribution. This model neglects any heat redistribution between the hot and cold hemispheres of the planet, resulting in a dipole brightness map where the hot side points to the substellar point at periastron. In this scenario, the emission of the cold hemisphere and planetary limb-darkening are negligible. The total luminosity of the planet is then that of a half-sphere blackbody with radius ${R}_{{\rm{P}}}$ and temperature T_hot. For the temperature of the hot hemisphere, we adopt two cases: first, our estimate of the equilibrium temperature for the case of orbitally averaged heat adjustment and $\beta =0.5$, hence ${T}_{\mathrm{hot}}={T}_{\mathrm{eq},\mathrm{avg}}=974\,{\rm{K}}$. Second, we consider the hottest temperature in Table 6 and assume ${T}_{\mathrm{hot}}=1519\,{\rm{K}}$. With a resulting peak-to-peak phase curve amplitude of ∼6 ppm in the cool case and ∼32 ppm in the hot case, a future detection of the phase curve or secondary eclipse might be challenging despite the expected precision of the James Webb Space Telescope (JWST; Beichman et al. 2014).

Measuring the Rossiter–McLaughlin effect (RM effect; McLaughlin 1924; Rossiter 1924) has proved a useful tool to measure the alignment of planetary orbits with the spin axis of host stars. The different proposed scenarios for the formation and migration theory of warm Jupiters differ in their predicted impact on the spin–orbit alignment (e.g., Triaud 2018). A detection of the RM effect could thus shed light on TIC 237913194b's enigmatic formation history. Analytical formulas exist to estimate the amplitude of its RV signature (e.g., Gaudi & Winn 2007, Equation (6)), however, the large impact parameter in the transit geometry of TIC 237913194b would entail a large error. Instead, we modeled the RV anomaly due to the RM effect with a velocity-weighted brightness map in starry using the median posterior values of the system's physical parameters (see Tables 2 and 5).

The resulting amplitude of the signal is ${K}_{\mathrm{RM}}\approx 10\,{\rm{m}}\,{{\rm{s}}}^{-1}$, which is just in the range of current state-of-the-art spectroscopic facilities.

In Figure 7 we compare TIC 237913194b with well-studied transiting exoplanets (Southworth 2011). ²¹ The left panel shows the periods and eccentricities of these planets (blue markers); our discovery is marked in red. We included planets with both mass and radius measurements that have constrained eccentricities (not only upper limits) and show those in the period range of 1–100 days. Marker sizes in the plot correspond to planet masses. Warm Jupiters with high eccentricities appear to be rare: 98% of this population have smaller eccentricities than TIC 237913194b, making it one of the most eccentric planets in this period range. It lies at the edge of a demographic feature that we discuss in the following.

Zoom In Zoom Out Reset image size

On close orbits, the planet occurrence rate ${\rm{d}}n/{\rm{d}}e$ shows a rapid drop with increasing eccentricity. With increasing period, the position of this ridge shifts to larger eccentricities. Through this, a triangular under-density of planets with high eccentricity on very short orbits emerges (upper left corner in Figure 7). While exoplanet detection sensitivities are expected to have a dependency on eccentricity, the effect is too small to account for the observed dearth of planets (Burke 2008). A plausible physical explanation would be tidal circularization (Adams & Laughlin 2006; Dawson & Johnson 2018). As discussed in Section 4.3, the strength of this mechanism is expected to scale inversely with orbital distance, which would explain the period-dependence of the distribution. The detection of planets close to this ridge can thus help constrain theories of tidal interaction and giant planet migration, which are crucial components for explaining planetary systems with close-in giant planets. Our discovery of TIC 237913194b adds to the small current sample of such planets.

In the right panel of Figure 7, we put our planet into the context of warm Jupiters with mass and radius measurements. Here, we include only planets from TEPCat with periods of 10–100 days, and color code them by equilibrium temperature. We further plot a theoretically predicted mass–radius relation for planets with a 10 ${M}_{\oplus }$ core (Fortney et al. 2007; blue line). TIC 237913194b is located close to this curve, which indicates that its bulk density is consistent with established structural models.

The large eccentricity of TIC 237913194b makes it a promising test bed to study the response of its atmosphere to external forcing (e.g., compare Carone et al. 2020). Atmospheres at the inferred temperature ranges are susceptible to a variety of chemical disequilibrium processes such as photochemistry and chemical quenching (e.g., Venot et al. 2012; Moses et al. 2013; Tsai et al. 2018; Kawashima & Ikoma 2019; Molaverdikhani et al. 2019a; Venot et al. 2020). While a thorough inspection of these processes is beyond the scope of this paper, we used the physical parameters constrained here to demonstrate the feasibility of atmospheric characterization. To this end, we used self-consistent models for cloud-free (Molaverdikhani et al. 2019b) and cloudy (Molaverdikhani et al. 2020) irradiated planetary atmospheres and calculated synthetic transmission and emission spectra using petitCODE (Mollière et al. 2015, 2017). For this fiducial model, we assumed solar composition, zero bond albedo, instantaneous thermal equilibrium, and β = 0.5. This resulted in an equilibrium temperature of the planet during transit and a secondary eclipse of 1070 K and 1340 K, respectively (see Table 6). Using the Pandexo package (Batalha et al. 2017), we predicted uncertainties for JWST observations of a single transit or secondary eclipse for the three observing modes: Mid-Infrared Instrument (MIRI) low-resolution spectroscopy (LRS), Near Infrared Imager and Slitless Spectrograph (NIRISS) single object slitless spectroscopy (SOSS), and NIRSpec-G395M.

We find relative magnitudes of the largest synthetic spectral features of ∼100 ppm in transmission and ∼1000 ppm in emission. While some emission features are well above the predicted noise floor, the largest transmission features are on the order of the predicted uncertainties for a single transit observation. However, since these two techniques probe different regions of the atmosphere and at different orbital phases, a joint analysis of the transmission and emission spectra may provide important clues on atmospheric dynamics and heat redistribution.

The peculiarly large eccentricity of TIC 237913194b could be an important lead in not only understanding the dynamical origin of the system, but also planet evolution theories in general. Possible origins of large warm Jupiter eccentricities include interaction with a massive companion through either scattering events (e.g., Rasio & Ford 1996), secular interactions (e.g., Kozai 1962; Lidov 1962; Petrovich & Tremaine 2016), or giant impacts (Frelikh et al. 2019); planet–disk interactions (e.g., Lubow 1991; Petrovich et al. 2019) or a combination of processes. The absence of evidence for an additional perturber that might sustainably excite TIC 237913194b's eccentricity or that could be the counterpart in a recent scattering event that makes it challenging to distinguish between these scenarios.

However, TIC 237913194b is subject to tidal dissipation through secular interaction with the host star (Goldreich & Soter 1966) and the rates of semimajor axis and eccentricity decay, given the current orbital parameters, can be estimated (Yoder & Peale 1981). A short orbit circularization timescale compared to the lifetime of the star would exclude the planet–disk interaction scenario and could provide an upper limit on the time that has passed since a hypothetical perturbation event. In this case, we would observe the system during high-eccentricity migration and TIC 237913194b would eventually become a hot Jupiter in a circular orbit.

Several caveats have to be considered when trying to trace the system back in time close to its primordial orbital configuration. First, the tidal evolution of a and e are strongly coupled, which may result in ambiguities. In addition, the tidal evolution strongly depends on the stellar and planetary tidal dissipation rates, typically parameterized by the dimensionless reduced tidal quality factors ${Q}_{{\rm{P}}}^{{\prime} }$ and ${Q}_{\star }^{{\prime} }$. Here, ${Q}^{{\prime} }=1.5\,Q/{k}_{2,}$, where k_2, is the Love number of the second order. Estimates of the planetary tidal quality factor range from ${Q}_{{\rm{P}}}^{{\prime} }={10}^{4}$ to ${Q}_{{\rm{P}}}^{{\prime} }={10}^{7}$ (e.g., Lainey et al. 2009; Hansen 2012; Lainey et al. 2020). The stellar tidal dissipation factor is even less well constrained but theoretical and observational works suggest ${Q}_{\star }^{{\prime} }\gtrsim {10}^{7}$ (e.g., Carone & Pätzold 2007; Hansen 2012; Damiani & Lanza 2015).

We studied the star–planet tides of TIC 237913194b using the EqTide ²² code (Barnes 2017), which calculates the tidal evolution of two bodies based on models by Ferraz-Mello et al. (2008) and a constant-phase-lag (CPL) model (Goldreich & Soter 1966; Cheng et al. 2014). For our tidal-torque test, we adopted a stellar rotational period of 30 days and an initial planetary rotational period of 0.5 day (i.e., similar to the solar and Jupiter rotational periods). We adopted tidal quality factors ${Q}_{{\rm{P}}}^{{\prime} }$ in the range of 3 × 10⁴–10⁶. For the primary, we chose a fixed value of ${Q}_{\star }^{{\prime} }$ = 10⁸ (e.g., Hansen 2010; Penev et al. 2012; Bonomo et al. 2017). We tested a large set of increased initial semimajor axes and eccentricities and integrated with EqTide. The results agree with the observed eccentricity and semimajor axis only for those samples that started a few percent above the observed values. This suggests that the orbital period of ∼15 days is too large for significant tidal circularization within the age of the system (∼5.5 Gyr, see Table 2) and TIC 237913194b's orbit has only slightly evolved from its primordial configuration. These results are in line with Barnes (2015), who showed that Jovians with periods longer than ∼8 days and a typical eccentricity of 0.3 do not experience significant orbital and eccentricity decay. While we cannot determine the origin of the high orbital eccentricity, we conclude that the planet we observe today is not a credible progenitor of a future hot Jupiter.